純無網格并行計算在傳熱方程數值模擬中的應用

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

(揚州大學 數學科學學院,江蘇 揚州 225002)

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

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)為熱傳導系數;函數u=u(x,t)是固體在熱傳導過程中t時刻、x處的溫度;Ω為求解區域.

2 基于FPM方法的并行算法

2.1 FPM方法

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

其中:ei為誤差;xi1,xi2,xi3是點xi的x,y,z分量;x1,x2,x3是點x的x,y,z分量;導數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.

函數ω為

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

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

(5)

求出相應的導數即可得到下一時間層的函數值.

2.2 MPI并行

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

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

3 數值算例

3.1 有解析解的三維算例

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

ut=κ(uxx+uyy+uzz),

初值條件為

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

邊值條件為

本文參數κ=0.1,對應該問題的解析解為

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

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

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

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

由表1可見：

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

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

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

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

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

通過例1及本文方法與FDM(有限差分)方法[8]的構造過程發現,本文方法不僅能精確可靠地模擬規則區域下的三維傳熱問題,較FDM方法還具有如下優點：

1) 程序實現相對簡單,特別對復雜區域,如求解圓柱形區域,FDM方法在程序上很難實現；

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

3) FPM方法容易推廣應用到非規則區域問題上的離散.

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

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

例2為體現本文方法在模擬非矩形復雜區域上溫度傳播問題時較FDM方法的優勢,考慮圓柱形區域且采用圓形粒子分布方式.

圓柱形區域上帶混合邊界的瞬態傳熱方程[9]為

初值條件為

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

Dirichlet邊界條件為

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

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

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

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

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

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

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

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

圖2 三維圓柱區域內的粒子分布情況(A)及圓柱形區域下沿z=1處截面極坐標方向上溫度的變化曲線(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 無解析解的三維算例

為進一步驗證本文并行算法的可靠性,下面對無解析解算例進行數值模擬,并與FDM方法求得的數值結果做對比.考慮求解區域為Ω: [0,1]×[0,1]×[0,1],帶有與時間有關的混合邊值條件的變系數瞬態傳熱方程[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]中的數值結果做對比,選取c(x,y,z)=1e3z,k(x,y,z)=5e3z,對應的k/c=5.取粒子數71×71×71,時間步長為dt=10-5,CPU數為24,計算到1 s,結果如圖3和圖4所示.圖3為3個不同時刻x=y=0.5截面上溫度沿z軸的變化曲線.由圖3可見,該并行算法模擬混合邊界條件變系數下瞬態傳熱方程是穩定可靠的.圖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

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

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

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

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