楊星月 楊榮奎 馮民富

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

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

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

收稿日期： ?2022-05-17

基金項(xiàng)目： ?國(guó)家自然科學(xué)基金（11971337）

作者簡(jiǎn)介： ??楊星月（1997-）， 女， 四川成都人， 碩士研究生， 主要研究方向?yàn)槲⒎址匠虜?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 引 言

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

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

在使用混合有限元方法時(shí)，通常需要強(qiáng)制地要求有限元空間組合滿足離散的Layzhenskaya-Brezzi-Babuska（簡(jiǎn)稱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)定性條件.

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

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

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

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

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

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

-εΔp+a·?

p+up=f， ?in Ω ，

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

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

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

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

v+ε?

p-ap=0， ?in Ω ，

?

·v+up=f， ?in Ω，

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

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

記 V=L 2 ??Ω ???d ， Q H 1 0 ?Ω ??.混合問(wèn)題的弱形式為： 求 （v，p）∈V×Q ，滿足

（v，w）+ε（?

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

-（v，?

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

上式可以進(jìn)一步寫為： 求 （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ù).

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

·a≥C 0≥0 ，問(wèn)題存在唯一解 ?［23］ .

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

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

p- 1 ε ap， 1 ε w-?

q+ 1 ε aq ??（5）

則問(wèn)題的穩(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 } .

對(duì)任意整數(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 上不超過(guò) k 次的Lagrange多項(xiàng)式.定義標(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 ，即對(duì) 任意v∈ L 2 （ Ω ） ?d ?，都能找到一個(gè) ?Π ?h v ∈H h 滿足如下正交性 ?［22］ ：

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

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

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

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

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

基于上述有限元空間的定義，穩(wěn)定化混合問(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）

在正式對(duì)離散穩(wěn)定化混合形式的穩(wěn)定性與收斂性進(jìn)行詳細(xì)分析前，我們先引入與雙線性形式 ?B（ ·，· ， ·，· ） ?相對(duì)應(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 ???對(duì)于式中的雙線性形式 B ·，· ?， 存在一個(gè)常數(shù) C ，使得對(duì)任意 ?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）

證明 ??首先令測(cè)試函數(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 .

對(duì)上式使用Cauchy-Schwarz與Young不等式，化簡(jiǎn)后可得

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）

接下來(lái)，取測(cè)試函數(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）

又因?yàn)?/p>

‖ 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 誤差估計(jì)

為給出離散穩(wěn)定化混合形式的誤差估計(jì)，我們引入Scott-Zhang插值算子.

引理3.1 ????［24］ ?Scott-Zhang插值算子 I h：H 1 ??Ω ???d→ ?H h 與 ?I ?h：Q→ ??Q 0 h 滿足如下誤差估計(jì)：

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 是離散問(wèn)題的解， ?（v，p）∈H ?k+1 ?（ Ω ） ?d×［H ?k+1 ??Ω ?∩H 1 0 ?Ω ?］ ， ?（v，p） 是混合問(wèn)題（5）的精確解. 則存在一個(gè)常數(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 的定義，運(yùn)用三角不等式并結(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 .

進(jìn)而有

‖（η 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） 是問(wèn)題（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)定項(xiàng)的相容性可得

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 ，對(duì)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ù)文獻(xiàn)［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）

下面我們將對(duì)（19）式中的每一項(xiàng)的誤差進(jìn)行分析.首先，有

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 .

對(duì)上式化簡(jiǎn)可得

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）

我們先估計(jì)式（24）中的第一項(xiàng) R 1 .對(duì) 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）

我們?cè)俟烙?jì)中的第二項(xiàng) R 2 .化簡(jiǎn)后有

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）

最后，我們估計(jì)中的第三項(xiàng) 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 ?非定常對(duì)流擴(kuò)散反應(yīng)方程的穩(wěn)定化混合元

考慮以下具有齊次Dirichlet型邊界條件的非定常對(duì)流擴(kuò)散反應(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 ?是時(shí)間區(qū)間， p 0（·） 是給定的初值. 類似地，我們引入通量 v -ε?

p+ap ，并將非定常對(duì)流擴(kuò)散反應(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í)數(shù).定義向下取整運(yùn)算， m=n ?表示， m 取小于等于 n 的最大整數(shù). 為了得到格式的全離散有限元格式，我們首先取時(shí)間步長(zhǎng) Δt=t n-t ?n-1 ?， ?n=1，2，…，N ，其中的 N=T/Δt ?是正常數(shù).在時(shí)刻 t=t n 處，我們對(duì) ??p ?t ?使用向后歐拉有限 差分逼近，并用 ?v n h，p n h ??v h ·，t n ，p h ·，t n ??來(lái)近似 ??v h，p h ?，則問(wèn)題的全離散格式為：

求 ?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 表示一個(gè)確定的時(shí)刻， f ?n+1 =f（·，t ?n+1 ） .

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

定義4.1 ???對(duì)于連續(xù)的時(shí)間 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=∞ 時(shí)，取標(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 ?，則有如下估計(jì)成立：

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

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

下文中我們僅對(duì)逼近誤差 （v ?^ ?h-v h，p ?^ ?h-p h） 進(jìn)行估計(jì).

引理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 ?，則有如下估計(jì)：

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

證明 ?在式中取時(shí)間 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）

選取測(cè)試函數(shù) ?w h，q h = ???n+1 ?v，???n+1 ?p ?. 我們對(duì)上式中的每一項(xiàng)進(jìn)行分析.

對(duì)于等式左邊的第一項(xiàng)，我們有如下估計(jì)：

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

進(jìn)而，結(jié)合范數(shù) ???·，· ???h 的定義式與投影算子 ?Π ?h 的正交性（8）可知式（38）的左邊項(xiàng)有如下估計(jì)：

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）

下面我們對(duì)式（38）的右邊項(xiàng)進(jìn)行分析.首先，經(jīng)計(jì)算并化簡(jiǎn)可得

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）兩邊同時(shí)乘以 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（ Ω ）） .

進(jìn)一步化簡(jiǎn)可得

‖???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） ，則如下誤差估計(jì)成立：

‖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ù)值算例

為簡(jiǎn)潔起見，記 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 ???對(duì)定常的對(duì)流擴(kuò)散反應(yīng)方程，選取空間區(qū)間 ?Ω ?= ?0，1 × 0，1 ?，并令 μ=0 ， a= ［1，2］ ?T . 方程的精確解為

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

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

例5.2 ???對(duì)于非定常的對(duì)流擴(kuò)散反應(yīng)方程，選取空間區(qū)間Ω= ?0，1 × 0，1 ?，時(shí)間區(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 ），取時(shí)間步長(zhǎng) Δt=h 2 ，擴(kuò)散系數(shù)為 ε= 10 ??-3 ?與 ε= 10 ??-5 ?時(shí)的誤差及收斂精度. 可以看到：穩(wěn)定化混合有限元格式在 H 1 -范數(shù)下變量 p 的誤差收斂率達(dá)到 O h k ?階，與第4節(jié)中的理論分析結(jié)果符合. 此外，在 L 2 -范數(shù)下，變量 p 的誤差收斂率達(dá)到 O h ?k+1 ??階.因此，該方法是有效的.

6 結(jié) 論

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

參考文獻(xiàn)：

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