MF-1飛行試驗(yàn)表面熱流辨識(shí)

周 宇, 楊慶濤, 王安齡, 袁先旭, 錢煒祺,*

(1.空氣動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室, 四川 綿陽 621000; 2.中國空氣動(dòng)力研究與發(fā)展中心, 四川 綿陽 621000)

0 引 言

熱傳導(dǎo)逆問題是相對(duì)于熱傳導(dǎo)正問題而言的。一般的正問題是給定其物理模型及其外形、物性參數(shù)等模型參數(shù)，分析模型對(duì)邊界或內(nèi)部擾動(dòng)的響應(yīng)；而熱傳導(dǎo)逆問題，是利用測(cè)得的系統(tǒng)響應(yīng)(即溫度)，來辨識(shí)邊界、內(nèi)部擾動(dòng)或模型的結(jié)構(gòu)參數(shù)(如外形、物性參數(shù)等)。辨識(shí)表面熱流是在已知傳熱模型結(jié)構(gòu)參數(shù)的情況下，由溫度測(cè)量結(jié)果辨識(shí)外表面熱流，是經(jīng)典的熱傳導(dǎo)逆問題,它可以為航空航天飛行器表面熱環(huán)境研究提供技術(shù)支撐[1]。美國20世紀(jì)60年代進(jìn)行Reentry-F飛行試驗(yàn)項(xiàng)目[2-3]利用近壁熱電偶測(cè)溫來計(jì)算表面熱流并由此判斷轉(zhuǎn)捩。X15/X17[4-5]則在金屬殼體內(nèi)安裝了熱電偶測(cè)量溫度來計(jì)算表面熱流。“FIRE”[6-8]飛行器采用了多層燒蝕防熱結(jié)構(gòu)，并在其中3層金屬薄殼中埋入多個(gè)熱電偶來測(cè)量溫度計(jì)算熱流。而國內(nèi)研究者也開展一些研究并應(yīng)用于地面試驗(yàn)，錢煒祺[9]等人開展了表面熱流辨識(shí)的原理性實(shí)驗(yàn)，張石玉[10]在高超聲速氣動(dòng)熱風(fēng)洞試驗(yàn)中利用試驗(yàn)件內(nèi)部溫度辨識(shí)了表面熱流，張聰[11]建立利用溫度計(jì)算超聲速燃燒室的內(nèi)壁面熱流方法并應(yīng)用于地面試驗(yàn)臺(tái)。

“MF-1”飛行試驗(yàn)是中國空氣動(dòng)力研究與發(fā)展中心首次開展的航天飛行試驗(yàn)，其目的是進(jìn)行技術(shù)驗(yàn)證并開展相關(guān)的空氣動(dòng)力學(xué)研究。飛行試驗(yàn)成功獲取了飛行器表面壓力和測(cè)點(diǎn)溫度數(shù)據(jù)。本文計(jì)劃利用這些測(cè)點(diǎn)溫度計(jì)算表面熱流。由溫度數(shù)據(jù)計(jì)算表面熱流關(guān)鍵在于建立與傳感器相匹配的辨識(shí)模型，考慮傳感器對(duì)結(jié)構(gòu)傳熱的影響。據(jù)此先后建立了一維、三維辨識(shí)模型，并利用氣動(dòng)熱環(huán)境數(shù)值計(jì)算所得的恢復(fù)焓考慮了傳感器表面溫度分布對(duì)壁面?zhèn)鳠釥顟B(tài)的影響，獲得了較為精確的表面熱流數(shù)據(jù)。辨識(shí)結(jié)果可以為確定邊界層流態(tài)變化(再層流化、轉(zhuǎn)捩)提供依據(jù)，并為更詳盡的氣動(dòng)熱環(huán)境數(shù)據(jù)對(duì)比提供數(shù)據(jù)支撐。

1 測(cè)溫結(jié)構(gòu)

MF-1是典型的頭錐、柱外型，具體結(jié)構(gòu)見文獻(xiàn)[12-13]。根據(jù)試驗(yàn)的目標(biāo)以及氣動(dòng)力熱、邊界層流態(tài)變化(再層流化、轉(zhuǎn)捩)計(jì)算結(jié)果，試驗(yàn)飛行器溫度測(cè)點(diǎn)總體布置方案如下：在測(cè)量錐艙和測(cè)量直艙沿流向布置了4條的測(cè)溫帶，分別位于0°～180°的4個(gè)子午線上，一共58個(gè)測(cè)點(diǎn)，主要集中在錐段以得到邊界層流態(tài)變化試驗(yàn)數(shù)據(jù)。考慮到結(jié)構(gòu)強(qiáng)度以及重心位置，各艙段的壁厚不同，其中錐艙前段不銹鋼壁厚為22 mm，測(cè)量錐艙后段不銹鋼壁厚為12 mm，測(cè)量直段的厚度為4.5 mm。由于殼體的熱傳導(dǎo)系數(shù)較低，直接用內(nèi)壁溫度來辨識(shí)飛行器外壁面熱流的精度較低。因此需要設(shè)計(jì)靈敏度較高即厚度較薄測(cè)溫部件，提高辨識(shí)精度。較薄測(cè)溫部件與旁邊較厚的殼體結(jié)構(gòu)之間必然存在相互傳熱，使得測(cè)溫部件傳熱相對(duì)復(fù)雜。而基于三維導(dǎo)熱模型進(jìn)行熱流辨識(shí)工作量大，也會(huì)引入邊界誤差。因此在設(shè)計(jì)測(cè)溫部件結(jié)構(gòu)時(shí)，主要考慮減小橫向?qū)嵝?yīng)，使測(cè)溫部件的導(dǎo)熱接近一維導(dǎo)熱，期望利用一維熱流辨識(shí)方法獲得較好的辨識(shí)效果。

為減少橫向?qū)嵝?yīng)[12]，結(jié)合結(jié)構(gòu)熱載荷設(shè)計(jì)與工藝的要求，將測(cè)溫部件結(jié)構(gòu)設(shè)計(jì)為整體厚2 mm，邊長40 mm的正方形平板，平板中間凸起直徑15 mm相對(duì)高1 mm的圓臺(tái)，而溫度測(cè)點(diǎn)在圓臺(tái)底面中心位置。對(duì)應(yīng)的在飛行器殼體外表面銑出與正方形平板匹配、深2 mm的方形槽，方形槽中心挖出直徑為30 mm的圓孔穿透金屬結(jié)構(gòu)層。裝配時(shí)將測(cè)溫薄板扣入方形槽，薄板四周與方形槽利用激光焊接成一體，熱電偶導(dǎo)線由方形槽圓孔中穿過導(dǎo)入飛行器內(nèi)部。結(jié)構(gòu)有限元模型(1/4模型)如圖1所示，圖中藍(lán)色部分為飛行器殼體，紅色部分為測(cè)溫薄壁結(jié)構(gòu)，上表面為飛行器外表面，即加熱面。

圖1 測(cè)量部件計(jì)算模型

下面利用有限元模型分析測(cè)溫結(jié)構(gòu)在沿彈道熱環(huán)境下的熱響應(yīng)。計(jì)算時(shí)假設(shè)加熱面的熱流均勻分布，加載的熱流為根據(jù)設(shè)計(jì)彈道利用數(shù)值方法計(jì)算得到的典型熱流，如圖2中的“Qexact”所示。圖3為計(jì)算所得溫度歷程。以此計(jì)算值采用一維表面熱流辨識(shí)方法辨識(shí)表面熱流。文獻(xiàn)[14-15]中有一維表面熱流辨識(shí)算法詳細(xì)的介紹，這里不再詳述。圖2中“Qestimated”為仿真辨識(shí)結(jié)果。飛行器殼體厚4.5 mm時(shí)辨識(shí)所得熱流峰值與真值相比誤差在1%左右，而殼體厚12 mm和22 mm時(shí)辨識(shí)所得熱流峰值誤差都在5%左右。由此可以看出：一維辨識(shí)結(jié)果峰值與真值偏差較小，可以為邊界層流態(tài)變化分析提供依據(jù)。但由于目前的計(jì)算模型沒有考慮飛行器殼體和測(cè)溫部件之間傳熱，辨識(shí)結(jié)果在20 s至100 s出現(xiàn)了明顯的負(fù)熱流，這與理論預(yù)期明顯不符。因此為了能夠獲得更精細(xì)的辨識(shí)結(jié)果，需要建立多維辨識(shí)模型，考慮測(cè)溫結(jié)構(gòu)和飛行器殼體之間的傳熱。

圖2 熱流辨識(shí)結(jié)果

圖3 有限元模型溫度計(jì)算結(jié)果

2 基于三維傳熱模型的表面熱流辨識(shí)原理

三維熱傳導(dǎo)控制方程：

(1)

接觸熱阻邊界條件：

初始條件：t=0時(shí)T=T0。

設(shè)觀測(cè)方程為：

t∈[0,tf];s=1,S

(2)

其中下標(biāo)“m”表示測(cè)量點(diǎn)，上標(biāo)s表示測(cè)點(diǎn)序號(hào)，S為測(cè)點(diǎn)總數(shù)，tf為末時(shí)刻值。

給出待辨識(shí)表面熱流初始值后,采用有限元方法求解溫度場(chǎng)，可得出各測(cè)點(diǎn)的溫度變化歷程。通過對(duì)比溫度計(jì)算結(jié)果和實(shí)測(cè)值的差異就可以對(duì)表面熱流進(jìn)行優(yōu)化辨識(shí)。即求解表面熱流以使得下式中的代價(jià)函數(shù)J極小：

(3)

采用拉格朗日乘數(shù)法后，代價(jià)函數(shù)可取為：

(4)

式中λ為拉格朗日乘子，通過變換可以得到當(dāng)λ滿足如下伴隨方程：

(5)

測(cè)點(diǎn)處：

s=1,S；

接觸熱阻邊界條件：

時(shí)域邊界條件：t=tf,λ(tf)= 0 。

可導(dǎo)出目標(biāo)函數(shù)對(duì)Q的導(dǎo)數(shù)為：

(6)

由伴隨方程求得梯度后可以利用牛頓類梯度法進(jìn)行迭代計(jì)算，在這里我們采用最為常用的共軛梯度法。

共軛梯度法的計(jì)算公式如下：

(7)

式中：

βn=

(8)

接觸熱阻邊界條件：

初始條件：t=0時(shí)U=0。

3 飛行試驗(yàn)的熱流辨識(shí)結(jié)果

3.1 基于一維模型的辨識(shí)結(jié)果

首先利用一維辨識(shí)模型分析了飛行試驗(yàn)溫度測(cè)量數(shù)據(jù)。圖4給出了頭錐段三個(gè)測(cè)點(diǎn)的溫升歷程和辨識(shí)所得的測(cè)點(diǎn)溫度，可以看到兩者符合很好，這說明了優(yōu)化算法的有效性。

圖4 飛行測(cè)量溫度與辨識(shí)所得溫度對(duì)比

圖5給出了利用θ=0°子午線上幾個(gè)測(cè)點(diǎn)辨識(shí)所得的熱流，F(xiàn)ly為一維模型辨識(shí)結(jié)果，Cal為工程數(shù)值計(jì)算結(jié)果，后標(biāo)TR表示邊界層再層流化、轉(zhuǎn)捩。由圖5可以看辨識(shí)結(jié)果與工程數(shù)值計(jì)算結(jié)果在大部分區(qū)域符合較好；但辨識(shí)結(jié)果在20 s～100 s出現(xiàn)了明顯的負(fù)熱流，而工程數(shù)值模擬結(jié)果為正熱流。辨識(shí)結(jié)果如此反常是因?yàn)橐痪S辨識(shí)模型沒有考慮測(cè)溫結(jié)構(gòu)和飛行器殼體之間的橫向傳熱。

圖6、圖7為不同時(shí)刻θ=0°子午線上熱流分布，Cal為數(shù)值模擬所得熱流，其中后標(biāo)CL、TL、TR表示數(shù)值模擬時(shí)氣體流動(dòng)狀態(tài)，分別對(duì)應(yīng)完全層流狀態(tài)、完全湍流狀態(tài)、層流-湍流之間轉(zhuǎn)化的狀態(tài)；Tw為數(shù)值模擬時(shí)所采用的壁溫；Fly為一維模型辨識(shí)結(jié)果。辨識(shí)結(jié)果結(jié)合數(shù)值模擬結(jié)果可以判斷實(shí)際飛行時(shí)的邊界層流態(tài)變化發(fā)生的時(shí)間與位置。特別是在圖6中的20 s和圖7中 220 s和222 s，辨識(shí)結(jié)果清晰地捕捉到了邊界層流態(tài)變化發(fā)生的區(qū)域。

圖5 一維模型辨識(shí)所得表面熱流

圖6 18 s～20 s熱流空間分布對(duì)比

圖7 218 s～222 s熱流空間分布對(duì)比

3.2 二維、三維辨識(shí)模型

為了能夠獲得更精細(xì)的辨識(shí)結(jié)果，需要建立能夠考慮測(cè)溫結(jié)構(gòu)和飛行器殼體之間傳熱的多維辨識(shí)模型。如圖1所示飛行器殼體選取的區(qū)域過大會(huì)導(dǎo)致計(jì)算量過大，但殼體選取的區(qū)域過小又會(huì)影響計(jì)算精度。因此局部模型精度還需通過仿真驗(yàn)證。

圖8為計(jì)算采用的有限元模型局部示意圖。圖中紅線所示為薄壁測(cè)溫結(jié)構(gòu)與飛行器殼體的激光焊接面;黑線所示為測(cè)溫結(jié)構(gòu)與殼體的無焊接接觸面，接觸面為內(nèi)圓外方的環(huán)形面，存在接觸熱阻。當(dāng)認(rèn)為接觸面完美接觸即接觸熱阻為0時(shí)，三維模型可退化為軸對(duì)稱模型,即二維模型。圖9為軸對(duì)稱模型，y軸為對(duì)稱軸，測(cè)點(diǎn)位于x=0處。其中二維模選取的殼體影響區(qū)域較大。二維、三維模型辨識(shí)表面熱流時(shí)均假設(shè)加熱面熱流均勻分布，其余各邊界絕熱。圖10給出了不同模型辨識(shí)所得飛行器直段測(cè)點(diǎn)表面熱流。圖中Qestimated(2D)標(biāo)示的是二維模型辨識(shí)所得熱流，Qestimated(3D，TCR=0)標(biāo)示的是接觸熱阻為0時(shí)三維模型辨識(shí)所得熱流。同時(shí)圖10還給出了二維加密網(wǎng)格情況下辨識(shí)的結(jié)果。通過對(duì)比可以看到三者符合較好，這驗(yàn)證了計(jì)算所使用的網(wǎng)格較為合適，網(wǎng)格帶來的誤差較小，而三維模型選取的傳熱影響區(qū)域也比較合適。

圖8 三維有限元模型

圖9 二維軸對(duì)稱有限元模型

圖10 二維三維模型辨識(shí)所得表面熱流

3.3 接觸熱阻對(duì)辨識(shí)結(jié)果的影響

圖8薄壁測(cè)溫結(jié)構(gòu)與飛行器殼體存在接觸熱阻，而接觸熱阻與接觸面粗糙度、壓力、溫度和材料有關(guān)。文獻(xiàn)[16]通過試驗(yàn)測(cè)量了不同壓強(qiáng)下不銹鋼接觸面的換熱系數(shù)(接觸熱阻系數(shù)的倒數(shù))。根據(jù)文獻(xiàn)的結(jié)果可以看到，接觸熱阻隨溫度和壓強(qiáng)升高而降低，在界面壓強(qiáng)為2.96 MPa時(shí)，換熱系數(shù)降至8000 W/m2K左右。而薄壁測(cè)溫結(jié)構(gòu)與飛行器殼體接觸面壓強(qiáng)大小不確定，因此在辨識(shí)中為了覆蓋接觸熱阻可能的取值范圍，分別考慮了接觸熱阻系數(shù)R為0，0.0001，0.001，0.01 m2K/W的情況，分析不同熱阻對(duì)辨識(shí)的影響，其中R為0.01對(duì)應(yīng)接觸熱阻的極大值。圖11給出了辨識(shí)結(jié)果的對(duì)比。比較幾組辨識(shí)結(jié)果可以看到，隨著接觸熱阻變大，熱流峰值不斷降低。三維的辨識(shí)結(jié)果沒有再出現(xiàn)反常的負(fù)熱流，與工程數(shù)值模擬比較接近，有效地改善了辨識(shí)結(jié)果。進(jìn)一步在圖12、圖13給出了不同模型辨識(shí)所得熱流峰值對(duì)比。

圖11 考慮接觸熱阻影響辨識(shí)所得的表面熱流

圖12 不同模型辨識(shí)所得熱流峰值對(duì)比(上升段)

圖13 不同模型辨識(shí)所得熱流峰值對(duì)比(下降段)

通過對(duì)比分析可以看到：在上升段試驗(yàn)窗口期，二維模型辨識(shí)結(jié)果峰值較高而三維模型辨識(shí)結(jié)果峰值較低。這是由于測(cè)溫結(jié)構(gòu)中間較厚的圓形區(qū)域和飛行器殼體之間的環(huán)形槽溫度較高，而二維模型的橫向傳熱較強(qiáng)，環(huán)形槽溫度不足以補(bǔ)償殼體的低溫，測(cè)點(diǎn)所在圓形區(qū)域會(huì)通過環(huán)形槽殼體傳熱，因此二維辨識(shí)結(jié)果峰值較高。三維模型中接觸熱阻阻礙了環(huán)形槽向飛行器殼體的傳熱，環(huán)形槽同時(shí)向測(cè)點(diǎn)所在圓形區(qū)域和飛行器殼體傳入熱量，此時(shí)三維模型辨識(shí)結(jié)果峰值較低。接觸熱阻為0和0.01 m2K/W分別對(duì)應(yīng)的接觸熱阻影響較小和較大的情況，其辨識(shí)結(jié)果也對(duì)應(yīng)了考慮接觸熱阻情況下的辨識(shí)上下界。

3.4 測(cè)點(diǎn)局部熱流空間分布對(duì)辨識(shí)結(jié)果的影響

上述辨識(shí)計(jì)算都假設(shè)了表面熱流均勻分布，而表面熱流與當(dāng)?shù)乇砻姹跍赜嘘P(guān)，測(cè)溫薄壁表面溫度明顯高于周圍飛行器的外殼，壁溫不同傳熱狀態(tài)亦不同。本節(jié)利用局部的表面熱流分布與表面溫度分布的關(guān)系，分析局部熱流空間分布對(duì)辨識(shí)結(jié)果的影響。

對(duì)于簡單外形利用雷諾比擬關(guān)系可以將表面熱流密度寫為：

q?St*ρ*ue(hr-hw)

其中hr為恢復(fù)焓，hw為壁面焓。hr可以根據(jù)彈道以及飛行器外形計(jì)算，計(jì)算精度也相對(duì)較高。這樣結(jié)合壁面溫度就可給出表面熱流密度的空間分布。將事先算好的氣體恢復(fù)焓代入傳熱方程，則式(1)邊界條件可以改寫為：

此時(shí)辨識(shí)表面熱流密度變?yōu)楸孀R(shí)傳熱系數(shù)C，并合理假設(shè)在計(jì)算區(qū)域C均勻分布，由此辨識(shí)所得熱流如圖14，圖中Qestimatedwithhr為考慮壁溫對(duì)熱流密度影響時(shí)辨識(shí)所得的表面熱流。通過對(duì)比可以看到,3.3節(jié)不考慮壁溫影響時(shí)辨識(shí)所得熱流稍稍高于考慮壁溫影響時(shí)辨識(shí)所得熱流，但兩者相差不大。圖15還給出了測(cè)溫薄壁結(jié)構(gòu)表面不同位置的熱流密度，x=0 mm、10 mm、20 mm分別對(duì)應(yīng)圖9中的位置。由圖9可以看到飛行器殼體表面熱流密度(接近于20 mm處)的結(jié)果明顯高于測(cè)溫薄壁結(jié)構(gòu)的表面熱流(0 mm處)的結(jié)果。這說明測(cè)溫部件感受到的熱流與飛行殼體的熱流略有不同。同時(shí)這也說明，在對(duì)比分析熱流的測(cè)量結(jié)果和理論計(jì)算結(jié)果時(shí)，需要注意其對(duì)應(yīng)的狀態(tài)。

圖14 考慮空間分布辨識(shí)所得表面熱流對(duì)比

圖15 同一測(cè)溫結(jié)構(gòu)不同位置表面熱流

4 小 結(jié)

本文根據(jù)MF-1飛行試驗(yàn)，利用表面熱流辨識(shí)方法處理了飛行試驗(yàn)的溫度測(cè)量數(shù)據(jù)，并將辨識(shí)所得的熱流與工程數(shù)值模擬預(yù)測(cè)的結(jié)果相比較。得出以下結(jié)論:

1)辨識(shí)結(jié)果與數(shù)值模擬結(jié)果比較接近，可以有效地判斷邊界層流態(tài)變化的發(fā)生。

2)合理簡化的二維、三維辨識(shí)模型可以有效模擬飛行器殼體與測(cè)溫薄壁結(jié)構(gòu)之間的橫向傳熱從而提高辨識(shí)精度。

3)由于缺乏相應(yīng)的先驗(yàn)知識(shí)，飛行器殼體與測(cè)溫薄壁之間的接觸熱阻不確定度較高。本文在辨識(shí)計(jì)算中分別考慮了接觸熱阻相對(duì)較小和較大的情況，給出了相應(yīng)的辨識(shí)結(jié)果及影響規(guī)律。

4)飛行器殼體與測(cè)溫薄壁表面溫度不一致會(huì)導(dǎo)致殼體與測(cè)溫薄壁表面熱流有明顯差別。本文利用熱環(huán)境數(shù)值計(jì)算的結(jié)果，在辨識(shí)模型中考慮了壁溫對(duì)氣動(dòng)傳熱的影響，給出了相對(duì)應(yīng)的表面熱流分布。

綜上所述，目前所采用溫度測(cè)量和分析處理方法滿足了此次飛行試驗(yàn)的要求，成功地獲取了真實(shí)飛行環(huán)境下氣動(dòng)熱數(shù)據(jù)，為邊界層流態(tài)變化等高超聲速關(guān)鍵氣動(dòng)問題研究提供了數(shù)據(jù)支撐。

致謝:中國空氣動(dòng)力研究與發(fā)展中心歐朝助理研究員、何坤副研究員、吉洪亮副研究員、張毅峰副研究員對(duì)本文工作亦有貢獻(xiàn)，在此表示感謝。