高分辨率格式在计算流体力学领域中已经得到了广泛的应用, 如数值模拟溃坝^[1-2]、环境污染^[3-4]、孤立水波^[5], 激波捕捉^[6]等实际问题.控制方程中的对流项看似简单, 但伪数值耗散和非物理振荡导致在数值求解上十分困难.计算流体力学中低阶格式(如一阶迎风格式)具有有界性、稳定性, 而且在间断附近不会产生数值震荡, 却会带来严重的数值耗散而导致局部的大梯度和间断抹平, 这是由于数值计算带来的而并非真实的物理现象^[7-8]；而高阶格式(如SOU格式, QUICK格式)虽然能够改进数值计算的准确性和减少数值耗散^[9], 但是在陡坡和间断附近却会产生数值振荡^[10-11].

因此, Harten^[12]提出了一种高分辨率格式, 即最小变差递减(total variation diminishing, TVD)格式.该方法能够克服上述计算格式的缺陷, 并继承了计算流体力学中的迎风特性.TVD格式通过引入流通量限制器实现在低阶格式与高阶格式之间的转换保证TVD特性.之后, 各种各样的流通量限制器(Superbee^[13], Van Leer^[14], Minmod^[12])被提出.Juntasaro等^[15]在结构网格下比较了不同限制器的性质.对于结构网格, 流通量限制器的因数r是一致的.但实际中, 特别是溃坝, 孤立波爬坡等问题, 计算域往往较为复杂的问题, 非结构网格能够快速、灵活地对计算域进行网格划分.而对于非结构网格, 通常无法直接定义迎风点, 给TVD格式的实施带来困难.

近年来, 一种源于无网格方法并采用梯度光滑技术的梯度光滑法(gradient smoothing method, GSM)被提出, 应用于求解流体问题^[16].由于有多种光滑函数与积分格式作为备选, 此方法具有更加优秀的通用性和灵活性.大量的数值实验证明了GSM能提供稳定、精确的数值解, 而且具有计算结果对网格畸变不敏感的性质^[17-19].GSM方法已经成功被用于求解非结构网格下的可压缩/不可压缩和黏性/无黏流动, 但目前尚未引入高分辨率计算格式.引入TVD格式将有助于避免应用GSM方法求解间断问题时所产生的数值振荡与耗散, 对于求解复杂工程实际问题非常重要.相对于经典的基于非结构网格的TVD格式, GSM方法中能够形成多种光滑域, 在光滑域上可以灵活、准确地定义迎风点来构造TVD格式.

在Li等^[20]构造非结构网格TVD格式的基础上, 本文提出了改进的TVD格式并用于非结构网格下求解对流方程.利用梯度光滑技术, 在3种不同的光滑域上分别定义迎风信息来构造TVD格式, 并将该方法应用到不同类型的对流方程, 进行数值实验.通过与现有的基于非结构网格TVD方法进行对比, 验证所提出的方法.

1 梯度光滑法

Liu等^[21]提出了广义梯度光滑技术, 并应用于固体力学中提出了基于G空间理论和双重弱形式的光滑点插值法.Liu等^[16]将梯度光滑技术应用于直接求解偏微分方程提出了应用于流体求解的梯度光滑法.梯度光滑法利用散布在计算域与域边界上的节点形成背景网格来表示该问题计算域及其边界, 场变量信息存储在背景网格的节点上.计算域在不同位置(如节点、边中点、单元中心点)的梯度可以在相应的光滑域上利用光滑技术求解.

1.1 梯度光滑技术

计算域Ω中任意一点X_i处场变量U的梯度可近似为

$ \nabla {U_i} \equiv \nabla U\left( {{X_i}} \right) \approx \int_\mathit{\Omega } {\nabla U\left( X \right)\hat \omega \left( {X - {X_i}} \right){\rm{d}}\mathit{\Omega }} \cdot $

(1)

式中：$ \nabla $为梯度算子, $ \begin{array}{*{20}{l}}{\hat \omega (X - {X_i})}\end{array} $为光滑函数；U(X)为场函数.

对式(1)分部积分并采用散度定理

$ \begin{array}{*{20}{l}}{\nabla U\left( {{X_i}} \right) \approx \int_\mathit{\Gamma } {U\left( X \right)\hat \omega \left( {X - {X_i}} \right)\mathit{\boldsymbol{n}}{\rm{d}}\mathit{\Gamma }} - }\\{\;\;\;\int_\Omega {U\left( X \right)\nabla \hat \omega \left( {X - {X_i}} \right){\rm{d}}} \mathit{\Omega } \cdot }\end{array} $

(2)

式中：Γ为域的边界；n为边界的外法向向量；如图 1所示.

光滑函数选取如下形式的分段函数：

$ \hat \omega \left( {X - {X_i}} \right) = \left\{ \begin{array}{l} 1/{A_i}, X \in \mathit{\Omega };\\ \;\;\;0\;{\rm{ }}\;\;\;, X \notin \mathit{\Omega } \cdot \end{array} \right.{\rm{ }} $

(3)

式中：A_i为光滑域Ω的面积.将式(3)代入式(2)中, 化简为

$ \nabla U\left( {{X_i}} \right) \approx \frac{1}{{{A_i}}}\int_\mathit{\Gamma } {U\left( X \right)\mathit{\boldsymbol{n}}{\rm{d}}\mathit{\Gamma }} \cdot $

(4)

1.2 光滑域的构造

在GSM方法中, 计算域及其边界用分布节点加以表示, 这些节点是场变量的载体.节点的分布一般为非均匀的, 节点分布密度取决计算精度的要求.通过连结这些节点来构成初始单元(背景网格).在GSM方法中提出了3种光滑域, 分别为：节点光滑域(node-based gradient smoothing domain, nGSD), 中点光滑域(midpoint-based gradient smoothing domain, mGSD)和中心点光滑域(centroid-based gradient smoothing domain, cGSD).节点光滑域的构成是通过连结相关单元中心点及其关联的边中点；中点光滑域是连结中点所在边的2个端点与边相邻的2个单元中点；中点光滑域就是计算域的初始单元, 3种光滑域如图 2所示.

1.3 空间梯度近似 1.3.1 节点的一阶梯度

计算节点一阶导数时采用一点积分格式, 并假设：

$ {U_{{c_k}}} = {U_{{c_{k - 1}}}} = {U_{{m_k}}} \cdot $

(5)

式中：U_ck, U_ck-1分别为单元中心点c_k, c_k-1的场变量值, U_mk为边中点m_k的场变量值.

在节点光滑域上应用式(4), 节点处的一阶梯度可近似为

$ \left. \begin{array}{l} \frac{{\partial {U_i}}}{{{\partial _x}}} \approx \frac{1}{{A_i^{nGSD}}}\sum\nolimits_{k = 1}^{{n_i}} {{{\left( {\Delta {\mathit{\boldsymbol{S}}_x}} \right)}_{i{j_k}}}{U_{{m_k}}}}, \\ \frac{{\partial {U_i}}}{{{\partial _y}}} \approx \frac{1}{{A_i^{nGSD}}}\sum\nolimits_{k = 1}^{{n_i}} {{{\left( {\Delta {\mathit{\boldsymbol{S}}_y}} \right)}_{i{j_k}}}{U_{{m_k}}}} \cdot \end{array} \right\} $

(6)

$ \left. \begin{array}{l} {\left( {\Delta {\mathit{\boldsymbol{S}}_x}} \right)_{i{j_k}}} = \left( {\Delta {\mathit{\boldsymbol{S}}_x}} \right)_{i{j_k}}^{\left( L \right)} + \left( {\Delta {\mathit{\boldsymbol{S}}_x}} \right)_{_{i{j_k}}}^{\left( R \right)}, \\ {\left( {\Delta {\mathit{\boldsymbol{S}}_y}} \right)_{i{j_k}}} = \left( {\Delta {\mathit{\boldsymbol{S}}_y}} \right)_{i{j_k}}^{\left( L \right)} + \left( {\Delta {\mathit{\boldsymbol{S}}_y}} \right)_{_{i{j_k}}}^{\left( R \right)} \cdot \end{array} \right\} $

(7)

式中：A_i^nGSD为节点光滑域的面积, n_i为与节点i相连的节点总数, (ΔS_x)_ijk和(ΔS_y)_ijk为边界m_kc_k的外法向向量, 上标(L)和(R)为原始单元的边ij_k相关联的两侧的光滑域边界, 如图 2(a)所示.

$ \;\;\;\left. \begin{array}{l} \;\;\;\;\left( {\Delta {\mathit{\boldsymbol{S}}_x}} \right)_{i{j_k}}^{\left( L \right)} = {r_{{m_k}{c_k}}}{\left( {{\mathit{\boldsymbol{n}}_x}} \right)_{{m_k}{c_k}}}\\ \;\;\;\;\left( {\Delta {\mathit{\boldsymbol{S}}_y}} \right)_{i{j_k}}^{\left( L \right)} = {r_{{m_k}{c_k}}}{\left( {{\mathit{\boldsymbol{n}}_y}} \right)_{{m_k}{c_k}}}\\ \;\left( {\Delta {\mathit{\boldsymbol{S}}_x}} \right)i{j_k}\left( R \right) = {r_{{m_k}{c_{k - 1}}}}{\left( {{\mathit{\boldsymbol{n}}_x}} \right)_{{m_k}{c_{k - 1}}}}\\ \;\left( {\Delta {\mathit{\boldsymbol{S}}_y}} \right)i{j_k}\left( R \right) = {r_{{m_k}{c_{k - 1}}}}{\left( {{\mathit{\boldsymbol{n}}_y}} \right)_{{m_k}{c_{k - 1}}}} \end{array} \right\} \cdot $

(8)

式中：r_mkck、r_mkck-1分别为光滑域边界m_kc_k和m_kc_k-1的长度, n_x和n_y为边界m_kc_k和m_kc_k-1的单位外法向向量.

边中点处的变量值利用线性插值计算,

$ {U_{{m_k}}} \approx \frac{{{U_i} + {U_j}_k}}{2} \cdot $

(9)

1.3.2 中点的一阶梯度中点梯度

可以在中点光滑域上利用式(4)来近似计算, 表达形式如下：

$ \begin{array}{l} \frac{{\partial {U_{{m_k}}}}}{{{\partial _x}}} \approx \left[{\frac{1}{2}{{\left( {\Delta \mathit{\boldsymbol{S}}_m^x} \right)}_{i{c_{k-1}}}}\left( {{U_i} + {U_{{c_{k-1}}}}} \right) + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}{\left( {\Delta \mathit{\boldsymbol{S}}_m^x} \right)_{{c_{k-1}}{j_k}}}\left( {{U_{_{{c_{k - 1}}}}} + {U_{_{{j_k}}}}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}{\left( {\Delta \mathit{\boldsymbol{S}}_m^x} \right)_{{j_k}{c_k}}}\left( {{U_{_{{j_k}}}} + {U_{{c_k}}}} \right) + \\ \left. {\;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}{{\left( {\Delta \mathit{\boldsymbol{S}}_m^x} \right)}_{{c_k}i}}\left( {{U_{{c_k}}} + {U_i}} \right)} \right]\frac{1}{{A_i^{mGSD}}} \cdot \end{array} $

(10)

$ \begin{array}{l} \frac{{\partial {U_{{m_k}}}}}{{{\partial _y}}} \approx \left[{\frac{1}{2}{{\left( {\Delta \mathit{\boldsymbol{S}}_m^y} \right)}_{i{c_{k-1}}}}\left( {{U_i} + {U_{{c_{k-1}}}}} \right) + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}{\left( {\Delta \mathit{\boldsymbol{S}}_m^y} \right)_{{c_{k-1}}{j_k}}}\left( {{U_{_{{c_{k - 1}}}}} + {U_{_{{j_k}}}}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}{\left( {\Delta \mathit{\boldsymbol{S}}_m^y} \right)_{{j_k}{c_k}}}\left( {{U_{_{{j_k}}}} + {U_{{c_k}}}} \right) + \\ \left. {\;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}{{\left( {\Delta \mathit{\boldsymbol{S}}_m^y} \right)}_{{c_k}i}}\left( {{U_{{c_k}}} + {U_i}} \right)} \right]\frac{1}{{A_i^{mGSD}}} \cdot \end{array} $

(11)

式中：U_i, U_jk为节点的场变量, A_i^mGSD为中点光滑域的面积, ΔS_m^x和ΔS_m^y为中点光滑域界面的外法向矢量, 计算方式与节点光滑域界面的外法向矢量相似.中心点处的场变量U_ck利用单元节点U_i, U_jk和U_jk+1的变量采用线性插值计算

$ {U_{{c_k}}} \approx \frac{{{U_i} + {U_{{j_k}}} + {U_{{j_{k + 1}}}}}}{3} \cdot $

(12)

1.3.3 中心点的一阶梯度

中心点一阶梯度的离散格式可近似为如下形式

$ \begin{array}{l} \frac{{\partial {U_{{c_k}}}}}{{{\partial _x}}} \approx \left[{\frac{1}{2}{{\left( {\Delta \mathit{\boldsymbol{S}}_c^x} \right)}_{i{j_k}}}\left( {{U_i} + {U_{{j_k}}}} \right) + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}{\left( {\Delta \mathit{\boldsymbol{S}}_c^x} \right)_{{j_k}{j_{k + 1}}}}\left( {{U_{_{{j_k}}}} + {U_{_{{j_{k + 1}}}}}} \right) + \\ \left. {\;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}{{\left( {\Delta \mathit{\boldsymbol{S}}_c^x} \right)}_{{j_{k + 1}}i}}\left( {{U_{{j_{k + 1}}}} + {U_i}} \right)} \right]\frac{1}{{A_i^{cGSD}}} \cdot \end{array} $

(13)

$ \begin{array}{l} \frac{{\partial {U_{{c_k}}}}}{{{\partial _y}}} \approx \left[{\frac{1}{2}{{\left( {\Delta \mathit{\boldsymbol{S}}_c^y} \right)}_{i{j_k}}}\left( {{U_i} + {U_{{j_k}}}} \right) + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}{\left( {\Delta \mathit{\boldsymbol{S}}_c^y} \right)_{{j_k}{j_{k + 1}}}}\left( {{U_{_{{j_k}}}} + {U_{_{{j_{k + 1}}}}}} \right) + \\ \left. {\;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}{{\left( {\Delta \mathit{\boldsymbol{S}}_c^y} \right)}_{{j_{k + 1}}i}}\left( {{U_{{j_{k + 1}}}} + {U_i}} \right)} \right]\frac{1}{{A_i^{cGSD}}} \cdot \end{array} $

(14)

式中:ΔS_c^x和ΔS_c^y为相应的中心点光滑域边界外法向向量, A_i^cGSD为中心点光滑域的面积.

2 基于非结构网格的TVD格式

光滑域边界上的流通量$ {\phi _i} + \frac{1}{2} $按照TVD格式可写成一阶迎风项和逆扩散与非线性流通量限制器的乘积之和, 结构网格下表达式如下^[23]：

$ {\phi _{i + \frac{1}{2}}} = {\phi _i} + \frac{1}{2}\psi \left( {{r_{i + \frac{1}{2}}}} \right)\left( {{\phi _{i + 1}} - {\phi _i}} \right), $

(15)

$ {r_{i + \frac{1}{2}}} = \frac{{{\phi _i} - {\phi _{i - 1}}}}{{{\phi _{i + 1}} - {\phi _i}}} \cdot $

(16)

式中:ψ(r)为流通量限制器.ψ(r)需要满足在r-ψ关系图中其图像落在二阶精度TVD单调区间内^[24].在本文中采用了3种常见的限制器：Superbee限制器ψ(r)=max (0, min (1, 2r), min (2, r)), Van Leer限制器$ \left( {\psi \left( r \right){\rm{ }} = \frac{{r + \left| r \right|}}{{1 + \left| r \right|}}} \right) $和Minmod限制器ψ(r)=max (0, min (1, r)).在非结构网格下, 变量如图 3所示.式(15)、(16)改写成

$ {\phi _f} = {\phi _C} + \frac{1}{2}\psi \left( r \right)\left( {{\phi _D} - {\phi _C}} \right) \cdot $

(17)

$ r = \frac{{{\phi _C} - {\phi _U}}}{{{\phi _D} - {\phi _C}}} \cdot $

(18)

式中：ϕ_D, ϕ_C和ϕ_U分别为下风点, 中心点和迎风点的变量值.可以看出在非结构网格上应用TVD格式, 迎风点U的变量信息无法直接得到.因此, 很多学者提出了不同的方法在非结构网格下构造TVD格式.

2.1 现有非结构网格的TVD格式

在构造非结构网格的TVD格式中, 主要在于计算迎风点处的变量信息和因数r.Darwish等^[25]假设在线段UD上, 变量梯度为$ \nabla {\phi _{\rm{C}}} $, 利用节点梯度$ \nabla {\phi _{\rm{C}}} $来外插计算迎风点变量信息并应用到任意网格.Li等^[20]指出了Darwish假设的不合理之处, 而采用更加准确的内插方式来计算ϕ_U, 得到了较好的计算结果.除了通过插值计算计算迎风变量, Zhang等^[26]在基于顶点的有限体积法下将TVD格式应用到非结构网格, 采用临近的节点近似得到迎风点变量.Hou等^[27]为消除在任意非结构网格下无法保证l_CD与所需计算的单元界面正交所带来的误差, 采用投影近似的方式对C、D进行重构得到新的C′和D′.Zhang等^[28]结合了Li和Hou的方法, 采用和Hou相同方式重构CD, 但在计算U点变量采用Li的内插值计算迎风点变量.由于GSM在计算梯度中依赖于构造的光滑域, 使得Hou的重构方法难以在GSM中直接实施.

可以看出, 非结构网格下迎风点的变量直接影响到整个计算结果的准确性.

Darwish等^[25]通过节点C的梯度及节点D值外插计算ϕ_U.

$ {\phi _U} = {\phi _D} - 2{\mathit{\boldsymbol{d}}_{CD}} \cdot {\left( {\nabla \phi } \right)_C} \cdot $

(19)

式中：$ {\left( {\nabla \phi } \right)_{\rm{C}}} $为中心点C的梯度, d_CD为线段CD的矢径.将上式代入式(18)中, 得到因数^[25]

$ {r_{{\rm{Darwish}}}} = \frac{{2{{(\nabla \phi )}_C} \cdot {\mathit{\boldsymbol{d}}_{CD}}}}{{{\phi _D} - {\phi _C}}} - 1 \cdot $

(20)

Li等^[20]在研究中发现Darwish方法对ϕ_U的插值计算并不准确.当变量ϕ呈抛物线分布时, 根据公式(19)计算上风点的变量ϕ_U是精确的；而当ϕ呈指数分布时, 会导致了对ϕ_U计算的误差较大.

与Darwish通过外插方式计算迎风点的变量值不同, Li提出了因数r的改进算法.假设l_CD=l_DU来确定上风点U在计算域中所位于的单元, 在单元内通过内插值的方式计算迎风点的变量值^[20], 表达式如下：

$ {\phi _U} = {\phi _{{U_r}}} - 2{\mathit{\boldsymbol{d}}_{U{U_r}}} \cdot {(\nabla \phi )_U}_r \cdot $

(21)

在Li关于r的计算格式中应用了更多的迎风信息, 将式(21)代入到公式(18)中, 可得

$ {r_{Li}} = - \frac{{2{\mathit{\boldsymbol{d}}_{U{U_r}}}\cdot{{(\nabla \phi )}_{{U_r}}}}}{{{\phi _D} - {\phi _C}}} + \frac{{{\phi _C} - {\phi _{{U_r}}}}}{{{\phi _D} - {\phi _C}}} \cdot $

(22)

式中: U_r为点U所在单元的中心点变量, $ {\left( {\nabla \phi } \right)_{{U_r}}} $为变量的梯度, d_UUr为点U到点U_r的矢径, 如图 4所示.在结构网格下, 式(22)与(18)是一致的.

与上述的2种方法通过梯度插值计算得到迎风点的变量信息不同, Barth等^[29]通过满足单调性准则的方式来构造TVD格式.Speikreijse单调性准则被改进来重构控制体上的 ϕ, 要求控制体上的变量ϕ必须在周围相邻单元ϕ的最大值与最小值范围内.构造出了不同表达形式的TVD格式并且被成功地应用到实际问题^{[30, 31]}.

Speikreijse单调性准则被改写成如下形式：

$ \begin{array}{l} {\rm{min}}\left( {{\phi _N}, {\phi _P}} \right) \le {R_P}\left( {{r_j}} \right) \le {\rm{max}}\left( {{\phi _N}, {\phi _P}} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\forall N \in {\rm{Neighbors}}\left( P \right) \cdot \end{array} $

(23)

式中：j为控制体P内的任意点, R为重构多项式,

$ {R_P}\left( {{r_j}} \right) = {\phi _P} + \psi {\left( {\nabla \phi } \right)_P} \cdot \left( {{\mathit{\boldsymbol{d}}_f} - {\mathit{\boldsymbol{d}}_P}} \right) \cdot $

(24)

式中：ψ为流通量限制器, d_f和d_P为坐标向量.控制体界f面积分点处的变量为

$ {R_P}\left( {{\phi _f}} \right) = {\phi _P} + \psi {\left( {\nabla \phi } \right)_P} \cdot \left( {{\mathit{\boldsymbol{d}}_f} - {\mathit{\boldsymbol{d}}_P}} \right) \cdot $

(25)

满足条件:min(ϕ_N, ϕ_P)≤R_P(r_f)≤max(ϕ_N, ϕ_P)

式(25)中的流通量限制器ψ表达式为^[29]

$ {\psi _f} = \left\{ \begin{array}{l} \mathit{\Phi }\left( {\frac{{({\rm{max}}\left( {{\phi _N}, {\phi _P}} \right) - {\phi _P}}}{{{{\left( {\nabla \phi } \right)}_P} \cdot \left( {{\mathit{\boldsymbol{d}}_f} - {\mathit{\boldsymbol{d}}_P}} \right)}}} \right), \;\;\;{\phi _f} > {\phi _P};\\ \mathit{\Phi }\left( {\frac{{{\phi _P} - {\rm{min}}\left( {{\phi _N}, {\phi _P}} \right)}}{{{{\left( {\nabla \phi } \right)}_P} \cdot \left( {{\mathit{\boldsymbol{d}}_f} - {\mathit{\boldsymbol{d}}_P}} \right)}}} \right)\;\;\;\;\;{\phi _f} < {\phi _P};\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;1, \;\;\;\;\;\;\;\;\;\;\;\;\;{\phi _f} = {\phi _P} \cdot \end{array} \right. \cdot $

(26)

式中：Φ(x)=min(x, 1).定义控制体上所有界面的ψ_f中的最小值作为流通量限制器ψ_P=min(ψ_f).

TVD格式应用于传统的有限体积法时, 变量值存储和计算都是在单元的中心点；而在GSM中, 存储和计算都是在节点上实施同时需要在初始单元上形成节点光滑域.因此, 应用到GSM中时都需要一定的修改.类似地, 针对于这一点Zhang等^[32]对比分析了非结构网格下不同TVD格式利用单元顶点方式和单元中心点方式有限体积法求解输运方程计的算结果.

2.2 基于光滑技术改进的非结构网格TVD格式

在Li方法基础上, 提出了不同对上风点U的插值计算格式.在3种不同的光滑域上计算迎风点上的变量, 计算步骤如下：

1)延长边DC到点U, 使l_CD=l_DU, 可得U点坐标(x_U, y_U), 由此确定U点在计算域中所在的单元.

2)在U点所在的单元分别根据建立的节点光滑域、中点光滑域、中心点光滑域提出nGSM、mGSM、cGSM3种方案来插值计算U点的变量值.

nGSM (node gradient smoothed method)：计算U点到单元3个节点的距离, 选择距U最近的节点N_min, 根据节点光滑域计算该节点梯度来插值计算ϕ_U, 如图 5(a)所示.

$ {\phi _U} = {\phi _{{N_{{\rm{mim}}}}}} + {{\boldsymbol{d}}_{{N_{{\rm{mim}}}}U}} \cdot {(\nabla \phi )_{{N_{{\rm{mim}}}}}} \cdot $

(27)

式中：d_NmimU为N_mimU的矢径.

mGSM (midpoint gradient smoothed method)：计算U点到单元3条边中点的距离, 选择距U最近的边中点M_min, 根据中点光滑域计算该中点梯度来插值计算ϕ_U, 如图 5(b)所示.

$ {\phi _U} = {\phi _{{M_{{\rm{mim}}}}}} + {{\boldsymbol{d}}_{{M_{{\rm{mim}}}}U}} \cdot {(\nabla \phi )_{{M_{{\rm{mim}}}}}} \cdot $

(28)

式中：d_MmimU为M_mimU的矢径.

cGSM (centre gradient smoothed method)：计算U点到单元中心点C₀距离, 根据中心点光滑域计算节点梯度来插值计算ϕ_U, 如图 5(c)所示.

$ {\phi _U} = {\phi _{_{{C_0}}}} + {\mathit{\boldsymbol{d}}_{{C_0}U}} \cdot {(\nabla \phi )_{{C_0}}} \cdot $

(29)

式中：d_C0U为C₀U的矢径.

(3)根据上面3种方案算出了迎风点ϕ_U, 利用式(17)、(18)计算光滑域边界上积分点的流通量, 利用式(6)在节点上建立节点光滑域进行空间离散.

可以看出在3种光滑域上插值计算迎风点ϕ_U方法是相似的.但由于在插值计算ϕ_U时所应用的光滑域不同, 在实例计算中的数值振荡与耗散也是不同的.在计算公式上, 步骤(2.3)中cGSM的ϕ_U的计算方法与Li的方法是基本相同的, 但2种方法计算$ {\left( {\nabla \phi } \right)_C} $的方法是不同的.在传统有限体积法中, 变量信息是存储在单元中心点上, 计算$ {\left( {\nabla \phi } \right)_C} $需要周围单元的信息；而在GSM中, 变量信息是存储在节点上, $ {\left( {\nabla \phi } \right)_C} $的计算仅需要本单元上的节点信息, 这就是2种方法的根本区别.如果将GSM中形成的光滑域看作初始单元的重叠网格, 利用步骤(2.1)nGSM方法计算的ϕ_U与Li的方法是相同的, 在计算结果中均以nGSM/Li标识.

当迎风点U落在计算域之外时, 计算线段UD与边界的交点$ \tilde U $来近似为迎风点, 并在边界上线性插值计算$ {\phi _{\tilde U}} $, 如图 6所示.

3 数值实验

为了验证本文提出的基于非结构网格的改进TVD格式, 通过经典算例来与应用Darwish, Li和BJ方法进行对比.算例中既包括间断问题和连续问题, 也包括稳态问题和瞬态问题.在数值计算中分别采用Superbee, Van Leer和Minmod3种常见的流通量限制器.数值计算的误差为

$ E = \frac{1}{N}\sqrt {\sum\nolimits_{i = 1}^N {{{\left( {\phi _i^n - \phi _i^e} \right)}^2}} } \cdot $

(30)

式中：N为计算域节点总数, ϕ_iⁿ为数值解, ϕ_i^e为解析解.

3.1 单阶梯对流

该问题可以检验不同TVD构造方案在横风方向上的数值耗散与振荡, 控制方程如下：

$ \frac{{\partial \phi }}{{\partial t}} + \frac{{\partial u\phi }}{{\partial x}} + \frac{{\partial v\phi }}{{\partial y}} = 0 \cdot $

(31)

式中：ϕ为输运变量, 速度场为V=(1, 1), 计算模型如图 7(a)所示.

计算域几何形状为1×1的正方形, 在计算域内任意布置7 481个节点、14 640个三角形单元, 如图 7(b)所示.

为了验证提出的改进TVD格式, 将数值结果与前人构造的TVD格式进行对比分析, 在计算中所有TVD格式均采用了Superbee流通量限制器(除BJ格式以外), 如图 8(a)所示.

图 8(a)表明Darwish和nGSM/Li的比值构造方案均在阶梯附近产生了较大的数值振荡.BJ的计算结果并未产生任何数值振荡, 但在横风方向上数值耗散最大的方法.Darwish方法在数值耗散不严重, 但存在数值振荡.基于梯度光滑技术提出的mGSM和cGSM 2种方案基本没有产生数值振荡, 在数值耗散方面也优于之前的方法.

对提出的3种不同的改进方案应用不同的流通量限制器, 计算结果如图 8(b)~(d)所示.Superbee限制器能够最大程度地保证结果的间断性, 最接近真实解, 而应用于nGSM时会产生较大的数值振荡；Minmod限制器会产生了最大的数值耗散；而Van Leer限制器则介于两者之间.此外, 提出的3种方案应用不同的限制器对流场完成初始化后, 收敛误差ε_con均能随着时间步数N快速地达到收敛, 如图 9所示.

3.2 双阶梯对流

此计算模型与单阶梯对流基本相同, 但是在计算域中具有双阶梯, 计算模型如图 10所示.网格布置与时间步长和上一个算例的参数相同.控制方程与上节相同, 边界条件如下：

$ \left. \begin{array}{l} \phi = 1\;\;\;x = 0, {\rm{ }}\;\;\;0 \le y \le 0.3;\\ \phi = 0\;\;\;\;x = 0, {\rm{ }}\;\;0.3 \le y \le 1 \cdot \end{array} \right\} $

(32)

双阶梯对流对于TVD格式提出了更为严格的要求, 进一步检验改进的TVD格式.与单阶梯对流算例相似, 将数值结果与前人构造的TVD格式进行对比分析, 如图 11(a)所示.

从图 11的计算结果中可以看出, 与单梯度对流的计算结果相似, 所有的方法中nGSM/Li在2个间断附件产生了最为严重的数值振荡, BJ方法的数值耗散最大, Darwish方法产生了数值振荡, mGSM在双阶梯的一侧产生了轻微的振荡, cGSM依然能够得到令人满意的结果.应用其他替他限制器代替Superbee限制器, nGSM和mGSM能够消除数值振荡.Minmod限制器会带来最大的数值耗散.计算收敛性如图 12所示.

3.3 风暴对流

为了验证改进TVD方法对于计算连续场问题的准确性, 以计算风暴形成旋转流场过程为例.初始状态(t=0)为场变量从计算域顶端正值逐渐过渡到低端的负值, 如图 13(a)所示.场变量随时间变化的解析解表达式如下：

$ \phi = - {\rm{tanh}}\left( {0.5\left[{y{\rm{cos}}\left( {\omega t} \right)-x{\rm{sin}}\left( {\omega t} \right)} \right]} \right) \cdot $

(33)

式中:速度场为V=(-ωy, ωx), 角速度为ω=Q/dmax (Q), 距风暴中心的距离为$ d = \sqrt {{x^2} + {y^2}} $, 切向流速为Q=sech²(d)tanh(d).

计算域为一个8 m×8 m的正方形.在计算域的中心处设置为细网格, 单元大小为0.05 m；而远离中心处选择粗网格, 单元大小为0.1 m, 如图 13(b)所示.计算时间步长为0.01 s, 初始条件由ϕ_t=0计算得到.

在旋转流场的作用下, 计算域中心处场变量的梯度随着时间推移而增大.当t=9 s时, 对比截面(y=0)上, 应用不同TVD方法的计算结果如图 14(a)所示.除了BJ方法, 与解析解对比其他方法均能够得到令人满意的计算结果.BJ方法产生了较大的数值耗散.其他的方法在远离中心处的梯度变化较大的位置造成了轻微的振荡, 随着靠近中心的位置存在一定的数值耗散.综合来看cGSM是最为接近解析解.应用不同限制器的计算结果如图 14(b)~(d)所示, 与之前的计算结果相似, Superbee限制器更加接近解析解, 而其他2种限制器却产生了不同程度的较大数值耗散.此例证明了本文基于光滑技术改进的非结构网格的TVD格式能够解决连续场变量的问题, 其中3种改进格式中cGSM依旧最为出色.

根据式(30)对本文提出的3种改进TVD格式对3个数值算例进行准确性对比, 如表 1所示(表中S, V, M分别表示Superbee, Van Leer, Minmod这3种流通量限制器).可以看出mGSM和cGSM这2种改进格式在所有的算例中都具有较高的精度, 其中cGSM改进格式采用Superbee流通量限制器时计算结果最为优秀.

4 非均匀网格验证

在提出的3种改进TVD格式中, cGSM格式在上面的3个算例中表现最为优秀, 能够准确模拟数值间断.对cGSM格式在相对非均匀的非结构网格下的计算精度进行验证.

以3.1节的单阶梯对流为例, cGSM格式采用Superbee限制器进行非均匀网格的验证及分析.本文选取了2种网格.网格A采用大小为0.15的相对均匀网格, 节点总数为5 023, 如图 15(a)所示.网格B将计算域分为3部分, 上部和下部处采用的网格尺寸为0.04, 而计算域中间的网格大小为0.01, 网格节点总数为4 527, 网格尺寸具有明显差别, 如图 15(b)所示.对比数值结果, 可以看出在相对非均匀网格B下计算结果并未产生任何数值振荡；但受网格非均匀性的影响, 在间断附近结果不如网格A结果光顺, 准确性略也差于网格A, 如图 16所示.总体上, cGSM在非均匀网格下依然能够保证在间断附近避免产生数值振荡.

5 结 语

根据梯度光滑技术, 在3种不同的光滑域上分别构造了非结构网格TVD格式, 灵活地定义迎风点, 讨论了3种构造方法的准确性.之前多数方法都是通过外插方式定义迎风点, 而本文首先确定迎风点所在的单元, 之后计算单元节点、边中点和单元中心点的梯度, 分别通过梯度插值计算迎风点的变量.在数值实验中对比了本文提出的基于nGSM、mGSM和cGSM构造的3种改进TVD格式和传统TVD格式在非结构网格下的计算结果.nGSM在3种格式中结果最差, 表现为采用Superbee限制器时在间断附近产生了数值振荡, 采用其他两种限制器时会带来较大的数值耗散.mGSM和cGSM利用Superbee限制器在计算无论间断问题还是连续问题中均优于之前的方法, 其中mGSM在能够基本不产生数值振荡或轻微数值的情况下依然保持较小的数值耗散, cGSM在不带来数值振荡的同时保证了间断解的锐利性.综合来看, cGSM要略优于mGSM.在应用Var Leer与Minmod限制器时, cGSM相对于其他的方法也能够取得较高的精度.通过对cGSM格式进行了网格非均匀性的验证, 发现对比相对均匀的非结构网格, 非均匀网格的计算结果精度略差, 但仍能够保证在间断处避免产生数值振荡.

根据以上的讨论可以看出, 通过合理的插值格式计算迎风点变量ϕ_U对于非结构网格的TVD格式非常重要.本文提出的TVD构造格式同样可用于基于非结构网格的其他高阶计算格式, 将在后续工作中深入研究.