«上一篇
文章快速检索     高级检索
下一篇»
  哈尔滨工程大学学报  2018, Vol. 39 Issue (3): 446-451  DOI: 10.11990/jheu.201608033
0

引用本文  

和文强, 秦国良, 包振忠, 等. 稳定化谱元方法求解二维稳态对流扩散方程[J]. 哈尔滨工程大学学报, 2018, 39(3): 446-451. DOI: 10.11990/jheu.201608033.
HE Wenqiang, QIN Guoliang, BAO Zhenzhong, et al. Stabilized spectral element method for solving a two-dimensional steady convection-diffusion equation[J]. Journal of Harbin Engineering University, 2018, 39(3): 446-451. DOI: 10.11990/jheu.201608033.

基金项目

国家自然科学基金项目(51776155)

通信作者

秦国良, E-mail:glqin@mail.xjtu.edu.cn

作者简介

和文强(1982-), 男, 博士研究生; 秦国良(1964-), 男, 教授, 博士生导师

文章历史

收稿日期:2016-08-17
网络出版日期:2017-12-15
稳定化谱元方法求解二维稳态对流扩散方程
和文强, 秦国良, 包振忠, 王亚洲    
西安交通大学 流体机械研究所, 陕西 西安 710049
摘要:为了求解对流项占优的对流扩散方程的非稳定性问题,本文提出了二维稳态对流扩散方程的一种高精度求解方法。该方法将一致逼近迎风方法和Chebyshev谱元方法相结合,形成了稳定化谱元方法。通过数值算例对该方法的可行性进行了验证,讨论了计算精度、收敛速度及边界层逼近和单元插值阶数的关系。研究表明:和谱元方法相比,稳定化谱元方法扩大了对流扩散方程的稳定求解区域;消除了计算区域内部及边界层附近的数值振荡,对内部和外部边界层都进行了很好的逼近;单元插值阶数的增加使计算精度及边界层的逼近效果都获得了迅速的提高。
关键词对流扩散方程    谱元法    稳定性    高精度    单元插值阶数    数值振荡    边界层    
Stabilized spectral element method for solving a two-dimensional steady convection-diffusion equation
HE Wenqiang, QIN Guoliang, BAO Zhenzhong, WANG Yazhou    
Fluid Machinery Institute, Xi'an Jiaotong University, Xi'an 710049, China
Abstract: A numerical solution method with high accuracy is proposed herein for the instability problems of a two-dimensional steady convection-dominated convection-diffusion equation. The new scheme is a stabilized spectral element method based on the consistent approximate upwind method and the Chebyshev spectral element method. The feasibility of the method is verified by numerical examples. Moreover, the relation among the calculation accuracy, convergence rate, boundary layer approximation, and the order of the element interpolation is discussed. The results show that compared with the spectral element method, the stabilized spectral element method greatly enlarges the stable solution region for the convection-diffusion equation. The numerical oscillations in the computational domain and near the boundary layers in the convection-diffusion problem are completely eliminated. In addition, the internal and external boundary layers are excellently approximated. The computational accuracy and the effect of approximation to the boundary layers are both rapidly enhanced with the increase of the element interpolation order.
Key words: convection-diffusion equation    spectral element method    stability    high accuracy    element interpolation order    numerical oscillation    boundary layer    

对流扩散方程是最为常见的运动学方程,可以用来描述大气、海洋、河流中的污染物分布,地下石油的模拟开采,多孔介质流动,流体中的热量传递等众多物理过程,因其重要性而备受关注[1-4]。由于实际问题的复杂性,对流扩散方程的精确解通常难以求出,因此该方程的求解十分依赖有效的数值方法。常见的有限容积法、有限差分法、有限元法等,都是求解对流扩散方程的主要方法[5-7],其中有限元法因其计算时的高效率及良好的适应性在对流扩散方程的研究中获得了广泛的应用。但有限元法求解对流项占优的对流扩散方程时,会在计算区域内产生剧烈的数值振荡,导致数值不稳定[8]。针对这种不稳定现象,国内外学者提出了稳定化有限元法[9-13],其中SUPG(streamline upwind Petrov-Galerkin)方法因其简单、有效的特点在实际问题的计算中得到了成功的应用。SUPG方法在有限元变分形式中增加了沿流线方向的人工黏性,通过适当的迎风函数的定义,所增加的人工黏性在提高稳定性的同时不破坏数值解的精度。然而对流项占优的对流扩散过程往往包含有边界层的产生,对于这类问题SUPG方法基本上消除了计算区域内的不稳定现象,而在边界层邻域内却依然存留有数值振荡。CAU(consistent approximate upwind)方法以SUPG方法为基础,同时又添加了边界层梯度方向的人工黏性[14-15],从而有效地抑制了边界层邻域内的数值振荡。但有限元法所采用的插值基函数使得它的数值解仅以代数阶的速度收敛。

谱元方法(spectral element method,SEM)结合了谱方法的高精度和有限元法处理复杂区域灵活性的特点[16-17],其采用在空间正交且无穷阶光滑的多项式作为插值基函数,当插值阶数增加时,数值解的收敛速度是指数阶的。然而谱元方法在求解对流扩散方程时也会出现剧烈的数值振荡,因此引入适当的稳定性措施,扩大谱元方法求解对流项占优的对流扩散方程的稳定求解域,保持数值解的精度,对于谱元方法显得非常的重要。

本文将CAU方法和Chebyshev谱元方法相结合形成稳定化谱元方法(stabilized spectral element method,SSEM)[18-19],并用其求解对流项占优且含有边界层的二维稳态对流扩散方程,数值验证该方法在扩大稳定求解域、维持计算精度方面的实用性,讨论插值阶数对计算误差、收敛速度及边界层逼近效果的影响。

1 对流扩散方程

二维稳态对流扩散方程为

$ \gamma \phi + \mathit{\boldsymbol{u}} \cdot \nabla \phi - \varepsilon \Delta \phi = f,\;\;\;\;\mathit{\boldsymbol{x}} \in \mathit{\Omega } $ (1)

相应的边界条件:

$ \phi = {g_D},\;\;\;\;\mathit{\boldsymbol{x}} \in {\mathit{\Gamma }_D};\;\;\;\;\frac{{\partial \phi }}{{\partial \mathit{\boldsymbol{n}}}} = {g_N},\;\;\;\;\mathit{\boldsymbol{x}} \in {\mathit{\Gamma }_N} $

式中:Ω为计算区域, ∂Ω=ΓDΓN, ΓDΓN=ϕ, ϕ为求解的未知量; γ为反应系数; u为速度矢量;ε为扩散系数;f为源项;gDgN分别为第一类和第二类边界条件;n为边界外法线单位矢量。

定义试探函数和检验函数空间为

$ U = \left\{ {\phi \left( {\rm{x}} \right):\phi \in {H^1}\left( \mathit{\Omega } \right),\phi \left| {_{{\mathit{\Gamma }_D}}} \right. = {\mathit{g}_D}} \right\} $ (2)
$ V = \left\{ {\eta \left( x \right):\eta \in {H^1}\left( \mathit{\Omega } \right),\eta \left| {_{{\mathit{\Gamma }_D}}} \right. = 0} \right\} $ (3)

式中:Hs(Ω)为通常的Sobolev空间。式(1)的Galerkin变分问题为:求ϕU,使得

$ {A_G}\left( {\phi ,\eta } \right) = {F_G}\left( \eta \right),\forall \eta \in V $ (4)

式中:

$ {A_G}\left( {\phi ,\eta } \right) = \varepsilon \left( {\nabla \phi ,\nabla \eta } \right) + \left( {\mathit{\boldsymbol{u}},\nabla \phi ,\eta } \right) + \left( {\gamma \phi ,\eta } \right) $
$ {F_G}\left( \eta \right) = \left( {f,\eta } \right) + {\left( {{g_N},\eta } \right)_{{\mathit{\Gamma }_N}}} $
$ {\left( {{g_N},\eta } \right)_{{\mathit{\Gamma }_N}}} = \int_{{\mathit{\Gamma }_N}} {{g_N}\eta {\rm{d}}\mathit{\Gamma }} $
2 CAU逼近解

将计算区域划分为互相不重叠的单元, $ \overline{\mathit{\Omega} }=\cup _{e=1}^{{{N}_{e}}}{{\overline{\mathit{\Omega} }}_{e}} $,∩e=1NeΩe=ϕNe为单元数,UhVhUV的有限元子空间:

$ {U^h} = \left\{ {{\phi ^h} \in {C^0}\left( \mathit{\Omega } \right);{\phi ^h}\left| {_{{\mathit{\Omega }_e}}} \right. \in P_e^k,{\phi ^h}\left| {_{{\mathit{\Gamma }_D}}} \right. = {\mathit{g}_D}} \right\} $ (5)
$ {V^h} = \left\{ {{\eta ^h} \in {C^0}\left( \mathit{\Omega } \right);{\eta ^h}\left| {_{{\mathit{\Omega }_e}}} \right. \in P_e^k,{\eta ^h}\left| {_{{\mathit{\Gamma }_D}}} \right. = 0} \right\} $ (6)

式中Pek为定义在Ωe上的k阶多项式空间。

2.1 SUPG方法

式(1)在子空间的SUPG变分问题为:求ϕhUh,使得

$ \begin{array}{*{20}{c}} {{A_G}\left( {{\boldsymbol{\phi} ^h},{\eta ^h}} \right) + {A_{{\rm{SUPG}}}}\left( {{\boldsymbol{\phi} ^h},{\eta ^h}} \right) = }\\ {{F_G}\left( {{\eta ^h}} \right) + {F_{{\rm{SUPG}}}}\left( {{\eta ^h}} \right),\forall {\eta ^h} \in {V^h}} \end{array} $ (7)

式中:

$ {A_{{\rm{SUPG}}}}\left( {{\boldsymbol{\phi} ^h},{\eta ^h}} \right) = \sum\limits_{e = 1}^{{N_e}} {\left( {{L_e}\left( {{\boldsymbol{\phi} ^h}} \right),\tau _e^s{\mathit{\boldsymbol{u}}_e} \cdot \nabla {\eta ^h}} \right)} \left| {_{{\mathit{\Omega }_e}}} \right. $
$ {F_{{\rm{SUPG}}}}\left( {{\eta ^h}} \right) = \sum\limits_{e = 1}^{{N_e}} {\left( {{f^h},\tau _e^s{\mathit{\boldsymbol{u}}_e} \cdot \nabla {\eta ^h}} \right)} \left| {_{{\mathit{\Omega }_e}}} \right. $
$ {L_e}\left( {{\phi ^h}} \right) = {\gamma _e}{\phi ^h} + {\mathit{\boldsymbol{u}}_e} \cdot \nabla {\boldsymbol{\phi} ^h} - \varepsilon \Delta {\boldsymbol{\phi} ^h} $
$ \tau _e^s \sim \min \left\{ {\frac{{{l_e}}}{{p_e^k\left\| {{\mathit{\boldsymbol{u}}_e}} \right\|}},\frac{{l_e^2}}{{{{\left( {p_e^k} \right)}^4}\varepsilon }}} \right\} $

式中:τes为流线迎风函数[20]le为流线方向的特征长度,pek为插值阶数。由式(7)可知,SUPG方法在Galerkin变分的检验函数ηh上增加了扰动τesue·∇ηh,式(1)中对流项关于扰动的内积(ue·∇ϕh, τesue·∇ηh)|Ωe为SUPG方法提供了额外的稳定性,而迎风函数τes的定义方法使高阶项关于扰动的内积(-εΔϕh, τesue·∇ηh)|Ωe,对数值解的精度没有破坏。

2.2 CAU方法

式(1)在子空间的CAU变分问题为:求ϕhUh,使得

$ \begin{array}{*{20}{c}} {{A_G}\left( {{\phi ^h},{\eta ^h}} \right) = {A_{{\rm{SUPG}}}}\left( {{\phi ^h},{\eta ^h}} \right) + {A_{{\rm{CAU}}}}\left( {{\phi ^h},{\eta ^h}} \right) = }\\ {{F_G}\left( {{\eta ^h}} \right) + {F_{{\rm{SUPG}}}}\left( {{\eta ^h}} \right),\;\;\;\;\;\;\forall {\eta ^h} \in {V^h}} \end{array} $ (8)

式中:

$ {A_{{\rm{CAU}}}}\left( {{\boldsymbol{\phi} ^h},{\eta ^h}} \right) = \sum\limits_{e = 1}^{{N_e}} {\left( {{R_e}\left( {{\boldsymbol{\phi} ^h}} \right),\tau _e^c{\mathit{\boldsymbol{v}}_e} \cdot \nabla {\eta ^h}} \right)} \left| {_{{\mathit{\Omega }_e}}} \right. $
$ {R_e}\left( {{\boldsymbol{\phi} ^h}} \right) = {\gamma _e}{\phi^h} + {\mathit{\boldsymbol{u}}_e} \cdot \nabla {\boldsymbol{\phi} ^h} - \varepsilon \Delta {\boldsymbol{\phi} ^h} - {f^h} $
$ {\mathit{\boldsymbol{v}}_e} = \frac{{{R_e}\left( {{\boldsymbol{\phi} ^h}} \right)}}{{{{\left( {\left\| {\nabla {\boldsymbol{\phi} ^h}} \right\| + {k_T}} \right)}^2}}}\nabla {\phi ^h} $
$ \tau _e^c = \tau _e^s\max \left\{ {\frac{{\left\| {{\mathit{\boldsymbol{u}}_e}} \right\|}}{{\left\| {{\mathit{\boldsymbol{v}}_e}} \right\|}} - {\alpha _e},0} \right\} $
$ {\alpha _e} = \max \left( {1,\frac{{{\mathit{\boldsymbol{u}}_e} \cdot \nabla {\boldsymbol{\phi} ^h}}}{{{R_e}\left( {{\boldsymbol{\phi} ^h}} \right)}}} \right) $

式中:Re(ϕh)为式(1)的残量,ve为边界层方向的辅助速度,τec为一致逼近迎风函数[21]αe为补偿系数,kT为常量,计算中取1.0。由式(8)可知,CAU方法在SUPG方法的基础上增加了各向同性的非线性人工黏性对边界层附近的数值振荡进行捕捉,同时迎风函数τec的定义方法保证所增加的人工黏性对数值解的精度没有影响。

2.3 简单线性化迭代

由于CAU方法在边界层方向增加了非线性形式的人工黏性,对流扩散方程也相应的转化为非线性方程,式(8)也应该采用迭代的方法进行求解。利用第k个迭代步的计算结果ϕkh对CAU人工黏性项进行线性化近似:

$ c\left( {\boldsymbol{\phi} _k^h;\boldsymbol{\phi} _{k + 1}^h,{\eta ^h}} \right) = \left( {\tau _{e,k}^c{K_{e,k}}\nabla \boldsymbol{\phi} _{k + 1}^h \cdot \nabla {\eta ^h}} \right)\left| {_{{\mathit{\Omega }_e}}} \right. $ (9)

式中:

$ \tau _{e,k}^c = \tau _{e,k}^c\left( {{\mathit{\boldsymbol{v}}_e}\left( {\boldsymbol{\phi} _k^h} \right)} \right) $
$ {K_{e,k}} = \frac{{R_e^2\left( {\boldsymbol{\phi} _k^h} \right)}}{{{{\left( {\left\| {\nabla \boldsymbol{\phi} _k^h} \right\| + {k_T}} \right)}^2}}} $

由给定的ϕkh,求ϕk+1hUh, 使得

$ \begin{array}{*{20}{c}} {{A_G}\left( {\boldsymbol{\phi} _{k + 1}^h,{\eta ^k}} \right) + {A_{{\rm{SUPG}}}}\left( {\boldsymbol{\phi} _{k + 1}^h,{\eta ^k}} \right) + }\\ {\sum\limits_{e = 1}^{{N_e}} {c\left( {\boldsymbol{\phi} _k^h;\boldsymbol{\phi} _{k + 1}^h,{\eta ^h}} \right)\left| {_{{\mathit{\Omega }_e}}} \right.} = {F_G}\left( {{\eta ^h}} \right) + }\\ {{F_{{\rm{SUPG}}}}\left( {{\eta ^k}} \right),\;\;\;\;\;\forall {\eta ^h} \in {V^h}} \end{array} $

且‖ϕk+1h-ϕkh‖ < σtolσtol为迭代精度,ϕ0h为采用SUPG方法计算一次的结果。

3 谱元空间离散

在标准单元内由Chebyshev多项式的极值点构成插值基函数,则求解变量的单元逼近形式为

$ \phi \left( {\xi ,\eta } \right) = \sum\limits_{i = 0}^{{I^d}} {{\phi ^i}{N^i}\left( {\xi ,\eta } \right) = \boldsymbol{\phi} {\mathit{\boldsymbol{N}}_e}} $ (10)

$ {N^i}\left( {\xi ,\eta } \right) = \sum\limits_{j = 0}^{{I^x}} {\sum\limits_{k = 0}^{{I^y}} {{h^j}\left( \xi \right){h^k}\left( \eta \right)} } $
$ {h^j}\left( \xi \right) = \frac{2}{{{I^x}}}\sum\limits_{m = 0}^{I_e^x} {\frac{1}{{{c_j}{c_m}}}{T_m}\left( {{\xi ^j}} \right){T_m}\left( \xi \right)} $
$ {h^k}\left( \eta \right) = \frac{2}{{{I^y}}}\sum\limits_{n = 0}^{{I^y}} {\frac{1}{{{c_k}{c_n}}}{T_n}\left( {{\eta ^k}} \right){T_n}\left( \eta \right)} $

式中:Id为插值节点数,ϕi为求解变量在i节点的值,Nii节点的插值基函数,ϕ为求解变量的节点值列向量,N为节点插值基函数列向量,Tm=cos(m arccos x),IxIyxy方向的插值阶数。在单元外hj(ξ)、hk(ξ)为0,在单元内hj(ξp)=δjphk(ξq)=δkqcm满足:

$ {c_m} = \left\{ \begin{array}{l} 2,\;\;\;m = 0,{I^x}\\ 1,\;\;\;m \ne 0,{I^x} \end{array} \right. $

标准单元通过映射转化到e单元,则式(8)的单元矩阵方程为

$ \begin{array}{*{20}{c}} {\left( {{\mathit{\boldsymbol{B}}_e} + {\mathit{\boldsymbol{C}}_e} + {\mathit{\boldsymbol{D}}_e} + \mathit{\boldsymbol{B}}_e^1 + \mathit{\boldsymbol{C}}_e^1 + \mathit{\boldsymbol{D}}_e^1 + D_e^2} \right){\boldsymbol{\phi} _e} = }\\ {\left( {{\mathit{\boldsymbol{E}}_e} + \mathit{\boldsymbol{C}}_e^2} \right){\mathit{\boldsymbol{f}}_e} + \left( {{\mathit{\boldsymbol{S}}_e} + \mathit{\boldsymbol{S}}_e^1} \right){\mathit{\boldsymbol{g}}_{e,N}}} \end{array} $ (11)

其中,

$ {\mathit{\boldsymbol{B}}_e} = \left( {{\gamma _e};{\mathit{\boldsymbol{N}}_e},{\mathit{\boldsymbol{N}}_e}} \right) $
$ \mathit{\boldsymbol{B}}_e^1 = \mathit{\boldsymbol{\tau }}_e^s\left( {{\gamma _e};{\mathit{\boldsymbol{N}}_e},{\mathit{\boldsymbol{u}}_e} \cdot \nabla {\mathit{\boldsymbol{N}}_e}} \right) $
$ {\mathit{\boldsymbol{C}}_e} = \left( {{\mathit{\boldsymbol{u}}_e} \cdot \nabla {\mathit{\boldsymbol{N}}_e},{\mathit{\boldsymbol{N}}_e}} \right) $
$ \mathit{\boldsymbol{C}}_e^1 = \mathit{\boldsymbol{\tau }}_e^s\left( {{\mathit{\boldsymbol{u}}_e} \cdot \nabla {\mathit{\boldsymbol{N}}_e},{\mathit{\boldsymbol{u}}_e} \cdot \nabla {\mathit{\boldsymbol{N}}_e}} \right) $
$ \mathit{\boldsymbol{C}}_e^2 = \tau _e^s\left( {{\mathit{\boldsymbol{u}}_e} \cdot \nabla {\mathit{\boldsymbol{N}}_e},{\mathit{\boldsymbol{N}}_e}} \right) $
$ {\mathit{\boldsymbol{D}}_e} = \varepsilon \left( {\nabla {\mathit{\boldsymbol{N}}_e},\nabla {\mathit{\boldsymbol{N}}_e}} \right) $
$ \mathit{\boldsymbol{D}}_e^1 = \tau _e^s\varepsilon \left( {\nabla {\mathit{\boldsymbol{N}}_e},\nabla \left( {{\mathit{\boldsymbol{u}}_e} \cdot \nabla {\mathit{\boldsymbol{N}}_e}} \right)} \right) $
$ \mathit{\boldsymbol{D}}_e^2 = \tau _e^c\left( {{\mathit{\boldsymbol{K}}_{e,k}};\nabla {\mathit{\boldsymbol{N}}_e},\nabla {\mathit{\boldsymbol{N}}_e}} \right) $
$ {\mathit{\boldsymbol{E}}_e} = \left( {{\mathit{\boldsymbol{N}}_e},{\mathit{\boldsymbol{N}}_e}} \right) $
$ {\mathit{\boldsymbol{S}}_e} = \varepsilon {\left( {{\mathit{\boldsymbol{N}}_e},{\mathit{\boldsymbol{N}}_e}} \right)_{{\mathit{\Gamma }_N}}} $
$ \mathit{\boldsymbol{S}}_e^1 = \tau _e^s\varepsilon {\left( {{\mathit{\boldsymbol{N}}_e},{\mathit{\boldsymbol{u}}_e} \cdot \nabla {\mathit{\boldsymbol{N}}_e}} \right)_{{\mathit{\Gamma }_N}}} $

式中:BeBe1CeCe1Ce2DeDe1De2EeSeSe1为单元矩阵,矩阵元素可由Chebyshev多项式的性质进行计算;fe为节点右端项列向量,ge, N为边界节点列向量。利用有限元中的矩阵合成方法,即可得式(1)总的离散方程:

$ \begin{array}{*{20}{c}} {\left( {\mathit{\boldsymbol{B}} + \mathit{\boldsymbol{C}} + \mathit{\boldsymbol{D}} + {\mathit{\boldsymbol{B}}^1} + {\mathit{\boldsymbol{C}}^1} + {\mathit{\boldsymbol{D}}^1} + {\mathit{\boldsymbol{D}}^2}} \right)\boldsymbol{\phi} = }\\ {\left( {\mathit{\boldsymbol{E}} + {\mathit{\boldsymbol{C}}^2}} \right)\mathit{\boldsymbol{f}} + \left( {\mathit{\boldsymbol{S}} + {\mathit{\boldsymbol{S}}^1}} \right){\mathit{\boldsymbol{g}}_N}} \end{array} $ (12)

利用式(9)的简单线性化对D2进行线性化处理,即可通过迭代方法求得离散方程(12)的数值解。

4 解析解及边界层算例验证 4.1 解析解算例

计算区域Ω=(0, 1)×(0, 1),给定式(1)的一个解析解为

$ \phi \left( {x,y} \right) = \sin \left( {{\rm{ \mathsf{ π} }}x} \right)\sin \left( {{\rm{ \mathsf{ π} }}y} \right) $ (13)

式(1)中u=[1 1]Tγ=0、ε=10-20,采用4种不同的网格划分4×4、8×8、16×16、32×32,图 1给出了单元插值阶数和误差L2范数的关系。当单元插值阶数增加时,可以看出误差迅速减小。

Download:
图 1 不同插值阶数的误差L2范数 Fig. 1 Error L2 norm for different element interpolation order

图 2给出了不同计算节点数的误差L2范数,可以看出相同的计算节点数时,较高的单元插值阶数获得了较小的计算误差。

Download:
图 2 不同节点数的误差L2范数 Fig. 2 Error L2 norm for different node number
4.2 45°斜坡算例

设定u=[1 0]Tγ=0、ε=10-10f=1,gD=0。所有的边界层算例取定计算区域Ω=(0, 1)×(0, 1),网格划分30×30,插值阶数p=3。斜坡算例为由入口到出口的45°斜坡,其在y侧边界存在抛物边界层,在出口边界存在指数边界层,如图 3所示。稳定化谱元方法准确地反映了指数和抛物边界层的发展,没有在边界层附近产生任何振荡。

Download:
图 3 45°斜坡3D Fig. 3 45° slope 3D

同时,为了进一步说明稳定化谱元方法在计算含有边界层问题时的稳定效果,图 4给出了谱元方法和稳定化谱元方法计算45°斜坡算例,在y=0.5的计算结果。由图 4可以看出,谱元方法的数值振荡充满了整个计算区域,而稳定化谱元方法则准确地描述了出口的指数边界层。

Download:
图 4 45°斜坡,y=0.5 Fig. 4 45° slope, y=0.5
4.3 内边界层算例

设定γ=0、ε=10-10f=0,边界条件:

$ \phi \left( {x,1} \right) = 1,\frac{{\partial \phi \left( {1,y} \right)}}{{\partial x}} = \frac{{\partial \phi \left( {x,0} \right)}}{{\partial y}} = 0, $
$ \begin{array}{*{20}{c}} {\phi \left( {0,y} \right) = }\\ {\left\{ {\begin{array}{*{20}{c}} \begin{array}{l} 0,\\ y - 0.6,\\ 18\left( {y - 0.65 + 0.05} \right),\\ y + 0.25,\\ 1, \end{array}&\begin{array}{l} 0 \le y \le 0.6\\ 0.6 \le y \le 0.65\\ 0.65 \le y \le 0.7\\ 0.7 \le y \le 0.75\\ 0.75 \le y \le 1 \end{array} \end{array}} \right.} \end{array} $

分别取两种不同的流场:倾斜流场u1=[1-1]T、旋转流场u2=[y-x]T,该算例会在计算区域内形成内部边界层。图 56分别给出了倾斜流场和旋转流场形成的倾斜和圆形内边界层,可以看出稳定化谱元方法同样对内边界层进行了准确的反映。

Download:
图 5 倾斜内边界层3D Fig. 5 Inclined internal boundary layer 3D
Download:
图 6 圆形内边界层3D Fig. 6 Circular internal boundary layer 3D

图 7给出了插值阶数对倾斜内边界层计算结果的影响,y=0.5,网格划分30×30。从图 7可以看出,内边界层的逼近效果随着插值阶数的增加迅速提高,且很快地趋于一致。

Download:
图 7 不同单元插值阶数的边界层逼近 Fig. 7 Boundary layer approximation for different element interpolation order
5 结论

1) 稳定化谱元方法在谱元离散时增加了流线及边界层方向的人工黏性,使得谱元方法在求解对流项占优的对流扩散问题时能够保持数值解的稳定,同时适当的迎风函数的定义使得人工黏性的添加没有破坏谱元方法原有的谱精度和谱收敛特性;

2) 人工黏性项的增加使得谱元方法能够求解含有指数边界层、抛物边界层和内边界层的对流扩散问题,当插值阶数增加时,边界层的逼近效果迅速提高;

3) 和谱元方法相比,稳定化谱元方法极大地扩大了对流扩散问题的稳定求解范围,完全消除了计算区域及边界层邻域内的数值振荡,获得了一致稳定的结果;

4) 研究工作为谱元方法在流动及和流动相关联的对流扩散问题中的应用奠定了一定的理论基础。

参考文献
[1]
AUGERAUND-VÉRON E, CHOQUET C, COMTE É. Optimal control for a groundwater pollution ruled by a convection-diffusion-reaction problem[J]. Journal of optimization theory and applications, 2017, 173(3): 941-966. DOI:10.1007/s10957-016-1017-8 (0)
[2]
LIU Yijian, LI Ziyu. Research on locating model of heavy metal pollutants source based on SFPI method and 2D convection-diffusion equation[J]. Environment and natural resources research, 2017, 7(2): 68-79. DOI:10.5539/enrr.v7n2p68 (0)
[3]
LI Lingyu, YIN Zhe. Numerical simulation of groundwater pollution problems based on convection diffusion equation[J]. American journal of computational matheatics, 2017, 7(3): 350-370. DOI:10.4236/ajcm.2017.73025 (0)
[4]
刘忠波, 房克照, 孙昭晨. 不同时空格式在求解污染物对流扩散方程中的应用[J]. 海洋技术, 2012, 31(1): 96-99.
LIU Z B, FANG K Z, SUN Z C. The application of different time & space schemes in pollutant convective-diffusive equation[J]. Marine technology, 2012, 31(1): 96-99. (0)
[5]
牛云霞, 李春光, 景何仿, 等. 求解对流扩散方程的高阶紧致Padé有限体积格式[J]. 高等学校计算数学学报, 2017, 39(1): 43-54.
NIU Y X, LI C G, JING H F, et al. A Padé compact high order finite volume scheme for solving convection-diffusion equation[J]. Journal of computational mathematics in colleges and universities, 2017, 39(1): 43-54. (0)
[6]
祁应楠, 武莉莉. 一维定常对流扩散反应方程的高精度紧致差分格式[J]. 华中师范大学学报(自然科学版), 2017, 51(1): 1-6.
QI Y N, WU L L. A high-order compact difference scheme for the 1D steady convection-diffusion-reaction equation[J]. Journal of central China normal university (Nat. Sci.), 2017, 51(1): 1-6. (0)
[7]
朱晓钢, 聂玉峰, 王俊刚, 等. 分数阶对流扩散方程的特征有限元方法[J]. 计算物理, 2017, 34(4): 417-424.
ZHU X G, NIE Y F, WANG J G, et al. A characteristic finite element method for fractional convection-diffusion equations[J]. Chinese journal of computational physics, 2017, 34(4): 417-424. (0)
[8]
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. (0)
[9]
BROOKS A N, HUGHES T J R. Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[J]. Computer methods in applied mechanics and ergineering, 1982, 32(1/2/3): 199-259. (0)
[10]
HUGHES T J R, FRANCA L P, HULBERT G M. A new finite element formulation for computational fluid dynamics:Ⅷ. The Galerkin least-squares method for advective-diffusive equations[J]. Computer methods in applied mechanices and engineering, 1989, 73(2): 173-189. DOI:10.1016/0045-7825(89)90111-4 (0)
[11]
FRANCA L P, VALENTIN F. On an improved unusual stabilized finite element method for the advective-reactive-diffusive equation[J]. Computer methods in applied mechanics and engineering, 2000, 190(13/14): 1785-1800. (0)
[12]
SANGALLI G. Global and local error analysis for the residual-free bubbles method applied to advection-dominated problems[J]. SIAM journal onnumerical analrsis, 2000, 38(5): 1496-1522. (0)
[13]
GANESAN S, SRIVASTAVA S. ALE-SUPG finite element method for convection-diffusion problems in time-dependent domains:conservative form[J]. Applied matheatics computer, 2017, 303: 128-145. DOI:10.1016/j.amc.2017.01.032 (0)
[14]
GALEÃO A C, DO CARMO E G D. A consistent approximate upwind Petrov-Galerkin method for convection-dominated problems[J]. Computer methods applied mechanics engrgineering, 1988, 68(1): 83-95. DOI:10.1016/0045-7825(88)90108-9 (0)
[15]
GALEÃO A C, ALMEIDA R C, MALTA S M C, et al. Finite element analysis of convection dominated reaction-diffusion problems[J]. Applied numerical matheatics, 2004, 48(2): 205-222. DOI:10.1016/j.apnum.2003.10.002 (0)
[16]
PATERA A T. A spectral element method for fluid dynamics:laminar flow in a channel expansion[J]. Journal of computer physics, 1984, 54(3): 468-488. DOI:10.1016/0021-9991(84)90128-1 (0)
[17]
RANJAN R, CHRONOPOULOS A T, FENG Y. Computational algorithms for solving spectral/hp stabilized incompressible flow problems[J]. Journal of matheatics research, 2016, 8(4): 21-29. DOI:10.5539/jmr.v8n4p21 (0)
[18]
秦国良, 徐忠. 谱元方法求解二维不可压缩Navier-Stokes方程[J]. 应用力学学报, 2000, 17(4): 20-25.
QIN Guoliang, XU Zhong. A spectral element method for incompressible Navier-Stokes equations[J]. Chinese journal of applied mechanics, 2000, 17(4): 20-25. (0)
[19]
陈雪江, 秦国良, 徐忠. 谱元法和高阶时间分裂法求解方腔顶盖驱动流[J]. 计算力学学报, 2002, 19(3): 281-285.
CHEN Xuejiang, QIN Guoliang, XU Zhong. Spectral element method and high order time splitting method for navies stokes equation[J]. Chinese journal of computational mechanics, 2002, 19(3): 281-285. (0)
[20]
LUBE G, RAPIN G. Residual-based stabilized higher-order FEM for advection dominated problems[J]. Computer methods in applied mechanics and engineering, 2006, 195: 4124-4138. DOI:10.1016/j.cma.2005.07.017 (0)
[21]
ALMEIDA R C, SILVA R S. A stable Petrov-Galerkin method for convection-dominated problems[J]. Computer methods in applied mechanics engineering, 1997, 140: 291-304. DOI:10.1016/S0045-7825(96)01108-5 (0)