純無網(wǎng)格并行計算在傳熱方程數(shù)值模擬中的應(yīng)用

胥 康,任恒飛,任金蓮,蔣 濤

(揚州大學(xué) 數(shù)學(xué)科學(xué)學(xué)院,江蘇 揚州 225002)

傳熱方程[1]是偏微分方程[2]的一個重要分支,目前求解傳熱方程的數(shù)值解法很多,有限點集法(FPM)[3]是其中重要的數(shù)值解法之一.通常這類問題的計算量很大,需要數(shù)億次的計算,因此如何提高計算效率,縮短求解時間,成為研究者們急需解決的問題.而高性能并行計算機的出現(xiàn)極大提高了大規(guī)模計算問題的計算效率,因此將并行計算運用到傳熱方程數(shù)值模擬中有一定的意義.目前,并行計算技術(shù)有很多,主要有MPI(massage passing interface)[4-5],OPENMP,CUDA,OPENGL,其中MPI是一種比較成熟高效的并行技術(shù).本文在FPM方法的基礎(chǔ)上,通過引入MPI并行計算,對三維傳熱方程進行數(shù)值求解,并分析并行效率,從而驗證研究傳熱方程時施加并行算法的必要性和重要性.

1 問題描述

傳熱方程的一般形式為

(1)

初始條件為

u(x,y,z,0)=φ(x,y,z), (x,y,z)∈Ω,

(2)

邊界條件為

u(x,y,z,t)=φ(x,y,z), (x,y,z,t)∈?Ω×(0,T],

(3)

其中:ki=ki(x)(i=1,2,…,n)為熱傳導(dǎo)系數(shù);函數(shù)u=u(x,t)是固體在熱傳導(dǎo)過程中t時刻、x處的溫度;Ω為求解區(qū)域.

2 基于FPM方法的并行算法

2.1 FPM方法

有限點集法屬于無網(wǎng)格方法[6],其思想是確定待求點的支持域,將支持域內(nèi)的每個點通過Taylor展開到三階導(dǎo)數(shù)得到關(guān)于導(dǎo)數(shù)的方程,再用最小二乘法使加權(quán)誤差最小,求得待求點處的各階導(dǎo)數(shù),最后迭代求出該點處的數(shù)值.設(shè)xi為點x附近的點(i=1,2,…,n.),函數(shù)u(x,t)，ui(t)表示u(x,t)在xi處、t時刻的函數(shù)值,則u(xi,t)在x點的三階Taylor展開式為

其中:ei為誤差;xi1,xi2,xi3是點xi的x,y,z分量;x1,x2,x3是點x的x,y,z分量;導(dǎo)數(shù)uk,ukl和uklj(k,l,j=1,2,3)可以通過最小二乘法求出.上述問題可寫成

en×1=Mn×19a19×1-bn×1,

(4)

其中

M中第i行為

Δxki=xik-xk, Δxkli=(xik-xk)(xil-xl),

Δxklji=(xik-xk)(xil-xl)(xij-xj),

a19×1=(u1,u2,u3,u11,u12,u13,u22,u23,u33,u111,u112,u113,u122,u123,u133,u222,u223,u233,u333)T,

bn×1=(u1-u,u2-u,u3-u,…,un-u)T,en×1=(e1,e2,…,en)T.

函數(shù)ω為

α為常數(shù),且α>0,取α=6.25.h決定x的支持域,即x為中心、h為半徑的一個球,記為p(x,h)={xi;i=1,2,…,n}.易推出

a=(MTWM)-1(MTW)b,

(5)

求出相應(yīng)的導(dǎo)數(shù)即可得到下一時間層的函數(shù)值.

2.2 MPI并行

在進行MPI計算時,使用若干個CPU以加快計算效率,這若干個CPU會運行一段相同的代碼.由于每個進程都有自己的進程號,因此可通過這些進程號決定不同進程執(zhí)行不同行為.

本文涉及的FPM算法,需要先確定支持域內(nèi)的相鄰粒子,這一步消耗的時間較多.為提高相鄰粒子搜索的計算效率,需考慮粒子搜索并行,即考慮將所有粒子分配在多個CPU上,同時進行相鄰粒子搜索并標(biāo)記.此外,粒子物理量的循環(huán)求解也需要實現(xiàn)并行,同樣將粒子分配給多個CPU同時進行求解,以提高計算效率.先后兩次并行為CPU分配粒子數(shù)相同,只需分配一次.因此,本文基于FPM方法的MPI并行算法主要體現(xiàn)在相鄰粒子搜索標(biāo)記和循環(huán)求解過程中,采用多個CPU計算以提高計算效率.

3 數(shù)值算例

3.1 有解析解的三維算例

例1為了分析該并行算法的并行效率及可靠性,先對有解析解算例進行數(shù)值模擬.考慮求解區(qū)域Ω: [0,1]×[0,1]×[0,1]中的常系數(shù)非穩(wěn)態(tài)傳熱問題[7],其方程為

ut=κ(uxx+uyy+uzz),

初值條件為

u(x,y,z,0)=sin(πx)+sin(πy)+sin(πz),

邊值條件為

本文參數(shù)κ=0.1,對應(yīng)該問題的解析解為

u(x,y,z,t)=[sin(πx)+sin(πy)+sin(πz)]e-κπ2t.

Dirichlet邊界條件易處理,可直接賦值.為體現(xiàn)數(shù)值模擬的準(zhǔn)確性,先取粒子數(shù)為61×61×61,時間步長為dt=10-4,CPU為24,圖1為數(shù)值模擬結(jié)果與解析解的對比曲線.由圖1可見,幾個不同時刻的數(shù)值結(jié)果均與解析解相符,表明該并行算法可靠.再取不同粒子數(shù),將數(shù)值解與解析解進行比較,得到最大誤差范數(shù)L∞,結(jié)果列于表1.

圖1 幾個不同時刻、不同位置處沿z方向的變化曲線Fig.1 Variation curves along z direction at several different times and locations

粒子數(shù)61×61×6181×81×81101×101×101誤差L∞0.000 0930.000 0750.000 049

由表1可見：

1) 最大誤差值隨著粒子數(shù)增加而減小；

2) 本文數(shù)值方法模擬常系數(shù)非穩(wěn)態(tài)傳熱問題時接近2.5階精度(由表1數(shù)據(jù)估計得到),進一步體現(xiàn)了本文算法的精確性.

為了考察并行運算對計算效率的影響,計算粒子數(shù)為61×61×61,時間步長為dt=10-4,運算到0.5 s時不同CPU數(shù)下的總消耗時間,結(jié)果列于表2.由表2可見,采用本文算法求解三維傳熱方程的計算量很大,因此考慮并行計算是非常必要的.

表2 粒子數(shù)為61×61×61時不同CPU數(shù)下運算到0.5 s時的消耗時間(s)

為了更好地體現(xiàn)并行計算的效率,表3列出了不同粒子數(shù)、不同CPU數(shù)下第一步(包含粒子搜索)所需的時間.由表3可見：當(dāng)CPU數(shù)不變時,計算時間隨著粒子數(shù)的增加而增加；當(dāng)粒子數(shù)不變時,計算效率隨著CPU數(shù)的增加而得到提高.表4列出了不同粒子數(shù)、不同CPU數(shù)下計算到1 s(除第一步)的平均消耗時間.由表4可見：當(dāng)CPU數(shù)不變時,計算時間隨著粒子數(shù)的增加而增加,且計算量增加比率與粒子數(shù)增加比率并不成線性正比關(guān)系；當(dāng)粒子數(shù)不變時,計算效率隨著CPU數(shù)的增加而得到提高,但計算效率的提高比率與CPU數(shù)增加比率也不成線性正比關(guān)系.這是因為在數(shù)值模擬過程中,CPU的計算時間受編程語言、網(wǎng)絡(luò)通信環(huán)境及高性能設(shè)備等因素的影響.

通過例1及本文方法與FDM(有限差分)方法[8]的構(gòu)造過程發(fā)現(xiàn),本文方法不僅能精確可靠地模擬規(guī)則區(qū)域下的三維傳熱問題,較FDM方法還具有如下優(yōu)點：

1) 程序?qū)崿F(xiàn)相對簡單,特別對復(fù)雜區(qū)域,如求解圓柱形區(qū)域,FDM方法在程序上很難實現(xiàn)；

2) 涉及線性方程組的計算時,FPM方法是局部的系數(shù)矩陣,FDM方法涉及的系數(shù)矩陣明顯大很多；

3) FPM方法容易推廣應(yīng)用到非規(guī)則區(qū)域問題上的離散.

表3 不同粒子數(shù)、不同CPU數(shù)下第一步的消耗時間(s)

表4 不同粒子數(shù)、不同CPU數(shù)下(除第一步)平均每步的消耗時間(s)

例2為體現(xiàn)本文方法在模擬非矩形復(fù)雜區(qū)域上溫度傳播問題時較FDM方法的優(yōu)勢,考慮圓柱形區(qū)域且采用圓形粒子分布方式.

圓柱形區(qū)域上帶混合邊界的瞬態(tài)傳熱方程[9]為

初值條件為

u(x,y,z,0)=0,

Dirichlet邊界條件為

u(r,t)=100,r=1(r為極坐標(biāo)),

Neumann邊界條件[6,10-11]為

u,z|z=0=u,z|z=2=0.

參數(shù)k/c=5.該算例的邊界條件是混合邊界條件,Dirichlet邊界條件可直接賦值,對于Neumann邊界條件:

可采用文獻[6]的處理方法.邊界點x處需添加一個方程:

0=u1(x,t)nx+u2(x,t)ny+u3(x,t)nz,

矩陣M和W相應(yīng)的要增加一行,其中nx,ny,nz為在邊界點x處單位法向量n的x,y,z分量,u1(x,t),u2(x,t),u3(x,t)函數(shù)關(guān)于x,y,z的偏導(dǎo)數(shù).

圖2(A)為三維圓柱區(qū)域及粒子的分布情況;圖2(B)為沿z=1處截面極坐標(biāo)方向上溫度的變化曲線.由圖2及FDM和FPM方法構(gòu)造過程可見：本文FPM-3D方法較FDM法容易求解非矩形區(qū)域熱傳導(dǎo)問題,且本文方法得到的結(jié)果與解析解相符；給出的粒子方法易實現(xiàn)帶混合邊界復(fù)雜區(qū)域傳熱問題的模擬,且計算結(jié)果可靠.

圖2 三維圓柱區(qū)域內(nèi)的粒子分布情況(A)及圓柱形區(qū)域下沿z=1處截面極坐標(biāo)方向上溫度的變化曲線(B)Fig.2 Particle distribution in three-dimensional cylindrical region (A) and variation curves of temperature along polar coordinate direction of cross-section at z=1 in cylindrical region (B)

3.2 無解析解的三維算例

為進一步驗證本文并行算法的可靠性,下面對無解析解算例進行數(shù)值模擬,并與FDM方法求得的數(shù)值結(jié)果做對比.考慮求解區(qū)域為Ω: [0,1]×[0,1]×[0,1],帶有與時間有關(guān)的混合邊值條件的變系數(shù)瞬態(tài)傳熱方程[9]:

c(x,y,z)ut=(k(x,y,z)u),

初值條件為

u(x,y,z,0)=0,

Dirichlet邊界條件為

u(x,y,1,t)=10t,

Neumann邊界條件為

u,x|x=0=u,x|x=1=u,y|y=0=u,y|y=1=u,z|z=0=0.

為方便與文獻[9]中的數(shù)值結(jié)果做對比,選取c(x,y,z)=1e3z,k(x,y,z)=5e3z,對應(yīng)的k/c=5.取粒子數(shù)71×71×71,時間步長為dt=10-5,CPU數(shù)為24,計算到1 s,結(jié)果如圖3和圖4所示.圖3為3個不同時刻x=y=0.5截面上溫度沿z軸的變化曲線.由圖3可見,該并行算法模擬混合邊界條件變系數(shù)下瞬態(tài)傳熱方程是穩(wěn)定可靠的.圖4為三維功能材料上的溫布分布.

圖3 不同時刻溫度沿z軸變化的曲線(x=y=0.5截面上)Fig.3 Variation curves of temperature along z axis at different times (x=y=0.5 cross section)

圖4 不同時刻三維功能材料上的溫度分布Fig.4 Temperature distribution on three-dimensional functional materials at different times

綜上所述,本文采用有限點集法的并行算法對熱傳導(dǎo)問題進行了求解,通過對有解析解傳熱問題的模擬,分析了并行計算的計算效率和可靠性,并把該并行算法用于求解變系數(shù)瞬態(tài)熱傳導(dǎo)方程中,可得以下結(jié)論：

1) 當(dāng)CPU不變時,計算時間隨著粒子數(shù)的增加而增加,且計算量增加比率與粒子數(shù)增加比率并不成線性正比關(guān)系；

2) 當(dāng)粒子數(shù)不變時,計算效率隨著CPU數(shù)的增加而得到了提高,但計算效率的提高比率與CPU數(shù)增加比率也不成線性正比關(guān)系；

3) 并行算法能可靠地求解無解析解的熱傳導(dǎo)方程.