一种用于对流-扩散问题的SUPG稳定化尺度边界有限元方法

《Computers & Mathematics with Applications》：A SUPG stabilized scaled boundary finite element method for convection-diffusion problems

　　M. Arrutselvi | V.S. Suvin | Ean Tat Ooi | Chongmin Song | Sundararajan Natarajan
机械工程系，印度理工学院马德拉斯分校，金奈，600036，泰米尔纳德邦，印度

摘要
在本文中，我们将基于流线迎风Petrov-Galerkin（SUPG）稳定化的缩放边界有限元方法（SBFEM）扩展到稳态对流-扩散问题。SBFEM提供了一个半解析框架，自然地支持具有悬挂节点的元素，并且不对用于空间离散化的元素的形状和质量施加限制，从而减轻了网格划分的负担。另一方面，SUPG在保持精度的同时确保了解的稳定性，并抑制了数值不稳定性。通过一系列复杂度不同的数值示例，展示了将SBFEM的计算效率和几何灵活性与SUPG稳定性的稳健性相结合的潜力。从系统的数值研究可以推断出，所提出的框架表现出改进的性能，并且结果以最佳的收敛速率收敛。

1. 引言
对流-扩散（CD）方程描述了一种物理现象，其中标量量（如热量、质量或化学浓度）的移动是由对流和扩散过程共同作用的结果。该方程包含两个主要项：对流项，涉及方向性传输；扩散项，涉及局部混合。CD方程在许多领域都有广泛的应用，包括环境工程[1]、[2]、流体力学[3]、[4]、生物医学工程[5]以及地球物理学[6]、[7]。典型的应用包括固体或流体中的热传递[8]、[9]、水或空气中的污染物扩散[10]、[11]以及半导体建模[12]、[13]。

传统的解决方法，如变量分离、相似变换或拉普拉斯变换，可以用来求解CD方程。尽管这些方法提供了有价值的见解，但它们仅适用于简单的几何形状和边界条件。然而，大多数实际场景涉及复杂的几何形状、非均质介质和非线性源项，在这些情况下寻找解析解是不可行的。已经开发并应用了多种数值方法来求解CD方程，其中包括有限元方法（FEM）[14]、[15]、虚拟元素方法（VEM）[16]、[17]、[18]、[19]、[20]、有限体积方法（FVM）[21]、[22]、[23]、有限差分方法（FDM）[24]、[25]、等几何分析（IGA）[26]、双网格方法[27]以及无网格方法[28]、[29]、[30]。当问题以对流为主时，上述数值方法面临重大挑战[31]。根据大量研究，使用经典方法求解CD方程得到的数值解往往不稳定，在陡峭的梯度附近会出现振荡或不准确的结果。这促使研究人员开发出能够减轻非物理数值振荡的数值技术。观察发现，虚假振荡的出现归因于：(a) 网格扭曲；(b) 网格不够细，无法捕捉到陡峭的梯度。研究者已经研究了几种自适应网格[32]、[33]、[34]。自适应网格技术的一个主要限制是需要预先了解分层区域的信息。此外，使用自适应网格细化技术并不会显著增加经典数值方法的收敛范围[35]。

稳定化策略通过在Galerkin公式中添加惩罚偏离流线解的行为的项来增强其效果。所使用的不同稳定化技术包括：Galerkin/最小二乘（GLS）方法[36]、[37]、伴随稳定化（ADJ）[36]、代数子网格尺度（ASGS）[38]、局部投影稳定化[39]以及流线迎风Petrov-Galerkin（SUPG）[40]、[41]、[42]。其中，SUPG因其简单性、有效性和一致性而被最广泛采用。SUPG方法减少了振荡，从而提高了结果的准确性。它通过适用于稳态和瞬态方程来证明其适应性。采用SUPG还可以降低计算成本，因为较粗的网格也能产生满意的结果[43]。

在数值方法中，VEM[18]、[19]、[20]受到了特别关注。这可能归因于以下几点：(a) 它放宽了网格划分的要求，因此允许使用任意形状和大小的元素；(b) 可以轻松扩展到高阶元素和三维空间；(c) 最重要的是，不需要显式的形状函数形式。另一方面，缩放边界有限元方法（SBFEM）最近也引起了人们的兴趣。SBFEM最初由Song和Wolf[44]、[45]针对无界域问题提出。根据特定应用，可以使用多边形元素[46]、自适应网格技术（如四叉树或八叉树网格[47]、[48]）或非常规网格方法[49]与SBFEM结合使用。有兴趣的读者可以参考关于SBFEM的综述文章[50]及其中的参考文献。

如上所述，SBFEM已广泛应用于弹性方程，类似于VEM，它也可以适应任意形状和大小的元素；然而，在文献中关于其在对流-扩散方程中的应用较少。SBFEM的一个特点是，它不需要稳定化项来计算双线性/线性形式，并且积分仅限于边界。所提出框架的主要目标是：
? 将SBFEM扩展到对流-扩散方程——通过采用缩放边界形状函数实现
? 结合SUPG稳定化来抑制对流主导问题中的虚假数值振荡

该框架已通过内部代码实现，并对其准确性进行了数值研究。此外，还研究了在不同空间离散化条件下（即变形的四边形元素、任意Voronoi多边形和变形的凸多边形元素）SUPG稳定化SBFEM框架的适用性和收敛特性。

本文的结构如下：第2节详细介绍了理论公式及SUPG稳定化方法。第3节阐述了使用缩放边界形状函数的SBFEM方法。第4节进行了数值分析，以评估所创建框架的收敛性和精度。该节包含了五个数值测试案例。

2. 模型问题
考虑一个二维有界域Ω?R2，其边界由?Ω表示。一个控制对流-扩散现象的偏微分方程具有简单的齐次Dirichlet边界条件，表达式为：
(1) ??Δu + v·?u = f ∈ Ω,
u = 0 ∣ ?Ω,
其中u = u(x, y)是感兴趣的变量，系数? = ?(x, y)是扩散率参数，v = [v1(x, y), v2(x, y)]是速度场，f = f(x, y)是线性源或汇函数。为了确保问题(1)的理论适定性，考虑以下规则性假设：? ∈ L∞(Ω)且?(x, y) ≥ ?0 ? 0，v ∈ [L∞(Ω)]2且?·v = 0（即流动是不可压缩的），f ∈ L2(Ω)。方程(1)的经典变分公式为：
寻找u ∈ H?1(Ω)，使得：
(2) ∫Ω ??u·?wx + ∫Ω(v·?u)wdx = ∫Ω fwdx，?w ∈ H?1(Ω)。
方程(2)有唯一解[51]。

设{Th}h > 0表示Ω的一系列星形多边形划分。Th中的多边形片段在SBFEM中称为S元素。对于S元素E ∈ Th，令hE表示元素直径，h定义为h := max{e ∈ Th | h(e) ≥ η/2}。缩放边界有限元空间Vh定义为：
(3) Vh := {? ∈ H?1(Ω) : ?|E ∈ H?(E) ∩ B?(?E) ∧ Δ? = 0 ∣ E ∈ Th}，
其中B?(?E)是边界?E上的连续分段线性多项式集合。在每个S元素上，拉普拉斯方程的解导致线性完备的缩放边界形状函数。形状函数及其实现的详细推导在第3节中给出。

Galerkin方法在数值解中会导致振荡，尤其是在对流现象占主导地位时。为了减轻振荡，提出了SUPG方法来稳定缩放边界有限元框架。在这种情况下，方程(2)的离散SUPG稳定化SBFEM公式为：
寻找uh ∈ Vh，使得：
(4) Asupg(uh, wh) = F(wh) ? wh ∈ Vh，
其中
Asupg(u, w) = ∑?∈Th {∫? ??u·?wx + ∫?(v·?u)wdx + ∫? τ?(e) (??Δu + v·?u)(v·?w)dx}，
F(w) = ∑?∈Th {∫? efwdx + ∫? τ? ef(v·?w)dx}，
τ?是依赖于hE和v的局部稳定化参数。在这项工作中，τE的选择（参考[52]）为：
τE := hE2/|v| (cosh(Pe) ? Pe^(-1))，
其中Pe := |v|hE|2?。
Péclet数Pe > 1表明问题以对流为主，而Pe < 1表明问题以扩散为主。

3. 缩放边界有限元方法概述
3.1. 坐标变换
在缩放边界坐标系中，坐标轴定义为径向（ξ）和周向（η）。选择缩放中心(x0, y0)作为坐标原点，使得整个边界从这个点/原点可见。对于每个一维线元素，坐标范围为0 ≤ ξ ≤ 1和?1 ≤ η ≤ 1，如图1所示。

考虑笛卡尔坐标系中的点P(x)，将其转换到缩放边界坐标系中：
(5) x(ξ, η) = x? + ξN(η)，
xb := x? + ξ，
y(ξ, η) = y? + ξN(η)，
yb := y? + ξ，
其中N(η) = [N?(η), N?(η), …, N_m(η)]是边界上线元素的形状函数，m是每个线元素的节点数。xb和yb表示线元素的节点坐标。微分体积dΩ和微分面积dΓ转换到缩放边界坐标系后表示为：
(6) dΩ = J_p(η)ξ dξ/dη，
dΓ = ξL(η)dη，
其中L(η) = x2 + y2，J_p是在边界定义的雅可比矩阵Jp的行列式：
(7) J_p = [xηy, yηx, yη]，
xη, η和yη, η分别是xη和yη的导数。

3.2. 缩放边界形状函数的推导
本节详细推导了缩放边界形状函数。这种形式简化了实现过程。该公式基于通用拉普拉斯方程[53]。形状函数方法非常灵活，已在多种应用中得到使用[54]、[55]、[56]。考虑没有源项的稳态条件下的标量场方程：
(8) ?·(??(x)) = 0，
其中?是标量场变量，? = [?/?x, ?/?y]?是微分算子。将微分算子转换到缩放边界坐标系后得到：
(9) ? = L??/?x + L??/?y = b?(η)?/?ξ + 1/ξ × b?(η)?/?η，
其中L? = [10]?，L? = [0, 1]?，b?(η) = 1/J_p(L?yη(η), η ? L?xη(η)，b?(η) = 1/J_p(L?xη(η) ? L?yη(η))。

通过将方程(8)的弱形式与任意权重函数w相乘并应用高斯散度定理，得到：
(10) ∫Ω χ(w) ?? dΩ = ∫Γ w??·n dΓ。

在缩放边界坐标系中，标量场变量和权重函数表示为：
(11) ?(ξ, η) = N(η)?h(ξ)，
w(ξ, η) = N(η)wh(ξ)，
其中?h(ξ)是满足控制微分方程的解析函数，wh(ξ)是任意解析函数。从弱形式方程(10)出发，并遵循传统的SBFEM方法，可以得到径向坐标系中的Euler-Cauchy方程[44]、[45]、[50]：
(12) E?ξ2?h(ξ) = 0，
ξ ξ + [E? + E???E?]ξ?h(ξ) ? E??h(ξ) = 0，
其中系数矩阵E?、E?和E?分别为：
(13a) E? = ∫?11 B?(η)?B?(η)Jp(η)dη，
(13b) E? = ∫?11 B?(η)?B?(η)Jp(η)dη，
(13c) E? = ∫?11 B?(η)?B?(η)Jp(η)dη，
B?(η)和B?(η)定义为：
B?(η) = b?(η)N(η)，
B?(η) = b?(η)N(η)，η。
方程(12)连同节点通量方程qp = E?ξ?h, ξ + E?T?h|ξ = 1通过标准SBFEM程序转换为二阶微分方程：
(14) ξ {?h(ξ)qp(ξ)}，
ξ = ?Z{?h(ξ)qp(ξ)}，
其中哈密顿矩阵Zp为：
Z = [(E?)?1(E?)??(E?)?1(E?)??E??E?(E?)?1]。

未知的?h(ξ)可以通过Z的特征值分解得到：
(15) Zφ = φS，
其中S是对角特征值矩阵，φ是相应的特征向量。特征值按从负到正的顺序排列，φ中的特征向量也进行了相应的调整[50]。对于n边形，仅选择前n个特征值及其对应的特征向量[50]、[57]。

注1：保留S的负特征值以确保在缩放中心ξ=0时解仍然是有限的。对于有界域，用Λ表示这些特征值的对角矩阵，因子ξ?Λ确保了在ξ=0时解是良定义的。特征值存储在Λ ? S(n)中。相应的特征向量被分成两部分，形成两个n × n矩阵并存储为[Ψ ? ? φ(n)。Ψ表示对应于标量场?h(ξ)的特征向量的前半部分，Φ表示对应于节点通量qp的特征向量。现在解?h(ξ)定义为：
(16) ?h(ξ) = Ψξ?Λc，
其中c是在边界(ξ=1)使用?h=?b得到的积分常数。现在，方程(11)重新表述为：
(17) ?h(ξ, η) = N(η)Ψ(e)ξ?Λψ?1?b = N(η)?b，
其中Ψ(e)是通过从多边形中选择对应于线元素e的行构建的。下一节将详细解释缩放边界形状函数的实现。

3.3. SBFEM实现
考虑第2节中的离散公式(4)，并将标量场和权重函数重写为缩放边界形状函数的形式：
(18) uh(ξ, η) = N(ξ, η)ub，
wh(ξ, η) = N(ξ, η)wb，
其中ub是标量场的节点值，wb是节点处的任意权重值。使用算子?求形状函数的导数为：
(19) ?N(ξ, η) = Bη(η)Bξ(ξ) =：
B(ξ, η)，
其中：
(20a) Bη = ?b?(η)N(η)Ψ(e)Λ + b?(η)N(η), ηψ(e)，
(20b) Bξ = ξ?Λ?IΨ?1，
其中I是单位矩阵。将方程6、(18)和19代入并应用边界条件后，对于任意的wb，双线性形式简化如下：
(21)
Asupg=∑E∈Th{∫01∫?11B(ξ,η)??B(ξ,η)Jp(η)ξdηdξub+∫01∫?11N(ξ,η)?vB(ξ,η)Jp(η)ξdηdξub+∫01∫?11τEB(ξ,η)?v?vB(ξ,η)Jp(η)ξdηdξub}=[Kd+Kc+Ks]ub，
其中Kd对应于扩散项的双线性项，Kc是对流项的双线性项，Ks是由于SUPG产生的双线性项，表示为
(22)
Kd=∑E∈ThΨ??XdΨ?1,
Kc=∑E∈ThΨ??XcΨ?1,
Ks=∑E∈ThΨ??XsΨ?1，
其中
(23)
Xd=∫01ξ?ΛYdξ?Λ?Idξ,
Yd=∫?11Bη(η)??Bη(η)Jp(η)dη,
Xc=∫01ξ?ΛYcξ?Λdξ,
Yc=Ψ(e)?∫?11N(η)?vBη(η)Jp(η)dη,
Xs=∫01ξ?Λ?IYsξ?Λdξ,
Ys=∫?11Bη(η)?v?τEvBη(η)Jp(η)dη。
Yd、Yc和Ys是通过沿η方向的高斯积分确定的。至于Xd、Xc和Xs的值，它们是通过分离变量积分得到，每个项可以表示为
(24)
Xdij=Ydij?Λii?Λjj,
Xcij=Ycij?Λii?Λjj+1,
Xsij=Ysij?Λii?Λjj。
注2：每当特征值对满足?Λii?Λjj=0时，Xd和Xs的ij分量被视为零，即Xdij=0和Xsij=0。对于Xcij中的分母项?Λii?Λjj+1，它始终满足条件?Λii?Λjj+1≥1，因为特征值Λii和Λjj是负的。

方程1中的线性项可以用缩放后的边界坐标表示为：
(25)
F=∑E∈Th{∫ΩN(ξ,η)?fdΩ+∫ΩτEB(ξ,η)?v?fdΩ}=Fb+Fs，
其中Fb是源项，Fs是由于SUPG稳定化而产生的项。对于Fb和Fs，源项f使用最小二乘法插值到缩放后的边界坐标，表示为f(ξ,η)=∑j=1rξjfj(η)，并在径向进行解析积分，在周向上使用高斯积分，表示为
(26)
Fb=∑E∈Th∑j=1rΨ??∫ξξ?Λ+(j+1)I(Ψ(e)?∫ηN(η)?fj(η)Jp(η)dη)dξ=∑E∈Th∑j=1rΨ??[?Λ+(j+2)I]?1(Ψ(e)?∫ηN(η)?fj(η)Jp(η)dη),
Fs=∑E∈Th∑j=1rΨ??∫ξξ?Λ+jI(∫ηBη(η)Tv?τEfj(η)Jp(η)dη)dξ=∑E∈Th∑j=1rΨ??[?Λ+(j+1)I]?1(∫ηBη(η)?v?τEfj(η)Jp(η)dη)，
其中r是最小二乘近似的阶数。

4. 数值示例
在本节中，我们数值研究了SUPG稳定化缩放边界有限元方法的准确性和收敛性质。为此，研究了五种不同复杂度的示例。此外，还展示了不同网格离散化的鲁棒性和准确性，例如变形的4节点四边形元素、规则的Voronoi多边形、变形的凸任意多边形元素和四叉树网格。与其他网格不同，四叉树网格有一个显著的特点，即四叉树网格中的单元类型是有限的（如果使用2:1的四叉树网格，则有16种不同的单元类型），无需重新计算网格中每个单元的形状函数和导数、刚度和耦合矩阵。这些可以预先计算并缩放到每个单元的大小[47]，[58]，从而减少了计算开销。为了用任意多边形离散化域，使用了Matlab工具PolyMesher，该工具可在http://paulino.ce.gatech.edu/software.html获取。为了研究数值解与精确解之间的误差（在适用的情况下），分别选择了L2范数和H1半范数中的误差eh,0和eh,1作为评估指标。还研究了网格尺寸h对Péclet数的影响。图2展示了本研究中使用的代表性网格。在Matlab中开发了内部代码，所有模拟都在配备64GB RAM、Intel Core i9-10900X × 20处理器（3.7 GHz）的计算机上执行。开发的代码可以从github仓库下载：SBFEM_CD。
下载：下载高分辨率图像（399KB）
下载：下载全尺寸图像
图2. 板的示意图离散化：(a) 变形四边形元素（64个元素），(b) 任意Voronoi多边形（50个元素）和(c) 变形凸任意多边形元素（168个元素）。

4.1. 示例1（平滑解 - 修补测试）
考虑域Ω=[0,1]×[0,1]，并使用上述网格对其进行离散化。对于这种情况，我们考虑一个制造解，u(x,y)=sinπxsinπy。源项f被确定为使假设的解决方案满足控制微分方程。假设的解决方案在所有边界上都得到强制，并且设v=[1,1]。考虑了两个不同的ε值，即ε=10^(-8)，代表两个极端情况。图3和图4分别显示了误差eh,0、eh,1的收敛性和Péclet数随网格尺寸的变化。可以推断，所提出的框架能够获得准确的结果，并具有最佳的收敛速率。此外，可以看出结果对使用的网格类型和ε的选择不敏感。

注3：为了评估方程4中的项，可以使用经典的多边形有限元方法（PFEM）。在PFEM中，未知场使用基于多边形元素的质心坐标（如Wachspress插值、均值坐标和Laplace插值）来近似。这些基函数是有理函数，通常需要大量的积分点来进行双线性形式的准确数值积分，从而增加计算成本。例如，在图2(c)所示的网格上，使用Wachspress插值的传统PFEM实现大约需要2秒，而所提出的方法在0.1秒内完成相同的计算。
下载：下载高分辨率图像（851KB）
下载：下载全尺寸图像
图3. 示例1：不同空间离散化下L2范数中误差的收敛性：(a) 变形正方形；(b) 随机Voronoi；以及(c) 变形Voronoi。还展示了扩散常数的影响。
下载：下载高分辨率图像（525KB）
下载：下载全尺寸图像
图4. 示例1：不同空间离散化下Péclet数随网格尺寸的变化：(a) 变形正方形；(b) 随机Voronoi；以及(c) 变形Voronoi。还展示了扩散常数的影响。

4.2. 示例2（具有内部层的解决方案）
接下来，为了展示所提出框架的鲁棒性，考虑一个精确解
(27)
u(x,y)=x(1?x)(1?y)[0.5?0.5tanh(α?xβ)]
对于模型问题（参见方程1）。上述假设的解决方案表现出内部层，参数α和β控制内部层的位置和厚度。对于这个示例，考虑一个正方形域Ω=[0,1]×[0,1]，并在所有边界上都强制应用精确解决方案。通过将方程27代入方程1并假设ε=10^(-9)和v=(1,0)来计算源项f。研究了两种不同的情况：(a) 情况A：α=0.5，β=0.1 和 (b) 情况B：α=0.5，β=0.01。图5显示了两种情况下的解析解和数值解的等高线图。结果显示，所提出的框架与解析解非常吻合。此外，为了测试准确性，进行了收敛性研究，结果分别显示在图6中。可以推断，当前框架对使用的空间离散化类型不敏感。图7显示了网格尺寸h对Péclet数的影响。正如预期的那样，Péclet数随网格尺寸的减小而减小，且所有类型的离散化行为都相同。
下载：下载高分辨率图像（661KB）
下载：下载全尺寸图像
图5. 对于h=0.0076的变形Voronoi网格，得到的解析解和数值解的等高线图。图(a,c,e)和(b,d,f)分别对应于情况A和B。
下载：下载高分辨率图像（987KB）
下载：下载全尺寸图像
图6. 不同网格下L2范数中误差的收敛性作为网格尺寸h的函数：(a) 变形正方形；(b) 随机Voronoi 和 (c) 变形Voronoi。
下载：下载高分辨率图像（660KB）
下载：下载全尺寸图像
图7. 不同网格下网格尺寸h对Péclet数的影响：(a) 变形正方形；(b) 随机Voronoi 和 (c) 变形Voronoi。还展示了扩散系数的影响。

4.3. 示例3
在这种情况下，考虑v(x)=(?2,1)和ε=10^(-6)。源项f根据以下假设的精确解来确定
(28)
u(x,y)=16(x?x^2)(y?y^2)[0.5+π?1arctan(200(0.252?(x?0.5)^2?(y?0.5)^2))。
上述解决方案具有圆形内部层，所提出的框架对于不同网格获得的结果显示在图8中。可以看出，无论使用哪种网格来划分区域，内部层都得到了准确的捕捉。图中显示了未知场u在L2范数中的收敛性和Péclet数随网格尺寸的变化。与之前的示例类似，所提出的框架获得了准确的结果，并且在L2范数中以最佳的收敛速率收敛。此外，Péclet数随网格尺寸的增加而减小，O(h)。
下载：下载高分辨率图像（197KB）
下载：下载全尺寸图像
图8. 不同网格离散化下的u(x, y)空间分布：(a) 解析解；(b) 变形正方形；(c) Voronoi 和 (d) 变形Voronoi网格。对于所有情况，网格尺寸h=0.076。
下载：下载高分辨率图像（634KB）
下载：下载全尺寸图像
图9. 比较：(a) 误差eh,0和eh,1的收敛性；以及(b) 不同网格类型的Péclet数。

4.4. 示例4（边界层测试）
考虑一个域Ω=[0,1]×[0,1]，以及源项f。对于这个示例，我们考虑以下解决方案：
(29)
u(x,y)=sin(πx^2)+sin(πy^2)(1?sin(πx^2))+e^(-1/ε)?e^(-1/(1?x))(1?y)/ε^(-1/ε)。
通过将方程29代入方程1来计算源项f。对于数值研究，我们考虑v=(1,1)和ε=10^(-6)。对于这种参数选择，方程29给出的解决方案在域的右侧和顶部边界附近表现出边界层。对于这个示例，我们考虑两种离散化：(a) 用4节点双线性四边形元素的均匀离散化；(b) 基于四叉树分解的局部自适应网格。图10显示了用于这项研究的两种代表性局部自适应网格。注意，这里的自适应网格是使用解决方案生成的。当前框架的结果与弱Galerkin方法（参见[59]中的表5.4）进行了比较。当前工作的结果还与使用均匀细化网格的WG方法和基于WG方法的自适应细化网格进行了比较。对于当前框架，只考虑了自适应细化。注意，传统FEM无法处理带有悬挂节点的元素，除非单独处理它们。图11显示了两种离散化下L2范数中误差的收敛性。可以看出，所提出框架的结果与弱Galerkin方法的结果非常吻合（见表1）。图12显示了解的等高线图。图12(b)和(c)还展示了SUPG稳定化的效果。从图12(d)和(e)可以看出，通过使用SUPG公式可以缓解整个计算域中的数值振荡。
下载：下载高分辨率图像（300KB）
下载：下载全尺寸图像
图10. 示例4：边界层测试 – 本例中使用的四叉树网格的示意图。
下载：下载高分辨率图像（300KB）
下载：下载全尺寸图像
图11. Weak Galerkin（WG）方法与SBFEM之间L2范数中相对误差的比较。域使用均匀正方形元素离散化，并使用四叉树分解。
表1. 示例4：在均匀正方形网格和自适应正方形网格上，Weak Galerkin（WG）方法与SBFEM之间L2范数中相对误差的比较。
均匀正方形网格
自适应正方形网格
WG方法 [59]
SBFEM
为空单元
WG方法 [59]
SBFEM
DOFe^h,0
eh,0
DOFe^h,0
25
66.14
90
e^(-18.34
09
e^(-16)
67
1.52
66
e^(-12.13
94
e^(-11)
02
42.87
22
e^(-15.79
24
e^(-11)
85
18.91
66
e^(-29.18
96
e^(-24)
06
91.50
66
e^(-13.04
22
e^(-15)
75
55.64
35
e^(-26.43)
81
e^(-21)
63
84
88.82
57
e^(-27.85)
81
e^(-21)
97
07
3.75
81
e^(-23.54)
38
e^(-26)
55
36
5.60
63
e^(-26.84)
24
e^(-27)
21
87
2.57
67
e^(-21.84)
09
e^(-2)
下载：下载高分辨率图像（647KB）
下载：下载全尺寸图像
图12. u(x, y)的等高线图：(a) 解析解；(b-c) 未使用SUPG稳定化的数值解；以及(d-e) 使用SUPG稳定化的数值解和误差。

4.5. 示例5（带障碍物的通道中的流动）
作为最后一个数值示例，我们考虑一个污染物（例如化学染料）通过充满水的通道中的恒定不可压缩流进行传输。在通道的不同位置放置了两根过滤管道，并假定它们在接触表面时会立即吸收任何污染物。污染物传输在通道的水平截面上进行建模。计算域定义为：Ω=[0,5]×[?0.5,0.5]，原点位于左边界的中点。管道的截面假设为椭圆形，由以下方程给出：
(30)
(x?1.5)^2
0.52+y^2
0.162=1,
以及
(x?3.5)^2
0.22+y^2
0.352=1。
假设流速是均匀的，并且与水平轴对齐。我们考虑低扩散情况，ε=10^(-5)，v=(1,0)，f=0。在?Ω的顶部和底部边界以及两个椭圆管道的边界上施加了均匀的Dirichlet边界条件，表示污染物在接触管道时立即被吸收。在左边界施加了抛物线形的Dirichlet流入条件，u(x,y)=?10(y?0.5)(y+0.5)，而在右侧（流出）边界施加了均匀的Neumann边界条件。由于存在内部几何特征，计算域不能使用结构化四边形网格进行离散化。相反，使用随机Voronoi多边形和特征自适应四叉树网格的组合对域进行离散化。图中展示了本研究中使用的不同网格类型。图14(a)显示了使用随机Voronoi网格获得的流场，清楚地表明污染物沿着直线轨迹传输。当化学染料接触到管道表面时，会迅速被吸收，从而使得过滤器流出的污染物浓度可以忽略不计。图14(b)比较了在距离左侧边界2.5处的不同网格配置下，污染物u随管厚变化的情况。结果表明，不同网格之间的结果非常一致，没有观察到虚假的振荡现象。这种一致性证实了所提出公式的数值稳定性、鲁棒性和对网格的独立性。下载：下载高分辨率图片（664KB）下载：下载全尺寸图片

图13. 通道内的流动：(a) 几何形状；(b-d) 不同网格类型的示意图，包括四叉树网格、体适性方形网格和随机Voronoi网格。下载：下载高分辨率图片（361KB）下载：下载全尺寸图片

图14. 有障碍物的通道内流动：(a) 等高线图示的流场；(b) 在x=2.5处不同网格类型下流场u的变化情况。

5. 结论

在这项工作中，通过将SBFEM与SUPG稳定性技术相结合，提出了一个用于求解对流-扩散问题的数值框架。该方法利用SBFEM在处理任意形状元素方面的灵活性，并结合SUPG稳定性技术有效抑制了对流主导问题中常见的数值不稳定性。通过一系列数值示例研究了该框架的收敛性能。从所有数值示例中可以看出：(a) 该方法能够得到准确的结果；(b) 对用于空间离散化的网格类型不太敏感；(c) 在L2范数和H1半范数下都以最优速率收敛。本研究采用了均匀分布的网格，但对于具有内部/边界层的问题（如示例2-4所示），这可能不是理想的选择。由适当的误差指标驱动的局部自适应细化值得进一步研究。将这种方法扩展到非线性、时变问题，以及在二维和三维空间中采用自适应网格策略以提高计算性能，是未来研究的方向。