反演航天器在軌瞬態外熱流的導熱反問題

宋馨，張有為，劉自軍

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

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

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

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

1 數學模型

本文所研究的是一維瞬態導熱問題，一維瞬態導熱方程為

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

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

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

初始條件為

式中:T0為初始溫度.

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

共軛梯度法的目標函數為

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

式中:b為迭代次數;dn的計算式為

γ由式(8)計算:

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

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

式(12)中

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

共軛梯度法的收斂目標是使

式中:δ為很小的正數.

為增強迭代過程的收斂性，與傳統的共軛梯度法迭代過程相比進行了兩處改進:

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

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

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

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

2 計算結果

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

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

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

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

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

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

2.1 數值試驗1

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

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

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

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

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

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

2.2 數值試驗2

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

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

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

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

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

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

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

2.3 小 結

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

3 結論

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

References)

［1］張天宇，董長虹.基于自適應反演法的導彈直/氣復合制導［J］.北京航空航天大學學報，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］.北京航空航天大學學報，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］程文龍，劉娜，鐘奇，等.衛星穩態熱模型參數修正方法研究［J］.宇航學報，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］程文龍，劉娜，李志，等.衛星熱模型蒙特卡羅混合算法的修正方法應用研究［J］.科學通報，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］張鏡洋，常海萍，王立國.小衛星瞬態熱分析模型修正方法［J］.中國空間科學技術，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］張中禮，李明海，胡紹全.外壁熱流響應計算的導熱反問題方法及其驗證［J］.強度與環境，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］薛齊文，楊海天，胡國俊.共軛梯度法求解瞬態傳熱組合邊界條件多宗量反問題［J］.應用基礎與工程科學學報，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］楊海天，胡國俊.共軛梯度法求解多宗量穩態傳熱反問題［J］.應用基礎與工程科學學報，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］王登剛，劉迎曦，李守巨，等.非線性二維穩態導熱反問題的一種數值解法［J］.西安交通大學學報，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］范春利，孫豐瑞，楊立.基于紅外側溫的試件內部缺陷的識別算法研究［J］.工程熱物理學報，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).