采用序列函數法在身管導熱反問題中的應用＊

周淵鍵 陳 剛 徐春雨

（1.解放軍陸軍軍官學院五系41隊 合肥 230031）（2.解放軍陸軍軍官學院機械工程教研室 合肥 230031）（3.解放軍陸軍軍官學院五系42隊 合肥 230031）

1 引言

導熱反問題從上個世紀60年代起就受到了人們的重視。在軍事科學、工業生產研究領域的許多場合都需要知道管內流動的內壁溫度。對于管外壁溫度，有許多種實驗測量方法；而對管內壁溫度的實驗測量比較困難。利用導熱反問題方法，可通過數值求解由管外壁溫度得到管內壁溫度。導熱反問題（IHCP）是相對于導熱正問題而言的。對于導熱正問題，物體的熱物性參數和邊界條件為已知，可以由此得出物體的溫度分布。而當物體的熱物性參數或部分邊界條件必須通過已知參數來確定時，就形成了導熱反問題。

導熱反問題在各個領域中有著廣泛的應用背景，但其非適定性、非線性、計算量大等特點，使得求解比較困難。導熱反問題作為一個優化問題，已有一些求解方法，如基于梯度的優化方法（高斯牛頓法、共扼梯度法、蟻群算法）等，但都存在一定的局限性。它們不能收斂到全局最優解，并且收斂情況與初始值有關［1］。

文中以對火炮發射時內膛熱交換過程的分析，基于導熱反問題的研究方法，根據Beck提出的序列函數法建立由結構內部一個或多個溫度測點的測量值計算內膛壁溫度與熱流密度的數學模型，并用FORTRAN語言編寫通用計算程序。

2 內壁熱交換分析及數學模型的確定

在火炮射擊過程中，由于身管內膛壁受到高溫高壓火藥燃氣及彈丸摩擦的瞬時熱作用。而熱量是通過三種基本方式傳遞，分別是導熱方式、對流換熱方式、熱輻射方式［2］。膛內熱交換過程雖然復雜，但是在發射過程結束時，通過內膛面傳入身管的熱量隨時間積累的總熱量是一定的。考察邊界條件的內涵，是指導熱物體在其邊界面上與外部環境之間在熱交換方面的聯系或相互作用。雖然導熱問題的三類邊界條件在數學表達形式上有差別，但無論導熱物體是以哪一種邊界條件形式與外部環境發生相互作用，熱交換結果在本質上卻是一致的，即熱量的轉移，在形式上均表現為通過某一換熱面積的熱流量（熱流密度）［3］。

導熱反問題是由導熱微分方程，初始條件，邊界條件，即導熱微分方程定解條件的三個組成部分再加上作為補充條件或輸入條件的內部溫度測點信息組成。寫成數學表達式為［4］

其中I、B、A分別為初始算子、邊界算子和附加算子，A（x）為面積因子，對于直角坐標系A（x）＝1，對于圓柱坐標系A（x）＝r；φ（t），φ（t），Y（x，t）分別為初始條件，邊界條件和附加條件。其中只有一個邊界條件φ（t），而另一邊界條件未知，即為待求未知熱流密度。

定義敏感系數X為

根據敏感系數X定義和導熱正問題微分方程，可得到關于X的微分方程［5］

由此可知，對線性及準線性問題［6］，X與qM無關，只與位置相關。若tM時刻節點k（位置xk）處熱電偶溫度測量值為YK，M（℃），則溫度可表達為

式中，T＊（k，M）為tM－1≤t≤tM期間熱流密度為q＊時節點k處溫度值，單位℃，Xk為邊界熱流對節點k處溫度的敏感系數。因此可采用如下計算步驟：任意假設一個迭代初值q＊。

由導熱微分方程計算T（x，t，tM－1，qM－1，qM）；對照導熱微分方程，由式（3）計算敏感系數 X（x，t，tM－1，qM－1）；利用實驗測得的Y（k，M）由式（4）求得qM；再一次求解導熱正問題離散方程，即可求出所有節點溫度（i＝1，2，3，4，5）。

對于計算中取r個未來時間步的情況，唯一不同的是，假設熱流密度q（t）＝q＊在tM－1≤t≤tM內均成立，如果在固體內部不同深度布置J個熱電偶，則可采用r個未來時間步的二乘計算方法。最小二乘誤差S定義為

對于任一假定熱流值q＊（tM－1≤t≤tM），根據溫度泰勒展開式

綜合式（5）和式（6）得到熱流計算表達式

定義測點k在第j未來時間步的加權系數為

以上求解過程為迭代過程，初始時刻給定的假定熱流值q＊與真實值之間的差異將影響初始段的熱流計算誤差，而其后所有時刻的假定熱流值均取上一時刻的最終收斂熱流。

3 數值實驗及結果分析

為驗證序列函數法的有效性，針對前面提出的問題，通過計算圖1所示鉻鋼炮管材料組成結構的內膛壁溫度、熱流密度來驗證上述數學計算模型的正確性。

為檢驗序列函數法的計算效果，本文給定了一組管內膛 壁 在 0℃、45℃、90℃、135℃、180℃節點溫度隨時間的變化，并作為內膛壁溫度的準實驗值，節點溫度變化如圖2所示。

計算檢驗過程分為求解導熱正問題和求解導熱反問題兩步進行，兩個過程都利用FORTRAN語言編制程序進行計算。首先求解導熱正問題，即通過求解導熱微分方程（1）得到管外壁溫度隨時間變化。身管外徑為20mm，管內徑為15mm。由于計算區域的對稱性，只計算半個圓周，采用四邊形網格，沿圓周向為13個節點，徑向為7個節點。身管內膛壁上未知的節點溫度由已知的管內膛壁5個節點溫度通過插值得到。

計算中取管外側流體溫度為T＝20℃，對流換熱系數h＝5w／（m·K），身管的熱擴散系數a＝4.5×10－6m2／s，導熱系數入＝17W／（m·K）。求解導熱微分方程（1），得到管外壁不同位置溫度隨時間的變化，并做為管外壁溫度的準實驗值。圖3為管內膛壁和管外壁0℃節點溫度變化比較。

圖1 炮管截面結構示意圖

圖2 管內膛壁不同位置溫度變化

圖3 管內膛壁與管外壁0℃位置溫度變化比

圖4為0℃位置計算得到的管內壁溫度與準實驗值比較。從圖中可以看到，計算值與準實驗值的曲線幾乎重合在一起，計算值與準實驗值相符很好。序列函數法不僅能準確地捕捉到管內壁溫度值的較小波動，對于30s～50s之間的較大波動也能精確地捕捉到。比較其它幾個節點計算值與準實驗值發現，序列函數法也準確捕捉到了其它位置的溫度波動，計算值與準實驗值很好地吻合。圖5為管內壁45℃位置序列函數法計算值與準實驗值的比較，計算值與準實驗值也符合得非常好。

圖4 管內膛壁0℃位置計算值與準實驗值比較

圖5 管內膛壁45℃位置計算值與準實驗值比較

4 結語

通過建立的數學模型是針對二維非線性瞬態導熱反問題的，根據Beck's序列函數思想，在計算中采用了多溫度測點、多未來時間步長的計算方法。針對算例中內膛壁熱流波形曲線變化情況，得到理論值和實際測量值誤差不大的理想結果。結果表明，該方法擁有較強的可行性和實用性。

文中建立的方法可以應用在火炮射擊特殊狀態下內膛壁溫度及熱流難以準確測量而外部溫度測量相對容易的情況下的結構內膛壁熱響應特性研究。該方法不需在內壁面布置熱電偶和熱流密度計，也不需要考慮熱傳遞方式的影響，只需在結構外部或外壁布置溫度熱電偶，簡化了測量難度，并提高了內膛壁熱流及溫度的精度。

［1］管屏，朱剛，馬良.生長競爭蟻群算法求解導熱反問題［J］.上海第二工業大學學報，2010（9）：20-21.

［2］陶文銓.傳熱學［M］.西安：西安工業大學出版社，2006，12：78-103.

［3］吳斌，夏偉，湯勇，等.射擊過程中熱影響及身管熱控制措施綜述［J］.兵工學報，2003，24（4）：525-529.

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

［5］J V Beck.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］楊冬，陳聽寬.導熱反問題方法在瞬態傳熱過程中的應用［J］.核動力工程，1997，18（6）：553-558.

［8］李明海.鋼－木組合結構的火燒熱響應模擬與導熱反問題在火燒試驗中的應用［D］.重慶：重慶大學，1998：42-60.

［9］薛齊文，楊海天.應用共軛梯度法求解非線性多宗量穩態熱傳導反問題［J］.計算力學學報，2005（22）：51-54.

［10］楊海天，薛齊文.兩級敏度分析求解非線性穩態多宗量熱傳導反問題［J］.工程熱物理學報，2003，24（3）：463-465.