3維非定常Navier-Stokes方程組高效全離散有限元方法研究進展?

何銀年，馮新龍

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

0 引言

3維非定常不可壓縮Navier-Stokes方程組（簡稱N-S方程組）[1?14]:

其中:u=(u1(x,t),u2(x,t),u3(x,t))表示速度向量,p=p(x,t)表示壓力,f=f(x,t)表示給定外力,u0(x)表示初始速度,ν ＞0 表示粘性系數及T ＞0表示大時間區間長度,Ω ?R3表示有界區域.3維非定常不可壓縮N-S方程組的初邊值問題是描述不可壓縮流體運動最一般的規律,是典型的非線性問題,在工程應用和非線性科學研究方面具有廣泛的應用背景和重要的科學意義.然而,由于人們對N-S方程組的非線性現象本質和解的特性認識有限,使得數值求解N-S方程組成為一種十分重要的研究手段.但數值求解N-S方程組面臨三大困難:即如何處理不可壓縮約束條件div u=0,強非線性(u·?)u和長時間區間[0,T]積分的問題[1?3].

對于依賴于空間時間變量(x,t)的非定常N-S方程組的數值求解,人們首先是設計空間變量離散化,考慮半離散解(uh(x,t),ph(x,t))所滿足的近似N-S方程組,然后再考慮時間變量離散化,即全離散解所滿足的有限維近似N-S方程組;或者相反先設計時間變量離散化,考慮半離散解(un(x),pn(x))所滿足的近似N-S方程組,然后再考慮空間變量離散化,即考慮全離散解所滿足的有限維近似N-S方程組.由于篇幅有限,本文僅考慮前一種情形,且僅考慮有限元的空間離散情形,不考慮有限差分方法、有限體積元方法、譜方法等其它方法的空間離散化情形.在設計空間變量離散化時,人們首先考慮的是如何克服不可壓縮約束條件div u=0的困難.通常采取的方法是:(1)設計滿足inf-sup條件的協調或非協調的有限元空間對(Xh,Mh)[1?8,13];(2)對一般的有限元空間對(Xh,Mh)構造穩定化的弱變分形式[15?18];(3)設計滿足不可壓縮約束條件div uh=0的有限元空間對(Xh,Mh)[19?22].考慮文章的篇幅限制,本文僅考慮滿足inf-sup條件的一階協調有限元空間對(Xh,Mh)情形.

非定常N-S方程組在經過空間變量離散化后,得到關于時間變量的一個非線性常微分方程組.因此為了得到非定常N-S方程組的數值解,還需要對時間變量實施有限差分離散化.為了克服長時間區間[0,T]積分的困難,需要設計一個大時間步長的有限差分數值方法.高階時間精度的離散差分格式可以使得時間離散步長取大一些.考慮到非定常N-S方程組解的正則性限制,人們通常選取時間二階精度的差分離散化格式.時間二階精度的離散格式有全隱格式(例如Crank-Noclson差分格式[7]),半隱格式(例如Crank-Noclson外推差分格式[23?25]) 和隱式/顯式差分格式(例如Crank-Noclson/Adams-Bashforth差分格式[26?27])等.眾所周知,盡管二階精度的全隱格式是無條件穩定的,但是對于每個時間層n,需要求解關于的非線性方程組,人們承受不了這個巨大的計算量耗費.其次時間二階精度的半隱離散格式也具有線性化和無條件穩定的優點,克服了非線性的困難,但是對于每個時間層n,求解關于的線性方程組時,需要求解變系數矩陣的大型代數方程組,即其系數矩陣仍然依賴前兩步的速度向量人們仍然難以承受這個大的計算量耗費.最后,用時間二階精度的隱式/顯式差分格式求解時,對每個時間層n,僅僅需要求解常系數矩陣的大型代數方程組.這樣既克服了長時間積分和非線性的困難,又花費較少的計算量耗費.因此,把空間離散的滿足inf-sup條件的一階協調有限元空間對(Xh,Mh)方法和時間離散的二階精度的隱式/顯式差分格式結合起來,就得到求解非定常N-S方程組的高效有限元方法.對于時間二階精度的隱式/顯式差分格式,除了Crank-Nicolson/Adams-Bashforth(CN/AB)差分格式,還有二階后差隱式/顯式差分格式或稱Gear外推差分格式[28?29]等可以被應用于求解非定常N-S方程組.

眾所周知,盡管時間二階精度的CN/AB全離散有限元方法吸引了眾多學者的極大興趣,然而在以往數值分析方面,許多學者認為即使在初值足夠光滑的條件下,該數值方法的CFL穩定性和收斂性條件都要求時間離散步長τ嚴格地依賴于空間離散尺度h.在文獻[3]中,Marion及Temam提出了下列的穩定性和收斂性條件:

其中d=2,3表示空間區域Ω的維數.最近,Tone[27]給出了下列的收斂性條件:

另外,在修正的CN/AB情形下（非線性項和壓力項都是顯式的）,Johnston及Liu[30]提出了下列的穩定性條件:

最近,He及Sun[26]證明了求解二維非定常N-S方程組時穩定性和收斂性條件是τ ≤C,即是幾乎無條件穩定和收斂的.此外,He及Sun[26]也證明了下列的收斂性結論:對每一時間步tn∈(0,T],全離散解滿足下列的收斂率

其中:σ(t)=min{1,t},κ是一個依賴于參數(ν,Ω,T,u0,f)的一般性常數.這里數值速度在L2范數意義下關于時空變量(x,t)達到了最優收斂階,數值速度和壓力在H1-L2范數意義下關于時間步長τ沒有達到二階最優收斂.

在最新的研究中,He,Zhang及Zou證明了求解三維非定常N-S方程組時穩定性和收斂性條件是τ ≤C,并且也證明了下列的收斂性結論:對每一時間步tn∈(0,T],全離散解滿足下列的收斂率

本文第一節主要介紹3維非定常N-S方程組有關數學描述和已有的存在性、唯一性和正則性結論.第二節主要介紹3維非定常N-S方程組有限元空間離散的有關數學描述和已有的收斂性和穩定性結論.第三節主要介紹3維非定常N-S方程組的時間二階精度的CN/AB全離散有限元方法的有關數學描述和最新收斂性結論.第四節給出本文主要結論.

1 3維非定常N-S方程組的有關數學描述及預備知識

令Ω ?R3是具有Lipschitz連續邊界的有界單連通區域,Banach空間Lp(Ω)是區域Ω上的p次Lebesgue可積函數全體,且賦予范數

為了書寫簡便,我們常常利用記號u(t)表示函數u(x,t).對于p=2,L2(Ω)是Hilbert空間且賦予內積和范數:

進一步,我們也引進下列Hilbert空間H1(Ω),H2(Ω):

以及下列內積和范數.

由Green公式,可以導出(1)的變分形式:求u ∈L∞(0,T;L2(Ω)3)∩L2(0,T;X) 及p ∈L2(0,T;M),使得對任意(v,q)∈X×M 滿足

其中:X=(H10(Ω))3,M=L20(Ω)={q ∈L2(Ω);∫Ωq(x)dx=0},以及

如果(u,p)滿足(4),被稱為非定常N-S方程組的弱解.此外,如果弱解(u,p)還滿足正則性u ∈L∞(0,T;X)∩L2(0,T;(H2(Ω))3∩X)及p ∈L2(0,T;H1(Ω)∩M),則(u,p)被稱為非定常N-S方程組的強解.目前,對于一般的數據(ν,Ω,u0,f,T),非定常N-S方程組的弱解的唯一性和強解的存在性還是一個未解的重大難題[1?3,8?12,14].

本文中我們關于問題(1)中已知信息u0及f(t)作下列假設.

假設(A0):初始速度u0∈H2(Ω)3∩X且滿足?·u0=0,以及外力項f滿足

這里κ及以下κ0被用來表示依賴于數據(ν,Ω,T,u0,f)的正常數.

進一步,如果3維區域Ω是一凸多面體或邊界?Ω是C1,1的,則下列的正則性假設成立.

假設(A1):如果f ∈(L2(Ω))3,則Stokes方程組：

允許有唯一解(u,p)滿足下列正則性：

其中c0=c0(Ω)是依賴于Ω的常數.

為了以后研究非定常N-S方程組解的正則性及其數值解的收斂性,我們需要引進下列Gadliardo-Nirenberg不等式[1?7]:

其中c1和γ0為依賴于Ω的常數.利用不等式(6)～(7),我們可以導出非線性算子B(w,u)的下列關系式:

其中N0=N0(Ω)是依賴于Ω的常數及

應用(8)并在(4)中取檢驗函數(v,q)=(u,p),我們得到下列的關系式:

為了進一步得到解的正則性結論,我們需要對弱解(u,p)做出進一步的假設.

假設(A2):假定問題(4)的弱解(u,p)滿足

利用假設(A0)～(A2)及估計式(8)～(10),可以導出(u,p)如下的正則性結果[4].

定理1如果假設(A0)～(A2) 成立,則問題(4)的弱解(u,p)滿足:

2 3維非定常N-S方程組的有限元離散化

為了方便,本文僅考慮3維有界的凸多面體區域Ω,對Ω進行擬一致正則的四面體網格剖分四面體單元K的直徑表示為hK,并定義網格尺度為引進有限元空間Xh?X及Mh?M滿足下列的一階逼近和inf-sup條件假設[4].

假設(A3):假定有限元空間對Xh×Mh滿足

(1) 對每一v ∈H2(Ω)3∩X及q ∈H1(Ω)∩M,存在逼近函數πhv ∈Xh及ρhq ∈Mh使得

成立及inf-sup條件:對每一qh∈Mh,存在vh∈Xh,vh/=0 使得

成立,其中c2及β是依賴于Ω的正常數.

根據文獻[2,32-33],存在一系列有限元空間對Xh×Mh滿足假設(A3),比如P2-P0元:

由于有限元空間對Xh×Mh滿足假設(A3),可以定義有限元空間

及離散Stokes算子Ah=-PhΔh,其中Ph是由(L2(Ω))3到Vh的L2-投影算子,滿足

以及離散拉普拉斯算子-Δh被定義為

由此可以定義離散的Sobolev范數

根據離散Stokes算子Ah的定義,結合(7),可得下列離散Gadliardo-Nirenberg不等式[4]:

于是可以由(14)推出非線性項B(wh,uh)滿足下列的有界性:

現在我們定義非定常N-S方程組(4)的有限元變分形式：對每一t ∈(0,T] 求(uh(t),ph(t))∈Xh×Mh使得對每個(vh,qh)∈Xh×Mh滿足

其中uh(0)=Phu0.在(17)中,取(vh,qh)=(uh,ph)及利用(8),我們可得到

再利用假設(A0)～(A3)和B(w,u)性質及定理1關于(u,p)的正則性結論,我們可以導出(uh,ph)關于(u,p)的下列誤差估計結果.

定理2如果假設(A0)～(A3)成立,則問題(17)的有限元解(uh,ph)滿足下列誤差估計:

由于假設(A2)～(A3),(18)～(19),我們可以證明有限元解(uh,ph)滿足下列穩定性[7]:

最后根據(15)～(20),我們可以證明對于每一t ∈(0,T],有限元解(uh,ph)滿足下列進一步的穩定性:

上述穩定性結果是進行全離散有限元解的基礎保障.

3 3維非定常N-S方程組的NS/AB全離散有限元解及其最優誤差估計

成立.利用不等式關系(8),(15)～(16),(24),定理2,有限元解(uh,ph)的穩定性以及歸納法原理,可以得到滿足誤差估計結果.

4 結論

求解三維非定常N-S方程組的CN/AB全離散有限元方法在每個時間步都是解線性方程組(Stokes方程組),對每一時間層n,其系數矩陣是固定不變的.在該數值格式中線性項用隱式格式離散以增加其格式的穩定性能,非線性項用顯式格式離散以保證格式的簡單性.此外,現有研究分析表明該方法是幾乎無條件穩定和收斂的,并且關于時間步長τ是最優階收斂的,即該方法不要求時間離散步長τ不隨著空間網格尺度h的縮小而影響時間步長τ的選取,從而允許選取大的時間步長τ,因此該方法也克服了求解3維非定常N-S方程組的強非線性和長時間區間積分的困難,另外由于使用滿足inf-sup條件的協調有限元空間對Xh×Mh克服了不可壓縮約束條件的困難.綜上所述,CN/AB全離散有限元方法是求解三維非定常N-S方程組的高效有限元方法.