對(duì)流擴(kuò)散方程的保正守恒特征有限體積法研究?

李昕姝， 付 凱

(中國(guó)海洋大學(xué)數(shù)學(xué)科學(xué)學(xué)院， 山東 青島 266100)

對(duì)流擴(kuò)散方程廣泛應(yīng)用于眾多科學(xué)技術(shù)領(lǐng)域，如大氣環(huán)境計(jì)算，油藏開(kāi)發(fā)數(shù)值模擬，地下水污染防控和金融計(jì)算等領(lǐng)域[1-5]。由于該方程通常不存在精確解，因此開(kāi)發(fā)其數(shù)值求解方法進(jìn)行計(jì)算模擬具有十分重要的實(shí)際應(yīng)用價(jià)值。

對(duì)流占優(yōu)擴(kuò)散方程的數(shù)值求解是一項(xiàng)極具挑戰(zhàn)性的工作。在網(wǎng)格尺寸和時(shí)間步長(zhǎng)不夠充分小時(shí)，傳統(tǒng)的有限差分、有限體積或有限元方法通常不能得到較好的結(jié)果，在陡峭前沿會(huì)出現(xiàn)虛假的非物理振蕩或過(guò)度的數(shù)值耗散。

為了克服求解對(duì)流占優(yōu)擴(kuò)散問(wèn)題的困難，學(xué)者們提出了許多基于特征線方法的數(shù)值格式。從物理角度，特征線方法可以有效地求解沿流線問(wèn)題，得到更為精確的近似結(jié)果；從計(jì)算角度，特征線方法可以使用較大的時(shí)間步長(zhǎng)，因此可大幅降低計(jì)算成本[6]。Douglas和Russell提出了用于解決一維對(duì)流占優(yōu)擴(kuò)散問(wèn)題的特征線修正方法[7]。隨后，Bermejo，Dawson和Hansbo等進(jìn)一步發(fā)展了時(shí)間一階特征線方法用于解決高維對(duì)流擴(kuò)散問(wèn)題[8-12]，F(xiàn)u，Liang和Rui等提出了時(shí)間二階數(shù)值格式[6,13-15]。

科學(xué)家和工程師們通常希望數(shù)值解能夠準(zhǔn)確地反映現(xiàn)實(shí)問(wèn)題的物理性質(zhì)，如保正性和守恒性等[15,16-21]。眾所周知，特征線方法需要在對(duì)前一時(shí)間層上追蹤點(diǎn)或追蹤單元處的解進(jìn)行近似。Colella和 Woodward關(guān)于對(duì)流方程的求解提出了一種分段拋物方法(PPM)，該方法滿足局部質(zhì)量守恒[22]。將特征線法與守恒插值結(jié)合，F(xiàn)u和Liang提出了求解對(duì)流擴(kuò)散問(wèn)題的高階守恒方法[6,13,15]，但未對(duì)保正性進(jìn)行研究。對(duì)于解的保正性和守恒性，通過(guò)引入限制器進(jìn)行修正是一種有效的實(shí)現(xiàn)方法。Liu，Zhang和Zhu等提出的局部限制器可以在保持格式二階精度的基礎(chǔ)上滿足保正性和質(zhì)量守恒[23-24]。

本文提出了對(duì)流占優(yōu)擴(kuò)散方程保正守恒特征有限體積法。通過(guò)將特征線方法和有限體積法相結(jié)合，對(duì)前一時(shí)間層追蹤單元上的解采用帶保正約束的二階精度守恒分段拋物插值近似，并運(yùn)用沿特征線取平均和局部限制器，得到時(shí)間二階保正守恒特征有限體積方法。我們給出了保正性和質(zhì)量守恒性的理論證明。對(duì)旋轉(zhuǎn)速度場(chǎng)中的Gaussian hill的輸運(yùn)問(wèn)題和陡峭前沿移動(dòng)問(wèn)題進(jìn)行了數(shù)值實(shí)驗(yàn)，驗(yàn)證了方法在時(shí)間和空間上的精度和保正性。通過(guò)計(jì)算離散質(zhì)量誤差，證明了格式具有質(zhì)量守恒。

1 數(shù)值格式

1.1 保正守恒特征有限體積法

考慮如下二維對(duì)流擴(kuò)散方程

(x,y,t)∈Ω×(0,T];

α▽u(x,y,t)·nΩ=0, (x,y,t)∈?Ω×(0,T];

u(x,y,0)=u0(x,y), (x,y)∈Ω。

(1)

對(duì)二維區(qū)域Ω=[ax,bx]×[ay,by]，給出如下剖分

ax=x1/2ay=y1/2

(2)

對(duì)于i=1,2,…,Nx，j=1,2,…,Ny，定義歐拉單元，單元中心和單元大小為

Ωi,j=[xi-1/2,xi+1/2]×[yj-1/2,yj+1/2],

(3)

(4)

(5)

令

(6)

對(duì)于時(shí)間離散，取均勻剖分。令Δt=T/M，tm=mΔt，m=0,1,2,…,M。

應(yīng)用特征線法[7]，令

ψ(x,y)=[1+β1(x,y)2+β2(x,y)2]1/2，

(7)

用τ表示與算子?tu+β1?xu+β2?yu相關(guān)聯(lián)的特征方向，

(8)

過(guò)tm時(shí)刻x=(x,y)點(diǎn)的特征線滿足常微分方程

(9)

(10)

每個(gè)Ωi,j在t時(shí)刻對(duì)應(yīng)的特征單元記為

t∈[tm-1,tm]},

(11)

在tm-1時(shí)刻的特征單元記為

(12)

定義時(shí)空單元

(13)

(14)

通過(guò)沿特征線取平均來(lái)處理擴(kuò)散項(xiàng)，給出如下時(shí)間二階近似：

(15)

(16)

(17)

對(duì)于在tm時(shí)刻的擴(kuò)散項(xiàng)，定義

v1(x,y,t)=-α1(x,y)ux(x,y,t)，
v2(x,y,t)=-α2(x,y)uy(x,y,t)。

(18)

則

(19)

其中

(20)

(21)

由于追蹤單元的邊界通常不與歐拉網(wǎng)格線平行，因此對(duì)tm-1時(shí)刻的擴(kuò)散項(xiàng)進(jìn)行高精度計(jì)算較為困難。為得到高精度質(zhì)量守恒數(shù)值解，應(yīng)用分片拋物函數(shù)在Ωi,j上近似Um-1來(lái)計(jì)算追蹤單元邊界擴(kuò)散通量[6]。該方案可保證tm-1時(shí)刻擴(kuò)散通量的計(jì)算精度及其在追蹤單元邊界的連續(xù)性。將其近似為

(22)

綜上可得二維對(duì)流占優(yōu)擴(kuò)散問(wèn)題(1)的時(shí)間二階格式

(23)

令

(24)

(25)

初始條件

(26)

邊界條件

(27)

1.2 守恒保正插值

(細(xì)黑線代表歐拉線，紅色線代表拉格朗日線，虛線代表中線。紅色單元是與灰色歐拉單元Ωi,j相對(duì)應(yīng)的拉格朗日追蹤單元 associated with gray Eulerian cell Ωi,j.)

第一步沿y=yj，j=1,…,Ny進(jìn)行計(jì)算。定義守恒插值 [Hm-1]p,j，p=1,…,Nx，

[Hm-1]p,j(x)=p-1/2,j+

(28)

滿足如下條件

(29)

p,j,L和p,j,R值由如下二階近似給出

(30)

(31)

[Hm-1]p,j表示在區(qū)間[xp-1/2,xp+1/2]上沿歐拉線y=yj方向的守恒插值分布。

U0,j=U1,j，
UNx+1,j=UNx,j。

(32)

(33)

(34)

滿足如下條件

(35)

(36)

1.3 PPM插值方法保正限制器

為使得所構(gòu)造出的 PPM 插值函數(shù)(見(jiàn)公式(28))滿足保正性，引入如下約束：

IfUmin<0，then

else

If p,j,R>p,j,L，then

else

end

end

end

1.4 守恒保正局部限制器

假設(shè)1速度β滿足如下條件

(37)

(38)

證明 對(duì)式(24)從i=1到Nx和j=1到Ny求和

由于在tm-1時(shí)刻擴(kuò)散通量在追蹤單元邊界的連續(xù)性，對(duì)于相鄰的兩個(gè)拉格朗日單元，通過(guò)一個(gè)邊界流入單元的量等于其他單元通過(guò)該邊界流出的量[6]，并注意邊界條件(公式(27))可得

即得證。

(39)

定義如下限制器:

(40)

其中

(41)

定理1限制器(40)具有如下性質(zhì)：

(42)

則

證明 ① 將式(41)代入式(40)，并注意M的定義，可直接得出保正性。

② 由式(42)，可得

2 格式性質(zhì)

(43)

其矩陣形式為

AmUm=U*，

(44)

其中

引理2[27-28]對(duì)于一個(gè)不可約矩陣A=(aij)N×N滿足aii>0(1≤i≤N)和aij≤0(1≤i,j≤N,i≠j)，如果A為行弱對(duì)角占優(yōu)，即

(45)

且至少有一個(gè)不等式是嚴(yán)格大于零，則矩陣A是一個(gè)M矩陣，A-1中的所有元素均為非負(fù)的。

定理2在假設(shè)1和Neumann邊界條件下，格式 T2-PC-BFVM(見(jiàn)公式(25))具有如下性質(zhì)：

(46)

② 質(zhì)量守恒:

(47)

則

注意邊界條件(27)，可得

因此

質(zhì)量守恒性得證。

3 數(shù)值算例

在本節(jié)中，通過(guò)數(shù)值算例驗(yàn)證了二維對(duì)流擴(kuò)散問(wèn)題時(shí)間二階保正守恒特征有限體積法(T2-PC-CFVM)的精度和性質(zhì)。定義tm時(shí)刻數(shù)值解的離散質(zhì)量誤差

(48)

用于評(píng)估數(shù)值格式的守恒性。

算例1首先考慮區(qū)域Ω=[-1.5,1.5]×[-1.5,1.5]上的Gaussian hill 旋轉(zhuǎn)問(wèn)題。速度場(chǎng)β=(4y,-4x)T；擴(kuò)散系數(shù)α1=α2=1×10-3。問(wèn)題的精確解為

(49)

其中

(x0,y0)=(-0.5,0)，σ=0.1。

(50)

初始條件由式(49)中t=0給出。

數(shù)值試驗(yàn)中將新格式與Douglas和Russell[7]提出的雙二次插值標(biāo)準(zhǔn)特征有限差分(Characteristic finite difference method with biquadratic interpolation，C-FDM-QI)進(jìn)行了比較，計(jì)算了兩種方法的收斂階。取T=π/4，擴(kuò)散系數(shù)α1=α2=1×10-3和均勻網(wǎng)格。選擇空間步長(zhǎng)h=2.5×10-3，以使得空間誤差足夠小，選擇時(shí)間網(wǎng)格點(diǎn)數(shù)M=20、24、28、32和36用于計(jì)算時(shí)間收斂階。表1列出了不同數(shù)值格式誤差的L2范數(shù)(E2)和時(shí)間階。從表1的結(jié)果中可以很明顯的看出，T2-PC-CFVM格式具有時(shí)間二階精度，而C-FDM-QI為一階精度。選擇時(shí)間網(wǎng)格數(shù)M=2 000和空間步長(zhǎng)h=1/50、1/60、1/70、1/80 和1/90，對(duì)空間階進(jìn)行了計(jì)算。如表2所示，兩種格式均為二階精度。

表1 格式C-FDM-QI和T2-PC-CFVM所得解的誤差、時(shí)間收斂階、最小值和質(zhì)量誤差

表2 格式C-FDM-QI和T2-PC-CFVM所得解的誤差、空間收斂階、最小值和質(zhì)量誤差

關(guān)于保正性，由表2列出不同結(jié)果的最小值可以看出T2-PC-CFVM得到非負(fù)結(jié)果,而C-FDM-QI會(huì)產(chǎn)生負(fù)值，例如，當(dāng)空間步長(zhǎng)為h=1/50的時(shí)，C-FDM-QI格式的最小值為-1.170 4×10-3，T2-PC-CFVM格式的最小值為0。圖2給出了在T=π/4時(shí)刻，應(yīng)用空間網(wǎng)格N=100和時(shí)間網(wǎng)格數(shù)M=100得到的精確解，以及C-FDM-QI和T2-PC-CFVM數(shù)值解的側(cè)視圖。可以看出格式T2-PC-CFVM的結(jié)果均為非負(fù)值，而格式C-FDM-QI產(chǎn)生了明顯的負(fù)值。

圖2 精確解和不同格式所得數(shù)值解的側(cè)視圖比較

同時(shí)還檢驗(yàn)了這兩種方法的守恒性。表1和表2展示了計(jì)算結(jié)果的質(zhì)量誤差。結(jié)果表明格式T2-PC-CFVM可以很好的保持守恒性，而格式 C-FDM-QI 的質(zhì)量誤差超過(guò)了3×10-4。圖3展示了格式C-FDM-QI和T2-PC-CFVM的計(jì)算結(jié)果離散總質(zhì)量隨時(shí)間的變化?？梢钥闯鲈赥=π/4時(shí)刻，與初始總質(zhì)量相比，C-FDM-QI丟失了大約20% 的質(zhì)量，而T2-PC-CFVM精確的保持了質(zhì)量守恒。

圖3 格式C-FDM-QI和T2-PC-CFVM所得解的離散總質(zhì)量隨時(shí)間的變化

算例2本算例中對(duì)陡峭前沿移動(dòng)問(wèn)題進(jìn)行了計(jì)算，計(jì)算區(qū)域?yàn)棣?[-1,1]2。初始條件為

(51)

速度場(chǎng)(β1(x,y),β2(x,y))=(-0.5,-0.5)，擴(kuò)散系數(shù)α1=α2=1×10-3。圖4展示了M=50，N=100，T=0.8時(shí)，不同數(shù)值格式 C-FDM-QI，T2-PC-CFVM和近似精確解的比較圖。取M=200和N=200的細(xì)網(wǎng)格，應(yīng)用格式T2-PC-CFVM得到近似精確解。由結(jié)果可以看出格式T2-PC-CFVM可以很好地近似陡峭前沿，得到的結(jié)果與參考精確解基本一致，T2-PC-CFVM的結(jié)果均為正值。C-FDM-QI 在陡坡前沿移動(dòng)附近產(chǎn)生強(qiáng)烈的震蕩，不能準(zhǔn)確地近似陡峭前沿且未能保正。

圖4 T=0.8 時(shí)刻的參考精確解和不同格式所得數(shù)值解(見(jiàn)算例2)

4 結(jié)語(yǔ)

本文給出了求解對(duì)流占優(yōu)擴(kuò)散方程的時(shí)空二階保正守恒格式。沿特征線平均得到時(shí)間二階格式。通過(guò)帶保正約束的守恒插值和局部限制器使得數(shù)值解為正值，保持質(zhì)量守恒且精度不變。數(shù)值實(shí)驗(yàn)證明了所提方法的性能，并從數(shù)值上證明了它們的守恒性和保正性。