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

時永興

(華電江蘇能源有限公司句容發電廠，江蘇 鎮江 212400)

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

時永興

(華電江蘇能源有限公司句容發電廠，江蘇 鎮江 212400)

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

對流擴散方程；譜元法；譜消去粘性法；穩定性；高雷諾數

0 引言

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

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

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

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

1 譜消去粘性法

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

00 ，

(1)

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

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

(2)

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

(3)

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

(4)

2 Chebyshev譜元離散

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

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

(5)

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

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

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

(6)

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

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

(7)

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

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

(f,v),?v∈H1，

(8)

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

(9)

(10)

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

(11)

(12)

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

3 時間離散

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

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

(13)

4 穩定性分析

4.1譜消去粘性譜元方法穩定性條件

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

(14)

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

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

(15)

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

Vn+1=EVn，n≥1 ，

(16)

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

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

4.2穩定性結果

圖1 不同時間步長的穩定性曲線

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

圖2 不同粘性振幅的穩定性曲線

圖3給出了不同截斷階數譜消去粘性譜元方法的求解穩定域。可以看出，對于2種時間步長，當截斷階數mN從Nx減小到2時，求解穩定域基本上保持不變；當截斷階數mN取0或1時，求解穩定域顯著增加，同樣說明較大的譜消去粘性項有利于求解穩定域的擴大。圖3中粘性振幅的設置為ε=5.0/Nx。

圖3 不同截斷階數的穩定性曲線

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數增大時，邊界層的厚度逐漸變小，速度變化劇烈，利用譜元方法進行計算時會產生數值振蕩。

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

(20)

(21)

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

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

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

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

譜元方法和譜消去粘性譜元方法在不同Re數的L2誤差如圖5所示。可以看出，譜元方法的L2誤差隨著Re數的增加逐步增大，特別是當Re大于臨界值時數值結果快速發散；譜消去粘性譜元方法的L2誤差隨Re數的增加基本上保持不變，即使在穩定區域較小的譜消去粘性項的增加仍有利于計算精度的提高。

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

5.2粘性參數對譜消去粘性譜元方法的影響

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

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

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

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

6 結論

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

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

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

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

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

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

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

[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]陶文銓.數值傳熱學[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—)，男，江蘇常州人，工程師，從事火力發電廠生產管理工作(E-mail:syx_jrhd@126.com)。