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

(1.新疆大學(xué) 數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，新疆 烏魯木齊 830017；2.西安交通大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，陜西 西安 710049)

0 引言

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

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

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

1 預(yù)備知識(shí)

在這一部分,我們將介紹一些曲面的預(yù)備知識(shí),比如曲面算子的一些符號(hào)和相應(yīng)的Sobolev空間以供后面使用.此外,還推導(dǎo)了曲面上對(duì)流-擴(kuò)散-反應(yīng)方程的混合弱形式.

1.1 曲面算子

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

對(duì)于充分光滑的函數(shù)f:Γ →R 在曲面上某一點(diǎn)x ∈Γ 的切向梯度可以定義為

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

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

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

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

1.2 Sobolev 空間

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

這將用于推導(dǎo)對(duì)流-擴(kuò)散-反應(yīng)方程的混合形式.

1.3 曲面對(duì)流-擴(kuò)散-反應(yīng)方程

我們所研究的問題是以下的對(duì)流-擴(kuò)散-反應(yīng)方程:在Γ上找到p,使得

其中:ε ＞0 是擴(kuò)散系數(shù),α ∈W1,∞(Γ)3是對(duì)流項(xiàng)且在Γ 上滿足??！う?0,μ＞0 是反應(yīng)系數(shù),f 是體積力且屬于L2(Γ).

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

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

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

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

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

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

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

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

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

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

其中:C?是一個(gè)正常數(shù).

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

2 穩(wěn)定混合有限元方法

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

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

其中:Ne和Np分別表示單元和節(jié)點(diǎn)的個(gè)數(shù),h=max{h1,h2,···,hNe} 是離散曲面網(wǎng)格的大小.對(duì)于離散曲面Γh上的任意點(diǎn)x,有光滑曲面Γ 上的點(diǎn)a(x) 與之對(duì)應(yīng),滿足以下的映射關(guān)系

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

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

其中:C?是一個(gè)正常數(shù).

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

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

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

標(biāo)量變量p 的離散子空間為:

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

本文在(17) 式的離散變分形式中加入了穩(wěn)定項(xiàng),此時(shí),求解的離散形式變?yōu)?找到(vdh,ph)∈Hh×Qh,使得

本文在(29) 式上加入了穩(wěn)定項(xiàng),此時(shí),求解的離散形式變?yōu)?找到(vh,ph)∈Hh×Qh,使得

其中:δ ＞0 是正的常數(shù).

此外,我們?cè)陔x散空間中定義算法一和算法二的范數(shù):

算法一范數(shù):

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

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

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

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

3 數(shù)值算例

在這部分,我們展示了幾個(gè)數(shù)值實(shí)驗(yàn)的結(jié)果,以說明所提方法的有效性.

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

球面隱函數(shù)為:

我們選擇連續(xù)的精確解如下,右端項(xiàng)f 由精確解給出

這個(gè)實(shí)驗(yàn)的目的是檢驗(yàn)方法的有效性.我們將取擴(kuò)散系數(shù)ε 和反應(yīng)系數(shù)c 等于1.為了計(jì)算誤差階,我們?nèi)【W(wǎng)格剖分和自由度分別為h=0.2,0.1,0.05,0.025和447,1 703,6 718,26 562,所得到的結(jié)果如圖1所示.從圖1中,我們可以看到變量的L2范數(shù)的收斂階是2階,我們分別運(yùn)用中間變量的定義將p 的H1范數(shù)的收斂階提升到2階,所得到的結(jié)果如圖1所示,這些結(jié)果與我們所做的理論分析一致.

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

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

曲面隱函數(shù)表達(dá)式為:

在這個(gè)實(shí)驗(yàn)中,我們所采用的精確解和對(duì)流項(xiàng)與例1相同,測(cè)試兩種混合形式的對(duì)流占優(yōu)問題.我們發(fā)現(xiàn)當(dāng)ε 較大時(shí),對(duì)流-擴(kuò)散-反應(yīng)方程具有光滑的精確解.然而,當(dāng)ε 取值較小時(shí),精確解的函數(shù)值是具有大梯度變化的對(duì)流占優(yōu)問題.

我們?cè)谌乔婢W(wǎng)格上求解這個(gè)問題,圖2是所得到的數(shù)值模擬結(jié)果.從圖2中可以發(fā)現(xiàn)當(dāng)ε=10?3時(shí),其解在心型區(qū)域的中心上有非物理震蕩現(xiàn)象.由于曲面的曲率不一樣,在曲率大的地方模擬的效果比曲率小的部分差一些,且誤差也較大.在這個(gè)數(shù)值實(shí)驗(yàn)中,我們將δ 設(shè)為1,并且采用相同的自由度求解此問題.由實(shí)驗(yàn)結(jié)果可知,我們的方法可以成功地捕獲解的大跳躍,模擬效果較好.對(duì)比這兩種方法可知算法二（運(yùn)用全局中間變量）較好.

圖2 當(dāng)ε=10?3時(shí),心型區(qū)域上兩種混合形式的比較

4 結(jié)論

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