藍(lán)林華，富明慧，劉祚秋

(中山大學(xué)應(yīng)用力學(xué)與工程系，廣東 廣州 510275)

熱傳導(dǎo)問題一直是熱門問題，隨著科學(xué)技術(shù)的發(fā)展，更多的領(lǐng)域和方向與之緊密相關(guān)，此問題的研究也越來越受到關(guān)注。在過去的幾十年里，人們進(jìn)行了大量的研究工作，提出了各種近似求解的方法(有限差分法、邊界元法、多重網(wǎng)格法、變分法，無網(wǎng)格法等)。劉繼軍等[1]針對(duì)二維熱傳導(dǎo)方程問題構(gòu)造了一種三層顯式格式并進(jìn)行了數(shù)值計(jì)算。Ochiai[2]采用邊界元方法討論了含熱源的穩(wěn)態(tài)熱傳導(dǎo)問題，同時(shí)用區(qū)域積分法改進(jìn)了邊界元方法并把它運(yùn)用到功能梯度材料的熱傳導(dǎo)問題中。李壽佛等[3]用無網(wǎng)格方法研究了二維、三溫?zé)醾鲗?dǎo)方程。陳光南等[4]用變分原理在不規(guī)則結(jié)構(gòu)網(wǎng)格上建立熱流通量形式的差分格式,研究了二維熱傳導(dǎo)方程的差分?jǐn)?shù)值模擬。Roy等[5]用多重尺度法計(jì)算了雙曲熱傳導(dǎo)問題。Chen等[6]結(jié)合直接法和伴隨矩陣法得到功能梯度材料穩(wěn)態(tài)和瞬態(tài)熱傳導(dǎo)問題的靈敏度方程并進(jìn)行了分析，隨后利用精細(xì)積分法求解了功能梯度材料的瞬態(tài)傳熱問題。精細(xì)積分法最初是應(yīng)用于結(jié)構(gòu)動(dòng)力分析的方法[7-8]，但由于具有精度高、穩(wěn)定性好等優(yōu)點(diǎn)，因此迅速被推廣到熱傳導(dǎo)問題[9-10]。本文針對(duì)單層或者多層介質(zhì)上的二維瞬態(tài)熱傳導(dǎo)方程開展研究，對(duì)空間進(jìn)行離散，并對(duì)空間離散后的拉普拉斯差分算子由塊三角矩陣轉(zhuǎn)化為塊對(duì)角形式，結(jié)合精細(xì)積分算法，建立了求解層合結(jié)構(gòu)瞬態(tài)熱傳導(dǎo)問題的一種有效方法。算例表明，它具有良好的數(shù)值穩(wěn)定性，而且精度高、收斂速度快，便于運(yùn)用。

1 單層介質(zhì)的二維瞬態(tài)熱傳導(dǎo)問題

考慮矩形區(qū)域上二維瞬態(tài)熱傳導(dǎo)方程的混合問題

(1)

其中φ(x,y,t)是所要求的溫度分布函數(shù)，a為熱擴(kuò)散系數(shù)，Ω為矩形區(qū)域內(nèi)部(lx×ly)，?Ω為矩形區(qū)域邊界，g(x,y,t)，φ(x,y)，φ0(x,y)都是已知的連續(xù)函數(shù)。

用圖1所示均勻網(wǎng)格將定解區(qū)域離散，陰影部分為一矩形小區(qū)域單元(h1×h2)，其中h1=lx/(m+1),h2=ly/(n+1)分別是x和y方向的步長(zhǎng)。

圖1 二維熱傳導(dǎo)矩形域的網(wǎng)格劃分

(2)

記

式中δij為Kronecker_δ符號(hào)。利用上述記號(hào)，式(2)可表達(dá)為

(4)

其中

(5)

式中In×n是n階單位矩陣，D是n階三對(duì)角矩陣。

容易求得，矩陣D的特征值λ及相應(yīng)的特征向量分別為[11]

其中k=1,…,n，η=π/(n+1)。

Q=(q1,q2,…,qn)=

由特征矩陣的性質(zhì)，有

Q-1CQ=Λ

(6)

由于Q是正交而且對(duì)稱的，有

Q-1=QT=Q

(7)

由式(6)，式(4)可化為

(8)

其中

(i=1,2,…，m)

(9)

式(8)可記為

(10)

其中

矩陣H具有塊三對(duì)角形式，這對(duì)于建立高效算法是十分有利的。

初始條件也相應(yīng)地化為

(11)

式(10)的解可表達(dá)為

(12)

2 層合結(jié)構(gòu)的瞬態(tài)熱傳導(dǎo)問題

圖2所示為一(k+1)層矩形層合結(jié)構(gòu)差分網(wǎng)格示意圖，沿x方向均勻離散，y方向分段均勻離散，黑粗線是每?jī)蓪硬牧现g的分界線。

圖2 層合結(jié)構(gòu)的差分網(wǎng)格示意圖

為方便推導(dǎo)，本文假設(shè)分界面的熱接觸良好，沒有接觸熱阻。其導(dǎo)熱微分方程(不包括分界線)可分段寫為

用與方程(2)類似的符號(hào)記法，方程(13)的差分形式可寫為

(14)

其中h0=l0/(m+1),hr=lr/(nr-nr-1) 分別是x和y方向的步長(zhǎng)。

對(duì)于第r條分界線上的某一差分點(diǎn)(xp,ynr)，取如圖2b陰影部分的一小單元h0×(hr+hr+1)/2， 根據(jù)Fourier定律和能量守恒定理，可以直接寫出其差分形式

(15)

記

(16)

式中δij為Kronecker_δ符號(hào)。

利用上述記號(hào)，式(15)、(16)可合并表達(dá)為

(17)

其中

(18)

式中Im是m階單位方陣，Cr、Cnr、Dr、Dnr均是m階三對(duì)角矩陣。類似于上一節(jié)，式(17)的系數(shù)矩陣的每一個(gè)子塊同樣可以化為對(duì)角矩陣，特別地，當(dāng)ar=ar+1，hr=hr+1時(shí)，式(18)可化為式(5)形式，即復(fù)合疊層材料變成單層介質(zhì)材料。同樣，上述方法同樣可推廣到三維的情況。

3 算 例

圖3 是用兩種方法得到的溫度場(chǎng)等值線分布圖，其中圖3(a)的在t=1 s后的結(jié)果，圖3(b)是t=10 s后的結(jié)果。此結(jié)果中有限元法的時(shí)間步長(zhǎng)取為Δt=0.005 s，可以看到，當(dāng)時(shí)間步長(zhǎng)取得合理時(shí)，兩種方法的結(jié)果吻合得很好。

圖3 溫度場(chǎng)等值線分布圖

表1 中心點(diǎn)在t=0.5 s的溫度值

Table 1 The temperature of center point att=0.5 s

時(shí)間步長(zhǎng)/s溫度值有限元法本文方法0.005 99.236 506 912 023 6099.211 297 459 871 740.01 99.261 546 525 305 5599.211 297 459 872 300.02 —99.211 297 459 872 610.1 —99.211 297 459 872 780.5 —99.211 297 459 873 08

表1 給出了中心點(diǎn)在t=0.5 s時(shí)的溫度值，分別用本文方法和有限元方法進(jìn)行了計(jì)算，兩種方法均采用相同的時(shí)間、空間網(wǎng)格劃分。數(shù)字的下劃線部分是兩種方法在不同時(shí)間步長(zhǎng)下得到結(jié)果的不同有效位數(shù)，從結(jié)果可以看出，對(duì)于不同的時(shí)間步長(zhǎng)，即使時(shí)間步長(zhǎng)等于時(shí)間t，本文方法的計(jì)算結(jié)果有長(zhǎng)達(dá)11位的相同有效位數(shù)，而有限元法的結(jié)果對(duì)時(shí)間步長(zhǎng)的變化情況比較敏感，特別是當(dāng)步長(zhǎng)增大到一定程度時(shí)還會(huì)出現(xiàn)數(shù)值不穩(wěn)定的情況(表中的虛橫線)。事實(shí)上，用有限元法計(jì)算二維瞬態(tài)熱傳導(dǎo)方程時(shí)其差分格式分顯式格式和隱式格式，隱式格式能夠保證收斂，但在高維的情況下其計(jì)算量很大，相當(dāng)于顯式格式計(jì)算量的平方陪，故對(duì)于高維的情況一般采用顯式格式(本文的有限元法也采用顯式格式)，但顯式格式為了保證收斂性，時(shí)間、空間的網(wǎng)格劃分須滿足一定條件

表2 有限元法在不同時(shí)間步長(zhǎng)下的穩(wěn)定性分析

表2 給出了有限元法在不同時(shí)間步長(zhǎng)下的穩(wěn)定性分析情況，其中紅色的大于0.5表示該時(shí)間步長(zhǎng)不滿足穩(wěn)定行的充要條件，顯然表2 的分析結(jié)果跟表1是一致的。

而對(duì)于本文采用的精細(xì)積分(PTI)方法，即使時(shí)間步長(zhǎng)取整個(gè)時(shí)間t，而精細(xì)積分算法中的dt=t/220將是一個(gè)很小的量(根據(jù)需要還可以取得更小)，它總能夠滿足上面的穩(wěn)定性條件。也就是說，本文的求解過程可以無條件地滿足算法的穩(wěn)定性要求。

4 結(jié) 論

本文針對(duì)單層或者多層介質(zhì)上的二維瞬態(tài)熱傳導(dǎo)方程，給出了對(duì)空間的詳細(xì)離散過程，并對(duì)空間離散后的拉普拉斯差分算子由塊三角矩陣轉(zhuǎn)化為塊對(duì)角形式，隨后結(jié)合已有的精細(xì)積分算法，建立了求解層合結(jié)構(gòu)瞬態(tài)熱傳導(dǎo)方程的一種有效方法。算例表明，本文方法不僅具有良好的精度(相同的有效位數(shù)長(zhǎng)達(dá)11位)，還可以無條件地滿足算法的穩(wěn)定性要求，并且該方法還可以推廣到三維的情況，因此可以成功地解決層合結(jié)構(gòu)的瞬態(tài)熱傳導(dǎo)問題，為國(guó)內(nèi)、外同類研究提供一定的參考。

參考文獻(xiàn)：

[1]劉繼軍. 二維熱傳導(dǎo)方程的三層顯式差分格式[J]. 應(yīng)用數(shù)學(xué)和力學(xué),2003,24(5):537-544.

[2]OCHIAI Y. Two-dimensional steady heat conduction in functionally gradient materials by triple-reciprocity boundary element method[J]. Eng Analysis with Boundary Elements,2004,28(12):1445-1453.

[3]李壽佛,張瑗,劉玉珍. 熱傳導(dǎo)方程的一類無網(wǎng)格方法[J]. 計(jì)算物理,2007,24(5):573-580.

[4]陳光南,張永慧. 基于變分原理的二維熱傳導(dǎo)方程差分格式[J]. 計(jì)算物理,2002,19(4):299-304.

[5]ROY S,VASUDEVA MURTHY A S,KUDUNATTI R B. A numerical method for the hyperbolic-heat conduction equation based on multiple scale technique[J]. Applied Numerical Mathematics,2009,59:1419-1430.

[6]GU Y X,CHEN B S,ZHANG H W,et al. A sensitivity analysis method for linear and nonlinear transient heat conduction with precise time integration[J]. Struct Multidisc Optim,2002,24:23-37.

[7]鐘萬勰. 一類非線性熱傳導(dǎo)方程的單點(diǎn)精細(xì)積分[J]. 自然科學(xué)進(jìn)展,1996,6(3):313-316.

[8]富明慧,林敬華. 精細(xì)積分法在非線性動(dòng)力學(xué)問題中的應(yīng)用[J]. 中山大學(xué)學(xué)報(bào)：自然科學(xué)版,2008,47(3):1-5.

[9]藍(lán)林華,富明慧,程正陽. 功能梯度材料瞬態(tài)熱傳導(dǎo)問題的降維精細(xì)積分法[J]. 固體力學(xué)學(xué)報(bào),2010,31(4): 406-410.

[10]藍(lán)林華,富明慧,高文樂. 功能梯度材料穩(wěn)態(tài)熱傳導(dǎo)方程的分層精細(xì)指數(shù)法[J]. 中山大學(xué)學(xué)報(bào)：自然科學(xué)版,2011,50(4): 1-6.

[11]李榮華,馮果沈. 微分方程數(shù)值解法[M]. 北京:人民教育出版社,1980.

[12]富明慧,劉祚秋,林敬華. 一種廣義精細(xì)積分法[J]. 力學(xué)學(xué)報(bào),2007,39(5):672-677.

[13]富明慧,梁華力. 一種改進(jìn)的精細(xì)-龍格庫塔法[J].中山大學(xué)學(xué)報(bào)：自然科學(xué)版,2009,48(5):1-5.

[14]戴嘉尊. 微分方程數(shù)值解法[M]. 南京:東南大學(xué)出版社,2002.