首页 理论教育数值方程算法与验证

数值方程算法与验证

【摘要】:Stone对全隐式差分方程提出了新的迭代解法——强隐式迭代解法,简称SIP方法。由于这种迭代方法利用解方程组对全部格点同时加以改进,迭代效率非常高,已经十分接近直接解法,因而也就更加有效。Trescott等所做的比较研究表明:对于复杂的实际问题,SIP方法明显优于其他的迭代方法。表5-3-1ADI算法时间为500天时浓度场计算结果表5-3-2SIP算法时间为500天时浓度场计算结果图5-3-6SIP算法实例验证图算例2作者对三维SIP算法也进行了编程验证,作为

Stone(1968)对全隐式差分方程提出了新的迭代解法——强隐式迭代解法,简称SIP方法。由于这种迭代方法利用解方程组对全部格点同时加以改进,迭代效率非常高,已经十分接近直接解法,因而也就更加有效。Trescott等(1977)所做的比较研究表明:对于复杂的实际问题,SIP方法明显优于其他的迭代方法。

1.强隐式差分格式

将一般的全隐式差分方程写为以下的形式(二维):

这里略去水头关于时间的上标n+1,与式(5-3-16)中每个系数相联系的各点位置见图5-3-2,共五个格点,称式(5-3-16)为五点格式。

对渗流区域内的每个格点都列出式(5-3-16)并利用给定的边界条件,就得到一线性方程组

暂时设渗流区域为一矩形,划分成N×M个差分网格(图5-3-3),格点号码按自左而右、由下而上的顺序编排,则由五点格式所确定的上式中的系数矩阵[K]具有下列形式:

图5-3-3 假设的矩形2D渗流区位置图

矩阵的第一行相当于对格点(1,1)列的方程,因为利用了边界条件,所以a1,1和c1,1都不出现,其他依此类推。可以看出,矩阵[K]的非零元素全部集中在五条对角线上。此外,它还是一个对称的矩形。矩阵[K]具有式(5-3-18)的形式是由格点号码的特殊编排方式所决定的。对于这样一个形状简单、高度稀疏的矩阵,在对它进行LU分解时,期望能分解成下列形式的[L]和[U]的乘积:

下三角矩阵[L]和上三角矩阵[U]包含的非零元素的位置类似把[K]沿主对角线分成两半。假若能做这样的分解,那么从[K]得到[L]和[U]的计算量一定是很小,因为[L]和[U]也是高度稀疏的。事实上对[K]分解得到[L]和[U]并不稀疏,使得计算量变得很大。

于是把问题倒过来,由式(5-3-19)和式(5-3-20)确定[L]和[U]的乘积,按照矩阵乘法的规则,有

式(5-3-21)的非零元素集中在七条对角线上,比[K]多了两条非零元素对角线。假若对每个格点列的差分方程不用五点格式而是在形式上用七点格式,如式(5-3-22):

与公式(5-3-22)的系数相关联的格点表示在图5-3-4中。由此形成的方程组设为

图5-3-4式(5-3-22)系数相关联的格点示意图

图5-3-4 式(5-3-22)系数相关联的格点示意图

其系数矩阵为

2.系数矩阵的近似因子分解

比较式(5-3-21)和式(5-3-24),可直接得到[^K]的元素[L]、[U]的元素之间的下列对应关系

另一方面,根据“七点”格式与“五点”格式之间的关系可以找到[^K]与[K]的关系。事实上,利用Taylor展开式并略去高阶项得

因为

所以由式(5-3-27)可得下面的近似关系

将上式代入到式(5-3-16)中,得到

这样就把式(5-3-22)近似地化成了式(5-3-16)的形式。与式(5-3-16)比较得

式中:ω为松弛因子,目的是为了加快收敛速度。

在上式前添加松弛因子ω(0<ω<1)。

计算是按照从左而右、自下而上的顺序进行,则式(5-3-31)中的δi,j-1、δi-1,j、ηi-1,j、ηi,j-1应作为已知,因为它们已经在左、下格的计算中被求出或者可根据边界条件定出。因此,借助式(5-3-32),对每个格点(i,j)都能由“五点”格式系数ai,j、bi,j、ci,j、di,j、ei,j算出[L]和[U]的非零元素αi,j、βi,j、γi,j、ηi,j和δi,j。矩阵[L]和[U]是[^K]的因子分解,不是[K]的因子分解,但[]和[K]很近似。到这里已经实现了对系数矩阵[K]的近似因子分解。

松弛因子ω的确定通过以下公式:

其中

SIP迭代对ω的值并不敏感,ω一般大于0.92,在[0.92,0.95]之间(Peric,1995)。

3.强隐式迭代解法

解方程组[K]H=F等价于解方程组

这一方程组左端的系数矩阵[^K]已有LU分解,且[L]和[U]均为稀疏矩阵,因此很容易求解,但因式(5-3-34)的右端也含有未知的水头H,所以只能用迭代解法,将式(5-3-34)改写为下列形式

式中 H——右上角的指数表示迭代的次数,两端同减Hm,得

或者

其中

迭代的过程是:已知Hm,可由式(5-3-39)决定Rm,以此作为式(5-3-40)的右端项,解方程组求得ξm+1,再由式(5-3-38)解出Hm+1,如此重复迭代直到收敛为止。收敛的准则采用下式

其中ε是足够小的正数。

上述迭代过程主要工作是解方程式(5-3-37),[^K]已有LU分解,所以在解方程时很容易求解,只需一个前推回代的过程。

4.前推回代工程

式(5-3-43)和式(5-3-44)为中间值V和ξm+1的表示方法,这样就能顺利的解出[L]和[U],从而求解整个方程。

实际含水层区域不一定是矩形,但一定可以采用矩形区域将其覆盖,用两组平行线划分成矩形网格,在这些矩形网格中可分成内格点、靠近边界的格点和外格点三类。这样就可以运用SIP方法求解。

5.SIP算法验证

算例1

在200m×200m的区域上(浓度为0),北部和西部隔水,东部和南部边界有稳定的浓度为1mg/L,且地下水为稳定流动,孔隙平均流速为0.005m/d,流动方向平行于X轴并与X方向相反。ΔX=ΔY=20m,Δt=5d。总时段为2000d。共400个步长,DL=10,DT=1(采用陈崇希、李国敏等编著的《地下水溶质运移理论及模型》“求矩形含水层中地下水溶质运移二维ADI交替隐式差分算例”)。

区域剖分分辨率为20m×20m,共100个网格。剖分网格与初始浓度场图见图5-3-5。

图5-3-5 有限差分网格剖分图及初始浓度场图

表5-3-1和表5-3-2是时间为500天时分别用ADI算法和SIP算法算出的区域浓度场分布结果。可以看出两结果稍有不同,但是仅在小数点后第三位开始略有不同,并不影响计算精度。这些差别也主要是计算方法不同引起的。图5-3-6为图化的输出结果。

表5-3-1 ADI算法时间为500天时浓度场计算结果

表5-3-2 SIP算法时间为500天时浓度场计算结果

图5-3-6 SIP算法实例验证图

算例2

作者对三维SIP算法也进行了编程验证,作为变密度三维渗流溶质运移模型开发的基础条件,算例如下。利用SIP算法计算三维拉普拉斯方程:将边长为1的立方体剖分成20×20×20=8000的网格,Δx=Δy=Δz=0.05;网格的初始值均为0;立方体的下边、左边和后边的边界条件均为零通量边界,即U(0,Y,Z)=U(X,0,Z)=U(X,Y,0)=0;立方体的上边、前边和右边的边界条件分别为U(X,Y,1)=X×Y,U(1,Y,Z)=Y×Z,U(X,1,Z)=X×Z,(来自J.H.Ferziger和M.Peric的著作《Computational Methods for Fluid Dynamics》,2002)

方程离散格式如下:

利用SIP迭代求解结果如图5-3-7和图5-3-8所示。

图5-3-7 SIP算法验证3D切片图(从左到右从上到下z为1~21)

注:横坐标为1~21,间隔为1;纵坐标为0~1,间隔为0.1

图5-3-7中从左到右从上到下21张图分别表示z=1,21时的切片图。图5-3-8表示的是其三维立体图[顶角格点坐标为U(20,20,20)],两图能清楚地反映出算法的有效性,即在立方体的最上端边长x=1、y=1、z=1时数值最大,其余内部格点值按照离最上端远近距离均匀分布。

图5-3-8 SIP算法3D三维验证图