陳國芳, 吳 丹, 呂俊良

(1. 吉林省教育學(xué)院 民族教育學(xué)院, 長春 130022; 2. 吉林大學(xué) 數(shù)學(xué)學(xué)院, 長春 130012)

0 引 言

考慮滲流方程：

其中:T<∞;m>1. 顯然, 方程(1)是一個非線性退化的拋物方程. 滲流方程常用于模擬物理上的多孔介質(zhì)或熱擴散等問題. 文獻[1]給出了標(biāo)準(zhǔn)有限元法求解該類問題的理論分析； 文獻[2]采用積分方程方法研究了退化的非線性拋物問題； 文獻[3]用移動網(wǎng)格方法研究了二維多孔介質(zhì)擴散問題； 文獻[4]給出了求解問題(1)的局部間斷Galerkin法(LDG); 文獻[5]提出了一種交錯的擬態(tài)差分法(MFDM)求解該類退化問題; 文獻[6-7]分別給出了改進的混合有限元法和混合有限體積元法求解退化非線性拋物問題. 由于節(jié)點型方法求解滲流方程時會出現(xiàn)非物理振蕩, 因此文獻[8]給出了3種后驗修補技術(shù); 文獻[9]建立了保正且守恒的二次元有限體積法. 文獻[5]研究表明, 當(dāng)使用有限體積法求解退化拋物問題時, 會產(chǎn)生數(shù)值曲線不能有效向前傳播的所謂“數(shù)值熱障”現(xiàn)象. 這主要是因為當(dāng)擴散系數(shù)等于零或趨于零時, 擴散系數(shù)在單元格上的調(diào)和平均也將等于零或趨于零, 因此數(shù)值求解時, 擴散過程不發(fā)生或幾乎不發(fā)生, 導(dǎo)致數(shù)值界面不會隨著時間的推移而向前傳播. 基于此, 本文考慮改進擴散系數(shù)的數(shù)值逼近方法, 取內(nèi)部邊兩側(cè)單元中心點處未知量的代數(shù)平均值作為非線性擴散系數(shù)中未知量的近似. 從而避免因一側(cè)未知量等于零或趨于零導(dǎo)致的“數(shù)值熱障”問題.

1 標(biāo)準(zhǔn)有限體積法

首先離散時間變量, 將時間[t0,T]做如下N等分:

t0

其中tn=t0+n(T-t0)/N. 記K(u)=mum-1, 使用向后Euler法離散時間導(dǎo)數(shù)項, 得半離散方程

(4)

其中: 時間步長Δt=(T-t0)/N;u(n)是u(x,t)在時刻tn的近似. 方程(4)可進一步改寫成

(5)

a=x0

記任意單元Ii=[xi-1,xi]的步長為hi(i=1,2,…,M), 剖分步長h=max{hi;i=1,2,…,M}.

1.1 中心型有限體積法

考慮方程(5), 在單元Ii+1=[xi,xi+1]上對方程(5)積分, 得

(6)

利用分部積分公式, 并對u(n+1)和u(n)采用中點近似, 有

(7)

(8)

(9)

將u(n+1)(xi+1)的表達式代入式(8), 可得流的近似為

顯然當(dāng)hi=hi+1時, 這里將導(dǎo)出擴散系數(shù)的調(diào)和平均作為非線性擴散系數(shù)的近似. 類似地, 還有

于是, 對于1≤i≤M-2, 格式為

即

(10)

類似地, 對于i=1, 格式為

(11)

對于i=M-1, 格式為

(12)

(A(n+1)+B)U(n+1)=F(n+1),n=0,1,…,N-1,

其中:

實際計算中, 需使用非線性迭代法, 如Picard迭代法. 因此線性方程組可寫成

(A(n+1,s)+B)U(n+1,s+1)=F(n+1),s=0,1,2,…,

其中:

這里

1.2 數(shù)值實例

問題(1)的精確解稱為Barenblatt解[10]:

其中u+=max{0,u}. 因此初值條件可取為g(x)=Bm(x,t0). 圖1為Barenblatt解的圖像. 使用1.1中推導(dǎo)的中心型有限體積法求解該問題的數(shù)值結(jié)果如圖2所示. 由圖2可見, 計算終止時刻數(shù)值波陣面未離開初始位置, 產(chǎn)生了數(shù)值振蕩. 因此, 標(biāo)準(zhǔn)單元中心型有限體積法不適合直接求解退化拋物問題.

圖1 當(dāng)t=0.1,1時問題(1)的Barenblatt解Fig.1 Barenblatt solutions of problem (1) when t=0.1,1

圖2 中心型有限體積法的數(shù)值結(jié)果Fig.2 Numerical results of cell-centered finite volume methods

2 修正的單元中心有限體積法

未知函數(shù)仍定義在單元中心, 但擴散系數(shù)的近似選取單元節(jié)點處的值. 對于非線性問題, 需要用單元中心點的值組合得到單元節(jié)點的值.

2.1 修正的中心型有限體積格式

考慮方程(7), 在xi+1點的流可用如下兩種方式近似:

(13)

(14)

將u(n+1)(xi+1)的表達式代入式(13)可得流的近似：

顯然這里導(dǎo)出擴散系數(shù)的近似不是調(diào)和平均值. 類似地, 還有

于是, 對于1≤i≤M-2, 格式為

即

(15)

類似地, 可得當(dāng)i=0時的格式為

(16)

當(dāng)i=M-1時的格式為

(17)

從而可得如下線性方程組:

(A(n+1)+IM×M)U(n+1)=F(n+1),n=0,1,…,N-1,

其中:

這里

實際計算中, 使用非線性迭代法所得的線性方程組可寫成

(A(n+1,s)+IM×M)U(n+1,s+1)=F(n+1),s=0,1,2,…,

其中:

這里

定理1格式(15)～(17)是保正的, 即若g(x)≥0, 則對每個n∈+和s∈+, 均有U(n+1,s+1)≥0.

2.2 修正的中心型有限體積法數(shù)值實例

考慮Barenblatt解問題, 使用上述推導(dǎo)的修正的中心型有限體積法求解, 數(shù)值結(jié)果如圖3所示. 由圖3可見, 新方法不僅沒有“數(shù)值熱障”現(xiàn)象, 還可以很好地捕捉波陣面且數(shù)值解保正.

圖3 修正的中心型有限體積法數(shù)值結(jié)果Fig.3 Numerical results of modified cell-centered finite volume methods

表1 幾種數(shù)值方法的L∞誤差比較

表2 幾種數(shù)值方法的L2誤差比較

表3 幾種數(shù)值方法保持能量守恒能力的比較

綜上, 針對使用標(biāo)準(zhǔn)有限體積法求解退化問題時, 會出現(xiàn)數(shù)值曲線不能向前傳播的所謂“數(shù)值熱障”現(xiàn)象, 本文提出了一種修正的有限體積法, 可有效解決該問題. 新方法中非線性擴散系數(shù)的取值為相鄰單元中心點處的代數(shù)平均, 從而避免了因擴散系數(shù)的調(diào)和平均近似帶來的缺陷.