劉子昂，陳 強，吳邵慶，李彥斌，費慶國

(1. 東南大學機械工程學院，南京 211189；2. 東南大學江蘇省空天機械裝備工程研究中心，南京 211189)

0 引 言

再入飛行器在進出大氣層時，其結構與外側空氣相互摩擦，表面將產(chǎn)生大量熱能，形成強烈的氣動加熱效應。結構表面的熱流會沿飛行器結構向內(nèi)傳遞，進而影響飛行器結構靜、動力學特性，嚴重時可導致結構發(fā)生強度破壞。因此，準確掌握再入飛行器結構在實際服役過程中所面臨的熱流載荷，對飛行器結構設計和評估至關重要。然而，直接測量服役狀態(tài)下再入飛行器結構表面的熱流載荷存在諸多難點：一方面，再入飛行器服役于嚴酷的力/熱/振動/噪聲等復雜環(huán)境，將熱流傳感器布置于飛行器外側容易遭受多場載荷環(huán)境影響；另一方面，在再入飛行器表面的熱防護結構中集成傳感器的難度較大。預埋測溫晶粒為一種對熱防護結構影響小的測量方式，但此方法只能記錄服役過程中的最高溫度，難以準確反映溫度隨時間變化的特征。近年來所發(fā)展的載荷反演方法通過在再入飛行器內(nèi)部布置溫度傳感器并測量其內(nèi)壁面溫度，進而求解傳熱學反問題(Inverse heat conduction problem, IHCP)得到外壁面熱載荷，成為獲取再入飛行器熱流載荷的有效手段。

傳熱學反問題作為一類典型的不適定性問題，對定解數(shù)據(jù)的擾動極為敏感，國內(nèi)外學者已通過各類數(shù)值方法開展了相關研究。錢煒祺等采用順序函數(shù)法對簡單平板結構表面瞬變熱流進行反演，以當前時刻后若干時刻的結果對當前時刻熱流進行反演，反演結果依賴于用于反演的時刻數(shù)量。Lu等基于共軛梯度法，利用彎管外側布置的溫度傳感器獲得的數(shù)據(jù)對薄壁彎管結構內(nèi)部流體溫度進行反演；在彎管內(nèi)外溫差不大的情況下共軛梯度法有著較好的反演結果。相比于順序函數(shù)法，共軛梯度法中沒有人為設定的參數(shù)，并且反演精度和計算速度與順序函數(shù)法相當，因此共軛梯度法有著更好的適用性。黃少君等采用隨機優(yōu)化算法對熱源位置進行反演。隨機優(yōu)化算法相比于梯度方法在跳出局部最優(yōu)解上有著一定的優(yōu)勢，但由于其計算量大、計算時間長，難以用于復雜結構上的熱流載荷反演。Lu等基于Tikhonov正則化方法以及共軛梯度法對三維瞬態(tài)熱傳導分析問題進行反演，通過正則化方法降低測量誤差對反演結果的影響。對于較厚的結構，Khajehpour等采用域分解法，將結構沿厚度方向劃分為若干個區(qū)域，分區(qū)域?qū)崃鬟M行識別，識別結果相比于對整個結構直接識別有更高的精度?，F(xiàn)有研究大多針對簡單的厚板、圓管等結構，采用一維傳熱模型進行熱載荷反演，僅能考慮結構沿厚度方向的熱量傳遞。而對于實際結構，熱載荷不僅會沿厚度方向傳遞，同樣會由于外表面溫度的非均勻性使得熱量沿面內(nèi)方向傳遞。再入飛行器尺寸較大，并且外形復雜。采用再入飛行器結構簡化傳熱模型時，其在實際工況下的計算結果將存在較大誤差。

此外，現(xiàn)有研究大多采用較為簡單的熱載荷進行理論研究和仿真分析，取得了較好的識別效果，但這些方法在復雜熱載荷工況下的穩(wěn)定性和收斂性有待研究。Nakamura等基于順序函數(shù)法以及奇異值分解法對飛行器再入階段表面熱流進行反演；通過在計算中引入奇異值分解法降低了測量噪聲對反演結果的影響，使得反演結果不會由于過度擬合從而偏離真實值。邵元培等針對飛行器中的燒蝕情況，考慮結構表面幾何域的變化反演其表面熱流。

對于再入飛行器而言，其服役過程中結構處于上千攝氏度的高溫中。此時材料的熱導率、比熱容相比于常溫下有著很大的差異，材料熱物性參數(shù)隨溫度的變化會引起傳熱的非線性特性。薛齊文等基于T形板結構,在考慮材料溫變熱物性參數(shù)的基礎上，用共軛梯度法對熱物性參數(shù)以及邊界條件同時進行反演，基于數(shù)值仿真分析的手段探討了測點數(shù)目、測量誤差和變量初值對反演結果的影響。Bergagi等基于共軛梯度法求解邊界上的溫度，考慮了材料的溫變熱物性參數(shù)對反演結果的影響，證明了此方法可以較好地反演薄壁結構的溫度邊界載荷。Xiong等結合共軛梯度法以及順序函數(shù)法的優(yōu)點，提出了順序共軛梯度法,當測量數(shù)據(jù)受較大噪聲影響時，此方法具有高于傳統(tǒng)順序函數(shù)法的計算精度以及高于傳統(tǒng)共軛梯度法的計算速度。

本文針對再入飛行器的熱流載荷反演問題開展研究。首先，介紹基于共軛梯度法的熱流載荷反演理論基礎。進而，針對一維結構開展熱流載荷反演的數(shù)值分析，討論材料溫變熱物性參數(shù)、溫度傳感器測量誤差等因素對反演結果的影響。然后，針對典型熱防護結構和返回艙結構,結合有限元方法開展熱流載荷反演的數(shù)值分析，并討論測點位置對反演結果的影響。

1 熱流載荷反演的理論基礎

再入飛行器在服役過程中面臨復雜的熱流載荷，為了更好地擬合熱流載荷變化，將整個時程均勻離散為若干個時程節(jié)點，其中各個節(jié)點之間的熱流載荷由兩個相鄰節(jié)點的熱流載荷差值獲得。

假定結構所受熱流時程曲線函數(shù)為()，分布規(guī)律的空間函數(shù)為()，則全域上熱流載荷為：

(,)=()()

(1)

式中:為結構熱流載荷所處邊界上的位置矢量。

將傳熱過程均勻劃分為個時間節(jié)點。在結構上選取個位置作為模擬布置傳感器的位置，所受熱流載荷邊界上選取個位置。以此個位置的溫度時程曲線為已知值，以個位置的熱流時程曲線為反演值，建立優(yōu)化目標函數(shù)：

(2)

式中：mea,,為第個傳感器測量點時刻的測量溫度，cal,,為第個傳感器測量點時刻的反演結果計算值。當小于極小值之后，認為其收斂，得到反演結果。

,為第個點的熱流時程曲線時刻的熱流值，則目標函數(shù)對未知參量,求偏導數(shù)，可得到目標函數(shù)的梯度為：

=1,2,…,;=1,2,…,

(3)

熱流時程曲線參量的迭代式為：

(,)+1=(,)-(,),

=1,2,…,;=1,2,…,

(4)

式中:表示迭代步數(shù)，(,)表示第次迭代計算中獲得第個點的熱流時程曲線時刻的熱流值，表示迭代步長，,表示對第個熱流時程曲線時刻熱流值的迭代搜索方向。,由下式獲得：

=1,2,…,;=1,2,…,

(5)

(6)

式中:為共軛系數(shù)；當=0時，=0。

計算中的迭代步長為：

=

(7)

2 熱流載荷反演驗證算例研究

為驗證熱流反演方法的準確性，并探討影響反演結果的因素。針對如圖1所示的一維傳熱模型，采用有限體積法(Finite volume method, FVM)進行瞬態(tài)傳熱正問題分析計算。以結構右側的仿真溫度結果，反演其左側的熱流載荷。

2.1 有限體積法瞬態(tài)傳熱計算

對于一維情況下采用有限體積法將總長度為的區(qū)域離散為長度為Δ的個控制體，如圖1所示，其中=0處為熱流入射面，=處為絕熱面，以此處的仿真結果模擬溫度傳感器的測量值。

圖1 有限體積法中的控制體離散Fig.1 Discretization of the control volume in finite volume method

依據(jù)能量守恒定律，對于每個控制體有：

(-)Δ=Δ

(8)

式中:為單個控制體的輸入熱流，為單個控制體的輸出熱流，為控制體面積，為材料密度，為材料比熱容，Δ為時間間隔。

對于控制體,,…,-1其輸入輸出熱流為：

(9)

式中:為材料熱導率，為第個控制體的溫度。

聯(lián)立式(8)和式(9)，可得：

(10)

定義式(11)所示的傅里葉數(shù)，式(10)可簡化為：

(11)

()(+1+-1-2)=1

(12)

依據(jù)實際工況設定第1個控制體與第個控制體為第二類邊界條件。控制體1邊界為外界熱流輸入面，控制體為絕熱邊界，則邊界條件為：

(13)

式中:為環(huán)境輸入熱流。

2.2 熱流載荷反演算例

對于如圖1所示的一維結構在=0處施加如式(14)所示熱流載荷：

()=-+300,0≤≤300

(14)

基于有限體積法計算熱傳導正問題獲取控制體的仿真值。以此仿真值模擬結構內(nèi)部的實際測量溫度值。通過共軛梯度法對外側熱流載荷進行反演，計算流程圖如圖2所示。分別考慮材料熱物性參數(shù)為常數(shù)和溫變參數(shù)時的熱流反演結果，如圖3所示。

反演熱流載荷的平均誤差為：

(15)

式中：為真實熱流值，為反演熱流值。

圖2 共軛梯度法反演熱流載荷流程圖Fig.2 Flow chart of heat flux identification using conjugate gradient method

由圖3可知，材料熱物性為常數(shù)時，反演熱流平均誤差為6.84 mW·mm；材料熱物性為溫變參數(shù)時，反演熱流平均誤差為3.88 mW·mm。當考慮材料熱物性參數(shù)為常數(shù)時，其反演誤差明顯高于考慮材料熱物性為溫變參數(shù)的結果。同時，這一誤差會隨著結構整體溫度的升高而增大。這是因為材料在高溫時的熱導率和比熱容與常溫時相比有著很大差異。因此很有必要在熱流載荷反演中考慮材料熱物性隨溫度的變化特征。

圖3 熱流載荷的反演結果Fig.3 Identification results of heat flux

為模擬溫度傳感器的測量誤差，在仿真溫度值的基礎上施加方差為0.5、均值為0的隨機誤差。依據(jù)3σ準則，其誤差所在區(qū)間為[-2.12 K, 2.12 K]。由于反演結果對測量誤差的敏感性較高，且容易產(chǎn)生過擬合的情況。因此，需要對加噪聲之后的數(shù)據(jù)進行降噪處理。采用降噪處理后的數(shù)據(jù)與未降噪處理的模擬溫度測量數(shù)據(jù)對熱流載荷進行反演，反演得到的熱流載荷如圖4所示。

圖4 測量數(shù)據(jù)含噪聲時熱流反演結果Fig.4 Identified results of heat flux with noisy measurement data

表1中列出了在模擬測量結果上考慮不同測量誤差情況下熱流反演結果的誤差。由表1和圖4可知，當考慮溫度傳感器測量誤差時，其反演結果會受到較大的影響，特別是這種誤差會導致反演的結果在真值附近出現(xiàn)較大的震蕩。為了減小這種由于測量誤差所造成的反演誤差，一方面需要采用降噪方法，降低測量誤差對反演結果的影響；另一方面需要制定合適的收斂準則，使其不會因為迭代次數(shù)過多產(chǎn)生對測量數(shù)據(jù)的過度擬合。

表1 測量數(shù)據(jù)含噪聲時熱流反演誤差Table 1 Error in identified heat flux with noisy measurement data

以測量噪聲方差為0.5，均值為0的反演結果為例，分析反演過程中目標函數(shù)與熱流誤差的收斂過程，如圖5所示。其中熱流反演平均誤差以式(15)計算。

圖5 反演結果隨迭代次數(shù)的變化Fig.5 Variation of the identified results with respect to the number of iterations

當考慮溫度傳感器測量誤差時，隨著迭代次數(shù)的增加,目標函數(shù)值不斷減小，但熱流反演結果的誤差卻不會一直減小。這說明目標函數(shù)值越小并不等同于反演結果越好。由于測量誤差的存在，會導致共軛梯度法計算過程中各個時刻的搜尋方向出現(xiàn)誤差，從而影響最終的反演結果。因此，在制定共軛梯度法的收斂準則時，不僅要設置目標函數(shù)達到的最小值，也需要考慮下降速率。當目標函數(shù)下降速率較低時需要停止計算，防止過擬合。

3 實際結構算例研究

為驗證熱流反演方法在再入飛行器上的運用，并探討影響反演結果的因素。本節(jié)首先針對夾芯式熱防護結構，反演二維結構表面均勻分布的熱流載荷；進而，針對典型再入飛行器返回艙，反演三維結構上所受熱流載荷。

3.1 熱流載荷

再入飛行器在服役過程中所受熱流載荷取決于飛行器的形狀、上升軌跡和再入條件等因素。參考文獻[21]中采用了如圖6中真實值曲線所示的一種兩級入軌飛行器在概念設計中采用的熱流載荷。為使熱流反演更貼近實際情況，在實際結構熱流反演計算中采用此熱流載荷作為結構的熱流載荷。

圖6 含測量誤差時熱流反演結果Fig.6 Identified results of heat flux with measurement error

3.2 夾芯式熱防護結構熱流載荷反演

本文所研究的夾芯式熱防護結構，其從外至內(nèi)分別包括復合材料面板結構、氣凝膠結構、復合材料面板結構、膠層以及內(nèi)側金屬底板。各層的厚度依次為2 mm、28 mm、0.5 mm、2 mm、5 mm，結構有限元模型如圖7所示。

圖7 夾芯式熱防護結構有限元模型Fig.7 Finite element model of the sandwich thermal-protection structure

在實際工程應用中，其結構內(nèi)部難以布置傳感器，因此只能在其金屬底板內(nèi)側布置傳感器。其外表面受到空間上均勻分布的熱流。以底部溫度傳感器獲得的溫度值，對外部熱流載荷進行反演，獲取結構外部施加的熱流載荷。

當不考慮傳感器的測量誤差時，以底部金屬板獲得的數(shù)值仿真結果作為已知條件，反演其頂部輸入熱流，此時的熱流反演平均誤差為0.362 mW·mm。當使用不含噪聲的溫度測量數(shù)據(jù)反演時，熱流的反演結果不會出現(xiàn)過擬合的情況，因此可以在設置收斂準則時設置一個較小的值，增加迭代次數(shù)。最終的反演結果能夠很好的反映實際熱流載荷。

在此基礎上，考慮測量中由于傳感器測量誤差所造成的數(shù)據(jù)波動。分別在原模擬測量溫度結果上加上方差分別為0.01、0.1,均值為0的正態(tài)分布隨機誤差。依據(jù)3σ準則，此兩工況下誤差分別在[-0.3 K, 0.3 K]與[-0.95 K, 0.95 K]之內(nèi)。將模擬測量值進行降噪處理前后的結果分別用于反演其外表面所受熱流，結果如圖6所示。在考慮測量誤差時，反演結果在整體的熱流趨勢上反演結果較好，但對于部分熱流峰值的反演存在誤差。

由表2可見，當不考慮傳感器測量誤差時，共軛梯度法能夠有效地反演對復雜結構的熱流載荷。而考慮傳感器測量誤差時，由于復合材料夾芯結構本身良好的隔熱能力,溫度傳感器獲取到的溫度測量值較低,此時測量誤差對于反演的結果有著很大的影響。因此，為了保證熱流反演結果的精度，在布置傳感器位置時應盡量選擇溫度較高的區(qū)域，以獲取較大的靈敏度。

表2 含測量誤差時熱流反演平均誤差Table 2 Average error in heat flux identification with measurement errors

3.3 再入飛行器結構熱流反演

返回艙結構在再入大氣過程中，與空氣摩擦產(chǎn)生大量熱量。本文以日本宇宙航空研究開發(fā)機構的PARTT返回艙結構為例展開研究，其結構如圖8所示，分為返回艙主體、熱防護結構、支撐結構和殼體四部分。并在結構內(nèi)側布置模擬傳感器測點A與B，分別位于熱防護結構的內(nèi)側與外側。

圖8 返回艙結構模型Fig.8 Structural model of the reentry capsule

對于復雜三維結構，由于其內(nèi)側各部分溫差較大。在此采用不同位置的溫度測量值作為輸入，討論測點位置對熱流反演結果的影響。分別采用測點A與B的模擬測量結果作為輸入，反演其所受熱流載荷。反演結果如圖9所示。當采用測點A的溫度值作為反演輸入數(shù)據(jù)時，熱流反演平均誤差為2.22 mW·mm;當采用測點B的溫度值作為反演輸入數(shù)據(jù)時，反演平均誤差為1.30 mW·mm。

圖9 測點不同時熱流反演結果Fig.9 Identified results of heat flux with different sensor arrangement

由圖9可見，當測點位置越靠近熱流邊界位置時，反演效果也越好。這是由于越靠近熱載荷的位置，其結構溫度越高，即靈敏度相對較高。此時誤差相對于有效測量值的比例較低，測量誤差對于反演結果的影響較小。因此在布置傳感器時，需要盡量使其位置靠近熱流邊界。

在測點B處的結果上加上正態(tài)分布的隨機誤差模擬傳感器的測量誤差。分別加方差為1與2的誤差，依據(jù)3σ準則,其誤差所在區(qū)間分別為[-3 K, 3 K]與[-4.24 K, 4.24 K]。反演誤差見表3，反演結果如圖10所示。在考慮較大的測量誤差下，熱流反演結果與真實值相差較小，說明該反演方法具有良好的抗噪性。

表3 測量誤差對熱流反演結果的影響Table 3 Influence of measured error on the heat flux identification results

圖10 不同測量誤差時熱流反演結果Fig.10 Identified results of heat flux with different measurement errors

4 結 論

本文針對再入飛行器結構開展熱流載荷反演問題研究，采用共軛梯度法求解傳熱學反問題。在已知熱流載荷作用區(qū)域和準確結構傳熱模型的基礎上，通過結構內(nèi)部區(qū)域布置的測溫傳感器反演了再入飛行器結構外部熱流載荷變化過程。在此基礎上探討了材料熱物性參數(shù)的溫變特性、測點位置、傳感器誤差等因素對反演結果的影響，主要研究結論為：

1)共軛梯度法能夠克服熱流反演問題中的不適定性，可以較好地反演典型熱防護結構和返回艙結構所受的熱流載荷；

2)對于結構內(nèi)部溫度差異較大的結構，其熱載荷反演過程中呈現(xiàn)出的非線性問題，通過修正靈敏度可以有效地提高識別精度。針對本文所研究的夾芯式熱防護結構，考慮材料熱物性為溫變參數(shù)時，熱流反演的平均誤差從6.84 mW·mm降至3.88 mW·mm；

3)當反演復雜結構的熱流載荷時，結構溫度測點位置應盡量靠近主要入射熱流位置，以提高計算中的靈敏度，從而提高反演精度。對于本文所研究的返回艙結構，當溫度測點布置于熱防護外側時相比于測點布置于熱防護結構內(nèi)部，熱流反演誤差從2.22 mW·mm下降至1.30 mW·mm；

4)為了降低測量誤差對反演結果的影響，需要采用合適的降噪手段并制定合適的收斂準則，從而防止反演結果出現(xiàn)過擬合。