錢煒祺，周宇，何開鋒，邵元培

（1．空氣動力學國家重點實驗室，四川 綿陽 621000；2．中國空氣動力研究與發(fā)展中心，四川 綿陽 621000）

考慮燒蝕情況下的表面熱流辨識

錢煒祺1，2，周宇1，2，何開鋒1，2，邵元培1，2

（1．空氣動力學國家重點實驗室，四川 綿陽 621000；2．中國空氣動力研究與發(fā)展中心，四川 綿陽 621000）

針對燒蝕傳熱問題，在熱解面模型的基礎上通過伴隨方程推導建立了基于測點溫度辨識表面熱流的方法，并進行了算例考核。結果表明：熱流辨識結果與真值符合較好，辨識結果與真值之間的偏差隨測量誤差的增加而增加；燒蝕后退量測量結果的誤差對辨識結果有較為顯著的影響。然后，將該辨識方法用于鈍頭型碳酚醛材料Narmco4028試件在陶瓷加熱風洞中的燒蝕試驗結果分析，結果表明辨識出的表面熱流與加熱功率基本符合，辨識方法是有效的，在工程實際中有較好的應用前景。

燒蝕；表面熱流辨識；熱解面模型

0 引 言

飛行器在再入大氣層過程中，高速氣流流過飛行器表面，由于氣體粘性的阻滯作用，會帶來氣動加熱問題，使飛行器表面溫度明顯升高，以美國阿波羅飛船指令艙和前蘇聯(lián)聯(lián)盟號飛船返回艙為例，再入時的最高表面溫度分別達到了2950 K和1900K，這會給飛行安全帶來影響。因此，再入飛行器需進行防熱系統(tǒng)設計。燒蝕防熱是目前主要采用的防熱系統(tǒng)設計方法，即通過防熱材料的熱解、燒蝕后退來吸收并帶走氣動加熱的熱量，從而達到熱防護的目的［1-2］。

表面熱流的辨識，或稱為表面熱流的反演，是指在飛行器表面溫度較高或表面有燒蝕后退的情況下，不適宜在表面直接安裝傳感器來進行熱流和溫度測量，只能在飛行器的防熱層內(nèi)某些點上安裝溫度傳感器進行溫度測量，然后根據(jù)這些點上的溫度測量歷程來反算飛行器的表面熱流。這類問題是一類典型的熱傳導逆問題IHCP（Inverse Heat Conduction Problem），是傳熱學研究的一個重要領域。目前，針對固定幾何域的熱傳導逆問題研究較多，但由于防熱層的燒蝕后退對應著防熱層幾何外形的變化，因而此時防熱層內(nèi)的熱傳導逆問題是一個變幾何域的熱傳導逆問題，研究工作相對較少，尤其是在防熱層燒蝕后退過程中還伴有材料熱解反應、熱解氣體流動等復雜物理化學過程，使得此時的逆問題更為復雜。文獻［3］中，Kanevce利用地面燒蝕試驗的溫度測量數(shù)據(jù)對酚醛材料的熱物性參數(shù)和材料熱解反應系數(shù)進行了辨識；文獻［4-5］建立了考慮防熱材料燒蝕情況下的表面熱流反演辨識算法，但文獻［4］采用的是一層的簡化燒蝕模型，并且通過給定燒蝕表面溫度來確定外表面位置，這樣處理會導致辨識結果出現(xiàn)非物理振蕩；文獻［5］采用了兩層燒蝕熱解面模型，建立了相應的表面熱流辨識算法，但未用實測數(shù)據(jù)進行驗證。本文在前期變幾何域表面熱流辨識方法［6］研究的基礎上，借鑒文獻［5］的方法，采用兩層燒蝕熱解面模型，實現(xiàn)了表面熱流辨識算法，并用典型試驗數(shù)據(jù)對辨識算法進行了驗證。

1 燒蝕模型

防熱材料的燒蝕是一個復雜的物理化學過程。本文針對目前工程上通常采用的碳化復合材料開展研究。碳化復合材料的燒蝕模型主要有熱解層模型、熱解面模型和簡化模型，本項工作中采用的是熱解面模型［1］，如圖1示，左邊界為材料受熱邊界，右邊界為防熱材料內(nèi)部邊界，即圖中的s3。熱解面模型是將防熱材料分為兩層：碳化層和原始材料層，熱解層則簡化為一個熱解面進行處理。圖1中給定碳化邊界和外邊界燒蝕后退規(guī)律，此時的傳熱區(qū)域包含兩個運動邊界：一個是表面燒蝕后退的邊界，即圖中的s1；一個是碳化邊界，即圖中的s2，對應的數(shù)學模型為：

圖1 熱解面模型Fig.1 Pyrogeneration-plane ablation model

本文采用變幾何域的有限控制體積法（ECV：Einite Control Volume method）［7］來對式（1）進行數(shù)值求解，式中的s1（t）和s2（t）可通過兩種方法來給出，一種是給定為測量結果；另一種是給定燒蝕碳化溫度后通過迭代計算來確定其位置。同時，此時由于計算域幾何外形隨時間變化，計算網(wǎng)格相應隨時間變化，因而在離散方程中需考慮由于控制體邊界運動而引起的控制體內(nèi)能變化項。該計算方法處理變幾何域傳熱問題的有效性在文獻［6］中已得到了有解析解的一維半無限區(qū)域燒蝕表面恒溫恒速后退算例的驗證。

2 考慮燒蝕情況下的表面熱流辨識方法

通常的燒蝕問題計算是在表面熱流Q（t）已知的情況下，求解式（1）得出材料表面的后退率和內(nèi)部溫度分布。而表面熱流辨識問題則是在表面熱流Q（t）未知的情況下，根據(jù)材料端面或內(nèi)部測點的溫度歷程來反演計算出表面熱流Q（t）的值。不失一般性，下面將測量點取在右端面處，建立由測量信息來辨識表面熱流Q（t）的方法。首先，利用拉格朗日（Lagrange）乘數(shù)法將辨識問題轉化為使選取表面熱流Q（t）值使如下目標函數(shù)達極小的優(yōu)化問題：

式中的λ為伴隨變量。對此式做分部積分，代入式（1）中的相應邊界條件后，再做變分，可知為使式（2）中的目標函數(shù)達極小，伴隨變量需滿足如下伴隨方程：

而目標函數(shù)對Q的梯度為：

在此基礎上即可采用梯度類優(yōu)化算法，如最速下降法和共軛梯度法［6，8］，來對Q（t）進行辨識計算，以最速下降法為例，優(yōu)化迭代計算公式為：

式中上標“l(fā)”和“l(fā)＋1”表示迭代步數(shù)；下標“i”表示時間方向上的離散；β為步長。在優(yōu)化計算過程中，為克服不適定性的影響，優(yōu)化計算設置了如下停止準則：J≤δ；δ＝σ2tf；σ為溫度測量結果的標準差。

3 驗證算例

下面進行算法驗證：設燒蝕材料的原始材料層ρ＝1300kg／m3，比熱Cp＝1507J／kgK，kv＝1W／（mK），初始厚度為L＝20mm；碳化后的ρ＝406．3kg／m3，比熱Cp＝2000J／kg K，kc＝3W／（m K）；裂解化學反應吸收熱ΔHp＝831500 J／kg；熱解面溫度為873K；裂解氣比熱Cg＝1300J／kg K。給定表面熱流如圖2中的“Qexact”，給定燒蝕端面s1（t）的后退歷程如圖3中的“s1without noise”，采用第1節(jié)中的數(shù)值方法對溫度場進行求解，s2（t）根據(jù)熱解面溫度通過迭代計算來確定，圖4示出了計算出的熱解面位置。

圖2 表面熱流給定值與辨識結果對比Fig.2 Comparison of exact and estimated value of heat flux

圖3 燒蝕后退量S1隨時間變化歷程Fig.3 Time history of ablated surface recession

圖4 熱解面時間變化歷程Fig.4 Time history of calculated pyrogeneration-plane

以此時計算出的內(nèi)壁測點溫度歷程為實測值，采用第2節(jié)介紹的方法對表面熱流進行辨識，圖2中的虛線“Qestimated（σ＝0K）”示出了此時的辨識結果，與給定值符合較好，但在末時刻有差異，這是由于傳熱存在一定的滯后性，末段時刻的表面熱流信息在內(nèi)壁溫度測量結果中還沒有體現(xiàn)出來，所以不能有效辨識。

接下來分析測量噪聲的影響，在計算出內(nèi)壁溫度歷程上疊加標準差σ＝6K（對應測量的相對誤差約為1%）的測量白噪聲來作為測量結果（圖5中“exp．”），對表面熱流進行辨識，圖2中的虛線“Qestimated（σ＝6K）”示出了此時的辨識結果，與給定值仍符合較好；圖5中的“Eitted”給出了利用辨識結果計算出的測點溫度歷程，和測量值符合較好。進一步將疊加的測量白噪聲標準差由σ＝6K增大到σ＝12K（對應2%的測量誤差），對表面熱流進行辨識，圖2中的虛線“Qestimated（σ＝12K）”示出了此時的辨識結果，可以看到：在相對誤差2%的情況下，表面熱流的辨識結果與給定值進一步偏離，但仍較好地反映出了表面熱流的變化趨勢與強度，辨識算法具有較好的魯棒性。

圖5 測點溫度對比Fig.5 Comparison of temperature measurement

在實際應用中，燒蝕后退量的測量通常有較大的測量誤差，因而在圖3中s1（t）上疊加σ＝0.11mm（相當于燒蝕后退量最大值的10%）的白噪聲作為測量噪聲，得到的測量結果如圖3中的“s1with noise（σ＝0.00011m）”示。利用不考慮s1（t）測量噪聲時的測點溫度計算值作為實測值（不疊加測量噪聲），對表面熱流進行辨識，得到的辨識結果與給定值的比較如圖6示，從圖中可以看到，由于燒蝕后退量測量結果存在明顯波動，表面熱流的辨識結果中也出現(xiàn)了明顯波動，二者之間存在較直接的關聯(lián)關系。

4 試驗算例分析

圖6 表面熱流給定值與辨識結果對比Fig.6 Comparison of exact and estimated value of heat flux

圖7 鈍頭型碳酚醛材料地面燒蝕試驗Fig.7 Ablation experiment of blunt Carbon-phenonic material

圖8 燒蝕表面、碳化邊界計算結果和表面溫度試驗結果Fig.8 Calculated front surface and interface of ablation and experimental surface temperature

圖9 表面熱流時間歷程辨識結果和加熱熱流值對比Fig.9 Comparison of estimated heat flux and heating power

下面通過一地面試驗的實例來對辨識方法的有效性做進一步分析驗證。鈍頭型碳酚醛材料（Narmco4028）試件在Langley中心11英寸陶瓷加熱風洞中的燒蝕試驗示意圖［9］如圖7所示，材料的幾何參數(shù)在圖中標注，來流中氧的質量分數(shù)為0.13，駐點壓力5.88atm，駐點焓值2.55MJ／kg。圖8為文獻［9］中給出的燒蝕表面、碳化邊界的后退歷程計算值（圖中分別記為＂front surface＂和＂interface＂）和表面溫度測量結果（圖中記為＂front surface temperature＂），其中燒蝕表面和碳化邊界在試驗結束時刻的計算結果與實測值吻合?；谶@組結果，利用前面建立的辨識算法，可辨識出表面加熱熱流的時間變化歷程，如圖9中“Qestimated”示。可以看到，熱流辨識結果和陶瓷加熱風洞的加熱熱流3.43MW／m2是基本符合的，偏差約20%。這一結果表明：本文所建立的考慮燒蝕情況下表面熱流辨識方法是基本有效的，經(jīng)過更多算例的檢驗校核后有望在工程實際中得到更進一步的應用。

5 小 結

本文針對燒蝕傳熱問題，在熱解面模型的基礎上，推導了伴隨方程，建立了基于測點溫度辨識表面熱流的方法，并對算法的有效性進行了算例考核。結果表明：在不考慮測量噪聲時，熱流辨識結果與真值符合較好；隨著測量誤差的增大，辨識結果與真值之間的偏差增大；燒蝕后退量測量結果的誤差對辨識結果有較為顯著的影響。此后，將該辨識方法用于鈍頭型碳酚醛材料Narmco4028試件在陶瓷加熱風洞中的燒蝕試驗結果的分析，結果表明辨識出的表面熱流與加熱功率基本符合，辨識方法是有效的，在工程實際中有較好的應用前景。

［1］姜貴慶，劉連元．高速氣流傳熱與燒蝕熱防護［M］．北京：國防工業(yè)出版社，2003．

［2］王希季，林華寶，李頤黎，等．航天器進入與返回技術（下）［M］．北京：宇航出版社，1991．

［3］KANEVCE L P，KANEVCE G H，ANGELEVSKI Z Z．Comparison of two kinds of experiments for estimation of thermal properties of ablative composite［R］．In：Inverse Problems in Engineering：Theory and Practice［C］．3rd Int．Conference on Inverse Problems in Engineering．Port Ludlow，WA，USA，1999．

［4］OLIVEIRA A P D，ORLANDE H R B．Estimation of the heat flux at the surface of ablating materials by using temperature and surface position measurements［J］．Inverse Problems in Science and Engineering，2004，12（5）：563-577．

［5］HAKKAKI-EARD A，KOWSARY E．Heat flux estimation in a charring ablator［J］．Numerical Heat Transfer，Part A，2008，53：543-560．

［6］邵元培，錢煒祺，周宇，等．變幾何域傳熱的表面熱流反演方法［J］．計算力學學報，2013，30（2）：296-301．

［7］AMAR A J，BLACKWELL B E，EDWARDSJ R．One-dimensional ablation using a full Newton’s method and finite control volume procedure［J］．Journal of Thermophysics and Heat transfer，2008，22（1）：71-82．

［8］錢煒祺，周宇，何開鋒，等．非線性熱傳導逆問題的表面熱流辨識方法［J］．空氣動力學學報，2012，30（2）：145-150．

［9］SUTTON K．An experimental study of a carbon-phenolic ablation material［R］．NASA TND-5930，1970．

Heat flux estimation for heat transfer problem with ablation

QIAN Weiqi1，2，ZHOU Yu1，2，HE Kaifeng1，2，SHAO Yuanpei1，2
（1．State Key Laboratory of Aerod ynamics，Mianyang Sichuan 621000，China；2．China Aerodynamics Research and Development Center，Mianyang Sichuan 621000，China）

Ablation is one of the main means of thermal protection for hypersonic flight vehicles．In this paper，based on pyrogeneration-plane ablation model and temperature measurement at some specific locations，the surface heat flux estimation method isdeveloped with adjoint equation analysis and validated by numerical examples．It is shown that the estimated heat flux agrees with the exact value well，and thedifference between the estimated heat flux and the exact value grows as the measurement noise increases．It is also revealed that the measured noise in the ablated surface recession value may have a significant influence on the estimated heat flux．Eurthermore，the estimation method is utilized to analyze the experimental ablation results of the blunt Carbon-phenonic material Narmco4028 in ceramic-heated tunnel，it is found that the estimated surface heat flux value is close to the heating power of the arc-heater．This result further verifies that the estimation method is effective and may have a bright prospect in engineering practice．

ablation；heat flux estimation；pyrogeneration-plane ablation model

V211.3

Adoi：10.7638／kqdlxxb2014.0038

0258-1825（2014）06-0772-05

2014-05-16；

2014-11-06

國家自然科學基金項目（11372338）；空氣動力學國家重點實驗室基金資助（JBKY11030903）

錢煒祺（1973-），男，江蘇無錫，研究員，主要從事空氣動力學、飛行力學與數(shù)理反問題研究．

錢煒祺，周宇，何開鋒，等．考慮燒蝕情況下的表面熱流辨識［J］．空氣動力學學報，2014，32（6）：772-776．

10．7638／kqdlxxb-2014．0038 QIAN W Q，ZHOU Y，HE K E，et al．Heat flux estimation for heat transfer problem with ablation［J］．ACTA Aerodynamica Sinica，2014，32（6）：772-776．