楊星月 楊榮奎 馮民富

摘要：本文針對對流擴散反應(yīng)方程提出了一個穩(wěn)定化混合有限元格式，該格式基于混合有限元法與最小二乘法的結(jié)合.在此格式中，由于最小二乘穩(wěn)定項的引入，有限元逼近空間的選取無需滿足經(jīng)典的Ladyzhenkaya-Babuska-Brezzi（LBB）穩(wěn)定性條件，從而對兩個變量的有限元逼近可以方便地使用等階有限元組合. 對于定常的對流擴散反應(yīng)方程，本文獲得了有限元的穩(wěn)定性，對誤差進行了估計， 并以數(shù)值算例驗證了理論分析和格式的有效性. 對于非定常的對流擴散反應(yīng)方程， 本文給出了有限元的誤差估計和數(shù)值算例.

關(guān)鍵詞：對流擴散反應(yīng)方程； 穩(wěn)定化方法； 混合有限元； LBB穩(wěn)定性條件

中圖分類號：O241.82???文獻標(biāo)識碼：A????DOI：10.19907/j.0490-6756.2023.051001

收稿日期： ?2022-05-17

基金項目： ?國家自然科學(xué)基金（11971337）

作者簡介： ??楊星月（1997-）， 女， 四川成都人， 碩士研究生， 主要研究方向為微分方程數(shù)值解. 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 引 言

對流擴散反應(yīng)方程可以刻畫大氣、水流中污染物濃度的分布，流體流動、傳熱以及溫度擴散等現(xiàn)象，是一類重要的數(shù)學(xué)模型.

對流擴散反應(yīng)方程的數(shù)值解法已有多種，如Raviart與Thomas最早提出的混合有限元方法 ?［1］ 將問題重寫為不同的一階系統(tǒng)， 并選取Raviart-Thomas-Nedelec有限元空間對變分格式進行離散化. 隨后， Babuska和Brezzi提出了混合有限元法的一般理論 ?［2，3］ . 與標(biāo)準(zhǔn)的有限元相比，混合有限元法引入了通量， 可以對未知函數(shù)的微分算子進行直接求解，而引入的通量也可以近似到與原始變量相同的精度 ?［4，5］ .

在使用混合有限元方法時，通常需要強制地要求有限元空間組合滿足離散的Layzhenskaya-Brezzi-Babuska（簡稱LBB）穩(wěn)定性條件 ?［4，6，7］ . 雖然已有多種有限元空間滿足LBB穩(wěn)定性條件，如Raviart-Thomas-Nedelec有限元空間 ?［1，8］ ， Brezzi-Douglas-Fortin-Marini有限元空間 ?［9］ ， Brezzi-Douglas-Marini有限元空間（2維） ?［5］ ， Chen-Douglas有限元空間 ?［10］ 等， 但也有部分被廣泛應(yīng)用的有限元空間不滿足LBB穩(wěn)定性條件.

在高維問題中，構(gòu)建滿足LBB條件的混合有限元空間并不簡單，計算復(fù)雜且通常還需要解決一個典型鞍點問題. 為處理這些問題，學(xué)術(shù)界提出了一些穩(wěn)定化的方法， 如Subgrid modeling法 ?［11-13］ ，流線擴散Petrov-Galerkin法（SUPG） ?［14，15］ 及最小二乘法 ?［16，17］ ，等. 其中，對于后者，Aziz等 ?［18］ 提出了橢圓方程的最小二乘有限元法的一般理論.該方法的優(yōu)點是可以將一個非自共軛的問題轉(zhuǎn)化為一個對稱正定問題，通過在混合有限元法中引進最小二乘穩(wěn)定項來解決被選有限元空間滿足LBB穩(wěn)定性條件的限制，從而使得有限元空間的選取更加靈活. 這一結(jié)果已被 Pehlivanovand等所證明 ?［16］ .

對于定常的反應(yīng)擴散方程的數(shù)值解法，已有多種混合有限元，如Fu等 ?［19］ 通過添加適當(dāng)?shù)臍埐铐椀綄ε蓟旌先跣问街袠?gòu)建了一種新的穩(wěn)定混合有限元法，并在加權(quán)范數(shù)意義下證明了格式的強制性和連續(xù)性；Masud等 ?［20］ 使用對偶混合有限元格式，但沒有對該方法的穩(wěn)定性與誤差估計進行分析，等.隨后，Barrenechea等 ?［21］ 在文獻［20］的基礎(chǔ)上通過添加恰當(dāng)?shù)姆€(wěn)定項得到了一種新的穩(wěn)定化格式，并分析了該格式的穩(wěn)定性與收斂性.

受到文獻［19-21］的啟發(fā)，本文構(gòu)造了一種新的穩(wěn)定化混合有限元，并將其用于求解對流擴散反應(yīng)方程. 該方法的主要思想是將混合有限元法與穩(wěn)定化法相結(jié)合，引入合適的最小二乘殘差項，使得所選混合有限元空間不必滿足LBB穩(wěn)定性條件，從而可以對兩個變量均采用廣泛使用的標(biāo)準(zhǔn)拉格朗日有限元. 對于定常對流擴散反應(yīng)方程，我們對該方法的穩(wěn)定性和誤差估計進行了分析. 對于非定常的對流擴散反應(yīng)方程，我們對時間項采用向后歐拉有限差分逼近，并通過數(shù)值模擬驗證了該穩(wěn)定化混合有限元方法的有效性.

后文的結(jié)構(gòu)如下.在第2節(jié)中，我們給出一些必要的預(yù)備知識，提出并分析新的穩(wěn)定化混合格式. 在第3節(jié)中，我們得到穩(wěn)定化格式的誤差估計. 在第4節(jié)中，我們將新的穩(wěn)定化方法拓展應(yīng)用于非定常的對流擴散反應(yīng)方程， 并給出誤差估計. 在第5節(jié)中，我們用兩個數(shù)值算例來驗證理論分析.最后，在第6節(jié)中我們將對主要結(jié)果進行總結(jié).

2 ?定常對流擴散反應(yīng)方程的穩(wěn)定化混合元

考慮以下具有齊次Dirichlet型邊界條件的定常對流擴散反應(yīng)方程：

-εΔp+a·?

p+up=f， ?in Ω ，

p=0， ?on ?Ω ?（1）

其中有界凸區(qū)域Ω ??R ??d（d=2，3） ，邊界 ?? Ω 是 ???Lipschitz 連續(xù)的， ε>0 為常數(shù)擴散系數(shù)， ?a∈W ?1，∞ ?（ Ω ） ?d ?為對流場，滿足 ?

·a=0， ?in Ω，常數(shù) μ>0 為反應(yīng)項系數(shù)， f 為已知源函數(shù). 我們引入獨立通量 v -ε?

p+ap ，并將方程重寫為如下的一階混合形式：

v+ε?

p-ap=0， ?in Ω ，

?

·v+up=f， ?in Ω，

p=0， ??on ?Ω ?（2）

為了得到問題的弱形式，我們引入Lebesgue可積函數(shù)空間 L p（ Ω ）（1≤p≤∞） ，其中對于 D ?Ω ?， ?L 2（D） ?上的內(nèi)積用 ??·，· ??D 表示. D 上的 m- 階Sobolev空間用 W ?m，p ??Ω ??表示，其上的范數(shù)與半范數(shù)分別表示為 ??· ???m，p，D ?和 ??· ???m，p，D ?. 當(dāng) p=2 時，記 H m D ?W ?m，2 ?D ?， ??· ???m，D = ?· ???m，p，D ?.當(dāng) ?D =Ω 時，下標(biāo)Ω將被省略. Sobolev空間的標(biāo)準(zhǔn)符號及其相應(yīng)的范數(shù)參見文獻［22］. 在后文中，不加特殊說明時， C 表示與 ε ， μ ， a 及 h 無關(guān)的正常數(shù).

記 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 均為充分光滑函數(shù).

對某個正常數(shù) C 0 ，滿足 μ- 1 2 ?

·a≥C 0≥0 ，問題存在唯一解 ?［23］ .

為引進適當(dāng)?shù)臍埐钭鳛榉€(wěn)定項，我們添加如下穩(wěn)定項：

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

p- 1 ε ap， 1 ε w-?

q+ 1 ε aq ??（5）

則問題的穩(wěn)定化混合弱形式為： 求 （v，p）∈V×Q ，使得

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

根據(jù) ?

·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）

設(shè) ??T h ???h>0 ?是區(qū)域Ω的擬一致三角形剖分，單元 T 的直徑為 h T ?diam （T） ， 且

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

對任意整數(shù) k k≥1 ?，通量 v 的協(xié)調(diào)有限元逼近空間為：

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

T∈T h}，

其中 P k T ?表示在 T 上不超過 k 次的Lagrange多項式.定義標(biāo)量變量 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）

且投影算子滿足如下的穩(wěn)定性與逼近性質(zhì).

引理2.1 ????［22］ ?存在一個常數(shù) C ，使得有如下誤差估計成立：

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

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

基于上述有限元空間的定義，穩(wěn)定化混合問題的離散形式為： 求 ?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）

在正式對離散穩(wěn)定化混合形式的穩(wěn)定性與收斂性進行詳細分析前，我們先引入與雙線性形式 ?B（ ·，· ， ·，· ） ?相對應(yīng)的網(wǎng)格范數(shù)

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

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

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

定理2.2 ???對于式中的雙線性形式 B ·，· ?， 存在一個常數(shù) 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）

證明 ??首先令測試函數(shù) ?w h，q h = v h，p h ?. 根據(jù)雙線性形式 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）

接下來，取測試函數(shù) ?w h，q h ?= ?w h，0 ?.同樣，根據(jù) B ?·，· ， ·，· ??的定義式，結(jié)合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中的結(jié)論，令 w ?⌒ ??h v h- Π ?h ap h ?，并設(shè) ?w h，q h ?= ?w ?⌒ ??h，0 ?.利用 ?Π ?h 的正交性，結(jié)合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）相加并定義正常數(shù) τ≤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 ???新的穩(wěn)定化混合有限元是相容的.因而它還滿足以下的Galerkin正交關(guān)系：

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

3 誤差估計

為給出離散穩(wěn)定化混合形式的誤差估計，我們引入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）的精確解. 則存在一個常數(shù) 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）

根據(jù)范數(shù) ???·，· ???h 的定義，運用三角不等式并結(jié)合引理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））.

根據(jù)式（18），結(jié)合穩(wěn)定項的相容性可得

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 應(yīng)用引理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）

根據(jù)文獻［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）

結(jié)合式（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）

根據(jù)以上分析，不妨令

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

結(jié)合三角不等式可知， ??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）

根據(jù)定義 ?θ v，θ p = ?v h-I hv ， p h- I ?hp ??，結(jié)合式（24）～（26），應(yīng)用定理2.1中的結(jié)論式（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 ?非定常對流擴散反應(yīng)方程的穩(wěn)定化混合元

考慮以下具有齊次Dirichlet型邊界條件的非定常對流擴散反應(yīng)方程：

? 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 ?是時間區(qū)間， p 0（·） 是給定的初值. 類似地，我們引入通量 v -ε?

p+ap ，并將非定常對流擴散反應(yīng)方程重寫為混合形式

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）

方程的穩(wěn)定化混合弱形式為：

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 為實數(shù).定義向下取整運算， m=n ?表示， m 取小于等于 n 的最大整數(shù). 為了得到格式的全離散有限元格式，我們首先取時間步長 Δt=t n-t ?n-1 ?， ?n=1，2，…，N ，其中的 N=T/Δt ?是正常數(shù).在時刻 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（·） 的近似，上標(biāo) n+1 表示一個確定的時刻， f ?n+1 =f（·，t ?n+1 ） .

記 ?L r 0，T;W ?m，p ??Ω ???為帶時間的Banach空間，其范數(shù)的定義如下.

定義4.1 ???對于連續(xù)的時間 t 及Banach空間 L r 0，T;W ?m，p ??Ω ???， ?p∈L r 0，T;W ?m，p ??Ω ???， ?r∈ 1，+∞ ?，定義范數(shù)

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

當(dāng) r=∞ 時，取標(biāo)準(zhǔn)定義下的 L ∞（0，T;W ?m，p （ Ω ）） 范數(shù).

假設(shè) ?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）

根據(jù)式（34），類似于與定理3.2可得如下結(jié)論.

引理4.2 ???令 （v，p） 與 ?v ?^ ?h，p ?^ ?h ?分別是式（32）與式（34）的解. ?t∈I ，設(shè) 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 ?.結(jié)合式（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）

選取測試函數(shù) ?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 .

進而，結(jié)合范數(shù) ???·，· ???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）的右邊項進行分析.首先，經(jīng)計算并化簡可得

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）

因此，根據(jù)式（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 ，并關(guān)于 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.

最后，結(jié)合離散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）

再結(jié)合式（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 數(shù)值算例

為簡潔起見，記 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 .

網(wǎng)格尺寸取 N=h ?-1 ?.

例5.1 ???對定常的對流擴散反應(yīng)方程，選取空間區(qū)間 ?Ω ?= ?0，1 × 0，1 ?，并令 μ=0 ， a= ［1，2］ ?T . 方程的精確解為

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

在數(shù)值計算中，源函數(shù) f 由如上的精確解確定.在表1中我們分別給出了選取 P 2-P 2 元（即 k=2 ），擴散系數(shù) ε= 10 ??-3 ?與 ε= 10 ??-5 ?時的誤差及收斂精度. 可以看到：新的穩(wěn)定化混合有限元格式在 H 1 范數(shù)下變量 p 的誤差收斂率達到 O h k ?階，與第3節(jié)中的理論分析結(jié)果符合. 此外，在 L 2 -范數(shù)下，變量 p 的誤差收斂率達到 O h ?k+1 ??階.數(shù)值結(jié)果還顯示；當(dāng)擴散項系數(shù)較小時，該方法仍然有效，能夠克服對流占優(yōu).

例5.2 ???對于非定常的對流擴散反應(yīng)方程，選取空間區(qū)間Ω= ?0，1 × 0，1 ?，時間區(qū)間 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 .

在我們的數(shù)值格式中，源函數(shù) f 可由如上的精確解確定.在表2中，我們分別給出了選取 P 1-P 1 元（即 k=1 ），取時間步長 Δt=h 2 ，擴散系數(shù)為 ε= 10 ??-3 ?與 ε= 10 ??-5 ?時的誤差及收斂精度. 可以看到：穩(wěn)定化混合有限元格式在 H 1 -范數(shù)下變量 p 的誤差收斂率達到 O h k ?階，與第4節(jié)中的理論分析結(jié)果符合. 此外，在 L 2 -范數(shù)下，變量 p 的誤差收斂率達到 O h ?k+1 ??階.因此，該方法是有效的.

6 結(jié) 論

本文對于對流擴散反應(yīng)方程提出并分析了一種新的穩(wěn)定化混合有限元，通過引入最小二乘穩(wěn)定項，該格式解決了在選取混合有限元空間時受LBB穩(wěn)定性條件限制的問題. 該穩(wěn)定化方法也可應(yīng)用于非定常的對流擴散反應(yīng)方程.數(shù)值算例說明了方法的有效性與可靠性.

參考文獻：

［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.