李 西， 羅加福， 馮民富

(1.四川大學(xué)數(shù)學(xué)學(xué)院, 成都 610064； 2.成都體育學(xué)院， 成都 610041)

1 引 言

粘性不可壓縮Navier-Stokes方程(NS方程)的混合有限元離散通常會(huì)出現(xiàn)兩個(gè)難點(diǎn):需滿足經(jīng)典的inf-sup條件和需要克服由高雷諾數(shù)(Reynolds number,Re)導(dǎo)致的數(shù)值震蕩.

為了解決上述兩個(gè)問(wèn)題, 大量專家學(xué)者提出了不同的穩(wěn)定化方法.比如，為了克服對(duì)流-擴(kuò)散方程中對(duì)流占優(yōu)時(shí)所引起的數(shù)值震蕩,Brooks和Hughes[1]提出了Streamline Upwind Petrov-Galerkin方法(SUPG),并將其推廣到不可壓縮非定常NS方程.隨后, Hughes 和Franca等[2]發(fā)現(xiàn)SUPG法不僅可以保持格式的相容性,不損失解的逼近性,還能增強(qiáng)離散解的穩(wěn)定性——可利用SUPG 法避免inf-sup條件從而得到了速度和壓力的等階有限元插值.這種用來(lái)穩(wěn)定離散壓力的Petrov-Galerkin變分又被稱為Pressure Stabilized Petrov-Galerkin(PSPG). 雖然SUPG/PSPG既可以做到避免inf-sup條件,又可以得到對(duì)高Re依然有效的先驗(yàn)誤差估計(jì)，但需在經(jīng)典的Galerkin變分格式中引入最小二乘項(xiàng),導(dǎo)致穩(wěn)定項(xiàng)中含有非物理耦合項(xiàng),從而在用高階元近似時(shí)穩(wěn)定化格式中需要計(jì)算二階導(dǎo)數(shù)[3-4].隨后,Becker和Braack提出的Local Projection Stabilization方法(LPS)[5]很好地解決了上述不足.與SUPG相比, LPS中添加的穩(wěn)定項(xiàng)中不含有速度-壓力的耦合項(xiàng),并且不需要計(jì)算二階導(dǎo)數(shù)，也可以避免inf-sup條件和克服對(duì)流占優(yōu).正因如此,LPS從提出到現(xiàn)在得到了充分地發(fā)展.從文獻(xiàn)[5]中首次采用LPS來(lái)處理Stokes問(wèn)題之后，LPS被推廣到了運(yùn)輸方程[6],Oseen方程的低階離散[7],NS方程[8-9]以及其它一些相關(guān)工作[10-12].

在處理由高Re帶來(lái)的數(shù)值震蕩時(shí),Burman和Fernandez等[13]采用了Bertoluzza在文獻(xiàn)[14]中證明了的離散commutator性質(zhì),得到的誤差估計(jì)隨著流體Re的增大依然有效,即該數(shù)值格式可以克服由高Re帶來(lái)的數(shù)值震蕩.之后, Chen和Feng等[9]以及Frutos和Garcia等人[12,15]將上述技巧用在NS方程的LPS方法中,同樣得到了對(duì)高Re依然有效的誤差估計(jì).

本文在文獻(xiàn)[9,11-12,16]基礎(chǔ)上,提出了一種基于非線性項(xiàng)的局部投影穩(wěn)定格式，其中時(shí)間半隱格式中的非線性局部投影項(xiàng)等價(jià)于文獻(xiàn)[16]中的SUPG-type穩(wěn)定項(xiàng)和文獻(xiàn)[11]中的高階term-by-term穩(wěn)定項(xiàng)(因而本文提出的非線性局部投影可看做是該兩類穩(wěn)定方法的非線性推廣)，對(duì)速度-壓力我們采用等階元(Pk,Pk)進(jìn)行逼近,并添加局部壓力梯度投影穩(wěn)定化項(xiàng)克服inf-sup條件.在分析離散解的先驗(yàn)誤差估計(jì)時(shí),我們采用Bertoluzza[14]提出的離散commutator性質(zhì)，得到的誤差估計(jì)右端項(xiàng)的常系數(shù)中不含有粘性系數(shù)的倒數(shù)1/ν,使得當(dāng)流體的Re增大(即流體的粘性系數(shù)ν減小)時(shí)該誤差估計(jì)依然有效.

2 半離散穩(wěn)定化有限元格式

Hm(Ω)=Wm,2(Ω),

‖·‖m,M=‖·‖Hm(M),

‖·‖k,∞,M=‖·‖Wk,∞(M).

當(dāng)M=Ω時(shí)，我們省略下標(biāo)M,即‖·‖m=‖·‖Hm(Ω),‖·‖k,∞=‖·‖Wk,∞(Ω), 其中k=1,2是一個(gè)整數(shù). 我們分別用Lp,Hs和Wm,p來(lái)表示相應(yīng)于Lp,Hs和Wm,p的向量值Sobolev 空間,以及

設(shè)X為Sobolev空間, 定義映射φ(x,t):[0,T]→X,

并簡(jiǎn)記‖φ‖Lp(X)=‖φ‖Lp(0,T;X),p=2或∞.

取時(shí)間區(qū)間I=[0,T], 其中T為一個(gè)固定的正常數(shù). 非定常不可壓縮粘性流體的流動(dòng)由以下非定常Navier-Stokes方程表示:

(1)

其中u=u(x,t)∈Rd,p=p(x,t)∈R,f=f(x,t)∈Rd分別表示不可壓縮粘性流體流動(dòng)的速度場(chǎng), 壓力場(chǎng)和外力場(chǎng). 設(shè)特征長(zhǎng)度L和特征速度U均為單位1. 于是該系統(tǒng)的運(yùn)動(dòng)粘度系數(shù)ν=Re-1, 其中Re為雷諾數(shù).

求(u,p)∈X×Q, 滿足

(?tu,v)+b(u;u,v)+ν(?u,?v)-(p,?·v)=

(2)

Xh(M)={vh|M:vh∈Vh,vh=0 on ΩM}.

由文獻(xiàn)[17], LPS方法誤差分析的關(guān)鍵在于存在一個(gè)插值算子j滿足如下最優(yōu)近似性質(zhì):

假設(shè)2.1令以下局部inf-sup條件成立：

β>0

(3)

由該局部inf-sup條件可以推導(dǎo)出插值算子ju和jp具有如下的正交性和近似性[17]:

引理2.2令假設(shè)2.1成立.則存在兩個(gè)插值算子ju:X→Xh,jp:Q→Qh滿足以下的正交性和近似性[18]:

(u-juu,vh)=0,?u∈X,?vh∈Dh,

(p-jpp,qh)=0,?p∈Q,?ph∈Dh,

‖u-juu‖0+h|u-juu|1≤

Chr+1‖u‖r+1,?u∈X∩Hr+1(Ω)，

‖p-jpp‖0+h|p-jpp|1≤

Chr+1‖p‖r+1,?p∈Q∩Hr+1(Ω)，

‖u-juu‖0,∞+h|u-juu|1,∞≤

Ch‖u‖1,∞,?u∈X∩W1,∞(Ω)

(4)

下面我們給出NS方程的非線性局部投影穩(wěn)定的空間半離散格式.?t∈(0,T), 求(uh(t),ph(t))∈Xh×Qh, 滿足,?(vh,qh),

(?tuh,vh)+ν(?uh,?vh)+b(uh;uh,vh)-

(ph,?·vh)+Sconv(uh;uh,vh)=f,vh，

(qh,?·uh)+Spres(ph,qh)=0

(5)

其中

Sconv(uh;uh,vh)=

其中α1,M和α2,M分別為相應(yīng)于對(duì)流項(xiàng)和壓力梯度項(xiàng)的局部投影穩(wěn)定參數(shù), 并且滿足以下假設(shè).

α1,M=α1hm1,α2,M=α2hm2，

其中m1,m2為穩(wěn)定參數(shù)的階,其具體數(shù)值在后文的誤差分析中給出.

3 含穩(wěn)定項(xiàng)的全離散格式

接下來(lái)我們對(duì)空間半離散格式(5)做時(shí)間上的有限差分.一般說(shuō)來(lái),時(shí)間差分上的全隱格式為無(wú)條件穩(wěn)定,但每層時(shí)間上需解一個(gè)非線性方程組，全顯格式在計(jì)算模擬時(shí)具有優(yōu)勢(shì),但為了格式的穩(wěn)定性需對(duì)時(shí)間步長(zhǎng)Δt有諸多限制.一種通常的做法是線性項(xiàng)采用隱式格式,非線性項(xiàng)采用顯示格式. 下面我們給出(5)的不同時(shí)間差分格式.

(6)

其中

注意到該格式中的非線性局部投影項(xiàng)實(shí)為線性型:

此時(shí)從程序?qū)崿F(xiàn)時(shí)做數(shù)值積分的角度可知, 該項(xiàng)等價(jià)于以下兩種線性化對(duì)流項(xiàng)的局部投影：

(1) 高階term-by-term穩(wěn)定項(xiàng)[11]

(2) SUPG-type的局部投影穩(wěn)定項(xiàng)[16]

(κM((uM·?)uh),κM((uM·?)vh))M.

其中uM為流場(chǎng)速度u的局部投影平均

(7)

(8)

其中

下文我們將給出格式2的穩(wěn)定性和收斂性分析.在前面的討論中我們知道,格式1中的非線性局部投影穩(wěn)定項(xiàng)可等價(jià)于文獻(xiàn)[16]中提出的SUPG-type的穩(wěn)定項(xiàng)和文獻(xiàn)[11]中的高階term-by-term穩(wěn)定項(xiàng).為簡(jiǎn)潔起見(jiàn),我們略去格式3中離散解的穩(wěn)定性和收斂性分析.

4 穩(wěn)定性與收斂性分析

本節(jié)我們給出全離散格式(7)的穩(wěn)定性與收斂性分析.為此我們需要一些必要的假設(shè)和不等式, 并且為了記法的簡(jiǎn)潔,我們定義

?(vh,qh)∈Xh×Qh.

在算子Π的定義中，曾假設(shè)Π具有局部穩(wěn)定性及近似性[18]，即

‖Π(v)‖0,M≤C‖v‖0,M.

注1有好幾類算子滿足上式給出的最優(yōu)近似性和局部穩(wěn)定性,如文獻(xiàn)[19]中的Scott-Zhang-like算子以及文獻(xiàn)[20]中的插值算子.

如文獻(xiàn)[12,21]中所述, 誤差分析要做到與1/ν一致的關(guān)鍵在于Bertoluzza所證明的離散commutator 性質(zhì), 即

引理4.2對(duì)于任意的u∈W1,∞(Ω),vh∈Vh，以下不等式成立:

‖(I-Π)(u·vh)‖0≤Ch‖u‖1,∞‖vh‖0.

(9)

再次利用Cauchy-Schwarz不等式得

利用上兩式,再將(7)式乘上2Δt,并從0加到n得

對(duì)任意的1+Δt

利用

C0(N+1)≥Δt?

定理4.4(離散壓力的穩(wěn)定性) 對(duì)任意的ph∈Qh, 存在一個(gè)與h,Δt和ν無(wú)關(guān)的正實(shí)數(shù)β, 使得

β‖ph‖0

(10)

證明 由連續(xù)inf-sup條件可知, 對(duì)任意ph∈Qh?Q, 存在v∈X使得

?·v=ph,|v|1≤C‖ph‖0.

由Green公式,(4)式中給出的算子ju的正交性和近似性以及Poincaré不等式可得

(?·(v-juv),ph)+(?·juv,ph)=

-(v-juv,κh?ph)+(?·juv,ph)≤

Ch|v|1·‖κh?ph‖0+(?·juv,ph).

又由ju:X→Xh有

定理得證.

定理4.5(離散速度的誤差估計(jì)) 令(u,p)∈X×Q為問(wèn)題(1)的解.我們假設(shè)?tf∈L∞(L2)及(u,p)有額外的正則性,即(u,p)∈C0(Hs+1)×C0(Hs+1).又令存在常數(shù)C, 使得

‖?tu‖L∞(Hs+1)+‖?ttu‖L∞(L2)+

‖u‖L∞(W1,∞)≤C.

(11)

證明 記

E1+E2+E3+E4+E5+E6.

由Young不等式可得

由

Cauchy-Schwarz不等式, Young不等式及(4)式可得

關(guān)于對(duì)流項(xiàng)之差，我們將其分成如下形式：

E21+E22+E23+E24+E25,

這里我們利用了b(u;v,v)=0. 由Cauchy-Schwarz不等式,u的正則性,(4)式,逆不等式以及Young不等式有

類似地,有

及

對(duì)于第二項(xiàng), 我們采用文獻(xiàn)[9]中的技巧,有

令

在方程式(2)中取時(shí)間t=tn+1后減去(7)式, 取測(cè)試函數(shù)為

(v,q)=(vh,qh)=

后有

對(duì)于上式(4)右端第一項(xiàng), 我們有

將(4)式及逆不等式應(yīng)用于上式第一項(xiàng)可得

對(duì)于最后一項(xiàng), 我們?cè)僖淮卫貌坏仁?/p>

‖juun‖Wk,∞≤‖un‖Wk,∞,k=0,1

以及(4)式, 逆不等式和算子Π的穩(wěn)定性,得

由逆不等式和(4)式,我們得到了最后一項(xiàng)的估計(jì)如下：

由引理4.2有

由以上估計(jì), 我們有

對(duì)于第三項(xiàng)E3, 運(yùn)用Cauchy-Schwarz不等式, Young不等式以及(4)式有

然后, 由分部積分, 算子ju的正交性, Cauchy-Schwarz不等式及Young不等式有

對(duì)于E5, 我們將其分成兩部分

由Young不等式可得

于是

C+ChM‖un+1‖1,M,∞+‖un+1‖1,M,∞.

對(duì)于剩下的部分, 由三角不等式及(4)式

從而由u∈L∞(W1,∞(Ω))可得

C(h2+1)h2s+m1‖u‖L∞(Hs+1).

于是

C(h2+1)h2s+m1‖u‖L∞(Hs+1).

類似可得

最后, 我們得到E6的估計(jì)

結(jié)合以上估計(jì)式,分別取m1,m2=0,2,再由逆不等式和范數(shù)‖·‖LPS的定義可得

在上式不等號(hào)左右同時(shí)乘上Δt,并且從0加到n,再由引理4.2可推出

最終, 由三角不等式和算子ju,jp的近似性, 定理得證.

定理4.6(離散壓力的誤差估計(jì)) 令速度誤差估計(jì)中的條件成立,且時(shí)間和空間的劃分具有相同的階, 即存在與h和Δt無(wú)關(guān)的常數(shù)C, 使得Δt≤Ch.于是我們可得如下的壓力誤差估計(jì)：

C(1+(1+ν2+Eu(1+h-d))/β2)Eu+

C/β2(1+h2+Eu)Su

(12)

其中我們令離散速度的穩(wěn)定性和誤差估計(jì)的右端項(xiàng)分別為

和

(Δt)2.

證明 我們采用速度誤差方程中相同的誤差分解, 在(2)式中令時(shí)間t=tn+1后減去(7)式, 并令測(cè)試函數(shù)(v,q)=(vh,qh)=(vh,0)可得壓力的誤差方程

采用估計(jì)E1時(shí)的類似技巧和Poincaré不等式可得

不同于E2的估計(jì), 這里我們采用文獻(xiàn)[9]中給出的方法來(lái)估計(jì)對(duì)流項(xiàng)之差

|T2|=|b(un+1;un+1,vh)-

對(duì)上式第二項(xiàng),我們采用分部積分得

T21+T22.

由逆不等式知

于是

由Cauchy-Schwarz不等式,‖·‖LPS的定義和速度誤差估計(jì)的運(yùn)用可以得到

至于T5,由離散速度的穩(wěn)定性和離散Cauchy-Schwarz不等式可得

即

最后, 我們有

利用Δt≤Ch,結(jié)合以上的估計(jì)式可得

最終,由三角不等式,(4)式以及離散壓力的穩(wěn)定性，定理得證.