对流扩散反应方程的一个稳定化混合有限元

杨星月 杨荣奎 冯民富

摘要：本文针对对流扩散反应方程提出了一个稳定化混合有限元格式，该格式基于混合有限元法与最小二乘法的结合.在此格式中，由于最小二乘稳定项的引入，有限元逼近空间的选取无需满足经典的Ladyzhenkaya-Babuska-Brezzi（LBB）稳定性条件，从而对两个变量的有限元逼近可以方便地使用等阶有限元组合. 对于定常的对流扩散反应方程，本文获得了有限元的稳定性，对误差进行了估计， 并以数值算例验证了理论分析和格式的有效性. 对于非定常的对流扩散反应方程， 本文给出了有限元的误差估计和数值算例.

关键词：对流扩散反应方程； 稳定化方法； 混合有限元； LBB稳定性条件

中图分类号：O241.82   文献标识码：A    DOI：10.19907/j.0490-6756.2023.051001

收稿日期：  2022-05-17

基金项目：  国家自然科学基金（11971337）

作者简介：   杨星月（1997-）， 女， 四川成都人， 硕士研究生， 主要研究方向为微分方程数值解. E-mail： sherylyoung@qq.com

A stabilized mixed finite element  for convection-diffusion-reaction equations

YANG Xing-Yue， YANG Rong-Kui， FENG Min-Fu

（School of Mathematics， Sichuan University， Chengdu 610064， China）

In this paper， we propose a stabilized finite element for the convection-diffusion-reaction equations. This finite element combines the mixed finite element method with the least-squares method. Due to the introduced least-square stability term， the selection of finite element spaces does not need to satisfy the classical Ladyzhenkaya-Babuska-Brezzi （LBB） stability condition. As a result， the finite element approximation of the two variables can conveniently use equal order finite elements. For the steady convection-diffusion-reaction equations， we obtain the stability and give the error estimate of the finite element and exemplify the theoretical analysis and reliability by numerical experiments. For the unsteady convection-diffusion-reaction equations， we estimate the error and give an example for the finite element.

Convection-diffusion-reaction equation; Stabilized method; Mixed finite element; LBB stability condition

（2010 MSC 65M60）

1 引 言

對流扩散反应方程可以刻画大气、水流中污染物浓度的分布，流体流动、传热以及温度扩散等现象，是一类重要的数学模型.

对流扩散反应方程的数值解法已有多种，如Raviart与Thomas最早提出的混合有限元方法  ［1］ 将问题重写为不同的一阶系统， 并选取Raviart-Thomas-Nedelec有限元空间对变分格式进行离散化. 随后， Babuska和Brezzi提出了混合有限元法的一般理论  ［2，3］ . 与标准的有限元相比，混合有限元法引入了通量， 可以对未知函数的微分算子进行直接求解，而引入的通量也可以近似到与原始变量相同的精度  ［4，5］ .

在使用混合有限元方法时，通常需要强制地要求有限元空间组合满足离散的Layzhenskaya-Brezzi-Babuska（简称LBB）稳定性条件  ［4，6，7］ . 虽然已有多种有限元空间满足LBB稳定性条件，如Raviart-Thomas-Nedelec有限元空间  ［1，8］ ， Brezzi-Douglas-Fortin-Marini有限元空间  ［9］ ， Brezzi-Douglas-Marini有限元空间（2维）  ［5］ ， Chen-Douglas有限元空间  ［10］ 等， 但也有部分被广泛应用的有限元空间不满足LBB稳定性条件.

在高维问题中，构建满足LBB条件的混合有限元空间并不简单，计算复杂且通常还需要解决一个典型鞍点问题. 为处理这些问题，学术界提出了一些稳定化的方法， 如Subgrid modeling法  ［11-13］ ，流线扩散Petrov-Galerkin法（SUPG）  ［14，15］ 及最小二乘法  ［16，17］ ，等. 其中，对于后者，Aziz等  ［18］ 提出了椭圆方程的最小二乘有限元法的一般理论.该方法的优点是可以将一个非自共轭的问题转化为一个对称正定问题，通过在混合有限元法中引进最小二乘稳定项来解决被选有限元空间满足LBB稳定性条件的限制，从而使得有限元空间的选取更加灵活. 这一结果已被 Pehlivanovand等所证明  ［16］ .

对于定常的反应扩散方程的数值解法，已有多种混合有限元，如Fu等  ［19］ 通过添加适当的残差项到对偶混合弱形式中构建了一种新的稳定混合有限元法，并在加权范数意义下证明了格式的强制性和连续性；Masud等  ［20］ 使用对偶混合有限元格式，但没有对该方法的稳定性与误差估计进行分析，等.随后，Barrenechea等  ［21］ 在文献［20］的基础上通过添加恰当的稳定项得到了一种新的稳定化格式，并分析了该格式的稳定性与收敛性.

受到文献［19-21］的启发，本文构造了一种新的稳定化混合有限元，并将其用于求解对流扩散反应方程. 该方法的主要思想是将混合有限元法与稳定化法相结合，引入合适的最小二乘残差项，使得所选混合有限元空间不必满足LBB稳定性条件，从而可以对两个变量均采用广泛使用的标准拉格朗日有限元. 对于定常对流扩散反应方程，我们对该方法的稳定性和误差估计进行了分析. 对于非定常的对流扩散反应方程，我们对时间项采用向后欧拉有限差分逼近，并通过数值模拟验证了该稳定化混合有限元方法的有效性.

后文的结构如下.在第2节中，我们给出一些必要的预备知识，提出并分析新的稳定化混合格式. 在第3节中，我们得到稳定化格式的误差估计. 在第4节中，我们将新的稳定化方法拓展应用于非定常的对流扩散反应方程， 并给出误差估计. 在第5节中，我们用两个数值算例来验证理论分析.最后，在第6节中我们将对主要结果进行总结.

2  定常对流扩散反应方程的稳定化混合元

考虑以下具有齐次Dirichlet型边界条件的定常对流扩散反应方程：

-εΔp+a·?

p+up=f，  in Ω ，

p=0，  on ?Ω  （1）

其中有界凸区域Ω   R   d（d=2，3） ，边界  ? Ω 是    Lipschitz 连续的， ε>0 为常数扩散系数，  a∈W  1，∞  （ Ω ）  d  为对流场，满足 ?

·a=0，  in Ω，常数 μ>0 为反应项系数， f 为已知源函数. 我们引入独立通量 v -ε?

p+ap ，并将方程重写为如下的一阶混合形式：

v+ε?

p-ap=0，  in Ω ，

?

·v+up=f，  in Ω，

p=0，   on ?Ω  （2）

为了得到问题的弱形式，我们引入Lebesgue可积函数空间 L p（ Ω ）（1≤p≤∞） ，其中对于 D  Ω  ，  L 2（D）  上的内积用   ·，·   D 表示. D 上的 m- 阶Sobolev空间用 W  m，p   Ω   表示，其上的范数与半范数分别表示为   ·    m，p，D  和   ·    m，p，D  . 当 p=2 时，记 H m D  W  m，2  D  ，   ·    m，D =  ·    m，p，D  .当  D =Ω 时，下标Ω将被省略. Sobolev空间的标准符号及其相应的范数参见文献［22］. 在后文中，不加特殊说明时， C 表示与 ε ， μ ， a 及 h 无关的正常数.

记 V=L 2   Ω    d ， Q H 1 0  Ω   .混合问题的弱形式为： 求 （v，p）∈V×Q ，满足

（v，w）+ε（?

p，w）-（ap，w）=0， w∈V，

-（v，?

q）+μ（p，q）=（f，q），  q∈Q   （3）

上式可以进一步写为： 求 （v，p）∈V×Q ，使得

1 ε （v，w）+（?

p，w）- 1 ε （ap，w）-（v，?

q）+

μ（p，q）=（f，q），  （w，q）∈V×Q  （4）

其中 a ， μ ， f 均为充分光滑函数.

对某个正常数 C 0 ，满足 μ- 1 2 ?

·a≥C 0≥0 ，问题存在唯一解  ［23］ .

为引进适当的残差作为稳定项，我们添加如下稳定项：

S （v，p），（w，q）     - ε 2   1 ε v+?

p- 1 ε ap， 1 ε w-?

q+ 1 ε aq   （5）

则问题的稳定化混合弱形式为： 求 （v，p）∈V×Q ，使得

B （v，p），（w，q） =（f，q），   （w，q）∈V×Q    （6）

根据 ?

·a=0 以及 v 的定义，整理可得双线性形式 B  ·，· ， ·，·   的定义为

B （v，p），（w，q）   1 2ε  v，w +

1 2  ?

p，w - 1 2ε  ap，w +μ（p，q）- 1 2ε （v，aq）+

1 2ε  ap，aq + ε 2  ?

p，?

q - 1 2  v，?

q   （7）

设   T h    h>0  是区域Ω的拟一致三角形剖分，单元 T 的直径为 h T  diam （T） ， 且

h  max {h T：   T∈T h } .

对任意整数 k k≥1  ，通量 v 的协调有限元逼近空间为：

H h { φ  h∈C 0   Ω     d：   φ  h   T∈P k  T   d，

T∈T h}，

其中 P k T  表示在 T 上不超过 k 次的Lagrange多项式.定义标量变量 p 的离散子空间 Q 0 h Q h∩H 1 0  Ω   ，其中

Q h：={q h∈C 0  Ω   ：  q h   T∈P k T  ，

T∈T h}.

定义 L 2 投影算子  Π  h：L 2   Ω    d→H h ，即对 任意v∈ L 2 （ Ω ）  d  ，都能找到一个  Π  h v ∈H h 满足如下正交性  ［22］ ：

v- Π  h（v），w h =0，  w h∈H h       （8）

且投影算子满足如下的稳定性与逼近性质.

引理2.1     ［22］  存在一个常数 C ，使得有如下误差估计成立：

Π  h（v）   0≤C  v   0，  v∈L 2 （ Ω ）  d，

v- Π  h（v）   0≤Ch  v   1，  v∈H 1 （ Ω ）  d.

基于上述有限元空间的定义，稳定化混合问题的离散形式为： 求  v h，p h ∈H h×Q 0 h ，使得

B  v h，p h ， w h，q h  =    f，q h ，   w h，q h ∈H h×Q 0 h    （9）

在正式对离散稳定化混合形式的稳定性与收敛性进行详细分析前，我们先引入与双线性形式  B（ ·，· ， ·，· ）  相对应的网格范数

（v，q）∈V×Q ，满足

（v，q）   2 h μ  q   2 0+ε  q   2 1+

1 ε   v- Π  h（aq）   2 0   （10）

定理2.2    对于式中的双线性形式 B ·，·  ， 存在一个常数 C ，使得对任意  v h，p h ∈H h×Q 0 h 有

sup    （w h，q h）∈ H  h×Q 0 h   B  v h，p h ， w h，q h      w h，q h    h ≥

C   v h，p h    h   （11）

证明   首先令测试函数  w h，q h = v h，p h  . 根据双线性形式 B  ·，· ， ·，·   的定义式可得

B  v h，p h ， v h，p h  = 1 2ε  v h，v h -    1 ε  ap h，v h +μ p h，p h + ε 2  ?

p h，?

p h +    1 2ε  ap h，ap h .

对上式使用Cauchy-Schwarz与Young不等式，化简后可得

B  v h，p h ， v h，p h  = 1 2ε ‖v n‖ 2 0-

1 ε  ap h，v h +μ‖p n‖ 2 0+ ε 2 |p n| 2 1+

1 2ε ‖p n‖ 2 0≥μ‖p n‖ 2 0+ ε 2 |p n| 2 0   （12）

接下来，取测试函数  w h，q h  =  w h，0  .同样，根据 B  ·，· ， ·，·   的定义式，结合Cauchy-Schwarz与Young不等式可得

B  v h，p h ， w h，0  = 1 2ε  v h-ap h，w h +

1 2  ?

p h，w h .

由引理2.1中的结论，令 w  ⌒   h v h- Π  h ap h  ，并设  w h，q h  =  w  ⌒   h，0  .利用  Π  h 的正交性，结合Cauchy-Schwarz与Young不等式可得

B  v h，p h ， w  ⌒   h，0  = 1 2ε （v h-ap h，v h-

Π  h（ap h））+ 1 2  ?

p h，v h- Π  h ap h  =

1 2ε ‖v h- Π （ap h）‖ 2 0+ 1 2  ?

p h，v h- Π  h ap h  ≥

1 4ε ‖v h- Π （ap h）‖ 2 0- ε 4 |p h| 2 1   （13）

将式（12）与式（13）相加并定义正常数 τ≤1 ，则有

B  v h，p h ， v h+τ w   ⌒   h，p h  ≥

μ‖p h‖ 2 0+ ε 2-τ  4 |p h| 2 1+ τ 4ε ‖v h-    Π （ap h）‖ 2 0≥ τ 4 ‖（v h，p h）‖ 2 h   （14）

又因为

‖ v h+τ w   ⌒   h，p h ‖ h≤‖（v h，p h）‖ h+

1+2τ ε  ‖v h- Π  h（ap h）‖ 0≤    1+ 1+2τ  ‖（v h，p h）‖ h，

由式（14）可得离散inf-sup条件

sup    （w h，q h）∈H h×Q 0 h   B  v h，p h ， w h，q h      w h，q h    h ≥

B  v h，p h ， v h+τ w   ⌒   h，p h      v h + τ w   ⌒   h，p h    h ≥

τ 4   1+ 1+2τ   ‖（v h，p h）‖ h.

定理证毕.

注1    新的稳定化混合有限元是相容的.因而它还满足以下的Galerkin正交关系：

B  v-v h，p-p h ， w h，q h  =0，     （w h，q h）∈H h×Q 0 h  （15）

3 误差估计

为给出离散稳定化混合形式的误差估计，我们引入Scott-Zhang插值算子.

引理3.1     ［24］  Scott-Zhang插值算子 I h：H 1   Ω    d→  H h 与  I  h：Q→   Q 0 h 满足如下误差估计：

p- I  hp   m≤Ch  k+1-m   p    k+1 ，

p∈H  k+1   Ω  ∩H 1 0  Ω .

v-I hv   m≤Ch  k+1-m   v    k+1 ，

v∈H  k+1    Ω    d.

定理3.2    令  v h，p h ∈H h×Q 0 h 是离散问题的解，  （v，p）∈H  k+1  （ Ω ）  d×［H  k+1   Ω  ∩H 1 0  Ω  ］ ，  （v，p） 是混合问题（5）的精确解. 则存在一个常数 C 使得

v-v h，p-p h    h≤   Ch k G 1  p    k+1 +G 2  v    k+1     （16）

成立，其中

G 1=C 1 1 +    a    0，∞ h ε   + μ   1 2  h，

G 2= C 1 ε  ， C 1= min     a    0，∞  μ   1 2   ， h  a    1，∞  μ   1 2    +ε   1 2  .

证明  记

（v-v h，p-p h）=（（v-I hv），（p- I  hp））-

（（v h-I hv），（p h- I  hp）） （η v，η p）-（θ v，θ p）   （17）

根据范数    ·，·    h 的定义，运用三角不等式并结合引理2.1可得

‖（η v，η p）‖ h≤ μ   1 2    η p   0+ε   1 2    η p   1+

1 ε   1 2     η v   0+C   a    0，∞  ε   1 2     η p   0 .

进而有

‖（η v，η p）‖ h≤

Ch k μ   1 2  h+   a    0，∞ h ε   1 2   +ε   1 2   ‖p‖  k+1 +

Ch  k+1  1 ε   1 2   ‖v‖  k+1    （18）

因 （v h，p h） 是问题（5）的解及

（v h，p h）=（v，p）+（θ v，θ p）-（η v，η p） ，

则 （θ v，θ p）∈H h×Q 0 h 满足如下方程：

B（（θ v，θ p），（w h，q h））=

（f，q）-B（（v，p），（w h，q h））+B（（η v，η p），

（w h，q h））.

根据式（18），结合稳定项的相容性可得

B  θ v，θ p ， w h，q h  =B  η v，η p ， w h，q h     （19）

用id表示恒等算子. 因 （w h，q h）∈H h×Q 0 h ，对id-Π  h 应用引理2.1可得

aq h- Π  h aq h    0≤C  a    0，∞   q h   0≤

C   a    0，∞  μ   1 2      w h，q h    h  （20）

根据文献［24］中的引理1.137（Bertoluzza），有

aq h- Π  h aq h    0≤Ch  a    1，∞   q h   0≤

C h  a    1，∞  μ   1 2      w h，q h    h  （21）

结合式（20）与式（21）可得

aq h- Π  h aq h    0≤C·

min     a    0，∞  μ   1 2   ， h  a    1，∞  μ   1 2       w h，q h    h  （22）

根据以上分析，不妨令

C 1  min     a    0，∞  μ   1 2   ， h  a    1，∞  μ   1 2    +ε   1 2  .

结合三角不等式可知，   w h，q h ∈H h×Q 0 h ，有

‖w h-aq h‖ 0≤‖w h-∏ h（aq h）‖ 0+

‖aq h-∏ h（aq h）‖ 0≤CC 1‖（w h，q h）‖ h  （23）

下面我们将对（19）式中的每一项的误差进行分析.首先，有

B  η  v ，η p ， w h，q h  = 1 2ε  η v-aη p，w h +

1 2  ?

η p，w h - 1 2ε  η v-aη p，aq h -

1 2  η v，?

q h +μ η p，q h + ε 2  ?

η p，?

q h .

对上式化简可得

B （η  v，η  p），（w  h，q  h） =

1 2   1 ε η  v- 1 ε aη  p+?

η  p，w  h-aq  h +

1 2  ε?

η  p-η  v-aη  p，?

q  h +μ η  p，q  h

∑ 3  i=1  R  i  （24）

我们先估计式（24）中的第一项 R 1 .对 R 1 使用Cauchy-Schwarz与Young不等式可得

R 1= 1 2   1 ε η v- 1 ε aη p+?

η p，w h-aq h ≤

1 ε η v- 1 ε aη p+?

η p  0‖w h-aq h‖ 0≤

CC 1h k  h ε   v    k+1 +    a    0，∞ h ε +1

p    k+1  ‖（w h，q h）‖ h   （25）

我们再估计中的第二项 R 2 .化简后有

R 2= 1 2  ε?

η p-aη p-η v，?

q h =

ε 2  ?

η p，?

q h - 1 2  aη p，?

q h - 1 2  η v，?

q h ≤

Ch k  ε   1 2  +   a    0，∞ h ε   1 2      p    k+1 +

h ε   1 2     v    k+1  ‖（w h，q h）‖ h  （26）

最后，我们估计中的第三项 R 3 .容易得到

R 3≤μ   1 2  ‖η p‖ 0‖（w h，q h）‖ h≤

Ch  k+1 μ   1 2  ‖p‖  k+1 ‖（w h，q h）‖ h  （27）

将式（25）～（27）代入式（24），整理后可得

B  η  v ，η p ， w h，q h  ≤

Ch k   C 1  a    0，∞ h ε +C 1+

μ   1 2  h   p    k+1  ‖（w h，q h）‖ h+

Ch k  C 1h ε  ‖v‖  k+1 ‖（w h，q h）‖ h  （28）

根据定义  θ v，θ p =  v h-I hv ， p h- I  hp   ，结合式（24）～（26），应用定理2.1中的结论式（11）可得

‖（θ v，θ p）‖ h≤

C  sup    （w h，q h）∈H h×Q 0 h   B  η v，η p ， w h，q h      w h，q h    h =

Ch k μ   1 2  h  p    k+1 +

C 1 1 +    a    0，∞ h ε  ‖p‖  k+1 +  C 1h ε   v    k+1     （29）

最后，由三角不等式可得

v-v h，p-p h    h≤   η v，η p    h+   θ v，θ p    h.

证毕.

4  非定常对流扩散反应方程的稳定化混合元

考虑以下具有齐次Dirichlet型边界条件的非定常对流扩散反应方程：

? tp-εΔp+a· p+up=f，  in Ω ×（0，t］，

p=0，  on  ? Ω ×（0，t］，

p（·，0）=p 0（·），  on Ω ，t=0    （30）

其中 T>0 ， I= 0，T  是时间区间， p 0（·） 是给定的初值. 类似地，我们引入通量 v -ε?

p+ap ，并将非定常对流扩散反应方程重写为混合形式

v+ε p-ap=0，  in Ω ×（0，t］，

? tp+ ·v+μp=f，  in Ω ×（0，t］，

p=0，  on  ? Ω ×（0，t］，

p（·，0）=p 0（·），  in Ω    （31）

方程的稳定化混合弱形式为：

t∈I ，求 （v，p）∈V×Q ，使得

（? tp，q）+B（（v，p），（w，q））=

（f，q）（w，q）∈v×Q，

p（·，0）=p 0（·）   （32）

其中 B（（v，p），（w，q）） 的定义与式相同.

令 m，n 为实数.定义向下取整运算， m=n  表示， m 取小于等于 n 的最大整数. 为了得到格式的全离散有限元格式，我们首先取时间步长 Δt=t n-t  n-1  ，  n=1，2，…，N ，其中的 N=T/Δt  是正常数.在时刻 t=t n 处，我们对  ?p ?t  使用向后欧拉有限 差分逼近，并用  v n h，p n h   v h ·，t n ，p h ·，t n   来近似   v h，p h  ，则问题的全离散格式为：

求  v  n+1  h，p  n+1  h ∈H h×Q 0 h ，使得

1 Δt （p  n+1  h-p n h，q h）+B（（v  n+1  h，p  n+1  h），（w h，q h））=

（f  n+1 ，q h），  （w h，q h）∈H h×Q 0 h，

p h（.，0）=p 0 h（.）， p 0 h∈Q 0 h     （33）

其中 p 0 h（·） 是 p 0（·） 的近似，上标 n+1 表示一个确定的时刻， f  n+1 =f（·，t  n+1 ） .

记  L r 0，T;W  m，p   Ω    为带时间的Banach空间，其范数的定义如下.

定义4.1    对于连续的时间 t 及Banach空间 L r 0，T;W  m，p   Ω    ，  p∈L r 0，T;W  m，p   Ω    ，  r∈ 1，+∞  ，定义范数

p     L  r（0，T;W   m，p （ Ω ）） =   ∫   T 0  p    r  m，p  d t      1 r  .

当 r=∞ 时，取标准定义下的 L ∞（0，T;W  m，p （ Ω ）） 范数.

假设  v  ^  h，p  ^  h ：  0，T  →H h×Q 0 h 满足

B   v-v  ^  h，p-p  ^  h ， w h，q h  =0，

w h，q h ∈H h×Q 0 h  （34）

根据式（34），类似于与定理3.2可得如下结论.

引理4.2    令 （v，p） 与  v  ^  h，p  ^  h  分别是式（32）与式（34）的解.  t∈I ，设 v ·，t ∈H  k+1  ，  p ·，t ∈ H 1 0  Ω  ∩H  k+1  ，则有如下估计成立：

v-v  ^  h，p-p  ^  h    h≤

Ch k G 1  v    k+1  + G 2  p    k+1    （35）

下文中我们仅对逼近误差 （v  ^  h-v h，p  ^  h-p h） 进行估计.

引理4.3    令  v  n+1  h，p  n+1  h  与  v  ^  h，p  ^  h  分别是式（33）与式（34）的解，且 p h（·，0）=p  ^  h ·，0  .若方程（32）的解 p∈［H 1 0，T;H  k+1   Ω   ×   H 2 0，T;L 2  Ω    ］，  v∈H 1 0，T;H  k+1    Ω    d  ，则有如下估计：

‖p  ^  h-p h‖  L ∞（0，T;L 2（ Ω ）） +   |‖（v  ^  h-v h，p  ^  h-p h）‖|  L 2（0，T;stab） ≤

Ch  k+1    p t    L 2（0，T;H  k+1 （ Ω ）） +

v t    L 2（0，T;H  k+1  （ Ω ）  d）  +   CΔt‖p  tt ‖  L 2（0，T;L 2（ Ω ））   （36）

其中

v，   p     2  L 2（0，T;stab）

∑  N-1   n=0  Δt       n+1  v，    n+1  p     2 h.

证明  在式中取时间 t=t  n+1  .记

v  n+1 -v  n+1  h  =  v  n+1 -v  ^   n+1  h +

v  ^   n+1  h-v  n+1  h  ζ  n+1  v+   n+1  v，

p  n+1 -p  n+1  h  =  p  n+1 -p  ^   n+1  h +

p  ^   n+1  h-p n h  ζ  n+1  p+   n+1  p  （37）

令式中的 （w，q）= w h，q h  .结合式（33）与式（34）可得：   w h，q h ∈H h×Q 0 h ，有

n+1  p-  n p Δt ，q h +B     n+1  v，   n+1  p ， w h，q h  =

- ? tp  n+1 - p  n+1 -p n Δt ，q h -  ζ  n+1  p-ζ n p Δt ，q h    （38）

选取测试函数  w h，q h =    n+1  v，   n+1  p  . 我们对上式中的每一项进行分析.

对于等式左边的第一项，我们有如下估计：

n+1  p-  n p Δt ，   n+1  p ≥ 1 2Δt       n+1  p   2-    n p   2 .

进而，结合范数    ·，·    h 的定义式与投影算子  Π  h 的正交性（8）可知式（38）的左边项有如下估计：

n+1  p-  n p Δt ，   n+1  p +

B     n+1  v，   n+1  p ，    n+1  v，   n+1  p  ≥

1 2Δt       n+1  p   2-    n p   2 +Ch  2  （39）

下面我们对式（38）的右边项进行分析.首先，经计算并化简可得

p   n+1 -p  n =  ∫    t   n+1    t  n p  t d t ≤ ∫    t   n+1    t  n  p  t  d t，

?  tp   n+1 - p   n+1 -p  n Δt  =  1 Δt  ∫    t   n+1    t  n  t-t  n p   tt  d t ≤

∫    t   n+1    t  n  p   tt   d t.

由Cauchy-Schwarz不等式可得

ζ   n+1  p-ζ  n p Δt ，    n+1  p ≤ 1 Δt  ∫    t   n+1    t  n  ζ   n+1   p，t  dt·     n+1  p ≤

1 Δt ‖ζ   p，t ‖  2     L 2 （t  n，t   n+1 ;L 2（ Ω ））   2+C‖    n+1  p‖ 2  （40）

? tp  n+1 - p  n+1 -p n Δt ，   n+1  p ≤

Δt  p  tt    2  L 2（t n，t  n+1 ;L 2（ Ω ）） +C     n+1  p   2  （41）

因此，根据式（39）～（41）有

1 2Δt       n+1  p   2-    n p   2 +C‖（   n+1  p，   n+1  p）‖ 2 h≤

1 Δt ‖ζ  p，t ‖  2   L 2（t n，t  n+1 ;L 2（ Ω ））   2+

Δt‖p  tt ‖   2    L 2（t n，t  n+1 ;L 2（ Ω ））   2+C‖   n+1  p‖ 2  （42）

由初始条件 p 0 h=p  ^  0 h ，即   0 p=0 ，在式（42）两边同时乘以 2Δt ，并关于 n 从0到 N-1 求和可得

∑  N-1   n=0         n+1  p   2-     n p   2 +

C∑  N-1   n=0  Δt       n+1  v，    n+1  p     2 h≤

C∑  N-1   n=0    ζ   p，t     2  L 2（t  n，t   n+1 ;L 2（ Ω ）） +C∑  N-1   n=0  Δt      n+1  p   2+

C∑  N-1   n=0    Δt   2  p   tt    2  L 2（t  n，t   n+1 ;L 2（ Ω ）） .

进一步化简可得

‖   N p‖ 2+C∑  N-1   n=0  Δt       n+1  v，    n+1  p     2 h≤

C‖ζ   p，t ‖    2    L 2（0，T;L 2（ Ω ）） +

C（Δt） 2‖p   tt ‖    2    L 2（0，T;L 2（ Ω ）） +

C∑  N-1   n=0  Δt      n+1  p   2.

最后，结合离散Gronwall引理  ［25］ 可得

‖   N p‖ 2+C∑  N-1   n=0  Δt       n+1  v，    n+1  p     2 h≤

C‖ζ   p，t ‖   2     L 2（0，T;L 2（ Ω ））   2+

C（Δt） 2‖p   tt ‖   2     L 2（0，T;L 2（ Ω ））   （43）

再结合式（43）与引理4.2，引理证毕.

定理4.4     令 （v，p） 与  v  n+1  h，p  n+1  h  分别是方程（32）与方程（33）的解.若引理4.2与引理4.3中的条件均成立，且 p∈L ∞（0，T;H 1 0（ Ω ）∩   H  k+1 （ Ω ））  ， v∈L ∞（0，T;H  k+1  （ Ω ）  d） ，则如下误差估计成立：

‖p-p h‖  L ∞（0，T;L 2（ Ω ）） +

‖（v-v h，p-p h）‖  L 2（0，T;stab） ≤

C（h k+Δt）   （44）

其中

v-v  h，p-p  h     2  L 2（0，T;stab）

∑  N-1   n=0  Δt   v   n+1 -v   n+1  h，p   n+1 -p   n+1  h     2 h.

5 数值算例

为简洁起见，记 e p=p-p h ，

E p=  max    1≤n≤N    p n-p n h   0 ，

E  ?

p =  max    1≤n≤N    ?

p n-?

p n h   0 .

网格尺寸取 N=h  -1  .

例5.1    对定常的对流扩散反应方程，选取空间区间  Ω  =  0，1 × 0，1  ，并令 μ=0 ， a= ［1，2］  T . 方程的精确解为

p x，y = sin  2 π x  sin  2 π y .

在数值计算中，源函数 f 由如上的精确解确定.在表1中我们分别给出了选取 P 2-P 2 元（即 k=2 ），扩散系数 ε= 10   -3  与 ε= 10   -5  时的误差及收敛精度. 可以看到：新的稳定化混合有限元格式在 H 1 范数下变量 p 的误差收敛率达到 O h k  阶，与第3节中的理论分析结果符合. 此外，在 L 2 -范数下，变量 p 的误差收敛率达到 O h  k+1   阶.数值结果还显示；当扩散项系数较小时，该方法仍然有效，能够克服对流占优.

例5.2    对于非定常的对流扩散反应方程，选取空间区间Ω=  0，1 × 0，1  ，时间区间 I= 0，1  ，即 T=1 .令 μ=1 ， a= ［y，-x］  T .方程的精确解为

p x，y =

1- cos t 100x 2 1-x 2 y 1-y  1-2y .

在我们的数值格式中，源函数 f 可由如上的精确解确定.在表2中，我们分别给出了选取 P 1-P 1 元（即 k=1 ），取时间步长 Δt=h 2 ，扩散系数为 ε= 10   -3  与 ε= 10   -5  时的误差及收敛精度. 可以看到：稳定化混合有限元格式在 H 1 -范数下变量 p 的误差收敛率达到 O h k  阶，与第4节中的理论分析结果符合. 此外，在 L 2 -范数下，变量 p 的误差收敛率达到 O h  k+1   阶.因此，该方法是有效的.

6 结 论

本文对于对流扩散反应方程提出并分析了一种新的稳定化混合有限元，通过引入最小二乘稳定项，该格式解决了在选取混合有限元空间时受LBB稳定性条件限制的问题. 该稳定化方法也可应用于非定常的对流扩散反应方程.数值算例说明了方法的有效性与可靠性.

参考文献：

［1］    Raviart P A， Thomas J M. Mathematical aspects of finite element methods ［M］. Berlin： Springer， 1977.

［2］  Russell  F T. Time stepping along characteristics with incomplete iteration for a Galerkin approximation of miscible displacement in porous media ［J］. SIAM J Numer Anal， 1985， 22： 970.

［3］  Demkowicz L， Brezzi F， Boffi D，  et   al . Mixed finiteelements， compatibility conditions， and applications ［M］. Berlin： Springer， 2008.

［4］  Brezzi F. On the existence， uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers ［J］. ESAIM-Math Model Num Anal， 1974， 8： 129.

［5］  Brezzi F， Douglas J， Marini L D. Two families of mixed finite elements for second order elliptic problems ［J］. Numer Math， 1985， 47： 217.

［6］  Ladyzhenskaya  O A. The mathematical theory of viscous incompressible flow ［M］. New York： Gordon and Breach， 1969.

［7］  Babuvka I， Rheinboldt W C. Error estimates for adaptive finite element computations ［J］. SIAM J Numer Anal， 1978， 15： 736.

［8］  Nedelec J C. Mixed finite elements in  R   3   ［J］. Num Math， 1980， 35： 315.

［9］  Brezzi  F， Fortin M， Marini L D，  et al . Efficient rectangular mixed finite elements in two and three space variables ［J］. ESAIM-Math Model Num Anal， 1987， 21： 581.

［10］  Chen  Z X， Douglas J. Prismatic mixed finite elements for second order elliptic problems ［J］. Calcolo， 1989， 26： 135.

［11］ Guermond J L. Subgrid stabilization of Galerkin approximations of monotone operators ［J］. Z Angew Math Mech， 1999， 79： 29.

［12］ Guermond J L. Stabilization of Galerkin approximations  of transport equations by subgrid modeling ［J］. ESAIM-Math Model Numer Anal， 1999， 33： 1293.

［13］ Layton W. A connection between subgrid scale eddy viscosity and mixed methods ［J］. Appl Math Comput， 2002， 133： 147

［14］ Brooks A N， Hughes T J R. Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations ［J］. Comput Method Appl M， 1982， 32： 199.

［15］ Johnson C， Navert U， Pitkaranta J. Finite element methods for linear hyperbolic problems ［J］. Comput Methods Appl Mech Eng， 1984， 45： 285.

［16］ Pehlivanovand  A L， Carey C F， Lazarov R D. Least-squares mixed finite elements for second-order elliptic problems ［J］. SIAM J Numer Anal， 1994， 31： 1368.

［17］ Cai Z Q， Lazarov R， Manteuel T A，  et al . First-order system least squares for second-order partial dierential equations： Part I ［J］. SIAM J Numer Anal， 1994， 31：1785.

［18］ Aziz  A K， Kellogg R B， Stephens A B. Least squares methods for elliptic systems ［J］. Math Comput， 1985， 44： 53.

［19］ Fu H F， Guo H， Hou J，  et al . A stabilized mixed finite element method for steady and unsteady reaction-diusion equations ［J］. Comput Method Appl M， 2016， 304： 102.

［20］ Masud A， Kwack J. A stabilized mixed finite element method for the first-order form of advection-diusion equation ［J］. Int J Numer Meth， 2008， 57： 1321.

［21］ Barrenechea  G R， Poza A H， Yorston H. A stabilised finite element method for the convection-diusion-reaction equation in mixed form ［J］. Comput Methods Appl Mech Eng， 2018， 339： 389.

［22］ Brenner S C， Scott L R. The mathematical theory of finite element methods ［M］. New York： Springer， 2008.

［23］ Ern A， Guermond J L. Theory and practice of finite elements ［M］. New York： Springer， 2004.

［24］ Burman  E， Fernández M A. Continuous interior penalty finite element method for the time-dependent Navier-Stokes equations： space discretization and convergence ［J］. Numer Math， 2007， 107： 39.

［25］ Heywood， John G， Rolf R. Finite element approximation of the nonstationary Navier-Stokes problem （I）， regularity of solutions and second-order error estimates for spatial discretization ［J］. SIAM J Numer Anal， 1982， 19： 275.