譜消去粘性譜元方法求解對(duì)流擴(kuò)散方程

時(shí)永興

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

譜消去粘性譜元方法求解對(duì)流擴(kuò)散方程

時(shí)永興

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

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

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

0 引言

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

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

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

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

1 譜消去粘性法

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

00 ，

(1)

相應(yīng)的初始條件和邊界條件：

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

(2)

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

(3)

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

(4)

2 Chebyshev譜元離散

采用譜元方法求解式(1)～式(3)即為尋找u∈H1，使得在時(shí)間區(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)可得，譜消去粘性項(xiàng)的弱形式為：

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

(6)

為了保持其系數(shù)矩陣的對(duì)稱(chēng)性，可改寫(xiě)為：

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

(7)

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

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

(f,v),?v∈H1，

(8)

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

(9)

(10)

在標(biāo)準(zhǔn)單元Λ內(nèi)對(duì)式(8)進(jìn)行離散，則可以得到第i個(gè)單元的譜元離散形式：

(11)

(12)

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

3 時(shí)間離散

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

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

(13)

4 穩(wěn)定性分析

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

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

(14)

由于僅考慮第1類(lèi)邊界條件，即在邊值的計(jì)算中不引入誤差，所有的誤差來(lái)自初始值的迭代過(guò)程。假定在初值的計(jì)算中引入了擾動(dòng)τ0，則擾動(dòng)τ滿足：

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

(15)

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

Vn+1=EVn，n≥1 ，

(16)

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

時(shí)永興：譜消去粘性譜元方法求解對(duì)流擴(kuò)散方程

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

圖1 不同時(shí)間步長(zhǎng)的穩(wěn)定性曲線

將粘性振幅設(shè)置為ε=α/Nx，通過(guò)改變自由參數(shù)α的值，來(lái)研究ε的大小對(duì)求解穩(wěn)定性的影響。圖2給出了不同粘性振幅譜消去粘性譜元方法的求解穩(wěn)定域。可以看出，對(duì)于較大的時(shí)間步長(zhǎng)，ε值的增加對(duì)求解穩(wěn)定域的擴(kuò)大作用有限；對(duì)于較小的時(shí)間步長(zhǎng)，求解穩(wěn)定域隨著ε值的增加不斷擴(kuò)大。說(shuō)明對(duì)于較小的時(shí)間步長(zhǎng)，較大的譜消去粘性項(xiàng)有利于求解穩(wěn)定域的擴(kuò)大。圖2中截?cái)嚯A數(shù)的設(shè)置為mN=Nx/2。

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

圖3給出了不同截?cái)嚯A數(shù)譜消去粘性譜元方法的求解穩(wěn)定域?？梢钥闯?，對(duì)于2種時(shí)間步長(zhǎng)，當(dāng)截?cái)嚯A數(shù)mN從Nx減小到2時(shí)，求解穩(wěn)定域基本上保持不變；當(dāng)截?cái)嚯A數(shù)mN取0或1時(shí)，求解穩(wěn)定域顯著增加，同樣說(shuō)明較大的譜消去粘性項(xiàng)有利于求解穩(wěn)定域的擴(kuò)大。圖3中粘性振幅的設(shè)置為ε=5.0/Nx。

圖3 不同截?cái)嚯A數(shù)的穩(wěn)定性曲線

5 算例及分析

當(dāng)f=0時(shí)，利用變量分離法可得式(1)的一個(gè)解析解為：

(17)

相應(yīng)的初始條件、邊界條件：

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

(18)

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

(19)

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

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

(20)

(21)

通過(guò)數(shù)值解和解析解的比較來(lái)驗(yàn)證譜消去粘性譜元方法求解對(duì)流擴(kuò)散方程的有效性。

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

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

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

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

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

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

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

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

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

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

6 結(jié)論

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

(1)譜元方法在求解含有高Re數(shù)的對(duì)流擴(kuò)散方程時(shí)穩(wěn)定區(qū)域十分有限，譜粘消去粘性項(xiàng)的增加有效地?cái)U(kuò)大了穩(wěn)定區(qū)域的范圍。

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

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

時(shí)永興：譜消去粘性譜元方法求解對(duì)流擴(kuò)散方程

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

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

[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ù)值傳熱學(xué)[M].2版.西安：西安交通大學(xué)出版社，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.

(本文責(zé)編：白銀雷)

O 357.1

A

1674-1951(2017)10-0025-05

2017-08-07；

2017-09-13

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