对流占优扩散方程广泛应用于河流污染, 空气污染, 组织渗透, 热盐环境等领域. 因此建立高效准确的计算方法求解对流占优扩散方程具有重要意义. 本文考虑如下对流占优扩散方程

其中 $ u(\boldsymbol{X},t) $ 表示物质浓度, $ T $ 为总时间, $ \varepsilon>0 $ 为扩散系数. $ \Omega\subseteq \boldsymbol{R}^d(d=2,3) $ 是一个有界多边形区域, 并且边界 $ \partial \Omega $ 是 Lipschitz 连续. $ \boldsymbol{\beta} $ 是速度矢量场满足无散度条件 $ \nabla \cdot \boldsymbol{\beta} = 0 $.

当扩散系数 $ \varepsilon $ 远小于对流系数 $ \boldsymbol{\beta} $ 时, 标准的 Galerkin 有限元和有限差分法在计算过程中会出现尖锐的锋面和过度的非物理振荡现象. 因此人们提出许多方法来减轻数值震荡现象, 如间断 Galerkin 法^[1-3], 特征有限元法^[4-6], 特征有限差分法^[9-11], 变分多尺度法^[12-14]和流线迎风 Petrov-Galerkin (SUPG)法^[11,12,13].

稳定化有限元方法既具有迎风方法的长处, 又能抑制非物理震荡现象, 具有较好的稳定性^[15]. Cengizci 等人^[13]将流线迎风 Petrov-Galerkin (SUPG) 稳定化有限元法结合激波捕捉机制求解对流占优扩散方程,并将该方法扩展到三维^[19]; Xie 等^[9]将变分多尺度稳定化方法与虚单元方法结合求解对流占优扩散问题, 证明了所提方法的有效性; Zhang等^[15]结合移动网格和变分多尺度有限元方法求解对流占优扩散方程, 提高了计算精度并减少数值震荡; Key 等^[11]提出了有限元和等几何稳定化方法求解对流扩散反应方程; Zhang 等^[18]提出一种自适应变分多尺度无单元 Galerkin 方法求解对流占优扩散方程. 容易看出, 变分多尺度方法在处理对流占优问题时表现出有效性和稳定性. 但它需要引进新的未知变量, 从而导致求解规模比传统有限元方法增大了许多, 使得该算法并不实用. 因此 Song 等^[19]提出了基于两局部高斯积分的数值方法来求解定常对流占优扩散方程. 这种方法利用局部单元上两局部高斯积分的残差来代替基于投影变分多尺度理论框架下的稳定化项, 其优点就是在不引入新的未知变量和不增加额外存储空间的情况下, 保持了经典变分多尺度方法的稳定性和处理对流占优问题的有效性.

本文进一步扩展了文献 [19] 方法的适用范围. 具体而言本文构造了一种变分多尺度稳定化有限元全离散加权 $ \theta $ 格式求解非定常对流占优扩散方程. 空间方向采用变分多尺度有限元方法进行离散, 利用局部单元上两局部高斯积分的残差来代替基于投影变分多尺度理论框架下的稳定化项构造全离散加权 $ \theta $ 格式. 同时分析了数值格式的误差估计, 得到 $ L^2 $ 范数意义下时空最优误差估计式. 最后通过数值算例验证该方法的有效性.

本文考虑二维情况. 设 $ W=H_{0}^{1}\left (\Omega \right) $, 方程 (1.1) 的变分问题为: 求 $ u\in W $ 满足

其中 $ A(u,v) = \varepsilon(\nabla u,\nabla v)+(\boldsymbol{\beta}\cdot \nabla u,v) $. 设 $ K_h=\{K\} $ 为 $ \Omega $ 的拟一致三角形剖分, 具有网格参数 $ h_i $ ($ K_h $ 中每个单元的最大直径), 并取 $ h=max\{h_i\} $. 选取如下有限元空间

其中 $ P_0 (K) $ 为 $ K $ 上的分段常数, $ P_1 (K) $ 为 $ K $ 上的线性多项式, $ B (K) =span\{b_K\} $ 为单元 $ K $ 上的函数空间. 对于每个 $ K\in K_h $ 有如下定义

其中 $ \lambda_{1} (x) $, $ \lambda_{2} (x) $ 和 $ \lambda_{3} (x) $ 表示三角单元 $ K $ 的面积坐标. 由面积坐标性质可知 $ \lambda_{i}\in H^1 (K) $ 且$ b_K\in H_0^1 (K) $. 考虑经典变分多尺度方法^[24]: $ \forall v_h\in V_h^b $, 求 $ u_h\in V_h^b $ 满足

其中 $ \delta $ 是取决于 $ h $ 的非负函数. $ P_H $ 为 $ \nabla u_h\rightarrow\boldsymbol{R}_0 $ 的 $ L^2 $ -正交投影.值得注意的是, 尽管方程 (2.2) 可以保持稳定性和有效性, 但在二维情况下, 它引入了两个额外的因变量, 增加了额外存储空间, 这种情况在三维情况下更为严重. 因此, 本文将给出一种新方法来有效减少这种额外的存储空间, 而不引入任何额外变量.

固定单元 $ K $, 设 $ u_h,v_h\in V_h^b(K)=P_1(K)\oplus B(K) $. 将其唯一分解为

记 $ K $ 的重心为 $ \bar q_K $, 本文给出满足一点高斯积分的新投影$ \varPi:\nabla V_h^b\rightarrow \boldsymbol{R}_0 $, 满足^[19]$ \varPi\nabla v_h|_K=\nabla v_h (\bar q_K). $ 又因为 $ u_b\in B(K) $, 在 $ \partial K $ 上为零^[20,21], 由散度定理

接着我们可以得到

其中 $ \iint_K $ 表示在 $ K $ 上的积分, $ \iint_{K, 1} $ 为 $ K $ 上的近似高斯积分, 它对一次多项式是精确的. 则等式 (2.3) 可以表示为 $ (\nabla u_h-\varPi\nabla u_h, \varPi \nabla v_h) =0. $ 因此方程 (2.2) 可以等价描述为: $ \forall v_h\in V_h^b $, 求 $ u_h\in V_h^b $ 满足

其中 $ A_h (u_h, v_h) = A(u_h,v_h) +G (u_h, v_h), $ 稳定化项 $ G (u_h, v_h) = (\delta\nabla u_h, \nabla v_h) - (\delta\varPi\nabla u_h, \varPi\nabla v_h) $ 具体表达式为: $ \forall u_h, v_h\in V_h^b $$ G (u_h, v_h) =\delta\sum_{K\in K_h} (\iint_{K, k}\nabla u_h\nabla v_h{\rm d}x{\rm d}y-\iint_{K, 1}\nabla u_h\nabla v_h{\rm d}x{\rm d}y), $其中 $ \iint_{K, i}g (x, y) $ 表示在 $ K $ 上的近似高斯积分, 对于 $ i (i=1\text{ 或 }k) $ 次多项式是精确成立的. 常数 $ \delta $ 作为 $ O (h) $ 的尺度, 用来稳定对流项. 这一稳定化方法避免了额外变量的引入, 显著降低了计算成本, 同时保持了传统变分多尺度方法的稳定性和有效性.

定理 2.1 对 $ \forall u_h,v_h\in V_h^b $, 存在与 $ h $ 无关的常数 $ C, c>0 $, 使得双线性形式 $ A_h(\cdot,\cdot) $ 是连续且强制的, 即

证 利用正交投影性质我们可以先对 $ G (u_h, v_h) $ 进行如下变形

根据正交投影性质我们可以得到恒等式 $ \lVert \nabla v_h\rVert_0^2=\lVert (I-\varPi) \nabla v_h\rVert_0^2 + \lVert \varPi \nabla v_h\rVert_0^2 $, 由 Cauchy-Schwarz 不等式可得

故算子 $ A_h(\cdot, \cdot) $ 是有界的. 接着考虑 $ A_h(u_h, u_h) $, 由分部积分有

由Poincaré不等式可得

将方程 (1.1) 中的时间 $ (0, T] $ 一致剖分步长为 $ \Delta t=t_n-t_{n-1} $ 的时间层单元, $ t^n=n\Delta t $, $ n = 1, 2, \dots, N $, $ N=T/\Delta t $. 我们采用有限差分法对方程 (2.4) 的时间尺度进行离散, 令 $ u_h^n $ 为$ u $ 在 $ t=t_n $ 处的近似值, 则方程 (2.4) 的全离散格式为

我们定义

其中 $ \theta $ 为具体格式的参数. 特别地, 当 $ \theta =1 $ 时该格式为向后欧拉格式; 当 $ \theta=\dfrac{1}{2} $ 时该格式为 Crank-Nicolson 格式.

接下来我们将对上述设计的对流占优扩散方程的数值格式进行理论分析. 在进行理论分析之前我们首先给出如下投影算子定义.

定义 2.1 对任意 $ u\in W $, 定义投影算子 $ \pi_h:W\rightarrow V_h^b $ 使

引理 2.2 对任意 $ u\in W $, 上述 $ \pi_hu $ 存在且唯一.

证 由引理 2.1 可知 $ A_h $ 在 $ V_h^b $ 上有界且强制. 设 $ F_u(v_h):=A(u,v_h) $ 为 $ V_h^b $ 上的有界线性泛函, 由 Lax-Milgram 定理可得 $ \pi_hu $ 存在且唯一.

引理 2.3 设 $ u\in W\cap H^{k+1}(\Omega) $, $ \pi_hu\in V_h^b $ 为定义 2.1 所给的, 则存在与 $ h $ 无关的常数 $ C>0 $ 使得

证 $\textbf{(1) $ H^1 $ 估计.}$ 由引理 2.1 可知, 离散双线性型 $ A_h(\cdot,\cdot) $ 在 $ V_h^b $ 上连续且强制. 因而对任意 $ w_h\in V_h^b $,

由 $ A_h $ 的连续性,

注意到 $ \|(I-\Pi)\nabla v_h\|_0\le \|\nabla v_h\|_0\le C\|v_h\|_1 $, 于是得到

为进一步利用稳定化项的结构, 记 $ Q(w_h):=\|u-w_h\|_1+\delta\,\|(I-\Pi)\nabla w_h\|_0, \ w_h=\tilde w_h+w_b, $ 其中 $ \tilde w_h\in P_1 $, $ w_b\in B $ 为泡空间部分, 并且这种分解在逐单元意义下唯一. 由 $ \varPi(\nabla \tilde w_h)=\nabla \tilde w_h,\ \varPi(\nabla w_b)=\nabla w_b(\bar q_K)=\frac{1}{9}(\nabla \lambda_1+\nabla \lambda_2+\nabla \lambda_3)=0, $ 可得

由 Poincaré 不等式可得

将 (2.12)-(2.13) 代入 $ Q(w_h) $, 并用三角不等式,

取 $ \delta\ge C_b $ 即得

这表明在该情况下去除泡部分不会增大 $ Q $. 同时因为 $ P_1\subset V_h^b $, 于是

据此我们可在 (2.11) 式中选取 $ w_h=I_h u\in P_1\subset V_h^b $. 注意 $ \nabla(I_hu) $ 在每个单元上为常向量, 从而 $ (I-\Pi)\nabla(I_hu)=0. $ 配合标准插值误差估计 $ \|u-I_hu\|_1\le C\,h^k\,\|u\|_{k+1}, $ 代入 (2.11)

$\textbf{(2) $ L^2 $ 估计.}$ 令 $ e:=u-\pi_h u $, 定义如下连续对偶问题: 求 $ \phi\in H_0^1(\Omega)\cap H^2(\Omega) $ 使得

并满足正则性 $ \|\phi\|_{2}\le C\|e\|_0 $ (参见文献 [22,定理 10.3.11]. 于是

由投影定义 $ A_h(\pi_h u,v_h)=A(u,v_h) $, 取 $ v_h=I_h\phi $ 且 $ A_h=A+G $, 有

又因为 $ \nabla I_h\phi $ 在每个单元上为常向量, $ \Pi\nabla I_h\phi=\nabla I_h\phi $, 故 $ G(\pi_h u, I_h\phi)=\delta\big((I-\varPi)\nabla\pi_h u,(I-\Pi)\nabla I_h\phi\big)=0, $ 得到 $ A(e,I_h\phi)=0 $. 从而

即 $ \lVert e\rVert_0\le C h\,\lVert e\rVert_1 $. 结合上一步的 $ H^1 $ 估计 $ \lVert e\rVert_1=\lVert u-\pi_h u\rVert_1\le C h^k\lVert u\rVert_{k+1} $, 有

定理 2.1 设方程 (1.1) 的解满足 $ u^n\in (W\cap H^{k+1}) $, 时间离散参数 $ \theta\in[\frac{1}{2}, 1] $, 则对于方程 (2.7) 的解 $ u_h^n $, 存在与 $ \Delta t $ 和 $ h $ 无关的常数 $ C $, 使得

证 设 $ u^n-u_h^n= (u^n-\pi_hu^n) - (u_h^n-\pi_hu^n) =\eta^n-\xi^n $. 我们可以得到如下误差方程

注意到 $ u_h^{n+1}-u_h^n=u^{n+1}-u^n- (\eta^{n+1}-\eta^n) + (\xi^{n+1}-\xi^n) $, 则我们有

取 $ v=\xi^{n+\theta} $, 带回方程 (2.16) 可得

上式两端同乘 $ \Delta t $, 从 $ i=0 $ 到 $ i=n-1 $ 累加可得

将上式记为 $ i_1+i_2=I_1+I_2 $. 由 $ (a, b) =\frac{1}{2} (\lVert a \rVert_0^2+\lVert b\rVert_0^2-\lVert a-b\rVert_0^2) $, 可得

注意到 $ \xi^0=0 $, 取 $ \theta\in[\frac{1}{2}, 1] $ 我们有

由引理 2.1 我们可以得到 $ i_2 $ 估计式

此外由 Cauchy-Schwarz 不等式成立

由 Young 不等式我们可以得到 $ I_2 $ 估计式

在对 $ I_1 $ 进行估计前我们先给出 $ u_t^{i+1} $ 和 $ u^{i+1} $ 的积分型余项的泰勒公式

因此由 Young 不等式我们可以得到 $ I_1 $ 估计式

结合方程 (2.18)-(2.22) 我们可以得到

根据离散形式的 Gronwall 引理可得

由引理 2.3 可得 $ \eta_t $ 的估计式

结合引理 2.3 和方程 (2.24) 和 (2.25), 利用结合三角不等式

定理2.1 可以得到证明.

我们选取如下精确解

考虑区域 $ \Omega= [0,1]^2 $, $ T=1 $, $ \Delta t=0.0001 $. $ f $ 是由精确解 $ u $ 确定的, $ \delta = 0.1h $. 考虑对流系数 $ \boldsymbol{\beta} = (1, 1) ^T $, 扩散系数 $ \varepsilon = 10^{-5} $. 表 1 和表 2 分别给出了向后欧拉和 Crank-Nicolson 格式在 $ P_{1b} $ 有限元空间采用变分多尺度方法得到的误差和空间收敛阶.

二维数值算例表明, 在一致网格进行等距的时间步长离散后可以获得理想结果, 两种格式在 $ L^2 $ 范数 $ H^1 $ 范数意义下收敛阶分别可以达到 2 阶和 1 阶.

此外取定参数 $ T=1 $, 扩散系数 $ \varepsilon = 10^{-5} $, 考虑对流系数 $ \boldsymbol{\beta} = (1, 1) ^T $, 网格剖分 $ h=1/512 $, 在表 3 给出了向后欧拉和 Crank-Nicolson 格式在 $ P_{1b} $ 有限元空间采用变分多尺度方法得到的误差和时间收敛阶.

从表 3 可以看出向后欧拉格式时间收敛阶可以达到 1 阶精度, 而 Crank-Nicolson 格式时间收敛阶可以达到 2 阶精度, 这与理论分析是一致的. 并且 Crank-Nicolson 格式可以比向后欧拉格式得到更小的误差.

另外我们在图 1 给出 $ h=1/32 $ 时, 向后欧拉和 Crank-Nicolson 格式的数值图像. 从图 1 我们可以看出这两种格式未出现数值振荡现象, 体现了这两种格式求解对流占优扩散方程的有效性.

高斯旋转问题是一个经典的对流占优扩散基准问题^[23,24]. 我们选定区域 $ \Omega =\{ (x, y) \in R^2:x^2+y^2\leq 1 \} $, $ \boldsymbol{\beta}= (-y, x) ^T $, $ f=0 $. 初始高斯条件以 $ (0.3, 0.3) $ 为中心, $ u_0\left ( x, y \right) =e^{-10\left[ \left ( x-0.3 \right) ^2+\left ( y-0.3 \right) ^2 \right]} $, 计算参数选择如下

图 2 给出了初始数值解图像, 图 3 给出了 $ T=2\pi $ 时数值图像.

从图像可以看出向后欧拉格式在 $ T=2\pi $ 时等值线向外围扩散, 出现了数值震荡现象. 而 Crank-Nicolson 格式同心圆结构保持完整, 与初始时刻几何特征高度一致, 说明 Crank-Nicolson 格式在旋转流场问题中可以有效减少数值震荡现象.

考虑在区域 $ \Omega = (0, 1) \times (0, 1) $ 上定义的源项 $ f=1 $, $ \boldsymbol{\beta} = (1, 0) ^T $, 边界条件 $ u=0 $, 扩散系数 $ \varepsilon=10^{-3} $, $ T=1 $, $ \Delta t=0.001 $. 我们将网格剖分成 $ 20\times20 $, 数值结果如图 4 所示. 从图像我们可以看出向后欧拉格式和 Crank-Nicolson 格式在 $ y=0 $ 和 $ y=1 $ 处出现抛物线层, 在 $ x=1 $ 处出现指数层, 并且存在 $ 45^{\circ} $ 斜率.

考虑区域 $ \Omega = (0, 1) \times (0, 1) $, $ \boldsymbol{\beta}= (1, 1) ^T $, 边界处 $ u=g $ 满足以下关系式

计算参数选取为: $ T = 2 $, $ \Delta t=0.01 $, $ \varepsilon =10^{-3} $. 在这种情况下, 解从 $ (0.3, 0) $ 开始形成一个与对流方向一致的内层, 并在 $ x = 1 $ 处显示一个指数边界层. 图 5 给出了向后欧拉和 Crank-Nicolson 格式的数值图像. 从图像可以看出向后欧拉格式和 Crank-Nicolson 格式都能很好的体现对流问题的内部层和指数外层.

在这个例子中, 我们考虑三维情况精确解 $ u (x,y,z,t) = e^t\sin (\pi x) \sin (\pi y) \sin (\pi z) $, $ f $ 是由精确解 $ u $ 决定的, 考虑区域 $ \Omega =\left[ 0, 1 \right] ^3 $. 计算参数选择如下

$ T=0.1, \Delta t = 0.001, \boldsymbol{\beta}= (1, 1, 1) ^T, \varepsilon=10^{-3}, \delta = 0.1h. $

表 4 给出了在时间步长 $ \Delta t $ 和空间步长 $ h $ 满足关系 $ \Delta t=h^2 $ 时向后欧拉在 $ P_{1b} $ 有限元空间采用变分多尺度方法得到的误差和收敛阶. 表 5 给出了在时间步长 $ \Delta t $ 和空间步长 $ h $ 满足关系 $ \Delta t=h $ 时 Crank-Nicolson 格式在 $ P_{1b} $ 有限元空间采用变分多尺度方法得到的误差和收敛阶.

对于向后欧拉方法, $ L^2 $ 范数误差的收敛阶趋近于理论 2 阶收敛, $ H^1 $ 范数误差收敛阶略高于 1 阶, 与理论分析结果基本一致. Crank-Nicolson 方法在 $ L^2 $ 范数和 $ H^1 $ 范数意义下空间收敛阶分别可以达到 2 阶和 1 阶, 这与理论分析结果一致.

在计算效率方面, Crank-Nicolson 方法在 $ 1/h=40 $ 时仅需要 239 秒就可以达到与向后欧拉相同的精度, 远低于向后欧拉的 7192 秒. 这是因为向后欧拉格式时间收敛阶是 1 阶精度, 而 Crank-Nicolson 格式可以达到 2 阶精度. 同时图 6 给出了在空间步长 $ h=1/16 $ 时 $ T=0.1 $ 的数值图像. 从图像可以看到该数值方法是可行且有效的.

本文提出了一种求解非定常对流占优扩散方程的稳定化有限元加权 $\theta$ 格式. 以两局部高斯积分替代 VMS 稳定项, 避免引入额外未知量与存储. 在此基础上构造投影算子 $\pi_h$, 证明了其存在唯一性与最优阶逼近, 得到全离散的 $L^2$ 时空误差估计. 二维与三维数值算例验证了理论结果, 并表明在相同精度下 Crank-Nicolson 格式具有更高的计算效率. 进一步我们可以将方法推广到非线性对流占优问题的数值求解中.