陳 剛，周淵鍵，吳 斌

(解放軍陸軍軍官學(xué)院 a.機(jī)械工程教研室;b.研究生管理大隊(duì)1 隊(duì);c.裝甲兵系，合肥 230031)

導(dǎo)熱反問題是相對(duì)導(dǎo)熱正問題而言的一個(gè)交叉學(xué)科，它涉及傳熱學(xué)、物理學(xué)、數(shù)學(xué)、計(jì)算機(jī)、實(shí)驗(yàn)技術(shù)等眾多的領(lǐng)域，在實(shí)際工程中有著廣泛的應(yīng)用前景［1-5］。

導(dǎo)熱反問題在各領(lǐng)域中有著廣泛的應(yīng)用背景，但其非適定性、非線性、計(jì)算量大等特點(diǎn)，使得求解比較困難［6-8］。導(dǎo)熱反問題作為一個(gè)優(yōu)化問題，已有一些求解方法，如基于梯度的優(yōu)化方法Levenberg-Marquardt 算法、遺傳算法、貝葉斯方法等［9-11］。

文中以對(duì)火炮發(fā)射時(shí)內(nèi)膛熱交換過程的分析，基于導(dǎo)熱反問題的研究方法，根據(jù)Levenberg-Marquardt 算法建立由結(jié)構(gòu)內(nèi)部一個(gè)或多個(gè)溫度測(cè)點(diǎn)的測(cè)量值計(jì)算內(nèi)膛壁溫度與熱流密度的數(shù)學(xué)模型。

1 導(dǎo)熱正問題

模擬三維水平火炮身管物理模型如圖1 所示，內(nèi)徑為50 mm，外徑70 mm，長(zhǎng)200 mm。本文中假設(shè)xoy 面為對(duì)稱面，水平身管上下熱工情況相同，因此只取水平身管的一側(cè)進(jìn)行討論。水平身管的內(nèi)外壁分別為第一類與第三類邊界條件，其控制方程及邊界條件為:

本文使用有限單元法求解導(dǎo)熱問題，采用八節(jié)點(diǎn)六面體等參單元離散模型，如圖1 所示。在內(nèi)壁邊界條件及身管導(dǎo)熱系數(shù)已知情況下，求解導(dǎo)熱正問題，獲得身管內(nèi)壁溫度場(chǎng)。式(1)～式(5)構(gòu)成三維火炮身管瞬態(tài)導(dǎo)熱初邊值定解問題。已知內(nèi)外壁邊界條件及初始溫度，即可計(jì)算出身管隨時(shí)間變化的溫度分布。

圖1 網(wǎng)格模型

2 導(dǎo)熱反問題

在求解導(dǎo)熱正問題時(shí)，內(nèi)外壁邊界條件已知，需要求解的未知量為外壁面的溫度;而在進(jìn)行反問題分析時(shí)，內(nèi)壁面的溫度分布為未知參量，需要根據(jù)外壁面的測(cè)量溫度，利用反演算法求得。

本文待反演的參量為身管內(nèi)壁面的溫度分布，屬于函數(shù)估計(jì)問題。Su 等借鑒有限單元法中單元形函數(shù)的思想，把反問題中的函數(shù)估計(jì)轉(zhuǎn)化為參數(shù)估計(jì)問題。本文亦用此法，把內(nèi)壁溫度場(chǎng)的反演轉(zhuǎn)化為內(nèi)壁有限結(jié)點(diǎn)溫度的反演。

在水平身管內(nèi)壁面的軸向方向，溫度值為

式中:t(x，θj)為在橫坐標(biāo)為x 的橫截面上ti(θj)為身管內(nèi)壁面圓周方向上θj角處的溫度值，為已知值，在反問題分析中為待反演參量;Ni(x)為身管內(nèi)壁軸向上的形函數(shù)，本文算例中采用線性插值函數(shù)，類似地，在水平身管內(nèi)壁面圓周方向，溫度值由式(7)確定

本文研制的是一種無機(jī)泡沫吸波材料，是利用碎玻璃或火山灰等為主要原材料，在其中添加電磁損耗物質(zhì)，熔融發(fā)泡而成的一種隱身材料，如圖7所示。

式中:t(xj，θ)為橫坐標(biāo)為xj的橫截面是θ 角處的溫度值，θ為變量;ti(xj)為管內(nèi)壁橫坐標(biāo)為xj的橫截面上的溫度值，為已知值，在反問題分析中為待反演參量;Ni(θ)為身管內(nèi)壁圓周上的形函數(shù)，本文算例采用線性插值函數(shù)。

利用式(6)和式(7)，反問題由函數(shù)估計(jì)問題轉(zhuǎn)化為參數(shù)估計(jì)問題。為了構(gòu)建導(dǎo)熱反問題的目標(biāo)函數(shù)，先把未知參量表示為向量形式

式中:pi表示內(nèi)壁面選定點(diǎn)處的溫度;N 為選定待反演未知參量的總數(shù)。

構(gòu)建反問題目標(biāo)函數(shù)

式中:T(Pm)為未知向量P 時(shí)計(jì)算導(dǎo)熱正問題得到的在測(cè)量點(diǎn)處的溫度值;Zm為測(cè)量溫度值;m = 1，2，…，M，M 為測(cè)點(diǎn)總數(shù)。

式(9)可寫為

式中:Y 為測(cè)量點(diǎn)計(jì)算溫度值與測(cè)量溫度值的向量差，且Ym=Tm-Zm。由此，導(dǎo)熱反問題即可視為一個(gè)優(yōu)化問題，即利用一定的反演算法，構(gòu)造出內(nèi)壁邊界條件，以使目標(biāo)函數(shù)J(P)達(dá)到最小。

對(duì)Y(P)進(jìn)行Tayor 展開，并省去二階以上高階項(xiàng)，代人式(11)并整理成迭代式得

式中:R 為Jacobian 矩陣。

引入修正系數(shù)a 調(diào)節(jié)收斂速度，得Levenberg-Marquar 算法迭代式

式中，I 為(n+1)階單位矩陣。

通過不斷地計(jì)算正問題，由式(14)求得新的反演參量，直到未知參量值P 滿足迭代終止條件式(16)

式中:ε 為較小的正數(shù)，如1 ×10-6。

3 試驗(yàn)驗(yàn)證

針對(duì)前面提出的問題，以直角坐標(biāo)系為例，通過計(jì)算鋼-耐火纖維-鋼三層材料組成結(jié)構(gòu)的外表面溫度、熱流密度來驗(yàn)證上述數(shù)學(xué)計(jì)算模型的正確性。試件結(jié)構(gòu)及尺寸如圖2所示，試件總厚度為80 mm，外殼是厚度為15 mm 的不銹鋼元件，內(nèi)殼不銹鋼元件厚度為15 mm，中間層為隔熱緩沖層，厚度為50 mm。

假設(shè)在距外壁面5 mm，10 mm，25 mm 處各布置有一對(duì)熱電偶測(cè)溫點(diǎn)，分別命名為1#，2#，3#測(cè)溫點(diǎn)。為了計(jì)算方便，其初始條件假設(shè)為室溫25 ℃，已知邊界條件假設(shè)為絕熱邊界。通過三角波形熱流變化模擬和驗(yàn)證導(dǎo)熱反問題在外壁熱流在該變化情況下的反演計(jì)算，并考察測(cè)溫點(diǎn)位置、數(shù)目及時(shí)間步長(zhǎng)對(duì)計(jì)算結(jié)果誤差的影響。

首先假設(shè)真實(shí)熱流變化變化曲線如圖2 所示，以圖2 中三角波形熱流作為初始條件，通過計(jì)算導(dǎo)熱正問題方程可得到結(jié)構(gòu)內(nèi)部溫度場(chǎng)，其中外壁面溫度變化曲線如圖3 所示;然后取3 個(gè)測(cè)溫點(diǎn)所在位置的溫度變化歷程作為反問題計(jì)算輸入條件，可以得到反問題計(jì)算熱流及外壁面溫度變化曲線如圖4 和圖5。

圖2 組合結(jié)構(gòu)和熱電偶分布

圖3 實(shí)際熱流變化曲線

圖4 正問題計(jì)算外壁面溫度曲線

圖5 反演熱流變化曲線

從圖2 ～5 可以看出，利用正問題計(jì)算得到的內(nèi)部測(cè)溫點(diǎn)溫度值，作為導(dǎo)熱反問題輸入條件的導(dǎo)熱反問題方法，對(duì)于三角波形變化熱流(也即緩變熱流)情況下的計(jì)算，無論是外壁熱流值還是外壁面溫度值都符合較好。

圖6 反問題計(jì)算外壁面溫度曲線

4 結(jié)束語(yǔ)

本文所討論的導(dǎo)熱反問題方法，適用于對(duì)航天器、火炮身管、核電及化工等大型工程設(shè)備進(jìn)行瞬態(tài)傳熱研究時(shí)計(jì)算介質(zhì)對(duì)容器內(nèi)壁的放熱系數(shù)，該方法不需要布置熱電偶和熱流密度計(jì)，而僅需要在容器外壁或內(nèi)壁布置熱電偶，大大簡(jiǎn)化了試驗(yàn)難度。通過計(jì)算模擬鋼-耐火纖維-鋼三層材料組成結(jié)構(gòu)的外表面溫度、熱流密度來驗(yàn)證上述數(shù)學(xué)計(jì)算模型的正確性。

［1］薛齊文，楊海天.應(yīng)用共軛梯度法求解非線性多宗量穩(wěn)態(tài)熱傳導(dǎo)反問題［J］.計(jì)算力學(xué)學(xué)報(bào)，2005(22):51-54.

［2］陶文銓.傳熱學(xué)［M］.西安:西安工業(yè)大學(xué)出版社，2006:78-103.

［3］吳斌，夏偉，湯勇，等.射擊過程中熱影響及身管熱控制措施綜述［J］.兵工學(xué)報(bào)，2003，24(4):525-529.

［4］Beck J V，Blackwell B，St C R.clair inverse heat conduction-Ill posed problems［M］. New York: Wiley - Interscience，1985(1):218-243.

［5］Beck J V.Methodology for comparison of inverse heat conduction methods［J］. J. of heat transfer，1988( 110): 30-37.

［6］Lawton B. Thermal-chemical erosion in gun barrels［J］.Wear，2001(251):827-838.

［7］楊冬，陳聽寬.導(dǎo)熱反問題方法在瞬態(tài)傳熱過程中的應(yīng)用［J］.核動(dòng)力工程，1997，18(6):553-558.

［8］李明海.鋼-木組合結(jié)構(gòu)的火燒熱響應(yīng)模擬與導(dǎo)熱反問題在火燒試驗(yàn)中的應(yīng)用［D］.重慶:重慶大學(xué)，1998:42-60.

［9］薛齊文，楊海天.應(yīng)用共軛梯度法求解非線性多宗量穩(wěn)態(tài)熱傳導(dǎo)反問題［J］.計(jì)算力學(xué)學(xué)報(bào)，2005(22):51-54.

［10］楊海天，薛齊文.兩級(jí)敏度分析求解非線性穩(wěn)態(tài)多宗量熱傳導(dǎo)反問題［J］. 工程熱物理學(xué)報(bào)，2003，24(3):463-465.

［11］王靖君，赫信鵬. 火炮概論［M］. 北京: 兵器工業(yè)出版社，1997.

［12］趙金輝，何忠波，傅建平，等.火炮發(fā)射過程中身管溫度場(chǎng)及彎曲度的有限元計(jì)算［J］. 火力與指揮控制，2011(5):106-109.