多尺度有限元法结合分层网格模拟二维奇异摄动的两端边界层问题

doi:10.3785/j.issn.1008-9497.2022.05.007

多尺度有限元法结合分层网格模拟二维奇异摄动的两端边界层问题

孙美玲^,^,¹^,², 江山^,^,²

1.南通职业大学数学教研室，江苏南通 226007

2.南通大学理学院，江苏南通 226019

Simulation of multiscale finite element method on graded meshes for two-dimensional singularly perturbed twin boundary layers problems

SUN Meiling^,^,¹^,², JIANG Shan^,^,²

1.Department of Mathematics，Nantong Vocational University，Nantong 226007，Jiangsu Province，China

2.School of Science，Nantong University，Nantong 226019，Jiangsu Province，China

通讯作者: ORCID：https：//orcid.org/0000-0001-7983-0012，E-mail：jiangshan@ntu.edu.cn.

收稿日期: 2021-08-24

基金资助:

南通市基础科学研究指令性项目. JC2021123
国家自然科学基金面上项目. 11771224

江苏省高校青蓝工程优秀骨干教师资助项目

Received: 2021-08-24

作者简介 About authors

孙美玲（1981—），ORCID：https：//orcid.org/0000-0003-0061-5155，女，博士，副教授，主要从事偏微分方程数值解及其应用研究，E-mail：sunmeiling81@163.com. , E-mail：sunmeiling81@163.com

摘要

通过摄动系数建立分层网格，用多尺度有限元法捕捉对流扩散方程的两端边界层，研究二维奇异摄动模型。基于分层网格并利用多尺度基函数刻画了边界层的微观信息，用有限的计算资源、较短的计算时间，得到了不依赖于摄动系数、一致稳定的模拟结果。

关键词： 奇异摄动 ; 自适应网格 ; 两端边界层 ; 多尺度有限元 ; 一致稳定

Abstract

To solve a two-dimensional singularly perturbed model, a multiscale finite element method on graded meshes built from the perturbed parameter is presented for capturing the twin boundary layers of convection-diffusion equations effectively. Based on the graded meshes, the multiscale basis functions are capable of subtly describing the microscopic information in the boundary layers. No wonder, it just costs a handful of computing resource and short time to achieve the accurate and efficient results, and the results are independent of the perturbed parameter with uniform stability.

Keywords： singular perturbation ; adaptive meshes ; twin boundary layers ; multiscale finite element ; uniform stability

PDF (2888KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

孙美玲, 江山. 多尺度有限元法结合分层网格模拟二维奇异摄动的两端边界层问题. 浙江大学学报(理学版)[J], 2022, 49(5): 564-569 doi:10.3785/j.issn.1008-9497.2022.05.007

SUN Meiling, JIANG Shan. Simulation of multiscale finite element method on graded meshes for two-dimensional singularly perturbed twin boundary layers problems. Journal of Zhejiang University(Science Edition)[J], 2022, 49(5): 564-569 doi:10.3785/j.issn.1008-9497.2022.05.007

0　引言

称奇异摄动方程最高阶导数项的系数 $ε$ 为摄动参数，奇异摄动方程在弹性力学、量子力学、动力学、生物种群、最优控制等领域应用广泛。由相关研究知，当 $ε$ 较大时，不会出现奇异摄动边界层，用传统方法便可处理。当 $ε$ 较小时，往往呈现伪振荡，即边界层现象，影响数值精度。苏煜城等^［1］系统介绍了数值方法，MILLER等^［2］的研究也流传很广。

文献［3-9］主要考虑一维奇异摄动的有效求解方法。KADALBAJOO等^［3］采用三次B-样条配点法在Shishkin网格上得到了关于最大模的收敛结果。但Shishkin分片等距网格存在一定局限性，粗略估计过渡点位置可能造成方法精度不高。GENG等^［4］用再生核方法模拟两端边界层现象；杨宇博等^［5］研究非对称内罚间断有限元法在分层网格中的一致收敛性，在一维情形下对左端附近的分层网格进行加密构建。受其启发，本文将其拓展为左右两端附近皆可用自适应迭代的网格，适用于二维情形下的 $x, y$ 方向。ZHENG等^［6］用混合有限差分格式处理拟线性边值问题；基于Bakhvalov-Shishkin网格，江山等^［7］、郑权等^［8］分别得到多尺度有限元法、混合差分法的奇异摄动高精度结果；CHENG^［9］利用局部间断有限元法有效模拟了双参数模型，并基于各范数给出了稳定性估计。

对二维奇异摄动的研究已取得一定成果，如FRANZ等^［10］利用单元边界稳定化技术，证明了高阶有限元格式总能得到理想的一致超收敛；BRDAR等^［11］面向二维对流反应扩散的双参数方程，基于Duran-Lombardi与Duran-Shishkin网格，用传统双线性有限元法得到了比Shishkin网格精度更高的能量范数误差；JIANG等^［12］提出了独立构造试探函数空间、检验函数空间的Petrov-Galerkin多尺度有限元基函数，在等级网格上自动消除共振误差，得到了精确、高效、稳定的一致收敛模拟；LI等^［13］针对反应扩散方程，基于分层网格证明了高阶Galerkin有限元的优化理论；XU等^［14］、CAI等^［15］分别给出了组合有限元、高阶有限元的鲁棒性分析与数值验证，但在二维情形下这类有限元模拟的计算代价较大。

本文针对二维奇异摄动中小参数导致的两端边界层问题，利用基于分层网格的多尺度有限元计算格式，实现较传统有限元计算格式数值精度更高、计算代价更小、运算时间更短的一致稳定结果。

1　模型

考虑二维情形的对流扩散方程

\{\begin{array}{l} - ε Δ u + b \cdot \nabla u + c u = f, 在 Ω \in {(0,1)}^{2}, \\ u = g, 在 \partial Ω, \end{array}

（1）

其中， $ε$ 为奇异摄动系数， $Δ$ 为拉普拉斯算子， $b = (b_{1}, b_{2})^{T}$ 为二维向量，与梯度算子 $\nabla = {(\frac{\partial}{\partial x}, \frac{\partial}{\partial y})}^{T}$ 做点积， $c$ 为变系数， $f$ 为右端项， $g$ 为区域边界 $\partial Ω$ 的边界条件， $u$ 为求解目标。

利用虚功原理，可知式（1）的变分形式是寻求 $u \in H_{}^{1} (Ω)$ ，使得

a (u, v) = (f, v), v \in H_{0}^{1} (Ω),

（2）

其中，双线性形式为

a (u, v) = \int_{Ω} (ε \nabla u \cdot \nabla v + b \cdot \nabla u \cdot v + c u v) d x d y

，

线性内积为 $(f, v) = \int_{Ω} f v d x d y$ ，

$H_{0}^{1}$ 为一阶可导且平方可积的齐次边界函数空间。

2　分层网格与解的分解

2.1　两端边界层的分层网格

为得到式（2）变分形式的有效近似，用有限维逼近无限维的思想进行区域离散。因常规的一致网格难以有效求解奇异摄动小参数问题，即使剖分数 $N$ 很大，其等距步长 $h = 1 / N$ 也无法满足 $h < ε$ ，故难以形成可靠的分辨率。对二维区域，先在 $x$ 方向形成适合左右两端边界层的优化分层网格，再类似处理 $y$ 方向的上下两端边界层。

分层（graded）网格^［5］由迭代格式生成，用摄动系数 $ε$ 和网格参数 $0 < h \leq 1$ 计算。 $x$ 方向的节点 $x_{i}$ 为

x_{i} = \{\begin{array}{l} 0, i = 0, \\ i h ε, i = 1, \\ (1 + h) x_{i - 1}, 2 \leq i \leq N - 1, \\ 1, i = N ， \end{array}

（3）

从而形成一端稠密、另一端稀疏的分层网格。为满足两端均有边界层的情况，将式（3）改进为

x_{i} = \{\begin{array}{l} 0, i = 0, \\ i h ε, i = 1, \\ (1 + h) x_{i - 1}, 2 \leq i \leq N_{0} - 1, \\ 0.5, i = N_{0}, \\ \begin{array}{l} 1 - (1 + h) (1 - x_{i + 1}), \\ N_{0} + 1 \leq i \leq N - 2, \end{array} \\ 1 - i h ε, i = N - 1, \\ 1, i = N, \end{array}

（4）

其中， $N_{0}$ 由式（3）的 $N$ 得到。式（4）在 $x$ 方向的剖分数约为式（3）的2倍。记第 $i$ 单元 $[x_{i - 1}, x_{i}]$ 的步长 $H_{i} = x_{i} - x_{i - 1}$ ，类似记 $y$ 方向第 $j$ 单元 $[y_{j - 1}, y_{j}]$ 的步长 $H_{j} = y_{j} - y_{j - 1}$ ，边界层宽度约为 $τ_{x}, τ_{y}$ 。这样可将二维区域更好地离散剖分，以形成适合两端均有边界层的情况，且稠密稀疏程度不同的上下左右分层网格均可进行数值计算。

下文将验证分层网格是一种能自适应地逼近边界层位置及宽度的优化网格，其剖分数 $N$ 并非简单地成倍增加，从而突破了一致网格和Shishkin网格剖分数偶数倍加密的局限，得到了更好的数值精度与稳定结果。

2.2　解的多尺度分解

依据解的多尺度性质，将其分解为若干部分之和，即

u = \sum_{i = 1}^{3} \sum_{j = 1}^{3} u_{i j},

（5）

其中， $u_{i j}$ 的下标与图1区域 $Ω$ 的子分块 $Ω_{i j}$ 相对应，则各解及其导数与摄动系数 $ε$ 相关，有

|\frac{\partial^{i + j} u_{22}}{\partial x^{i} \partial y^{j}}| \leq C,

|\frac{\partial^{i + j} u_{12}}{\partial x^{i} \partial y^{j}}| \leq C e^{- \frac{b_{2} y}{ε}} ε^{- j},

|\frac{\partial^{i + j} u_{32}}{\partial x^{i} \partial y^{j}}| \leq C e^{\frac{b_{2} (y - 1)}{ε}} ε^{- j},

|\frac{\partial^{i + j} u_{21}}{\partial x^{i} \partial y^{j}}| \leq C e^{- \frac{b_{1} x}{ε}} ε^{- i},

|\frac{\partial^{i + j} u_{23}}{\partial x^{i} \partial y^{j}}| \leq C e^{\frac{b_{1} (x - 1)}{ε}} ε^{- i},

|\frac{\partial^{i + j} u_{11}}{\partial x^{i} \partial y^{j}}| \leq C e^{- \frac{b_{1} x}{ε}} e^{- \frac{b_{2} y}{ε}} ε^{- (i + j)},

|\frac{\partial^{i + j} u_{13}}{\partial x^{i} \partial y^{j}}| \leq C e^{\frac{b_{1} (x - 1)}{ε}} e^{- \frac{b_{2} y}{ε}} ε^{- (i + j)},

|\frac{\partial^{i + j} u_{31}}{\partial x^{i} \partial y^{j}}| \leq C e^{- \frac{b_{1} x}{ε}} e^{\frac{b_{2} (y - 1)}{ε}} ε^{- (i + j)},

|\frac{\partial^{i + j} u_{33}}{\partial x^{i} \partial y^{j}}| \leq C e^{\frac{b_{1} (x - 1)}{ε}} e^{\frac{b_{2} (y - 1)}{ε}} ε^{- (i + j)}

。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 区域 $Ω$ 的子分块

Fig.1 Sub-domains of domain $Ω$

用式（4）对两端边界层的分层网格进行区域离散，再采用多尺度有限元计算格式，使其更好地逼近子分块上多尺度解式（5）的局部形态。

3　有限元法与多尺度有限元法

3.1　有限元的变分原理

传统有限元法（FEM）通过分片多项式构造基底以形成有限维函数空间。如选定一组基 $ψ_{0}, ψ_{1},$ $ψ_{2}, \dots, ψ_{n}$ ，记有限元空间 $V^{h} = s p a n {ψ_{i}}_{i = 0}^{n} \subset H_{0}^{1}$ ，其变分形式对应为寻求 $u_{g} \in V^{h}$ ，使得

a (u_{g}, v) = (f, v), v \in V^{h},

（6）

其中， $u_{g}$ 为Galerkin有限元解。利用等参变换 $ξ = \frac{x - x_{i - 1}}{H_{i}}$ ， $η = \frac{y - y_{j - 1}}{H_{j}}$ ，4个局部双线性基函数依次为 $\hat{ψ_{1}} = (1 - ξ) (1 - η)$ ， $\hat{ψ_{2}} = ξ (1 - η)$ ， $\hat{ψ_{3}} = ξ η$ ， $\hat{ψ_{4}} = (1 - ξ) η$ ，因此局部有限元的解 $\hat{u_{g}} = \sum_{i = 1}^{4} \hat{u_{i}} \hat{ψ_{i}}$ ，其中 $\hat{u_{i}}$ 为局部节点离散值。通过局部节点和全局节点的对应数据结构，形成总刚度矩阵和总载荷向量，求解代数方程组，得到有限元解 $u_{g}$ 。

当摄动系数 $ε$ 很小时，传统方法无法进行高效处理，即使 $x, y$ 方向的剖分数 $N_{F E M}$ 取很大，二维问题仍达 $O (N_{F E M}^{2})$ 量级，计算消耗巨大，难以得到可靠的数值精度，存在计算局限。

3.2　多尺度有限元的变分原理

不同于传统有限元，多尺度有限元法（MsFEM）在构造有限维空间时，不采用显式多项式函数，而采用基于与原问题相同的微分算子在粗风格单元中求解非显式基函数。在每个粗网格单元 $K$ 中用有限元法求对应的齐次子问题：

\{\begin{array}{l} - ε Δ ϕ_{i} + b \cdot \nabla ϕ_{i} + c ϕ_{i} = 0, 在 K, \\ ϕ_{i} = θ_{i}, 在 \partial K ， \end{array}

（7）

其中， $θ_{i}$ 为单元边界上的线性边界条件。可知，式（1）是直接用传统有限元法在非常密的二维网格上求解有限元，而式（7）是间接用有限元法在较粗糙的二维网格上求多尺度基函数，其子单元剖分数为 $M$ 。

记多尺度有限元空间为 $U^{h} = s p a n {ϕ_{i}}_{i = 0}^{n}$ $\subset$ $H_{0}^{1}$ ，用对应的变分形式寻求 $u_{h} \in U^{h}$ ，使得

a (u_{h}, v) = (f, v), v \in U^{h},

（8）

其中， $u_{h}$ 为多尺度有限元解。

显然，传统有限元空间 $V^{h} = s p a n {ψ_{i}}_{i = 0}^{n}$ 与多尺度有限元空间 $U^{h} = s p a n {ϕ_{i}}_{i = 0}^{n}$ 不同，两者由不同的基底张成。此外，传统有限元法需在非常密的二维网格 $O (N_{F E M}^{2})$ 上求解方程组

A_{f} u_{g} = b_{f},

（9）

其中， $A_{f}, b_{f}$ 分别表示细网格上的总刚度矩阵和总右端载荷。多尺度有限元法是在较粗的 $O (N_{M s F E M}^{2})$ （其 $x, y$ 方向剖分数 $N_{M s F E M}^{} ≪ N_{F E M}^{}$ ）上求解方程组

A_{c} u_{h} = b_{c},

（10）

其中， $A_{c}, b_{c}$ 分别表示粗网格上的总刚度矩阵和总右端载荷，得到多尺度解 $u_{h}$ 。可以验证，多尺度有限元法的离散化多尺度基函数有效地刻画了边界层的微观信息，计算代价更小，计算精度有时较传统有限元法更高，对处理奇异摄动边界层问题具有模拟优势。

3.3　多尺度有限元解的误差估计

当真解 $u$ 用式（8）中的多尺度有限元解 $u_{h}$ 近似时，可用两者之间的能量范数

{‖u - u_{h}‖}_{E n e r g y}^{2} = {‖u - u_{h}‖}_{L^{2}}^{2} + ε {‖u - u_{h}‖}_{H^{1}}^{2}

（11）

度量误差，结合分层网格得到收敛结果。

定理1

{‖u - u_{h}‖}_{E n e r g y}^{} \leq [C \frac{1}{N^{α}} + ε^{\frac{1}{2}} {(\frac{l n N}{N})}^{α + 1} m a x {|ω^{'}|}^{p + \frac{1}{2}}],

（12）

其中， $α = p + 1$ 为收敛阶，当 $p = 1$ 时为线性基函数，

ω = (C_{1} e^{- \frac{x_{i}}{α ε}} + C_{2} e^{\frac{x_{i} - 1}{α ε}}) (C_{3} e^{- \frac{y_{j}}{α ε}} + C_{4} e^{\frac{y_{j} - 1}{α ε}})

为二维分层网格的生成函数。

4　数值验证

用已有文献算例和程序结果验证相应方法的精度和效率，本文仅讨论当 $ε$ 很小时产生奇异摄动边界层求解困境的情况。将分层网格上与传统有限元、多尺度有限元对应的结果分别记作FEM（G）和MsFEM（G），用真解、近似解和误差的三维图示、范数值分析度量实际模拟效果。

令 $b = {(- 2, - 2)}^{T}$ ，变系数 $c = e^{x + y}$ ，式（1）的真解来自文献［13］，为

\begin{array}{l} u (x, y) = (1 - e^{- \frac{x}{ε}}) (1 - e^{\frac{x - 1}{ε}}) \times \\ (1 - e^{- \frac{y}{ε}}) (1 - e^{\frac{y - 1}{ε}}) 。 \end{array}

（13）

当 $ε = 10^{- 1} 和 10^{- 5}$ 时的真解见图2。可知，当 $ε$ 较大时，解光滑且无边界层现象；当 $ε$ 很小时，边界层出现了上下左右四边界。分别用传统有限元法与多尺度有限元法求解，2种方法结合自适应的分层网格均能很好地模拟 $ε = 10^{- 5}$ 时的真解，见图3。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 当 $ε$ = $10^{- 1} 和 10^{- 5}$ 时的真解

Fig.2 Exact solutions when $ε = 10^{- 1}$ and $10^{- 5}$ ，respectively

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 当 $ε = 10^{- 5}$ 时传统有限元法与多尺度有限元法的解

Fig.3 The solutions of FEM（G） and MsFEM（G） when $ε = 10^{- 5}$

为更清晰地展现相应方法的精确性与稳定性，通过网格加密的方法观察数值变化。表格格式与文献［5］的表1与表2一致，区别在于文献［5］处理的是一维问题，行数较多、单方向剖分数较小，本文研究的是二维问题，因受算力限制，行数较少、单方向两端边界层的剖分数较大。由表1知，无论摄动参数如何选取，依据网格参数 $h$ 的递减 $ε$ ，由迭代式（4）自适应生成两端疏密不同的分层网格，用于离散化计算。横向看，表1中单方向剖分数 $N$ 、范数误差均微增，纵向看，其范数误差随 $h$ 递减呈稳定收敛。表1为用传统有限元法求解二维问题，若网格剖分数 $N_{F E M}^{}$ 较大，计算消耗很大，继续剖分将超出单机的运行限定，具有一定的局限性。表2采用的是多尺度有限元法，仅在较粗的分层网格 $N_{M s F E M}^{}$ 上计算误差的能量范数，求式（7）时其子单元剖分数 $M = 4$ ，对应行的精度略逊于表1，但其计算消耗小、时间短，继续剖分能得到更精确的结果，其收敛结果与理论估计式（12）一致。另需指出，多尺度有限元法的精度在剖分数 $N_{M s F E M}^{}$ =222，262，300时较传统有限元法在剖分数 $N_{F E M}^{}$ =240，288，328时更高。当然此优化结果也有计算消耗，主要用以刻画奇异摄动问题的边界层微观属性。

表1 不同 $ε$ 下传统有限元法在分层网格上的误差能量范数

Table 1 FEM（G） with different parameters for errors of energy norm

h	$ε = 10^{- 5}$		$ε = 10^{- 6}$		$ε = 10^{- 7}$
h	$N_{F E M}^{}$	误差	$N_{F E M}^{}$	误差	$N_{F E M}^{}$	误差
$2^{+ 0}$	136	4.510×10^-3	160	4.664×10^-3	192	4.748×10^-3
$2^{- 1}$	240	1.310×10^-3	288	1.362×10^-3	328	1.402×10^-3
$2^{- 2}$	448	3.870×10^-4	536	4.054×10^-4	616	4.189×10^-4
$2^{- 3}$	888	1.083×10^-4	1 048	1.142×10^-4	1 200	1.186×10^-4

新窗口打开| 下载CSV

表2 不同 $ε$ 下多尺度有限元法在分层网格上的误差能量范数

Table 2 MsFEM（G） with different parameters for errors of energy norm

h	$ε = 10^{- 5}$		$ε = 10^{- 6}$		$ε = 10^{- 7}$
h	$N_{M s F E M}^{}$	误差	$N_{M s F E M}^{}$	误差	$N_{M s F E M}^{}$	误差
$2^{+ 0}$	34	8.866×10^-2	40	9.084×10^-2	48	9.083×10^-2
$2^{- 1}$	60	2.270×10^-2	72	2.303×10^-2	82	2.329×10^-2
$2^{- 2}$	112	4.686×10^-3	134	4.938×10^-3	154	4.955×10^-3
$2^{- 3}$	222	9.228×10^-4	262	9.565×10^-4	300	9.627×10^-4

新窗口打开| 下载CSV

当 $ε = 10^{- 5}$ 时，2种数值方法在大致相当的剖分数 $N_{F E M}^{2} = 136^{2}$ 和 $N_{M s F E M}^{2} = 112^{2}$ 下的对应误差直观比较见图4。可以验证，传统有限元法与多尺度有限元法结合分层网格均有效逼近了 $x$ 轴、 $y$ 轴、 $x = 1$ 、 $y = 1$ 的四边边界层，两者的误差上限均有效控制在1%内，展示了很好的模拟效果。再对应给出剖分数 $N_{F E M}^{2} = 240^{2}$ 和 $N_{M s F E M}^{2} = 222^{2}$ 的误差比较，两者上限分别为 $3 \times 10^{- 3}$ 和 $8 \times 10^{- 3}$ ，且后者在全局范围误差实现了整体消磨，边界层误差不如前者明显，但在区域中部表现出一定的振荡，如图5所示。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 当 $ε = 10^{- 5}$ 时传统有限元法于 $N_{F E M}^{} = 136$ 和多尺度有限元法于 $N_{M s F E M}^{} = 112$ 的误差

Fig.4 Errors of FEM（G） on $N_{F E M}^{} = 136$ and MsFEM（G） on $N_{M s F E M}^{} = 112$ when $ε = 10^{- 5}$

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 当 $ε = 10^{- 5}$ 时传统有限元法于 $N_{F E M}^{} = 240$ 和多尺度有限元法于 $N_{M s F E M}^{} = 222$ 的误差

Fig.5 Errors of FEM（G） on $N_{F E M}^{} = 240$ and MsFEM（G） on $N_{M s F E M}^{} = 222$ when $ε = 10^{- 5}$

在上述数值精度与稳定分析的基础上，考虑2种数值方法所需的运行时间和效率，表3给出了当 $ε = 10^{- 5}$ 时台式机Intel Core i9 CPU 3.7 GHz运行相应程序所需的CPU时间，可见在较密二维网格上用传统有限元法所需的CPU时间是同一行较粗二维网格用多尺度有限元法的近10倍，显然多尺度有限元法的计算效率更高。进一步，图6为摄动参数取更小（ $10^{- 6}$ 与 $10^{- 7}$ ）时相应方法的CPU时间与剖分数的对数比例关系，再次证实了多尺度有限元法的计算代价更小、计算效率更高。

表3 当 $ε = 10^{- 5}$ 时传统有限元法与多尺度有限元法的CPU时间

Table 3 FEM（G） and MsFEM（G）'s CPU time when $ε = 10^{- 5}$

$N_{F E M}^{2}$	FEM（G）CPU时间/s	$N_{M s F E M}^{2}$	MsFEM（G）CPU时间/s
$136^{2}$	4.1	$34^{2}$	3.3
$240^{2}$	49	$60^{2}$	12
$448^{2}$	659	$112^{2}$	84
$888^{2}$	10 277	$222^{2}$	1 106

新窗口打开| 下载CSV

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 $ε$ = $10^{- 6}$ 和 $10^{- 7}$ 时2种方法的剖分数 $N$ 与CPU时间的log-log图示

Fig.6 Two methods' log-log on partition $N$ and CPU time when $ε = 10^{- 6}$ and $10^{- 7}$

综上所述，多尺度有限元法只需在较粗分层网格上进行计算，消耗的计算资源较少，且能保证稳定收敛的有效精度，因此，多尺度有限元法在高维奇异摄动问题求解中具有广阔的应用前景。

5　结束语

基于自适应的分层网格生成机制，主要利用多尺度有限元法处理奇异摄动的二维对流扩散变系数方程。用分层迭代精确逼近边界层位置与宽度，结合多尺度计算格式有效捕捉了两端边界层的微观信息，实现了不依赖摄动系数的精确高效模拟结果，充分展现了多尺度有限元法结合分层网格求解高维奇异摄动问题的一致稳定性和优势。

http://dx.doi.org/10.3785/j.issn.1008-9497.2022.05.007

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

苏煜城，吴启光. 奇异摄动问题数值方法引论［M］. 重庆：重庆出版社， 1991.

[本文引用: 1]

Y C

， WU

Q G

. An Introduction to Numerical Methods for the Singular Perturbation Problems［M］. Chongqing： Chongqing Publishing House， 1991.

[本文引用: 1]

[2]

MILLER

J J H

， O’RIORDAN

， SHISHKIN

G I

. Fitted Numerical Methods for Singular Perturbation Problems （Revised Edition）［M］. Singapore： World Scientific， 2012. doi:10.1142/8410

[本文引用: 1]

[3]

KADALBAJOO

M K

， GUPTA

A parameter uniform B-spline collocation method for solving singularly perturbed turning point problem having twin boundary layers

［J］. International Journal of Computer Mathematics， 2010， 87（14）： 3218-3235. doi:10.1080/00207160902980492

[本文引用: 2]

[4]

GENG

F Z

， QIAN

S P

Reproducing kernel method for singularly perturbed turning point problems having twin boundary layers

［J］. Applied Mathematics Letters， 2013， 26（10）： 998-1004. DOI：10.1016/j.aml.2013.05.006

[本文引用: 1]

[5]

杨宇博，祝鹏，尹云辉.

分层网格上奇异摄动问题的一致NIPG分析

［J］. 计算数学， 2014， 36（4）： 437-448. doi:10.12286/jssx.2014.4.437

[本文引用: 4]

YANG

Y B

， ZHU

， YIN

Y H

Uniform analysis of the NIPG method on graded meshes for singularly perturbed convection-diffusion problems

［J］. Mathematica Numerica Sinica， 2014， 36（4）： 437-448. doi:10.12286/jssx.2014.4.437

[本文引用: 4]

[6]

ZHENG

， LI

X Z

， GAO

Uniformly convergent hybrid schemes for solutions and derivatives in quasilinear singularly perturbed BVPs

［J］. Applied Numerical Mathematics， 2015， 91： 46-59. DOI：10. 1016/j.apnum.2014.12.010

[本文引用: 1]

[7]

江山，孙美玲.

多尺度有限元结合Bakhvalov-Shishkin网格法高效处理边界层问题

［J］. 浙江大学学报（理学版）， 2015， 42（2）： 142-146. DOI：10. 3785/j.issn.1008-9497.2015.02.004

[本文引用: 1]

JIANG

， SUN

M L

Combining the multiscale finite element and Bakhvalov-Shishkin grid to solve the boundary layer problems

［J］. Journal of Zhejiang University （Science Edition）， 2015， 42（2）： 142-146. DOI：10.3785/j.issn.1008-9497.2015.02.004

[本文引用: 1]

[8]

郑权，刘颖，刘忠礼.

奇异摄动问题在修正的Bakhvalov-Shishkin网格上的混合差分格式

［J］. 浙江大学学报（理学版）， 2020， 47（4）： 460-468. DOI：10. 3785/j.issn.1008-9497.2020.04.009

[本文引用: 1]

ZHENG

， LIU

Z L

The hybrid finite difference schemes on the modified Bakhvalov-Shishkin mesh for the singularly perturbed problem

［J］. Journal of Zhejiang University （Science Edition）， 2020， 47（4）： 460-468. DOI：10.3785/j.issn.1008-9497.2020. 04.009

[本文引用: 1]

[9]

CHENG

On the local discontinuous Galerkin method for singularly perturbed problem with two parameters

［J］. Journal of Computational and Applied Mathematics， 2021， 392： 113485. DOI：10. 1016/j.cam.2021.113485

[本文引用: 2]

[10]

FRANZ

， LINß

， ROOS

H G

， et al.

Uniform superconvergence of a finite element method with edge stabilization for convection-diffusion problems

［J］. Journal of Computational Mathematics， 2010， 28（1）： 32-44. DOI：10.4208/jcm.2009.09-m1005

[本文引用: 1]

[11]

BRDAR

， ZARIN

， TEOFANOV

A singularly perturbed problem with two parameters in two dimensions on graded meshes

［J］. Computers and Mathematics with Applications， 2016， 72（10）： 2582-2603. DOI：10.1016/j.camwa.2016.09.021

[本文引用: 1]

[12]

JIANG

， PRESHO

， HUANG

Y Q

An adapted Petrov-Galerkin multi-scale finite element method for singularly perturbed reaction-diffusion problems

［J］. International Journal of Computer Mathematics， 2016， 93（7）： 1200-1211. DOI：10.1080/00207160. 2015.1041935

[本文引用: 1]

[13]

Z W

， WU

， XU

Y S

High order Galerkin methods with graded meshes for two-dimensional reaction-diffusion problems

［J］. International Journal of Numerical Analysis and Modeling， 2016， 13（3）： 319-343.

[本文引用: 2]

[14]

S P

， DENG

W B

， WU

H J

A combined finite element method for elliptic problems posted in domains with rough boundaries

［J］. Journal of Computational and Applied Mathematics， 2018， 336： 235-248. DOI：10.1016/j.cam.2017.12.049

[本文引用: 1]

[15]

CAI

D F

， CAI

Z Q

， ZHANG

Robust equilibrated a posteriori error estimator for higher order finite element approximations to diffusion problems

［J］. Numerische Mathematik， 2020， 144： 1-21. DOI：10.1007/s00211-019-01075-1

[本文引用: 1]

1991

... 称奇异摄动方程最高阶导数项的系数

ε

为摄动参数，奇异摄动方程在弹性力学、量子力学、动力学、生物种群、最优控制等领域应用广泛.由相关研究知，当

ε

较大时，不会出现奇异摄动边界层，用传统方法便可处理.当

ε

较小时，往往呈现伪振荡，即边界层现象，影响数值精度.苏煜城等^［1］系统介绍了数值方法，MILLER等^［2］的研究也流传很广. ...

1991

... 称奇异摄动方程最高阶导数项的系数

ε

为摄动参数，奇异摄动方程在弹性力学、量子力学、动力学、生物种群、最优控制等领域应用广泛.由相关研究知，当

ε

较大时，不会出现奇异摄动边界层，用传统方法便可处理.当

ε

较小时，往往呈现伪振荡，即边界层现象，影响数值精度.苏煜城等^［1］系统介绍了数值方法，MILLER等^［2］的研究也流传很广. ...

2012

... 称奇异摄动方程最高阶导数项的系数

ε

为摄动参数，奇异摄动方程在弹性力学、量子力学、动力学、生物种群、最优控制等领域应用广泛.由相关研究知，当

ε

较大时，不会出现奇异摄动边界层，用传统方法便可处理.当

ε

较小时，往往呈现伪振荡，即边界层现象，影响数值精度.苏煜城等^［1］系统介绍了数值方法，MILLER等^［2］的研究也流传很广. ...

A parameter uniform B-spline collocation method for solving singularly perturbed turning point problem having twin boundary layers

2010

... 文献［3-9］主要考虑一维奇异摄动的有效求解方法.KADALBAJOO等^［3］采用三次B-样条配点法在Shishkin网格上得到了关于最大模的收敛结果.但Shishkin分片等距网格存在一定局限性，粗略估计过渡点位置可能造成方法精度不高.GENG等^［4］用再生核方法模拟两端边界层现象；杨宇博等^［5］研究非对称内罚间断有限元法在分层网格中的一致收敛性，在一维情形下对左端附近的分层网格进行加密构建.受其启发，本文将其拓展为左右两端附近皆可用自适应迭代的网格，适用于二维情形下的

x, y

方向.ZHENG等^［6］用混合有限差分格式处理拟线性边值问题；基于Bakhvalov-Shishkin网格，江山等^［7］、郑权等^［8］分别得到多尺度有限元法、混合差分法的奇异摄动高精度结果；CHENG^［9］利用局部间断有限元法有效模拟了双参数模型，并基于各范数给出了稳定性估计. ...

... ［3］采用三次B-样条配点法在Shishkin网格上得到了关于最大模的收敛结果.但Shishkin分片等距网格存在一定局限性，粗略估计过渡点位置可能造成方法精度不高.GENG等^［4］用再生核方法模拟两端边界层现象；杨宇博等^［5］研究非对称内罚间断有限元法在分层网格中的一致收敛性，在一维情形下对左端附近的分层网格进行加密构建.受其启发，本文将其拓展为左右两端附近皆可用自适应迭代的网格，适用于二维情形下的

x, y

Reproducing kernel method for singularly perturbed turning point problems having twin boundary layers

2013

x, y

分层网格上奇异摄动问题的一致NIPG分析

2014

x, y

... 分层（graded）网格^［5］由迭代格式生成，用摄动系数

ε

和网格参数

0 < h \leq 1

计算.

x

方向的节点

x_{i}

为 ...

... 为更清晰地展现相应方法的精确性与稳定性，通过网格加密的方法观察数值变化.表格格式与文献［5］的表1与表2一致，区别在于文献［5］处理的是一维问题，行数较多、单方向剖分数较小，本文研究的是二维问题，因受算力限制，行数较少、单方向两端边界层的剖分数较大.由表1知，无论摄动参数如何选取，依据网格参数

h

的递减

ε

，由迭代式（4）自适应生成两端疏密不同的分层网格，用于离散化计算.横向看，表1中单方向剖分数

N

、范数误差均微增，纵向看，其范数误差随

h

递减呈稳定收敛.表1为用传统有限元法求解二维问题，若网格剖分数

N_{F E M}^{}

较大，计算消耗很大，继续剖分将超出单机的运行限定，具有一定的局限性.表2采用的是多尺度有限元法，仅在较粗的分层网格

N_{M s F E M}^{}

上计算误差的能量范数，求式（7）时其子单元剖分数

M = 4

，对应行的精度略逊于表1，但其计算消耗小、时间短，继续剖分能得到更精确的结果，其收敛结果与理论估计式（12）一致.另需指出，多尺度有限元法的精度在剖分数

N_{M s F E M}^{}

=222，262，300时较传统有限元法在剖分数

N_{F E M}^{}

=240，288，328时更高.当然此优化结果也有计算消耗，主要用以刻画奇异摄动问题的边界层微观属性. ...

... 一致，区别在于文献［5］处理的是一维问题，行数较多、单方向剖分数较小，本文研究的是二维问题，因受算力限制，行数较少、单方向两端边界层的剖分数较大.由表1知，无论摄动参数如何选取，依据网格参数

h

的递减

ε

，由迭代式（4）自适应生成两端疏密不同的分层网格，用于离散化计算.横向看，表1中单方向剖分数

N

、范数误差均微增，纵向看，其范数误差随

h

递减呈稳定收敛.表1为用传统有限元法求解二维问题，若网格剖分数

N_{F E M}^{}

较大，计算消耗很大，继续剖分将超出单机的运行限定，具有一定的局限性.表2采用的是多尺度有限元法，仅在较粗的分层网格

N_{M s F E M}^{}

上计算误差的能量范数，求式（7）时其子单元剖分数

M = 4

N_{M s F E M}^{}

=222，262，300时较传统有限元法在剖分数

N_{F E M}^{}

=240，288，328时更高.当然此优化结果也有计算消耗，主要用以刻画奇异摄动问题的边界层微观属性. ...

分层网格上奇异摄动问题的一致NIPG分析

2014

x, y

... 分层（graded）网格^［5］由迭代格式生成，用摄动系数

ε

和网格参数

0 < h \leq 1

计算.

x

方向的节点

x_{i}

为 ...

h

的递减

ε

，由迭代式（4）自适应生成两端疏密不同的分层网格，用于离散化计算.横向看，表1中单方向剖分数

N

、范数误差均微增，纵向看，其范数误差随

h

递减呈稳定收敛.表1为用传统有限元法求解二维问题，若网格剖分数

N_{F E M}^{}

较大，计算消耗很大，继续剖分将超出单机的运行限定，具有一定的局限性.表2采用的是多尺度有限元法，仅在较粗的分层网格

N_{M s F E M}^{}

上计算误差的能量范数，求式（7）时其子单元剖分数

M = 4

N_{M s F E M}^{}

=222，262，300时较传统有限元法在剖分数

N_{F E M}^{}

=240，288，328时更高.当然此优化结果也有计算消耗，主要用以刻画奇异摄动问题的边界层微观属性. ...

h

的递减

ε

，由迭代式（4）自适应生成两端疏密不同的分层网格，用于离散化计算.横向看，表1中单方向剖分数

N

、范数误差均微增，纵向看，其范数误差随

h

递减呈稳定收敛.表1为用传统有限元法求解二维问题，若网格剖分数

N_{F E M}^{}

较大，计算消耗很大，继续剖分将超出单机的运行限定，具有一定的局限性.表2采用的是多尺度有限元法，仅在较粗的分层网格

N_{M s F E M}^{}

上计算误差的能量范数，求式（7）时其子单元剖分数

M = 4

N_{M s F E M}^{}

=222，262，300时较传统有限元法在剖分数

N_{F E M}^{}

=240，288，328时更高.当然此优化结果也有计算消耗，主要用以刻画奇异摄动问题的边界层微观属性. ...

Uniformly convergent hybrid schemes for solutions and derivatives in quasilinear singularly perturbed BVPs

2015

x, y

多尺度有限元结合Bakhvalov-Shishkin网格法高效处理边界层问题

2015

x, y

多尺度有限元结合Bakhvalov-Shishkin网格法高效处理边界层问题

2015

x, y

奇异摄动问题在修正的Bakhvalov-Shishkin网格上的混合差分格式

2020

x, y

奇异摄动问题在修正的Bakhvalov-Shishkin网格上的混合差分格式

2020

x, y

On the local discontinuous Galerkin method for singularly perturbed problem with two parameters

2021

x, y

... ［9］利用局部间断有限元法有效模拟了双参数模型，并基于各范数给出了稳定性估计. ...

Uniform superconvergence of a finite element method with edge stabilization for convection-diffusion problems

2010

... 对二维奇异摄动的研究已取得一定成果，如FRANZ等^［10］利用单元边界稳定化技术，证明了高阶有限元格式总能得到理想的一致超收敛；BRDAR等^［11］面向二维对流反应扩散的双参数方程，基于Duran-Lombardi与Duran-Shishkin网格，用传统双线性有限元法得到了比Shishkin网格精度更高的能量范数误差；JIANG等^［12］提出了独立构造试探函数空间、检验函数空间的Petrov-Galerkin多尺度有限元基函数，在等级网格上自动消除共振误差，得到了精确、高效、稳定的一致收敛模拟；LI等^［13］针对反应扩散方程，基于分层网格证明了高阶Galerkin有限元的优化理论；XU等^［14］、CAI等^［15］分别给出了组合有限元、高阶有限元的鲁棒性分析与数值验证，但在二维情形下这类有限元模拟的计算代价较大. ...

A singularly perturbed problem with two parameters in two dimensions on graded meshes

2016

An adapted Petrov-Galerkin multi-scale finite element method for singularly perturbed reaction-diffusion problems

2016

High order Galerkin methods with graded meshes for two-dimensional reaction-diffusion problems

2016

... 令

b = {(- 2, - 2)}^{T}

，变系数

c = e^{x + y}

，式（1）的真解来自文献［13］，为 ...

A combined finite element method for elliptic problems posted in domains with rough boundaries

2018

Robust equilibrated a posteriori error estimator for higher order finite element approximations to diffusion problems

2020

〈

〉