考慮約束的材料熱傳導(dǎo)系數(shù)辨識方法

周 宇錢煒祺邵元培顏維旭

(1.空氣動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室,綿陽 621000；2.中國空氣動力研究與發(fā)展中心 計(jì)算空氣動力研究所,綿陽 621000；3.中國運(yùn)載火箭技術(shù)研究院,北京 100076)

0 引 言

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

1 辨識材料熱傳導(dǎo)系數(shù)的數(shù)學(xué)模型

確定材料熱傳導(dǎo)系數(shù)的實(shí)驗(yàn)原理圖如圖1所示,假設(shè)上下邊界為絕熱邊界條件后,傳熱問題可以簡化為一維傳熱問題。假設(shè)在材料左右邊界可以進(jìn)行溫度和熱流的測量,在材料內(nèi)部可以進(jìn)行溫度的測量,則當(dāng)邊界條件和觀測條件為以下兩種情況時(shí),可以進(jìn)行材料熱傳導(dǎo)系數(shù)的辨識:

圖1 辨識材料熱傳導(dǎo)系數(shù)示意圖Fig.1 Sketch of thermal conductivity estimation

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

2)材料內(nèi)部進(jìn)行了測點(diǎn)溫度測量的情況:此時(shí)在左右邊界分別給定溫度或熱流測量結(jié)果作為邊界條件,將內(nèi)部測點(diǎn)的溫度作為觀測條件。

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

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

圖2 隨溫度變化熱傳導(dǎo)系數(shù)示意圖Fig.2 Sketch of temperature-dependant thermal conductivity

式(1)中的熱傳導(dǎo)系數(shù)值是溫度的函數(shù),若函數(shù)形式未知,則可將溫度的可能取值范圍分段,熱傳導(dǎo)系數(shù)值在各段取線性函數(shù),如圖2。設(shè)在整個溫度區(qū)間取M個節(jié)點(diǎn),對應(yīng)溫度為T i(i＝1,M),則k(T)可寫為:

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

其中,t＝[0,t f]表示溫度測量的時(shí)間段,w i為各測點(diǎn)溫度偏差在目標(biāo)函數(shù)中占的權(quán)值。由于k針對溫度區(qū)間進(jìn)行了離散,因此式(3)中的函數(shù)優(yōu)化問題轉(zhuǎn)化為了針對參數(shù)k i(i＝1,N)的參數(shù)優(yōu)化問題,下面介紹其優(yōu)化策略。

2 辨識材料熱傳導(dǎo)系數(shù)的遺傳算法

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

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

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

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

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

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

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

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

表1 熱傳導(dǎo)系數(shù)的真值和辨識結(jié)果Table 1 Exact and estimated values of const thermal conductivity

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

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

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

表2 材料熱傳導(dǎo)系數(shù)隨溫度變化的函數(shù)及辨識結(jié)果Table 2 Specified temperature-dependant thermal conductivity and estimation results

(1)當(dāng)σ＝0時(shí),辨識結(jié)果和給定值符合較好,相對誤差＜5‰,驗(yàn)證了辨識方法的有效性。

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

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

圖5 熱傳導(dǎo)系數(shù)辨識結(jié)果Fig.5 Estimated value of thermal conductivity

3 考慮物理約束的熱傳導(dǎo)系數(shù)辨識算法

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

3.1 辨識增量的方法

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

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

3.2 罰函數(shù)方法

由于材料的熱傳導(dǎo)系數(shù)值隨溫度增加不減小,這就相當(dāng)于對辨識結(jié)果施加了如下約束:

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

對應(yīng)到遺傳算法的實(shí)際計(jì)算中,樣本目標(biāo)函數(shù)的計(jì)算公式需要由式(3)變?yōu)橄率?

其中ε為罰因子,ε值的選取對辨識結(jié)果有較大影響,當(dāng)ε值較小時(shí),罰函數(shù)的作用未充分顯現(xiàn)；而當(dāng)ε值較大時(shí),在優(yōu)化計(jì)算早期會出現(xiàn)很多樣本都因不滿足約束而被淘汰的情況,影響優(yōu)化計(jì)算效率,因此實(shí)際使用中需要根據(jù)優(yōu)化結(jié)果是否滿足約束的情況來對ε值進(jìn)行調(diào)整。

3.3 算例與結(jié)果分析

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

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

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

(3)兩種考慮約束的算法得出的辨識結(jié)果基本一致,但從方法構(gòu)建上來看,辨識增量的方法推導(dǎo)較為簡便,沒有自由參數(shù)；而罰函數(shù)方法中罰因子ε值的選取有一定的經(jīng)驗(yàn)成分。

圖6 考慮約束情況下的熱傳導(dǎo)系數(shù)辨識結(jié)果Fig.6 Estimated value of thermal conductivity with constraints taken into account

4 工程應(yīng)用約束的實(shí)現(xiàn)

下面進(jìn)一步討論工程實(shí)際應(yīng)用中“溫度較高時(shí)計(jì)算溫度計(jì)算預(yù)測結(jié)果能略高于溫度實(shí)測值”這一約束的實(shí)現(xiàn)。以第三節(jié)中的辨識增量的方法為基礎(chǔ),將待辨識參數(shù)記為＝{k1,Δk2,Δk3,Δk4,Δk5},此時(shí)的優(yōu)化目標(biāo)函數(shù)取為:

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

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

圖8 測點(diǎn)溫度值對比Fig.8 Comparison of temperatures at measurement points

5 結(jié) 論

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