二階混合有限體積法求解Navier-Stokes 方程的穩定性及誤差估計

張杰華, 韓明華

(凱里學院理學院，貴州 凱里 556011)

1 引言

Navier-Stokes(N-S)方程是計算流體力學(CFD)中的基本方程，而局部守恒性是CFD 數值方法發展的一個重要準則.有限體積法(FVM)最大的特點就是能在其控制體積上保持某種物理量的局部守恒，同時繼承了有限元法和有限差分法的優點，因此FVM 已成為解決CFD 問題的一種流行數值方法.

近年來FVM 的研究非常活躍.文獻[1,2]分析和構造了一類橢圓問題的二階和三階FVM.它們根據三角網格上的某些幾何要求，建立了離散雙線性型一致局部橢圓性的充要條件，為研究各類高階FVM 方程的穩定性提供了一種通用的方法.文獻[3]利用兩套網格—速度場的粗網格和壓力場的細網格，在三角形網格上建立了求解N-S 方程的一階FVM，并證明了其穩定性.文獻[4]利用增加穩定項的方法使其N-S 方程對應的一階FVM 方程滿足穩定性，從而得到三角形網格上的最優誤差估計.另外還有許多其它關于研究FVM 求解N-S 方程穩定性理論的工作，如文獻[5–13]及其參考文獻.然而，這些工作大多是局限于在三角形網格上采用分片線性多項式函數逼近速度場的低階FVM，其主要思想是利用增加穩定項的技術使N-S 方程的一階FVM 方程滿足穩定性條件.這在一定程度上導致了求解FVM 方程的計算量增大及其穩定性理論分析的難度，并且其理論結果只有關于速度項的一階H1模誤差估計.到目前為止，很少出現利用二次多項式函數逼近N-S 方程速度場的高階FVM 的相關理論研究工作.

本文考慮利用二階混合FVM 求解定常N-S 方程，即對速度場的試驗函數空間取為分層二次多項式函數的有限元空間，檢驗函數空間由對偶剖分網格上的特征函數及二次多項式函數所組成，壓力場的試驗函數空間和檢驗函數空間均取為線性有限元空間.利用二階混合FVM 的離散雙線性型與其對應的有限元法離散雙線性型之間的等價關系，本文證明了求解N-S 方程二階混合FVM 的穩定性，進而得到了二階混合FVM 關于速度項和壓力項的最優誤差估計收斂階.

本文有如下幾個特點.首先，本文構建的混合FVM 的對偶剖分網格與線性FVM 的相同，從而其理論分析與其技術的實現都較一般傳統FVM 簡單.另外，本文不用增加穩定項，不需要粗細兩套三角形網格，即可得到求解N-S 方程的二階混合FVM 的穩定性及其最優誤差估計.其次，對于N-S 方程中的非線性項，本文直接在FVM 的控制體積上對其進行離散化，從而利用其FVM 的非線性項與其對應的有限元法的非線性項之間的關系，使得二階FVM 的理論分析及其計算技術的實現更加方便.此外，本文首次證明了二階FVM 求解N-S 方程關于速度場的H1模和壓力場的L2模的最優誤差估計.不僅本文的理論結果與其對應的有限元法結論相同，而且其證明方法具有一般性，可被應用于更一般的高階及高維有限體積法理論分析的情形.

第2 節將介紹一些相關記號和N-S 方程的變分問題.第3 節建立求解N-S 方程的二階混合FVM.第4 節將證明此二階FVM 的穩定性及其該方法的最優階誤差估計.第5 節將給出數值算例來驗證本文結論的正確性.

2 Navier-Stokes 方程

記? 為R2中的一個多邊形區域，?? 為其邊界，設

并用//·//0表示L2(?)空間的范數，用|·|m和//·//m分別表示Hm(?)空間的半范數和全范數.設C,c(或含下標)表示一類與網格尺寸無關的正常數，與文中出現的(μ,?,f)有關.

本文考慮如下定常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]

本文假設原方程(1)總存在唯一解，即μ與f 滿足如下條件[15]

3 有限體積法

本節將建立求解N-S 方程(1)的二階混合FVM 方程.下面依次介紹FVM 的三角形網格剖分及其試驗函數空間，對偶網格剖分及其檢驗函數空間，最后給出求解N-S 方程(1)的FVM 方程.

設T:={K}為? 上的擬一致三角形網格剖分[16,17]，即存在一個與網格尺寸h

無關的正常數θinf，使得三角形單元K ∈T的最小內角θmin,K滿足θmin,Kθinf，且存在常數γ> 0，對所有的K ∈T滿足γh2meas(K)h2.記T:={T}為? 上的擬一致三角形剖分族.設取N-S 方程(1)中的速度場與壓力場的試驗函數空間分別為如下的有限元函數空間

其中Pr(D)表示在區域D上的r次多項式函數空間.具體地，取速度場的試驗函數空間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}的內部節點的集合，則?P表示在節點P ∈N處的線性插值基函數，即?P ∈P1(K)，且若P=P′，則?P(P′)=1，若P ?=P′，則?P(P′)=0.壓力場的試驗函數空間及檢驗函數空間均取為QT，其基函數集為

ΨT:={?P:P ∈N}.

下面介紹二階混合FVM 的對偶剖分及其關于速度項的檢驗函數空間[1,2,18].分別用直線連結三角形的重心與各邊的中點，則將每個三角形K ∈T分解成三個區域.把包含三角形頂點的多邊形區域稱之為一個控制體積K?，則所有的這些控制體積K?就形成了區域? 的一個對偶剖分T ?，即T ?:={K?}.顯然，二階混合FVM 的對偶剖分與線性FVM 的相同，因此其理論分析較傳統Lagrange 二階有限體積法簡單.取速度場的檢驗函數空間為UT ?滿足

UT ?=V0,T ?⊕W2,T,

這里V0,T ?:= span{χP:P ∈N}為控制體積K?上的分片常數函數空間，χP為其頂點P ∈N處的控制體積的特征函數.顯然dimUT=dimUT ?.下面引入從UT到UT ?的一個線性映射.對任意的

設ΠT ?:UT →UT ?滿足

這里cP,cP,P′為向量常數.

利用從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)的第一式加上連續性方程后在每個控制體積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)的關系，可以簡化FVM(6)的理論分析.

方程(6)即為本文所分析的FVM 方程.注意到UT ?包含控制體積上的特征函數，因此FVM(6)關于速度項的數值解具有局部守恒性.

4 穩定性及誤差估計分析

本節主要分析FVM(6)在三角形網格上的穩定性及其誤差估計.注意到FVM(6)中含壓力項的離散雙線性型與其對應的有限元法離散雙線性型具有等價的關系，以及FVM(6)的橢圓雙線性型滿足強制性，從而在粘度μ和外力f 滿足一定的條件下(與有限元法類似)，利用Brouwers 不動點定理，可以得到FVM(6)的穩定性及其方程解的唯一性.一旦FVM(6)的穩定性成立，即可推出其關于速度場的H1(?)模與壓強場的L2(?)模的最優誤差估計.

首先，我們引入一些結論.關于線性映射ΠT ?: UT →UT ?，由文獻[18,19]可知如下性質成立.

引理1對任意的u∈UT，有

若T為擬一致剖分，則對任意的u∈UT，有

由文獻[14]可知，對任意的(v,q)∈(UT,QT)，雙線性型dT(v,q)具有如下估計.

引理2若T為擬一致剖分，則對任意的(v,q)∈(UT,QT), T ∈T，存在常數C，使得

|dT(v,q)|C|v|1//q//0.

另外，對任意的q ∈QT，存在常數c滿足

結合引理1、引理2，由文獻[19]可知：

引理3對任意的(v,q)∈(UT,QT)，有

結合不等式(10)，則對任意的q ∈QT，存在常數c滿足

引理4若T為擬一致剖分，則對所有的(v,p)∈(UT,L20(?)∩H1(?)), T ∈T，存在常數C滿足

引理5若T為擬一致剖分，θmin,K> 14.18?, K ∈T，則對所有的u∈(?)∩(?),v∈UT, T ∈T，存在常數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，存在常數C, c5滿足

有了上述準備，下面證明FVM(6)的穩定性及解的唯一性.記//·//m,p為Soblev 空間Wm,p(?)的范數，如下結論將會被用到.

定理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]證明結論.記

對任意的(v,q)∈(UT,QT)，設映射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(·,·,·)的定義及關系式(3)可知，對任意的w,u,v∈UT，有

下面估計等式(23)左邊的第三項.注意到u∈UT ?W1,∞(?)，利用Cauchy-Schwartz 不等式及估計(9)式和(17)式，則可得

注意到(w,·)∈BM，即

因此利用不等式(28)和條件(19)式，可將估計(27)式化為

對于等式(23)的右端項，利用Cauchy-Schwartz 不等式及范數估計(8),(16)可得

從而將上述關系式(24),(26),(29),(30)分別代入等式(23)，整理可得

再結合條件(18)式，因此有

下面考慮//p//0的估計.利用上述類似的方法，注意到w,u,v∈UT，則由估計(25),(27)及關系式(4)可得

取

則由估計(32), (33)及(28)可得

另外，由引理6 中的第二個不等式，結合等式(22)，可知對任意的(u,p)∈(UT,QT)有

將估計(30),(34)代入上式可知，對任意的p ∈QT，有

結合估計(31)可得

從而PT(BM)?BM.因此由Brouwers 不動點定理可知FVM(6)至少存在一組解.

下面證明FVM(6)解的唯一性.設(u,p)∈(UT,QT)和(ˉu,ˉp)∈(UT,QT)為FVM(6)的兩組解，即對任意的(v,q)∈(UT,QT)，有

易知u,均滿足估計(31).利用上面兩式作差可得

這里用到了

設(v,q) = (u?,p ?ˉp) := (e,η)∈(UT,QT)，則由估計(24)的結論可知，等式(35)左端第一項滿足

另外，由關系式(25)和(26)可知，等式(35)左端第二與第三項之和滿足

上式結合估計(27)和(4)得

再由估計(31)和條件(19)可知，等式(35)左端第二與第三項之和滿足

利用估計(36)和(38)，整理等式(35)，可得

再由條件(20)的第二式可得

因此p=.

綜上所述，FVM(6)在條件(20)下存在唯一解.

此定理說明，當網格尺寸h加密到一定的粗細時，即條件(18)成立時，FVM(6)存在數值解.若μ與f 滿足條件(20)的第二式，則當網格尺寸h進一步加密時，FVM(6)具有唯一解.結合假設條件(5)與條件(20)可知，在原問題(1)存在唯一解的條件下，要使FVM(6)也有唯一解，那么μ與f 須滿足

注意到由//f//?的定義及范數估計(16)可知

一般地可取4C3> 2c4(c4)2，從而可知當μ與f 滿足條件(20)時，條件(39)一定成立.即，在FVM(6)存在唯一解的條件下，原問題(1)也一定有解.

有了FVM(6)的穩定性，下面將考慮FVM(6)關于速度場的H1(?)模與壓強場的L2(?)模的誤差估計.為此，需要一個強于(20)的假設條件

其中f∈L2(?).即要求μ與f 滿足

這里仍取4C3> 2c4(c4)2.從而利用上述FVM(6)的穩定性結及插值誤差估計，可得如下定理.

定理2設(u,p)∈((?)∩H3(?),(?)∩H2(?))為原方程(1)的真解，(uT,pT)∈(UT,QT)為FVM(6)的數值解.若T為擬一致剖分，θmin,K> 14.18?, K ∈T，則當h> 0 滿足條件(18)，且f∈L2(?)和μ> 0 滿足(40)時，存在常數C使得對所有的T ∈T，有

為了證明定理2，下面首先介紹兩個引理.

引理7如果定理2 的所有條件成立，則存在常數C滿足這里ξh:=uT ?u∈UT, ηh:=pT ?p ∈QT，其中,為插值算子.

證明 此引理的證明主要利用離散雙線性型aT(·,ΠT ?·)的連續性和橢圓性，離散雙線性型bT(ΠT ?·,·)的連續性，另外還用到插值誤差估計.

在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(?))，并利用范數估計(16)可知

即

其中f∈L2(?).結合估計(53)和(54)，利用插值誤差估計可得

下面考慮等式(49)右端的第四項的估計.事實上，利用估計(21)的結論

及插值誤差估計，等式(49)右邊的第四項滿足如下估計

因此將估計(50),(51),(52),(55),(56)代入等式(49)右端，可得

接下來考慮等式(49)左端項的估計.類似于不等式(37)的推導，可得

這里

注意到c42C3，利用條件(18)和(19)，則由不等式(58)可得

因此，結合估計(15)的結論，可知等式(49)的左端滿足如下估計

其中

考慮到條件(20)和(40)，可知N(f)>0.因此由不等式(59)可得

再結合估計(57)和(60)，可得

即證結論(43).

引理8如果定理2 的所有條件成立，則存在常數C滿足

證明 此結論的證明利用雙線性型bT(ΠT ?·,·)的inf-sup 條件，雙線性型aT(·,ΠT ?·)的連續性，以及插值誤差估計.

首先，由等式(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，有

結合(69)和(70)易得

即得估計(61).

定理2的證明：利用引理7 和引理8 的結論，得到|u?uT|1和//p ?pT //0的估計.

首先由估計(43)和(61)易知

為了方便，記|||(u,p)|||:=h2(//u//3+//p//2)，則利用Young 不等式可得

此處ε>0 為一任意常數.取ε=C，則利用不等式(72)，估計(71)可化為

即

因此，利用三角不等式，插值誤差估計和估計(73)可得

另一方面，結合估計(73)和(61)，利用三角不等式及插值誤差估計可得

因此由估計(74)和(75)即得結論(42).

此定理說明，當μ和f 滿足條件(40)或(41)時，FVM(6)存在唯一解且具有最優誤差估計.進一步分析可以發現條件(41)意味著C1C3C6>//f//0，即要求//f//0不宜過大.另外，在實際問題中粘度往往非常小，從而條件(41)的第一個不等式可恒成立，也就是條件(40)的第二個不等式可恒成立.因此只需考慮條件(40)的第一個不等式，即條件(20)的第二式.綜合前述的穩定性，可知當μ和f 滿足條件(20)時，原問題(1)有解，且FVM(6)存在唯一解及其最優誤差估計.

此定理的證明具有一般性.這種方法可以直接應用到更一般的高階及高維情形的誤差估計證明.

5 算例

本節將利用數值算例來檢驗FVM(6)的有效性.

步驟1取迭代初值:=0(n=0)代入方程(76)，求解對應的線性方程(即Stokes方程)可得.顯然滿足條件(21).

為非線性方程FVM(6)的解.

下面檢驗FVM(6)在三角形網格上的誤差估計.考慮將正方形區域? := (0,1)×(0,1)進行一致三角形網格剖分，網格大小為h=1/N.具體地，先將? 分解成N×N個矩形網格，然后用對角線將每個正方形分成兩個三角形.兩個算例如下.

例1 設N-S 方程(1)的真解為u=(u1,u2)和p，其中

例2設N-S 方程(1)的真解為u=(u1,u2)和p，其中

例1 和例2 中的解均滿足條件u|??=0,∫?p=0.

為使N-S 方程(1)中的μ和f 滿足條件(20)的第二式，此處取μ=1.表1 和表2 分別列出了例1 和例2 的數值計算結果，其中|u?uT|1, //p ?pT //0分別為速度場的H1模和壓強場的L2模誤差，r(u), r(p)分別表示速度場與壓強場所對應的誤差收斂階.數值結果表明本文的計算方法正確，且當網格逐漸加密(滿足條件(20)的第一式)時，FVM(6)具有速度項的最優H1模和壓強項的最優L2模的誤差收斂階O(h2)，顯然這與本文的理論結果相一致.從而也驗證了求解N-S 方程(1)的二階混合FVM(6)與相對應的有限元法具有相同的誤差收斂階.

表1 例1 的誤差估計和收斂階

表2 例2 的誤差估計和收斂階

下面考慮在正方域? = (0,1)×(0,1)上利用FVM(6)對經典lid-driven cavity 問題進行數值模擬.

例3 設N-S 方程(1)中的μ=1,f =0.真解u=(u1,u2)滿足如下邊界條件

顯然，N-S 方程(1)中的μ和f 滿足條件(20).利用GMSH 軟件自動生成區域? 的非結構三角形網格剖分，見圖1.根據FVM(6)的數值計算結果，利用TECPLOT 軟件畫出了其速度場的流線圖，見圖2.從流線圖中可以看出，對于具有粘性(系數μ= 1)的方腔驅動流體，當邊界以一定的速度移動時，方腔內會產生一個漩渦流，且漩渦中心位于方腔的中間位置.流線的疏密可反映速度的大小，因此可見靠近方腔驅動壁面上的速度要大于其它地方的速度.總之，模似結果顯示出了良好的計算效果.

圖1 三角形網格剖分

圖2 速度場流線圖

6 結論

由定理1 和定理2 可知，當μ和f 滿足適當的條件時(類似于但略強于有限元法)，例如條件(20)，則求解N-S 方程(1)的FVM(6)存在唯一解及其最優誤差估計.在有限體積法中，μ和f 所需滿足條件略強于對應的有限元法，這主要是由有限體積法的橢圓離散雙線性型的不對稱性所引起的，因此這是不可避免的.不過有限體積法可以保持物理量的局部守恒性(這是傳統有限元法不具備的)，且計算量比有限元法(尤其是間斷有限元法)小得多，這正是利用有限體積法處理CFD 問題的優勢所在.