CN102880589B - 一种检测网格品质的有限差分模拟绕复杂构型流动的方法 - Google Patents

Abstract

本发明公开了一种检测网格品质的有限差分模拟绕复杂构型流动的方法，它包括：步骤一、计算网格的生成、检测和调整；步骤二、控制微分方程的离散和求解；步骤三、后置处理；其重点在于：对整个流动计算区域生成复杂多块对接结构网格；网格变换雅可比的数学定义式、对称形式和对称守恒形式采用线性有限差分算子离散后计算得到的网格变换雅可比分别记为V¹、V²和V³；根据所述V¹、V²和V³来确定网格品质的检测指标；根据所述网格品质的检测指标对生成的计算网格进行检测，确定出那些未达到检测指标的网格单元的具体位置；对所述的未达到检测指标的网格单元进行调整后重新检测，通过这一可重复的过程直到该调整后的网格单元全部达到检测指标。

Description

一种检测网格品质的有限差分模拟绕复杂构型流动的方法

技术领域

本发明涉及一种复杂计算区域物理特性的数值模拟方法，特别涉及一种检测绕复杂构型计算区域网格生成品质的有限差分模拟复杂流动的方法。

背景技术

尽管风洞试验到目前仍是预测各种航天航空飞行器气动特性的重要手段，但是数值模拟技术在气动设计中发挥着越来越重要的作用。近年来，随着计算机浮点运算能力的不断提升和数值计算方法的逐渐完善，人们越来越青睐于高精度、高分辨率的数值计算方法。现有大量的研究事实表明：基于传统二阶精度有限体积方法离散微分方程得到微分方程数值解的数值模拟方法并不能很好的满足实际工程问题对计算精度的需求，尤其是战斗机在作大攻角飞行时产生的大范围分离流动，二阶精度的有限体积方法往往不能给出令人满意的数值模拟结果，需要采用高精度、高分辨率的数值方法来进行模拟。战斗机在作大攻角飞行时产生的大范围分离流动是一种绕复杂构型流动。绕复杂构型流动是指空气或水等流体绕过飞机、导弹等各种真实航空飞行器复杂构型或水下飞行器的复杂流动。

当在多块对接结构网格中采用高阶精度有限差分方法求解时，如何保证数值模拟结果的高阶精度特性显得尤为重要。众所周知，要想得到高阶精度的数值解除了需要采用高阶精度的数值离散方法外，对网格生成质量的要求也不同于传统的二阶精度数值方法。需要特别指出的是，生成的网格高阶精度数值方法能够顺利进行计算和高阶精度数值方法能够得到较好的数值解是两个完全不同的概念。生成的网格高阶精度数值方法能够顺利进行计算是首先要满足的，在此基础上才有可能采用高阶精度数值方法得到较好的数值解。二者对网格生成质量的要求也不尽相同。显然，如果想要利用高阶精度数值方法得到较好的数值解，其对网格质量的要求要比仅仅是高阶精度数值方法能够顺利进行计算的要求高。那么，如何判定生成的网格是否能够在高阶精度数值方法的离散下得到较好的数值解呢？是否存在某种可以量化的网格质量的指标能够很好的指出网格的质量问题呢？此类问题鲜有国内外文献提及。事实上绕复杂构型的计算网格通常是由几百万个甚至数亿个的网格单元组成，很多情况下初始生成的计算网格不能够在高阶精度数值方法的离散下得到较好的数值解仅仅是由于极少数(如几个或者几十个)品质差的网格单元的存在造成的，而这些极少数品质差的网格单元又不太容易从大量的计算网格单元中被识别出来。本发明公开的网格品质检测方法能够有效地检测出计算网格的质量问题，确定出品质不达标网格单元的具体位置，通过简单地调整以保障网格品质，有利于模拟绕复杂构型流动时采用高阶精度数值方法顺利得到较好的数值解。

发明内容

本发明的目的是解决运用高阶精度有限差分方法模拟绕复杂构型流动时，如何保障网格生成品质以获得较好的或者是满足工程实用精度的数值解的问题。

为了达到上述目的，本发明采用的技术方案如下：

一种检测网格品质的有限差分模拟绕复杂构型流动的方法，它包括：

步骤一、计算网格的生成、检测和调整；

步骤二、控制微分方程的离散和求解；

步骤三、后置处理；

其特征在于在所述步骤一中1)、对整个流动计算区域生成复杂多块对接结构网格，2)、网格变换雅可比的数学定义式、对称形式和对称守恒形式采用线性有限差分算子离散后计算得到的网格变换雅可比分别记为V¹、V²和V³，3)、根据所述V¹、V²和V³来确定网格品质的检测指标V_N和V_M，所述V_N为所述V¹、V²和V³的最小值；所述V_M为V¹²、V²³和V¹³的最大值，其中V¹²、V²³和V¹³由所述V¹、V²和V³按V¹²＝∣V¹‐V²∣/V²、V²³＝∣V²‐V³∣/V³、V¹³＝∣V¹‐V³∣/V³确定，4)、根据所述网格品质的检测指标对生成的计算网格进行检测，确定出那些未达到检测指标的网格单元的具体位置，5)、对所述的未达到检测指标的网格单元进行调整后重新检测，通过这一可重复的过程直到该调整后的网格单元全部达到检测指标；在所述步骤二中对整个流动计算区域采用高阶精度有限差分方法离散控制微分方程，并对其进行求解；在所述步骤三中对数值模拟流场分析其对应的物理意义。

本发明具有如下技术效果：

1)由于网格检测只采用了网格变换雅可比一个变量，检测过程方便可行，检测标准明确易于判断；

2)由于网格检测会采用网格变换雅可比的不同计算形式，每种网格变换雅可比的计算形式涉及到的网格坐标点不完全相同，能够有效地对网格生成的品质进行较全面的检测；

3)由于采用本发明中的网格检测技术能够比较全面的检测网格生成的品质，在通过网格检测后，采用高阶精度的数值方法能够顺利地得到较好的高阶精度的数值解。

此外，本发明所述网格变换雅克比的三种计算形式均存在多种表现形式，这些多种表现形式在线性有限差分算子离散下等同。

附图说明

图1为具体实施方式的网格变换雅可比的计算网格示意图

图2为具体实施方式的网格变换雅可比的数学定义式(7)计算得到的体积示意图

图3为具体实施方式的网格变换雅可比的对称形式(9)计算得到的体积示意图(与ξ方向矢量面积相关的体积部分)

图4为具体实施方式的网格变换雅可比的对称形式(9)计算得到的体积示意图(与η方向矢量面积相关的体积部分)

图5为具体实施方式的网格变换雅可比的对称形式(9)计算得到的体积示意图(与ζ方向矢量面积相关的体积部分)

图6为具体实施方式的网格变换雅可比的守恒对称形式(16)计算得到的体积示意图(ξ方向左边的体积部分)

图7为具体实施方式的网格变换雅可比的守恒对称形式(16)计算得到的体积示意图(ξ方向右边的体积部分)

图8为具体实施方式的网格变换雅可比的守恒对称形式(16)计算得到的体积示意图(η方向左边的体积部分)

图9为具体实施方式的网格变换雅可比的守恒对称形式(16)计算得到的体积示意图(η方向右边的体积部分)

图10为具体实施方式的网格变换雅可比的守恒对称形式(16)计算得到的体积示意图(ζ方向左边的体积部分)

图11为具体实施方式的网格变换雅可比的守恒对称形式(16)计算得到的体积示意图(ζ方向右边的体积部分)

具体实施方式

结合附图以如下直角坐标系下的微分方程为例说明：

$\frac{\∂Q}{\∂t}+\frac{\∂E}{\∂x}+\frac{\∂F}{\∂y}+\frac{\∂G}{\∂z}=0---\left(1\right)$

其中Q为绕复杂构型流动中求解的物理变量，E、F和G均为关于Q的函数。方程(1)在多块对接结构网格中进行有限差分离散时，需将其变换到计算坐标系下，建立计算坐标系(τ,ξ,η,ζ)与直角坐标系(t,x,y,z)之间一一对应的变换关系为：

$\left\{\begin{array}{c}\mathrm{\τ}=t\\ \mathrm{\ξ}=\mathrm{\ξ}(t,x,y,z)\\ \mathrm{\η}=\mathrm{\η}(t,x,y,z)\\ \mathrm{\ζ}=\mathrm{\ζ}(t,x,y,z)\end{array}\right.---\left(2\right)$

则方程(1)在计算坐标系下的表现形式为：

$\frac{\∂\hat{Q}}{\∂\mathrm{\τ}}+\frac{\∂\hat{E}}{\∂\mathrm{\ξ}}+\frac{\∂\hat{F}}{\∂\mathrm{\η}}+\frac{\∂\hat{G}}{\∂\mathrm{\ζ}}=0---\left(3\right)$

其中：

$\begin{array}{c}\hat{Q}=J^{-1}Q\\ \hat{E}={\widehat{\mathrm{\ξ}}}_{t}Q+{\widehat{\mathrm{\ξ}}}_{x}E+{\widehat{\mathrm{\ξ}}}_{y}F+{\widehat{\mathrm{\ξ}}}_{z}G\\ \hat{F}={\widehat{\mathrm{\η}}}_{t}Q+{\widehat{\mathrm{\η}}}_{x}E+{\widehat{\mathrm{\η}}}_{y}F+{\widehat{\mathrm{\η}}}_{z}G\\ \hat{G}={\widehat{\mathrm{\ζ}}}_{t}Q+{\widehat{\mathrm{\ζ}}}_{x}E+{\widehat{\mathrm{\ζ}}}_{y}F+{\widehat{\mathrm{\ζ}}}_{z}G\end{array}---\left(4\right)$

静止网格中网格导数的数学定义式为：

$\begin{array}{c}\left\{\begin{array}{c}{\widehat{\mathrm{\ξ}}}_x=J^{-1}{\mathrm{\ξ}}_x=y_{\mathrm{\η}}{z}_{\mathrm{\ζ}}-y_{\mathrm{\ζ}}{z}_{\mathrm{\η}}\\ {\widehat{\mathrm{\ξ}}}_y=J^{-1}{\mathrm{\ξ}}_y=x_{\mathrm{\ζ}}{z}_{\mathrm{\η}}-x_{\mathrm{\η}}{z}_{\mathrm{\ζ}}\\ {\widehat{\mathrm{\ξ}}}_z=J^{-1}{\mathrm{\ξ}}_z=x_{\mathrm{\η}}{y}_{\mathrm{\ζ}}-x_{\mathrm{\ζ}}{y}_{\mathrm{\η}}\end{array}\right.\\ \left\{\begin{array}{c}{\widehat{\mathrm{\η}}}_x=J^{-1}{\mathrm{\η}}_x=y_{\mathrm{\ζ}}{z}_{\mathrm{\ξ}}-y_{\mathrm{\ξ}}{z}_{\mathrm{\ζ}}\\ {\widehat{\mathrm{\η}}}_y=J^{-1}{\mathrm{\η}}_y=x_{\mathrm{\ξ}}{z}_{\mathrm{\ζ}}-x_{\mathrm{\ζ}}{z}_{\mathrm{\ξ}}\\ {\widehat{\mathrm{\η}}}_z=J^{-1}{\mathrm{\η}}_z=x_{\mathrm{\ζ}}{y}_{\mathrm{\ξ}}-x_{\mathrm{\ξ}}{y}_{\mathrm{\ζ}}\end{array}\right.\\ \left\{\begin{array}{c}{\widehat{\mathrm{\ζ}}}_x=J^{-1}{\mathrm{\ζ}}_x=y_{\mathrm{\ξ}}{z}_{\mathrm{\η}}-y_{\mathrm{\η}}{z}_{\mathrm{\ξ}}\\ {\widehat{\mathrm{\ζ}}}_y=J^{-1}{\mathrm{\ζ}}_y=x_{\mathrm{\η}}{z}_{\mathrm{\ξ}}-x_{\mathrm{\ξ}}{z}_{\mathrm{\η}}\\ {\widehat{\mathrm{\ζ}}}_z=J^{-1}{\mathrm{\ζ}}_z=x_{\mathrm{\ξ}}{y}_{\mathrm{\η}}-x_{\mathrm{\η}}{y}_{\mathrm{\ξ}}\end{array}\right.\end{array}---\left(5\right)$

其中下标表示偏导数，如x_ξ表示坐标x对计算坐标ξ方向的偏导数。

为了能够严格满足几何守恒律，同时又能准确地反映计算网格的几何特性，需要采用网格导数的对称守恒计算形式：

$\begin{array}{c}\left\{\begin{array}{c}{\widehat{\mathrm{\ξ}}}_x=\frac{1}{2}\[{\left({\mathrm{zy}}_{\mathrm{\η}}\right)}_{\mathrm{\ζ}}+{\left({\mathrm{yz}}_{\mathrm{\ζ}}\right)}_{\mathrm{\η}}-{\left({\mathrm{zy}}_{\mathrm{\ζ}}\right)}_{\mathrm{\η}}-{\left({\mathrm{yz}}_{\mathrm{\η}}\right)}_{\mathrm{\ζ}}\]\\ {\widehat{\mathrm{\ξ}}}_y=\frac{1}{2}\[{\left({\mathrm{xz}}_{\mathrm{\η}}\right)}_{\mathrm{\ζ}}+{\left({\mathrm{zx}}_{\mathrm{\ζ}}\right)}_{\mathrm{\η}}-{\left({\mathrm{xz}}_{\mathrm{\ζ}}\right)}_{\mathrm{\η}}-{\left({\mathrm{zx}}_{\mathrm{\η}}\right)}_{\mathrm{\ζ}}\]\\ {\widehat{\mathrm{\ξ}}}_z=\frac{1}{2}\[{\left({\mathrm{yx}}_{\mathrm{\η}}\right)}_{\mathrm{\ζ}}+{\left({\mathrm{xy}}_{\mathrm{\ζ}}\right)}_{\mathrm{\η}}-{\left({\mathrm{yx}}_{\mathrm{\ζ}}\right)}_{\mathrm{\η}}-{\left({\mathrm{xy}}_{\mathrm{\η}}\right)}_{\mathrm{\ζ}}\]\end{array}\right.\\ \left\{\begin{array}{c}{\widehat{\mathrm{\η}}}_x=\frac{1}{2}\[{\left({\mathrm{zy}}_{\mathrm{\ζ}}\right)}_{\mathrm{\ξ}}+{\left({\mathrm{yz}}_{\mathrm{\ξ}}\right)}_{\mathrm{\ζ}}-{\left({\mathrm{zy}}_{\mathrm{\ξ}}\right)}_{\mathrm{\ζ}}-{\left({\mathrm{yz}}_{\mathrm{\ζ}}\right)}_{\mathrm{\ξ}}\]\\ {\widehat{\mathrm{\η}}}_y=\frac{1}{2}\[{\left({\mathrm{xz}}_{\mathrm{\ζ}}\right)}_{\mathrm{\ξ}}+{\left({\mathrm{zx}}_{\mathrm{\ξ}}\right)}_{\mathrm{\ζ}}-{\left({\mathrm{xz}}_{\mathrm{\ξ}}\right)}_{\mathrm{\ζ}}-{\left({\mathrm{zx}}_{\mathrm{\ζ}}\right)}_{\mathrm{\ξ}}\]\\ {\widehat{\mathrm{\η}}}_z=\frac{1}{2}\[{\left({\mathrm{yx}}_{\mathrm{\ζ}}\right)}_{\mathrm{\ξ}}+{\left({\mathrm{xy}}_{\mathrm{\ξ}}\right)}_{\mathrm{\ζ}}-{\left({\mathrm{yx}}_{\mathrm{\ξ}}\right)}_{\mathrm{\ζ}}-{\left({\mathrm{xy}}_{\mathrm{\ζ}}\right)}_{\mathrm{\ξ}}\]\end{array}\right.\\ \left\{\begin{array}{c}{\widehat{\mathrm{\ζ}}}_x=\frac{1}{2}\[{\left({\mathrm{zy}}_{\mathrm{\ξ}}\right)}_{\mathrm{\η}}+{\left({\mathrm{yz}}_{\mathrm{\η}}\right)}_{\mathrm{\ξ}}-{\left({\mathrm{zy}}_{\mathrm{\η}}\right)}_{\mathrm{\ξ}}-{\left({\mathrm{yz}}_{\mathrm{\ξ}}\right)}_{\mathrm{\η}}\]\\ {\widehat{\mathrm{\ζ}}}_y=\frac{1}{2}\[{\left({\mathrm{xz}}_{\mathrm{\ξ}}\right)}_{\mathrm{\η}}+{\left({\mathrm{xz}}_{\mathrm{\η}}\right)}_{\mathrm{\ξ}}-{\left({\mathrm{xz}}_{\mathrm{\η}}\right)}_{\mathrm{\ξ}}-{\left({\mathrm{zx}}_{\mathrm{\ξ}}\right)}_{\mathrm{\η}}\]\\ {\widehat{\mathrm{\ζ}}}_z=\frac{1}{2}\[{\left({\mathrm{yx}}_{\mathrm{\ξ}}\right)}_{\mathrm{\η}}+{\left({\mathrm{xy}}_{\mathrm{\η}}\right)}_{\mathrm{\ξ}}-{\left({\mathrm{yx}}_{\mathrm{\η}}\right)}_{\mathrm{\ξ}}-{\left({\mathrm{xy}}_{\mathrm{\ξ}}\right)}_{\mathrm{\η}}\]\end{array}\right.\end{array}---\left(6\right)$

以及网格变换雅可比的数学定义式为：

$J^{-1}=\left|\begin{array}{ccc}{x}_{\mathrm{\ξ}}& x_{\mathrm{\η}}& x_{\mathrm{\ζ}}\\ y_{\mathrm{\ξ}}& y_{\mathrm{\η}}& y_{\mathrm{\ζ}}\\ z_{\mathrm{\ξ}}& z_{\mathrm{\η}}& z_{\mathrm{\ζ}}\end{array}\right|=x_{\mathrm{\ξ}}{y}_{\mathrm{\η}}{z}_{\mathrm{\ζ}}-x_{\mathrm{\ξ}}{y}_{\mathrm{\ζ}}{z}_{\mathrm{\η}}+x_{\mathrm{\η}}{y}_{\mathrm{\ζ}}{z}_{\mathrm{\ξ}}-x_{\mathrm{\η}}{y}_{\mathrm{\ξ}}{z}_{\mathrm{\ζ}}+x_{\mathrm{\ζ}}{y}_{\mathrm{\ξ}}{z}_{\mathrm{\η}}-x_{\mathrm{\ζ}}{y}_{\mathrm{\η}}{z}_{\mathrm{\ξ}}---\left(7\right)$

网格变换雅可比的数学定义式(7)的几何本质是计算网格单元的体积。考虑到体积的几何定义，任意单元的体积由其面积围成，可由其面积唯一确定。将网格变换雅可比的数学定义式(7)表示为面积的形式：

$\begin{array}{c}{J}^{-1}=x_{\mathrm{\ξ}}{\widehat{\mathrm{\ξ}}}_x+x_{\mathrm{\η}}{\widehat{\mathrm{\η}}}_x+x_{\mathrm{\ζ}}{\widehat{\mathrm{\ζ}}}_x\stackrel{\mathrm{\Δ}}{=}{J}^{-1}\left(x\right)\\ J^{-1}=y_{\mathrm{\ξ}}{\widehat{\mathrm{\ξ}}}_y+y_{\mathrm{\η}}{\widehat{\mathrm{\η}}}_y+y_{\mathrm{\ζ}}{\widehat{\mathrm{\ζ}}}_y\stackrel{\mathrm{\Δ}}{=}{J}^{-1}\left(y\right)\\ J^{-1}=z_{\mathrm{\ξ}}{\widehat{\mathrm{\ξ}}}_z+z_{\mathrm{\η}}{\widehat{\mathrm{\η}}}_z+z_{\mathrm{\ζ}}{\widehat{\mathrm{\ζ}}}_z\stackrel{\mathrm{\Δ}}{=}{J}^{-1}\left(z\right)\end{array}---\left(8\right)$

式(8)在进行有限差分离散时并不具备体积的几何意义，将式(8)改写为对称形式：

$J^{-1}=\frac{1}{3}\left(x_{\mathrm{\ξ}}{\widehat{\mathrm{\ξ}}}_x+y_{\mathrm{\ξ}}{\widehat{\mathrm{\ξ}}}_y+z_{\mathrm{\ξ}}{\widehat{\mathrm{\ξ}}}_z+x_{\mathrm{\η}}{\widehat{\mathrm{\η}}}_x+y_{\mathrm{\η}}{\widehat{\mathrm{\η}}}_y+z_{\mathrm{\η}}{\widehat{\mathrm{\η}}}_z+x_{\mathrm{\ζ}}{\widehat{\mathrm{\ζ}}}_x+y_{\mathrm{\ζ}}{\widehat{\mathrm{\ζ}}}_y+z_{\mathrm{\ζ}}{\widehat{\mathrm{\ζ}}}_z\right)---\left(9\right)$

将式(8)改写为守恒形式：

$\begin{array}{c}{J}^{-1}\left(x\right)={\left(x{\widehat{\mathrm{\ξ}}}_x\right)}_{\mathrm{\ξ}}+{\left(x{\widehat{\mathrm{\η}}}_x\right)}_{\mathrm{\η}}+{\left(x{\widehat{\mathrm{\ζ}}}_x\right)}_{\mathrm{\ζ}}-{\mathrm{xI}}_x\\ J^{-1}\left(y\right)={\left(y{\widehat{\mathrm{\ξ}}}_y\right)}_{\mathrm{\ξ}}+{\left(y{\widehat{\mathrm{\η}}}_y\right)}_{\mathrm{\η}}+{\left(y{\widehat{\mathrm{\ζ}}}_y\right)}_{\mathrm{\ζ}}-{\mathrm{yI}}_y\\ J^{-1}\left(z\right)={\left(z{\widehat{\mathrm{\ξ}}}_z\right)}_{\mathrm{\ξ}}+{\left(z{\widehat{\mathrm{\η}}}_z\right)}_{\mathrm{\η}}+{\left(z{\widehat{\mathrm{\ζ}}}_z\right)}_{\mathrm{\ζ}}-{\mathrm{zI}}_z\end{array}---\left(10\right)$

其中I_x、I_y和I_z为几何守恒律余项：

$\begin{array}{c}{I}_x=\frac{\∂{\widehat{\mathrm{\ξ}}}_x}{\∂\mathrm{\ξ}}+\frac{\∂{\widehat{\mathrm{\η}}}_x}{\∂\mathrm{\η}}+\frac{\∂{\widehat{\mathrm{\ζ}}}_x}{\∂\mathrm{\ζ}}\\ I_y=\frac{\∂{\widehat{\mathrm{\ξ}}}_y}{\∂\mathrm{\ξ}}+\frac{\∂{\widehat{\mathrm{\η}}}_y}{\∂\mathrm{\η}}+\frac{\∂{\widehat{\mathrm{\ζ}}}_y}{\∂\mathrm{\ζ}}\\ I_z=\frac{\∂{\widehat{\mathrm{\ξ}}}_z}{\∂\mathrm{\ξ}}+\frac{\∂{\widehat{\mathrm{\η}}}_z}{\∂\mathrm{\η}}+\frac{\∂{\widehat{\mathrm{\ζ}}}_z}{\∂\mathrm{\ζ}}\end{array}---\left(11\right)$

在采用有限差分方法离散时，将网格导数的对称守恒形式(6)中的内层差分算子记为δ³，外层差分算子记为δ²，将式(11)中的差分算子记为δ¹，则网格变换雅可比的守恒计算形式(10)和(15)以及其对称守恒计算形式(16)中的差分算子均为δ¹。为了能够严格满足几何守恒律，δ¹、δ²和δ³需均为线性有限差分算子。所谓线性有限差分算子，是指差分算子δ具有如下的性质：

对于任意的常数a、b及变量φ、均有下式成立：

以及求导顺序可交换性质，即：

δ_ξ(δ_ηφ)＝δ_η(δ_ξφ)(13)

几何守恒律的严格满足要求所述线性有限差分算子满足如下的约束条件：

${\mathrm{\δ}}_{\mathrm{\ξ}}^1={\mathrm{\δ}}_{\mathrm{\ξ}}^2={\mathrm{\δ}}_{\mathrm{\ξ}}^3,{\mathrm{\δ}}_{\mathrm{\η}}^1={\mathrm{\δ}}_{\mathrm{\η}}^2={\mathrm{\δ}}_{\mathrm{\η}}^3,{\mathrm{\δ}}_{\mathrm{\ζ}}^1={\mathrm{\δ}}_{\mathrm{\ζ}}^2={\mathrm{\δ}}_{\mathrm{\ζ}}^3---\left(14\right)$

在几何守恒律严格满足的条件下：I_x＝I_y＝I_z＝0，将式(10)写为：

$\begin{array}{c}{J}^{-1}\left(x\right)={\left(x{\widehat{\mathrm{\ξ}}}_x\right)}_{\mathrm{\ξ}}+{\left(x{\widehat{\mathrm{\η}}}_x\right)}_{\mathrm{\η}}+{\left(x{\widehat{\mathrm{\ζ}}}_x\right)}_{\mathrm{\ζ}}\\ J^{-1}\left(y\right)={\left(y{\widehat{\mathrm{\ξ}}}_y\right)}_{\mathrm{\ξ}}+{\left(y{\widehat{\mathrm{\η}}}_y\right)}_{\mathrm{\η}}+{\left(y{\widehat{\mathrm{\ζ}}}_y\right)}_{\mathrm{\ζ}}\\ J^{-1}\left(z\right)={\left(z{\widehat{\mathrm{\ξ}}}_z\right)}_{\mathrm{\ξ}}+{\left(z{\widehat{\mathrm{\η}}}_z\right)}_{\mathrm{\η}}+{\left(z{\widehat{\mathrm{\ζ}}}_z\right)}_{\mathrm{\ζ}}\end{array}---\left(15\right)$

同理，式(15)在进行有限差分离散时并不具备体积的几何意义，将式(15)改写网格变换雅可比的对称守恒形式：

$J^{-1}=\frac{1}{3}\[{\left(x{\widehat{\mathrm{\ξ}}}_x+y{\widehat{\mathrm{\ξ}}}_y+z{\widehat{\mathrm{\ξ}}}_z\right)}_{\mathrm{\ξ}}+{\left(x{\widehat{\mathrm{\η}}}_x+y{\widehat{\mathrm{\η}}}_y+z{\widehat{\mathrm{\η}}}_z\right)}_{\mathrm{\η}}+{\left(x{\widehat{\mathrm{\ζ}}}_x+y{\widehat{\mathrm{\ζ}}}_y+z{\widehat{\mathrm{\ζ}}}_z\right)}_{\mathrm{\ζ}}\]---\left(16\right)$

网格变换雅可比的几何本质是网格单元的体积。在对网格变换雅可比不同的计算形式进行有限差分离散时，网格变换雅可比的数学定义式(7)、网格变换雅可比的对称形式(9)、网格变换雅可比的对称守恒形式(16)均能保持其网格单元体积的几何特性，而表达式(8)和守恒形式(15)均不能保持其应有的几何特性，这正是需要将式(8)与式(15)分别写改写其对称形式(9)和式(16)的原因。

当采用线性有限差分算子离散时，不妨记由网格变换雅可比的数学定义式(7)计算得到的网格变换雅可比为V¹；记由网格变换雅可比的对称形式(9)计算得到的网格变换雅可比为V²；记由网格变换雅可比的对称守恒形式(16)计算得到的网格变换雅可比为V³。如果计算采用的网格示意图如图1所示，则当采用二阶中心差分格式离散时，由网格变换雅可比的数学定义式(7)计算得到点O处的网格变换雅可比V¹所用网格点坐标信息的示意图如图2所示；由网格变换雅可比的对称形式(9)计算得到点O处的网格变换雅可比V²所用网格点坐标信息的示意图分别为图3、图4和图5中所示的三个部分体积相加得到；由网格变换雅可比的对称守恒形式(16)计算得到点O处的网格变换雅可比V³所用网格点坐标信息的示意图如图6、图7、图8、图9、图10和图11中所示的六个部分体积相加得到。特别值得注意的是，采用二阶中心差分格式离散后，点O处的网格变换雅可比V¹并不完全等于图2所示的体积，而是图2所示体积的四分之三；V²也不完全等于图3、图4和图5所示的三部分体积之和，而是图3、图4和图5所示的三部分体积之和的八分之一；同理V³也不完全等于图6至图11中所示的六部分体积之和，而是图6至图11中所示的六部分体积之和的八分之一。

容易看出，图6至图11中所有六部分体积之和刚好为图1中所示包含八个网格单元的大网格单元的体积。由于V¹、V²和V³在数学定义上是完全等价的，其计算得到的体积乘以八倍后均能在一定程度上代表图1中所示由八个网格单元组成的大网格单元的体积，其中V¹乘以八倍后相当于以图2中所示体积的六倍来代表图1所示的大网格单元的体积(V¹为图2中所示体积的四分之三)，而V²和V³在乘以八倍以后均可以直接表示图1所示的大网格单元的体积。

由图2至图11中所示的不同体积示意图可以看出，V¹、V²和V³计算体积时采用的网格点坐标信息不完全相同，其中V¹所采用的网格点坐标信息最少，V²其次，V³采用的网格点坐标信息最多，V³用到了图1中所示构成大网格单元的所有网格点的坐标信息。因此V¹、V²和V³之间的差异刚好可以体现网格点坐标之间相对位置的差异。举例说明：图2中所示体积的六倍(计算体积V¹的八倍)与图1中所示的大网格体积(计算体积V³的八倍)之间的差异正好反映了图1中组成大网格单元的八个网格单元与平行六面体之间的差异，亦即图1中封闭大网格单元的所有四边形与平行四边形的差异。换言之，如果图1中所示组成大网格单元的八个网格单元均为平行六面体，则图2中所示体积的六倍刚好与图1中所示大网格单元的体积相等。如果图1中所示组成大网格单元的八个网格单元均不是平行六面体或者不全是平行六面体，则图2中所示体积的六倍与图1中所示大网格单元的体积一般不相等。

同理，图3至图5中所示的体积V²所采用的网格点坐标信息比V¹多，但比V³少，其分别与V¹和V³的差别亦可以分别反映出图1中围成大网格单元的某些面与平行四边形的差异。如果图1中围成大网格单元的所有的四边形均为平行四边形，即计算所采用的网格品质很好(最理想的计算网格)，则V¹、V²和V³均相等。鉴于此，本发明设计如下的网格品质检测指标：

$\begin{array}{c}\begin{array}{cc}{V}_N=\min \left(V_{i,j,k}^1,V_{i,j,k}^2,V_{i,j,k}^3\right)& foralli,j,k\end{array}\\ \begin{array}{cc}{V}_M=\max \left(V_{i,j,k}^{12},V_{i,j,k}^{23},V_{i,j,k}^{13}\right)& foralli,j,k\end{array}\\ V^{12}=\frac{\left|V^1-V^2\right|}{V^2}\×100\%,V^{23}=\frac{\left|V^2-V^3\right|}{V^3}\×100\%,V^{13}=\frac{\left|V^1-V^3\right|}{V^3}\×100\%\end{array}---\left(17\right)$

即V_N为网格变换雅可比三种计算形式V¹、V²和V³在所有计算网格上的最小值；V_M为V¹²、V²³和V¹³各网格变换雅可比计算分量之间相对误差的最大值。值得注意的是，尽管上述对网格变换雅克比的分析都是基于二阶精度的中心差分格式而言的，对于高阶精度的有限差分格式同样适用，其对网格品质的检测同样能够应用于高阶精度的有限差分方法中。

对于任意的多块对接结构网格而言，首先要保证其体积非负，即V_N＞0。在采用高阶精度的数值方法离散计算域时如果想要得到较好的高阶精度数值解，则网格品质还需满足另一要求：V_M≤ε。事实上指标V_M反映了网格单元畸变的程度，是所采用的计算网格本身的属性，网格品质越好则其网格单元越规则，V_M值就会越小；指标ε反映了数值计算方法能够得到较好数值模拟结果所允许的网格单元畸变的程度，是所采用数值方法计算能力的体现，ε值越小表明数值计算方法对网格生成品质的要求越高。根据实际工程应用的经验，对于高阶精度(及二阶精度)数值计算方法我们给出ε＝10％。

通过检测确定出未达到检测指标网格单元的具体位置，对所述的未达到检测指标的网格单元进行调整后重新检测，通过这一可重复的过程直到调整后的网格单元全部达到检测指标；这样能够有效地保证网格品质满足高阶精度有限差分方法计算的要求。则在采用高阶精度有限差分方法离散求解时，能够顺利地得到较好的高阶精度数值解。对数值模拟结果进行后置处理能够实现对绕复杂构型计算域内物理变量Q的较为理想的高阶精度数值模拟。

网格变换雅可比的三种计算形式均有不同的表现形式，如将V¹、V²和V³的计算表达式可以分别改写为矢量的形式：

网格变换雅可比的数学定义式(7)的矢量形式：

$J^{-1}=\left({\stackrel{\→}{r}}_{\mathrm{\ξ}}\×{\stackrel{\→}{r}}_{\mathrm{\η}}\right)\·{\stackrel{\→}{r}}_{\mathrm{\ζ}}=\left({\stackrel{\→}{r}}_{\mathrm{\η}}\×{\stackrel{\→}{r}}_{\mathrm{\ζ}}\right)\·{\stackrel{\→}{r}}_{\mathrm{\ξ}}=\left({\stackrel{\→}{r}}_{\mathrm{\ζ}}\×{\stackrel{\→}{r}}_{\mathrm{\ξ}}\right)\·{\stackrel{\→}{r}}_{\mathrm{\η}}---\left(18\right)$

网格变换雅可比的对称形式(9)的矢量形式：

$J^{-1}=\frac{1}{3}\left({\stackrel{\→}{r}}_{\mathrm{\ξ}}\·{\stackrel{\→}{S}}^{\mathrm{\ξ}}+{\stackrel{\→}{r}}_{\mathrm{\η}}\·{\stackrel{\→}{S}}^{\mathrm{\η}}+{\stackrel{\→}{r}}_{\mathrm{\ζ}}\·{\stackrel{\→}{S}}^{\mathrm{\ζ}}\right)---\left(19\right)$

网格变换雅可比的对称守恒形式(16)的矢量形式：

$J^{-1}=\frac{1}{3}\[{\left(\stackrel{\→}{r}\·{\stackrel{\→}{S}}^{\mathrm{\ξ}}\right)}_{\mathrm{\ξ}}+{\left(\stackrel{\→}{r}\·{\stackrel{\→}{S}}^{\mathrm{\η}}\right)}_{\mathrm{\η}}+{\left(\stackrel{\→}{r}\·{\stackrel{\→}{S}}^{\mathrm{\ζ}}\right)}_{\mathrm{\ζ}}\]---\left(20\right)$

其中：

$\begin{array}{c}{\stackrel{\→}{S}}^{\mathrm{\ξ}}={\widehat{\mathrm{\ξ}}}_x\stackrel{\→}{i}+{\widehat{\mathrm{\ξ}}}_y\stackrel{\→}{j}+{\widehat{\mathrm{\ξ}}}_x\stackrel{\→}{k}\\ {\stackrel{\→}{S}}^{\mathrm{\η}}={\widehat{\mathrm{\η}}}_x\stackrel{\→}{i}+{\widehat{\mathrm{\η}}}_y\stackrel{\→}{j}+{\widehat{\mathrm{\η}}}_z\stackrel{\→}{k}\\ {\stackrel{\→}{S}}^{\mathrm{\ζ}}={\widehat{\mathrm{\ζ}}}_x\stackrel{\→}{i}+{\widehat{\mathrm{\ζ}}}_y\stackrel{\→}{j}+{\widehat{\mathrm{\ζ}}}_z\stackrel{\→}{k}\end{array}$

显然，式(18)、式(19)和式(20)分别为网格变换雅可比的数学定义式(7)、网格变换雅可比的对称形式(9)和网格变换雅可比的对称守恒形式(16)的另一种表现形式，其本质分别与式(7)、式(9)和式(16)相同，这些多种表现形式在所述线性有限差分算子离散下等同，体现了网格变换雅可比计算的唯一性；另一方面网格变换雅可比的多种表现形式又便于在不同物理模型中的应用。

Claims (3)

1.一种检测网格品质的有限差分模拟绕复杂构型流动的方法，它包括：

步骤一、计算网格的生成、检测和调整；

步骤二、控制微分方程的离散和求解；

步骤三、后置处理；

其特征在于在所述步骤一中1)、对整个流动计算区域生成复杂多块对接结构网格，2)、网格变换雅可比的数学定义式、对称形式和对称守恒形式采用线性有限差分算子离散后计算得到的网格变换雅可比分别记为V¹、V²和V³，3)、根据所述V¹、V²和V³来确定网格品质的检测指标V_N和V_M，所述V_N为所述V¹、V²和V³的最小值；所述V_M为V¹²、V²³和V¹³的最大值，其中V¹²、V²³和V¹³由所述V¹、V²和V³按V¹²＝∣V¹‐V²∣/V²、V²³＝∣V²‐V³∣/V³、V¹³＝∣V¹‐V³∣/V³确定，4)、根据所述网格品质的检测指标对生成的计算网格进行检测，确定出那些未达到检测指标的网格单元的具体位置，5)、对所述的未达到检测指标的网格单元进行调整后重新检测，通过这一可重复的过程直到该调整后的网格单元全部达到检测指标；在所述步骤二中对整个流动计算区域采用高阶精度有限差分方法离散控制微分方程，并对其进行求解；在所述步骤三中对数值模拟流场分析其对应的物理意义。

2.根据权利要求1所述的检测网格品质的有限差分模拟绕复杂构型流动的方法，其特征在于采用高阶精度有限差分方法离散时所述网格品质的检测指标为同时满足V_N＞0和V_M≤10％。

3.根据权利要求1或2所述的检测网格品质的有限差分模拟绕复杂构型流动的方法，其特征在于所述网格变换雅可比的三种计算形式均有多种表现形式，这些多种表现形式在所述线性有限差分算子离散下等同。