基于網格自適應的飛行器防熱材料熱傳導系數辨識

章 勝，周 宇，錢煒祺，何開鋒

(1. 空氣動力學國家重點實驗室，綿陽 621000；2. 中國空氣動力研究與發展中心計算空氣動力研究所，綿陽 621000)

0 引 言

高超聲速飛行器在飛行過程中會遇到復雜的氣動加熱現象，其對防隔熱材料技術提出了嚴苛的要求，發展滿足設計熱性能要求的材料是飛行器完成既定任務的前提。氣動熱辨識是利用飛行器表面或防熱層內部的溫度測量結果來反演飛行器表面熱流、防熱材料熱物性參數等氣動熱相關物理量的技術，其研究包括飛行器表面熱流時間變化歷程的辨識、防熱材料熱物性參數隨溫度變化規律的辨識與防熱材料熱源項的辨識等[1]，這些研究對促進高超聲速飛行器的發展具有重要意義[2]。

氣動熱辨識問題又可稱為熱傳導反演問題，其是熱傳導物理問題的逆問題，它在數學上的描述是針對偏微分方程系統的優化問題。近年來，國外學者在氣動熱辨識技術的理論研究、數值計算與實驗驗證等方面開展了大量深入的研究工作并在工程實際中得到了應用[3]。材料熱傳導系數的辨識是典型的熱傳導逆問題，熱傳導系數一般是溫度的連續函數[4]，其辨識問題是無限維的約束泛函極值問題。對于此類問題，一般是通過參數化離散手段將其轉化為有限維優化問題進行求解。目前，工程上材料熱傳導系數辨識中通常假設熱傳導系數為常值或隨溫度的變化規律確定，在此前提下進一步采用穩態法、激光法、靈敏度法等方法進行參數反演[1-3]。為了更準確地辨識熱傳導系數，獲得其隨溫度的變化規律，學者們將熱傳導系數值按溫度區間分段離散，采用折半法[5]、復變量求導法[6]、遺傳算法[7]等方法辨識各溫度段的數值。由于反演問題是數學上的不適定問題，當測量誤差較大或反演數學模型和真實物理模型存在明顯差異時，辨識結果中會出現較大的數值振蕩，甚至可能出現熱傳導系數值小于零的違背物理實際的情況。因此，在熱傳導系數辨識中有必要引入相應的物理約束，從而得到更可信的結果。文獻[8]在反演材料熱傳導系數的遺傳算法基礎上，考慮熱傳導系數值大于零、熱傳導系數值隨溫度增加而增加的物理約束，建立了熱傳導系數辨識方法并進行了算例考核。

從數據學習的角度講，辨識問題也可以納入機器學習的范疇，由統計學習理論可知：測量樣本固定時，模型的復雜度越高，則模型的學習能力也越強，但過于復雜的模型通常會使得泛化性能降低[9]。對于觀測數據比較有限同時噪聲特性不清楚的辨識問題，研究中需要考慮在數據擬合誤差與模型泛化能力之間進行折中，避免將干擾噪聲誤認為物理量變化規律的“過擬合”問題。另一方面，對于參數化離散得到的非線性規劃問題(Non-linear programming, NLP)，不同于常用的迭代優化求解方法[10-11]，基于動態優化方程(Dynamic optimization equation, DOE)的方法是一種將NLP轉化為初值問題(Initial-value problem, IVP)加以求解的方法[12-19]。雖然該類方法的研究最早可以溯源到二十世紀40年代[12]，但是其并沒有得到廣泛的關注，由于未知參數的尋優過程可以通過一定的動力學方程描述，其可能更適宜于在線參數辨識與控制器設計(作為增廣狀態)。針對考慮物理約束的高超聲速飛行器防熱材料熱傳導系數辨識問題，本文將采用多項式參數化建模技術，借鑒機器學習領域中防止“過擬合”的樣本驗證手段設置網格自適應調節準則，對基于DOE的優化方法在氣動熱辨識領域的應用進行探索，實現材料熱傳導系數的準確辨識。

1 材料熱傳導系數辨識模型

高超聲速飛行器熱防護層一般為多層多孔隙結構，在一定的合理假設下，防熱層的三維熱傳導問題可以近似為一維問題進行研究[20-21]。為確定防熱材料的熱傳導物性，工程中可對其不同方向的特性分別加以研究，確定材料熱傳導系數的實驗原理圖如圖1所示，在材料上下邊界為絕熱邊界條件的假設下，材料特定方向的傳熱問題可以簡化為一維傳熱問題，當邊界條件和觀測條件為下述兩種情形之一時，可以進行材料熱傳導系數的反演[8]。

1)材料內部無測點溫度測量的情況：左右邊界分別給定溫度或熱流量，將其中一個邊界上除邊界條件給定的狀態量外的另一狀態量作為觀測量。

2)材料內部有測點溫度測量的情況：左右邊界分別給定溫度或熱流量，將內部測點的溫度作為觀測量。

圖1 反演材料熱傳導系數的實驗原理示意圖Fig.1 Sketch of the experimental principle for thermal conductivity inversion

本文針對某種典型防熱材料，對第二種情形下的材料熱傳導系數辨識進行研究，在材料左右邊界給定溫度邊界條件，在材料內部選取一個或多個測點進行溫度測量來對材料的熱物性參數進行反演。材料的傳熱方程為

(1)

式中:T(x,t)為溫度;ρ為材料密度;Cp為材料比熱;k(T)為材料的熱傳導系數，它是溫度T的函數。方程的邊界條件為

(2)

初始條件為

t=0∶T(x,0)=Tt0

(3)

觀測方程為

(4)

基于邊界條件式(2)與初始條件式(3)，給定熱傳導系數k(T)，可以采用有限差分法[22]實現熱傳導正問題的數值求解。對于逆問題，此時熱傳導系數k(T)未知，需要通過觀測方程式(4)中的信息對其進行辨識，在數學上該逆問題可表示為尋求k(T)使得如下性能指標達到極小

(5)

0≤k(T)

(6)

(7)

綜上，熱傳導系數辨識問題可整理為如下約束泛函極值問題：

s.t.

x=0∶T(0,t)=Tx0

x=L∶T(L,t)=TxL

t=0∶T(x,0)=Tt0

0≤k(T)

(8)

2 辨識方法

2.1 參數化離散

參數化離散是求解形如式(8)類型的約束泛函極值問題的有效方法。所有的數值方法都是通過參數化離散手段將無窮維問題轉化為有限維問題加以解決的。在眾多離散方法中，基于多項式的參數化離散是一種重要的手段，其理論基礎是Weierstrass于1885年提出的多項式逼近定理[23]。事實上，Weie-rstrass定理正是當前流行的求解最優控制問題的偽譜方法的理論基礎[24-25]。下面首先給出該定理。

(9)

式中:上標^指代建模結果;φi(T)為Lagrange插值多項式，為

(10)

顯然φi(T)滿足

(11)

相應的，θi(i=1,2,…,n)的物理意義為溫度Ti(i=1,2,…,n)處的熱傳導系數值，分析可知熱傳導系數模型式(9)為n-1階多項式。

另一方面，針對物理約束式(6)和式(7)，在溫度區間[Ta,Tb]取N個節點(同樣T1=Ta,TN=Tb，但N?n)，則物理約束可以表示為

(12)

采用這種方式，即使對于熱傳導參數物性復雜變化的情形，文中方法也易于處理。

s.t.

x=0∶T(0,t)=Tx0

x=L∶T(L,t)=TxL

t=0∶T(x,0)=Tt0

0≤k(T1)

k(Tj)≤k(Tj+1)j=1,2,…,N-1

(13)

2.2 自適應網格迭代算法

為實現對未知連續變量的準確求解，數值計算中通常需要對離散網格進行加密、實施多次參數化離散以提高解的精度[26-27]。為改善熱傳導系數辨識結果，文中采用自適應網格迭代算法，通過調節熱傳導系數模型式(9)的階次實現對真實熱傳導系數的逼近。熱傳導系數辨識問題基于對溫度信號的有限含噪觀測，通過使重構溫度信息與真實信號間誤差最小建立熱傳導系數辨識模型。在溫度測量噪聲特性未知的條件下，該誤差的期望值(即期望風險)無法準確求解，根據統計學習理論中常用的經驗風險最小化原則[28]，采用經驗風險(文中即式(5)定義的優化指標J)作為對期望風險的估計并使之最小。但是，經驗風險最小化原則是建立在測量信息足夠充分的前提之下，當測量樣本比較有限時(尤其是通過飛行試驗獲得的數據)，經驗風險最小不一定代表期望風險最小，反而復雜的學習模型(在本文中，相應于多項式模型式(9)階次增加)會導致結構風險，即模型泛化性能下降，產生“過擬合”問題。為解決這一問題，參考機器學習領域的典型做法：采用一定時域TV上的測量樣本對辨識結果進行檢驗，定義“過擬合”檢驗指標

(14)

網格自適應迭代算法具體實現為：

為提高計算效率，從第1次網格迭代開始，NLP問題式(13)中參數θ通過上一次的辨識結果進行初始化，即

(15)

2.3 基于DOE的NLP求解方法

針對參數化離散得到的NLP，學者們已發展了大量的求解方法，如序列二次規劃方法[10]、內點方法[11]等，這些方法通過離散迭代可以對NLP進行有效求解，但它們的實現往往比較繁瑣。基于DOE的優化方法是一種將NLP轉化為IVP加以求解的方法，這類方法具有簡潔的動力學公式表達，便于理解與使用。對于無約束的優化問題，學者們提出了梯度動力學方程[13]與連續牛頓方程[14]等DOE。為解決包含等式約束與不等式約束的一般NLP，文獻[15-16]建立了求解含等式約束NLP的DOE。通過松弛參數技術，文獻[17-18]將不等式約束轉化為等式約束，給出了可處理不等式約束的DOE。松弛參數的引入增加了DOE的維度，對于含有大量不等式約束的問題,這會引入過多的待求參數，增加問題復雜度。文獻[19]利用不等式約束的動力學特性直接對不等式約束進行處理，推導的DOE中避免了松弛參數的引入。

對于動力學優化方法得到的IVP，可以通過成熟的常微分方程積分方法進行計算，易于實現數值求解。尤其DOE對最優解的求解有理論保證，能嚴格保證解的全局收斂性，而且還可以得到經典最優性條件中Lagrange乘子與Karush-Kuhn-Tucker(KKT)乘子的解析表達式，直觀地揭示了乘子參數與優化參數之間的理論關系[19]。

考慮一般形式的NLP，性能指標為

minJ=f(θ)

(16)

受到如下的不等式與等式約束：

g(θ)≤0

(17)

h(θ)=0

(18)

式中:θ∈Rn為待優化參數;f:Rn→R為一階導數連續的標量函數;g:Rn→Rp、h:Rn→Rq為一階偏導數連續的矢量函數。為建立求解NLP的DOE，文獻[19]定義了如下達可行性動力學優化問題，這是一個二次規劃問題，為

s.t.

(19)

式中：fθ與(gi)θ分別為函數f與gi對參數θ的梯度(列)矢量，hθ代表Jacobi矩陣，Kθ、Kh與kg分別為正定矩陣與正常數，τ為描述參數演化的虛擬時間維度，指標集Ⅱ為

Ⅱ={i|gi≥0,i=1,2,…,p}

(20)

定理2[19].對于上述定義的NLP，給定任意初始條件θτ=0，求解如下DOE

(21)

式中:Lagrange乘子參數πE∈Rq與KKT乘子參數πI∈Rp通過下式確定

(22)

(23)

(24)

(25)

i=1,2,…,n

(26)

式(26)中需要確定溫度T對θi(i=1,2,…,n)的梯度，它們可以通過求解下述靈敏度方程獲得

i=1,2,…,n

(27)

相應的邊界條件為

(28)

初始條件為

(29)

對于該靈敏度偏微分方程組，同樣可以采用有限差分方法進行計算。對于NLP問題式(13)中的N個約束，由于熱傳導系數對參數θ的梯度為

(30)

基于上式，容易求得約束函數對參數θ的梯度。

3 算例分析

算例中研究的防熱材料的參數為ρ=1000 kg·m-3，Cp=1000 J·kg-1·K-1，材料長為L=10 mm，左右邊界的溫度條件如圖2中“x=0 mm”,“x=10 mm”所示，為

初始條件為

T(x,0)=600 K

測點數為m=2個，它們分別位于“x=3 mm”和“x=6 mm”位置。假設材料的熱傳導系數隨溫度的變化規律為

基于熱傳導方程可計算出兩個測點的溫度歷程，圖2給出了無測量噪聲情況下“x=3 mm”和“x=6 mm”處的溫度變化曲線。

圖2 溫度邊界條件及測點理想溫度曲線Fig.2 Temperature boundary conditions and noiseless measurement curves

3.1 無測量噪聲情形

采用文中的辨識方法，第4次網格迭代后“過擬合”檢驗指標增加，因此最終辨識熱傳導系數采用第3次網格迭代后的計算結果，表1給出了網格迭代過程中的具體結果，從表1中可以看到，最終辨識的熱傳導系數對應的優化與檢驗指標值均很小，這說明其與真實值很接近。

表1 無測量噪聲情形的辨識結果Table 1 Identified results without measurement noise

圖3給出了初始網格下NLP參數動力學發展曲線，可以看到在設置的Kθ下參數收斂很快，在τ=1 s時結果已趨于穩定，這說明了本文發展的DOE的高效性。圖4給出了第4次迭代網格下NLP參數辨識的過程曲線，最終結果與初始值幾乎完全相同，表明文中基于上一步辨識結果設置參數初值的方式可以提供良好的初始值。根據表1的辨識結果，圖5給出了初始網格下得到的熱傳導系數、最終辨識的熱傳導系數與真實系數的比較曲線，可以看到辨識結果很好地滿足了實際物理約束，但線性近似與真實曲線之間還是存在較明顯的差異，而最終辨識結果與真實曲線幾乎完全重合，最大誤差僅為3.5020×10-5W·m-1·K-1。

圖3 無測量噪聲情形初始網格參數辨識曲線Fig.3 Identification curve for the parameters of the initial mesh without measurement noise

圖4 無測量噪聲情形第4次迭代網格參數辨識曲線Fig.4 Identification curve for the parameters of the 4th refined mesh without measurement noise

圖5 無測量噪聲情形辨識熱傳導系數與真實系數比較Fig.5 Comparison between the identified thermal conductivity without measurement noise and the true quantity

3.2 含測量噪聲情形

進一步考慮測量中存在噪聲的情形，假設噪聲服從正態分布，均值為0，標準差為σ=4 K，辨識計算條件與無噪聲情形完全相同。表2給出了熱傳導系數辨識過程中的具體計算結果，其中第2次網格迭代后的結果對應的“過擬合”檢驗指標最小，說明此時達到了“經驗風險”與“結構風險”的最佳折中。注意不同于無測量噪聲情形，由于噪聲的影響，辨識結果對應的指標值并不接近于0值。

圖6給出了辨識的熱傳導系數與真實系數的比較曲線，從圖中看到存在噪聲時候辨識結果與真實曲線之間有一定的細微差別(最大誤差為0.0165 W·m-1·K-1)，但是其仍然很好地反映了真實物理量的變化規律。為進行比較，圖7給出了第3次迭代網格下的辨識結果，表2中的檢驗指標顯示此時結果出現了“過擬合”現象，從圖中也看到辨識的熱傳導系數與真值間的誤差開始增大(最大誤差為0.0744 W·m-1·K-1)。通過對比，表明本文自適應網格迭代算法“學習”到了更接近于真實物理量的結果。圖8給出了基于辨識的熱傳導系數重構的溫度狀態與測量量的對比結果，可以看到重構結果較好地擬合了測量信息。

表2 含測量噪聲情形的辨識結果Table 2 Identified results with measurement noise

圖6 有測量噪聲情形辨識熱傳導系數與真實系數比較Fig.6 Comparison between the identified thermal conductivity with measurement noise and the true quantity

圖7 有測量噪聲情形“過擬合”熱傳導系數與真實系數比較Fig.7 Comparison between the over-fitted thermal conductivity with measurement noise and the true quantity

圖8 溫度測量結果與重構曲線比較Fig.8 Comparison between the temperature measurement and the reconstructed results

4 結 論

考慮物理約束的高超聲速飛行器防熱材料熱傳導系數辨識是針對偏微分方程系統的約束泛函極值問題，本文對該問題進行了有效求解，辨識結果較準確地反映了熱傳導系數的變化規律。研究表明本文采用的多項式建模手段可以較好地模擬材料熱傳導系數的物理特性，發展的動態優化方程能有效求解最優參數，對于測量數據比較有限的問題，借鑒“機器學習”思想，提出的基于“過擬合”指標的網格自適應迭代方法能實現模型準確性與泛化性能的較好折中。

由于熱傳導系數的模型結構是通過“學習”確定，本文方法相較于之前確定模型結構的辨識方法更具靈活性。基于動態優化方程的非線性規劃方法在工程上的應用較新穎，本文工作也證明了這種方法的有效性。在本文中，三維傳熱問題是簡化為一維問題加以求解，為更準確地獲得材料的熱傳導物性，下一步將針對多維情形的問題開展研究。