考慮約束的材料熱傳導系數辨識方法

周 宇錢煒祺邵元培顏維旭

(1.空氣動力學國家重點實驗室,綿陽 621000；2.中國空氣動力研究與發展中心 計算空氣動力研究所,綿陽 621000；3.中國運載火箭技術研究院,北京 100076)

0 引 言

高速飛行器艙內熱環境、航空發動機和火箭發動機燃燒室承熱能力的有效預測,一方面取決于飛行器表面熱流和燃燒室壁面熱流的準確模擬,另一方面也取決于防隔熱材料熱物性參數的準確確定。目前,工程上確定材料熱傳導系數的方法大多假設材料的熱傳導系數為常數或是隨溫度變化的已知函數,采用穩態法、激光法、靈敏度法等方法進行參數辨識[1-8]。對于熱傳導系數隨溫度變化規律未知的情況,文獻[9-12]將材料的熱傳導系數值按溫度區間分段離散,并通過材料內部或邊界溫度場的測量結果采用折半法、復變量求導法、遺傳算法等方法來辨識出各溫度段的熱傳導系數值。但是,在工程實際應用中發現,由于辨識問題是數學上的不適定問題,當測量誤差較大或是計算數學模型和真實物理模型存在一定差異時,在辨識結果中可能會出現比較大的數值振蕩,甚至出現某些溫度下熱傳導系數值小于0的非物理情況。因此,有必要在辨識過程中對熱傳導系數值引入一些物理約束,例如,熱傳導系數值應大于0；對碳纖維增強樹脂基這一類防隔熱材料而言,材料的熱傳導系數值隨溫度增加而增加等。同時,在實際工程應用中,使用內置的熱電偶直接測量低熱導率材料內部溫度時,熱電偶的傳熱能力遠高于材料會導致材料測點附近的溫度場發生顯著變化,測得溫度往往低于實際溫度。例如文獻[13]在利用熱電偶測量燃燒固體溫度時,分析了熱電偶傳熱導致的溫度損失。為保證在防熱系統設計時留有較充足的設計余量,工程上通常要求溫度較高時的溫度計算預測結果能略高于溫度實測值。因此,本文在辨識材料熱傳導系數的遺傳算法基礎上進一步考慮上述這些物理約束和工程約束,建立了考慮約束的材料熱傳導系數辨識方法,并進行了算例考核分析。

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

確定材料熱傳導系數的實驗原理圖如圖1所示,假設上下邊界為絕熱邊界條件后,傳熱問題可以簡化為一維傳熱問題。假設在材料左右邊界可以進行溫度和熱流的測量,在材料內部可以進行溫度的測量,則當邊界條件和觀測條件為以下兩種情況時,可以進行材料熱傳導系數的辨識:

圖1 辨識材料熱傳導系數示意圖Fig.1 Sketch of thermal conductivity estimation

1)沒有材料內部測點溫度測量的情況:此時在左右邊界分別給定溫度或熱流測量結果作為邊界條件,同時將其中一個邊界上除邊界條件狀態量外的另一狀態量作為觀測量。

2)材料內部進行了測點溫度測量的情況:此時在左右邊界分別給定溫度或熱流測量結果作為邊界條件,將內部測點的溫度作為觀測條件。

不失一般性,本文針對第二種情況進行分析,即在材料左右邊界給定溫度邊界條件,在材料內部選取一個或多個測點進行溫度測量后來對材料的熱物性參數進行辨識。此時的傳熱方程為:

其中,ρ為材料密度,c p為材料比熱,為溫度測量值；v(t)為測量噪聲。

圖2 隨溫度變化熱傳導系數示意圖Fig.2 Sketch of temperature-dependant thermal conductivity

式(1)中的熱傳導系數值是溫度的函數,若函數形式未知,則可將溫度的可能取值范圍分段,熱傳導系數值在各段取線性函數,如圖2。設在整個溫度區間取M個節點,對應溫度為T i(i＝1,M),則k(T)可寫為:

可見,如果給定熱傳導系數在各溫度區間內的值k i(i＝1,M),則可以采用有限控制體積法[1]來進行熱傳導正問題的數值求解。然而就逆問題而言,熱傳導系數k(T)未知,需要由觀測方程式(2)中的信息來辨識,因此該逆問題可轉化為求合適的k(T)使如下目標函數達極小的優化問題:

其中,t＝[0,t f]表示溫度測量的時間段,w i為各測點溫度偏差在目標函數中占的權值。由于k針對溫度區間進行了離散,因此式(3)中的函數優化問題轉化為了針對參數k i(i＝1,N)的參數優化問題,下面介紹其優化策略。

2 辨識材料熱傳導系數的遺傳算法

對于上述優化問題,可以采用的優化方法主要分為兩類:一類是梯度類優化算法,如最速下降法和共軛梯度法；另一類是遺傳算法。由于梯度類優化算法對初值選取比較敏感,尤其是當M值較大(即溫度區間分段較密)時,優化計算易于陷入局部極值。因此,本文優化計算采用具有較強魯棒性和全局優化能力的遺傳算法[14]。其基本思路是隨機生成一個由待優化參數組合構成個體的種群,計算所有個體對應的目標函數,并轉換為適應值F的形式,如:

以使最小的目標函數對應于最大的適應值,再通過適當的選擇、復制、進化機制,使適應值較大的個體能在優化過程中被保留并不斷進化,類似于生物界中的“適者生存”,這樣經過若干代進化后最終得到的適應值最大的個體即為優化問題的最優解。對于材料熱傳導系數的辨識,J為式(3)中的形式,其具體算法如下:

(1)首先是編碼,可以采用二進制編碼或十進制編碼,將待優化參數k i(i＝1,M)分別對應成二進制數或十進制浮點數序列,組合在一起形成個體的染色體；

(2)初始化:利用隨機方法產生N個個體,形成初始種群,利用每個個體對應的k i(i＝1,M)進行熱傳導正問題計算,利用式(3)和式(4)得到各個體的J和適應值；

(3)個體選擇、復制操作:對初始種群,根據適應值決定的概率分布,采用輪盤賭選擇法來重新生成N個個體；

(4)對步驟(3)生成的N個個體進行單點交叉、均勻變異操作,產生新的種群；

(5)對新種群中的N個體進行熱傳導正問題計算,找出新種群中適應值最大的個體,判斷計算是否收斂,若收斂,則停止,否則返回步驟(3)。

下面給出兩個算例,第一個算例對熱傳導系數為常值,左端受常值熱流Q右端絕熱的一維平板傳熱問題進行分析,這一問題的邊界溫度和內點溫度都存在解析解,以這些解析解加上標準差為σ的白噪聲作為測量值,對熱傳導系數進行辨識。當L＝1,c p＝1,Q＝1,測點取在x＝0.5時分析了三種情況,一是σ＝0,即不考慮測量噪聲；二是σ取為x＝0.5處最大溫度的1%；三是σ取為x＝0.5處最大溫度的5%。表1中給出了辨識結果,從表中可以看到,即使在測量噪聲比較大的情況下,辨識結果與給定值較為符合。這一結果表明熱傳導正問題的計算方法和遺傳算法本身都是有效的。

表1 熱傳導系數的真值和辨識結果Table 1 Exact and estimated values of const thermal conductivity

第二個算例的參數設置為:材料的ρ＝1000 kg/m3,c p＝1000 J/kg·K,L＝0.01 m,T0＝600 K,材料左右邊界的溫度條件如圖3中“x＝0”,“x＝10mm”示,兩個測點分別位于x＝3 mm和x＝6 mm位置。設材料熱傳導系數隨溫度變化的函數為式(2)中的形式,T i和k i給定值列于表2的前兩行。由此可計算出兩個測點的溫度歷程,如圖3中“x＝3 mm”和“x＝6 mm”。

圖3 溫度邊界條件及測點溫度Fig.3 Temperatures of boundary condition and measurement points

下面分析逆問題,將兩個測點溫度計算值分別疊加標準差σ＝0、4 K和10 K的白噪聲來做為測量值(示于圖4),采用遺傳算法來對k i(i＝1,5)進行辨識,其中σ＝0表示不疊加白噪聲,用于驗證算法。遺傳算法采用二進制編碼,種群樣本數取為50,交叉概率為0.8,變異概率為0.05。表2和圖5中給出了此時的辨識結果,從結果中可以看到:

表2 材料熱傳導系數隨溫度變化的函數及辨識結果Table 2 Specified temperature-dependant thermal conductivity and estimation results

(1)當σ＝0時,辨識結果和給定值符合較好,相對誤差＜5‰,驗證了辨識方法的有效性。

(2)考慮測量噪聲后,辨識結果與給定值有較大偏差,辨識結果中出現了比較大的數值振蕩,對應較低溫度的k1值大于k2、k3和k4值,與真實物理情況存在一定差異。因此,下面考慮在辨識算法中加入約束來進行計算分析。

圖4 考慮測量噪聲后的測點溫度值Fig.4 Temperatures with measurement noise

圖5 熱傳導系數辨識結果Fig.5 Estimated value of thermal conductivity

3 考慮物理約束的熱傳導系數辨識算法

首先在辨識算法中考慮一項物理約束即碳纖維增強樹脂基復合材料的熱傳導系數值隨溫度增加不減小。對于上面約束條件的施加方法,考慮如下兩種方法。

3.1 辨識增量的方法

這一方法的基本思想是從某個k i(如k2)開始,不辨識熱傳導系數的絕對量,而是辨識其增量:

并且將增量值的取值范圍限定為大于或等于0的數,這樣就保證了從k2開始即滿足“熱傳導系數隨溫度增加而不減”的約束條件。

3.2 罰函數方法

由于材料的熱傳導系數值隨溫度增加不減小,這就相當于對辨識結果施加了如下約束:

為了在遺傳算法中引入此約束,需要采用罰函數的方法[15],其基本思想是:在進化過程中,若種群中某樣本個體對應的優化變量不滿足上述約束,則在其目標函數值上增加一正值作為“懲罰”,從而減小其適應值,使其在后期的交叉選擇過程中被選中的概率減小；而其它滿足約束的樣本個體則不施加這一“懲罰”,在后期的交叉選擇過程中被選中的概率較大；這樣不斷進化,可得到既有較大適應值又滿足約束的樣本作為辨識結果。

對應到遺傳算法的實際計算中,樣本目標函數的計算公式需要由式(3)變為下式:

其中ε為罰因子,ε值的選取對辨識結果有較大影響,當ε值較小時,罰函數的作用未充分顯現；而當ε值較大時,在優化計算早期會出現很多樣本都因不滿足約束而被淘汰的情況,影響優化計算效率,因此實際使用中需要根據優化結果是否滿足約束的情況來對ε值進行調整。

3.3 算例與結果分析

對應第2節中算例考慮測量噪聲的情況,分別采用辨識增量的方法和罰函數方法來進行辨識計算,辨識增量的方法中,5個優化變量為:k1,Δk2,Δk3,Δk4,Δk5；罰函數方法中優化變量仍為k i(i＝1,5),罰因子ε值取為100。表2和圖6中的“GA(With IM)”和“GA(With PM)”給出了此時的辨識結果。從結果中可以看到:

(1)在兩組不同測量噪聲水平下的辨識結果都滿足熱傳導系數值大于0和熱傳導系數值隨溫度增加不減小的約束條件,上節所建立的算法是有效的。

(2)在兩組不同測量噪聲水平下,隨著測量誤差的增加,辨識結果與給定值的偏差增大,但整體來看,辨識結果與給定值基本符合,在溫度測量標準差σ＝10 K情況(對x＝6mm測點而言,該標準差近似對應于溫升的6%)下,辨識結果與給定值的偏差小于10%。

(3)兩種考慮約束的算法得出的辨識結果基本一致,但從方法構建上來看,辨識增量的方法推導較為簡便,沒有自由參數；而罰函數方法中罰因子ε值的選取有一定的經驗成分。

圖6 考慮約束情況下的熱傳導系數辨識結果Fig.6 Estimated value of thermal conductivity with constraints taken into account

4 工程應用約束的實現

下面進一步討論工程實際應用中“溫度較高時計算溫度計算預測結果能略高于溫度實測值”這一約束的實現。以第三節中的辨識增量的方法為基礎,將待辨識參數記為＝{k1,Δk2,Δk3,Δk4,Δk5},此時的優化目標函數取為:

其中,N為時間方向上的離散點數；εil和ηil均為罰因子；T C為需要施加約束的溫度閾值,即溫度超過T C后才需要施加“計算預測結果略高于實測值”的約束；G1、G2為大于1的正數。

采用這一算法對第2節中考慮σ＝4 K測量噪聲的算例進行計算,取T C＝700 K,G1＝10000,G2＝500,此時辨識出的結果如表2的最后一行和圖7中的“GA(With IM and Penalty)”示,圖8給出了利用這一辨識結果計算出的測點溫度。從中可以看到:溫度為700 K和1000 K的熱傳導系數辨識結果顯著高于給定值而溫度為800 K的熱傳導系數辨識結果稍低于給定值,這使得利用辨識結果計算出的溫度值大于測量值,滿足工程約束要求。

圖8 測點溫度值對比Fig.8 Comparison of temperatures at measurement points

5 結 論

本文通過辨識增量和引入罰函數的方法,在辨識材料熱傳導系數的遺傳算法中實現了一項物理約束:材料熱傳導系數值隨溫度增加不減小,和一項工程應用約束:溫度較高時的仿真計算結果應高于實測值。算例計算結果表明:辨識計算結果與真值整體符合較好,考慮約束后,在測量噪聲較大的情況下,熱傳導系數辨識結果中沒有出現非物理的振蕩；測點溫度計算預測結果在高溫條件下達到了高于實測值的設計條件。辨識方法是有效的,在工程實際中有良好的應用前景,目前已在某型飛行器防熱材料的熱物性參數地面標定試驗中得到了初步的工程應用,下步將針對方法中罰因子的合理選取開展進一步的深入分析。