反演航天器在軌瞬態(tài)外熱流的導(dǎo)熱反問題

宋馨，張有為，劉自軍

(北京空間飛行器總體設(shè)計部，北京 100094)

獲得航天器在軌飛行過程中的外熱流數(shù)據(jù)對于研究熱控涂層在軌退化規(guī)律、各種空間因素對熱控產(chǎn)品的影響以及航天器姿軌控發(fā)動機羽流熱效應(yīng)都有非常重要的意義，然而中國航天領(lǐng)域在這方面做的工作非常少，主要原因是由于在航天器上安裝測量外熱流的熱流計裝置進(jìn)行直接測量存在很多困難，并且要消耗寶貴的重量、空間、功率等航天器資源.因此可以利用反問題方法來克服直接測量的困難，得到滿足一定精度要求的結(jié)果.

反問題方法是根據(jù)可觀測的量值反演系統(tǒng)變化規(guī)律或參數(shù)影響規(guī)律的數(shù)學(xué)物理方法，航空航天領(lǐng)域中在導(dǎo)彈制導(dǎo)［1］、火箭推進(jìn)系統(tǒng)故障［2］等方面已經(jīng)得到良好的應(yīng)用.導(dǎo)熱反問題是反問題方法研究的一個分支，一般是利用研究對象的測量溫度，通過一定反演算法計算得到熱物性參數(shù)、邊界條件等未知量.程文龍等［3-4］提出了衛(wèi)星熱分析模型不確定參數(shù)進(jìn)行分層修正的反問題模型并取得了非常好的反演結(jié)果;楊滬寧和鐘奇［5］建立了基于蒙特卡羅法反演修正航天器熱模型參數(shù)的反問題模型;張鏡洋等［6］研究了基于蒙特卡羅法的小衛(wèi)星瞬態(tài)熱分析模型參數(shù)分層修正方法;張中禮等研究了由內(nèi)壁面溫度反演火箭外壁熱流的導(dǎo)熱反問題［7］;Lin 和 Chen 等［8-9］分別求解了平行板通道和肋片的導(dǎo)熱反問題;Huang等［10-11］求解了三維的熱傳導(dǎo)反問題;薛奇天和楊海天等求解了多宗量的熱傳導(dǎo)反問題［12-14］，得到了較滿意的結(jié)果;王登剛等［15］研究了非線性二維穩(wěn)態(tài)導(dǎo)熱反問題的數(shù)值求解方法;范春利等［16］求解熱傳導(dǎo)反問題，識別試件內(nèi)壁的缺陷.

本文采用導(dǎo)熱反問題方法，利用現(xiàn)有的航天器在軌遙測溫度數(shù)據(jù)，通過反演計算就能夠得到滿足一定精度要求的航天器在軌外熱流數(shù)據(jù).由于導(dǎo)熱反問題的不適定性及非線性，使導(dǎo)熱反問題的求解遠(yuǎn)比導(dǎo)熱正問題復(fù)雜.反演航天器在軌外熱流的導(dǎo)熱反問題需要在邊界條件中引入表征輻射熱流的4次方非線性項，大大增加了求解的難度，目前國內(nèi)外還缺乏這方面的研究工作.本文研究了利用航天器設(shè)備在軌遙測溫度值反演出航天器在軌瞬態(tài)外熱流的導(dǎo)熱反問題方法.首先推導(dǎo)了反演航天器在軌外熱流的導(dǎo)熱反問題數(shù)學(xué)模型，采用共軛梯度法求解導(dǎo)熱反問題并改進(jìn)了共軛梯度法的迭代過程以增加其抗不適定性用于求解;在此基礎(chǔ)上采用FORTRAN語言編寫通用計算程序，構(gòu)造了兩組數(shù)值試驗用于檢驗反演算法的效果和計算程序的正確性.

1 數(shù)學(xué)模型

本文所研究的是一維瞬態(tài)導(dǎo)熱問題，一維瞬態(tài)導(dǎo)熱方程為

式中:T為溫度;τ為時間;k為熱傳導(dǎo)系數(shù);ρ為密度;cp為比熱容.內(nèi)側(cè)即朝向航天器內(nèi)部的一側(cè)(x=0)為絕熱邊界條件:

外側(cè)即朝向外界空間環(huán)境的一側(cè)(x=L)為第3類邊界條件:

式中:q為研究對象的吸收外熱流值，包括吸收的太陽輻射、天體的紅外輻射以及天體反照的太陽輻射;ε為半球向發(fā)射率;σ為斯蒂芬-波爾茲曼常數(shù).

初始條件為

式中:T0為初始溫度.

本文研究的導(dǎo)熱反問題是通過測量溫度反演出q值.從式(3)中可以看到，由于邊界條件中引入了表征輻射熱流的4次方非線性項，大大增加了導(dǎo)熱反問題的不適定性.導(dǎo)熱反問題的求解方法主要包括共軛梯度法、最大熵法、正則化算法等，與其他幾種算法相比共軛梯度法在迭代過程中具有一定的抗不適定性，目前在國內(nèi)外導(dǎo)熱反問題研究領(lǐng)域仍然是最常用的方法，本文采用共軛梯度法進(jìn)行導(dǎo)熱反問題的求解.

共軛梯度法的目標(biāo)函數(shù)為

式中:Tcal，n為溫度計算值;Tmea，n為溫度在軌遙測值;n為不同時間點.吸收外熱流值q的迭代式為

式中:b為迭代次數(shù);dn的計算式為

γ由式(8)計算:

式中:N為總時間點數(shù);迭代步長β為

邊界條件式(2)和式(3)對qn求微分得到:

式(12)中

靈敏度方程的初始條件由式(4)對qn求微分得到:

共軛梯度法的收斂目標(biāo)是使

式中:δ為很小的正數(shù).

為增強迭代過程的收斂性，與傳統(tǒng)的共軛梯度法迭代過程相比進(jìn)行了兩處改進(jìn):

2)從物理概念出發(fā)，航天器在軌吸收外熱流不小于0，因此式(6)整理為

共軛梯度法的求解步驟如下:

1)求解導(dǎo)熱方程(1)得到溫度計算值Tcal，n.

3)根據(jù)溫度計算值和在軌遙測值，檢查收斂目標(biāo)式(16)是否達(dá)到;如果達(dá)到收斂目標(biāo)，則停止迭代，否則，轉(zhuǎn)入第4)步.

2 計算結(jié)果

為檢驗共軛梯度法的計算效果，本文設(shè)計了兩組數(shù)值試驗，每組數(shù)值試驗檢驗步驟如下:

1)給出一組隨時間變化的吸收外熱流值qmea.

2)以qmea為邊界條件求解導(dǎo)熱方程(1)，得到的結(jié)果作為溫度測量值Tmea.

3)Tmea作為導(dǎo)熱反問題的輸入條件，采用共軛梯度法反演吸收外熱流值qcgm.

4)比較qmea和qcgm，檢驗反演算法的效果.

研究對象為高度30mm的鋁合金圓柱體，圓柱體除一個端面朝向空間環(huán)境通過輻射交換熱量外，其余各面均為絕熱邊界，因此沿軸向可視為一維導(dǎo)熱;計算網(wǎng)格為沿軸向劃分5個節(jié)點;熱物理參數(shù)取值為密度 ρ=2 700 kg/m3，比熱容 cp=900 J/(kg·K)，熱傳導(dǎo)系數(shù)k=120W/(m·K)，輻射邊界半球向發(fā)射率ε=0.6.

2.1 數(shù)值試驗1

數(shù)值試驗1給出1組吸收外熱流值qmea如圖1所示，數(shù)值試驗1給出的吸收外熱流曲線能夠代表目前多數(shù)地球軌道航天器和深空探測航天器在軌吸收外熱流變化趨勢.通過求解導(dǎo)熱正問題(式(1))得到溫度測量值Tmea并作為導(dǎo)熱反問題的輸入條件，采用共軛梯度法反演吸收外熱流值qcgm.共軛梯度法迭代200次后J(qbn)＜0.001滿足收斂條件，查看導(dǎo)熱反問題反演出的溫度值Tcgm與測量溫度值Tmea可以看到兩者符合的很好，見圖2.

圖1 數(shù)值試驗1的吸收外熱流曲線Fig.1 Absorbed external heat flux curve of numerical experiment1

圖2 數(shù)值試驗1的溫度反演值與測量值比較Fig.2 Compare between inversion temperature and measuring data of numerical experiment1

圖3和表1分別為導(dǎo)熱反問題反演出的吸收外熱流值qcgm與真實值qmea的比較.從圖3中看到，qcgm與qmea的曲線幾乎重合在一起;由表1中的數(shù)據(jù)qcgm與qmea偏差值在－23～19W/m2之間;除了0值區(qū)域外，最大相對偏差為2.0%，共軛梯度法很好的反演出了吸收外熱流值.

圖3 數(shù)值試驗1的吸收外熱流反演值與真實值比較Fig.3 Compare between inversion absorbed external heat flux and real value of numerical experiment1

表1 數(shù)值試驗1的吸收外熱流反演值與真實值比較Tabel1 Com pare between inversion absorbed external heat flux and real value of numerical experiment 1

2.2 數(shù)值試驗2

數(shù)值試驗2給出1組吸收外熱流值qmea如圖4所示，數(shù)值試驗2的目的是為了檢驗階躍突變情況下(如航天器姿軌控發(fā)動機點火工作)的共軛梯度法反演效果.通過求解導(dǎo)熱正問題(式(1))得到溫度測量值Tmea并作為導(dǎo)熱反問題的輸入條件，采用共軛梯度法反演吸收外熱流值qcgm.共軛梯度法迭代 350次后，增大迭代次數(shù)后)在0.62～0.63之間波動不再減小.查看導(dǎo)熱反問題反演出的溫度值Tcgm與測量溫度值Tmea如圖5所示，從圖中可以看到，Tcgm與Tmea的偏差主要出現(xiàn)在曲線拐點附近，其他區(qū)域符合較好.

圖4 數(shù)值試驗2的吸收外熱流曲線Fig.4 Absorbed external heat flux curve of numerical experiment2

圖6和表2分別為導(dǎo)熱反問題反演出的吸收外熱流值qcgm與真實值qmea的比較.

圖5 數(shù)值試驗2的溫度反演值與測量值比較Fig.5 Compare between inversion temperature and measuring data of numerical experiment 2

圖6 數(shù)值試驗2的吸收外熱流反演值與真實值比較Fig.6 Compare between inversion absorbed external heat flux and real value of numerical experiment2

表2 數(shù)值試驗2的吸收外熱流反演值與真實值比較Table 2 Compare between inversion absorbed external heat flux and real value of numerical experiment 2

從圖6中看到除了在qmea出現(xiàn)階躍變化位置外其他區(qū)域共軛梯度法反演結(jié)果較好，最大相對偏差為2.9%;在階躍變化處有一個時間點的相對偏差達(dá)到了31.6%，在實用中需要對階躍位置的反演結(jié)果進(jìn)行分析處理.

2.3 小 結(jié)

通過兩組數(shù)值試驗對共軛梯度法的效果進(jìn)行了檢驗，數(shù)值試驗1中共軛梯度法很好地反演出了吸收外熱流值;數(shù)值試驗2中除了階躍變化位置其他區(qū)域能夠得到較好的反演結(jié)果，階躍位置的反演結(jié)果需要根據(jù)階躍位置以外區(qū)域的反演結(jié)果進(jìn)行分析處理.數(shù)值試驗1給出的吸收外熱流曲線能夠代表大多數(shù)地球軌道航天器以及深空探測航天器在軌吸收外熱流變化趨勢，因此本文研究的方法適用于多數(shù)航天器.

3 結(jié)論

本文研究了利用航天器設(shè)備在軌遙測溫度值反演出航天器在軌瞬態(tài)外熱流的導(dǎo)熱反問題方法.首先推導(dǎo)了導(dǎo)熱反問題數(shù)學(xué)模型，采用共軛梯度法求解導(dǎo)熱反問題并從物理概念角度改進(jìn)了共軛梯度法的迭代過程以增加其抗不適應(yīng)性.然后根據(jù)大多數(shù)地球軌道航天器以及深空探測航天器在軌吸收外熱流的特點，構(gòu)造了兩組數(shù)值試驗對共軛梯度法的反演效果進(jìn)行了檢驗.計算結(jié)果表明本文研究的方法能夠很好地反演出目前大多數(shù)地球軌道航天器以及深空探測航天器在軌吸收外熱流，對于階躍變化的吸收外熱流情況在對反演結(jié)果進(jìn)行分析處理后也能夠得到較好的反演結(jié)果.

References)

［1］張?zhí)煊睿L虹.基于自適應(yīng)反演法的導(dǎo)彈直/氣復(fù)合制導(dǎo)［J］.北京航空航天大學(xué)學(xué)報，2013，39(7):902-906.Zhang TY，Dong CH.Compound control system design based on adaptive backstepping theory［J］.Journal of Beijing University of Aeronautics and Astronautics，2013，39(7):902-906(in Chinese).

［2］楊爾輔，張振鵬，劉國球.一種推進(jìn)系統(tǒng)故障診斷反問題模型與算法［J］.北京航空航天大學(xué)學(xué)報，1999，25(6):684-687.Yang E F，Zhang Z P，Liu G Q.Model and algorithm of inverse problems on fault diagnosis for propulsion systems［J］.Journal of Beijing University of Aeronautics and Astronautics，1999，25(6):684-687(in Chinese).

［3］程文龍，劉娜，鐘奇，等.衛(wèi)星穩(wěn)態(tài)熱模型參數(shù)修正方法研究［J］.宇航學(xué)報，2010，31(1):270 －275.Cheng W L，Liu N，Zhong Q，et al.Study on parameters correction method of steady-state thermal model for spacecraft［J］.Journal of Astronautics，2010，31(1):270 －275(in Chinese).

［4］程文龍，劉娜，李志，等.衛(wèi)星熱模型蒙特卡羅混合算法的修正方法應(yīng)用研究［J］.科學(xué)通報，2010，55(20):2056-2061.Cheng W L，Liu N，Li Z，et al.Application study of a correction method for a spacecraft thermal model with a Monte-Carlo hybrid algorithm［J］.Chinese Science Bulletin，2010，55(20):2056-2061(in Chinese).

［5］楊滬寧，鐘奇.航天器熱模型蒙特卡羅法修正論述［J］.航天器工程，2009，18(3):53-58.Yang H N，Zhong Q.Monte-Carlo method for thermal model correction of spacecraft［J］.Spacecraft Engineering，2009，18(3):53-58(in Chinese).

［6］張鏡洋，常海萍，王立國.小衛(wèi)星瞬態(tài)熱分析模型修正方法［J］.中國空間科學(xué)技術(shù)，2013，33(4):24-30.Zhang JY，Chang H P，Wang L G.Correction method for transient thermal analysis model of small satellite［J］.Chinese Space Science and Technology，2013，33(4):24-30(in Chinese).

［7］張中禮，李明海，胡紹全.外壁熱流響應(yīng)計算的導(dǎo)熱反問題方法及其驗證［J］.強度與環(huán)境，2009，36(4):54-59.Zhang Z L，Li M H，Hu S Q.Nonlinear transient inverse heat conduction problems method of calculating the boundary heat response［J］.Structure ＆ Environment Engineering，2009，36(4):54-59(in Chinese).

［8］ Lin D T，Yan W M，Li H Y.Inverse problem of unsteady conjugated forced convection in parallel plate channels［J］.International Journal of Heat and Mass Transfer，2008，51(5-6):993-1002.

［9］ Chen U C，Cheng W J，Hsu JC.Two-dimensional inverse problem in estimating heat flux of pin fins［J］.International Communication of Heat and Mass Transfer，2001，28(6):793-801.

［10］ Huang C H，Wang S P.A three-dimensional inverse heat conduction problem in estimated surface heat flux by conjugate gradient method［J］.International Journal of Heat and Mass Transfer，1999，42(18):3387-3403.

［11］ Huang C H，Chen W C.A three-dimensional inverse forced convection problem in estimating surface heat flux by conjugate gradient method［J］.International Journal of Heat and Mass Transfer，2000，43(17):3171-3181.

［12］薛齊文，楊海天，胡國俊.共軛梯度法求解瞬態(tài)傳熱組合邊界條件多宗量反問題［J］.應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報，2004，12(2):113-120.Xue QW，Yang H T，Hu G J.Solving inverse heat conduction problems with multi-variables of boundary conditions in transient-state via conjugate gradient method［J］.Journal of Basic Science and Engineering，2004，12(2):113-120(in Chinese).

［13］ Xue Q W，Yan H T.A conjugate gradient method for the hyperbolic inverse heat conduction problem with multi-variables［J］.Chinese Journal of Computational Physics，2005，22(5):417-424.

［14］楊海天，胡國俊.共軛梯度法求解多宗量穩(wěn)態(tài)傳熱反問題［J］.應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報，2002，10(2):174-181.Yang H T，Hu G J.Solving inverse heat conduction problems with multi-variables in steady-state via conjugate gradient method［J］.Journal of Basic Science and Engineering，2002，10(2):174-181(in Chinese).

［15］王登剛，劉迎曦，李守巨，等.非線性二維穩(wěn)態(tài)導(dǎo)熱反問題的一種數(shù)值解法［J］.西安交通大學(xué)學(xué)報，2000，34(1):49-52.Wang D G，Liu Y X，Li S J，et a1.Two dimensional numerical solution for nonlinear inverse steady heat conduction［J］.Journal of Xi’an Jiaotong University，2000，34(1):49-52(in Chinese

［16］范春利，孫豐瑞，楊立.基于紅外側(cè)溫的試件內(nèi)部缺陷的識別算法研究［J］.工程熱物理學(xué)報，2007，28(2):304-306.Fan C L，Sun FR，Yang L.An algorithm study on identification of subsurface defect based on thermographic temperature measurement［J］.Journal of Engineering Thermophysics，2007，28(2):304-306(in Chinese).