宋馨，張有為，劉自軍

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

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

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

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

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

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

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

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

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

初始條件為

式中:T0為初始溫度.

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

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

式中:Tcal，n為溫度計(jì)算值;Tmea，n為溫度在軌遙測(cè)值;n為不同時(shí)間點(diǎn).吸收外熱流值q的迭代式為

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

γ由式(8)計(jì)算:

式中:N為總時(shí)間點(diǎn)數(shù);迭代步長(zhǎng)β為

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

式(12)中

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

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

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

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

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

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

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

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

2 計(jì)算結(jié)果

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

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

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

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

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

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

2.1 數(shù)值試驗(yàn)1

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

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

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

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

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

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

2.2 數(shù)值試驗(yàn)2

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

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

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

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

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

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

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

2.3 小 結(jié)

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

3 結(jié)論

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

References)

［1］張?zhí)煊?，董長(zhǎng)虹.基于自適應(yīng)反演法的導(dǎo)彈直/氣復(fù)合制導(dǎo)［J］.北京航空航天大學(xué)學(xué)報(bào)，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］楊爾輔，張振鵬，劉國(guó)球.一種推進(jìn)系統(tǒng)故障診斷反問(wèn)題模型與算法［J］.北京航空航天大學(xué)學(xué)報(bào)，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é)報(bào)，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é)通報(bào)，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］張鏡洋，常海萍，王立國(guó).小衛(wèi)星瞬態(tài)熱分析模型修正方法［J］.中國(guó)空間科學(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)計(jì)算的導(dǎo)熱反問(wèn)題方法及其驗(yàn)證［J］.強(qiáng)度與環(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］薛齊文，楊海天，胡國(guó)俊.共軛梯度法求解瞬態(tài)傳熱組合邊界條件多宗量反問(wèn)題［J］.應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào)，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］楊海天，胡國(guó)俊.共軛梯度法求解多宗量穩(wěn)態(tài)傳熱反問(wèn)題［J］.應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào)，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)熱反問(wèn)題的一種數(shù)值解法［J］.西安交通大學(xué)學(xué)報(bào)，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)部缺陷的識(shí)別算法研究［J］.工程熱物理學(xué)報(bào)，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).