多尺度有限元法在分层网格上的边界层数值模拟

孙美玲, 刘 雪, 江 山*

(1. 南通大学理学院, 江苏 南通 226019; 2. 南通职业大学数学教研室, 江苏 南通 226007)

奇异摄动问题广泛存在于物理工程、化学反应、流体力学和生态种群等领域[1]. 近年来, 针对摄动参数ε很小时存在的边界层现象,主要采取基于有限元法(finite element method, FEM)的数值方法求解. 例如, Duran等[2]基于分层网格采用经典有限元法, Zhu等[3]基于Shishkin网格采用高阶间断有限元法, 杨宇博等[4]基于分层网格采用非对称内罚间断有限元法, Chen等[5]利用奇异值分解采用间断有限元法, 分别对一维对流扩散方程进行求解; Kumar[6]针对两端出现边界层的扩散模型, 应用五次B-样条配点法获得一致稳定性结果; Brdar等[7]在分层网格上采用经典有限元法求解对流项和扩散项含双参数的二维摄动方程, 获得了较精确的数值结果; Teofanov等[8]采用流线扩散有限元法求解双摄动参数ε1,ε2的扩散模型. 然而, 上述工作均在经典有限元格式框架下进行, 计算代价偏高. 多尺度有限元法(multiscale finite element method, MsFEM)在经典有限元理论和格式的基础上通过多尺度基函数将微观性态代入宏观刻画, 不仅能节约计算资源, 而且可保证计算精度, 为有效求解奇异摄动边界层问题提供了新的思路和技术[9-11]. Song等[10]提出组合有限元和多尺度有限元的计算格式, 高效模拟了带奇异性的多尺度摄动问题; Kopteva[12]给出了奇异摄动对流扩散方程基于各向异性网格的能量范数误差估计; Cheng[13]使用局部间断有限元法求解非线性方程. 本课题组在前期工作[14]中利用多尺度基函数在Shishkin网格上求解反应扩散方程, 由于Shishkin网格依赖摄动参数ε和事先选定剖分数N来估计过渡点, 难以精确反映边界层的确切位置,故数值精度和收敛速度不够.因分层网格可依靠摄动系数ε和选取网格参数0

1 问题的提出

考虑二维对流扩散方程[11-12]

2 自适应分层网格的离散

通常情况下,若对求解区域Ω进行一致剖分, 即在单维x或y方向各取剖分数N, 当其取值充分大时, 二维区域[0,1]2的等距步长h=1/N→0.然而,当摄动系数ε很小时难以满足h<ε, 无法有效处理小参数奇异摄动问题.

分层网格以x方向的节点xi为例,

将二维区域离散化以自适应地逼近边界层位置.其剖分数N须满足xN-1<1, (1+h)xN-1≥1, 并确保执行时最后一个单元[xN-1, 1]的步长不小于前一个单元[xN-2,xN-1]的步长, 否则去除节点xN-1.记x方向第i单元[xi-1,xi]的步长hi=xi-xi-1.类似地, 记y方向第j单元[yj-1,yj]的步长hj=xj-xj-1.

分层网格能够自动捕捉边界层的准确位置, 其网格剖分数N并非简单地呈倍数增加,在加密过程中直接由迭代公式(3)随机生成所需的剖分数N.

3 本文方法

3.1 经典有限元法

a(ug,v)=(f,v), ∀v∈Vh.

(4)

Afug=bf,

(5)

其中Af,bf分别为总刚度矩阵和总载荷向量, 从而解得有限元解ug.

3.2 多尺度有限元法

基于分层网格(3), 离散区域在x,y方向均已形成粗网格剖分数NMsFEM, 在每个粗网格单元K求解齐次子问题

其中φi为多尺度有限元的基函数.给定边界条件li, 当i=j时φi(xj,yj)=1; 当i≠j时φi(xj,yj)=0, 且φi在边界∂K呈线性变化.再采用经典有限元法取子剖分数M来求解子问题(6)的多尺度基函数φi.值得注意的是, 由于问题(6)是在局部粗单元K再进行网格剖分求解多尺度基函数, 故对其在粗尺度上进行离散和运算所消耗的计算资源较少, 且子问题(6)自身的局部微观信息通过多尺度基函数有效代入了全局尺度,基于分层网格可进一步有效刻画边界层的局部摄动信息.

a(uh,v)=(f,v), ∀v∈Uh.

(7)

Acuh=bc,

(8)

其中Ac,bc分别为相对粗网格上的总刚度矩阵和总载荷向量, 从而得到多尺度有限元解uh.

在程序执行时, 多尺度有限元法的分层网格粗剖分数NMsFEM=N由式(3)迭代所得并直接用于二维问题, 而经典有限元法的细剖分数NFEM在每个分层粗单元还须乘以子剖分数M, 即NFEM=NMsFEMM, 本文取M=4.

4 数值算例

数值模拟平台为MATLAB 2020, Intel Core i9-10900K CPU 3.7 GHz, 32 GB内存. 定义真解u与数值解ua(本文为经典有限元解ug或多尺度有限元解uh)之间误差的L2范数和H1范数为

对式(1)取b=(0,0),c=2,令真解

u=(x-1)(y-1)(exp(-x/ε)-1)(exp(-y/ε)-1),

(11)

图1 ε分别取10-1, 10-4时的真解

图2 ε=10-4时, 分层网格上经典有限元法(a)和多尺度有限元法(b)的数值解

图3 ε=10-4时, 分层网格上经典有限元法(a)和多尺度有限元法(b)的离散节点误差

随着分层网格的不断剖分加密, 2种方法在分层网格上的L2范数、H1范数及其收敛情况如表1～2所示.由表1～2可见:分层网格的单方向粗剖分数NMsFEM与NFEM大约呈2倍加密, 这是因为NMsFEM是由式(3)迭代计算所得, 故其值既可能是偶数, 也可能是奇数, 并非如一致网格或Shishkin网格一般仅呈偶数倍加密; 多尺度有限元法面向L2范数有二阶收敛, 面向H1范数有二阶以上的超收敛, 收敛速度快且求解精度高.通过本文方法求解二维奇异摄动的边界层问题,最终获得了高精度、高效率和快收敛的数值模拟.

表1 当ε=10-4时, 经典有限元法在分层网格上的L2、H1范数误差和收敛阶数

表2 当ε=10-4时, 多尺度有限元法在分层网格上的L2、H1范数误差和收敛阶数