王小娟, 賈宏恩

(太原理工大學(xué) 數(shù)學(xué)學(xué)院, 山西 晉中030600 )

0 引言

考慮不可壓縮的Brinkman-Forchheimer方程:

其中：?是Rn(n=2或3)中有界連通并具有Lipschitz連續(xù)的區(qū)域, 其邊界為??,u是流體速度,p是壓力,ν>0是粘度系數(shù),f是體力, Forchheimer項(xiàng)α|u|r?2u(2 ≤r<∞,α>0)是由Forchheimers定律導(dǎo)出的, 用來(lái)描述各種物理情況, 如多孔介質(zhì)流動(dòng)、摩擦效應(yīng)以及某些耗散機(jī)制等.

方程(1)具有非線性和單調(diào)性. 近年來(lái), 一些學(xué)者致力于Brinkman-Forchheimer方程數(shù)值分析的研究. 如文獻(xiàn)[1]研究了Brinkman-Forchheimer方程的收斂性和連續(xù)依賴性; 文獻(xiàn)[2]給出了求解這一問(wèn)題的混合有限元方法,證明了弱解的存在唯一性并進(jìn)行了誤差估計(jì); 文獻(xiàn)[3]討論了混合有限元方法的超收斂性, 也可參考文獻(xiàn)[4-7].

文獻(xiàn)[8]提出了求解非線性問(wèn)題的兩重網(wǎng)格法, 該方法先在粗網(wǎng)格上求解一個(gè)非線性方程作為細(xì)網(wǎng)格的迭代初值近似,然后在細(xì)網(wǎng)格求解線性問(wèn)題,與傳統(tǒng)的方法相比可以節(jié)省大量的計(jì)算時(shí)間,已被廣泛應(yīng)用于求解其他非線性方程[9?14]. 文獻(xiàn)[15]提出求解Navier-Stokes方程的回溯兩重網(wǎng)格法.

本文將采用回溯兩重網(wǎng)格法求解問(wèn)題(1). 該方法需要先在粗網(wǎng)上求解非線性方程,再在細(xì)網(wǎng)格上求解Stokes問(wèn)題, 然后在粗網(wǎng)格上線性校正. 最后, 通過(guò)數(shù)值實(shí)驗(yàn)驗(yàn)證算法的有效性及穩(wěn)定性.

1 符號(hào)和準(zhǔn)備工作

Wk,p表示Sobolev空間,其中范數(shù)和半范為‖·‖k,p和|·|k,p.L2(?)上的內(nèi)積空間定義為:

引進(jìn)如下的雙線性、擬線性和三線性形式:

并且有以下不等式成立[4?7]:

引理1[15]假設(shè)?的邊界滿足強(qiáng)Lipschitz條件, 下面三線性形式的不等式成立:

其中：當(dāng)n=2時(shí),s=?>0(任意小); 當(dāng)n=3時(shí),

方程(1)的弱形式: 求(u,p)∈(X,M) 使得

問(wèn)題(3)的解的存在性和唯一性, 詳見(jiàn)文獻(xiàn)[7].

2 兩重網(wǎng)格法

(A1) 對(duì)于?(v,q)∈(Hk+1(?))2×Hk(?),k≥1, 存在(Iηv,Jηq)∈Xη×Mη使得

‖v?Iηv‖0≤Cηk+1‖v‖k+1, ‖v?Iηv‖1≤Cηk‖v‖k+1, ‖q?Jηq‖0≤Cηk‖q‖k, 其中C是與η無(wú)關(guān)的正常數(shù).

(A2) LBB條件: 存在與η無(wú)關(guān)的正常數(shù)β使得, 其中η=H或h.

方程(3)的混合有限元格式是: 求(uη,pη)∈Xη×Mη使得

問(wèn)題(4)的存在性、唯一性及收斂性分析, 參見(jiàn)文獻(xiàn)[7].

引理2[6?7,15]設(shè)(u,p),(uη,pη)(η=h或H)分別是方程(3)和(4)的解,對(duì)于有:

下面給出兩重網(wǎng)格法:

步驟1: 求(uH,pH)∈(XH,MH), 使得對(duì)于?(v,q)∈(XH,MH)有

步驟2: 求(uh,ph)∈(Xh,Mh), 使得對(duì)于?(v,q)∈(Xh,Mh)有

步驟3: 求(eH,εH)∈(XH,MH), 使得對(duì)于?(v,q)∈(XH,MH)有

其中：u?=uh+eH, p?=ph+εH.

為了便于分析, 將內(nèi)積空間定義為:Yη=(Xη,Mη)(η=H或h)和Y=(X,M), 定義范數(shù)為, 并引入連續(xù)雙線性形式AH:Y×Y→R和BH:Y×Y→R使得

則兩重網(wǎng)格法的步驟2和步驟3可以改寫(xiě)為:AH和BH在Y×Y上連續(xù)且依賴于|uH|1. 當(dāng)H足夠小, 可以證明|uH|1≤|uH?u|1+|u|1≤CH(|u|2+|p|1)+|u|1≤C(u,p). 并且AH和BH滿足[15]:

引入Galerkin投影(Q,R): Y→YH滿足BH[(w,t);(v?Q(v,q),q?R(v,q))]=0,?(w,t)∈YH,(v,q)∈Y.

為了導(dǎo)出v?Q(v,q)的L2估計(jì), 令(φ,?)∈Y滿足

因?yàn)閡是非奇異解, 所以(φ,?)存在唯一. 若(10)是H2正則的, 則對(duì)?f∈(L2(?))n, 解

引理3若u是非奇異解, 則(Q,R)滿足:

證明(I) 利用(9)式和BH的連續(xù)性得：

然后利用三角不等式證明.

(II) 令(φf(shuō),?f)為(10)式的解, 對(duì)于右邊f(xié)=v?Q(v,q)∈(L2(?))n. 通過(guò)計(jì)算

設(shè)w=v?Q(v,q),t=q?R(v,q), 然后利用投影性質(zhì), (2)式, 引理1, 引理2和嵌入定理得：

再利用上述不等式和插值不等式得:

(III) 當(dāng)(v,q)∈Y∩((H2(?))n×H1(?))時(shí), 利用(IHv,JHv)∈YH和(8)式得：

然后, 利用三角不等式和插值的性質(zhì)得：

再重復(fù)以上步驟并利用(11)式得:

最后, 利用插值不等式和(12)式證明結(jié)論.

3 能量誤差估計(jì)

定理1假設(shè)XH?Xh?X,MH?Mh?M滿足逼近性假設(shè)及LBB條件. 記(u,p)為(3)式的解且則有:

其中:u?=uh+eH,p?=ph+εH,C為常數(shù). 當(dāng)n=2時(shí),s=?>0(任意小); 當(dāng)n=3時(shí),

證明首先, 考慮步驟2之后的誤差, 誤差滿足:

由(13)式及AH的連續(xù)性得:

根據(jù)三角不等式, (A1), (2)式, (8)式, (14)式, 嵌入定理, 引理1和引理2得:

利用(3)式和(6)式得:

然后, 由(3)式, (7)式和(Q,R)的定義得:

(17)式減去(16)式, 且u?=uh+eH,p?=ph+εH得:

再由三角不等式, (A1)和(9)式得：

記

然后, 利用(2)式, 嵌入定理, 引理1, 引理2和引理3分別估計(jì)Sj(j=1,···,11):

最后, 將(19)式代入(18)式證明定理.

4 數(shù)值算例

在這部分進(jìn)行數(shù)值實(shí)驗(yàn). 第一個(gè)例子是用光滑精確解驗(yàn)證算法的收斂階. 第二個(gè)例子是用方腔驅(qū)動(dòng)流證明算法的穩(wěn)定性. 另外, 我們將比較標(biāo)準(zhǔn)有限元算法和兩重網(wǎng)格法的CPU時(shí)間. 所有的算法都是用有限元軟件FreeFem++[16]實(shí)現(xiàn).

4.1 收斂測(cè)驗(yàn)

問(wèn)題(1)的精確的解:

其中f是將精確解代入(1)式所得.

取?=[0,1]×[0,1],h=1/49,1/64,1/81,1/100,1/121, 采用P2?P1元時(shí),取r=3,α=1,ν=0.1; 采用P1b?P1元時(shí),取r=4,α=10,ν=1. 數(shù)值結(jié)果分別見(jiàn)表1～表6.

表1 利用P2 ?P1元的標(biāo)準(zhǔn)有限元算法的數(shù)值結(jié)果Tab 1 The numerical results of the standard finite element method using P2 ?P1 element

表2 利用P2 ?P1元的兩重網(wǎng)格法的數(shù)值結(jié)果Tab 2 The numerical results of the two-grid method using P2 ?P1 element

表3 利用P2 ?P1元的兩種方法的CPU時(shí)間(單位: s)Tab 3 The comparison of the CPU times (unit: s) for two methods using P2 ?P1 element

表4 利用P1b?P1元的標(biāo)準(zhǔn)有限元算法的數(shù)值結(jié)果Tab 4 The numerical results of the standard finite element method using P1b?P1 element

表5 利用P1b?P1元的兩重網(wǎng)格法的數(shù)值結(jié)果Tab 5 The numerical results of the two-grid method using P1b?P1 element

表6 利用P1b?P1元的兩種方法的CPU時(shí)間(單位: s)Tab 6 The comparison of the CPU times (unit: s) for two methods using P1b?P1 element

由表1、表2得: 采用P2?P1元時(shí)所有的收斂階完全符合理論分析;由表4、表5知: 采用P1b?P1元時(shí)速度的收斂階可以達(dá)到理論值, 而壓力的收斂階比理論分析高. 此外, 兩種方法的收斂階幾乎相同, 但兩重網(wǎng)格法比標(biāo)準(zhǔn)有限元方法節(jié)省了60%以上的CPU時(shí)間, 見(jiàn)表3和表6, 所以回溯兩重網(wǎng)格法更有效.

4.2 方腔驅(qū)動(dòng)流

在區(qū)域?=[0,1]×[0,1]上進(jìn)行計(jì)算, 在沒(méi)有其他體力的情況下, 流動(dòng)是由施加在頂部邊界上的切向速度場(chǎng)驅(qū)動(dòng)的. 在??上施加速度的法向分量為零, 切向分量除沿頂部的邊界設(shè)為1其余為零. 采用P2?P1元取h=1/25及不同的r,α,ν比較兩種方法的CPU時(shí)間, 見(jiàn)表7. 然后取r=3,α=10,ν=1和r=4,α=10,ν=0.1分別描述兩種方法的流線和壓力等高線, 如圖1～圖4所示.

表7 兩種方法的CPU(時(shí)間: s)比較Tab 7 The comparison of the CPU times (unit: s) for two methods

通過(guò)圖1～圖4和表7可以觀察到: 這兩種方法的流線和壓力等值線幾乎相同, 但兩重網(wǎng)格法比標(biāo)準(zhǔn)有限元方法節(jié)省時(shí)間. 因此, 回溯兩重網(wǎng)格法是穩(wěn)定的并且比標(biāo)準(zhǔn)有限元方法更有效.

圖1 標(biāo)準(zhǔn)有限元法的流線（左）和壓力等值線（右）Fig 1 The streams (left) and pressure contours (right) of the standard finite element method

圖2 兩重網(wǎng)格法的流線（左）和壓力等值線（右）Fig 2 The streams (left) and pressure contours (right) of the two-grid method

圖3 標(biāo)準(zhǔn)有限元法的流線（左）和壓力等值線（右）Fig 3 The streams (left) and pressure contours (right) of the standard finite element method

圖4 兩重網(wǎng)格法的流線（左）和壓力等值線（右）Fig 4 The streams (left) and pressure contours (right) of the two-grid method

5 結(jié)論

本文討論了不可壓縮的Brinkman-Forchheimer方程的回溯兩重網(wǎng)格法. 通過(guò)兩個(gè)數(shù)值實(shí)驗(yàn)對(duì)理論分析進(jìn)行驗(yàn)證，可知該算法是穩(wěn)定的, 兩種方法具有相同的收斂階, 但是在保持相同的收斂階時(shí)兩重網(wǎng)格法可以節(jié)省大量的時(shí)間, 因此回溯兩重網(wǎng)格法更有效.