張甜甜,許文文,郭紅,杜明洋,余昌彪
(齊魯工業大學(山東省科學院) 數學與統計學院,山東 濟南 250301)
變系數對流擴散方程對研究自然界和實際工程中的很多物理現象具有重要意義,且在眾多學科中的應用極其廣泛,如流體力學、氣體力學等[1-2]。在研究此類問題時,擴散項的計算相對簡單,對流方程的計算通常是研究重點。
一階變系數對流方程在自然科學領域的應用背景較廣泛。王國棟等[3]將迎風格式應用到一階變系數對流方程的推廣模型交通流中,并通過新的格式模擬相關實例;Chen 等[4]通過耗散譜方法求解一階變系數雙曲方程(對流方程是最簡單的雙曲方程),通過傅里葉耗散譜方法討論周期問題,通過Legendre耗散譜方法討論邊界問題,且列出嚴格的誤差估計式; Aguirre等[5]在以Hermit函數為正交基的Herbert空間里開展對一階變系數雙曲方程周期邊界問題的研究,通過理論分析給出了收斂階 。關于對流方程的傳統數值解差分格式有迎風格式、Lax-Friedrichs格式、Lax-Wendroff格式、Beam-Warming格式和蛙跳格式。其中,迎風格式的基本思想是用向前差商或向后差商來代替空間偏導數,其關于時間和空間都是一階的,且是條件穩定的,算法在計算機上便于實現,因此在實際計算中受到廣泛重視,并得到了一些好的方法和技巧。本文我們將采用迎風格式逼近一階變系數對流方程并探索格式的穩定性條件[6-9]。
常系數問題差分格式的穩定性分析一般采用傅里葉方法和直接方法[9-11],但是變系數問題的穩定性分析相對復雜,通常不采用上述兩種方法,目前國內外關于一階變系數對流方程穩定性的相關研究較少。對于變系數方程差分格式穩定性問題,能量不等式方法是一個嚴格且很有技巧的方法。用能量不等式方法討論差分格式穩定性是從穩定性的定義出發,通過一系列估計式完成的,這個方法是偏微分方程中常用的能量方法的離散模擬[9]。本文將采用能量不等式方法分析一階變系數對流方程迎風格式的穩定性條件,分別就變系數的正負取值情況并結合凍結系數法驗證網格比條件,進而推導出穩定性條件,最后進行數值模擬驗證理論分析的正確性。
對于簡單的一階線性變系數對流方程的初值問題:
(1.1)
如果a(x,t)對x和t都是一次連續可微的,那么a(x,t)光滑變化,與常系數情形類似,(1.1)式的特征線滿足的方程為:
(1.2)
令x=x(t,x0)和u(x,t)分別是方程(1.1)和方程(1.2)的解,則:
(1.3)
于是,方程(1.1)的解沿特征線為常數。此時特征線為曲線,且有:
u(x,t)=u0(x),x=x(t,x0)。
(1.4)
因此,我們將已學常系數方程推導的差分格式推廣到變系數方程(1.1)。
設初值問題(1.1)的解區域為[0,l]×[0,T],將此區域沿x軸和t軸方向進行矩形網格剖分,其中空間步長為Δx=h=l/J,時間步長為Δt=τ,網格點記為(xj,tn),網格線可寫作:
xj=jΔx=jh,j=0,1,2,…,J,
tn=nΔt=nτ,n=0,1,2,…。
由于變系數對流方程(1.1)的迎風差分格式是對流方程關于空間偏導數用在特征線方向一側的單邊差商來代替的,且其系數a(x,t)符號是變化的,因此迎風格式可以寫成如下兩種形式:
(2.1)
(2.2)

|λj(G(τ,k))|≤1+Mτ,j=1,2,…,p,
其中|λj(G(τ,k))|表示增廣矩陣G(τ,k)的特征值,M為常數[9]。此條件稱為von Neumann條件。



(3.1)

vn+1eikjh=vneikjh-aλ(vneikjh-vneik(j-1)h),
(3.2)
兩邊同時消去公因子eikjh得:
vn+1=vn[1-aλ(1-e-ikh)],
(3.3)
所以增長因子為:
(3.4)
則有
|G(τ,k)|2=[1-aλ(1-coskh)]2+a2λ2sin2kh
(3.5)

接下來我們分情況討論迎風格式穩定性:

(3.6)

(3.7)
網格比滿足條件:
(3.8)
根據基本不等式a2+b2≥2ab并進一步化簡有:
用h乘上式兩邊并對j求和,記離散范數:
(3.9)
那么有:
(3.10)

(3.11)

(3.12)
從而有
(3.13)
重復使用上式有:
(3.14)


(3.15)



同理可知網格比滿足如下條件:
(3.16)
根據基本不等式a2+b2≥2ab并進一步化簡有:


用h乘上式兩邊并對j求和,記離散范數:
(3.17)
那么有:
(3.18)

(3.19)

(3.20)
從而有
(3.21)
重復使用上式有:
(3.22)

現在通過一個簡單的數值例子驗證用能量不等式的方法分析變系數對流方程的初值問題的迎風格式的穩定性條件。
對于初值問題:

u(x,t)=x2et。
(4.1)
取空間步長h=0.01, 時間步長τ=0.01,則λ=1,將迎風格式計算到tn=0.1時,計算得到初值問題的迎風格式的數值解與解析解,然后將數值解與解析解進行對比判斷此初值問題的穩定性,參考文獻[12],對隨機選取的部分輸出結果進行對比,見表1。

表1 數值解與解析解
從表中數據可知,迎風格式的數值解與解析解的誤差很小,即可認為此初值問題是穩定的。
圖1和圖2為輸出迎風格式的數值解圖像與解析解圖像,x軸代表空間方向的長度,y軸代表時間方向的長度,圖1的z軸是所研究初值問題的迎風格式關于x軸和y軸的數值解u(x,t),圖2中的z軸是所研究初值問題關于x軸和y軸的解析解z(x,t)。從迎風格式的數值解與解析解圖像可以直觀地看出此初值問題的迎風格式是穩定的。

圖1 迎風格式數值解Fig.1 Upwind scheme numerical solution

圖2 解析解Fig.2 Analytical solution
本文采用迎風格式逼近變系數對流方程,然后通過凍結系數法求出一階變系數對流方程的迎風格式穩定需要滿足的網格比條件。再對其迎風差分格式通過能量不等式的方法并結合凍結系數法得出的網格比條件進行討論得出其穩定性條件。最后通過一個數值算例對比其解析解與數值解,驗證穩定性條件的正確性。