張杰華, 韓明華
(凱里學(xué)院理學(xué)院,貴州 凱里 556011)
Navier-Stokes(N-S)方程是計算流體力學(xué)(CFD)中的基本方程,而局部守恒性是CFD 數(shù)值方法發(fā)展的一個重要準(zhǔn)則.有限體積法(FVM)最大的特點就是能在其控制體積上保持某種物理量的局部守恒,同時繼承了有限元法和有限差分法的優(yōu)點,因此FVM 已成為解決CFD 問題的一種流行數(shù)值方法.
近年來FVM 的研究非?;钴S.文獻[1,2]分析和構(gòu)造了一類橢圓問題的二階和三階FVM.它們根據(jù)三角網(wǎng)格上的某些幾何要求,建立了離散雙線性型一致局部橢圓性的充要條件,為研究各類高階FVM 方程的穩(wěn)定性提供了一種通用的方法.文獻[3]利用兩套網(wǎng)格—速度場的粗網(wǎng)格和壓力場的細(xì)網(wǎng)格,在三角形網(wǎng)格上建立了求解N-S 方程的一階FVM,并證明了其穩(wěn)定性.文獻[4]利用增加穩(wěn)定項的方法使其N-S 方程對應(yīng)的一階FVM 方程滿足穩(wěn)定性,從而得到三角形網(wǎng)格上的最優(yōu)誤差估計.另外還有許多其它關(guān)于研究FVM 求解N-S 方程穩(wěn)定性理論的工作,如文獻[5–13]及其參考文獻.然而,這些工作大多是局限于在三角形網(wǎng)格上采用分片線性多項式函數(shù)逼近速度場的低階FVM,其主要思想是利用增加穩(wěn)定項的技術(shù)使N-S 方程的一階FVM 方程滿足穩(wěn)定性條件.這在一定程度上導(dǎo)致了求解FVM 方程的計算量增大及其穩(wěn)定性理論分析的難度,并且其理論結(jié)果只有關(guān)于速度項的一階H1模誤差估計.到目前為止,很少出現(xiàn)利用二次多項式函數(shù)逼近N-S 方程速度場的高階FVM 的相關(guān)理論研究工作.
本文考慮利用二階混合FVM 求解定常N-S 方程,即對速度場的試驗函數(shù)空間取為分層二次多項式函數(shù)的有限元空間,檢驗函數(shù)空間由對偶剖分網(wǎng)格上的特征函數(shù)及二次多項式函數(shù)所組成,壓力場的試驗函數(shù)空間和檢驗函數(shù)空間均取為線性有限元空間.利用二階混合FVM 的離散雙線性型與其對應(yīng)的有限元法離散雙線性型之間的等價關(guān)系,本文證明了求解N-S 方程二階混合FVM 的穩(wěn)定性,進而得到了二階混合FVM 關(guān)于速度項和壓力項的最優(yōu)誤差估計收斂階.
本文有如下幾個特點.首先,本文構(gòu)建的混合FVM 的對偶剖分網(wǎng)格與線性FVM 的相同,從而其理論分析與其技術(shù)的實現(xiàn)都較一般傳統(tǒng)FVM 簡單.另外,本文不用增加穩(wěn)定項,不需要粗細(xì)兩套三角形網(wǎng)格,即可得到求解N-S 方程的二階混合FVM 的穩(wěn)定性及其最優(yōu)誤差估計.其次,對于N-S 方程中的非線性項,本文直接在FVM 的控制體積上對其進行離散化,從而利用其FVM 的非線性項與其對應(yīng)的有限元法的非線性項之間的關(guān)系,使得二階FVM 的理論分析及其計算技術(shù)的實現(xiàn)更加方便.此外,本文首次證明了二階FVM 求解N-S 方程關(guān)于速度場的H1模和壓力場的L2模的最優(yōu)誤差估計.不僅本文的理論結(jié)果與其對應(yīng)的有限元法結(jié)論相同,而且其證明方法具有一般性,可被應(yīng)用于更一般的高階及高維有限體積法理論分析的情形.
第2 節(jié)將介紹一些相關(guān)記號和N-S 方程的變分問題.第3 節(jié)建立求解N-S 方程的二階混合FVM.第4 節(jié)將證明此二階FVM 的穩(wěn)定性及其該方法的最優(yōu)階誤差估計.第5 節(jié)將給出數(shù)值算例來驗證本文結(jié)論的正確性.
記? 為R2中的一個多邊形區(qū)域,?? 為其邊界,設(shè)
并用//·//0表示L2(?)空間的范數(shù),用|·|m和//·//m分別表示Hm(?)空間的半范數(shù)和全范數(shù).設(shè)C,c(或含下標(biāo))表示一類與網(wǎng)格尺寸無關(guān)的正常數(shù),與文中出現(xiàn)的(μ,?,f)有關(guān).
本文考慮如下定常N-S 方程:求u := (u1,u2)∈(?)∩H2(?)和p ∈(?)∩H1(?)滿足
此處,f := (f1,f2)為作用于流體上的外力,μ> 0 為流體粘度,u 為所求速度場,p為所求壓力場.
方程(1)的變分形式為:求u∈(?)∩H2(?), p ∈(?)∩H1(?)滿足
此處
利用Green 公式和方程(1)的第二式及其邊界條件可知,對任意的u,v,w∈(?)∩H2(?),三線性型c(·,·,·)滿足
另外,對任意的u,v,w∈(?)∩H2(?),三線性型c(·,·,·)還滿足[14,15]
本文假設(shè)原方程(1)總存在唯一解,即μ與f 滿足如下條件[15]
本節(jié)將建立求解N-S 方程(1)的二階混合FVM 方程.下面依次介紹FVM 的三角形網(wǎng)格剖分及其試驗函數(shù)空間,對偶網(wǎng)格剖分及其檢驗函數(shù)空間,最后給出求解N-S 方程(1)的FVM 方程.
設(shè)T:={K}為? 上的擬一致三角形網(wǎng)格剖分[16,17],即存在一個與網(wǎng)格尺寸h
無關(guān)的正常數(shù)θinf,使得三角形單元K ∈T的最小內(nèi)角θmin,K滿足θmin,Kθinf,且存在常數(shù)γ> 0,對所有的K ∈T滿足γh2meas(K)h2.記T:={T}為? 上的擬一致三角形剖分族.設(shè)取N-S 方程(1)中的速度場與壓力場的試驗函數(shù)空間分別為如下的有限元函數(shù)空間
其中Pr(D)表示在區(qū)域D上的r次多項式函數(shù)空間.具體地,取速度場的試驗函數(shù)空間UT為[1,2,18]
UT=V1,T ⊕W2,T,
此處
V1,T:=span{?P:P ∈N},W2,T:=span{?P?P′:P,P′∈N}.
記N表示剖分T={K}的內(nèi)部節(jié)點的集合,則?P表示在節(jié)點P ∈N處的線性插值基函數(shù),即?P ∈P1(K),且若P=P′,則?P(P′)=1,若P ?=P′,則?P(P′)=0.壓力場的試驗函數(shù)空間及檢驗函數(shù)空間均取為QT,其基函數(shù)集為
ΨT:={?P:P ∈N}.
下面介紹二階混合FVM 的對偶剖分及其關(guān)于速度項的檢驗函數(shù)空間[1,2,18].分別用直線連結(jié)三角形的重心與各邊的中點,則將每個三角形K ∈T分解成三個區(qū)域.把包含三角形頂點的多邊形區(qū)域稱之為一個控制體積K?,則所有的這些控制體積K?就形成了區(qū)域? 的一個對偶剖分T ?,即T ?:={K?}.顯然,二階混合FVM 的對偶剖分與線性FVM 的相同,因此其理論分析較傳統(tǒng)Lagrange 二階有限體積法簡單.取速度場的檢驗函數(shù)空間為UT ?滿足
UT ?=V0,T ?⊕W2,T,
這里V0,T ?:= span{χP:P ∈N}為控制體積K?上的分片常數(shù)函數(shù)空間,χP為其頂點P ∈N處的控制體積的特征函數(shù).顯然dimUT=dimUT ?.下面引入從UT到UT ?的一個線性映射.對任意的
設(shè)ΠT ?:UT →UT ?滿足
這里cP,cP,P′為向量常數(shù).
利用從UT到UT ?的線性映射ΠT ?,建立如下求解N-S 方程(1)的FVM 方程:求(u,p)∈(UT,QT),使得對所有的(v,q)∈(UT,QT)滿足
這里
其中n 為?K?上的單位外法向量.此處的非線性項cT(u,w,ΠT ?v),u,w,v∈UT,是利用方程(1)的第一式加上連續(xù)性方程后在每個控制體積K?∈T ?上進行積分求和而得到.特別地,若用v∈UT代替非線性項cT(u,w,ΠT ?v)中的ΠT ?v∈UT ?,則有
cT(u,w,v)=c(u,w,v), ?u,w,v∈UT.
因此利用非線性項cT(u,w,ΠT ?v)與c(u,w,v)的關(guān)系,可以簡化FVM(6)的理論分析.
方程(6)即為本文所分析的FVM 方程.注意到UT ?包含控制體積上的特征函數(shù),因此FVM(6)關(guān)于速度項的數(shù)值解具有局部守恒性.
本節(jié)主要分析FVM(6)在三角形網(wǎng)格上的穩(wěn)定性及其誤差估計.注意到FVM(6)中含壓力項的離散雙線性型與其對應(yīng)的有限元法離散雙線性型具有等價的關(guān)系,以及FVM(6)的橢圓雙線性型滿足強制性,從而在粘度μ和外力f 滿足一定的條件下(與有限元法類似),利用Brouwers 不動點定理,可以得到FVM(6)的穩(wěn)定性及其方程解的唯一性.一旦FVM(6)的穩(wěn)定性成立,即可推出其關(guān)于速度場的H1(?)模與壓強場的L2(?)模的最優(yōu)誤差估計.
首先,我們引入一些結(jié)論.關(guān)于線性映射ΠT ?: UT →UT ?,由文獻[18,19]可知如下性質(zhì)成立.
引理1對任意的u∈UT,有
若T為擬一致剖分,則對任意的u∈UT,有
由文獻[14]可知,對任意的(v,q)∈(UT,QT),雙線性型dT(v,q)具有如下估計.
引理2若T為擬一致剖分,則對任意的(v,q)∈(UT,QT), T ∈T,存在常數(shù)C,使得
|dT(v,q)|C|v|1//q//0.
另外,對任意的q ∈QT,存在常數(shù)c滿足
結(jié)合引理1、引理2,由文獻[19]可知:
引理3對任意的(v,q)∈(UT,QT),有
結(jié)合不等式(10),則對任意的q ∈QT,存在常數(shù)c滿足
引理4若T為擬一致剖分,則對所有的(v,p)∈(UT,L20(?)∩H1(?)), T ∈T,存在常數(shù)C滿足
引理5若T為擬一致剖分,θmin,K> 14.18?, K ∈T,則對所有的u∈(?)∩(?),v∈UT, T ∈T,存在常數(shù)C, c4滿足
這里
記C(·,·)為空間(UT,QT)×(UT,QT)上的雙線性型泛函,即令
C((u,p),(v,q)):=aT(u,ΠT ?v)+bT(ΠT ?v,p)+dT(u,q).
類似于文獻[19]中Theorem 7 的證明,利用引理2、引理4、引理5 及Babuska-Lax-Milgram 定理[20],可得如下估計.
引理6若T為擬一致剖分,θmin,K>14.18?, K ∈T,則對所有(u,p)∈(UT,QT),T ∈T,存在常數(shù)C, c5滿足
有了上述準(zhǔn)備,下面證明FVM(6)的穩(wěn)定性及解的唯一性.記//·//m,p為Soblev 空間Wm,p(?)的范數(shù),如下結(jié)論將會被用到.
定理1若T為擬一致剖分,θmin,K>14.18?, K ∈T,且對所有的h>0 滿足
此處
則FVM(6)至少存在一組解(u,p)∈(UT,QT).如果μ>0,f∈L2(?), h>0,且滿足
則FVM(6)存在唯一解,且滿足
證明 這里主要利用Brouwers 不動點定理[21]證明結(jié)論.記
對任意的(v,q)∈(UT,QT),設(shè)映射PT:(UT,QT)→(UT,QT)滿足
此處PT(w,p):=(PTw,p)=(u,p)∈(UT,QT)以及(w,·)∈BM.
先證FVM(6)解的存在性,即證PT(BM)?BM.首先考慮|u|1的估計.在等式(22)中取(v,q)=(u,p)∈(UT,QT),則有
利用引理3 中的等式(11),則上式可簡化為
下面考慮等式(23)中的每一項.對于第一項,由估計(15)可知,對任意的u∈UT,有
對于第二項,由cT(·,·,·)的定義及關(guān)系式(3)可知,對任意的w,u,v∈UT,有
下面估計等式(23)左邊的第三項.注意到u∈UT ?W1,∞(?),利用Cauchy-Schwartz 不等式及估計(9)式和(17)式,則可得
注意到(w,·)∈BM,即
因此利用不等式(28)和條件(19)式,可將估計(27)式化為
對于等式(23)的右端項,利用Cauchy-Schwartz 不等式及范數(shù)估計(8),(16)可得
從而將上述關(guān)系式(24),(26),(29),(30)分別代入等式(23),整理可得
再結(jié)合條件(18)式,因此有
下面考慮//p//0的估計.利用上述類似的方法,注意到w,u,v∈UT,則由估計(25),(27)及關(guān)系式(4)可得
取
則由估計(32), (33)及(28)可得
另外,由引理6 中的第二個不等式,結(jié)合等式(22),可知對任意的(u,p)∈(UT,QT)有
將估計(30),(34)代入上式可知,對任意的p ∈QT,有
結(jié)合估計(31)可得
從而PT(BM)?BM.因此由Brouwers 不動點定理可知FVM(6)至少存在一組解.
下面證明FVM(6)解的唯一性.設(shè)(u,p)∈(UT,QT)和(ˉu,ˉp)∈(UT,QT)為FVM(6)的兩組解,即對任意的(v,q)∈(UT,QT),有
易知u,均滿足估計(31).利用上面兩式作差可得
這里用到了
設(shè)(v,q) = (u?,p ?ˉp) := (e,η)∈(UT,QT),則由估計(24)的結(jié)論可知,等式(35)左端第一項滿足
另外,由關(guān)系式(25)和(26)可知,等式(35)左端第二與第三項之和滿足
上式結(jié)合估計(27)和(4)得
再由估計(31)和條件(19)可知,等式(35)左端第二與第三項之和滿足
利用估計(36)和(38),整理等式(35),可得
再由條件(20)的第二式可得
因此p=.
綜上所述,F(xiàn)VM(6)在條件(20)下存在唯一解.
此定理說明,當(dāng)網(wǎng)格尺寸h加密到一定的粗細(xì)時,即條件(18)成立時,F(xiàn)VM(6)存在數(shù)值解.若μ與f 滿足條件(20)的第二式,則當(dāng)網(wǎng)格尺寸h進一步加密時,F(xiàn)VM(6)具有唯一解.結(jié)合假設(shè)條件(5)與條件(20)可知,在原問題(1)存在唯一解的條件下,要使FVM(6)也有唯一解,那么μ與f 須滿足
注意到由//f//?的定義及范數(shù)估計(16)可知
一般地可取4C3> 2c4(c4)2,從而可知當(dāng)μ與f 滿足條件(20)時,條件(39)一定成立.即,在FVM(6)存在唯一解的條件下,原問題(1)也一定有解.
有了FVM(6)的穩(wěn)定性,下面將考慮FVM(6)關(guān)于速度場的H1(?)模與壓強場的L2(?)模的誤差估計.為此,需要一個強于(20)的假設(shè)條件
其中f∈L2(?).即要求μ與f 滿足
這里仍取4C3> 2c4(c4)2.從而利用上述FVM(6)的穩(wěn)定性結(jié)及插值誤差估計,可得如下定理.
定理2設(shè)(u,p)∈((?)∩H3(?),(?)∩H2(?))為原方程(1)的真解,(uT,pT)∈(UT,QT)為FVM(6)的數(shù)值解.若T為擬一致剖分,θmin,K> 14.18?, K ∈T,則當(dāng)h> 0 滿足條件(18),且f∈L2(?)和μ> 0 滿足(40)時,存在常數(shù)C使得對所有的T ∈T,有
為了證明定理2,下面首先介紹兩個引理.
引理7如果定理2 的所有條件成立,則存在常數(shù)C滿足這里ξh:=uT ?u∈UT, ηh:=pT ?p ∈QT,其中,為插值算子.
證明 此引理的證明主要利用離散雙線性型aT(·,ΠT ?·)的連續(xù)性和橢圓性,離散雙線性型bT(ΠT ?·,·)的連續(xù)性,另外還用到插值誤差估計.
在FVM(6)中令(u,p):=(uT,pT)∈(UT,QT),則有
利用(44)式減去(45)式,可得
及
這里ξ:=u?u∈(?)∩H2(?), η:=p ?p ∈(?)∩H1(?).
在等式(46)中取v=ξh ∈UT,在等式(47)中取q=ηh ∈QT,則由等式(11)可知
因此由等式(46)和(47),可得
注意到上式的最后兩項可展開為
再將其回代到等式(48)中可得
下面考慮等式(49)右端的每一項估計.首先,注意到ξ ∈∩H2?∩,因此由估計(14)及插值誤差估計可知,對任意的ξ ∈(?)∩H2(?), ξh ∈UT,有
再由估計(13)及插值誤差估計可知,對任意的ξh ∈UT, η ∈(?)∩H1(?),有
利用Cauchy-Schwartz 不等式及插值誤差估計,可知對任意的ξ ∈(?)∩H2(?), ηh ∈QT,有
下面考慮等式(49)右端的第三項估計.利用估計(32)和(33)可知,等式(49)右端的第三項滿足
因為u∈(?)∩H3(?)為方程(1)的解,因此必滿足變分方程(2).在方程(2)中取(v,q):=(u,p)∈((?)∩H2(?),(?)∩H1(?)),并利用范數(shù)估計(16)可知
即
其中f∈L2(?).結(jié)合估計(53)和(54),利用插值誤差估計可得
下面考慮等式(49)右端的第四項的估計.事實上,利用估計(21)的結(jié)論
及插值誤差估計,等式(49)右邊的第四項滿足如下估計
因此將估計(50),(51),(52),(55),(56)代入等式(49)右端,可得
接下來考慮等式(49)左端項的估計.類似于不等式(37)的推導(dǎo),可得
這里
注意到c42C3,利用條件(18)和(19),則由不等式(58)可得
因此,結(jié)合估計(15)的結(jié)論,可知等式(49)的左端滿足如下估計
其中
考慮到條件(20)和(40),可知N(f)>0.因此由不等式(59)可得
再結(jié)合估計(57)和(60),可得
即證結(jié)論(43).
引理8如果定理2 的所有條件成立,則存在常數(shù)C滿足
證明 此結(jié)論的證明利用雙線性型bT(ΠT ?·,·)的inf-sup 條件,雙線性型aT(·,ΠT ?·)的連續(xù)性,以及插值誤差估計.
首先,由等式(46)可知
其中
下面考慮等式(63)右端各項的估計.注意到
則由估計(32)和(33)知
將以上估計(65)代入等式(63)的右端可得
由估計(13)可知,對任意的(v,η)∈(UT,(?)∩H1(?)),有
由估計(14)可知,對任意的(ξ ?ξh)∈∩H2?∩,有
此處用到了有限元逆估計|ξh|2Ch?1|ξh|1, ξh ∈UT.
將以上估計(66),(67),(68)代入等式(62)的右端,并利用插值誤差估計可得
又由引理3 中的估計(12)可知,對任意的ηh=pT ?p ∈QT,有
結(jié)合(69)和(70)易得
即得估計(61).
定理2的證明:利用引理7 和引理8 的結(jié)論,得到|u?uT|1和//p ?pT //0的估計.
首先由估計(43)和(61)易知
為了方便,記|||(u,p)|||:=h2(//u//3+//p//2),則利用Young 不等式可得
此處ε>0 為一任意常數(shù).取ε=C,則利用不等式(72),估計(71)可化為
即
因此,利用三角不等式,插值誤差估計和估計(73)可得
另一方面,結(jié)合估計(73)和(61),利用三角不等式及插值誤差估計可得
因此由估計(74)和(75)即得結(jié)論(42).
此定理說明,當(dāng)μ和f 滿足條件(40)或(41)時,F(xiàn)VM(6)存在唯一解且具有最優(yōu)誤差估計.進一步分析可以發(fā)現(xiàn)條件(41)意味著C1C3C6>//f//0,即要求//f//0不宜過大.另外,在實際問題中粘度往往非常小,從而條件(41)的第一個不等式可恒成立,也就是條件(40)的第二個不等式可恒成立.因此只需考慮條件(40)的第一個不等式,即條件(20)的第二式.綜合前述的穩(wěn)定性,可知當(dāng)μ和f 滿足條件(20)時,原問題(1)有解,且FVM(6)存在唯一解及其最優(yōu)誤差估計.
此定理的證明具有一般性.這種方法可以直接應(yīng)用到更一般的高階及高維情形的誤差估計證明.
本節(jié)將利用數(shù)值算例來檢驗FVM(6)的有效性.
步驟1取迭代初值:=0(n=0)代入方程(76),求解對應(yīng)的線性方程(即Stokes方程)可得.顯然滿足條件(21).
為非線性方程FVM(6)的解.
下面檢驗FVM(6)在三角形網(wǎng)格上的誤差估計.考慮將正方形區(qū)域? := (0,1)×(0,1)進行一致三角形網(wǎng)格剖分,網(wǎng)格大小為h=1/N.具體地,先將? 分解成N×N個矩形網(wǎng)格,然后用對角線將每個正方形分成兩個三角形.兩個算例如下.
例1 設(shè)N-S 方程(1)的真解為u=(u1,u2)和p,其中
例2設(shè)N-S 方程(1)的真解為u=(u1,u2)和p,其中
例1 和例2 中的解均滿足條件u|??=0,∫?p=0.
為使N-S 方程(1)中的μ和f 滿足條件(20)的第二式,此處取μ=1.表1 和表2 分別列出了例1 和例2 的數(shù)值計算結(jié)果,其中|u?uT|1, //p ?pT //0分別為速度場的H1模和壓強場的L2模誤差,r(u), r(p)分別表示速度場與壓強場所對應(yīng)的誤差收斂階.數(shù)值結(jié)果表明本文的計算方法正確,且當(dāng)網(wǎng)格逐漸加密(滿足條件(20)的第一式)時,F(xiàn)VM(6)具有速度項的最優(yōu)H1模和壓強項的最優(yōu)L2模的誤差收斂階O(h2),顯然這與本文的理論結(jié)果相一致.從而也驗證了求解N-S 方程(1)的二階混合FVM(6)與相對應(yīng)的有限元法具有相同的誤差收斂階.
表1 例1 的誤差估計和收斂階
表2 例2 的誤差估計和收斂階
下面考慮在正方域? = (0,1)×(0,1)上利用FVM(6)對經(jīng)典lid-driven cavity 問題進行數(shù)值模擬.
例3 設(shè)N-S 方程(1)中的μ=1,f =0.真解u=(u1,u2)滿足如下邊界條件
顯然,N-S 方程(1)中的μ和f 滿足條件(20).利用GMSH 軟件自動生成區(qū)域? 的非結(jié)構(gòu)三角形網(wǎng)格剖分,見圖1.根據(jù)FVM(6)的數(shù)值計算結(jié)果,利用TECPLOT 軟件畫出了其速度場的流線圖,見圖2.從流線圖中可以看出,對于具有粘性(系數(shù)μ= 1)的方腔驅(qū)動流體,當(dāng)邊界以一定的速度移動時,方腔內(nèi)會產(chǎn)生一個漩渦流,且漩渦中心位于方腔的中間位置.流線的疏密可反映速度的大小,因此可見靠近方腔驅(qū)動壁面上的速度要大于其它地方的速度.總之,模似結(jié)果顯示出了良好的計算效果.
圖1 三角形網(wǎng)格剖分
圖2 速度場流線圖
由定理1 和定理2 可知,當(dāng)μ和f 滿足適當(dāng)?shù)臈l件時(類似于但略強于有限元法),例如條件(20),則求解N-S 方程(1)的FVM(6)存在唯一解及其最優(yōu)誤差估計.在有限體積法中,μ和f 所需滿足條件略強于對應(yīng)的有限元法,這主要是由有限體積法的橢圓離散雙線性型的不對稱性所引起的,因此這是不可避免的.不過有限體積法可以保持物理量的局部守恒性(這是傳統(tǒng)有限元法不具備的),且計算量比有限元法(尤其是間斷有限元法)小得多,這正是利用有限體積法處理CFD 問題的優(yōu)勢所在.