扩散方程九点格式中节点未知量的一种新的插值算法

董 成,邬吉明

(1.中国工程物理研究院研究生院,北京 100088)

(2.北京应用物理与计算数学研究所,北京 100088)

1 引言

考虑扩散方程边值问题

其中Ω⊂R2是平面上的多边形区域,Λ是扩散张量,f是源项,∂Ω=¯ΓD∪¯ΓN是Ω的边界,n为单位外法向量,gD和gN分别是Dirichlet和Neumann边值.

上述扩散问题有着广泛的实际应用背景,例如在辐射流体力学、辐射磁流体力学、油藏模拟等应用中,扩散过程与其他物理过程相耦合.我们常常需要在网格强扭曲、强间断、强各项异性、强非线性等极端条件下求解扩散方程,这是一件具有挑战性的工作.有限体积方法是求解扩散问题最常用的方法之一,它具有局部守恒、简单容易实现等良好性质.多年来,许多科研工作者致力于扩散方程的有限体积方法研究,提出了众多的离散格式,如九点格式(NPS)[1]、多点通量逼近格式(MPFA)[2]、支撑算子格式[3]、广义差分格式[4]等.按照未知量的类型,这些格式可大致分为单元中心格式、节点格式、杂交格式、混合格式等等,详细的最新研究进展可参见文献[5–9].

我国学者李德元教授等人在二维网格上基于积分插值法构造了一个单元中心型有限体积格式[1],由于该格式在结构四边形网格上有九点模板,故常称其为九点格式,它是若干辐射流体力学程序的基本格式[10,11].由于单元中心格式在一个单元上只有一个未知量,在构造离散格式时需要引入辅助未知量来提高精度.九点格式的辅助未知量定义在网格节点处,如何用单元中心未知量逼近节点辅助未知量是九点格式研究中的一个重要内容.一个理想的二维节点插值算法应当具有简单、保正、不依赖于网格拓扑、不依赖于间断线的位置、有二阶精度、易于向三维推广等性质.目前已知的节点插值算法,如算术平均插值[1]、逆距离加权插值[12]、最小二乘插值[13]、显式加权算法[14]等,均只能满足部分要求.由于缺乏容易向三维推广的二阶插值算法,严重制约了九点格式在三维问题中的应用.

本文在多点通量逼近的边辅助未知量插值算法的基础上,通过一个特殊的极限技巧获得了一个新的节点插值算法,并在给定假设下证明了该算法中局部线性系统的可解性.本文插值算法的一个重要之处在于容易向三维推广,且满足除保正性以外的其他要求.全文的推导基于所谓的线性精确准则[15],即当扩散张量关于网格是分片常数、解析解关于网格是分片线性函数时,算法推导的每一步都是精确成立的.为了和通常的等式相区别,用符号’表示相关等式是在线性精确意义下成立的.

2 九点格式的构造

最初的九点格式是在扩散系数为标量的情形下推导的,但我们容易将其推广到扩散系数为张量的情形[14],为保持本文内容的完整性,我们简要介绍一下扩散张量情形下九点格式的一种构造算法．首先将Ω剖分为互不重叠的多边形网格,K表示其中的一个单元,σ表示单元边,K 的所有边组成的集合记为EK,单元K 关于边σ的单位外法向量记为nK,σ.记流向量为F=−Λ∇u,这里假定Λ关于网格是分片常数,并用ΛK表示Λ在单元K 上的限制.单元K 通过边σ的流记为FK,σ,即

在单元K上对方程(1.1)两端积分并利用散度定理,有

于是有限体积法解扩散问题就归结为构造流FK,σ的离散表达式.

在图1中,K和L表示两个相邻单元,它们的单元中心分别是OK和OL,σ是它们的公共边,其顶点为A 和B,nK,σ,nL,σ分别是K,L关于σ的单位外法向量.dK,σ,dL,σ分别表示OK,OL到σ的距离,uA,uB,uK,uL表示u在点A,B,OK,OL处的近似值.对于向量v=(a,b)T,记v◦=(b,−a)T.此外,为叙述简单起见,引入以下记号(见图1)

首先,有如下的向量分解

其中

图1:离散流构造模板

|σ|表示σ的长度.其次,在单元K 上,根据文献[14]的引理2.1,

注意到 |σ|nK,σ=t◦,t=aK− bK,·t= −|σ|dK,σ,将 (2.4) 和 (2.5) 式代入 (2.1) 式,经过简单的计算后就可以得到

类似地,在单元L上,有

通常要求通过边σ的流满足连续条件,即

由此得到

将(2.6)和(2.7)式代入(2.9)式的右端,利用aK−bK=aL−bL=t,aK−aL=bK−bL=s,并经过简单化简就得到内部边σ上流的离散表达式

其中

容易看出,上述离散流涉及单元中心未知量uK和uL以及节点辅助未知量uA和uB,要得到纯单元中心格式,需要用单元中心量对节点辅助未知量进行插值.

3 节点辅助未知量的插值算法

如图2(a),假设Q0是一内部节点,它周围的单元Ωk(1≤k≤n)按逆时针顺序排列.Q0Pk和Q0Pk+1是Ωk的以Q0为顶点的边,Ok为其中心.假设Tk是边Q0Pk上的一点,由下式确定

其中τ∈(0,1).令u0,uk,分别表示Q0,Ok,Tk处的未知量,Λk表示Λ在Ωk上所取的分片常数值,nk表示Ωk关于Q0Pk的单位外法向,见图2(b).本文中k将被视为以n为周期的指标,例如当k=n+1和0时,有Pn+1=P1,T0=Tn.

图2:内部节点Q0处的局部结构和记号

3.1 边辅助未知量的插值算法

在多点通量逼近格式中,边未知量¯uk可以通过单元中心未知量uk表示.为叙述简洁起见,引入以下记号

其中ek表示n阶单位矩阵的第k个列向量.显然,Tk是一个n×2的矩阵,并且

在单元Ωk的边Q0Pk和Q0Pk+1上,定义如下离散流

注意到对于 Q0Pk(Q0Pk+1),有于是

其中

这里记号(τ)表示相关的量是τ的函数,在不引起歧义的情况下将其省略.利用(3.3)式,可将(3.4)式改写为

再利用边Q0Pk+1上的流连续条件,即

或者其等价形式

就得到

其中

通过求解(3.8)式,就能够用单元中心未知量来表示边未知量,即=M−1NU.容易看出,只有当M可逆时,上述边未知量插值算法才是可行的.自多点通量逼近算法提出来的二十多年里,人们还没有在数值计算中遇到过M奇异的情形,但这一点至今仍缺乏严格的理论证明.

3.2 新的节点未知量插值算法

其中Sk,1(Sk,2)表示4Q0PkOk(4Q0OkPk+1)的面积.于是(3.5)式可以写为

其中

再利用(3.9)式,可得

其中

需要注意,由于(3.11)式的分母中含有τ,Ml和Nl依然是τ的函数.根据(3.11)式,发现

故有

其中1和0分别表示所有分量为1和0的向量,O表示零矩阵.利用(3.12)和(3.14)式,可证(3.8)式等价于

注意到解沿着Q0Pk(1≤k≤n)是切向连续可微的,通过泰勒展开,可以得到

其中 Γ =(γ1,···,γn)T是与 τ无关的常向量,O(τ2)表示所有分量为 O(τ2)的向量. 将(3.16)式代入(3.15)式并且利用(3.14)式中的第一个等式,得到

最后将以上方程两端同时除以τ,并令τ→0,就得到最终的方程

其中M1(0),M2(0)和N2(0)表示M1,M2和N2取τ=0的值.将M1(0)的第一列替换为M2(0)1其余列不变得到一个新的矩阵,记为,则(3.17)式可以写为与τ无关的线性系统

求解上述局部线性系统就得到一个新的节点插值算法,即

这里需要特别强调,虽然上述算法基于多点通量逼近的边辅助未知量插值,但它与(3.8)式中矩阵M的可逆性没有关系.

4 局部线性系统(3.18)的可解性

虽然(3.8)式的可解性目前尚无理论证明,但局部线性系统(3.18)式的可解性在一定条件下是可以严格证明的.本节的主要结果依赖以下两个假设:

在扩散系数是标量的情况下,假设(H1)和(H2)的几何意义是明显的,此时,两个假设变为

图3:假设(H1)和(H2)的几何意义

引理4.1在假设(H1)下,M2(0)1的所有元素都是负的,即

证 由M2(0)的定义,有

其中δij表示Kronecker delta函数.利用(3.11)式并且注意到τ=0,进一步得到

再结合假设(H1),立即得到(4.1)式.

引理4.2记A=(aij)为n×n的矩阵,满足

其中ak,bk(1≤k≤n)是正数.则

证 根据(4.2)式,有

分裂An−1的最后一列,可得

对于An−1的第一部分,通过将最后一列加到第(n−2)列,再将新的第(n−2)列加到第(n−3)列,依此类推,最终可得到一个上三角型行列式,进而有

用同样的方法可以得到

注意到A1=a1+bn>0以及ai,bi>0,通过数学归纳法最终得到An−1>0,从而>0.类似地,可以证明>0.当i 6=1,n时,记aii对应的余子式为,并令

其中Ik是k阶单位矩阵.容易验证是形如(4.4)式的三对角矩阵.于是

定理4.3 在假设(H1)和(H2)的条件下,(3.18)式中的系数矩阵 eM是非奇异的.

证 由M1(0)的定义和(3.13)式,有

由(3.11)式和假设(H2),可以得到

令A=(aij)n×n= −M1(0).根据引理4.2,>0(1≤k≤n).由 eM 的定义,并利用引理4.1,可得

5 数值实验

首先,定义如下的L2误差和流误差

其中M为单元集合,eK表示单元K中心的误差,Nσ表示边σ的相邻单元数,L表示σ的除K 之外的另一相邻单元（如果存在的话）,当σ在边界上时,我们约定dL,σ=eL=0.两个网格层之间的收敛率Rα(α=u,q)通过以下公式获得

其中h1,h2代表两个网格层的网格尺寸,Eα(h1),Eα(h2)为对应的离散误差.

用BICGSTAB方法求解离散线性系统,并取收敛准则为εlin=1.0E-15.离散流采用公式(2.10),节点插值算法分别采用文献[14]中的显式加权算法2和本文的极限加权算法,对应的格式简记为NPS-EW2和NPS-LW.采用的计算网格见图4,且所有的计算均在双精度下完成.

图4:计算网格

5.1 实验1:各向异性问题

考虑Ω=[0,1]2上的具有全Dirichlet边界的扩散问题(1.1).扩散张量和精确解如下

这个数值实验是文献[5]中的一个测试的一个简单变形.

从实验结果来看,NPS-LW在各网格上与NPS-EW2的实验结果接近,这里仅给出前者的计算结果.表1给出了NPS-LW在各个网格上数值解和流的误差,图5以对数折线图的形式展示了收敛速度.可以看出NPS-LW在这些网格上的数值解的误差随着网格加密以趋近2阶的速度收敛.

表1:NPS-LW在实验1中的数值结果

图5:NPS-LW在实验1中的误差曲线

5.2 实验2:泊松问题

在[0,1]2上求解泊松方程−∆u=4.0,采用Dirichlet边界条件,解析解取为

本算例十分简单,其设计目的是测试新的节点加权算法与已有显式加权算法在数值表现上的差异.NPS-LW的数值结果在Mesh1,Mesh2和Mesh4上与NPS-EW2非常接近,但在Mesh3上明显优于NPS-EW2,在表2中进行了对比.

表2:NPS-LW与NPS-EW2在Mesh3上的数值结果对比

5.3 实验3:间断系数问题

在Ω=[0,1]2上求解扩散问题(1.1)–(1.2),扩散张量取为

I为二阶单位矩阵,精确解取为

这里采用Mesh1和Mesh4.对于Mesh1,所有位于x=0.5上的节点在x方向上都不扰动.数值计算结果见表3.从该表中可以看到,数值解误差收敛速度均趋于2阶,流的误差接近1阶,均为最优阶,表明新的节点加权算法也能很好地适应间断系数问题.

表3:实验3 NPS-LW在Mesh1,4上的数值结果