曲面上對流-擴散-反應方程的兩種穩定化混合有限元方法的數值比較?

金孟晴，馮新龍?，何銀年,2

(1.新疆大學 數學與系統科學學院，新疆 烏魯木齊 830017；2.西安交通大學 數學與統計學院，陜西 西安 710049)

0 引言

對流-擴散-反應方程是基本的動力學方程,它描述了三種不同的物理現象：對流、擴散和反應[1?2].對流擴散現象包括河流污染、空氣污染、污染物在核廢料污染中的分布等.隨著數學和工程研究的發展,人們需要使用更加復雜的數學模型來描述自然界中的一些物理現象,這就使得研究者們將注意力集中在曲面問題的建立與求解上.因此,曲面偏微分方程的重要性越來越突出.

我們知道求解曲面橢圓方程的有限元方法可以追溯到1988年[3],作者用三角剖分來離散曲面,此后相繼出現求解曲面方程的方法.在流體力學中,一些文獻介紹和總結了兩相不可壓流的數值模擬方法.在生物學領域,一些生物學家用曲面對流和擴散現象來描述生物膜上的物質轉移[4?5]、種群遷移以及細菌在生物膜上的聚集現象[6?7].在大多數實際應用中,都需要求解曲面對流-擴散-反應方程.然而,用解析的方法在曲面上求解此類方程是一個巨大的挑戰,尤其是對流占優問題的相關求解.盡管一些學者對對流占優問題的數值逼近做了大量的工作,但仍然存在一個問題,即如何找到一種精確、穩定、快速的數值方法對其進行模擬.

對流占優問題,使方程有了雙曲方程的性質,這讓求解變得更加困難,用一般的方法,比如標準的有限元方法、有限差分方法會導致所得到的解產生非物理震蕩.只有當網格剖分的足夠細時,解才會相對穩定.但對于高維問題,網格剖分較細會導致大的計算量和存儲空間.因此,我們希望設計一種算法,即使網格尺寸稍大,也能克服上述問題.在本文中,我們采用混合的一階形式來近似對流-擴散-反應方程，其中兩個變量使用有限元對(P1-P1),根據平面微分方程中運用的穩定化思想[3,8?13],文中使用混合穩定化方法消除對流占優情況所產生的非物理震蕩.最后進行了數值實驗,來驗證兩種算法的有效性并得到了相應的結果.

1 預備知識

在這一部分,我們將介紹一些曲面的預備知識,比如曲面算子的一些符號和相應的Sobolev空間以供后面使用.此外,還推導了曲面上對流-擴散-反應方程的混合弱形式.

1.1 曲面算子

假設Γ ?R3是一個開區域,Γ ?Ω 是一個緊致曲面,可以用一個函數的零水平集φ(x)∈C2(Γ)定義,使得

對于充分光滑的函數f:Γ →R 在曲面上某一點x ∈Γ 的切向梯度可以定義為

這里,P(x)=I-n(x)?n(x),?表示平面上標準的梯度算子,其中曲面法向n定義為

函數φ 滿足?φ(x)/=0.

同樣的,曲面上的散度算子F(x)=(F1(x),F2(x),F3(x))∈R3可以被定義為

自然的,Laplace-Beltrami 算子可以定義為

1.2 Sobolev 空間

假設Γ 是屬于C2的一個曲面,Lp(Γ) 是可測量函數f:Γ →R 的空間,則我們將Sobolev 空間Wr,p(Γ) 定義為

這將用于推導對流-擴散-反應方程的混合形式.

1.3 曲面對流-擴散-反應方程

我們所研究的問題是以下的對流-擴散-反應方程:在Γ上找到p,使得

其中:ε ＞0 是擴散系數,α ∈W1,∞(Γ)3是對流項且在Γ 上滿足??！う?0,μ＞0 是反應系數,f 是體積力且屬于L2(Γ).

引理1[14]假設α 在曲面Γ 上是給定的無散度速度場,即滿足??！う?0.此外,如果α 和μ 滿足下式,則方程(12) 有弱解

我們思考(12) 的變分形式:找到p ∈H1(Γ),使得

定理1[14]假設Γ ∈C2且p ∈H2(Γ),我們得到

算法一:定義擴散通量為vd:=-ε?Γp,因此在Γ 上(12) 變成

(16) 的弱解是:找到(vd,p)∈V×Q:=H(??！?Γ)×L2(Γ),使得

對所有的(w,q)∈V×Q.

算法二:我們定義總通量為v:=-ε?Γp+αp,因此在Γ 上(12) 變成

(18) 的弱解是:找到(v,p)∈V ×Q:=H(??！?Γ)×L2(Γ),使得

對所有的(w,q)∈V ×Q.

引理2[15]若α 和μ 滿足(13) 且方程(12) 有唯一的弱解,使得

其中:C?是一個正常數.

注1:使用Lax-Milgram 引理,可以證明(12)具有唯一的弱解p ∈H1(Γ).因此,問題(16)或(18)的任何一個解都是(12) 的弱解,反之亦然.

2 穩定混合有限元方法

在這個部分,將介紹一種穩定混合有限元方法.首先,我們用一致的分片三角形來逼近光滑曲面Γ.因此,由近似曲面Γh引入的幾何誤差一般為O(h2),其中h 是分片三角形Γh的網格尺度.

我們選擇一系列的點{x1,x2,···,xNp} 生成一個不重疊的三角形單元{e1,e2,···,eNe},由這些單元近似光滑曲面所組成的曲面Γh被定義為:

其中:Ne和Np分別表示單元和節點的個數,h=max{h1,h2,···,hNe} 是離散曲面網格的大小.對于離散曲面Γh上的任意點x,有光滑曲面Γ 上的點a(x) 與之對應,滿足以下的映射關系

其中:d(x) 是一個定向距離函數并且滿足|d(x)|=dist(x,Γ),?d(x)=n(a(x)) 和|?d(x)|=1,詳細介紹可參閱文獻[16].

引理3[17]光滑曲面Γ 的定向距離函數d(x) 滿足以下估計

其中:C?是一個正常數.

對于曲面Γ 上定義的函數，我們可以將它投影到近似曲面Γh上

定理2[3]假設v 是v?l∈H1(Γ) 在Γh上的投影,則有下列范數等價

對于上面定義的有限元空間,我們引入了基于曲面Γh的有限元空間,其中通量變量v 的離散子空間為:

標量變量p 的離散子空間為:

我們可以得到(19) 的離散變分形式如下:找到(vh,ph)∈Hh×Qh,使得

本文在(17) 式的離散變分形式中加入了穩定項,此時,求解的離散形式變為:找到(vdh,ph)∈Hh×Qh,使得

本文在(29) 式上加入了穩定項,此時,求解的離散形式變為:找到(vh,ph)∈Hh×Qh,使得

其中:δ ＞0 是正的常數.

此外,我們在離散空間中定義算法一和算法二的范數:

算法一范數:

定理3設A(·,·)Γh是(30) 中的雙線性形式.存在一個不依賴于ε,μ 和h 的正常數C,使得

其中:(wh,qh)∈Hh×Qh.

證明 用A(·,·)Γh的定義,Cauchy-Schwarz 和Young 不等式,我們得到

注2:離散形式(32) 的穩定性在[15]中已被證明.

3 數值算例

在這部分,我們展示了幾個數值實驗的結果,以說明所提方法的有效性.

例1在球面上考慮問題(12).

球面隱函數為:

我們選擇連續的精確解如下,右端項f 由精確解給出

這個實驗的目的是檢驗方法的有效性.我們將取擴散系數ε 和反應系數c 等于1.為了計算誤差階,我們取網格剖分和自由度分別為h=0.2,0.1,0.05,0.025和447,1 703,6 718,26 562,所得到的結果如圖1所示.從圖1中,我們可以看到變量的L2范數的收斂階是2階,我們分別運用中間變量的定義將p 的H1范數的收斂階提升到2階,所得到的結果如圖1所示,這些結果與我們所做的理論分析一致.

圖1 比較兩種方法中p 和v 的誤差階

例2在流行上,對方程(12),考慮擴散系數是ε=10?3的解間斷問題.

曲面隱函數表達式為:

在這個實驗中,我們所采用的精確解和對流項與例1相同,測試兩種混合形式的對流占優問題.我們發現當ε 較大時,對流-擴散-反應方程具有光滑的精確解.然而,當ε 取值較小時,精確解的函數值是具有大梯度變化的對流占優問題.

我們在三角曲面網格上求解這個問題,圖2是所得到的數值模擬結果.從圖2中可以發現當ε=10?3時,其解在心型區域的中心上有非物理震蕩現象.由于曲面的曲率不一樣,在曲率大的地方模擬的效果比曲率小的部分差一些,且誤差也較大.在這個數值實驗中,我們將δ 設為1,并且采用相同的自由度求解此問題.由實驗結果可知,我們的方法可以成功地捕獲解的大跳躍,模擬效果較好.對比這兩種方法可知算法二（運用全局中間變量）較好.

圖2 當ε=10?3時,心型區域上兩種混合形式的比較

4 結論

本文拓展并分析了兩種穩定的混合有限元方法來求解曲面對流-擴散-反應方程.在給出穩定性分析外,還給出了一些數值實例,驗證了已知理論的結果.采用這兩種方法不僅可以有效地避免非物理震蕩問題,而且提高了?Γp 的精度.在文中還可以探討高斯曲率、平均曲率等對其方程求解時的影響,今后將是我們的研究重點.