纯无网格并行计算在传热方程数值模拟中的应用

胥 康,任恒飞,任金莲,蒋 涛

(扬州大学 数学科学学院,江苏 扬州 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) 并行算法能可靠地求解无解析解的热传导方程.