谱元方法求解对流扩散方程及其稳定性分析

和文强,秦国良

(西安交通大学流体机械研究所,710049,西安)



谱元方法求解对流扩散方程及其稳定性分析

和文强,秦国良

(西安交通大学流体机械研究所,710049,西安)

探讨了一维对流扩散方程的一种高精度数值解法,该解法在空间上采用了Chebyshev谱元方法,在时间上结合了半隐式Adams方法。通过数值算例验证了解法的可行性,利用特征分析法得到了对流扩散方程谱元求解时不同离散形式的稳定性条件,并对数值求解的稳定性进行了预测。通过时间步长、网格划分、对流项和黏性项插值阶数的影响研究表明:耦合Chebyshev谱元方法和半隐式Adams方法在求解对流扩散方程时能够获得高精度的数值解;时间离散时Adams方法的黏性项采用一阶插值形式、对流项采用二阶插值形式,在未增加计算量的同时能够获得较大的稳定区域和较高的计算精度。

对流扩散方程;谱元法;稳定性;半隐式Adams方法

对流扩散方程是描述流体流动和传热、污染物输移和扩散、电化学反应等的典型方程之一,利用数值方法进行求解时,由于对流项的存在,求解过程可能会出现不稳定现象,因此数值方法的稳定性问题一直是对流扩散方程研究的热点。当对流扩散方程的对流项占优时,应用传统的有限差分法和有限元法求解会出现数值振荡[1-2],Gottlieb等在利用拟谱方法求解时也发现了不稳定现象[3],而Mofid等通过引入罚方法来处理边界条件能够部分或完全避免这种不稳定现象[4]。

谱元方法[5]结合了谱方法的高精度和有限元法的灵活性特点,使得对流扩散问题求解有望获得高精度的数值解。目前,谱元方法应用的主要障碍是高雷诺数Re下的稳定性问题。因为谱元方法在单元上需对求解变量进行谱近似,当Re增大、黏性耗散作用减小时,任意微小扰动都会引起单元数值解不稳定且最终扩展到整个求解区域。为了扩大Re范围,一些稳定性措施[6-7]引入到了谱元方法之中。

本文在前期谱元方法[8]研究的基础上,采用特征分析法系统地分析了Chebyshev谱元方法的半离散形式和Chebyshev谱元方法耦合半隐式Adams时间离散方法的全离散形式的稳定性条件,讨论了半隐式Adams时间离散方法中对流项和黏性项在不同插值阶数、网格划分、时间步长时对稳定性的影响,探索了最优离散格式,试图扩大稳定区域,提高计算精度和计算效率,最后结合具体算例对Chebyshev谱元方法耦合半隐式Adams方法在求解对流扩散方程时的有效性进行了验证。

1 对流扩散方程

一维对流扩散方程为

(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)

插值函数为

(7)

式中:Tm、Tn为Chebyshev多项式;Nx为单元网格数;cm满足

(8)

在标准单元内对式(4)进行离散得

(9)

(10)

式(10)为u关于时间的一阶线性常微分方程组。

3 时间离散

式(10)在时间上的逐步积分途径一般有显式和隐式两种。显式积分计算效率较高,需要的计算机内存相对较小,但计算精度较低,条件稳定;隐式积分的精度高,稳定性较好,但每一个积分步都要求解一个线性方程组,不便于直接计算。本文时间项采用半隐式Adams方法进行离散,对流项采用Je阶显式Adams-Bashforth形式,黏性项采用Ji阶隐式Adams-Monlton形式[9],这样便融合了显式方法耦合简单和隐式方法稳定性好的特点。式(10)的全离散形式为

(11)

式中:βj为Adams-Bashforth系数[10];γk为Adams-Monlton系数[10]。

4 谱元方法的稳定性分析

4.1 半离散方程的稳定性

将式(10)两侧左乘B的逆矩阵可得

(12)

由于仅考虑第一类边界条件,所以在边值计算中不引入误差,而假定在初值的计算中引入了扰动ε0,扰动ε满足

(13)

式中:ε=(ε1,ε2,…,εn-1);L为D的内配置点n-1阶子矩阵。

根据Lyapunov稳定性理论[11],式(12)稳定的充要条件为矩阵L的一切特征根都有非正实部,且每个有0实部的特征根对应矩阵L的约当标准形中的一维约当块。也就是说,如果存在一个有正实部的特征根,式(12)对于任何时间离散都是不稳定的。矩阵L的特征值实部最大值与Re、单元数、单元网格数的关系如图1所示。

图1 特征值实部最大值与单元数、单元网格数的关系

对于不同的Re和单元网格数,图1中2种单元划分的空间谱元离散始终满足半离散方程稳定性要求,其他单元划分的情形也可以进行验证。

4.2 全离散方程的稳定性

4.2.1 单步法全离散的稳定性条件 采用半隐式Adams法时间离散,对流项、黏性项的插值阶数均为1,即n+1时刻的变量仅与n时刻的变量相关时,式(11)可表示为

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

(14)

式中:φ、φf为时间离散函数。

扰动ε的全离散形式为

εn+1=Lεn

(15)

Sousa分析了形如式(14)的全离散格式的稳定性条件[12],即:当L为正规矩阵时,其谱半径ρ(L)≤1为稳定的充分条件;当L为非正规矩阵时,其谱半径ρ(L)≤1仅为稳定的必要条件。可以验证,当采用的Chebyshev函数为单元基函数时,L为正规矩阵,所以单步谱元法稳定的充分条件为ρ(L)≤1。

4.2.2 多步法全离散的稳定性条件 采用半隐式Adams法时间离散,对流项或黏性项的插值阶数大于1,即n+1时刻的变量与前n-k+1时刻的变量相关时,式(11)可表示为

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

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

(16)

式中:φj、φf为时间离散函数;Dn-k+j为un-k+j的系数矩阵。

扰动ε的全离散形式为

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

(17)

式中:Lj为φj(ΔtDn-k+j)的内配置点n-1阶子矩阵。

Dorsselaer等分析了形如式(16)的全离散格式的稳定性条件[13],定义扰动

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

满足

Vn+1=EVn,n≥k

(18)

式中:I为单位矩阵;0为零矩阵。

多步谱元法稳定的充分条件为ρ(E)≤1。

4.3 不同形式全离散方程的稳定区域

通过选用4种不同的网格划分方式,分别表示为3×5(Nm=3,Nx=5)、5×7(Nm=5,Nx=7)、10×6(Nm=10,Nx=6)、10×10(Nm=10,Nx=10),时间(计算时间与特征时间的比值)步长Δt∈[0.001,0.1],来比较不同网格划分、时间步长和时间离散方式的求解稳定性。

对于半隐式Adams法时间离散,当对流项采用一阶插值形式,黏性项分别采用一阶、二阶插值形式时,式(10)的全离散形式可分别表示为

(19)

(20)

将式(19)、式(20)分别表示为式(14)、式(16)的形式,即可由单步法和多步法的稳定性条件对求解稳定性进行判定。对流项一阶插值,黏性项一阶、二阶插值时的求解稳定域如图2所示。从图2可以看出,黏性项二阶插值的稳定域和一阶插值基本相同,网格划分对稳定域的影响很小,上临界Re随着时间步长的减小而逐步增大。

图2 黏性项不同插值阶数的稳定性曲线

对于半隐式Adams法时间离散,对流项采用一阶插值形式,黏性项采用三阶插值形式时,式(10)的全离散形式可表示为

(21)

将式(21)表示为式(16)的形式,利用多步法稳定性条件对求解稳定性进行判定。对流项一阶插值、黏性项三阶插值时的求解稳定域如图3所示。从图3可以看出,黏性项三阶插值时,不仅存在上临界Re,也存在下临界Re,只有当Re位于上临界Re和下临界Re之间时求解才是稳定的。上临界Re随时间步长的增大而减小,下临界Re随时间步长的增大而增大,同时求解稳定域随网格节点数的增加而减小。黏性项三阶插值对求解的稳定性的影响很大,且显著减小了稳定域的范围。

图3 黏性项三阶插值的稳定性曲线

对于半隐式Adams法时间离散,当黏性项采用一阶插值形式,对流项分别采用二阶、三阶插值形式时,式(10)的全离散形式可表示为

(22)

(23)

将式(22)、式(23)表示为式(16)的形式,同样利用多步法稳定性条件对求解稳定性进行判定。黏性项一阶插值,对流项一阶、二阶、三阶插值时的求解稳定域如图4所示。从图4可以看出:对流项二阶插值存在临界时间步长,上临界Re在时间步长大于临界值时有所减小,在时间步长小于临界值时大幅度增加;对流项三阶插值也存在临界时间步长,上临界Re在时间步长大于临界值时比二阶插值进一步减小,在时间步长小于临界值时迅速增加到无穷大。网格节点数对临界时间步长的影响也很大,随节点数的增加,临界时间步长减小,对于10×10的网格划分,二阶、三阶插值的临界时间步长均小于0.001。

图4 对流项不同插值阶数的稳定性曲线

5 计算实例和数值结果

(24)

统一选取时间步长Δt=0.001,比较3单元5网格和10单元10网格2种网格划分,并给出t=1.0时的计算误差。

5.1 单步法全离散数值结果

单步法全离散的数值结果如图5所示。Re小于上临界Re时,2种网格划分都取得了正确的数值结果,但由于时间离散阶数较低,所以数值结果的精度不高。Re大于上临界Re时,10单元10网格的数值结果快速发散,3单元5网格的数值结果虽有发散的趋势,但速度非常平缓,说明节点数越多,不稳定作用越强。

图5 单步谱元法的计算误差

5.2 多步法全离散数值结果

对流项采用一阶插值,黏性项采用二阶、三阶插值时的数值结果如图6所示。对于黏性项二阶插值,当Re小于上临界Re时,计算误差和一阶插值相差很小;当Re大于上临界Re且节点数较多时,数值结果的发散速率远大于节点数较少时的发散速率。黏性项三阶插值的计算误差和二阶插值基本相同,当Re大于上临界Re或小于下临界Re且节点数较多时,数值结果快速发散。

图6 黏性项二阶、三阶插值的计算误差

黏性项采用一阶插值,对流项采用二阶、三阶插值时的数值结果如图7所示。对流项采用二阶插值时,计算误差相较一阶插值大大减小,同时Re较大时的误差下降更多,Re大于上临界Re且节点数较多时,数值结果快速发散。对流项三阶插值的计算精度与二阶插值基本相同。

图7 对流项二阶、三阶插值的计算误差

6 结 论

本文采用Chebyshev谱元方法耦合半隐式Adams时间离散方法求解了一维对流扩散方程,并对求解过程中方程的半离散和全离散形式的稳定性条件进行了讨论,分析了半隐式Adams时间离散对流项、黏性项在不同插值阶数、网格划分和时间步长时对求解稳定性及计算精度的影响,得出如下结论。

(1)对流扩散方程的空间谱元离散,对于任意的Re,网格划分都满足半离散方程的稳定性条件。

(2)随着对流项插值阶数的增加,求解稳定域得到了扩大,计算精度大大提高。对流项二阶插值、时间步长小于临界值时,上临界Re能取得较大的值;时间步长大于临界值时,上临界Re减小幅度减缓,适宜在实际数值计算中应用。

(3)随着黏性项插值阶数的增加,求解稳定域缩小。三阶插值时出现了下临界Re,稳定域大大减小。

今后的工作可以从以下两方面展开:在方程中通过增加谱消除人工黏性项,使得谱元方法能够求解的稳定区域扩大;探索高精度的时间处理方法,如时空耦合的最小二乘谱元方法。

[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]. 应用力学学报, 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] 陈雪江, 秦国良, 徐忠. 谱元法和高阶时间分裂法求解方腔顶盖驱动流 [J]. 计算力学学报, 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]张家忠. 非线性动力系统的运动稳定性、分岔理论及其应用 [M]. 西安: 西安交通大学出版社, 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.

[本刊相关文献链接]

耿艳辉,秦国良.Chebyshev谱元法求解含吸收边界的二维均匀稳定流场的声传播.2012,46(3):100-106.[doi:10.7652/xjtuxb201203018]

詹飞龙,张楚华.对流扩散方程的拟谱格式稳定性分析及扩稳方法.2012,46(3):113-118.[doi:10.7652/xjtux2012030 20]

孙旭,张家忠,黄科峰.基于弹簧近似的非结构化网格自适应处理方法.2010,44(9):104-108.[doi:10.7652/xjtuxb2010 09021]

章争荣,张湘伟.对流扩散方程的数值流形格式及其稳定性分析.2010,44(9):104-108.[doi:10.7652/xjtuxb201009021]

张荣欣,秦国良.用切比雪夫谱元方法求解均匀流场中的声传播问题.2009,43(7):120-124.[doi:10.7652/xjtuxb2009 07026]

朱昌允,秦国良,徐忠.谱元方法求解波动方程时显式与隐式差分方法的比较.2008,42(9):1142-1145.[doi:10.7652/xjtuxb200809017]

朱昌允,秦国良,徐忠.谱元方法求解波动方程及影响其数值精度的相关因素.2008,42(1):56-59.[doi:10.7652/xjtuxb 200801013]

许靖,秦国良,朱昌允.谱元方法求解含吸收边界的声学问题.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。

和文强(1982—),男,博士生;秦国良(通信作者),男,教授,博士生导师。

国家重点基础研究发展规划资助项目(2012CB026004)。

时间:2014-10-31

10.7652/xjtuxb201501001

O357.1

A

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

网络出版地址:http:∥www.cnki.net/kcms/detail/61.1069.T.20141031.1643.019.html