和文強(qiáng),秦國良

(西安交通大學(xué)流體機(jī)械研究所,710049,西安)

?

譜元方法求解對流擴(kuò)散方程及其穩(wěn)定性分析

和文強(qiáng),秦國良

(西安交通大學(xué)流體機(jī)械研究所,710049,西安)

探討了一維對流擴(kuò)散方程的一種高精度數(shù)值解法,該解法在空間上采用了Chebyshev譜元方法,在時(shí)間上結(jié)合了半隱式Adams方法。通過數(shù)值算例驗(yàn)證了解法的可行性,利用特征分析法得到了對流擴(kuò)散方程譜元求解時(shí)不同離散形式的穩(wěn)定性條件,并對數(shù)值求解的穩(wěn)定性進(jìn)行了預(yù)測。通過時(shí)間步長、網(wǎng)格劃分、對流項(xiàng)和黏性項(xiàng)插值階數(shù)的影響研究表明:耦合Chebyshev譜元方法和半隱式Adams方法在求解對流擴(kuò)散方程時(shí)能夠獲得高精度的數(shù)值解;時(shí)間離散時(shí)Adams方法的黏性項(xiàng)采用一階插值形式、對流項(xiàng)采用二階插值形式,在未增加計(jì)算量的同時(shí)能夠獲得較大的穩(wěn)定區(qū)域和較高的計(jì)算精度。

對流擴(kuò)散方程;譜元法;穩(wěn)定性;半隱式Adams方法

對流擴(kuò)散方程是描述流體流動和傳熱、污染物輸移和擴(kuò)散、電化學(xué)反應(yīng)等的典型方程之一,利用數(shù)值方法進(jìn)行求解時(shí),由于對流項(xiàng)的存在,求解過程可能會出現(xiàn)不穩(wěn)定現(xiàn)象,因此數(shù)值方法的穩(wěn)定性問題一直是對流擴(kuò)散方程研究的熱點(diǎn)。當(dāng)對流擴(kuò)散方程的對流項(xiàng)占優(yōu)時(shí),應(yīng)用傳統(tǒng)的有限差分法和有限元法求解會出現(xiàn)數(shù)值振蕩[1-2],Gottlieb等在利用擬譜方法求解時(shí)也發(fā)現(xiàn)了不穩(wěn)定現(xiàn)象[3],而Mofid等通過引入罰方法來處理邊界條件能夠部分或完全避免這種不穩(wěn)定現(xiàn)象[4]。

譜元方法[5]結(jié)合了譜方法的高精度和有限元法的靈活性特點(diǎn),使得對流擴(kuò)散問題求解有望獲得高精度的數(shù)值解。目前,譜元方法應(yīng)用的主要障礙是高雷諾數(shù)Re下的穩(wěn)定性問題。因?yàn)樽V元方法在單元上需對求解變量進(jìn)行譜近似,當(dāng)Re增大、黏性耗散作用減小時(shí),任意微小擾動都會引起單元數(shù)值解不穩(wěn)定且最終擴(kuò)展到整個(gè)求解區(qū)域。為了擴(kuò)大Re范圍,一些穩(wěn)定性措施[6-7]引入到了譜元方法之中。

本文在前期譜元方法[8]研究的基礎(chǔ)上,采用特征分析法系統(tǒng)地分析了Chebyshev譜元方法的半離散形式和Chebyshev譜元方法耦合半隱式Adams時(shí)間離散方法的全離散形式的穩(wěn)定性條件,討論了半隱式Adams時(shí)間離散方法中對流項(xiàng)和黏性項(xiàng)在不同插值階數(shù)、網(wǎng)格劃分、時(shí)間步長時(shí)對穩(wěn)定性的影響,探索了最優(yōu)離散格式,試圖擴(kuò)大穩(wěn)定區(qū)域,提高計(jì)算精度和計(jì)算效率,最后結(jié)合具體算例對Chebyshev譜元方法耦合半隱式Adams方法在求解對流擴(kuò)散方程時(shí)的有效性進(jìn)行了驗(yàn)證。

1 對流擴(kuò)散方程

一維對流擴(kuò)散方程為

(1)

u(x,0)=u0(x)

(2)

u(0,t)=g1(t);u(1,t)=g2(t)

(3)

2 譜元方法

僅考慮Re影響,式(1)的Galerkin變分形式為

(4)

原微分方程的變分問題:求u∈H1使得

B(u,v)=f(v), ?v∈H1

(5)

(6)

插值函數(shù)為

(7)

式中:Tm、Tn為Chebyshev多項(xiàng)式;Nx為單元網(wǎng)格數(shù);cm滿足

(8)

在標(biāo)準(zhǔn)單元內(nèi)對式(4)進(jìn)行離散得

(9)

(10)

式(10)為u關(guān)于時(shí)間的一階線性常微分方程組。

3 時(shí)間離散

式(10)在時(shí)間上的逐步積分途徑一般有顯式和隱式兩種。顯式積分計(jì)算效率較高,需要的計(jì)算機(jī)內(nèi)存相對較小,但計(jì)算精度較低,條件穩(wěn)定;隱式積分的精度高,穩(wěn)定性較好,但每一個(gè)積分步都要求解一個(gè)線性方程組,不便于直接計(jì)算。本文時(shí)間項(xiàng)采用半隱式Adams方法進(jìn)行離散,對流項(xiàng)采用Je階顯式Adams-Bashforth形式,黏性項(xiàng)采用Ji階隱式Adams-Monlton形式[9],這樣便融合了顯式方法耦合簡單和隱式方法穩(wěn)定性好的特點(diǎn)。式(10)的全離散形式為

(11)

式中:βj為Adams-Bashforth系數(shù)[10];γk為Adams-Monlton系數(shù)[10]。

4 譜元方法的穩(wěn)定性分析

4.1 半離散方程的穩(wěn)定性

將式(10)兩側(cè)左乘B的逆矩陣可得

(12)

由于僅考慮第一類邊界條件,所以在邊值計(jì)算中不引入誤差,而假定在初值的計(jì)算中引入了擾動ε0,擾動ε滿足

(13)

式中:ε=(ε1,ε2,…,εn-1);L為D的內(nèi)配置點(diǎn)n-1階子矩陣。

根據(jù)Lyapunov穩(wěn)定性理論[11],式(12)穩(wěn)定的充要條件為矩陣L的一切特征根都有非正實(shí)部,且每個(gè)有0實(shí)部的特征根對應(yīng)矩陣L的約當(dāng)標(biāo)準(zhǔn)形中的一維約當(dāng)塊。也就是說,如果存在一個(gè)有正實(shí)部的特征根,式(12)對于任何時(shí)間離散都是不穩(wěn)定的。矩陣L的特征值實(shí)部最大值與Re、單元數(shù)、單元網(wǎng)格數(shù)的關(guān)系如圖1所示。

圖1 特征值實(shí)部最大值與單元數(shù)、單元網(wǎng)格數(shù)的關(guān)系

對于不同的Re和單元網(wǎng)格數(shù),圖1中2種單元劃分的空間譜元離散始終滿足半離散方程穩(wěn)定性要求,其他單元劃分的情形也可以進(jìn)行驗(yàn)證。

4.2 全離散方程的穩(wěn)定性

4.2.1 單步法全離散的穩(wěn)定性條件 采用半隱式Adams法時(shí)間離散,對流項(xiàng)、黏性項(xiàng)的插值階數(shù)均為1,即n+1時(shí)刻的變量僅與n時(shí)刻的變量相關(guān)時(shí),式(11)可表示為

un+1=φ(ΔtC)un+φf(ΔtB)fn+1

(14)

式中:φ、φf為時(shí)間離散函數(shù)。

擾動ε的全離散形式為

εn+1=Lεn

(15)

Sousa分析了形如式(14)的全離散格式的穩(wěn)定性條件[12],即:當(dāng)L為正規(guī)矩陣時(shí),其譜半徑ρ(L)≤1為穩(wěn)定的充分條件;當(dāng)L為非正規(guī)矩陣時(shí),其譜半徑ρ(L)≤1僅為穩(wěn)定的必要條件?？梢则?yàn)證,當(dāng)采用的Chebyshev函數(shù)為單元基函數(shù)時(shí),L為正規(guī)矩陣,所以單步譜元法穩(wěn)定的充分條件為ρ(L)≤1。

4.2.2 多步法全離散的穩(wěn)定性條件 采用半隱式Adams法時(shí)間離散,對流項(xiàng)或黏性項(xiàng)的插值階數(shù)大于1,即n+1時(shí)刻的變量與前n-k+1時(shí)刻的變量相關(guān)時(shí),式(11)可表示為

φj(ΔtDn-k+j)un-k+j+φf(ΔtB)fn+1

n=k,k+1,k+2,…

(16)

式中:φj、φf為時(shí)間離散函數(shù);Dn-k+j為un-k+j的系數(shù)矩陣。

擾動ε的全離散形式為

n=k,k+1,k+2,…

(17)

式中:Lj為φj(ΔtDn-k+j)的內(nèi)配置點(diǎn)n-1階子矩陣。

Dorsselaer等分析了形如式(16)的全離散格式的穩(wěn)定性條件[13],定義擾動

Vn+1=((εn+1)T,(εn)T,…,(εn-k+1)T)T

滿足

Vn+1=EVn,n≥k

(18)

式中:I為單位矩陣;0為零矩陣。

多步譜元法穩(wěn)定的充分條件為ρ(E)≤1。

4.3 不同形式全離散方程的穩(wěn)定區(qū)域

通過選用4種不同的網(wǎng)格劃分方式,分別表示為3×5(Nm=3,Nx=5)、5×7(Nm=5,Nx=7)、10×6(Nm=10,Nx=6)、10×10(Nm=10,Nx=10),時(shí)間(計(jì)算時(shí)間與特征時(shí)間的比值)步長Δt∈[0.001,0.1],來比較不同網(wǎng)格劃分、時(shí)間步長和時(shí)間離散方式的求解穩(wěn)定性。

對于半隱式Adams法時(shí)間離散,當(dāng)對流項(xiàng)采用一階插值形式,黏性項(xiàng)分別采用一階、二階插值形式時(shí),式(10)的全離散形式可分別表示為

(19)

(20)

將式(19)、式(20)分別表示為式(14)、式(16)的形式,即可由單步法和多步法的穩(wěn)定性條件對求解穩(wěn)定性進(jìn)行判定。對流項(xiàng)一階插值,黏性項(xiàng)一階、二階插值時(shí)的求解穩(wěn)定域如圖2所示。從圖2可以看出,黏性項(xiàng)二階插值的穩(wěn)定域和一階插值基本相同,網(wǎng)格劃分對穩(wěn)定域的影響很小,上臨界Re隨著時(shí)間步長的減小而逐步增大。

圖2 黏性項(xiàng)不同插值階數(shù)的穩(wěn)定性曲線

對于半隱式Adams法時(shí)間離散,對流項(xiàng)采用一階插值形式,黏性項(xiàng)采用三階插值形式時(shí),式(10)的全離散形式可表示為

(21)

將式(21)表示為式(16)的形式,利用多步法穩(wěn)定性條件對求解穩(wěn)定性進(jìn)行判定。對流項(xiàng)一階插值、黏性項(xiàng)三階插值時(shí)的求解穩(wěn)定域如圖3所示。從圖3可以看出,黏性項(xiàng)三階插值時(shí),不僅存在上臨界Re,也存在下臨界Re,只有當(dāng)Re位于上臨界Re和下臨界Re之間時(shí)求解才是穩(wěn)定的。上臨界Re隨時(shí)間步長的增大而減小,下臨界Re隨時(shí)間步長的增大而增大,同時(shí)求解穩(wěn)定域隨網(wǎng)格節(jié)點(diǎn)數(shù)的增加而減小。黏性項(xiàng)三階插值對求解的穩(wěn)定性的影響很大,且顯著減小了穩(wěn)定域的范圍。

圖3 黏性項(xiàng)三階插值的穩(wěn)定性曲線

對于半隱式Adams法時(shí)間離散,當(dāng)黏性項(xiàng)采用一階插值形式,對流項(xiàng)分別采用二階、三階插值形式時(shí),式(10)的全離散形式可表示為

(22)

(23)

將式(22)、式(23)表示為式(16)的形式,同樣利用多步法穩(wěn)定性條件對求解穩(wěn)定性進(jìn)行判定。黏性項(xiàng)一階插值,對流項(xiàng)一階、二階、三階插值時(shí)的求解穩(wěn)定域如圖4所示。從圖4可以看出:對流項(xiàng)二階插值存在臨界時(shí)間步長,上臨界Re在時(shí)間步長大于臨界值時(shí)有所減小,在時(shí)間步長小于臨界值時(shí)大幅度增加;對流項(xiàng)三階插值也存在臨界時(shí)間步長,上臨界Re在時(shí)間步長大于臨界值時(shí)比二階插值進(jìn)一步減小,在時(shí)間步長小于臨界值時(shí)迅速增加到無窮大。網(wǎng)格節(jié)點(diǎn)數(shù)對臨界時(shí)間步長的影響也很大,隨節(jié)點(diǎn)數(shù)的增加,臨界時(shí)間步長減小,對于10×10的網(wǎng)格劃分,二階、三階插值的臨界時(shí)間步長均小于0.001。

圖4 對流項(xiàng)不同插值階數(shù)的穩(wěn)定性曲線

5 計(jì)算實(shí)例和數(shù)值結(jié)果

(24)

統(tǒng)一選取時(shí)間步長Δt=0.001,比較3單元5網(wǎng)格和10單元10網(wǎng)格2種網(wǎng)格劃分,并給出t=1.0時(shí)的計(jì)算誤差。

5.1 單步法全離散數(shù)值結(jié)果

單步法全離散的數(shù)值結(jié)果如圖5所示。Re小于上臨界Re時(shí),2種網(wǎng)格劃分都取得了正確的數(shù)值結(jié)果,但由于時(shí)間離散階數(shù)較低,所以數(shù)值結(jié)果的精度不高。Re大于上臨界Re時(shí),10單元10網(wǎng)格的數(shù)值結(jié)果快速發(fā)散,3單元5網(wǎng)格的數(shù)值結(jié)果雖有發(fā)散的趨勢,但速度非常平緩,說明節(jié)點(diǎn)數(shù)越多,不穩(wěn)定作用越強(qiáng)。

圖5 單步譜元法的計(jì)算誤差

5.2 多步法全離散數(shù)值結(jié)果

對流項(xiàng)采用一階插值,黏性項(xiàng)采用二階、三階插值時(shí)的數(shù)值結(jié)果如圖6所示。對于黏性項(xiàng)二階插值,當(dāng)Re小于上臨界Re時(shí),計(jì)算誤差和一階插值相差很小;當(dāng)Re大于上臨界Re且節(jié)點(diǎn)數(shù)較多時(shí),數(shù)值結(jié)果的發(fā)散速率遠(yuǎn)大于節(jié)點(diǎn)數(shù)較少時(shí)的發(fā)散速率。黏性項(xiàng)三階插值的計(jì)算誤差和二階插值基本相同,當(dāng)Re大于上臨界Re或小于下臨界Re且節(jié)點(diǎn)數(shù)較多時(shí),數(shù)值結(jié)果快速發(fā)散。

圖6 黏性項(xiàng)二階、三階插值的計(jì)算誤差

黏性項(xiàng)采用一階插值,對流項(xiàng)采用二階、三階插值時(shí)的數(shù)值結(jié)果如圖7所示。對流項(xiàng)采用二階插值時(shí),計(jì)算誤差相較一階插值大大減小,同時(shí)Re較大時(shí)的誤差下降更多,Re大于上臨界Re且節(jié)點(diǎn)數(shù)較多時(shí),數(shù)值結(jié)果快速發(fā)散。對流項(xiàng)三階插值的計(jì)算精度與二階插值基本相同。

圖7 對流項(xiàng)二階、三階插值的計(jì)算誤差

6 結(jié) 論

本文采用Chebyshev譜元方法耦合半隱式Adams時(shí)間離散方法求解了一維對流擴(kuò)散方程,并對求解過程中方程的半離散和全離散形式的穩(wěn)定性條件進(jìn)行了討論,分析了半隱式Adams時(shí)間離散對流項(xiàng)、黏性項(xiàng)在不同插值階數(shù)、網(wǎng)格劃分和時(shí)間步長時(shí)對求解穩(wěn)定性及計(jì)算精度的影響,得出如下結(jié)論。

(1)對流擴(kuò)散方程的空間譜元離散,對于任意的Re,網(wǎng)格劃分都滿足半離散方程的穩(wěn)定性條件。

(2)隨著對流項(xiàng)插值階數(shù)的增加,求解穩(wěn)定域得到了擴(kuò)大,計(jì)算精度大大提高。對流項(xiàng)二階插值、時(shí)間步長小于臨界值時(shí),上臨界Re能取得較大的值;時(shí)間步長大于臨界值時(shí),上臨界Re減小幅度減緩,適宜在實(shí)際數(shù)值計(jì)算中應(yīng)用。

(3)隨著黏性項(xiàng)插值階數(shù)的增加,求解穩(wěn)定域縮小。三階插值時(shí)出現(xiàn)了下臨界Re,穩(wěn)定域大大減小。

今后的工作可以從以下兩方面展開:在方程中通過增加譜消除人工黏性項(xiàng),使得譜元方法能夠求解的穩(wěn)定區(qū)域擴(kuò)大;探索高精度的時(shí)間處理方法,如時(shí)空耦合的最小二乘譜元方法。

[1] GE L X, ZHANG J. High accuracy iterative solution of convection diffusion equation with boundary layers on nonuniform grids [J]. Journal of Computational Physics, 2001, 171(2): 560-578.

[2] ZIENKIEWICZ O C, HEINRICH J C. The finite element method and convection problems in fluid mechanics [M]. New York, USA: John Wiley & Sons Inc., 1978: 1-22.

[3] GOTTLIEB D, ORSZAG S A. Numerical analysis of spectral method: theory and application [M]. Philadelphia, PA, USA: SIAM, 1977: 139-143.

[4] MOFID A, PEYRET R. Stability of the Chebyshev collocation approximation to the advection-diffusion equation [J]. Computers Fluids, 1993, 22(4/5): 453-465.

[5] PETERA A T. A spectral element method for fluid dynamics: laminar flow in a channel expansion [J]. Journal of Computational Physics, 1984, 54(3): 468-488.

[6] FISCHER P, MULLEN J. Filter-based stabilization of spectral element methods [J]. Comptes Rendus de L’Académie des Sciences, 2001, 332(3): 265-270.

[7] XU C J, PSAQUETTI R. Stabilized spectral element computations of high Reynolds number incompressible flows [J]. Journal of Computational Physics, 2004, 196(2): 680-704.

[8] 秦國良, 徐忠. 譜元方法求解二維不可壓縮Navier-Stokes方程 [J]. 應(yīng)用力學(xué)學(xué)報(bào), 2000, 17(4): 20-26. QIN Guoliang, XU Zhong. A spectral element method for incompressible Navier-Stokes equations [J]. Chinese Journal of Applied Mechanics, 2000, 17(4): 20-26.

[9] 陳雪江, 秦國良, 徐忠. 譜元法和高階時(shí)間分裂法求解方腔頂蓋驅(qū)動流 [J]. 計(jì)算力學(xué)學(xué)報(bào), 2002, 19(3): 281-285. CHEN Xuejiang, QIN Guoliang, XU Zhong. Spectral element method and high order splitting method for Navies-Stokes equation [J]. Chinese Journal of Computational Mechanics, 2002, 19(3): 281-285.

[10]WILLIAM GEAR C. Numerical initial value problem in ordinary differential equations [M]. Englewood Cliffs, NJ, USA: Prentice-Hall Inc., 1973: 147-158.

[11]張家忠. 非線性動力系統(tǒng)的運(yùn)動穩(wěn)定性、分岔理論及其應(yīng)用 [M]. 西安: 西安交通大學(xué)出版社, 2010: 12-19.

[12]ERCILIA S. The controversial stability analysis [J]. Applied Mathematics and Computation, 2003, 145(2/3): 777-794.

[13]VAN DORSSELAER J L M, HUNDSDORFER W. Stability estimates based on numerical ranges with an application to a spectral method [J]. BIT Numerical Mathematics, 1994, 34(2): 228-238.

[本刊相關(guān)文獻(xiàn)鏈接]

耿艷輝,秦國良.Chebyshev譜元法求解含吸收邊界的二維均勻穩(wěn)定流場的聲傳播.2012,46(3):100-106.[doi:10.7652/xjtuxb201203018]

詹飛龍,張楚華.對流擴(kuò)散方程的擬譜格式穩(wěn)定性分析及擴(kuò)穩(wěn)方法.2012,46(3):113-118.[doi:10.7652/xjtux2012030 20]

孫旭,張家忠,黃科峰.基于彈簧近似的非結(jié)構(gòu)化網(wǎng)格自適應(yīng)處理方法.2010,44(9):104-108.[doi:10.7652/xjtuxb2010 09021]

章爭榮,張湘?zhèn)?對流擴(kuò)散方程的數(shù)值流形格式及其穩(wěn)定性分析.2010,44(9):104-108.[doi:10.7652/xjtuxb201009021]

張榮欣,秦國良.用切比雪夫譜元方法求解均勻流場中的聲傳播問題.2009,43(7):120-124.[doi:10.7652/xjtuxb2009 07026]

朱昌允,秦國良,徐忠.譜元方法求解波動方程時(shí)顯式與隱式差分方法的比較.2008,42(9):1142-1145.[doi:10.7652/xjtuxb200809017]

朱昌允,秦國良,徐忠.譜元方法求解波動方程及影響其數(shù)值精度的相關(guān)因素.2008,42(1):56-59.[doi:10.7652/xjtuxb 200801013]

許靖,秦國良,朱昌允.譜元方法求解含吸收邊界的聲學(xué)問題.2007,41(7):875-878.[doi:10.7652/xjtuxb200707027]

(編輯 苗凌)

Spectral Element Method for Convection-Diffusion Equation with Stability Analysis

HE Wenqiang,QIN Guoliang

(Institute of Fluid Machinery, Xi’an Jiaotong University, Xi’an 710049, China)

The Chebyshev spectral element method combining with the semi-implicit Adams method is presented for solving the one-dimensional convection-diffusion equation, and the feasibility is verified by a numerical example. The stability condition of different discrete forms of spectral element method is deduced with character analysis, and the influences of time step, grid partitioning, interpolation order of convective and viscous terms are discussed. It is demonstrated that numerical solution with high accuracy can be gained with the coupled spectral element and semi-implicit Adams method for the convection-diffusion equation. Larger stability domain and higher accuracy can be achieved without extra-calculation as first-order viscous terms and second-order convective terms are interpolated in Adams method.

convection-diffusion equation; spectral element method; stability; semi-implicit Adams method

2014-04-21。

和文強(qiáng)(1982—),男,博士生;秦國良(通信作者),男,教授,博士生導(dǎo)師。

國家重點(diǎn)基礎(chǔ)研究發(fā)展規(guī)劃資助項(xiàng)目(2012CB026004)。

時(shí)間:2014-10-31

10.7652/xjtuxb201501001

O357.1

A

0253-987X(2015)01-0006-01

網(wǎng)絡(luò)出版地址:http:∥www.cnki.net/kcms/detail/61.1069.T.20141031.1643.019.html