王曉飛，李興國(guó)，任紅萍

(太原科技大學(xué) 應(yīng)用科學(xué)學(xué)院,太原 030024)

瞬態(tài)熱傳導(dǎo)問(wèn)題是工程上的一類(lèi)很重要的問(wèn)題,但一般這類(lèi)溫度場(chǎng)問(wèn)題的解析解是很難找到的,所以提出來(lái)了數(shù)值解。Belytschko等人根據(jù)改進(jìn)的擴(kuò)散單元法,提出了EFG方法[1]。由于無(wú)網(wǎng)格方法的計(jì)算量大,程玉民等人提出了新的構(gòu)造形函數(shù)的方法。同時(shí)將復(fù)變量的相關(guān)知識(shí)引入了MLS法,就是我們所說(shuō)的復(fù)變量MLS法[2-5]。隨后,插值型復(fù)變量MLS法被提出來(lái)。由于相比其它方法來(lái)說(shuō),插值型復(fù)變量MLS法的試函數(shù)中的待定系數(shù)少,所以二維問(wèn)題的復(fù)變量EFG方法可選取較少的節(jié)點(diǎn)。根據(jù)上述所提出的方法,在任何的場(chǎng)點(diǎn),它的影響域中所包含的節(jié)點(diǎn)的最小數(shù)量就會(huì)大大地減少,然后可以極大地減少在求解域中所要選擇的節(jié)點(diǎn)數(shù)量。和基于MLS法的無(wú)網(wǎng)格方法相比來(lái)說(shuō),在具有相同的節(jié)點(diǎn)分布和相近的計(jì)算精度時(shí),插值型復(fù)變量MLS法有較好的計(jì)算精度以及計(jì)算速度。陳麗和程玉民等提出了通過(guò)復(fù)變量重構(gòu)核粒子法來(lái)求解二維彈性力學(xué)的方法[6]。Chen和Cheng[7]采用復(fù)變量重構(gòu)核粒子的方法來(lái)求解了瞬態(tài)熱傳導(dǎo)問(wèn)題,它的優(yōu)點(diǎn)是二維問(wèn)題可以用一維的基函數(shù)來(lái)解決。Sigh等用無(wú)網(wǎng)格方法分析了半無(wú)限體的非穩(wěn)態(tài)傳熱導(dǎo)[8]。任紅萍等[9-10]通過(guò)插值型MLS法給出了插值型CVEFG方法以及改進(jìn)的邊界無(wú)單元法。張贊[11]等用改進(jìn)的Galerkin方法解決了三維瞬態(tài)熱傳導(dǎo)問(wèn)題,并且討論了尺度參數(shù)和節(jié)點(diǎn)數(shù)量等因素和數(shù)值解的關(guān)系。張國(guó)達(dá)[12]等研究了瞬態(tài)熱傳導(dǎo)問(wèn)題的插值型無(wú)單元伽遼金方法及誤差。

插值型復(fù)變量MLS法構(gòu)成了插值型復(fù)變量無(wú)單元伽遼金方法的基本原理。在此基礎(chǔ)上,本文以插值型復(fù)變量MLS法作為基礎(chǔ),由伽遼金積分弱形式,給出了二維瞬態(tài)熱傳導(dǎo)的第一類(lèi)邊值問(wèn)題的插值型CVEFG方法,并給出了相應(yīng)的推導(dǎo)過(guò)程。為了證明此方法的有效性,分別對(duì)兩個(gè)瞬態(tài)熱傳導(dǎo)的第一類(lèi)邊值問(wèn)題通過(guò)使用插值型復(fù)變量EFG方法進(jìn)行了數(shù)值求解。

1 控制方程的伽遼金積分弱形式

通過(guò)用插值型復(fù)變量EFG法推導(dǎo)二維瞬態(tài)熱傳導(dǎo)問(wèn)題的方程。通常,導(dǎo)熱率在空間上變化的各向同性材料的二維瞬態(tài)熱傳導(dǎo)問(wèn)題的控制方程可以表示為:

(1)

初始條件為:

T(x,y,0)=T0(x,y)∈Γ1

(2)

邊界條件為:

(3)

(4)

方程(1),(4)的等效積分的弱形式:

(5)

泛函Π(T)可表示為:

(6)

令δΠ=0，可得:

(7)

其中δ是變分運(yùn)算符,L是一個(gè)微分算子矩陣,

(8)

(9)

由于本文的形函數(shù)是由插值型復(fù)變量MLS法來(lái)構(gòu)造的,故而形函數(shù)符合Delta的特性,因此在求解方程的過(guò)程中可以直接施加本質(zhì)邊界條件，不需額外進(jìn)行處理。

2 積分弱形式的離散化

首先把空間域Ω離散為M個(gè)節(jié)點(diǎn)ZI,ZI=1,2,3…M,節(jié)點(diǎn)ZI的溫度為:

TI(t)=T(zI,t)

(10)

域內(nèi)任意場(chǎng)點(diǎn)在某時(shí)刻t的溫度T(z,t)采用影響域覆蓋場(chǎng)點(diǎn)的節(jié)點(diǎn)溫度進(jìn)行逼近。由于在任意時(shí)刻T(z,t)和TI(t)都為標(biāo)量,由插值型復(fù)變量MLS法的逼近函數(shù)表達(dá)式,可以得出任意場(chǎng)點(diǎn)z在時(shí)刻t的溫度可表示為:

T(z,t)=Re[Φ(z)T(t)]=

(11)

其中n為影響域中覆蓋z的節(jié)點(diǎn)的數(shù)目。

Φ(z)=(Φ1(z),Φ2(z),…Φn(z))=

(12)

T(t)=(T1(t),T2(t)，…，Tn(t))T

(13)

(14)

(15)

(16)

(17)

(18)

W(z)=

(19)

vT(z)=(v(z-z1),v(z-z2)…v(z-zn))

(20)

(21)

(22)

由式(8)和(11)可得:

(23)

其中：

B(z)=(B1(z),B2(z)…Bn(z))

(24)

(25)

將(12)和(24)代入得:

(26)

其中:

(27)

方程(26)的第一項(xiàng)為:

(28)

(29)

C=[CIJ]是M×M的矩陣

(30)

方程(26)的第二項(xiàng)為:

其中:

(32)

K=[KIJ]是M×M的矩陣

(33)

方程(26)的第三項(xiàng)為:

δTT·F(1)

(34)

(35)

(36)

方程(26)的第四項(xiàng)為:

(37)

(38)

(39)

將式(28),(31),(34),(37)代入方程(26)得到:

(40)

由δTT的任意性，可得：

(41)

其中:

F=F(1)+F(2)

(42)

將方程(7) 離散為僅存在以時(shí)間為獨(dú)立變量的常微分方程(41).

對(duì)時(shí)間域的離散:

(43)

當(dāng)取θ=0時(shí),方程(43)能寫(xiě)為:

(44)

(45)

其中：

Tn+1=T((n+1)Δt),Tn=T(nΔt)

(46)

Fn+1=F((n+1)Δt),Fn=F(nΔt)

(47)

以上內(nèi)容為用插值型復(fù)變量EFG方法來(lái)求解二維瞬態(tài)熱傳導(dǎo)問(wèn)題的過(guò)程。

3 算法實(shí)施流程

由二維瞬態(tài)熱傳導(dǎo)的第一類(lèi)邊值問(wèn)題,它的插值型復(fù)變量EFG法的步驟如下:

(1)輸入已知的數(shù)據(jù),包括求解域的幾何特性,時(shí)間的步長(zhǎng)和材料常數(shù);

(2)基于已知的第一類(lèi)邊值問(wèn)題,首先建立它的坐標(biāo)系,其次在求解域Ω和它的Γ上分布數(shù)量為M個(gè)的節(jié)點(diǎn),對(duì)節(jié)點(diǎn)進(jìn)行標(biāo)號(hào),建立節(jié)點(diǎn)信息,建立相應(yīng)的節(jié)點(diǎn)坐標(biāo);

(3)生成背景積分網(wǎng)格,并將其帶入數(shù)值積分;

(4)建立邊界積分單元的相關(guān)信息;

(5)Gauss積分點(diǎn)由背景的網(wǎng)格以及邊界生成,并且計(jì)算對(duì)應(yīng)的Gauss積分點(diǎn)的信息(包括Gauss積分點(diǎn)的坐標(biāo)、積分權(quán)重和相應(yīng)的Jacobi行列式|J|);

(6)形成熱容矩陣C、熱傳導(dǎo)矩陣K和F的第一項(xiàng)F(1);

A.循環(huán)每個(gè)背景積分網(wǎng)格;

1)循環(huán)高斯積分點(diǎn)(背景網(wǎng)格內(nèi));

2)若高斯積分點(diǎn)在Ω內(nèi),則進(jìn)行步驟(3)到(6)否則直接進(jìn)行步驟(6);

3)確定其影響域覆蓋的當(dāng)前高斯積分點(diǎn)zQ的節(jié)點(diǎn);

4)計(jì)算影響域覆蓋當(dāng)前的高斯積分點(diǎn)zQ的所有節(jié)點(diǎn)zI處的ΦI(zQ)及其導(dǎo)數(shù)的值;

5)計(jì)算當(dāng)前高斯積分點(diǎn)zQ對(duì)矩陣C、K和列向量F(1)的貢獻(xiàn);

6)結(jié)束循環(huán);

B.結(jié)束背景積分網(wǎng)格的循環(huán)。

(7)在式(38)的基礎(chǔ)上計(jì)算列向量F(2).

(8)將所計(jì)算出的F(1)和上述步驟所得的F(2)進(jìn)行相加,得出F.

(9)對(duì)于已知初值條件的瞬態(tài)熱傳導(dǎo)的第一類(lèi)邊值問(wèn)題，將初值條件代入遞推的關(guān)系式(41)來(lái)近似地給出另一組數(shù)據(jù),可以得到時(shí)間域內(nèi)的各個(gè)時(shí)間的場(chǎng)函數(shù)的數(shù)值T；

(10)由上述步驟得出的T和式(11)可以擬合出Ω內(nèi)任一場(chǎng)點(diǎn)在任意時(shí)刻的溫度數(shù)值解；

(11)輸出溫度值.

4 數(shù)值算例

下面利用剛才建立的二維瞬態(tài)熱傳導(dǎo)問(wèn)題的插值型復(fù)變量EFG法,對(duì)下面的兩個(gè)二維瞬態(tài)熱傳導(dǎo)的第一類(lèi)邊值問(wèn)題進(jìn)行數(shù)值求解。

4.1 算例1

控制方程是：

u,t=u,11+u,22+(1+t2)u+

(2π2-t2-2)e-tsin(πx)cos(πy)x∈Ω,

t>0

(48)

邊界條件為:

u(0,y,t)=u(1,y,t)=0

(49)

u(x,0,t)=-u(x,1,t)=e-tsin(πx)

(50)

初始條件為:

u(x,y,0)=sin(πx)cos(πy)

(51)

其中Ω=[0,1]×[0,1]

該問(wèn)題對(duì)應(yīng)的解析解為:

u(x,y,0)=e-tsin(πx)cos(πy)

(52)

當(dāng)使用插值型復(fù)變量EFG方法進(jìn)行計(jì)算時(shí),區(qū)域內(nèi)的節(jié)點(diǎn)分布選擇為11×11個(gè)節(jié)點(diǎn)。在影響域中的權(quán)函數(shù)的比例參數(shù)選取dmax=2.0.時(shí)間步長(zhǎng)選取Δt=0.001s.圖1為t=0.01時(shí),x=0.5的數(shù)值解和解析解。

圖1 當(dāng)x=0.5的溫度分布

4.2 算例2

控制方程是：

(53)

邊界條件為:

T(0,y,t)=e-2tsiny

(54)

T(1,y,t)=e-2tsin(1+y)

(55)

T(x,0,t)=e-2tsinx

(56)

T(x,1,t)=e-2tsin(1+x)

(57)

初始條件為:

T(x,y,0)=sin(x+y)

(58)

解析解為:

T(x,y,t)=e-2tsin(x+y)

(59)

當(dāng)使用插值型復(fù)變量EFG方法進(jìn)行計(jì)算時(shí),區(qū)域內(nèi)的節(jié)點(diǎn)分布選擇為11×11個(gè)節(jié)點(diǎn)。在影響域中的權(quán)函數(shù)的影響域比例參數(shù)選取dmax=1.5.時(shí)間步長(zhǎng)選取Δt=0.001s.圖2和圖3分別為t=0.01時(shí),x=0.4和y=0.7的數(shù)值解和解析解。

圖2 t=0.01時(shí)，x=0.4的溫度分布

圖3 t=0.01時(shí)，y=0.7的溫度分布

5 結(jié)論

本文通過(guò)使用插值型復(fù)變量MLS法,推導(dǎo)了求解二維瞬態(tài)熱傳導(dǎo)的第一類(lèi)邊值問(wèn)題的基本過(guò)程。由上述兩個(gè)算例得出,本文提出的瞬態(tài)熱傳導(dǎo)問(wèn)題的插值型復(fù)變量EFG方法具有較好的擬合效果。