对流扩散方程是最为常见的运动学方程，可以用来描述大气、海洋、河流中的污染物分布，地下石油的模拟开采，多孔介质流动，流体中的热量传递等众多物理过程，因其重要性而备受关注^[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为源项；g_D、g_N分别为第一类和第二类边界条件；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)

式中：H^s(Ω)为通常的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=1^N_eΩ_e=ϕ，N_e为单元数，U^h、V^h为U、V的有限元子空间：

$ {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)

式中P_e^k为定义在Ω_e上的k阶多项式空间。

2.1 SUPG方法

式(1)在子空间的SUPG变分问题为：求ϕ^h∈U^h，使得

$ \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\} $

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

2.2 CAU方法

式(1)在子空间的CAU变分问题为：求ϕ^h∈U^h，使得

$ \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) $

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

2.3 简单线性化迭代

由于CAU方法在边界层方向增加了非线性形式的人工黏性，对流扩散方程也相应的转化为非线性方程，式(8)也应该采用迭代的方法进行求解。利用第k个迭代步的计算结果ϕ_k^h对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}}} $

由给定的ϕ_k^h，求ϕ_k+1^h∈U^h, 使得

$ \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+1^h-ϕ_k^h‖ < σ_tol，σ_tol为迭代精度，ϕ₀^h为采用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)} $

式中：I^d为插值节点数，ϕⁱ为求解变量在i节点的值，Nⁱ为i节点的插值基函数，ϕ为求解变量的节点值列向量，N为节点插值基函数列向量，T_m=cos(m arccos x)，I^x、I^y为x、y方向的插值阶数。在单元外h^j(ξ)、h^k(ξ)为0，在单元内h^j(ξ^p)=δ^jp、h^k(ξ^q)=δ^kq。c_m满足：

$ {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}}} $

式中：B_e、B_e¹、C_e、C_e¹、C_e²、D_e、D_e¹、D_e²、E_e、S_e、S_e¹为单元矩阵，矩阵元素可由Chebyshev多项式的性质进行计算；f_e为节点右端项列向量，g_{e, 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)的简单线性化对D²进行线性化处理，即可通过迭代方法求得离散方程(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给出了单元插值阶数和误差L₂范数的关系。当单元插值阶数增加时，可以看出误差迅速减小。

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

4.2 45°斜坡算例

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

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

4.3 内边界层算例

设定γ=0、ε=10^-10、f=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} $

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

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

5 结论

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

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

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

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