瞬態(tài)熱傳導(dǎo)問題的精細(xì)積分邊界元法

周楓林， 袁小涵， 欽宇， 潘先云， 余江鴻

(湖南工業(yè)大學(xué)機(jī)械工程學(xué)院， 株洲 412007)

瞬態(tài)熱傳導(dǎo)問題存在于各種工程領(lǐng)域，例如機(jī)械、水利、化工、能源、航天航空等，并隨著瞬態(tài)熱傳導(dǎo)問題的廣泛應(yīng)用，傳熱學(xué)也在迅速發(fā)展，從剛開始的理論與實(shí)驗(yàn)結(jié)合的分析方法到各式各樣的數(shù)值分析方法。在實(shí)際問題中，由于結(jié)構(gòu)形狀以及邊界條件的復(fù)雜性，很難獲得溫度分布的解析解[1]，數(shù)值分析方法已成為有效解決問題的方法之一。

根據(jù)傅里葉熱傳導(dǎo)定律，瞬態(tài)熱傳導(dǎo)問題通常由拋物線偏微分方程(PDF)控制。瞬態(tài)熱傳導(dǎo)問題最常用的數(shù)值方法是有限差分法(finite difference method，F(xiàn)DM)[2]，因?yàn)樗谇蠼馑矐B(tài)熱傳導(dǎo)問題中顯得非常自然。除了FDM外，還有發(fā)展了許多求解瞬態(tài)熱傳導(dǎo)問題的數(shù)值計(jì)算方法，如有限體積法(finite volume method，F(xiàn)VM)[3]、有限元法(finite element method，F(xiàn)EM)[4]、無網(wǎng)格法(meshless method，MM)[5-9]和邊界元法(boundary element method，BEM)[10]，其中，邊界元法是一種只涉及邊界網(wǎng)格化的數(shù)值計(jì)算方法，優(yōu)點(diǎn)是可節(jié)省計(jì)算時(shí)間和存儲(chǔ)量，計(jì)算精度比較高，邊界區(qū)域的形狀可以任意，求解無限大區(qū)域的熱傳導(dǎo)問題比較理想[11]，但是因?yàn)閱栴}的復(fù)雜程度，一般找不到問題的基本解，這時(shí)邊界積分方程中會(huì)出現(xiàn)與內(nèi)部熱源等相關(guān)域積分。為避免域積分的直接計(jì)算，Yang等[12-13]、Feng等[14]在邊界元法中引入了徑向積分方法，將域積分轉(zhuǎn)化為邊界積分。周楓林等[15]將雙互易邊界元法運(yùn)用于標(biāo)量波傳播問題。Guo等[16]采用了三重互易法來避免瞬態(tài)熱傳導(dǎo)問題中的域積分。對(duì)于雙互易法的精度問題，黃誠(chéng)斌等[17]探討了域內(nèi)點(diǎn)的布置位置對(duì)計(jì)算精度的影響，并提出了可讓精度得到進(jìn)一步提高的一種自適應(yīng)布點(diǎn)方法。

在邊界元法的許多計(jì)算運(yùn)用中，通常采用有限差分格式對(duì)時(shí)間進(jìn)行離散，此時(shí)計(jì)算結(jié)果的精度與時(shí)間步長(zhǎng)的選取有著密切的聯(lián)系[18]，雖然有限差分邊界元法是條件穩(wěn)定的，但是在大時(shí)間步長(zhǎng)上，數(shù)值結(jié)果往往不穩(wěn)定。為了避免解決上述問題，Yu等[19]將精細(xì)時(shí)域展開法與雙互易邊界元法相耦合，避免了直接計(jì)算域積分；Yao等[20]將精細(xì)積分方法[21]運(yùn)用于瞬態(tài)熱傳導(dǎo)的邊界元模型中；Yao等[22]采用精細(xì)積分邊界元法解決了雙重相變問題；陳豪龍等[23]將雙互易邊界元法和精細(xì)積分方法用于求解二維瞬態(tài)熱傳導(dǎo)問題。

對(duì)于雙互易精細(xì)積分方法在三維熱傳導(dǎo)問題上還有待深入研究。為此，提出一種雙互易法(dual reciprocity method，DRM)與精細(xì)積分法(precise integration method，PIM)相耦合的方法解決含內(nèi)部熱源的三維瞬態(tài)熱傳導(dǎo)問題。首先將瞬態(tài)熱傳導(dǎo)問題視為準(zhǔn)穩(wěn)態(tài)問題，將溫度對(duì)時(shí)間的導(dǎo)數(shù)項(xiàng)視為等效熱源。通過DRM將熱源的域積分轉(zhuǎn)化為邊界積分。對(duì)空間變量進(jìn)行離散后，轉(zhuǎn)化為關(guān)于域內(nèi)溫度分布的初值問題，最后運(yùn)用PIM解決該初值問題，在計(jì)算了域內(nèi)節(jié)點(diǎn)的溫度后，再通過邊界積分方程計(jì)算邊界的溫度和熱通量。通過邊界元方法對(duì)結(jié)構(gòu)進(jìn)行熱傳導(dǎo)分析，給工程結(jié)構(gòu)的改進(jìn)提供一種數(shù)據(jù)支持，并促進(jìn)邊界元分析方法在熱傳導(dǎo)問題方面的應(yīng)用。

1 瞬態(tài)熱傳導(dǎo)的雙互易邊界元法

1.1 控制方程及其邊界條件

具有內(nèi)部熱源的各向同性介質(zhì)瞬態(tài)常系數(shù)熱傳導(dǎo)問題可表示為

(1)

1.2 邊界積分方程

對(duì)于瞬態(tài)熱傳導(dǎo)問題邊界元法，控制方程的邊界化是關(guān)鍵步驟，借助加權(quán)余量法和三維穩(wěn)態(tài)熱傳導(dǎo)問題的基本解，將控制方程轉(zhuǎn)換為邊界積分方程，可表示為

(2)

對(duì)于三維問題有

(3)

1.3 徑向基函數(shù)逼近

雙互易邊界元法處理瞬態(tài)熱傳導(dǎo)問題的關(guān)鍵在于對(duì)相關(guān)項(xiàng)的逼近，用徑向基函數(shù)(radial basis function,RBF)對(duì)已知函數(shù)(熱源密度函數(shù))和未知函數(shù)(與時(shí)間相關(guān)的溫度函數(shù))進(jìn)行近似表示。取徑向基函數(shù)fi(x)為

(4)

式(4)中：rdis為RBF插值點(diǎn)到源點(diǎn)的距離；s為形狀參數(shù)；i為內(nèi)部節(jié)點(diǎn)和邊界節(jié)點(diǎn)數(shù)之和。

(5)

故之前所提到的已知函數(shù)和未知函數(shù)可表示為

(6)

(7)

式中：Q(x，t)為已知的與位置相關(guān)的熱源密度函數(shù)；αi(t)、βi(t)為未知系數(shù)。

系數(shù)向量β與α可表示為

(8)

F-1Q(t)=β(t)

(9)

并將式(5)、式(6)和式(7)代入式(2)可得

(10)

令：

(11)

然后與推導(dǎo)式(2)相類似，運(yùn)用分部積分法將式(10)中的域積分轉(zhuǎn)化為邊界積分，可表示為

(12)

1.4 邊界離散

式(12)中只涉及邊界積分，所以在邊界離散后，可得到方程的離散形式為

(13)

式中：j為邊界單元個(gè)數(shù)；m為單元中的節(jié)點(diǎn)數(shù)，采用八節(jié)點(diǎn)二次面單元，m=8；xm為單元中節(jié)點(diǎn)的三維空間坐標(biāo)；Nm為八節(jié)點(diǎn)面單元的形函數(shù)。

然后在每個(gè)邊界的插值點(diǎn)上配置場(chǎng)點(diǎn)，將上述方程寫為矩陣形式

(14)

δijC(yi)，yi∈Γ

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

將式(8)和式(9)代入式(14)可得

(24)

特別的是，由于Q(x,t)一般為已知的與位置相關(guān)的熱源密度函數(shù)，所以向量Q(t)為已知向量。

1.5 形成微分方程組

式(15)為半離散系統(tǒng)，將其改寫為

(25)

假設(shè)邊界類型不隨時(shí)間改變，可將邊界條件的表達(dá)式改寫為

c1u(x,t)+c2q(x,t)=c3，x∈Γ

(26)

式(26)中：c3為邊界條件；c1和c2為常數(shù)值，其值可由邊界類型來確定，則有

(27)

所有邊界節(jié)點(diǎn)上的邊界條件可統(tǒng)一表示為

(28)

然后合并式(25)和式(28)可得

(29)

則邊界溫度和邊界熱通量可表示為

(30)

令

(31)

(32)

則式(30)簡(jiǎn)寫為

(33)

然后考慮內(nèi)部點(diǎn)，則有

(34)

(35)

(36)

(37)

(38)

將式(33)代入式(34)得

(39)

令

(40)

(41)

式(39)可簡(jiǎn)寫為

(42)

式(42)為一階常系數(shù)微分方程組。

2 精細(xì)積分方法

對(duì)于式(42)的一階常系數(shù)微分方程組，可用已知的初始條件求解。

u(x,t0)=u0(x), ?x∈Ω

(43)

式(42)的通解形式為

(44)

式(44)中：τ為積分時(shí)間變量。

(45)

式(45)中：I為單位矩陣。

式(33)稱為矩陣指數(shù)函數(shù)，為使得計(jì)算成本最小，采用逐步計(jì)算的方法求解式(44)。

(46)

因而，式(46)可化簡(jiǎn)為

(47)

這是一個(gè)代數(shù)推進(jìn)方程，依據(jù)式(43)的初始條件，通過該方程計(jì)算，可以逐步得到域內(nèi)溫度，在計(jì)算完域內(nèi)溫度后，通過式(33)和式(39)可計(jì)算出邊界上的量。

矩陣指數(shù)函數(shù)的計(jì)算精確度決定了整個(gè)系統(tǒng)的一個(gè)計(jì)算精度，由于矩陣指數(shù)函數(shù)涉及大量的矩陣乘法，所以計(jì)算非常耗時(shí)，為了準(zhǔn)確計(jì)算該矩陣指數(shù)函數(shù)，采取一種精細(xì)計(jì)算方法。

式(45)實(shí)際上可以看作是泰勒展開式，其收斂速度在很大程度上取決于時(shí)間步長(zhǎng)。理論上，級(jí)數(shù)收斂于任意時(shí)間步長(zhǎng)，然而，與較大的時(shí)間步長(zhǎng)相比，較小的時(shí)間步長(zhǎng)下，式(45)可以用很少的截?cái)囗?xiàng)來計(jì)算。在精細(xì)積分方法中，矩陣指數(shù)函數(shù)中的時(shí)間步長(zhǎng)被細(xì)分為

(48)

式(48)中：M為細(xì)分參數(shù)。

運(yùn)用泰勒級(jí)數(shù)展開，可表示為

=I+Er

(49)

=(I+Er)(I+Er)

=I+2Er+ErEr

(50)

然后重新賦值Er，使得

Er=2Er+ErEr

(51)

可得

=(I+Er)(I+Er)

=I+2Er+Er×Er

(52)

然經(jīng)過M次重復(fù)計(jì)算和重新賦值Er后，最終可得

eAΔt=I+Er

(53)

在整個(gè)計(jì)算過程中，Er的計(jì)算涉及了p個(gè)矩陣乘法，eAΔt的計(jì)算涉及了M個(gè)矩陣乘法，特別的是，假如時(shí)間步長(zhǎng)為常量，那么eAΔt只會(huì)被計(jì)算一次。并且在時(shí)間步推進(jìn)中，只涉及矩陣的乘法運(yùn)算。

3 數(shù)值算例

為驗(yàn)證所提出方法的有效性、準(zhǔn)確性和穩(wěn)定性，給出3個(gè)數(shù)值算例，運(yùn)用雙互易精細(xì)積分方法獲得的計(jì)算結(jié)果與解析解進(jìn)行對(duì)比。用于評(píng)估的相對(duì)誤差的表達(dá)式為

(54)

在這3個(gè)算例中，對(duì)于雙互易方法，使用的近似函數(shù)為前面所提到的式(6)，形狀參數(shù)s為1，對(duì)于精細(xì)積分方法中的細(xì)分參數(shù)M和截?cái)鄥?shù)p，3個(gè)案例都取10和6；材料的熱導(dǎo)率、熱容和密度分別給定為1 W/(m·C)、1 /(kg·C)和1 kg/m-3。

3.1 狄利克雷邊界條件下的環(huán)面結(jié)構(gòu)

在這個(gè)數(shù)值算例中，分析圖1所示的環(huán)形結(jié)構(gòu)上的瞬態(tài)熱傳導(dǎo)問題，圓環(huán)的外徑和內(nèi)徑分別為26 m和14 m。

狄利克雷邊界條件為

(55)

熱源項(xiàng)為

Q(x,t)=2 W/m-3

(56)

初始條件為

(57)

解析解為

(58)

在整個(gè)計(jì)算過程中，如圖1所示，總共劃分了750個(gè)矩形單元，2 381個(gè)邊界節(jié)點(diǎn)，1 473個(gè)徑向基函數(shù)插值點(diǎn)。

考慮0～50 s的時(shí)間段，并且選取了6個(gè)時(shí)間步長(zhǎng)用于計(jì)算，即0.1、0.2、0.5、1.0、2.5、5 s，圖2展示了在不同的時(shí)間步長(zhǎng)下，內(nèi)部點(diǎn)計(jì)算溫度與解析解的相對(duì)誤差，可以看出，隨著時(shí)間的變化，相對(duì)誤差趨近于0，逐漸收斂，并且在時(shí)間步比較小的情況下，即使是接近初始時(shí)刻，計(jì)算結(jié)果也有很高的精確度。在本案例中，最佳時(shí)間步長(zhǎng)為最小的時(shí)間步0.1 s，并且在同一時(shí)刻，相對(duì)誤差大小根據(jù)時(shí)間步長(zhǎng)大小排序，越小的時(shí)間步，相對(duì)誤差越小。

同時(shí)，對(duì)該模型進(jìn)行有限元分析，運(yùn)用軟件為workbench，計(jì)算模型如圖3所示，設(shè)置單元大小為0.5 m，節(jié)點(diǎn)數(shù)為169 388，單元數(shù)為119 334。

此時(shí)考慮一般情況，即解析解不易求得的情況，給定狄利克雷邊界條件為

u(x,t)=60 ℃，x∈Γ

(59)

圖1 圓環(huán)邊界元模型Fig.1 Ring boundary element model

圖2 溫度的相對(duì)誤差隨時(shí)間變化情況Fig.2 Variation of relative error of temperature along time

熱源項(xiàng)為

Q(x,t)=1 W/m3

(60)

初始條件為

u(x,t0)=u0(x)=0 ℃，x∈Ω

(61)

在兩個(gè)模型中，取同樣的內(nèi)部點(diǎn)A(7.927，-7.444，0.634 8)，運(yùn)用相同的時(shí)間步長(zhǎng)0.1 s，考慮0～10 s的時(shí)間段，將邊界元溫度計(jì)算結(jié)果與有限元分析溫度結(jié)果進(jìn)行對(duì)比，如圖4所示。

圖3 圓環(huán)的有限元模型Fig.3 Finite element model of ring

圖4 PI-DRBEM與FEM結(jié)果比較Fig.4 Comparison of results between PI-DRBEM and FEM

圖4表明，雙重互易法和精細(xì)積分法耦合擁有很高的精度，論證了該方法的有效性和準(zhǔn)確性。

3.2 混合條件下的山型結(jié)構(gòu)

此算例考慮在混合邊界條件下，一個(gè)山型結(jié)構(gòu)的瞬態(tài)熱傳導(dǎo)問題，三維模型及其尺寸如圖5所示。邊界劃分如圖6所示。

圖5 山型結(jié)構(gòu)幾何示意圖Fig.5 Geometric diagram of mountain structure

Γ1和Γ2為諾依曼型邊界圖6 邊界劃分幾何示意圖Fig.6 Geometric diagram of boundary division

邊界條件為

(62)

q(x,t)=62x2+16t，x1,x2,x3∈Γ1

(63)

其余邊界Γ3為狄利克雷型邊界，邊界條件為

8t(x1+2x2+1.5x3)

(64)

熱源項(xiàng)為

Q(x,t)=2x1+4x2+3x3

(65)

初始條件為

u(x,t0)=u0(x)

(66)

解析解為

1.5x3)，x1,x2,x3∈Ω

(67)

在此算例的計(jì)算中，總共劃分了500個(gè)連續(xù)的隨機(jī)四邊形單元，1 920個(gè)邊界節(jié)點(diǎn)，164個(gè)徑向基函數(shù)插值點(diǎn)。使用不同的時(shí)間步長(zhǎng)來驗(yàn)證方法對(duì)時(shí)間變量的收斂性,相對(duì)誤差結(jié)果如圖7所示。

在本算例中，考慮同樣的6個(gè)時(shí)間步長(zhǎng)，在混合邊界條件下，時(shí)間步長(zhǎng)趨近于0時(shí)，計(jì)算結(jié)果從一開始就具有很高的精度。稍大的時(shí)間步長(zhǎng)中，盡管剛開始誤差偏大，但隨著時(shí)間的推移，相對(duì)誤差快速減小?？梢钥闯?，雙互易精細(xì)積分法的精確性和有效性。

3.3 狄利克雷邊界條件下的彎管結(jié)構(gòu)

此算例考慮空心彎管結(jié)構(gòu)，模型結(jié)構(gòu)及尺寸如圖8所示。所有表面為狄利克雷邊界，其條件為

圖7 溫度的相對(duì)誤差隨時(shí)間變化情況Fig.7 Variation of relative error of temperature along time

R為半徑圖8 彎管幾何結(jié)構(gòu)示意圖Fig.8 Geometric diagram of elbow

(68)

熱源項(xiàng)、初始條件、解析解與算例2保持一致。在整個(gè)計(jì)算過程中，總共劃分了448個(gè)矩形單元，1 718個(gè)邊界節(jié)點(diǎn)，326個(gè)徑向基函數(shù)插值點(diǎn)，插值點(diǎn)分布如圖9所示。

圖9 彎管RBF插值點(diǎn)分布Fig.9 Distribution of RBF interpolation points of elbow

在本算例中，同樣考慮0～50 s時(shí)間段，結(jié)果如圖10所示，可以看出，相對(duì)誤差在任何一個(gè)時(shí)間步中都是隨著時(shí)間推移而減小，并且時(shí)間步長(zhǎng)取0.1 s時(shí)，內(nèi)部點(diǎn)熱通量計(jì)算結(jié)果的相對(duì)誤差最小，并且在初始時(shí)刻附近也具有較小的相對(duì)誤差。有趣的是，內(nèi)部點(diǎn)熱通量的計(jì)算精度，并不是隨著時(shí)間步長(zhǎng)減小而變高，在選取的6個(gè)時(shí)間步中，在取中間大小的時(shí)間步1.0 s時(shí)，相對(duì)誤差在任何時(shí)刻都是最大的。

圖10 熱通量隨時(shí)間的相對(duì)誤差變化Fig.10 Relative error variation of heat flux along time

4 結(jié)論

提出了求解含有內(nèi)部熱源的三維瞬態(tài)熱傳導(dǎo)問題的雙互易精細(xì)積分方法。在處理該熱學(xué)問題中，首先運(yùn)用雙互易法處理域積分，再利用精細(xì)積分方法處理所得到的初始條件問題。得出如下結(jié)論。

(1)雙互易邊界元法和精細(xì)積分方法的耦合具有較高的精度和穩(wěn)定性。

(2)對(duì)于不同的結(jié)構(gòu)體和邊界條件，并不是時(shí)間步長(zhǎng)越小越好。存在時(shí)間步長(zhǎng)減小，但是計(jì)算結(jié)果的相對(duì)誤差卻增大的情況，但通常情況下，時(shí)間步長(zhǎng)越接近零，越有可能獲得高精度結(jié)果。

(3)在無解析解的情況下，該方法依舊有很高的精確性。

不足是該方法的計(jì)算效率有待進(jìn)一步提高，特別的是，在第一個(gè)時(shí)間步內(nèi)，需要進(jìn)行大量的矩陣與矩陣的乘法運(yùn)算及矩陣求逆運(yùn)算。而在后續(xù)時(shí)間步中，只需進(jìn)行已有矩陣與向量的乘法運(yùn)算，容易實(shí)現(xiàn)并行計(jì)算。