CN102831304B - 一种运用几何守恒的有限差分模拟绕复杂构型流动的方法 - Google Patents

Abstract

本发明公开了一种运用几何守恒的有限差分模拟绕复杂构型流动的方法，它包括：步骤一、计算网格的生成；步骤二、控制微分方程的离散和求解；步骤三、后置处理；主要是在所述步骤二中需满足：1、采用网格导数的对称守恒计算形式，使得网格导数的计算在保证几何守恒律严格满足的同时又能体现其网格单元矢量面积的几何本质，2、采用网格变换雅可比的对称守恒计算形式，使得离散的网格变换雅可比能够体现网格单元的体积由其矢量面积包围而成的几何本质，3、离散微分方程求导的空间差分算子记为δ¹，守恒网格导数的外层差分算子记为δ²，守恒网格导数的内层差分算子记为δ³，则所述的三个差分算子均为线性差分算子，且它们在各个计算坐标方向上分别相等。本发明大大增强了高阶精度有限差分方法在绕复杂构型流动中的计算稳定性和适应能力。

Description

一种运用几何守恒的有限差分模拟绕复杂构型流动的方法

技术领域

本发明涉及一种复杂计算区域物理特性的数值模拟方法，特别涉及一种绕复杂构型流动的高阶精度的有限差分模拟方法。

背景技术

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

在现有的数值离散方法中，二阶精度的有限体积方法因为其良好的计算稳定性而得到了较大规模的应用；而古老的有限差分方法由于计算稳定性不好，尤其是在复杂多块对接的结构网格中的计算稳定性较差，而限制了其在实际工程问题中的大规模应用；有限元方法的计算量过大也不利于实际工程问题的大规模应用。但是如果考虑到三维空间中高精度、高分辨率数值计算方法的实现过程，有限差分方法由于可以采用逐维离散求导的方式达到高精度、高分辨率，而有限体积方法和有限元方法由于需要采用多维重构才能达到高精度，其计算量比相同精度的有限差分方法大一个量级，甚至更多。在现有条件下，考虑到计算资源的限制条件，基于有限差分离散的高精度数值方法成了解决实际工程问题中绕复杂构型流动数值模拟的最佳技术途径。

计算的稳定性和适应能力是指数值模拟方法对网格构建质量的依赖度。有限体积方法的计算稳定性好，其主要原因在于有限体积方法在计算网格导数时，网格面积和体积都是基于计算网格的几何意义计算的，能够严格满足几何守恒律(意即空间任一封闭的单元体其矢量面积和应为零，在计算中则表现为控制方程自由流守恒的物理特性)，同时又能较好的反映计算网格的几何特性。而有限差分方法在复杂多块对接的结构网格中的计算稳定性较差的一个重要原因就在于其几何守恒律不容易满足，尤其是高阶精度的有限差分方法其几何守恒律更不容易满足。尽管邓小刚等人在《Journal of ComputationalPhysics》Volume 230.Issue 4，20February 2011，Pages 1100-1115提出的守恒网格导数方法能够严格满足几何守恒律，增强高精度有限差分方法对复杂多块对接结构网格的适应能力，但仍存在以下问题：1、守恒网格导数方法只讨论了一种形式的守恒网格导数，且不具备网格导数应当具备的计算网格矢量面积的几何特性，2、网格导数的计算形式未能实质上唯一确定，导致不同的网格导数的计算形式对计算结果的影响也难以评估；3、守恒网格导数方法中对网格导数计算过程中守恒网格导数内层差分算子的选取没有给出充分有力的证明，存在歧义，容易产生误解，导致计算过程不能被正确的实施；4、守恒网格导数方法对网格变换雅可比的计算形式没有提及，在实际计算过程中会存在多种形式的网格变换雅可比，其对计算结果的影响也难以评判。以上问题导致目前的高精度有限差分方法模拟绕复杂构型流动时仍亟待提高计算稳定性和适应能力。

发明内容

本发明的目的是解决运用高精度有限差分方法模拟绕复杂构型流动时的计算稳定性和适应能力问题。

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

一种运用几何守恒的有限差分模拟绕复杂构型流动的方法，它包括：

步骤一、计算网格的生成；

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

步骤三、后置处理；

其特征在于在所述步骤一中，对整个流动计算区域生成复杂多块对接结构网格；在所述步骤二中1、采用网格导数的对称守恒计算形式，使得网格导数的计算在保证几何守恒律严格满足的同时又能体现其矢量面积的几何本质，2、采用网格变换雅可比的对称守恒计算形式，使得离散后的网格变换雅可比能够体现网格单元的体积由其矢量面积包围而成的几何本质，3、离散微分方程求导的空间差分算子记为δ¹，守恒网格导数的外层差分算子记为δ²，守恒网格导数的内层差分算子记为δ³，则所述的三个差分算子均为线性差分算子，且它们在各个计算坐标方向上分别相等；在所述步骤三中对数值模拟流场分析其对应的物理意义。

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

1)由于网格导数的计算采用了其对称守恒形式，所以能够保证几何守恒律的严格满足；

2)在多块对接结构网格中采用线性有限差分算子离散后，网格导数可以表示为网格单元的矢量面积或者是网格单元矢量面积的线性组合的形式；

3)在多块对接结构网格中采用线性有限差分算子离散后，网格变换雅可比可以表示为网格单元体积或者是网格单元体积的线性组合形式；

4)在多块对接结构网格中采用线性有限差分算子离散后，网格变换雅可比由其矢量面积围成且也由其矢量面积计算得出；

5)由于网格导数采用了微分算子的形式，所以能够直接采用高精度、高分辨率的有限差分算子逐维离散网格导数中的微分算子而得到高精度、高分辨率的数值模拟结果；

6)由于网格变换雅可比采用了微分算子的形式，所以能够直接采用高精度、高分辨率的有限差分算子逐维离散网格变换雅可比中的微分算子而得到高精度、高分辨率的数值模拟结果；

7)在较小计算量的情况下，大大提高了高精度、高分辨率有限差分方法在绕复杂构型流动中数值模拟的稳定性和适应能力。

此外，本发明所述网格导数的对称守恒计算形式存在多种表现形式，这些多种表现形式在所述线性差分算子离散下等同。

本发明所述网格变换雅克比的对称守恒计算形式存在多种表现形式，这些多种表现形式在所述线性差分算子离散下等同。

附图说明

图1为具体实施方式的网格导数的计算网格示意图

图2为具体实施方式的对称守恒网格导数(8)计算得到的矢量面积的示意图(二阶中心差分格式)

图3为具体实施方式的原始网格导数(5)计算得到的矢量面积的示意图(二阶中心差分格式)

图4为类有限体积方法计算得到的矢量面积的示意图(二阶中心差分格式)

图5为具体实施方式的对称守恒网格导数(8)计算得到的矢量面积的示意图(一阶前差格式)

图6为具体实施方式的对称守恒网格导数(8)计算得到的矢量面积的示意图(一阶后差格式)

图7为具体实施方式的对称守恒网格导数(8)计算得到的矢量面积的示意图(方向：一阶前差格式；方向：二阶中心差分格式)

图8为具体实施方式的对称守恒网格导数(8)计算得到的矢量面积的示意图(方向：二阶中心差分格式；方向：一阶后差格式)

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

具体实施方式

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

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

其中：

$\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---\left(4\right)$

$\hat{G}={\widehat{\mathrm{\ζ}}}_{t}Q+{\widehat{\mathrm{\ζ}}}_{x}E+{\widehat{\mathrm{\ζ}}}_{y}F+{\widehat{\mathrm{\ζ}}}_{z}G$

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

$\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(5\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.$

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

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

J^-1＝x_ξy_ηz_ζ-x_ξy_ζz_η+x_ηy_ζz_ξ-x_ηy_ξz_ζ+x_ζy_ξz_η-x_ζy_ηz_ξ    (6)

几何守恒律的严格满足要求I_x＝I_y＝I_z＝0，其中：

$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{\ζ}}---\left(7\right)$

$I_z=\frac{\∂{\widehat{\mathrm{\ξ}}}_z}{\∂\mathrm{\ξ}}+\frac{\∂{\widehat{\mathrm{\η}}}_z}{\∂\mathrm{\η}}+\frac{\∂{\widehat{\mathrm{\ζ}}}_z}{\∂\mathrm{\ζ}}$

为了能够严格满足几何守恒律，同时又能准确地反映计算网格的几何特性，本发明给出如下的网格导数的对称守恒的计算形式：

$\left\{\begin{array}{c}{\widehat{\mathrm{\ξ}}}_x=\frac{1}{2}[{\left({\mathrm{zy}}_{\mathrm{\η}}\right)}_{\mathrm{\ζ}}+{\left(y{z}_{\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(y{z}_{\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(y{z}_{\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(8\right)$

由网格导数的数学定义式(5)可以看出，离散的网格导数的几何本质是网格单元的矢量面积，式(5)中的各个几何量实为网格单元矢量面积在直角坐标系下的分量。本发明中给出的网格导数的计算形式(8)采用合适的线性有限差分算子离散时亦可以保持其网格单元矢量面积的几何特性。

网格变换雅可比的数学定义式(6)的几何本质是网格单元的体积。考虑到体积的几何定义，任意单元的体积由其矢量面积围成，也可由其矢量面积唯一确定。鉴于此，本发明给出的网格变换雅可比的对称守恒计算形式：

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

网格变换雅可比的几何本质是网格单元的体积。在对网格变换雅可比的对称守恒计算形式采用线性有限差分算子离散时，如果同时采用对称守恒的网格导数，则式(9)计算得到的体积刚好为图9中的整个体积。

本发明中所指的有限差分算子δ仅限线性差分算子，线性差分算子δ具有如下特性：

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

进一步的推导可以证明：线性差分算子还满足求导顺序可交换性质，即：

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

将对称守恒形式(8)中的网格导数的外层差分算子记为δ²，内层差分算子记为δ³，将式(8)中的微分导数表述为差分算子的形式：

$\left\{\begin{array}{c}{\widehat{\mathrm{\ξ}}}_x=\frac{1}{2}[{\mathrm{\δ}}_{\mathrm{\ζ}}^2\left(z{\mathrm{\δ}}_{\mathrm{\η}}^{3}y\right)+{\mathrm{\δ}}_{\mathrm{\η}}^2\left(y{\mathrm{\δ}}_{\mathrm{\ζ}}^{3}z\right)-{\mathrm{\δ}}_{\mathrm{\η}}^2\left(z{\mathrm{\δ}}_{\mathrm{\ζ}}^{3}y\right)-{\mathrm{\δ}}_{\mathrm{\ζ}}^2\left(y{\mathrm{\δ}}_{\mathrm{\η}}^{3}z\right)]\\ {\widehat{\mathrm{\ξ}}}_y=\frac{1}{2}[{\mathrm{\δ}}_{\mathrm{\ζ}}^2\left(x{\mathrm{\δ}}_{\mathrm{\η}}^{3}z\right)+{\mathrm{\δ}}_{\mathrm{\η}}^2\left(z{\mathrm{\δ}}_{\mathrm{\ζ}}^{3}x\right)-{\mathrm{\δ}}_{\mathrm{\η}}^2\left(x{\mathrm{\δ}}_{\mathrm{\ζ}}^{3}z\right)-{\mathrm{\δ}}_{\mathrm{\ζ}}^2\left(z{\mathrm{\δ}}_{\mathrm{\η}}^{3}x\right)]\\ {\widehat{\mathrm{\ξ}}}_z=\frac{1}{2}[{\mathrm{\δ}}_{\mathrm{\ζ}}^2\left(y{\mathrm{\δ}}_{\mathrm{\η}}^{3}x\right)+{\mathrm{\δ}}_{\mathrm{\η}}^2\left(x{\mathrm{\δ}}_{\mathrm{\ζ}}^{3}y\right)-{\mathrm{\δ}}_{\mathrm{\η}}^2\left(y{\mathrm{\δ}}_{\mathrm{\ζ}}^{3}x\right)-{\mathrm{\δ}}_{\mathrm{\ζ}}^2\left(x{\mathrm{\δ}}_{\mathrm{\η}}^{3}y\right)]\end{array}\right.$

$\left\{\begin{array}{c}{\widehat{\mathrm{\η}}}_x=\frac{1}{2}[{\mathrm{\δ}}_{\mathrm{\ξ}}^2\left(z{\mathrm{\δ}}_{\mathrm{\ζ}}^{3}y\right)+{\mathrm{\δ}}_{\mathrm{\ζ}}^2\left(y{\mathrm{\δ}}_{\mathrm{\ξ}}^{3}z\right)-{\mathrm{\δ}}_{\mathrm{\ζ}}^2\left(z{\mathrm{\δ}}_{\mathrm{\ξ}}^{3}y\right)-{\mathrm{\δ}}_{\mathrm{\ξ}}^2\left(y{\mathrm{\δ}}_{\mathrm{\ζ}}^{3}z\right)]\\ {\widehat{\mathrm{\η}}}_y=\frac{1}{2}[{\mathrm{\δ}}_{\mathrm{\ξ}}^2\left(x{\mathrm{\δ}}_{\mathrm{\ζ}}^{3}z\right)+{\mathrm{\δ}}_{\mathrm{\ζ}}^2\left(z{\mathrm{\δ}}_{\mathrm{\ξ}}^{3}x\right)-{\mathrm{\δ}}_{\mathrm{\ζ}}^2\left(x{\mathrm{\δ}}_{\mathrm{\ξ}}^{3}z\right)-{\mathrm{\δ}}_{\mathrm{\ξ}}^2\left(z{\mathrm{\δ}}_{\mathrm{\ζ}}^{3}x\right)]\\ {\widehat{\mathrm{\η}}}_z=\frac{1}{2}[{\mathrm{\δ}}_{\mathrm{\ξ}}^2\left(y{\mathrm{\δ}}_{\mathrm{\ζ}}^{3}x\right)+{\mathrm{\δ}}_{\mathrm{\ζ}}^2\left(x{\mathrm{\δ}}_{\mathrm{\ξ}}^{3}y\right)-{\mathrm{\δ}}_{\mathrm{\ζ}}^2\left(y{\mathrm{\δ}}_{\mathrm{\ξ}}^{3}x\right)-{\mathrm{\δ}}_{\mathrm{\ξ}}^2\left(x{\mathrm{\δ}}_{\mathrm{\ζ}}^{3}y\right)]\end{array}\right.$

$\left\{\begin{array}{c}{\widehat{\mathrm{\ζ}}}_x=\frac{1}{2}[{\mathrm{\δ}}_{\mathrm{\η}}^2\left(z{\mathrm{\δ}}_{\mathrm{\ξ}}^{3}y\right)+{\mathrm{\δ}}_{\mathrm{\ξ}}^2\left(y{\mathrm{\δ}}_{\mathrm{\η}}^{3}z\right)-{\mathrm{\δ}}_{\mathrm{\ξ}}^2\left(z{\mathrm{\δ}}_{\mathrm{\η}}^{3}y\right)-{\mathrm{\δ}}_{\mathrm{\η}}^2\left(y{\mathrm{\δ}}_{\mathrm{\ξ}}^{3}z\right)]\\ {\widehat{\mathrm{\ζ}}}_y=\frac{1}{2}[{\mathrm{\δ}}_{\mathrm{\η}}^2\left(x{\mathrm{\δ}}_{\mathrm{\ξ}}^{3}z\right)+{\mathrm{\δ}}_{\mathrm{\ξ}}^2\left(z{\mathrm{\δ}}_{\mathrm{\η}}^{3}x\right)-{\mathrm{\δ}}_{\mathrm{\ξ}}^2\left(x{\mathrm{\δ}}_{\mathrm{\η}}^{3}z\right)-{\mathrm{\δ}}_{\mathrm{\η}}^2\left(z{\mathrm{\δ}}_{\mathrm{\ξ}}^{3}x\right)]\\ {\widehat{\mathrm{\ζ}}}_z=\frac{1}{2}[{\mathrm{\δ}}_{\mathrm{\η}}^2\left(y{\mathrm{\δ}}_{\mathrm{\ξ}}^{3}x\right)+{\mathrm{\δ}}_{\mathrm{\ξ}}^2\left(x{\mathrm{\δ}}_{\mathrm{\η}}^{3}y\right)-{\mathrm{\δ}}_{\mathrm{\ξ}}^2\left(y{\mathrm{\δ}}_{\mathrm{\η}}^{3}x\right)-{\mathrm{\δ}}_{\mathrm{\η}}^2\left(x{\mathrm{\δ}}_{\mathrm{\ξ}}^{3}y\right)]\end{array}\right.---\left(12\right)$

在本发明中，当差分算子满足如下条件时：

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

则式(12)在线性有限差分算子离散后总是可以表示为网格单元矢量面积的形式。以图1所示的计算网格示意图为例说明：如果δ²和δ³在计算坐标系ξ和η方向均为二阶中心差分，则式(12)计算得到的矢量面积如图2所示(图2至图8中的矢量面积均指阴影部分的面积)。采用相同的差分算子，则网格导数的数学定义式(5)计算得到的矢量面积如图3所示。而类有限体积方法给出的矢量面积如图4所示。当差分算子δ²和δ³不为中心差分时，如δ²和δ³在计算坐标系ξ和η方向均为一阶向前差分格式时，则式(12)计算得到的矢量面积如图5所示；同理图6则为δ²和δ³在计算坐标系ξ和η方向均为一阶向后差分格式计算的矢量面积。注意到式(13)中并不要求不同计算坐标系方向的差分算子要一样，当两个坐标方向的差分算子不一样时，如δ²和δ³在计算坐标系ξ方向均为一阶向前格式，而在计算坐标系η方向均为二阶中心差分格式，其矢量面积如图7所示，图8则表示另外一种情况下的矢量面积：δ²和δ³在计算坐标系ξ方向均为二阶中心差分格式，而在计算坐标系η方向均为一阶向后差分格式。

在对方程(3)进行有限差分离散时，将方程(3)中的空间差分算子记为δ¹，则有自由流守恒(几何守恒律)要求：

$I_x={\mathrm{\δ}}_{\mathrm{\ξ}}^1\left({\widehat{\mathrm{\ξ}}}_x\right)+{\mathrm{\δ}}_{\mathrm{\η}}^1\left({\widehat{\mathrm{\η}}}_x\right)+{\mathrm{\δ}}_{\mathrm{\ζ}}^1\left({\widehat{\mathrm{\ζ}}}_x\right)$

$I_y={\mathrm{\δ}}_{\mathrm{\ξ}}^1\left({\widehat{\mathrm{\ξ}}}_y\right)+{\mathrm{\δ}}_{\mathrm{\η}}^1\left({\widehat{\mathrm{\η}}}_y\right)+{\mathrm{\δ}}_{\mathrm{\ζ}}^1\left({\widehat{\mathrm{\ζ}}}_y\right)---\left(14\right)$

$I_z={\mathrm{\δ}}_{\mathrm{\ξ}}^1\left({\widehat{\mathrm{\ξ}}}_z\right)+{\mathrm{\δ}}_{\mathrm{\η}}^1\left({\widehat{\mathrm{\η}}}_z\right)+{\mathrm{\δ}}_{\mathrm{\ζ}}^1\left({\widehat{\mathrm{\ζ}}}_z\right)$

如果式(14)中的网格导数采用如式(8)中的守恒计算形式，则依照线性有限差分算子的固有特性(11)，容易给出严格满足几何守恒律的充分条件为：

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

结合式(13)和式(15)，本发明给出的线性有限差分算子应满足：

${\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(16\right)$

由图2和图9可以看出网格导数的对称守恒形式离散后即为网格单元的矢量面积，网格变换雅可比完全由其矢量面积围成且也由其矢量面积计算得出，网格导数和网格变换雅可比各自的几何本质得到了充分体现，有利于增强有限差分方法的数值计算稳定性和适应能力。在较小计算量的情况下，大大提高了高精度、高分辨率有限差分方法在绕复杂构型流动中数值模拟的稳定性和适应能力。对数值模拟结果进行后置处理能够实现对计算流域内物理变量Q较好的数值模拟。

网格导数的对称守恒计算形式(8)还可以简单写为：

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

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

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

以及其对应的矢量运算形式：

${\stackrel{\→}{S}}^{\mathrm{\ξ}}={\widehat{\mathrm{\ξ}}}_x\stackrel{\→}{i}+{\widehat{\mathrm{\ξ}}}_y\stackrel{\→}{j}+{\widehat{\mathrm{\ξ}}}_x\stackrel{\→}{k}=\frac{1}{2}[{({\stackrel{\→}{r}}_{\mathrm{\η}}\×\stackrel{\→}{r})}_{\mathrm{\ζ}}-{({\stackrel{\→}{r}}_{\mathrm{\ζ}}\×\stackrel{\→}{r})}_{\mathrm{\η}}]$

${\stackrel{\→}{S}}^{\mathrm{\η}}={\widehat{\mathrm{\η}}}_x\stackrel{\→}{i}+{\widehat{\mathrm{\η}}}_y\stackrel{\→}{j}+{\widehat{\mathrm{\η}}}_z\stackrel{\→}{k}=\frac{1}{2}[{({\stackrel{\→}{r}}_{\mathrm{\ζ}}\×\stackrel{\→}{r})}_{\mathrm{\ξ}}-{({\stackrel{\→}{r}}_{\mathrm{\ξ}}\×\stackrel{\→}{r})}_{\mathrm{\ζ}}]---\left(18\right)$

${\stackrel{\→}{S}}^{\mathrm{\ζ}}={\widehat{\mathrm{\ζ}}}_x\stackrel{\→}{i}+{\widehat{\mathrm{\ζ}}}_y\stackrel{\→}{j}+{\widehat{\mathrm{\ζ}}}_z\stackrel{\→}{k}=\frac{1}{2}[{({\stackrel{\→}{r}}_{\mathrm{\ξ}}\×\stackrel{\→}{r})}_{\mathrm{\η}}-{({\stackrel{\→}{r}}_{\mathrm{\η}}\×\stackrel{\→}{r})}_{\mathrm{\ξ}}]$

显然式(17)和(18)在采用上述线性差分算子离散后均与式(8)的离散形式(12)等同。同理，式(19)与式(9)网格变换雅可比在采用上述线性差分算子离散后等同。

因此，在本发明中，可以说网格导数的对称守恒计算形式实质上是唯一确定的，网格变换雅可比的对称守恒形式实质上也是唯一确定的，这样就大大提高了数值模拟计算的稳定性和适应能力；网格导数的对称守恒计算形式和网格变换雅克比的对称守恒计算形式存在多种表现形式有利于不同物理模型中的应用。

在唯一确定了网格导数和网格变换雅可比的计算形式，消除了计算过程中的网格导数和网格变换雅可比的不确定性后，对数值模拟结果进行后置处理能够实现对计算流域内物理变量Q的较为理想的数值模拟。

Claims (3)

1.一种运用几何守恒的有限差分模拟绕复杂构型流动的方法，它包括：

步骤一、计算网格的生成；

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

步骤三、后置处理；

其特征在于在所述步骤一中，对整个流动计算区域生成复杂多块对接结构网格；在所述步骤二中1、采用网格导数的对称守恒计算形式，使得网格导数的计算在保证几何守恒律严格满足的同时又能体现其矢量面积的几何本质，2、采用网格变换雅可比的对称守恒计算形式，使得离散后的网格变换雅可比能够体现网格单元的体积由其矢量面积包围而成的几何本质，3、离散微分方程求导的空间差分算子记为δ¹，守恒网格导数的外层差分算子记为δ²，守恒网格导数的内层差分算子记为δ³，则所述的三个差分算子均为线性差分算子，且它们在各个计算坐标方向上分别相等；在所述步骤三中对数值模拟流场分析其对应的物理意义。

2.根据权利要求1所述的运用几何守恒的有限差分模拟绕复杂构型流动的方法，其特征在于所述网格导数的对称守恒计算形式存在多种表现形式，这些多种表现形式在所述线性差分算子离散下等同。

3.根据权利要求1所述的运用几何守恒的有限差分模拟绕复杂构型流动的方法，其特征在于所述网格变换雅克比的对称守恒计算形式存在多种表现形式，这些多种表现形式在所述线性差分算子离散下等同。