對流擴散方程的保正守恒特征有限體積法研究?

李昕姝， 付 凱

(中國海洋大學數學科學學院， 山東 青島 266100)

對流擴散方程廣泛應用于眾多科學技術領域，如大氣環境計算，油藏開發數值模擬，地下水污染防控和金融計算等領域[1-5]。由于該方程通常不存在精確解，因此開發其數值求解方法進行計算模擬具有十分重要的實際應用價值。

對流占優擴散方程的數值求解是一項極具挑戰性的工作。在網格尺寸和時間步長不夠充分小時，傳統的有限差分、有限體積或有限元方法通常不能得到較好的結果，在陡峭前沿會出現虛假的非物理振蕩或過度的數值耗散。

為了克服求解對流占優擴散問題的困難，學者們提出了許多基于特征線方法的數值格式。從物理角度，特征線方法可以有效地求解沿流線問題，得到更為精確的近似結果；從計算角度，特征線方法可以使用較大的時間步長，因此可大幅降低計算成本[6]。Douglas和Russell提出了用于解決一維對流占優擴散問題的特征線修正方法[7]。隨后，Bermejo，Dawson和Hansbo等進一步發展了時間一階特征線方法用于解決高維對流擴散問題[8-12]，Fu，Liang和Rui等提出了時間二階數值格式[6,13-15]。

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

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

1 數值格式

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

考慮如下二維對流擴散方程

(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)

對二維區域Ω=[ax,bx]×[ay,by]，給出如下剖分

ax=x1/2ay=y1/2

(2)

對于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)

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

應用特征線法[7]，令

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

(7)

用τ表示與算子?tu+β1?xu+β2?yu相關聯的特征方向，

(8)

過tm時刻x=(x,y)點的特征線滿足常微分方程

(9)

(10)

每個Ωi,j在t時刻對應的特征單元記為

t∈[tm-1,tm]},

(11)

在tm-1時刻的特征單元記為

(12)

定義時空單元

(13)

(14)

通過沿特征線取平均來處理擴散項，給出如下時間二階近似：

(15)

(16)

(17)

對于在tm時刻的擴散項，定義

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)

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

(22)

綜上可得二維對流占優擴散問題(1)的時間二階格式

(23)

令

(24)

(25)

初始條件

(26)

邊界條件

(27)

1.2 守恒保正插值

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

第一步沿y=yj，j=1,…,Ny進行計算。定義守恒插值 [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表示在區間[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插值方法保正限制器

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

IfUmin<0，then

else

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

else

end

end

end

1.4 守恒保正局部限制器

假設1速度β滿足如下條件

(37)

(38)

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

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

即得證。

(39)

定義如下限制器:

(40)

其中

(41)

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

(42)

則

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

② 由式(42)，可得

2 格式性質

(43)

其矩陣形式為

AmUm=U*，

(44)

其中

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

(45)

且至少有一個不等式是嚴格大于零，則矩陣A是一個M矩陣，A-1中的所有元素均為非負的。

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

(46)

② 質量守恒:

(47)

則

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

因此

質量守恒性得證。

3 數值算例

在本節中，通過數值算例驗證了二維對流擴散問題時間二階保正守恒特征有限體積法(T2-PC-CFVM)的精度和性質。定義tm時刻數值解的離散質量誤差

(48)

用于評估數值格式的守恒性。

算例1首先考慮區域Ω=[-1.5,1.5]×[-1.5,1.5]上的Gaussian hill 旋轉問題。速度場β=(4y,-4x)T；擴散系數α1=α2=1×10-3。問題的精確解為

(49)

其中

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

(50)

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

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

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

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

關于保正性，由表2列出不同結果的最小值可以看出T2-PC-CFVM得到非負結果,而C-FDM-QI會產生負值，例如，當空間步長為h=1/50的時，C-FDM-QI格式的最小值為-1.170 4×10-3，T2-PC-CFVM格式的最小值為0。圖2給出了在T=π/4時刻，應用空間網格N=100和時間網格數M=100得到的精確解，以及C-FDM-QI和T2-PC-CFVM數值解的側視圖?？梢钥闯龈袷絋2-PC-CFVM的結果均為非負值，而格式C-FDM-QI產生了明顯的負值。

圖2 精確解和不同格式所得數值解的側視圖比較

同時還檢驗了這兩種方法的守恒性。表1和表2展示了計算結果的質量誤差。結果表明格式T2-PC-CFVM可以很好的保持守恒性，而格式 C-FDM-QI 的質量誤差超過了3×10-4。圖3展示了格式C-FDM-QI和T2-PC-CFVM的計算結果離散總質量隨時間的變化。可以看出在T=π/4時刻，與初始總質量相比，C-FDM-QI丟失了大約20% 的質量，而T2-PC-CFVM精確的保持了質量守恒。

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

算例2本算例中對陡峭前沿移動問題進行了計算，計算區域為Ω=[-1,1]2。初始條件為

(51)

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

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

4 結語

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