孔令发, 董义道, 刘伟^,2)国防科技大学空天科学学院, 长沙410073

THE INFLUENCE OF GLOBAL-DIRECTION STENCIL ON GRADIENT AND HIGH-ORDER DERIVATIVES RECONSTRUCTION OF UNSTRUCTURED FINITE VOLUME METHODS ¹⁾

Kong Lingfa, Dong Yidao, Liu Wei^,2)College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China

通讯作者: 2)刘伟, 教授, 主要研究方向: 计算流体力学与应用. E-mail:fishfather6525@sina.com

收稿日期:2020-03-25接受日期:2020-03-25网络出版日期:2020-09-18

基金资助:

1)国家重大资助项目.GJXM92579

Received:2020-03-25Accepted:2020-03-25Online:2020-09-18
作者简介 About authors


摘要
模板选择方式对非结构有限体积方法的计算准确性会产生显著影响. 在之前的工作中, 基于局部方向模板存在的问题, 我们探索了一种更加简单有效的全局方向模板选择方法, 并将其应用于二阶精度非结构有限体积求解器. 基于该方法找到的模板单元均沿着壁面法向与流向, 可有效捕捉流场变化, 反映流动的各向异性, 并且模板选择过程脱离了对网格拓扑的依赖, 避免了局部方向模板选择方法中复杂的阵面推进与方向判断过程, 克服了在大压缩比三角形网格上模板单元偏离壁面法向的现象, 同时在二阶精度求解器上得到了较高的计算精度与计算准确性. 为了进一步验证全局方向模板在高阶精度非结构有限体积方法中应用的可行性, 本文初步测试了该模板对变量梯度及高阶导数重构的影响. 经检验, 在不同类型的网格上, 采用全局方向模板得到的变量梯度与高阶导数误差明显低于局部方向模板, 同时也低于共点模板的计算误差. 此外, 在高斯积分点处由全局方向模板得到的变量点值与导数误差同样在三种模板中最低. 因此该模板选择方法在非结构有限体积梯度与高阶导数重构方面具有较好的数值表现, 具备在高阶精度非结构有限体积求解器中应用并推广的可行性.
关键词： 高阶精度非结构有限体积方法;$k$-exact重构算法;模板选择方法;全局方向;局部方向

Abstract
The accuracy of unstructured finite volume methods is greatly influenced by different stencils. In previous work, based on the existing problems of local-direction stencil, we explored a more concise global-direction stencil selection method for the second-order unstructured finite volume solver, and stencil cells selected by this novel stencil selection method are always along the boundary normal and circumferential directions even on grids with high aspect ratio. As a result, the variation of flowfield is effectively captured, and flow anisotropy are well reflected. In addition, the novel method is topology-independent, since global directions are determined by the flowfield, while the local directions are strongly coupled with the grid. Therefore, the complex process of advancing front as well as local directions estimation are completely avoided in the novel stencil selection method, and the phenomenon that stencil cells deviate from the boundary normal vector is effectively eliminated on high-aspect-ratio triangular grids. What's more, a better computational accuracy and lower truncation errors on the second-order accurate finite volume solver are obtained by the employment of global-direction stencil. In order to further test the effectiveness of global-direction stencil on high-order unstructured finite volume methods, we will preliminarily utilize this stencil to test the effect of gradient and high-order derivatives reconstruction. After verification, computational errors of global-direction stencil are lower than that of local-direction stencil, and also lower than that of commonly used vertex-neighbor stencil on different grid types. Besides, errors of variable and derivatives at the Gauss point obtained by global-direction stencil are also the lowest among three methods we tested. Therefore, the global-direction stencil is well performed on gradient as well as high-order derivatives reconstruction, and it is feasible to extend this novel stencil selection method to high-order unstructured finite volume methods.
Keywords：high-order unstructured finite volume methods;$k$-exact reconstruction algorithm;stencil selection method;global direction;local direction


PDF (6125KB)元数据多维度评价相关文章导出EndNote|Ris|Bibtex收藏本文
本文引用格式
孔令发, 董义道, 刘伟. 全局方向模板对非结构有限体积梯度与高阶导数重构的影响 ¹⁾. 力学学报[J], 2020, 52(5): 1334-1349 DOI:10.6052/0459-1879-20-093
Kong Lingfa, Dong Yidao, Liu Wei. THE INFLUENCE OF GLOBAL-DIRECTION STENCIL ON GRADIENT AND HIGH-ORDER DERIVATIVES RECONSTRUCTION OF UNSTRUCTURED FINITE VOLUME METHODS ¹⁾. Chinese Journal of Theoretical and Applied Mechanics[J], 2020, 52(5): 1334-1349 DOI:10.6052/0459-1879-20-093

引言

非结构网格是CFD数值模拟较常用的一类网格, 其生成方便快捷, 能够在较少的人为干预下自动生成^[1-2], 这大幅度缩减了网格的生成周期; 同时非结构网格具备较高的灵活性, 其能快速高效地处理网格运动、网格变形以及网格加密等问题, 并对不同外形具有较强的适应性^[3-4]. 因此, 相比多块结构化网格, 非结构网格的使用带来了极大的便利. 但当前制约非结构网格在复杂外形上应用的主要原因是, 其很难保证计算效率与求解精度^[5-7]. 因此, 为了充分利用非结构网格自动化的生成能力, ****们不断尝试改进离散算法来突破精度损失以及稳定性恶化的瓶颈, 实现网格生成自动化和数值模拟精准化的统一.

当前常用的非结构网格算法是二阶精度有限体积方法^[8-13], 该方法被用于很多知名软件与程序库, 如: Ansys's Fluent (http://www.ansys.com), 以及开源软件OpenFOAM (http://openfoam.org) 等^[14-17]. 在众多不同的有限体积离散方法中, 单元中心型^[8-9,12]的非结构有限体积方法应用最为广泛, 该方法的离散过程相比单元顶点型更加简便^[18-19]. 但不论对于单元中心型还是单元顶点型的有限体积离散方法, 梯度重构是影响求解器求解精度的重要过程^[20-24]. 对于二阶精度有限体积方法, 目前较为常见的梯度重构方法包括格林高斯与最小二乘方法两类^[16]. 但为了提高计算效率与对流场的分辨率, 一些****尝试将非结构有限体积方法向高精度推广. 关于高精度非结构有限体积算法的构造, Barth等^[25-26]基于最小二乘方法, 构造了$k$-exact 重构算法. 此后, Olliver-Gooch^[27-33]基于$k$-exact 重构, 发展了一系列高精度非结构有限体积离散方法. 一般而言, 随着重构精度的增加, 所需要的网格单元数量, 即模板个数也会随之增多, 并且采用不同的模板选择方法, 对梯度与高阶导数重构乃至最终流场计算结果会产生重要影响.

考虑到非结构网格不具备结构网格的笛卡尔索引, 所以在梯度重构的过程中只能依赖于一定的准则来寻找模板单元^[34]. 目前较为常用的模板选择方法是选择与中心单元共面或共点的网格单元来构建模板集合. 这两种模板的数量可控, 但选择方式均依赖于固定的网格拓扑关系, 无法在模板选择过程中体现出流动的变化特征. 如边界层流动, 该流动沿着壁面法向的变化剧烈, 因此在重构过程中, 我们更希望模板单元能够尽可能多地体现沿壁面法向的流动信息. 考虑到这个因素, 2018年Xiong等^[35]提出了基于局部方向的模板选择方法, 该方法在各向同性以及压缩比较小的三角形网格上可近似沿着壁面法向与切向来扩充网格单元, 能够较为准确地反映流场变化, 并且模板单元的扩充严格依赖于两个局部方向, 具备类似于结构网格沿着$i$, $j$两个索引方向求解的特点, 并试图通过引入这种结构化特性^[34]来改善基于非结构网格算法的数值表现. 经过分析, 这种模板选择方式虽然初步考虑了流场的特征, 但其单元局部方向的确定严格依赖于阵面推进过程(advancing front technique)^[36-37], 该过程在方向判断时较为复杂. 其次, 在高度各向异性的三角形网格上, 基于此方法找到的模板单元与壁面法向间的偏差较大, 这将对梯度重构以及流场计算过程产生不利影响.

基于局部方向模板存在的问题, 在前期的工作中, 针对二阶精度有限体积方法探索了一种新的模板选择方法, 并将其命名为全局方向模板. 该方法不再依赖于每个网格单元独立的局部方向, 而是求解域的两个全局方向, 即壁面法向与切向, 并且模板单元严格沿着全局方向依次扩充. 该方法的优势在于, 找到的模板单元均沿着两个相互垂直的全局方向, 其空间延展性较好, 具备更加明显的结构化特征^[34]. 更重要的是, 全局方向模板选择方法不再依赖于网格拓扑约束, 即使在大压缩比三角形网格上也能够较为准确地反映流动的各向异性. 此外, 相比局部方向模板选择方法, 该方法更为简便, 避免了复杂的阵面推进与方向判断过程. 经算例验证, 采用全局方向模板得到的计算误差相比局部方向模板更低, 并且也接近甚至低于常用的共点模板.

为了进一步验证全局方向模板在高阶精度有限体积方法中应用的可行性, 在此工作中, 我们将初步测试采用全局方向模板对梯度以及高阶导数重构准确性的影响. 文章的总体结构如下: 第1节主要介绍了高阶精度非结构有限体积离散方法以及$k$-exact 重构算法的基本过程; 第2节回顾了几种常用的模板选择方式, 包括共点模板、共面模板, 同时重述与分析了局部方向模板选择方法的基本思路, 并给出全局方向模板选择方法的具体过程和在各向异性三角形网格上的表现效果; 在文章的第3节, 为了验证全局方向模板在高阶精度有限体积方法上应用的可行性, 我们在不同大压缩比三角形网格上设置典型数值算例, 来比较几种不同模板对梯度与高阶导数重构精度及准确性的影响, 并在此基础上对比了由不同模板得到的高斯点处的变量点值与导数误差; 第4节给出了本文结论与下一步工作展望.

为验证全局方向模板在高阶精度有限体积方法上应用的可行性, 本文初步测试了该模板对变量梯度以及高阶导数重构的影响, 以期为高精度非结构有限体积离散提供一种更为简单有效的模板选择方式.

1 非结构有限体积方法

1.1 控制方程及其离散

针对无黏流动, 积分型欧拉控制方程为

(1) $ \begin{eqnarray} \frac{\partial }{{\partial t}}\iiint\limits_V {{ u}{\rm d}V + \iiint\limits_V {\nabla \cdot {{ F}_{\rm c}}({ u})}{\rm d}V = 0} \label{eq1} \end{eqnarray}$
基于高斯积分公式, 可将方程(1)中的散度体积分转化为环面积分

(2) $ \begin{eqnarray} \frac{\partial }{{\partial t}}\iiint\limits_V { u}{\rm d}V + \mathop{{\int\!\!\!\!\!\int}\mkern-21mu \bigcirc}\limits_{\partial V} {{{F}_{\rm c}}({ u}) \cdot { n}{\rm d}S} = 0 \label{eq2} \end{eqnarray}$
式中, ${ u}$为守恒变量, ${ F}_{\rm c}({ u})$为对流通量, ${ n}$为控制体边界面外法向矢量, $V$与$\partial V$分别代表积分域与积分域边界. 对于高阶有限体积方法, 该控制方程可被离散为

(3) $ \begin{eqnarray} \frac{\partial }{{\partial t}}{\left( {{\bar{ u}}V} \right)_i} + {\left( {\sum\limits_{j = 1}^{{N_{\rm f}}} {\sum\limits_{k = 1}^{{N_{\rm G}}} {{\omega _k}{{ {\varPhi }}_{jk}} \cdot {{{ n}}_j}{S_j}} } } \right)_i} = 0 \label{eq3} \end{eqnarray}$
对于控制体$i$, $N_{\rm f}$代表单元面数量. 由于高阶有限体积离散时, 仅构造单元面中心点处的通量无法达到二阶以上精度, 因此在每个单元面处随着精度要求的不同, 所需要的高斯积分点的数量也不同. $N_{\rm G}$代表当前单元面处高斯积分点的数量, $\omega_k$与$ {\varPhi}_{jk}$分别为每个高斯积分点处的权系数与数值通量, ${ n}_j$与$S_j$为单元面处的单位外法矢量与单元面面积(二维时为单元边长). 图1给出了有限体积离散与通量构造的基本过程.

图1

新窗口打开|下载原图ZIP|生成PPT
图1非结构有限体积离散与通量构造示意图

Fig.1Unstructured finite volume discretization and flux construction



如图1所示, 为采用近似黎曼方法^[30-31]构造单元面某高斯积分点处的数值通量, 需要首先构造单元面处的左右状态矢量.

(4) $ \begin{eqnarray}\label{eq4} \left.\begin{array}{l} {{ u}_{\rm GL}^n} = {{ u}_i} + {\left. {\dfrac{{\partial { u}}}{{\partial x}}} \right|_i}\left( {{x_i} - {x_{\rm G}}} \right) + {\left. {\dfrac{{\partial { u}}}{{\partial y}}} \right|_i}\left( {{y_i} - {y_{\rm G}}} \right)+\\ \qquad\dfrac{1}{2}{\left. {\dfrac{{{\partial ^2}{ u}}}{{\partial {x^2}}}} \right|_i}{\left( {{x_i} - {x_{\rm G}}} \right)^2} + \dfrac{{{\partial ^2}{ u}}}{{\partial x\partial y}}\left( {{x_i} - {x_{\rm G}}} \right)\left( {{y_i} - {y_{\rm G}}} \right)+\\ \qquad\dfrac{1}{2}{\left. {\dfrac{{{\partial ^2}{ u}}}{{\partial {y^2}}}} \right|_i}{\left( {{y_i} - {y_{\rm G}}} \right)^2} + \cdots\\ {{ u}_{\rm GR}^n} = {{ u}_j} + {\left. {\dfrac{{\partial { u}}}{{\partial x}}} \right|_j}\left( {{x_j} - {x_{\rm G}}} \right) + {\left. {\dfrac{{\partial { u}}}{{\partial y}}} \right|_j}\left( {{y_j} - {y_{\rm G}}} \right) +\\ \qquad\dfrac{1}{2}{\left. {\dfrac{{{\partial ^2}{ u}}}{{\partial {x^2}}}} \right|_j}{\left( {{x_j} - {x_{\rm G}}} \right)^2} + \dfrac{{{\partial ^2}{ u}}}{{\partial x\partial y}}\left( {{x_j} - {x_{\rm G}}} \right)\left( {{y_j} - {y_{\rm G}}} \right)+\\ \qquad\dfrac{1}{2}{\left. {\dfrac{{{\partial ^2}{ u}}}{{\partial {y^2}}}} \right|_i}{\left( {{y_j} - {y_{\rm G}}} \right)^2} + \cdots \end{array}\right\}\quad \end{eqnarray}$
式中, ${ u}_{\rm GL}^n$与${ u}_{\rm GR}^n$分别代表当前时间步高斯积分点处的左右状态矢量, $x_{\rm G}$, $y_{\rm G}$; $x_i$, $y_i$与$x_j$, $y_j$分别代表高斯积分点以及单元$i$ 与单元$j$的中心点坐标值. 在得到左右状态矢量的基础上, 可采用Roe格式^[38-39]构造数值通量, 其基本表达式如下

(5) $ \begin{eqnarray} {{ {\varPhi}}_{{c_k}}} = \frac{1}{2}\left[ {{{ {F}}_n}({{ u}}_{\rm L}^n) + {{{ F}}_n}({{ u}}_{\rm R}^n)} \right] - \frac{1}{2}\left| {{\bar{ A}}({{ u}}_{\rm L}^n,{{ u}}_{\rm R}^n)} \right|({{ u}}_{\rm L}^n - {{ u}}_{\rm R}^n) \label{eq5}\quad \end{eqnarray}$
式中$\left| {{\bar{ A}}({{ u}}_{\rm L}^n,{{ u}}_{\rm R}^n)} \right|$代表经Roe平均后的雅可比矩阵.

1.2 $k$-exact 重构算法

由于本文研究全局方向模板在高阶精度非结构有限体积方法中的应用, 因此采用$k$-exact 方法来实现梯度与高阶导数重构. 此过程需要利用中心控制体单元$i$的模板单元来构造重构函数${R_i}({{ x}} - {{{ x}}_i})$, 其包含单元$i$的平均值, 并且任意多项式都可被函数${R_i}({{ x}} - {{{ x}}_i})$在不大于$k$阶精度的条件下精确重构, 即

(6) $ \begin{eqnarray} {R_i}({{ x}} - {{{ x}}_i}) - {{ u}}({{ x}}) = O(\Delta {x^{k + 1}}) \label{eq6} \end{eqnarray}$
但在重构过程中, 为达到重构精度, 始终需保证重构函数与待重构变量在控制体单元上的几何平均值相等

(7) $ \begin{eqnarray} \frac{1}{{{A_i}}}\int_{{V_i}} {{R_i}({{ x}} - {{{ x}}_i}){\rm d}A = } \frac{1}{{{A_i}}}\int_{{V_i}} {{{ u}}({{ x}}){\rm d}A} \equiv {{\bar{ u}}_i} \label{eq7} \end{eqnarray}$
式中$V$代表控制体. 对于$k$-exact 重构, 其重构函数的形式为

(8) $ \begin{eqnarray}\label{eq8} &&{R_i}({{ x}} - {{{ x}}_i}) = {\left. {{ u}} \right|_{{{{ x}}_i}}} + {\left. {\dfrac{{\partial {{ u}}}}{{\partial x}}} \right|_{{{{ x}}_i}}}(x - {x_i}) + {\left. {\dfrac{{\partial {{ u}}}}{{\partial y}}} \right|_{{{{ x}}_i}}}(y - {y_i})+ \\&&\qquad \dfrac{1}{2}{\left. {\dfrac{{{\partial ^2}{{ u}}}}{{\partial {x^2}}}} \right|_{{{{ x}}_i}}}{\left( {x - {x_i}} \right)^2} + {\left. {\dfrac{{{\partial ^2}{{ u}}}}{{\partial x\partial y}}} \right|_{{{{ x}}_i}}}(x - {x_i})(y - {y_i})+ \\&&\qquad \dfrac{1}{2}{\left. {\frac{{{\partial ^2}{{ u}}}}{{\partial {y^2}}}} \right|_{{{{ x}}_i}}}{\left( {y - {y_i}} \right)^2} + \cdots \end{eqnarray}$
将方程(8)等式两边同时在控制体单元$i$上积分可得

(9) $ \begin{eqnarray}\label{eq9} && \dfrac{1}{{{A_i}}}\int_{{V_i}} {{R_i}({{ x}} - {{{ x}}_i}){\rm d}A} \equiv {{{\bar{ u}}}_i} = {\left. {{ u}} \right|_{{{{ x}}_i}}} + {\left. {\dfrac{{\partial {{ u}}}}{{\partial x}}} \right|_{{{{ x}}_i}}}{{\bar x}_i} + {\left. {\dfrac{{\partial {{ u}}}}{{\partial y}}} \right|_{{{{ x}}_i}}}{{\bar y}_i}+ \\&&\qquad {\left. {\dfrac{{{\partial ^2}{{ u}}}}{{\partial {x^2}}}} \right|_{{{{ x}}_i}}}{\overline {{x^2}} _i} + {\left. {\dfrac{{{\partial ^2}{{ u}}}}{{\partial x\partial y}}} \right|_{{{{ x}}_i}}}{\overline {xy} _i} + {\left. {\dfrac{{{\partial ^2}{{ u}}}}{{\partial {y^2}}}} \right|_{{{{ x}}_i}}}{\overline {{y^2}} _i} + \cdots \end{eqnarray}$
式中, $\overline {{x^n}{y^m}} $为几何量在控制体上的体积分, Ollivier-Gooch在文献^[32] 中将其命名为Moments, 该积分量的具体表达式为

(10) $ \begin{eqnarray}\label{eq10} \overline {{x^n}{y^m}}_i = \frac{1}{{{A_i}}}\int_{C{V_i}} {{{\left( {x - {x_i}} \right)}^n}{{\left( {y - {y_i}} \right)}^m}{\rm d}A} \end{eqnarray}$
在此基础上, 将重构函数${R_i}({{ x}} - {{{ x}}_i})$在单元$j$ (即中心单元$i$的模板单元上积分可得

(11) $ \begin{eqnarray}\label{eq11} && \dfrac{1}{{{A_j}}}\int_{{V_j}} {{R_i}({{ x}} - {{{ x}}_i})} {\rm d}A = \\&&\qquad{\left. {{ u}} \right|_{{{{ x}}_i}}} + {\left. {\dfrac{{\partial {{ u}}}}{{\partial x}}} \right|_{{{{ x}}_i}}}\left[ {\dfrac{1}{{{A_j}}}\int_{{V_j}} {(x - {x_i})} {\rm d}A} \right] + \\&&\qquad {\left. {\dfrac{{\partial {{ u}}}}{{\partial y}}} \right|_{{{{ x}}_i}}}\left[ {\dfrac{1}{{{A_j}}}\int_{{V_j}} {(y - {y_i})} {\rm d}A} \right]+ \\&&\qquad {\left. {\dfrac{{{\partial ^2}{{ u}}}}{{\partial {x^2}}}} \right|_{{{{ x}}_i}}}\left[ {\dfrac{1}{{{A_j}}}\int_{{V_j}} {\dfrac{{{{(x - {x_i})}^2}}}{2}} {\rm d}A} \right] + \\&&\qquad {\left. {\dfrac{{{\partial ^2}{{ u}}}}{{\partial x\partial y}}} \right|_{{{{ x}}_i}}}\left[ {\dfrac{1}{{{A_j}}}\int_{{V_j}} {(x - {x_i})(y - {y_i})} {\rm d}A} \right]+ \\&&\qquad {\left. {\dfrac{{{\partial ^2}{{ u}}}}{{\partial {y^2}}}} \right|_{{{{ x}}_i}}}\left[ {\dfrac{1}{{{A_j}}}\int_{{V_j}} {\dfrac{{{{(y - {y_i})}^2}}}{2}} {\rm d}A} \right] + \cdots \end{eqnarray}$
经简化, 方程(11)最终可表示为

(12) $ \begin{eqnarray}\label{eq12} && \dfrac{1}{{{A_j}}}\int_{{V_j}} {{R_j}({{ x}} - {{{ x}}_i})} {\rm d}A = {\left. {{ u}} \right|_{{{{ x}}_i}}} + {\left. {\dfrac{{\partial {{ u}}}}{{\partial x}}} \right|_{{{{ x}}_i}}}{{\hat x}_{ij}} + {\left. {\dfrac{{\partial {{ u}}}}{{\partial y}}} \right|_{{{{ x}}_i}}}{{\hat y}_{ij}} + \\&&\qquad \dfrac{1}{2}{\left. {\dfrac{{{\partial ^2}{{ u}}}}{{\partial {x^2}}}} \right|_{{{{ x}}_i}}}{\widehat {{x^2}}_{ij}} + {\left. {\dfrac{{{\partial ^2}{{ u}}}}{{\partial x\partial y}}} \right|_{{{{ x}}_i}}}{\widehat {xy}_{ij}} + \\&&\qquad \dfrac{1}{2}{\left. {\dfrac{{{\partial ^2}{{ u}}}}{{\partial {y^2}}}} \right|_{{{{ x}}_i}}}{\widehat {{y^2}}_{ij}} + \cdots \end{eqnarray}$
式中${\widehat {{x^n}{y^m}}_{ij}}$同样代表几何量, 其计算表达式为

(13) $ \begin{eqnarray}\label{eq13} && \widehat{{x^n}{y^m}}_{ij} \equiv \dfrac{1}{{{A_j}}}\int_{{V_j}} \left\{ \left[ \left( {x - {x_j}} \right) + \left( {{x_j} - {x_i}} \right) \right]^n \right. \left[ \left( {y - {y_j}} \right) \right. + \\&&\qquad \left. \left. \left( {{y_j} - {y_i}} \right) \right]^m \right\}{\rm d}A = \sum\limits_{l = 0}^m \sum\limits_{k = 0}^n \dfrac{{m!}}{{l!(m - l)!}}\dfrac{{n!}}{{k!(n - k)!}}\cdot \\&&\qquad {{({x_j} - {x_i})}^k}{{({y_j} - {y_i})}^l}{{\overline {{x^{n - k}}{y^{m - l}}} }_j} \end{eqnarray}$
参照方程(12), 对于任意模板单元$j$, 均可构造一个重构多项式, 现将所有多项式整合可得到如方程(14)所示地矩阵方程组, 通过求解该方程组, 便可得到单元$i$的变量梯度及其高阶导数

(14) $ \left( \begin{matrix} 1 & {\bar x} & {\bar y} & {\overline {{x^2}} } & {\overline {xy} } & {\overline {{y^2}} } & \cdots \\ {{\omega_{i1}}}&{{\omega_{i1}}{{\widehat x}_{i1}}}&{{\omega_{i1}}{{\widehat y}_{i1}}}&{{\omega_{i1}}{{\widehat {{x^2}}}_{i1}}}&{{\omega_{i1}}{{\widehat {xy}}_{i1}}}&{{\omega_{i1}}{{\widehat {{y^2}}}_{i1}}}& \cdots \\ {{\omega_{i2}}}&{{\omega_{i2}}{{\widehat x}_{i2}}}&{{\omega_{i2}}{{\widehat y}_{i2}}}&{{\omega_{i2}}{{\widehat {{x^2}}}_{i2}}}&{{\omega_{i2}}{{\widehat {xy}}_{i2}}}&{{\omega_{i2}}{{\widehat {{y^2}}}_{i2}}}& \cdots \\ {{\omega_{i3}}}&{{\omega_{i3}}{{\widehat x}_{i3}}}&{{\omega_{i3}}{{\widehat y}_{i3}}}&{{\omega_{i3}}{{\widehat {{x^2}}}_{i3}}}&{{\omega_{i3}}{{\widehat {xy}}_{i3}}}&{{\omega_{i3}}{{\widehat {{y^2}}}_{i3}}}& \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \\ {{\omega_{iN}}}&{{\omega_{iN}}{{\widehat x}_{iN}}}&{{\omega_{iN}}{{\widehat y}_{iN}}}&{{\omega_{iN}}{{\widehat {{x^2}}}_{iN}}}&{{\omega_{iN}}{{\widehat {xy}}_{iN}}}&{{\omega_{iN}}{{\widehat {{y^2}}}_{iN}}}& \cdots \end{matrix} \right) \cdot \\ \left( \begin{matrix} { u} \\ \frac{{\partial \boldsymbol{u}}}{{\partial x}} \\ \frac{{\partial \boldsymbol{u}}}{{\partial y}} \\ {\frac{1}{2}\frac{{{\partial ^2}\boldsymbol{u}}}{{\partial {x^2}}}} \\ \vdots\\ {\frac{1}{2}\dfrac{{{\partial ^2}\boldsymbol{u}}}{{\partial {y^2}}}}\\ \end{matrix} \right) = \left( \begin{matrix} {{{\boldsymbol{\bar u}}_i}} \\ {{{\omega_{i1}\boldsymbol{\bar u}}_1}} \\ {{{\omega_{i2}\boldsymbol{\bar u}}_2}} \\ {{{\omega_{31}\boldsymbol{\bar u}}_3}} \\ \vdots \\ {{{\omega_{iN}\boldsymbol{\bar u}}_N}} \end{matrix} \right) $
式中每一个重构方程前面都乘以权系数$\omega_{ij}$

(15) $ \begin{eqnarray} {\omega_{ij}} = \dfrac{1}{{\left| {{{{ x}}_j} - {{{ x}}_i}} \right|}} \label{eq15} \end{eqnarray}$
该系数采用参考点之间距离倒数的形式, 用来强调距离中心单元较近的模板单元对重构结果的影响.

2 模板选择方法

首先介绍了两种常用的模板选择方法, 其次讨论了基于局部方向的模板选择方法, 并在大压缩比三角形网格上分析了这种模板存在的问题. 最后, 在局部方向模板缺陷的基础上, 给出了全局方向模板在不同压缩比三角形网格上的表现效果, 以进一步突显全局方向模板的优势.

2.1 两种常用的模板选择方法

较常用的模板选择方式有两种, 分别为共面模板与共点模板. 顾名思义, 共面模板由与中心单元共单元面的网格单元组成, 共点模板则由所有与中心单元共顶点的网格单元组成. 如图2所示, 两种方法的模板单元数量由模板层数控制. 如对于共点模板, 其第一层模板是所有与中心单元共顶点的网格单元, 而第二层模板则是所有与第一层模板共顶点的网格单元.

图2

新窗口打开|下载原图ZIP|生成PPT
图2共面与共点模板

Fig.2Face-neighbor and vertex-neighbor stencils



如图所示, 共点与共面模板均依赖于固定的网格拓扑关系, 针对不同流动, 这两种模板选择方式无法较好反映流场特性, 例如边界层型流动所表现出的流动各向异性.

2.2 局部方向模板选择方法

针对上述两种模板选择方式存在的问题, Xiong等^[35]在2018年提出了一种具有分维特征的结构化模板构造方法. 所谓结构化^[34], 是指模板选择沿着两个相互垂直的方向进行. 其次在确定局部方向的过程中, 能够尽量保证找到的两个局部方向沿着壁面法向与流向, 以此来捕捉流动的变化.

对于局部方向的确定, 首先要将网格单元分为四边形与三角形两种情况来考虑. 如图3所示, 如果网格单元是四边形, 两个局部方向则为两组对边中点的连线; 但三角形单元由于不存在两组对边, 因此Xiong等^[35]将三角形的其中一个顶点视作由四边形的一条边退化而成, 并将这个点称为三角形单元的退化点. 进而三角形网格单元的局部方向由退化点与对边中点的连线(其被命名为第一局部方向), 以及另外两条边中点的连线(第二局部方向)组成. 从图3可以看出, 对于一个三角形单元而言, 其退化点的位置有三种不同的可能性, 因此三角形网格单元的局部方向也存在三种可能的组合方式.

图3

新窗口打开|下载原图ZIP|生成PPT
图3三角形单元局部方向的可能组合方式

Fig.3Possible combinations of local directions on triangular cell



采用阵面推进技术^[36-37]可最终确定退化点的位置以及局部方向的最终组合方式, 由于其过程繁琐, 在此不再阐述. 经过阵面推进后, 可以找到所有网格单元的局部方向组合, 如图4所示. 图中黑色线段代表三角形网格单元及壁面, 红色线段表示单元内的局部方向.

图4

新窗口打开|下载原图ZIP|生成PPT
图4三角形网格单元的局部方向

Fig.4Local directions of triangular cell



在确定所有单元局部方向的基础上, 可沿着两个局部方向来扩充模板单元. 但对于三角形单元, 沿着局部方向确定模板单元还需要分为按点寻找与按边寻找两类, 该过程涉及局部方向的传递与夹角的判断, 由于较为繁琐, 本文不再重述其详细判断规则, 具体过程见文献^[35].

图5展示的是在不同网格压缩比的三角形网格上沿局部方向找到的模板单元效果. 如图5(a) 所示, 当网格压缩较小时, 第一局部方向接近壁面法向, 沿着两个局部方向找到的模板单元体现出结构化分维特征; 此外, 模板的空间延展性更好, 能够较为准确地反映流场变化. 但该方法仍存在一些问题, 如图5(b)所示, 当网格压缩比较大时, 模板呈现"畸形"现象.

图5

新窗口打开|下载原图ZIP|生成PPT
图5在不同压缩比三角形网格上的局部方向模板

Fig.5Local-direction stencil cells on triangular grids with different aspect ratios



出现这一问题的原因在于针对高度各项异性的大压缩比三角形网格, 基于"阵面推进"方法得到的局部方向严重偏离壁面法向, 并且两个局部方向之间的夹角较小, 导致模板在空间的延展性较差, 无法正确反映局部流动变化.

鉴于上述局部方向模板存在的问题, 迫切需要发展一种新的模板选择方式, 并具备能脱离网格拓扑约束与准确反映流动变化特征的能力.

2.3 全局方向模板选择方法

2.1与2.2节回顾了两种常用的模板选择方法以及局部方向模板的基本思路, 讨论了其在不同网格上的表现效果. 经分析, 局部方向模板虽体现了结构化特征, 在压缩比较小的三角形网格上能够较为准确地反映流场变化, 但随着网格压缩比的增大, 其第一局部方向逐渐偏离壁面法向, 导致方向性的失效. 据此, 在之前的工作中, 我们设计了一种基于全局方向的模板选择方式, 并将其应用于二阶精度非结构有限体积方法. 所谓全局方向, 是指对于边界层型各向异性流动, 首先可以确定两个特征方向, 其中一个方向沿着壁面法向, 也就是流场变化剧烈的方向, 通常沿该方向的网格尺度较小; 另外一个方向沿着流向, 对应的网格尺度较大. 所有的模板单元均沿着两个全局方向依次扩充, 因此避免了局部方向模板选择方法中的阵面推进过程, 并且具备准确反映流场特征的能力.

具体而言, 首先对于控制体单元, 过单元中心作两条平行于壁面法向与切向的直线, 其次在给定的单元集合中找到与两条直线相交的网格单元来构建新的模板. 通常, 该单元集合由与中心单元共点的网格单元构成. 图6所反映的是全局方向模板在不同压缩比网格上的效果, 图中红色与黑色虚线分别代表全局方向以及网格边界面法向, 绿色、蓝色、橙色以及黄色的网格单元分别代表中心单元、第一层与第二层共点单元, 以及全局方向模板单元. 如图6所示, 所有的模板单元均沿着两个正交的方向, 更加直观地体现出结构化分维特征; 其次, 不论网格是否存在大压缩比, 全局方向模板始终沿着壁面法向与切向. 因此, 相比局部方向模板, 全局方向模板可更加准确地准确捕捉到流动的各向异性特征.

图6

新窗口打开|下载原图ZIP|生成PPT
图6在不同压缩比三角形网格下的全局方向模板

Fig.6Global-direction stencil cells on triangular grids with different aspect ratios



根据图中所反映的现象可以看出, 基于该方法确定的模板单元能够更加准确地反映流场特征, 并且随着共点单元层数的增加, 模板单元沿两个全局方向依次增加, 因此模板数量可得到较为准确的控制. 此外, 经算例验证, 采用全局方向模板可有效提高二阶精度非结构有限体积方法的计算准确性. 为了进一步推广全局方向模板在高精度方法中的应用, 在下一节我们将通过典型的数值算例来初步检验该模板对变量梯度及高阶导数重构的影响.

3 梯度和高阶导数重构精度测试

为了检验全局方向模板在高阶精度非结构有限体积方法中的数值表现, 本节在不同压缩比的三角形网格上设置了一个典型的数值算例^[33]来检验采用该模板得到的变量梯度与高阶导数的重构准确性, 并将其与局部方向模板以及共点模板的重构效果对照, 进而验证新模板选择方法的有效性, 及其在高精度方法中应用并推广的可行性. 待测函数的形式如下

(16) $ \begin{eqnarray} u(x,y) = {(2 + \alpha x + \beta y)^4} \label{eq16} \end{eqnarray} $
该函数连续可导, 且给定的求解域范围是$x,y \in [ - 0.5,0.5] \times [ - 5{{\rm e}^{ - AR}},5{{\rm e}^{ - AR}}]$, 其中$AR$表示网格的压缩比, $\alpha$与$\beta$为两个实数, 其之间的关系为: $\beta / \alpha = AR$. 在网格压缩比较大时, 此函数沿着$y$方向(对于直线边界网格即壁面法向)的变化剧烈, 而在大压缩比情况下, 沿着$x$方向变化较小, 这种变化特征与边界层型流动所表现出的现象相似.

此外, 采用的计算网格为由背景四边形剖分得到的三角形网格. 为了体现非结构网格的随机特性, 并降低计算对网格的要求, 如图7所示, 我们分别测试了规则剖分与随机剖分两种网格拓扑, 并对两种拓扑下的规则网格添加随机扰动, 以充分对比不同模板的计算效果差异.

图7

新窗口打开|下载原图ZIP|生成PPT
图7计算网格示意图

Fig.7Schematic diagram of computational grids



为更加直观地体现在高度各向异性网格下, 不同模板对梯度与高阶导数重构的影响, 文章中给出了压缩比$AR = 10^3$时由4种网格分别得到的计算结果, 并在每种网格拓扑下, 设置了由疏到密的5套网格. 5套疏密程度不同的背景四边形网格在$x$与$y$方向网格量的分布如表1 所示.

Table 1
表1
表1背景四边形网格的网格量分布
Table 1The distribution of background quadrilateral grids

新窗口打开|下载CSV

对于高阶精度非结构有限体积方法, 重构过程的待求量除变量梯度外, 还有其高阶导数. 本算例测试了在三阶精度下的重构效果, 并采用变量梯度以及高阶导数的$L_2$误差(平均误差)来进行误差收敛性分析. 算例的解析解已知, 通过对函数(16)求导可得

(17) $ \begin{eqnarray}\label{eq17} \left.\begin{array}{l} \dfrac{{\partial u}}{{\partial x}} = 4\alpha {(2 + \alpha x + \beta y)^3} \hfill \\ \dfrac{{\partial u}}{{\partial y}} = 4\beta {(2 + \alpha x + \beta y)^3} \hfill \\ \dfrac{{{\partial ^2}u}}{{\partial {x^2}}} = 12{\alpha ^2}{\left( {2 + \alpha x + \beta y} \right)^2} \hfill \\ \dfrac{{{\partial ^2}u}}{{\partial x\partial y}} = 12\alpha \beta {\left( {2 + \alpha x + \beta y} \right)^2} \hfill \\ \dfrac{{{\partial ^2}u}}{{\partial {y^2}}} = 12{\beta ^2}{\left( {2 + \alpha x + \beta y} \right)^2} \hfill \\ \end{array}\right\} \end{eqnarray} $
考虑到重构单元梯度与高阶导数的最终目的是计算单元面高斯积分点处的变量点值与导数值, 以构造对流通量与黏性通量. 因此, 本节在测试单元梯度与高阶导数的基础上, 为了进一步验证不同模板的重构结果对实际计算的影响, 还分别统计了在不同网格上, 单元面高斯积分点处的变量点值与导数误差.

计算变量点值与导数时, 由单元中心插值到高斯积分点即可, 如公式(18)所示

(18) $ \begin{eqnarray}\label{eq18} \left. \begin{array}{l} {u_{{G_k}}} = {u_i} + {\left. {\dfrac{{\partial u}}{{\partial x}}} \right|_i}\left( {{x_i} - {x_{{G_k}}}} \right) + {\left. {\dfrac{{\partial u}}{{\partial y}}} \right|_i}\left( {{y_i} - {y_{{G_k}}}} \right)+\\ \qquad \dfrac{1}{2}{\left. {\dfrac{{{\partial ^2}u}}{{\partial {x^2}}}} \right|_i}{\left( {{x_i} - {x_{{G_k}}}} \right)^2} +\\ \qquad {\left. {\dfrac{{{\partial ^2}u}}{{\partial x\partial y}}} \right|_i}\left( {{x_i} - {x_{{G_k}}}} \right)\left( {{y_i} - {y_{{G_k}}}} \right) + \\ \qquad \dfrac{1}{2}{\left. {\dfrac{{{\partial ^2}u}}{{\partial {y^2}}}} \right|_i}{\left( {{y_i} - {y_{{G_k}}}} \right)^2} + \cdots \\ {\left. {\dfrac{{\partial u}}{{\partial x}}} \right|_{{G_k}}} = \dfrac{{\partial \left( {\dfrac{{\partial u}}{{\partial x}}} \right)}}{{\partial x}}\left( {{x_i} - {x_{{G_k}}}} \right) + \dfrac{{\partial \left( {\dfrac{{\partial u}}{{\partial y}}} \right)}}{{\partial x}}\left( {{y_i} - {y_{{G_k}}}} \right)\\ {\left. {\dfrac{{\partial u}}{{\partial y}}} \right|_{{G_k}}} = \dfrac{{\partial \left( {\dfrac{{\partial u}}{{\partial y}}} \right)}}{{\partial x}}\left( {{x_i} - {x_{{G_k}}}} \right) + \dfrac{{\partial \left( {\dfrac{{\partial u}}{{\partial y}}} \right)}}{{\partial y}}\left( {{y_i} - {y_{{G_k}}}} \right) \\ \end{array} \right\} \end{eqnarray}$
式中, 下标$G_k$与下标$i$分别代表高斯点与单元中心.

3.1 规则网格下的数值计算结果

图8和图9给出的计算结果分别代表在规则网格上, 使用3种模板得到的两个一阶导数${{\partial u}}/{{\partial x}}$, ${{\partial u}}/{{\partial y}}$, 以及高阶导数${{{\partial ^2}u}}/{{\partial {x^2}}}$, ${{{\partial ^2}u}}/{{\partial x\partial y}}$和${{{\partial ^2}u}}/{{\partial {y^2}}}$的$L_2$误差. 如果该函数的解具有三阶精度, 则其一阶导数(即梯度)应具备二阶精度, 二阶导数应具有一阶精度. 同时, 在图表中为简化对不同模板的描述, 使用V-Stencil、 L-Stencil以及G-Stencil来分别表示共点模板、局部方向模板以及全局方向模板, Slope2与Slope1代表二阶与一阶精度标准线. 此外, 表2所列举的数据为: 在规则密网格上, 由3种模板得到的梯度以及高阶导数$L_2$误差的具体数值.

图8

新窗口打开|下载原图ZIP|生成PPT
图8规则网格梯度$L_2$误差

Fig.8$L_2$ norm of gradient errors on regular grids

图9

新窗口打开|下载原图ZIP|生成PPT
图9规则网格上高阶导数的$L_2$误差

Fig.9$L_2$ norm of high-order derivative errors on regular grids



Table 2
表2
表2密网格上梯度与高阶导数的$L_2$误差
Table 2$L_2$ errors of derivatives on the vfin grid

新窗口打开|下载CSV

从图8反映的结果可以看出, 3种模板选择方法下梯度的计算结果均达到了设计精度, 并且采用全局方向模板得到的${{\partial u}}/{{\partial x}}$与${{\partial u}}/{{\partial y}}$误差与共点模板接近. 然而局部方向模板的数值表现较差, 采用该模板选择方法得到的计算误差在3种方法中最高. 此外, 我们统计了高阶导数的重构误差, 其结果如下.

从图9可以看出, 局部方向模板的计算误差在3种方法中偏高, ${{{\partial ^2}u}}/{{\partial x\partial y}}$和${{{\partial ^2}u}}/{{\partial {y^2}}}$两项与常用的共点模板差距较大. 此外, 结合表2中的数据可以看出, 全局方向模板的${{{\partial ^2}u}}/{{\partial {y^2}}}$误差与共点模板接近, 但稍高于共点模板, 而${{{\partial ^2}u}}/{{\partial {x^2}}}$与${{{\partial ^2}u}}/{{\partial x\partial y}}$误差在3种模板中最低. 因此, 当采用规则剖分的三角形网格时, 全局方向模板在梯度与高阶导数重构方面具有较好的数值表现.

此外, 为了进一步检验重构结果对实际计算的影响, 我们统计了高斯积分点处的变量导数与点值误差. 统计结果如图10、图11以及表3所示.

图10

新窗口打开|下载原图ZIP|生成PPT
图10高斯积分点处梯度$L_2$误差

Fig.10$L_2$ norm of gradient errors at Gauss point

图11

新窗口打开|下载原图ZIP|生成PPT
图11高斯积分点变量点值$L_2$误差

Fig.11$L_2$ norm of variable errors at Gauss point



Table 3
表3
表3密网格上高斯点处的变量点值与导数的$L_2$ 误差
Table 3$L_2$ errors of variable and derivatives at Gauss points

新窗口打开|下载CSV

从图10与图11以及表3可以看出, 在高斯积分点处, 采用全局方向模板得到的变量导数与点值误差与共点模板接近, 但均低于局部方向模板. 因此, 全局方向模板对单元梯度与高阶导数重构的优势, 在高斯积分点处得到了较好地验证.

3.2 扰动网格下的数值计算结果

在规则三角形网格的基础上, 我们对其添加随机节点扰动, 以测试不同模板对随机扰动网格的适应性. 图12与图13分别列举了梯度与高阶导数的$L_2$误差.

图12

新窗口打开|下载原图ZIP|生成PPT
图12扰动网格梯度$L_2$误差

Fig.12$L_2$ norm of gradient errors on perturbed grids

图13

新窗口打开|下载原图ZIP|生成PPT
图13扰动网格上高阶导数的$L_2$误差

Fig.13$L_2$ norm of high-order derivative errors on perturbed grids



从图12的结果可以看出, 由全局方向模板得到的梯度误差明显低于局部方向模板, 并且与共点模板的计算结果接近. 此外, 高阶导数的重构结果如下.

从图13可以看出, 在扰动网格上高阶导数的误差变化规律与规则网格相似, 全局方向模板对三个高阶导数的重构均有较好的数值表现. 此外, 高斯积分点处的变量导数与点值误差统计结果如图14和图15所示.

图14

新窗口打开|下载原图ZIP|生成PPT
图14高斯积分点处梯度$L_2$误差

Fig.14$L_2$ norm of gradient errors at Gauss point

图15

新窗口打开|下载原图ZIP|生成PPT
图15高斯积分点变量点值$L_2$误差

Fig.15$L_2$ norm of variable errors at Gauss point



从图中可以看出, 当网格存在扰动时, 高斯积分点处的误差变化规律同样与规则网格接近. 此外, 结合对单元梯度与高阶导数的重构结果, 可以充分说明全局方向模板对扰动网格同样具有较好的适应性.

3.3 随机剖分扰动网格下的数值计算结果

从3.2节可以看出, 对网格添加随机节点扰动后, 误差的变化规律与规则网格相似. 为简化结果分析, 本节给出了对随机剖分网格添加节点扰动后的数值结果. 图16与图17分别给出了梯度与高阶导数的$L_2$误差. 表4中列出的数据为密网格上, 各阶导数的具体误差值.

图16

新窗口打开|下载原图ZIP|生成PPT
图16随机剖分扰动网格梯度$L_2$误差

Fig.16$L_2$ norm of gradient errors on randomly-diag perturbed grids

图17

新窗口打开|下载原图ZIP|生成PPT
图17随机剖分扰动网格上高阶导数的$L_2$误差

Fig.17$L_2$ norm of high-order derivative errors on randomly-diag perturbed grids



Table 4
表4
表4密网格上梯度与高阶导数的$L_2$误差
Table 4$L_2$ errors of derivatives on the vfin grid

新窗口打开|下载CSV

从图中反映的结果可以明显看出, 由全局方向模板得到的两个一阶导数均达到了二阶精度, 但由局部方向模板得到的${{\partial u}}/{{\partial x}}$项的计算精度低于上述两种模板, 并且其计算误差明显偏高. 在此基础上, 还统计了高阶导数的重构误差.

从图17可以看出, 由全局方向模板得到的${{{\partial ^2}u}}/{{\partial {y^2}}}$ 与共点模板接近, 并稍高于共点模板, 而${{{\partial ^2}u}}/{{\partial {x^2}}}$ 与${{{\partial ^2}u}}/{{\partial x\partial y}}$误差在3种模板中最低. 此外, 从2.3 节对全局方向模板的描述以及图6中可以看出, 相比共点模板, 全局方向模板可有效减少用于重构过程的模板数量. 因此, 该模板的使用不仅显著改善了高阶精度非结构有限体积方法的计算准确性, 而且可提高计算效率.

为检验在当前网格拓扑上, 不同模板的单元梯度与高阶导数重构结果对实际计算的影响, 我们同样统计了高斯积分点处变量导数与点值误差, 统计结果如图18和图19所示. 此外, 表5中列出的数据为3种模板在密网格上的具体误差值.

图18

新窗口打开|下载原图ZIP|生成PPT
图18高斯积分点处梯度$L_2$误差

Fig.18$L_2$ norm of gradient errors at Gauss point

图19

新窗口打开|下载原图ZIP|生成PPT
图19高斯积分点变量点值$L_2$误差

Fig.19$L_2$ norm of variable errors at Gauss point



Table 5
表5
表5密网格上高斯点处的变量点值与导数的$L_2$ 误差
Table 5$L_2$ errors of variable and derivatives at Gauss points

新窗口打开|下载CSV

从统计结果来看, 由全局方向模板得到的高斯点处的变量导数与点值误差与共点模板接近, 并稍低于共点模板. 因此, 全局方向模板对梯度与高阶导数重构的优势得到了较好地验证.

反观局部方向模板, 由2.2节的分析可知, 在高度各向异性的三角形网格上, 由于局部方向与壁面法向间存在较大偏差, 导致模板单元无法较好地反映流动的各向异性特征. 从重构结果与高斯积分点处的误差统计也可以看出, 局部方向模板的计算误差较高. 因此, 当采用高度各向异性的三角形网格时, 不宜使用这种模板选择方式.

4 结论

本文通过对比不同模板得到的变量梯度与高阶导数重构结果, 验证了全局方向模板在高阶精度非结构有限体积方法中应用的可行性. 首先从操作层面而言, 相比局部方向模板, 全局方向模选择方法的过程简便, 避免了阵面推进与局部方向传递过程, 显著降低了模板选择这个前处理过程的复杂程度; 相比常用的共点模板选择方法, 新方法只需要在确定共点单元的基础上, 找到其中与两个全局方向相交的网格单元即可, 这个过程并未引入任何复杂性.

其次从计算结果来看, 不论网格拓扑是规则剖分还是随机剖分, 采用全局方向模板得到的梯度与高阶导数误差明显低于局部方向模板, 同时在计算高阶导数时, 全局方向模板得到的$({{{\partial ^2}u}})/({{\partial {x^2}}})$与$({{{\partial ^2}u}})/({{\partial x\partial y}})$误差值相比常用的共点模板更低. 并且相比共点模板, 全局方向模板的使用可有效减少重构过程所需的模板单元数量, 在一定程度上可节约CPU时间, 提高计算效率.

此外, 为检验不同模板的重构结果对实际计算的影响, 本文还统计了高斯积分点处变量点值与一阶导数的误差. 经检验, 全局方向模板得到的变量点值及导数误差与共点模板接近, 并稍低于共点模板, 而局部方向模板的误差值相对偏高, 不利于对流通量与黏性通量的构造.

综上, 全局方向模板在梯度与高阶导数重构过程数值表现良好, 并且重构优势在高斯积分点处得到了较好地验证, 具有在高阶精度非结构有限体积求解器上应用并推广的可行性. 接下来的工作将从两个方面开展, 首先从应用层面, 将进一步测试其在高精度非结构有限体积方法中的数值表现. 其次从数值分析层面, 在不同类型的网格上对重构矩阵的条件数与奇异性分析, 以及相应的模板优化工作同样也是下一步研究的重点.

参考文献 View Option 原文顺序
文献年度倒序
文中引用次数倒序
被引期刊影响因子

[1] Mavriplis DJ. Unstructured grid techniques
Annual Review of Fluid Mechanics, 1997,29:473-547
DOIURL [本文引用: 1]

[2] 曾扬兵, 沈孟育, 王保国 等. N-S方程在非结构网格下的求解
力学学报, 1996,28(6):641-650
URL [本文引用: 1]
在Ｒｏｅ的矢通量差分分裂的基础上，吸收了ＮＮＤ格式的优点，提出了一种非结构网格下求解Ｅｕｌｅｒ方程和Ｎ－Ｓ方程的高分辨率高精度迎风格式．这种格式具有捕捉强激波和滑移线的良好性能．在时间方向上采用了显式和隐式两种解法．文中还给出了自适应技术．最后，成功地完成了ＧＡＭＭ超音速前台阶绕流、二维平板无粘激波反射、三维Ｈｏｂｓｏｎ叶栅流动、ＶＫＩ叶栅流动、Ｃ３Ｘ叶栅流动的数值模拟，得到了满意的结果
( Zeng Yangbing, Shen Mengyu, Wang Baoguo, et al. Solution of N-S equations on unstructured grid
Chinese Journal of Theoretical and Applied Mechanics, 1996,28(6):641-650 (in Chinese))
URL [本文引用: 1]
在Ｒｏｅ的矢通量差分分裂的基础上，吸收了ＮＮＤ格式的优点，提出了一种非结构网格下求解Ｅｕｌｅｒ方程和Ｎ－Ｓ方程的高分辨率高精度迎风格式．这种格式具有捕捉强激波和滑移线的良好性能．在时间方向上采用了显式和隐式两种解法．文中还给出了自适应技术．最后，成功地完成了ＧＡＭＭ超音速前台阶绕流、二维平板无粘激波反射、三维Ｈｏｂｓｏｎ叶栅流动、ＶＫＩ叶栅流动、Ｃ３Ｘ叶栅流动的数值模拟，得到了满意的结果

[3] 王年华, 李明, 张来平. 非结构网格二阶有限体积法中黏性通量离散格式精度分析与改进
力学学报, 2018,50(3):527-537
[本文引用: 1]

( Wang Nianhua, Li Ming, Zhang Laiping. Accuracy analysis and improvement of viscous flux schemes in unstructured second-order finite-volume discretization
Chinese Journal of Theoretical and Applied Mechanics, 2018,50(3):527-537 (in Chinese))
[本文引用: 1]

[4] Dufresne Y, Moureau V, Lartigue G, et al. A massively parallel CFD/DEM approach for reactive gas-solid flows in complex geometries using unstructured meshes
Computers & Fluids, 2020,198:104402
[本文引用: 1]

[5] 王年华, 常兴华, 马戎 等. HyperFlow软件非结构网格亚跨声速湍流模拟的验证与确认
力学学报, 2019,51(3):813-825
[本文引用: 1]

( Wang Nianhua, Chang Xinghua, Ma Rong, et al. Verification and validation of HyperFlow solver for subsonic and transonic turbulent flow simulations on unstructured/hybrid grids
Chinese Journal of Theoretical and Applied Mechanics, 2019,51(3):813-825 (in Chinese))
[本文引用: 1]

[6] Park MA, Balan A, Anderson WK, et al. Verification of unstructured grid adaptation components
// AIAA Scitech 2019 Forum, San Diego, 2019-1-7-11, 2019: 1723


[7] 邵帅, 李明, 王年华 等. 基于非结构/混合网格模拟黏性流的高阶精度DDG/FV 混合方法
力学学报, 2018,50(6):1470-1482
[本文引用: 1]

( Shao Shuai, Li Ming, Wang Nianhua, et al. High-order DDG/FV hybrid method for viscous flow simulation based on unstructured / hybrid grid
Chinese Journal of Theoretical and Applied Mechanics, 2018,50(6):1470-1482 (in Chinese))
[本文引用: 1]

[8] Diskin B, Thomas JL, Nielsen EJ, et al. Comparison of node-centered and cell-centered unstructured finite-volume discretizations: Viscous fluxes
AIAA Journal, 2010,48:1326-1338
[本文引用: 2]

[9] Diskin B, Thomas JL. Comparison of node-centered and cell-centered unstructured finite-volume discretizations: Inviscid fluxes
AIAA Journal, 2011,49:836-854
DOIURL [本文引用: 1]

[10] Jameson A, Mavriplis DJ. Finite volume solution of the two-dimensional Euler equations on a regular triangular mesh
AIAA Journal, 1986,24:611-618


[11] Caraeni D, Hill DC. Unstructured-grid third-order finite volume discretization using a multistep quadratic data-reconstruction method
AIAA Journal, 2010,48:808-817


[12] Schwoppe A, Diskin B. Accuracy of the cell-centered grid metric in the DLR TAU-Code
Notes on Numerical Fluid Mechanics and Multidisciplinary Design, 2013,121:429-437
[本文引用: 1]

[13] 贺立新, 张来平, 张涵信. 间断Galerkin有限元和有限体积
力学学报, 2007,23(1):15-22
[本文引用: 1]

( He Lixin, Zhang Laiping, Zhang Hanxin. Discontinuous Galerkin finite element and finite volume
Chinese Journal of Theoretical and Applied Mechanics, 2007,23(1):15-22 (in Chinese))
[本文引用: 1]

[14] Mavriplis DJ. Multigrid solution of the two-dimensional Euler equations on unstructured triangular meshes
AIAA Journal, 1988,26:824-831
[本文引用: 1]

[15] 赵辉, 张耀冰, 陈江涛 等. 非结构网格体心梯度求解方法的精度分析
空气动力学学报, 2019,37(5):844-854


( Zhao Hui, Zhang Yaobing, Chen Jiangtao, et al. Accuracy analysis of the method for solving the cell-centered gradient on unstructured grid
Acta Aerodynamica Sinica, 2019,37(5):844-854 (in Chinese))


[16] Shima E, Kitamura K, Haga T. Green-gauss/weighted-least-squares hybrid gradient reconstruction for arbitrary polyhedra unstructured grids
AIAA Journal, 2013,51:2740-2747
[本文引用: 1]

[17] 王亚辉, 刘伟, 袁礼 等. 求解二维Euler方程有限单元边插值的降维重构算法
气体物理, 2019,3:34-41
[本文引用: 1]

( Wang Yahui, Liu Wei, Yuan Li, et al. A dimension reduction reconstruction algorithm for solving two-dimensional Euler equations with finite element edge interpolation
Physics of Gases, 2019,3:34-41 (in Chinese))
[本文引用: 1]

[18] Moukalled F, Mangani L, Darwish M. The Finite Volume Method in Computational Fluid Dynamics
Berlin: Springer, 2016: 46-49
[本文引用: 1]

[19] 张帆. 非结构网格有限体积法的空间离散算法研究. [博士论文]
大连: 大连理工大学, 2017
[本文引用: 1]

( Zhang Fan. Research on spatial discretization algorithm of unstructured mesh finite volume method. [PhD Thesis]
Dalian: Dalian University of Technology, 2017 (in Chinese))
[本文引用: 1]

[20] 王年华, 张来平, 赵钟 等. 基于制造解的非结构二阶有限体积离散格式的精度测试与验证
力学学报, 2017,49(3):627-637
[本文引用: 1]

( Wang Nianhua, Zhang Laiping, Zhao Zhong, et al. Accuracy verification of unstructured second-order finite volume discretization schemes based on the method of manufactured solutions
Chinese Journal of Theoretical and Applied Mechanics, 2017,49(3):627-637 (in Chinese))
[本文引用: 1]

[21] Sokolova I, Bastisya MG, Hajibeygi H. Multiscale finite volume method for finite-volume-based simulation of poroelasticity
Journal of Computational Physics, 2019,379:309-324.


[22] White JA, Nishikawa H, Baurle RA. Weighted least-squares cell-average gradient construction methods for The VULCAN-CFD second-order accurate unstructured grid cell-centered finite-volume solver
// AIAA Scitech 2019 Forum, San Diego, 2019-1-7-11, 2019: 0127


[23] Nishikawa H, White JA. Face-and cell-averaged nodal-gradient approach to cell-centered finite-volume method on mixed grids
// AIAA Scitech 2020 Forum, Orlando, 2020-1-6-10, 2020: 1787


[24] Nishikawa H. Efficient gradient stencils for robust implicit finite-volume solver convergence on distorted grids
Journal of Computational Physics, 2019,386:486-501
[本文引用: 1]

[25] Barth TJ, Jaespersen DC. The design and application of upwind schemes on unstructured meshes
// 27th Aerospace Sciences Meeting, Nevada, 1989-1-9-12, 1989: 0366
[本文引用: 1]

[26] Barth TJ. A 3-D upwind Euler solver for unstructured meshes
// 10th Computational Fluid Dynamics Conference. Hawaii, 1991-1-24-27, 1991: 1548
[本文引用: 1]

[27] Jalali A, Ollivier-Gooch CF. Higher-order unstructured finite volume RANS solution of turbulent compressible flows
Computers & Fluids, 2017,143:32-47
[本文引用: 1]

[28] Ollivier-Gooch CF. Analysis of unstructured meshes from GMGW-1/HiLiftPW-3//2018 AIAA Aerospace Sciences Meeting, Florida
2018-1-8-10, 2018: 0132


[29] Jalali A, Ollivier-Gooch CF. Higher-order unstructured finite volume methods for turbulent aerodynamic flows
// 22nd AIAA Computational Fluid Dynamics Conference, Dallas, 2015-1-22-26, 2015: 2284


[30] Jalali A, Ollivier-Gooch CF. Higher-order least-squares reconstruction for turbulent aerodynamic flows
// 11th World Congress on Computational Mechanics, Spain, 2014-7-20-25, 2014: 20-25
[本文引用: 1]

[31] Zangeneh R, Ollivier-Gooch CF. Stability analysis and improvement of the solution reconstruction for cell-centered finite volume methods on unstructured meshes
Journal of Computational Physics, 2019,393:375-405
[本文引用: 1]

[32] Ollivier-Gooch CF, Nejat A, Michalak K. Obtaining and verifying high-order unstructured finite volume solutions to the euler equations
AIAA Journal, 2009,47:2105-2120
DOIURL [本文引用: 1]

[33] Jalali A, Ollivier-Gooch CF. Higher-order finite volume solution reconstruction on highly anisotropic meshes
// 21st AIAA Computational Fluid Dynamics Conference, California, 2013-6-24-27, 2013: 2565
[本文引用: 2]

[34] Sozer E, Brehm C, Kiris CC. Gradient calculation methods on arbitrary polyhedral unstructured meshes for cell-centered cfd solvers
// 52nd Aerospace Science Meeting, Washington DC, 2014-1-13-17, 2014: 1440
[本文引用: 4]

[35] Xiong M, Deng XG, Gao X, et al. A novel stencil selection method for the gradient reconstruction on unstructured grid based on OpenFOAM
Computers & Fluids, 2018,172:426-442
[本文引用: 4]

[36] Mavriplis DJ. An advancing front Delaunay triangulation algorithm designed for robustness
Journal of Computational Physics, 1995,117:90-101
DOIURL [本文引用: 2]

[37] Luo SH. A new mesh generation scheme for arbitrary planar domains
International Journal for Numerical Methods in Engineering, 1985,21:1403-1426
DOIURL [本文引用: 2]

[38] Li XS, Ren XD, Gu CW, et al. Shock-stable roe scheme combining entropy fix and rotated Riemann solver
AIAA Journal, 2020,58:779-786
DOIURL [本文引用: 1]

[39] Roe PL. Characteristic-based schemes for the Euler equations
Annual Review of Fluid Mechanics, 1986,18:337-365
[本文引用: 1]