時永興

(華電江蘇能源有限公司句容發(fā)電廠，江蘇 鎮(zhèn)江 212400)

譜消去粘性譜元方法求解對流擴散方程

時永興

(華電江蘇能源有限公司句容發(fā)電廠，江蘇 鎮(zhèn)江 212400)

針對譜元方法求解高雷諾數(shù)下對流擴散方程的穩(wěn)定性問題，采用Chebyshev譜元方法結(jié)合譜消去粘性法求解一維對流擴散方程。利用特征分析法預測了數(shù)值方法的求解穩(wěn)定性，通過數(shù)值算例驗證了該解法的可行性，討論了譜消去粘性參數(shù)對求解穩(wěn)定性及數(shù)值精度的影響。結(jié)果表明：和譜元方法相比，譜消去粘性譜元方法求解對流擴散方程的穩(wěn)定區(qū)域有了明顯的擴大，在高雷諾數(shù)時能夠獲得具有較高精度的數(shù)值解；較大的譜消去粘性項有利于穩(wěn)定區(qū)域的擴大，而在計算穩(wěn)定的條件下較小的粘性項有利于數(shù)值精度的提高，所以適當?shù)卦O置粘性項的大小，在保證計算穩(wěn)定的同時提高數(shù)值精度。

對流擴散方程；譜元法；譜消去粘性法；穩(wěn)定性；高雷諾數(shù)

0 引言

對流擴散方程一般用來描述濃度、溫度、渦量等物理量的轉(zhuǎn)移擴散過程，在環(huán)境、化工、氣象、水利等多個工程領域中都有著廣泛的應用，因其重要性而備受關注。由于一般情況下難以求得解析解，高精度的數(shù)值解對于準確地理解物理過程中的對流擴散現(xiàn)象就有著重要的意義。

譜元方法[1]SEM (spectral element method)作為一種靈活的高精度數(shù)值方法，結(jié)合了譜方法的高精度和有限元法的復雜區(qū)域處理能力，非常適宜應用于對流擴散方程的數(shù)值求解。但當方程中含有很高的雷諾數(shù)Re時，應用譜元方法進行求解會出現(xiàn)傳統(tǒng)的有限差分法[2]、有限元法[3]、譜方法[4]所遇到的穩(wěn)定性問題，即對于一定的時間步長和網(wǎng)格劃分存在Re數(shù)區(qū)間，只有當計算位于該Re數(shù)區(qū)間時，求解才是穩(wěn)定的。因此，擴大穩(wěn)定區(qū)域的范圍成為譜元方法求解對流擴散方程的關鍵問題。

譜消去粘性法SVV(spectral vanishing viscosity)是由Tadmar[5]在利用Fourier譜方法求解雙曲型守恒方程時首先提出的，Maday[6]等將其擴展到求解非周期問題的Legendre譜方法獲得了守恒方程穩(wěn)定的數(shù)值解；Heinrichs[7]利用譜消去粘性譜方法求解了含有邊界層的穩(wěn)態(tài)對流擴散方程，消除了邊界層所造成的數(shù)值振蕩；Karamanos和Karniadakis[8]首次將譜消去粘性法和譜元方法相結(jié)合，成功地計算了管道湍流；Xu和Pasquetti[9]將譜消去粘性法應用于配置點譜元方法，推廣了其求解高維復雜區(qū)域問題的形式，并獲得了高Re數(shù)下圓柱繞流的穩(wěn)定數(shù)值解；Kirby和Sherwin[10]利用矩陣形式分析了譜消去粘性法和譜元方法的結(jié)合方式，對三維三角形管道湍流進行了求解。

本文在前期譜消去粘性譜元法研究的基礎上，采用矩陣分析法提出求解過程的穩(wěn)定性條件，系統(tǒng)地研究譜消去粘性譜元法求解對流擴散方程的穩(wěn)定性問題及其擴大穩(wěn)定域的效果。討論譜消去粘性參數(shù)對穩(wěn)定性和計算精度的影響，通過數(shù)值算例對譜消去粘性譜元法求解高Re數(shù)下對流擴散方程的有效性進行驗證。

1 譜消去粘性法

譜消去粘性法求解對流擴散方程時，在方程中加入人工粘性項，穩(wěn)定化的一維對流擴散方程形式為：

00 ，

(1)

相應的初始條件和邊界條件：

u(x,0)=u0(x),0

(2)

u(0,t)=g1(t),u(1,t)=g2(t)，t>0

(3)

ε，Q依賴于用式(1)求解的多項式的階數(shù)N，ε通常取ε≈N-1，Q在一維標準區(qū)域(-1，1)內(nèi)的定義為：

(4)

2 Chebyshev譜元離散

采用譜元方法求解式(1)～式(3)即為尋找u∈H1，使得在時間區(qū)間[0,1]上滿足：

ε(Qu,v)=(f,v),?v∈H1，

(5)

在邊界上滿足u(0,t)=g1(t)，u(1,t)=g2(t)，其中(·,·)表示定義在[0,1]上的L2內(nèi)積。

由式(5)可得，譜消去粘性項的弱形式為：

Svv(u,v)=ε(Qu,v) ，

(6)

為了保持其系數(shù)矩陣的對稱性，可改寫為：

Svv(u,v)=ε(Q1/2u,Q1/2v) 。

(7)

由Chebyshev多項式的性質(zhì)可知，式(6)、式(7)所表示的譜消去粘性項的弱形式并不完全相同，但式(7)保證了系數(shù)矩陣的對稱性，可在計算中近似地采用。

實際計算時，將式(5)中的粘性項和譜消去粘性項相合并，則可以得到：

(f,v),?v∈H1，

(8)

式中：S=I+εReQ，I為單位矩陣。

(9)

(10)

在標準單元Λ內(nèi)對式(8)進行離散，則可以得到第i個單元的譜元離散形式：

(11)

(12)

初始條件為：u(x,0)=u0(x)。

3 時間離散

在時域上求解微分方程組一般有顯式和隱式兩種積分法。隱式積分法絕對穩(wěn)定，但每一個時間步都需要求解線性方程組，顯式積分法計算高效，需要的計算機內(nèi)存小，但一般條件穩(wěn)定。本文采用半隱式中心差分格式對時間項進行離散，對流項采用顯式處理，粘性項采用隱式處理，融合了顯式方法計算簡單和隱式方法穩(wěn)定性好的特點，式(12)的全離散形式為：

c1un+1=c2un+c3un-1+Bfn，

(13)

4 穩(wěn)定性分析

4.1譜消去粘性譜元方法穩(wěn)定性條件

將式(13)兩側(cè)左乘c1的逆矩陣得：

(14)

由于僅考慮第1類邊界條件，即在邊值的計算中不引入誤差，所有的誤差來自初始值的迭代過程。假定在初值的計算中引入了擾動τ0，則擾動τ滿足：

τn+1=M1τn+M2τn-1，

(15)

Dorsselaer等[11]分析了形如式(15)的全離散格式的穩(wěn)定性條件，定義擾動Vn+1=((τn+1)T，(τn)T)T，其滿足：

Vn+1=EVn，n≥1 ，

(16)

可以驗證當采用Chebyshev函數(shù)為單元基函數(shù)時，擾動系數(shù)矩陣E為正規(guī)矩陣。譜消去粘性譜元法求解穩(wěn)定的充分條件為：矩陣E的譜半徑ρ(E)≤1。

時永興：譜消去粘性譜元方法求解對流擴散方程

4.2穩(wěn)定性結(jié)果

圖1 不同時間步長的穩(wěn)定性曲線

將粘性振幅設置為ε=α/Nx，通過改變自由參數(shù)α的值，來研究ε的大小對求解穩(wěn)定性的影響。圖2給出了不同粘性振幅譜消去粘性譜元方法的求解穩(wěn)定域?？梢钥闯?，對于較大的時間步長，ε值的增加對求解穩(wěn)定域的擴大作用有限；對于較小的時間步長，求解穩(wěn)定域隨著ε值的增加不斷擴大。說明對于較小的時間步長，較大的譜消去粘性項有利于求解穩(wěn)定域的擴大。圖2中截斷階數(shù)的設置為mN=Nx/2。

圖2 不同粘性振幅的穩(wěn)定性曲線

圖3給出了不同截斷階數(shù)譜消去粘性譜元方法的求解穩(wěn)定域?？梢钥闯觯瑢τ?種時間步長，當截斷階數(shù)mN從Nx減小到2時，求解穩(wěn)定域基本上保持不變；當截斷階數(shù)mN取0或1時，求解穩(wěn)定域顯著增加，同樣說明較大的譜消去粘性項有利于求解穩(wěn)定域的擴大。圖3中粘性振幅的設置為ε=5.0/Nx。

圖3 不同截斷階數(shù)的穩(wěn)定性曲線

5 算例及分析

當f=0時，利用變量分離法可得式(1)的一個解析解為：

(17)

相應的初始條件、邊界條件：

u(x,0)=0， 0≤x≤1 ，

(18)

u(0，t)=0，u(1,t)=1 ，t>0 ，

(19)

該算例在邊界x=1處有一個流動邊界層的存在，當Re數(shù)增大時，邊界層的厚度逐漸變小，速度變化劇烈，利用譜元方法進行計算時會產(chǎn)生數(shù)值振蕩。

本文采用ea和L2誤差來考察n時間層上數(shù)值解的準確性，兩種誤差的定義為：

(20)

(21)

通過數(shù)值解和解析解的比較來驗證譜消去粘性譜元方法求解對流擴散方程的有效性。

5.1譜元方法和譜消去粘性譜元方法的精度比較

選取網(wǎng)格劃分Nm=5，Nx=7，時間步長Δt=0.01，計算時間t=1.0，粘性參數(shù)ε=0.012/Nx，mN=Nx/2。由穩(wěn)定性條件可以計算出對于所選定的時間步長和網(wǎng)格劃分譜元方法求解的臨界Re=557。譜元方法和譜消去粘性譜元方法在Re=570求解時，不同節(jié)點的絕對誤差如圖4所示?？梢钥闯?，在靠近邊界層的區(qū)域譜元方法的數(shù)值結(jié)果已經(jīng)開始發(fā)散，而譜消去粘性譜元方法的數(shù)值結(jié)果仍然和解析解相吻合。

圖4 譜元方法和譜消去粘性譜元方法的絕對誤差

譜元方法和譜消去粘性譜元方法在不同Re數(shù)的L2誤差如圖5所示?？梢钥闯?，譜元方法的L2誤差隨著Re數(shù)的增加逐步增大，特別是當Re大于臨界值時數(shù)值結(jié)果快速發(fā)散；譜消去粘性譜元方法的L2誤差隨Re數(shù)的增加基本上保持不變，即使在穩(wěn)定區(qū)域較小的譜消去粘性項的增加仍有利于計算精度的提高。

圖5 譜元方法和譜消去粘性譜元方法的L2誤差

5.2粘性參數(shù)對譜消去粘性譜元方法的影響

圖6給出了譜消去粘性譜元方法不同ε值的L2誤差。選取2種不同的網(wǎng)格劃分Nm=4、Nx=4，Nm=8、Nx=8，研究粘性參數(shù)的變化對兩種網(wǎng)格計算精度的影響。網(wǎng)格劃分Nm=4，Nx=4的時間步長為Δt=0.05，相應的譜元求解臨界Re=100，計算Re=150；網(wǎng)格劃分Nm=8，Nx=8的時間步長為Δt=0.005，相應的譜元求解臨界Re=874，計算Re=1 000。計算時間t=1.0，截斷階數(shù)mN=Nx/2。可以看出，在保證計算穩(wěn)定的條件下，較小的粘性振幅有利于計算精度的提高，同時對于稀疏的網(wǎng)格，保證其計算穩(wěn)定所需要的譜消去粘性項較大，即粘性振幅較大。

圖6 譜消去粘性譜元方法不同ε的L2誤差

圖7給出了譜消去粘性譜元方法不同mN值的L2誤差。選取網(wǎng)格劃分Nm=4、Nx=6，Nm=6、Nx=6，研究截斷階數(shù)的變化對兩種網(wǎng)格計算精度的影響。網(wǎng)格劃分Nm=4、Nx=6的時間步長為Δt=0.02，相應的譜元求解臨界Re=270，計算Re=350；網(wǎng)格劃分Nm=6、Nx=6的時間步長為Δt=0.01，相應的譜元求解臨界Re=500，計算Re=600。計算時間t=1.0，粘性振幅ε=0.1/Nx?？梢钥闯?，截斷階數(shù)的變化對于計算精度的影響不是十分明顯，稀疏的網(wǎng)格需要較小的截斷階數(shù)以保證計算的穩(wěn)定。對于Nm=4、Nx=6的網(wǎng)格其計算在mN=Nx-1時已經(jīng)發(fā)散，而對于Nm=6、Nx=6的網(wǎng)格，其所需的譜消去粘性項較小，在mN=Nx-1時計算誤差反而達到最小。

圖7 譜消去粘性譜元方法不同mN的L2誤差

6 結(jié)論

本文采用Chebyshev譜元方法耦合譜消去粘性法求解了一維對流擴散方程，對該方法的求解穩(wěn)定性和計算精度進行了分析，討論了譜消去粘性參數(shù)對穩(wěn)定性及計算精度的影響，并與譜元方法的求解結(jié)果進行了對比，可以得出如下結(jié)論。

(1)譜元方法在求解含有高Re數(shù)的對流擴散方程時穩(wěn)定區(qū)域十分有限，譜粘消去粘性項的增加有效地擴大了穩(wěn)定區(qū)域的范圍。

(2)在高Re數(shù)時譜消去粘性項保證了譜元求解的數(shù)值精度，即使在穩(wěn)定區(qū)域較小的譜消去粘性項的增加仍然有利于數(shù)值精度的提高。

(3)粘性參數(shù)的設置對于譜消去粘性譜元方法的穩(wěn)定區(qū)域的擴大和計算精度的提高都有著非常大

時永興：譜消去粘性譜元方法求解對流擴散方程

的影響。較大的譜消去粘性項有利于穩(wěn)定域的擴大，但在保證計算穩(wěn)定的條件下，較小的粘性項又有利于精度的提高。

今后的工作可以從以下兩個方面展開：應用譜元方法求解非線性對流擴散方程，并對求解過程的穩(wěn)定性和影響因素進行分析；將譜消去粘性譜元方法拓展至三維非定常流動。

[1]PETERA A T.A spectral element method for fluid dynamics:laminar flow in a channel expansion[J].J comp phys,1984,54(3):468-488.

[2]陶文銓.數(shù)值傳熱學[M].2版.西安：西安交通大學出版社，2001:48-186.

[3]MOHAMED S A,MOHSMED N A,ABDEL GAWAD A F,et al.A modified diffusion coefficient technique for the convection diffusion equation[J].Appl math comput，2013,219:9317-9330.

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

[5]TADMOR E.Convergence of spectral methods for nonlinear conservation laws[J].SIAM J numer anal,1989,26(1):30-44.

[6]MADAY Y,OULD KABER S M,TADMOR E.Legendre pseudospectral viscosity method for nonlinear conservation laws[J].SIAM J numer anal,1993,30(2):321-342.

[7]HEINRICHS W.Spectral viscosity for convection dominated flow[J].J sci comput,1994,9(2):137-148.

[8]KARAMANOS G S,KARNIADAKIS G E.A spectral vanishing viscosity method for large-eddy simulation[J].J comput phys,2000,163:22-50.

[9]XU C J,PASQUETTI R.Stabilized spectral element computations of high Reynolds number incompressible flows[J].J comput phys,2004,196:680-704.

[10]KIRBY R M,SHERWIN S J.Stabilization of spectral/hp element methods through spectral vanishing viscosity:Application to fluid mechanics modeling[J].Comput methods appl mech.engrg,2006,195:3128-3144.

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

(本文責編：白銀雷)

O 357.1

A

1674-1951(2017)10-0025-05

2017-08-07；

2017-09-13

時永興(1972—)，男，江蘇常州人，工程師，從事火力發(fā)電廠生產(chǎn)管理工作(E-mail:syx_jrhd@126.com)。