采用序列函數(shù)法在身管導(dǎo)熱反問(wèn)題中的應(yīng)用＊

周淵鍵 陳 剛 徐春雨

（1.解放軍陸軍軍官學(xué)院五系41隊(duì) 合肥 230031）（2.解放軍陸軍軍官學(xué)院機(jī)械工程教研室 合肥 230031）（3.解放軍陸軍軍官學(xué)院五系42隊(duì) 合肥 230031）

1 引言

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

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

文中以對(duì)火炮發(fā)射時(shí)內(nèi)膛熱交換過(guò)程的分析，基于導(dǎo)熱反問(wèn)題的研究方法，根據(jù)Beck提出的序列函數(shù)法建立由結(jié)構(gòu)內(nèi)部一個(gè)或多個(gè)溫度測(cè)點(diǎn)的測(cè)量值計(jì)算內(nèi)膛壁溫度與熱流密度的數(shù)學(xué)模型，并用FORTRAN語(yǔ)言編寫通用計(jì)算程序。

2 內(nèi)壁熱交換分析及數(shù)學(xué)模型的確定

在火炮射擊過(guò)程中，由于身管內(nèi)膛壁受到高溫高壓火藥燃?xì)饧皬椡枘Σ恋乃矔r(shí)熱作用。而熱量是通過(guò)三種基本方式傳遞，分別是導(dǎo)熱方式、對(duì)流換熱方式、熱輻射方式［2］。膛內(nèi)熱交換過(guò)程雖然復(fù)雜，但是在發(fā)射過(guò)程結(jié)束時(shí)，通過(guò)內(nèi)膛面?zhèn)魅肷砉艿臒崃侩S時(shí)間積累的總熱量是一定的。考察邊界條件的內(nèi)涵，是指導(dǎo)熱物體在其邊界面上與外部環(huán)境之間在熱交換方面的聯(lián)系或相互作用。雖然導(dǎo)熱問(wèn)題的三類邊界條件在數(shù)學(xué)表達(dá)形式上有差別，但無(wú)論導(dǎo)熱物體是以哪一種邊界條件形式與外部環(huán)境發(fā)生相互作用，熱交換結(jié)果在本質(zhì)上卻是一致的，即熱量的轉(zhuǎn)移，在形式上均表現(xiàn)為通過(guò)某一換熱面積的熱流量（熱流密度）［3］。

導(dǎo)熱反問(wèn)題是由導(dǎo)熱微分方程，初始條件，邊界條件，即導(dǎo)熱微分方程定解條件的三個(gè)組成部分再加上作為補(bǔ)充條件或輸入條件的內(nèi)部溫度測(cè)點(diǎn)信息組成。寫成數(shù)學(xué)表達(dá)式為［4］

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

定義敏感系數(shù)X為

根據(jù)敏感系數(shù)X定義和導(dǎo)熱正問(wèn)題微分方程，可得到關(guān)于X的微分方程［5］

由此可知，對(duì)線性及準(zhǔn)線性問(wèn)題［6］，X與qM無(wú)關(guān)，只與位置相關(guān)。若tM時(shí)刻節(jié)點(diǎn)k（位置xk）處熱電偶溫度測(cè)量值為YK，M（℃），則溫度可表達(dá)為

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

由導(dǎo)熱微分方程計(jì)算T（x，t，tM－1，qM－1，qM）；對(duì)照導(dǎo)熱微分方程，由式（3）計(jì)算敏感系數(shù) X（x，t，tM－1，qM－1）；利用實(shí)驗(yàn)測(cè)得的Y（k，M）由式（4）求得qM；再一次求解導(dǎo)熱正問(wèn)題離散方程，即可求出所有節(jié)點(diǎn)溫度（i＝1，2，3，4，5）。

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

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

綜合式（5）和式（6）得到熱流計(jì)算表達(dá)式

定義測(cè)點(diǎn)k在第j未來(lái)時(shí)間步的加權(quán)系數(shù)為

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

3 數(shù)值實(shí)驗(yàn)及結(jié)果分析

為驗(yàn)證序列函數(shù)法的有效性，針對(duì)前面提出的問(wèn)題，通過(guò)計(jì)算圖1所示鉻鋼炮管材料組成結(jié)構(gòu)的內(nèi)膛壁溫度、熱流密度來(lái)驗(yàn)證上述數(shù)學(xué)計(jì)算模型的正確性。

為檢驗(yàn)序列函數(shù)法的計(jì)算效果，本文給定了一組管內(nèi)膛 壁 在 0℃、45℃、90℃、135℃、180℃節(jié)點(diǎn)溫度隨時(shí)間的變化，并作為內(nèi)膛壁溫度的準(zhǔn)實(shí)驗(yàn)值，節(jié)點(diǎn)溫度變化如圖2所示。

計(jì)算檢驗(yàn)過(guò)程分為求解導(dǎo)熱正問(wèn)題和求解導(dǎo)熱反問(wèn)題兩步進(jìn)行，兩個(gè)過(guò)程都利用FORTRAN語(yǔ)言編制程序進(jìn)行計(jì)算。首先求解導(dǎo)熱正問(wèn)題，即通過(guò)求解導(dǎo)熱微分方程（1）得到管外壁溫度隨時(shí)間變化。身管外徑為20mm，管內(nèi)徑為15mm。由于計(jì)算區(qū)域的對(duì)稱性，只計(jì)算半個(gè)圓周，采用四邊形網(wǎng)格，沿圓周向?yàn)?3個(gè)節(jié)點(diǎn)，徑向?yàn)?個(gè)節(jié)點(diǎn)。身管內(nèi)膛壁上未知的節(jié)點(diǎn)溫度由已知的管內(nèi)膛壁5個(gè)節(jié)點(diǎn)溫度通過(guò)插值得到。

計(jì)算中取管外側(cè)流體溫度為T＝20℃，對(duì)流換熱系數(shù)h＝5w／（m·K），身管的熱擴(kuò)散系數(shù)a＝4.5×10－6m2／s，導(dǎo)熱系數(shù)入＝17W／（m·K）。求解導(dǎo)熱微分方程（1），得到管外壁不同位置溫度隨時(shí)間的變化，并做為管外壁溫度的準(zhǔn)實(shí)驗(yàn)值。圖3為管內(nèi)膛壁和管外壁0℃節(jié)點(diǎn)溫度變化比較。

圖1 炮管截面結(jié)構(gòu)示意圖

圖2 管內(nèi)膛壁不同位置溫度變化

圖3 管內(nèi)膛壁與管外壁0℃位置溫度變化比

圖4為0℃位置計(jì)算得到的管內(nèi)壁溫度與準(zhǔn)實(shí)驗(yàn)值比較。從圖中可以看到，計(jì)算值與準(zhǔn)實(shí)驗(yàn)值的曲線幾乎重合在一起，計(jì)算值與準(zhǔn)實(shí)驗(yàn)值相符很好。序列函數(shù)法不僅能準(zhǔn)確地捕捉到管內(nèi)壁溫度值的較小波動(dòng)，對(duì)于30s～50s之間的較大波動(dòng)也能精確地捕捉到。比較其它幾個(gè)節(jié)點(diǎn)計(jì)算值與準(zhǔn)實(shí)驗(yàn)值發(fā)現(xiàn)，序列函數(shù)法也準(zhǔn)確捕捉到了其它位置的溫度波動(dòng)，計(jì)算值與準(zhǔn)實(shí)驗(yàn)值很好地吻合。圖5為管內(nèi)壁45℃位置序列函數(shù)法計(jì)算值與準(zhǔn)實(shí)驗(yàn)值的比較，計(jì)算值與準(zhǔn)實(shí)驗(yàn)值也符合得非常好。

圖4 管內(nèi)膛壁0℃位置計(jì)算值與準(zhǔn)實(shí)驗(yàn)值比較

圖5 管內(nèi)膛壁45℃位置計(jì)算值與準(zhǔn)實(shí)驗(yàn)值比較

4 結(jié)語(yǔ)

通過(guò)建立的數(shù)學(xué)模型是針對(duì)二維非線性瞬態(tài)導(dǎo)熱反問(wèn)題的，根據(jù)Beck's序列函數(shù)思想，在計(jì)算中采用了多溫度測(cè)點(diǎn)、多未來(lái)時(shí)間步長(zhǎng)的計(jì)算方法。針對(duì)算例中內(nèi)膛壁熱流波形曲線變化情況，得到理論值和實(shí)際測(cè)量值誤差不大的理想結(jié)果。結(jié)果表明，該方法擁有較強(qiáng)的可行性和實(shí)用性。

文中建立的方法可以應(yīng)用在火炮射擊特殊狀態(tài)下內(nèi)膛壁溫度及熱流難以準(zhǔn)確測(cè)量而外部溫度測(cè)量相對(duì)容易的情況下的結(jié)構(gòu)內(nèi)膛壁熱響應(yīng)特性研究。該方法不需在內(nèi)壁面布置熱電偶和熱流密度計(jì)，也不需要考慮熱傳遞方式的影響，只需在結(jié)構(gòu)外部或外壁布置溫度熱電偶，簡(jiǎn)化了測(cè)量難度，并提高了內(nèi)膛壁熱流及溫度的精度。

［1］管屏，朱剛，馬良.生長(zhǎng)競(jìng)爭(zhēng)蟻群算法求解導(dǎo)熱反問(wèn)題［J］.上海第二工業(yè)大學(xué)學(xué)報(bào)，2010（9）：20-21.

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

［3］吳斌，夏偉，湯勇，等.射擊過(guò)程中熱影響及身管熱控制措施綜述［J］.兵工學(xué)報(bào)，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］楊冬，陳聽寬.導(dǎo)熱反問(wèn)題方法在瞬態(tài)傳熱過(guò)程中的應(yīng)用［J］.核動(dòng)力工程，1997，18（6）：553-558.

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

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

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