李慶華，馮子超，陳莘莘，孔祥祿

（1.華東交通大學(xué)土木建筑學(xué)院，江西 南昌330013； 2.中國(guó)鐵路哈爾濱局集團(tuán)有限公司，黑龍江 哈爾濱150001）

材料參數(shù)隨溫度變化引起的非線性熱傳導(dǎo)問(wèn)題廣泛存在于航空航天、生物醫(yī)學(xué)工程和冶金工程等諸多實(shí)際工程領(lǐng)域[1-3]。 由于非線性熱傳導(dǎo)問(wèn)題的復(fù)雜性，很難得到解析解，因而研究其數(shù)值解十分必要。 Thakur 等[4]求解了二維非線性瞬態(tài)熱傳導(dǎo)問(wèn)題。 王峰等[5]采用基于Kriging 插值的MLPG 法分析了二維非線性穩(wěn)態(tài)和瞬態(tài)熱傳導(dǎo)問(wèn)題。 Yang 等[6]采用徑向積分邊界元法分析導(dǎo)熱系數(shù)隨溫度變化的瞬態(tài)非線性熱傳導(dǎo)問(wèn)題。Cui 等[7]采用單元微分法分析多維瞬態(tài)非線性熱傳導(dǎo)問(wèn)題。 朱強(qiáng)華等[8]提出了一種基于特征正交分解（POD）和有限元法的瞬態(tài)非線性熱傳導(dǎo)問(wèn)題的模型降階快速分析方法。Mierzwiczak 等[9]采用Kirchhoff 變換對(duì)導(dǎo)熱系數(shù)隨溫度變化的穩(wěn)態(tài)非線性熱傳導(dǎo)問(wèn)題進(jìn)行線性化，并采用奇異邊界法進(jìn)行分析。

作為一種半解析的數(shù)值方法，比例邊界有限元法[10]只需對(duì)計(jì)算域的環(huán)向邊界離散為單元，并且在未離散的徑向可以解析求解。 相對(duì)于邊界元法，這種方法不需要基本解，也不涉及奇異積分的數(shù)值困難。 目前，比例邊界有限元法已被用于求解斷裂力學(xué)問(wèn)題[11-12]、彈性動(dòng)力學(xué)問(wèn)題[13-14]、波導(dǎo)本征問(wèn)題[15]、土-結(jié)構(gòu)動(dòng)力相互作用分析[16]和熱傳導(dǎo)問(wèn)題[17]。 但目前尚未見(jiàn)到采用比例邊界有限元法進(jìn)行非線性熱傳導(dǎo)問(wèn)題的工作。

鑒于比例邊界有限元法的顯著優(yōu)勢(shì)，本文將其與Kirchhoff 變換[18-19]相結(jié)合，提出了求解穩(wěn)態(tài)非線性熱傳導(dǎo)問(wèn)題的一種新方法。 利用Kirchhoff 變換，將非線性的偏微分控制方程轉(zhuǎn)化為線性方程，然后采用比例邊界有限元法求解，并借助Kirchhoff 反變換求得溫度場(chǎng)。 最后，通過(guò)兩個(gè)典型算例驗(yàn)證了本文所提方法的有效性。

1 穩(wěn)態(tài)非線性熱傳導(dǎo)問(wèn)題

考慮一個(gè)二維穩(wěn)態(tài)熱傳導(dǎo)問(wèn)題， 其計(jì)算域?yàn)棣福?邊界為Γ。 當(dāng)熱傳導(dǎo)問(wèn)題的導(dǎo)熱系數(shù)隨溫度變化，且無(wú)內(nèi)熱源，則相應(yīng)的控制方程和邊界條件可寫(xiě)為

式中：k（T）為隨溫度變化的熱傳導(dǎo)系數(shù)；T 為溫度；n 為邊界外法向向量；Γ1為Dirichlet 邊界；Tˉ為其上給定的溫度；Γ2為Neumann 邊界；q 為其上給定的熱流密度。

為了消除式（1）的非線性，本文采用Kirchhoff變換[18－19]。 引入下式定義與溫度相關(guān)的新變量

2 比例邊界有限元法

考慮如圖1 所示的比例邊界坐標(biāo)系，相似中心為點(diǎn)O，要求從該點(diǎn)可以看到邊界上任意一點(diǎn)。 對(duì)于有限域問(wèn)題，徑向坐標(biāo)ξ 的變化范圍從0 到1，即0≤ξ≤1。 當(dāng)相似中心點(diǎn)與直角坐標(biāo)系原點(diǎn)一致時(shí)，計(jì)算域內(nèi)任意點(diǎn)的笛卡爾坐標(biāo)（x^，y^）可寫(xiě)為

圖1 比例邊界坐標(biāo)系統(tǒng)Fig.1 The scaled boundary coordinate system

式中：（x（s），y（s））分別為邊界上任意點(diǎn)的笛卡爾坐標(biāo)。

在比例邊界坐標(biāo)系下，梯度算子可表示為[14－15]

式中：N（s）為邊界結(jié)點(diǎn)形函數(shù)。

將式（12）和式（15）代入式（14），然后對(duì)含δθ（ξ），ξ項(xiàng)做分部積分，并考慮δθ（ξ）T任意性整理可得

式（16）為二階常微分方程組，采用矩陣函數(shù)[20]法可將其解表示為

式中：Sn為與負(fù)特征值對(duì)應(yīng)的Schur 標(biāo)準(zhǔn)型；Ψθ為Schur 分解中與θ 變量相關(guān)的模態(tài)；c 為可由邊界條件確定的積分常數(shù)。

將式（24）代入式（15）可求得域內(nèi)各點(diǎn)的變量θ（ξ，s）值后，由Kirchhoff 反變換式（6）可進(jìn)一步求得相應(yīng)各點(diǎn)的溫度T（ξ，s）。

3 數(shù)值算例

為了驗(yàn)證所提方法的有效性，本文對(duì)一些典型的穩(wěn)態(tài)非線性熱傳導(dǎo)問(wèn)題進(jìn)行分析求解，并與相應(yīng)的解析解進(jìn)行了對(duì)比。

3.1 算例1

考慮一個(gè)邊長(zhǎng)為1 的方板的穩(wěn)態(tài)非線性熱傳導(dǎo)問(wèn)題，如圖2 所示。 板的上端和下端為絕熱邊界，板左端的溫度保持為T(mén)l=300 K，右端的溫度保持為T(mén)r=400 K。

圖2 算例1 和算例2 的模型Fig.2 The model of example 1 and example 2

導(dǎo)熱系數(shù)與溫度的關(guān)系為

計(jì)算時(shí)取k1=-2 和k2=0.01，將相似中心取在方板中心，并將邊界離散成16 個(gè)二次單元。 采用本文方法，可以計(jì)算得到m 取不同值時(shí)y=0.5 處的溫度分布，如圖3 所示。 為了方便比較，圖3 還給出了相應(yīng)的解析解。 由圖3 可知，本文方法計(jì)算精度很高，與解析解之間的最大相對(duì)誤差只有2.02×10-10%。

圖3 算例1 中y=0.5 處的溫度分布Fig.3 Temperature distributions at y=0.5 for example 1

3.2 算例2

算例2 同算例1 一樣，只是導(dǎo)熱系數(shù)與溫度的關(guān)系假設(shè)為

計(jì)算時(shí)取k1＝-2，將相似中心取在方板中心，并將邊界離散成16 個(gè)二次單元。 圖4 給出了k2取不同值時(shí)y=0.5 處溫度分布的數(shù)值解和解析解的比較。 顯然，本文數(shù)值解與解析解吻合非常好（最大相對(duì)誤差為2.32×10-10%）， 這進(jìn)一步驗(yàn)證了本文所提方法的有效性。

圖4 算例2 中y=0.5 處的溫度分布Fig.4 Temperature distributions at y=0.5 for example 2

4 結(jié)論

本文首次采用比例邊界有限元法，建立了穩(wěn)態(tài)非線性熱傳導(dǎo)問(wèn)題分析的整套求解算法。 由本文分析和算例求解結(jié)果可以看出。

1） 不同于迭代解法，對(duì)于導(dǎo)熱系數(shù)隨溫度變化的穩(wěn)態(tài)非線性熱傳導(dǎo)問(wèn)題應(yīng)用Kirchhoff 變換可轉(zhuǎn)化為線性熱傳導(dǎo)問(wèn)題進(jìn)行求解，進(jìn)而可以避免迭代計(jì)算。

2） 作為一種半解析的數(shù)值方法，比例邊界有限元法具有精度高、 計(jì)算量小和數(shù)學(xué)處理方便等優(yōu)點(diǎn)，在拉普拉斯邊值問(wèn)題的求解中同有限元和邊界元相比更具優(yōu)勢(shì)。

3） 數(shù)值算例驗(yàn)證了本文方法在求解穩(wěn)態(tài)非線性熱傳導(dǎo)問(wèn)題時(shí)具有較高的計(jì)算精度和計(jì)算效率。