CN103914602A - 可压缩旋流场的数值模拟方法 - Google Patents

Abstract

本发明是一种模拟可压缩旋流场的数值方法。其主要特征是在计算网格边界上加入一个涡量耗散补偿项通量，补偿因为涡量的耗散方向和涡量的模的最大变化梯度方向的分离而造成的涡量过度耗散。这个涡量耗散补偿项应该与人工耗散项等价、作用相反。该方法被称为涡量耗散的补偿技术。这个被加入涡量耗散补偿项是在对流项加入人工耗散项以后加入的，不影响原来数值解的稳定性策略、是由数值粘性引起的、方向是沿着涡量的耗散项与涡量的模的最大梯度变化方向的叉乘方向。

Description

可压缩旋流场的数值模拟方法

技术领域

本发明涉及计算流体力学（CFD：Computational Fluid Dynamics）领域中的一种数值模拟方法，具体是一种模拟可压缩旋流场的数值方法。

背景技术

计算流体力学综合了流体力学、应用数学、计算机科学，是一门应用性极强的学科。流体力学问题的数值模拟以其低成本、直观性强的优势，在流体流动的机理探索、工业产品设计等各个相关领域占据重要地位。计算流体力学面临的最大的问题和挑战之一即是如何提高数值模拟的精度，忠实地表现流体流动的特性。

影响流体力学问题的数值模拟的精度的重要因素之一是：当使用数值方法求解流体控制方程，即欧拉（Euler）方程或者纳维尔-斯托克斯(Navier-Stokes)方程时，会产生数值耗散(numerical diffusion)，造成数值解的误差。例如数值方法中对控制方程中的对流项的各种空间离散方法（如中心差分、迎风差分）、时间离散方法（如显示时间积分、隐式时间积分）、湍流模型（如双方程模型、大涡模拟）的使用，以及计算网格的正交性都会产生不同程度的数值耗散。此外，数值模拟中还经常要使用一种人工数值耗散(artificial diffusion)技术，其目的是通过适当降低计算精度而获得稳定的数值解。例如流场中的不连续(discontinuity )界面的捕捉，如激波（shock）的捕捉，即是依靠加入适量的人工耗散项，以避免二阶以上精度的数值解在流动变量梯度较大的地方出现数值振荡现象。数值耗散可以理解为是流场中的一种能量损失，这种能量损失在某种程度上使得数值模拟结果不能忠实的体现流体的流动特性，降低了计算精度。先进的数值方法应该是在保证获得稳定性的数值解的前提下，将数值耗散减至最小。

数值耗散对流场的数值模拟结果最明显的影响体现在对流动变量的间断界面的捕捉。如前所述，激波是强间断界面，跨过激波， 速度、密度，压力均不连续，对其捕捉必须加入一定的人工数值耗散，以抵消数值解在不连续处的振荡。因为激波前后存在熵增，即能量的损失，所以，通过加入人工数值耗散捕捉激波具有合理的物理意义。但是，过大的数值耗散会使数值解获得的激波界面变得模糊，降低了数值解对激波强度和空间位置的预测精度。流场中存在另一类流动不连续现象，即接触不连续（contact discontinuity）。相对激波而言，接触不连续的特点是弱不连续，不连续界面是滑移线（slip-line）,跨过滑移线，没有质量传递、密度和切向速度不连续，但径向速度和压力连续。一类旋流场（Vortex-dominated Flows）中具有最典型的接触不连续界面的流动现象。例如，具有三角机翼的飞行器，具有较大展翼比，在大攻角飞行时周围的流场是一明显的旋涡运动为主。其他例如直升机旋翼后方的流场、涡轮发动机转子周围的流场，等等，同样会遇到类似的现象。计算机的数值模拟中对于这种接触不连续的捕捉更加困难，因为数值方法中的数值耗散即使很小也会使弱不连续界面变得模糊，降低了数值解对流场的预测精度，这也是旋流场的数值模拟技术成为CFD领域的重大挑战的原因。

描述可压缩粘性流体运动的控制方程包括连续方程和动量方程，即Navier-Stokes方程，分别由以下两式给出，

$\frac{\∂\mathrm{\ρ}}{\∂t}+\▿\·\left(\mathrm{\ρ}\stackrel{\→}{V}\right)=0,---\left(1\right)$

$\frac{\∂\stackrel{\→}{V}}{\∂t}+(\stackrel{\→}{V}\·\▿)\stackrel{\→}{V}=-\frac{1}{\mathrm{\ρ}}\▿p+\frac{1}{3}\frac{\mathrm{\μ}}{\mathrm{\ρ}}\▿(\▿\·\stackrel{\→}{V})+\frac{\mathrm{\μ}}{\mathrm{\ρ}}{\▿}^2\stackrel{\→}{V}+\stackrel{\→}{B},---\left(2\right)$

式中，是速度矢量，，包含直角坐标系下的三维i,j,k分量 u，v，w；ρ、p、t分别是密度、压力和时间。式中算子表示为；符号·代表内积计算。动量方程（2）中的代表压缩性引起的剪切动力粘性项、代表正压运动粘性项（是拉普拉斯算子）、代表体积力项。对于一类以旋涡运动为主的流场，常引入描述旋流流动的变量—涡量（vorticity）对流场进行描述。涡量用表示，在直角坐标系下有。按照涡量的定义有

$\stackrel{\→}{\mathrm{\ω}}=\▿\×\stackrel{\→}{V}\left|\begin{array}{ccc}i& j& k\\ \frac{\∂}{\∂x}& \frac{\∂}{\∂y}& \frac{\∂}{\∂z}\\ u& v& w\end{array}\right|=(\frac{\∂w}{\∂y}-\frac{\∂v}{\∂z})i+(\frac{\∂u}{\∂z}-\frac{\∂w}{\∂x})j+(\frac{\∂v}{\∂x}-\frac{\∂u}{\∂y})k,---\left(3\right)$

式中符号×代表叉乘运算。

对方程（2）等号两边作叉乘运算，并利用方程（1）和涡量的定义，可以推导出描述可压缩粘性旋涡流场的涡量方程，

$\frac{\∂\stackrel{\→}{\mathrm{\ω}}}{\∂t}+(\stackrel{\→}{V}\·\▿)\stackrel{\→}{\mathrm{\ω}}=(\stackrel{\→}{\mathrm{\ω}}\·\▿)\stackrel{\→}{V}-\stackrel{\→}{\mathrm{\ω}}(\▿\·\stackrel{\→}{V})+v{\▿}^2\mathrm{\ω}+\frac{\▿\mathrm{\ρ}\×\▿p}{{\mathrm{\ρ}}^2}+\▿\×(\frac{1}{3}v\▿(\▿\·\stackrel{\→}{V}))+\▿\·\stackrel{\→}{B}.---\left(4\right)$

上式使用了可压缩流的动力粘度v，定义为v=μ/ρ。公式（4）中等号左边表示涡量的变化率；右边第一项表示涡量的三维变形产生的涡量；第二项表示可压缩性对涡量的影响；第三项和第五项表示粘性产生的涡量耗散；第四项表示压力梯度和密度梯度不平行时产生的涡量；第六项表示体积产生涡量的作用。公式（4）对于二维情况下，；对于正压流场，；对于无粘流，且，该式可以被简化。

对于涡量场的影响最显著的是涡量的耗散项，主要是因为流体粘性作用引起的，同样的道理，数值方法中需要加入的人工耗散是以粘性形式起作用，直接导致了流体粘性引起的涡量耗散以外的涡量耗散。流体的粘性作用，包括人工粘性的作用具有椭圆型特征，即粘性作用在空间是各向同性的。

为了提高旋流场的数值模拟的精度，一类方法是加密计算网格，在更加细小的空间尺度内求解流体控制方程。加密计算网格首先会使计算量加大，增加计算成本。此外，数值计算的误差随着计算网格的增加会不断积累，在一定程度上造成相反的效果。另一类方法是在流场中采用物理模型来增加流场中涡量的强度。例如在流场中加入点涡模型，可以人为地增加涡量，或者在流场局部直接求解涡量方程（公式4），以减小涡量的输运过程中的数值耗散。但是，这些方法在应用上仍受到一定限制。点涡模型是在预先明确旋涡发生位置的前提下才能使用，仅适合一些简单的流动现象。除了二维不可压缩正压流场，涡量方程的求解比与欲求解的Euler、Navier-Stokes方程更为复杂。最后一类方法是开发具有最小数值耗散的空间离散算法，例如，对流项离散方法ENO、WENO方法，主要是针对求解直角坐标系下控制方程（1）和（2）的求解过程中使用的。但是，这些数值方法的数值耗散依然过大，会降低弱不连续的捕捉精度，使不连续界面变得模糊。

总之，为了更加精确地模拟可压缩流的旋涡运动，需要进一步改进对于流场中诸如旋涡运动一类的接触不连续的捕捉机制。一种途径是按照公式（4）给出的涡量变化的原理，在流场中补偿人工耗散项所产生的涡量耗散，同时在数值解求解构成中将补偿的涡量耗散进行高精度空间离散，更精确地模拟流场中的旋涡运动，形成一种新的可压缩流的旋涡运动的数值模拟技术。

发明内容

本发明涉及计算流体力学领域中一种模拟可压缩旋流场的数值方法，具体是一种在流场中补偿人工耗散项产生的涡量耗散，同时在数值解求解构成中将补偿的涡量耗散进行高精度空间离散，更精确地模拟可压缩流的旋涡运动的数值模拟方法—涡量耗散的补偿技术。

本发明采用的技术方案：

可压缩流体流动的控制方程是一组偏微分方程，包括连续、动量和能量方程。无粘流的动量方程又被称为Euler方程；粘性流动量方程又被称为Navier-Stokes方程。工程应用上，流体流动的控制方程的数值解的求解常采用有限体积法（Finite Volume Method）将控制方程在预先生成计算网格上进行离散化。有限体积法需要将控制方程写成守恒积分形式，每一个计算网格单元被是为控制体，控制体内流场守恒变量随时间的变化是通过在每一个特定时间间隔内，通过控制体表面（即计算网格边界）上的守恒变量的通量值的和来确定。这也是的流体控制方程数值解的基本原理。上述通过计算网格界面的通量值包括对流通量项、人工耗散项和粘性项。各种不同的空间离散方法，如中心差分的代表JST方法、考虑偏微分方程组的特征的代表VanLeer、CUSP、AUSM、Roe等方法，都是用来逼近在计算网格边界上的对流项和人工耗散项通量值。粘性项的空间离散，因为粘性项的椭圆形特性，均采用中心差分。

如前所述，空间离散的数值方法中需要加入一定的人工数值耗散。其加入的形式可以是直接的，如JST、CUSP、Roe或者是隐含式的，如VanLeer、AUSM方法。人工耗散在使偏微分方程组获得数值解稳定的作用也降低了旋涡结构的捕捉精度，构成了数值模拟误差的一个来源。对于旋流场，加入的人工耗散使得流场中的涡量过度耗散。如前所述，对于涡量场的影响最显著的是涡量的耗散项，所以要在旋流场中减小因加入人工耗散而造成的涡量耗散。一个使得旋流场的涡量耗散保持最小化的原则是：涡量的耗散方向和涡量的模的最大变化梯度方向保持一致。而人工耗散的加入使得这两个方向分离。因此，本发明采取的解决方法是，在计算网格边界上加入一个涡量耗散补偿项通量，补偿因为涡量的耗散方向和涡量的模的最大变化梯度方向的分离而造成的涡量过度耗散。这个涡量耗散补偿项应该与人工耗散项等价、作用相反。该方法被称为涡量耗散的补偿技术。其原理图见图1。图中表示，通过计算网格的界面，通常是由两个通量，即对流项通量（包含粘性项通量）和加入的人工耗散项，通过加入与人工耗散项相反的作用、等价的涡量耗散补偿项通量，可以抵消部分数值耗散。涡量耗散的补偿技术有以下三个特点：

1.这个加入涡量耗散补偿项是在对流项加入人工耗散项以后加入的，所以不影响原来数值解的稳定性策略；

2.这个被加入的涡量耗散补偿项不是由有流体粘性引起的，而由数值粘性引起的；

3.涡量的耗散项与涡量的的模的最大梯度变化的方向一致，则不需要涡量耗散补偿项。

按照量纲分析，涡量耗散补偿项是密度和单位涡量耗散补偿项的乘积，定义为

$\stackrel{\→}{f}={\stackrel{\→}{n}}_{\mathrm{\ω}}\×\left(v_{\mathrm{\ω}}\left({\▿}^2\stackrel{\→}{\mathrm{\ω}}\right)R_c^2\right),---\left(5\right)$

式中，代表涡量的模的梯度最大变化的方向；v_ω代表数值粘性系数，被称为涡运动的补偿粘性系数，与流体的运动粘度有相同的量纲，是一个标量；R_c代表补偿涡量的特征半径。式中符号×代表叉乘运算。当地涡量的模的梯度最大变化的方向表示为

${\stackrel{\→}{n}}_{\mathrm{\ω}}=\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|},---\left(6\right)$

其中，φ是涡量的模；是涡量的模的梯度；是涡量的模的梯度的模，分别定义为

$\mathrm{\φ}=\left|\stackrel{\→}{\mathrm{\ω}}\right|=\sqrt{{\mathrm{\ω}}_x^2+{\mathrm{\ω}}_y^2+{\mathrm{\ω}}_z^2},---\left(7\right)$

$\▿\mathrm{\φ}=\frac{\∂\mathrm{\φ}}{\∂x}i+\frac{\∂\mathrm{\φ}}{\∂y}j+\frac{\∂\mathrm{\φ}}{\∂z}k,---\left(8\right)$

$|\▿\mathrm{\φ}|=\sqrt{{\left(\frac{\∂\mathrm{\φ}}{\∂x}\right)}^2+{\left(\frac{\∂\mathrm{\φ}}{\∂y}\right)}^2+{\left(\frac{\∂\mathrm{\φ}}{\∂z}\right)}^2}.---\left(9\right)$

涡量的模的梯度的模表示了最大涡量的变化程度。

数值粘性系数的确定

因为人工耗散项是以流体的粘性耗散的形式体现的，对于动量方程具有速度的拉普拉斯算子（）的量纲，对具有速度的拉普拉斯算子的量纲人工耗散项取旋度（作运算），则有

$\▿\×\left({\▿}^2\stackrel{\→}{V}\right)={\▿}^2\stackrel{\→}{\mathrm{\ω}}+(\▿\left({\▿}^2\right))\×\stackrel{\→}{V}.---\left(10\right)$

如果舍去式中右边的三阶导数项，速度的拉普拉斯算子的旋度完全变成涡量的拉普拉斯算子，既是涡量方程（4）中的耗散项，说明了公式（2）中的人工耗散项与公式（3）中的涡量耗散项的等价性。所以涡量补偿也就是以涡量的形式补偿了数值方法中的人工耗散，消除了这部分误差，提高数值模拟的精度。

各种空间离散的数值方法中的人工耗散项均需要用一个具有速度量纲的参数作为人工耗散的尺度放大系数。例如，在JST方法中使用谱半径（spectral radii）作为人工耗散的尺度放大系数。谱半径是指偏微分方程在计算网格边界处的最大特征值。对流项的谱半径Λ是

Λ=|V|+c，（11）

式中c是当地的声速；V是抗变速度，定义为，是计算网格界面的法线方向向量。如果对具有速度量纲的人工耗散的放大尺度系数取旋度，可得到具有涡量的量纲的参数（见公式3）。这个放大系数可以用涡量的形式表示，再用来放大涡量耗散项，而且这个系数具有流体运动粘度的量纲，因此，该系数可以被认为是一个数值粘性引起的，被称为涡运动补偿粘性系数。

因为涡运动补偿是在空间不同方向进行的，涡运动补偿粘性系数向量定义为

${\stackrel{\→}{v}}_{\mathrm{\ω}}=\frac{{\stackrel{\→}{R}}_{\mathrm{\ω}}^2}{\left|\stackrel{\→}{\mathrm{\ω}}\right|},---\left(12\right)$

其中，补偿涡量的半径向量，定义成

${\stackrel{\→}{R}}_{\mathrm{\ω}}=R_c\stackrel{\→}{k}.---\left(14\right)$

上式和公式（5）中，补偿涡量的特征半径R_c在二维和三维空间分别定义为

$R_c=\frac{1}{2}{\mathrm{\Ω}}^{1/3},---\left(13\right)$

$R_c=\frac{1}{2}{\mathrm{\Ω}}^{1/3}.---\left(14\right)$

补偿涡量的特征半径的几何意义解释：在计算网格上的补偿涡量一定有一个与之等效的环量，而环量是围绕该网格内某一封闭区域上的速度的线积分。该封闭区域的大小用公式（13）或（14） 表示，计算网格单元的空间尺度的二分之一作路径形成的区域的半径，是涡量作用在该计算网格内的平均空间范围。

公式（14）中，涡量在最大作用的方向系数向量，。其作用是为了在流体的粘性层内减小涡量补偿，以及考虑到计算网格的大变形率，需要根据涡量向量的各向异性程度进行补偿涡量的半径向量的坐标方向的各向异性化，则

$k^I=1+\max \left[\right|\frac{{\mathrm{\ω}}_y}{{\mathrm{\ω}}_x}|,|\frac{{\mathrm{\ω}}_z}{{\mathrm{\ω}}_x}\left|\right],---\left(15\right)$

$k^J=1+\max \left[\right|\frac{{\mathrm{\ω}}_x}{{\mathrm{\ω}}_y}|,|\frac{{\mathrm{\ω}}_z}{{\mathrm{\ω}}_y}\left|\right],---\left(16\right)$

$k^K=1+\max \left[\right|\frac{{\mathrm{\ω}}_x}{{\mathrm{\ω}}_z}|,|\frac{{\mathrm{\ω}}_y}{{\mathrm{\ω}}_z}\left|\right].---\left(17\right)$

最终，在计算网格界面上的抗变涡运动补偿粘性系数v_ω定义为

$v_{\mathrm{\ω}}\≡{\stackrel{\→}{v}}_{\mathrm{\ω}}\·\stackrel{\→}{n}=\frac{1}{\left|\stackrel{\→}{\mathrm{\ω}}\right|}[{\left(R_{\mathrm{\ω}}^I\right)}^2{n}_x+{\left(R_{\mathrm{\ω}}^J\right)}^2{n}_y+{\left(R_{\mathrm{\ω}}^K\right)}^2{n}_z].---\left(18\right)$

按照以上推导，图2给出计算涡量耗散补偿项的步骤流程。

公式（5）给出的涡量耗散补偿项可以在三维直角坐标系下展开，其形式表示为

${\stackrel{\→}{F}}_{\mathrm{\ω}}=\frac{{\mathrm{\ρv}}_{\mathrm{\ω}}{R}_c^2}{|\▿\mathrm{\φ}|}\left[\begin{array}{c}0\\ 0\\ \frac{\∂\mathrm{\φ}}{\∂z}\frac{{\∂}^2{\mathrm{\ω}}_z}{\∂x^2}-\frac{\∂\mathrm{\φ}}{\∂x}\frac{{\∂}^2{\mathrm{\ω}}_z}{\∂z^2}\\ \frac{\∂\mathrm{\φ}}{\∂x}\frac{{\∂}^2{\mathrm{\ω}}_z}{\∂y^2}-\frac{\∂\mathrm{\φ}}{\∂y}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂x}^2}\\ v(\frac{\∂\mathrm{\φ}}{\∂z}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂x}^2}-\frac{\∂\mathrm{\φ}}{\∂x}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂z}^2})+w(\frac{\∂\mathrm{\φ}}{\∂x}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂y}^2}-\frac{\∂\mathrm{\φ}}{\∂y}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂x}^2})\end{array}\right]n_x+$

$\frac{{\mathrm{\ρv}}_{\mathrm{\ω}}{R}_c^2}{|\▿\mathrm{\φ}|}\left[\begin{array}{c}0\\ \frac{\∂\mathrm{\φ}}{\∂\mathrm{\∂}}\frac{{\∂}^2{\mathrm{\ω}}_z}{\∂z^2}-\frac{\∂\mathrm{\φ}}{\∂z}\frac{{\∂}^2{\mathrm{\ω}}_z}{\∂y^2}\\ 0\\ \frac{\∂\mathrm{\φ}}{\∂x}\frac{{\∂}^2{\mathrm{\ω}}_z}{\∂y^2}-\frac{\∂\mathrm{\φ}}{\∂y}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂x}^2}\\ u(\frac{\∂\mathrm{\φ}}{\∂y}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂z}^2}-\frac{\∂\mathrm{\φ}}{\∂z}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂y}^2})+w(\frac{\∂\mathrm{\φ}}{\∂x}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂y}^2}-\frac{\∂\mathrm{\φ}}{\∂y}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂x}^2})\end{array}\right]n_y+$

$\frac{{\mathrm{\ρv}}_{\mathrm{\ω}}{R}_c^2}{|\▿\mathrm{\φ}|}\left[\begin{array}{c}0\\ \frac{\∂\mathrm{\φ}}{\∂y}\frac{{\∂}^2{\mathrm{\ω}}_z}{\∂z^2}-\frac{\∂\mathrm{\φ}}{\∂z}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂y}^2}\\ \frac{\∂\mathrm{\φ}}{\∂z}\frac{{\∂}^2{\mathrm{\ω}}_z}{\∂x^2}-\frac{\∂\mathrm{\φ}}{\∂x}\frac{{\∂}^2{\mathrm{\ω}}_z}{\∂z^2}\\ 0\\ u(\frac{\∂\mathrm{\φ}}{\∂y}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂z}^2}-\frac{\∂\mathrm{\φ}}{\∂z}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂y}^2})+v(\frac{\∂\mathrm{\φ}}{\∂z}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂x}^2}-\frac{\∂\mathrm{\φ}}{\∂x}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂z}^2})\end{array}\right]n_z,---\left(19\right)$

其中，第一个元素用于连续方程、第二至四个元素用于三维直角坐标系三个方向上动量方程的、第五个元素用于能量方程。

本发明的优点：

本发明提出的用于提高模拟可压缩流的旋涡运动的数值方法，不改变原来空间离散的策略，同时可以在数值解中将以向量形式补偿的涡量耗散进行高精度空间离散，更精确地用模拟流场中的旋涡运动。该方法可以方便地在现有数值模拟的源程序中完成，功能的实现成本低，具有较高的工业实用价值，是一种新的可压缩流的旋涡运动的数值模拟技术。

附图说明

图1 涡量耗散的补偿技术的原理图

图中，1:计算网格的界面、2:对流项通量（包含粘性项通量）、3：人工耗散项通量、4：涡量耗散补偿项通量。

图2 涡量耗散补偿项计算步骤的流程图

图3计算空间示意图

图4 数值方法的验证

具体实施方式

以下以一个具体实施例子进一步说明本发明提出的的工作原理。本发明提出的是在流场中补偿人工耗散项产生的涡量耗散，同时在数值解求解构成中将补偿的涡量耗散进行高精度空间离散，更精确地用模拟可压缩流的旋涡运动的数值模拟方法—涡量耗散的补偿技术。

具体实施例子是三维可压缩流流经正方体的流场的数值模拟。图3是计算空间示意图，正方体变长D,整体计算域在正方体下游为10 倍边长，其他空间距离正方体取5倍边长。计算网格采用正六面体结构化计算网格。计算采用有限体积法空间离散，共200x150x150个计算网格。每个计算网格为一个控制体单元，其体积为Ω，控制体外表面由六个面元dS组成，流动变量保存在计算网格中心。

首先将描述可压缩无粘流控制方程，即Euler方程，写成守恒的积分形式,

$\frac{\∂}{\∂t}\underset{\mathrm{\Ω}}{\∫}\stackrel{\→}{W}\mathrm{d\Ω}+\underset{\∂\mathrm{\Ω}}{\∫}({\stackrel{\→}{F}}_c-{\stackrel{\→}{F}}_{\mathrm{\ω}})\mathrm{dS}=\underset{\mathrm{\Ω}}{\∫}\mathrm{Qd\Ω},---\left(20\right)$

式中，守恒变量和对流项向量分别为

$\stackrel{\→}{W}=\left[\begin{array}{c}\mathrm{\ρ}\\ \mathrm{\ρu}\\ \mathrm{\ρv}\\ \mathrm{\ρw}\\ \mathrm{\ρE}\end{array}\right],{\stackrel{\→}{F}}_c\left[\begin{array}{c}\mathrm{\ρV}\\ \mathrm{\ρuV}+n_{x}p\\ \mathrm{\ρvV}+n_{y}p\\ \mathrm{\ρwV}+n_{z}p\\ \mathrm{\ρHV}\end{array}\right],---\left(21\right)$

其中，时间、密度、压力、三维速度向量分别用t、ρ、p、表示，u,v,w是在三维直角坐标系的i,j,k三个方向的分量；抗变速度V是速度向量与网格边界的外法线向量的点乘，即

$V=\stackrel{\→}{V}\·\stackrel{\→}{n}=n_{x}u+n_{y}v+n_{z}w;---\left(22\right)$

总内能E和总焓H分别为

$E=\frac{1}{\mathrm{\γ}-1}\frac{p}{\mathrm{\ρ}}+\frac{{\left|\right|\stackrel{\→}{V}\left|\right|}_2}{2};---\left(23\right)$

$H=E+\frac{p}{\mathrm{\ρ}}.---\left(24\right)$

不计体积力，则源项。涡量补偿项向量由公式（19） 给出。

数值解的求解采用中心差分加上人工耗散的格式，又称JST方法。其中，在计算网格边界上的守恒变量取相邻两个计算网格内的守恒变量的算术平均，同时加入人工耗散项以避免数值解的奇偶振荡现象。控制方程（20）的空间离散形式为

$\frac{\∂\stackrel{\→}{W}}{\∂t}=-\frac{1}{\mathrm{\Ω}}\underset{m=1}{\overset{6}{\mathrm{\Σ}}}[{\stackrel{\→}{F}}_c\left({\stackrel{\→}{W}}_m\right)\mathrm{\Δ}{S}_m-{\stackrel{\→}{D}}_m-{\stackrel{\→}{F}}_m{\mathrm{\ΔS}}_m],---\left(25\right)$

式中变量的下标m表示是在计算网格边界上。以三维直角坐标系中的i-方向为例，在计算网格i和i+1之间的边界上，即i+1/2的位置上，

${\stackrel{\→}{W}}_{i+1/2}=\frac{1}{2}({\stackrel{\→}{W}}_i+{\stackrel{\→}{W}}_{i+1});---\left(26\right)$

带入公式（），求得对流项通量。人工粘性项为

${\stackrel{\→}{D}}_{i+1/2}={\mathrm{\Λ}}_{i+1/2}[{\mathrm{\ϵ}}_{i+1/2}^2({\stackrel{\→}{W}}_{i+1}-{\stackrel{\→}{W}}_i)-{\mathrm{\ϵ}}_{i+1/2}^4({\stackrel{\→}{W}}_{i+2}-3{\stackrel{\→}{W}}_{i+1}+3{\stackrel{\→}{W}}_i-{\stackrel{\→}{W}}_{i-1})],---\left(27\right)$

其中，

${\mathrm{\ϵ}}_{i+1/2}^2=k_2\max (v_{i+1},v_i)$ ；

${\mathrm{\ϵ}}_{i+1/2}^4=\max [0,(k_4-{\mathrm{\ϵ}}_{i+1/2}^2)],---\left(28\right)$

其中，系数 k₂=0.5； k₄=1/96。一个压力传感器可以在流场的大梯度区域或平滑区域启动或关闭人工粘性项的1阶或2阶项。人工粘性项需要用计算网格的谱半径Λ_i+1/2自适应调整，以保证计算精度。

在计算网格边界m上的涡量补偿项向量的计算需要利用辅助网格。在原来计算网格的各个边界为中心，构造另外一套计算网格。辅助网格的中心是原计算网格的边界中心。辅助网格的建立是为了是利用格林公式来计算涡量补偿项向量中的一阶和二阶导数。

计算网格边界上的涡量补偿项向量的计算过程：

1.按照图2给出计算抗变涡运动补偿粘性系数的步骤流程

（1）求辅助网格内的涡量表示，在直角坐标系下有，其中的速度一阶导数项采用格林公式获得；

（2）按照公式（7）（8）计算是涡量的模φ和涡量的模的梯度；

（3）按照公式（13）求辅助网格的补偿涡量的特征半径R_c；

（4）按照公式（15）—（17）求方向系数向量，；

（5）按照公式（14）计算补偿涡量的半径向量；

（6）按照公式（12）计算涡运动补偿粘性系数向量

（7）按照公式（18）计算网格界面上的抗变涡运动补偿粘性系数v_ω。

2.按照公式（19）计算涡量耗散补偿项。

时间离散采用混合的五步龙格-库塔 (Hybrid Multi-stage Runge-Kutta) 显式时间积分方法，利用当地时间和残差光顺方法加速收敛。表达式为

${\stackrel{\→}{W}}_i^{\left(0\right)}={\stackrel{\→}{W}}_i^n$

$\begin{array}{c}{\stackrel{\→}{W}}_i^{\left(1\right)}={\stackrel{\→}{W}}_i^{\left(0\right)}-{\mathrm{\α}}_1\frac{\mathrm{\Δt}}{\mathrm{\Ω}}{({\stackrel{\→}{R}}_c^{\left(0\right)}-{\stackrel{\→}{R}}_d^{\left(0\right)})}_i\\ {\stackrel{\→}{W}}_i^{\left(2\right)}={\stackrel{\→}{W}}_i^{\left(0\right)}-{\mathrm{\α}}_2\frac{\mathrm{\Δt}}{\mathrm{\Ω}}{({\stackrel{\→}{R}}_c^{\left(1\right)}-{\stackrel{\→}{R}}_d^{\left(0\right)})}_i\\ {\stackrel{\→}{W}}_i^{\left(3\right)}={\stackrel{\→}{W}}_i^{\left(0\right)}-{\mathrm{\α}}_3\frac{\mathrm{\Δt}}{\mathrm{\Ω}}{({\stackrel{\→}{R}}_c^{\left(2\right)}-{\stackrel{\→}{R}}_d^{\left(\mathrm{2,0}\right)})}_i\\ {\stackrel{\→}{W}}_i^{\left(4\right)}={\stackrel{\→}{W}}_i^{\left(0\right)}-{\mathrm{\α}}_4\frac{\mathrm{\Δt}}{\mathrm{\Ω}}{({\stackrel{\→}{R}}_c^{\left(3\right)}-{\stackrel{\→}{R}}_d^{\left(\mathrm{2,0}\right)})}_i\\ {\stackrel{\→}{W}}_i^{\left(5\right)}={\stackrel{\→}{W}}_i^{\left(0\right)}-{\mathrm{\α}}_5\frac{\mathrm{\Δt}}{\mathrm{\Ω}}{({\stackrel{\→}{R}}_c^{\left(4\right)}-{\stackrel{\→}{R}}_d^{\left(\mathrm{4,2}\right)})}_i\end{array},---\left(29\right)$

式中，变量的上标表示积分步骤数；为公式（25）中的右手项，而

${\stackrel{\→}{R}}_d^{\left(\mathrm{2,0}\right)}={\mathrm{\β}}_3{\stackrel{\→}{R}}_d^{\left(2\right)}+(1-{\mathrm{\β}}_3){\stackrel{\→}{R}}_d^{\left(0\right)};---\left(30\right)$

${\stackrel{\→}{R}}_d^{\left(\mathrm{4,2}\right)}={\mathrm{\β}}_5{\stackrel{\→}{R}}_d^{\left(4\right)}+(1-{\mathrm{\β}}_5){\stackrel{\→}{R}}_d^{\left(\mathrm{2,0}\right)}.---\left(31\right)$

每一步的系数和CFL数σ在表中列出

远场边界取特征边界条件，对于壁面，采用粘性流的壁面无滑移条件，表示为

$\stackrel{\→}{V}=0.---\left(32\right)$

在一定的来流雷诺数下，流体流经非流线型物体，如正方形、圆形，在其后方的尾迹中有固定的流动模式，表现形式是沿着流动方向出现成对的、旋转方向相反的漩涡，又称卡门涡街（Karman vortex street）。尾迹中的流动是以旋涡为主的运动，流动过程的数值模拟需要对漩涡的弱不连续界面进行捕捉。

具体实施案例中，流体的物性取理想可压缩空气，来流速度取马赫数0.3。图4 给出了本发明给出的数值方法的验证的结果，具体是沿流动方向上正方体中心剖面上的、正方体后方近处的涡量场的等值线图。等值线的细密表示涡量的精度水平较高。图4（a）是作为参考的、采用湍流的大涡模拟技术获得的结果，其细节是公知的。图4（b）采用的是本发明给出的涡量补偿技术在Euler方程下获得的结果。其计算没有考虑流体粘性和湍流的作用，但是在固体壁面采用了粘性流体的壁面无滑移条件。结果表明，计算结果比LES的结果能获得更加细密的涡量场的等值线，表明本发明提出涡量补偿技术的能更加准确地模拟可压缩流的旋流场的运动。

Claims (7)

1.一种模拟可压缩旋流场的数值方法，其特征在于在计算网格的界面加入涡量耗散补偿项，其等于流体密度ρ乘以单位质量涡量耗散补偿项向量，的表达形式是

${\stackrel{\→}{f}}_{\mathrm{\ω}}={\stackrel{\→}{n}}_{\mathrm{\ω}}\×\left(v_{\mathrm{\ω}}\left({\▿}^2\stackrel{\→}{\mathrm{\ω}}\right)R_c^2\right)$ ；

其中，代表涡量的模的梯度最大变化的方向、v_ω代表涡运动补偿粘性系数，与流体的运动粘度有相同的量纲，是一个标量、R_c代表补偿涡量的特征半径；涡量用表示，在直角坐标系下有；符号×代表叉乘运算；是拉普拉斯算子。

2.根据权利要求1所述的一种可压缩旋流场的数值模拟方法，其特征在于所述的涡量的模的梯度最大变化的方向表示为，其中，φ是涡量的模、是涡量的模的梯度、是涡量的模的梯度的模，分别定义为

$\mathrm{\φ}=\left|\stackrel{\→}{\mathrm{\ω}}\right|=\sqrt{{\mathrm{\ω}}_x^2+{\mathrm{\ω}}_y^2+{\mathrm{\ω}}_z^2}$ 、 $\▿\mathrm{\φ}=\frac{\∂\mathrm{\φ}}{\∂x}i+\frac{\∂\mathrm{\φ}}{\∂y}j+\frac{\∂\mathrm{\φ}}{\∂z}k$ 、 $|\▿\mathrm{\φ}|=\sqrt{{\left(\frac{\∂\mathrm{\φ}}{\∂x}\right)}^2+{\left(\frac{\∂\mathrm{\φ}}{\∂y}\right)}^2+{\left(\frac{\∂\mathrm{\φ}}{\∂z}\right)}^2}$ 。

3.根据权利要求1所述的一种可压缩旋流场的数值模拟方法，其特征在于所述的补偿涡量的特征半径R_c在二维和三维空间分别定义为和，其中Ω在二维和三维空间分别代表计算网格面积和体积。

4.根据权利要求1所述的一种可压缩旋流场的数值模拟方法，其特征在于所述的抗变涡运动补偿粘性系数v_ω定义为涡运动补偿粘性系数向量与计算网格界面的法线方向向量的点乘，即，其中，符号·代表点乘计算。

5.根据权利要求4所述的涡运动补偿粘性系数向量，其特征在于，定义为补偿涡量的半径向量与涡量的模值比，即。

6.根据权利要求5所述的补偿涡量的半径向量，其特征在于，定义成补偿涡量的特征半径R_c与涡量在最大作用的方向系数向量的乘积，即，其中，

$k^I=1+\max \left[\right|\frac{{\mathrm{\ω}}_y}{{\mathrm{\ω}}_x}|,|\frac{{\mathrm{\ω}}_z}{{\mathrm{\ω}}_x}\left|\right],---\left(15\right)$ ， $k^J=1+\max \left[\right|\frac{{\mathrm{\ω}}_x}{{\mathrm{\ω}}_y}|,|\frac{{\mathrm{\ω}}_z}{{\mathrm{\ω}}_y}\left|\right],---\left(16\right)$ ， $k^K=1+\max \left[\right|\frac{{\mathrm{\ω}}_x}{{\mathrm{\ω}}_z}|,|\frac{{\mathrm{\ω}}_y}{{\mathrm{\ω}}_z}\left|\right].---\left(17\right)$ 。

7.一种模拟可压缩旋流场的数值方法，其特征在于在计算网格的界面加入涡量耗散补偿项，其在三维直角坐标系下展开形式表示为

${\stackrel{\→}{F}}_{\mathrm{\ω}}=\frac{{\mathrm{\ρv}}_{\mathrm{\ω}}{R}_c^2}{|\▿\mathrm{\φ}|}\left[\begin{array}{c}0\\ 0\\ \frac{\∂\mathrm{\φ}}{\∂z}\frac{{\∂}^2{\mathrm{\ω}}_z}{\∂x^2}-\frac{\∂\mathrm{\φ}}{\∂x}\frac{{\∂}^2{\mathrm{\ω}}_z}{\∂z^2}\\ \frac{\∂\mathrm{\φ}}{\∂x}\frac{{\∂}^2{\mathrm{\ω}}_z}{\∂y^2}-\frac{\∂\mathrm{\φ}}{\∂y}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂x}^2}\\ v(\frac{\∂\mathrm{\φ}}{\∂z}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂x}^2}-\frac{\∂\mathrm{\φ}}{\∂x}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂z}^2})+w(\frac{\∂\mathrm{\φ}}{\∂x}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂y}^2}-\frac{\∂\mathrm{\φ}}{\∂y}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂x}^2})\end{array}\right]n_x+$

$\frac{{\mathrm{\ρv}}_{\mathrm{\ω}}{R}_c^2}{|\▿\mathrm{\φ}|}\left[\begin{array}{c}0\\ \frac{\∂\mathrm{\φ}}{\∂\mathrm{\∂}}\frac{{\∂}^2{\mathrm{\ω}}_z}{\∂z^2}-\frac{\∂\mathrm{\φ}}{\∂z}\frac{{\∂}^2{\mathrm{\ω}}_z}{\∂y^2}\\ 0\\ \frac{\∂\mathrm{\φ}}{\∂x}\frac{{\∂}^2{\mathrm{\ω}}_z}{\∂y^2}-\frac{\∂\mathrm{\φ}}{\∂y}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂x}^2}\\ u(\frac{\∂\mathrm{\φ}}{\∂y}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂z}^2}-\frac{\∂\mathrm{\φ}}{\∂z}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂y}^2})+w(\frac{\∂\mathrm{\φ}}{\∂x}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂y}^2}-\frac{\∂\mathrm{\φ}}{\∂y}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂x}^2})\end{array}\right]n_y+$

$\frac{{\mathrm{\ρv}}_{\mathrm{\ω}}{R}_c^2}{|\▿\mathrm{\φ}|}\left[\begin{array}{c}0\\ \frac{\∂\mathrm{\φ}}{\∂y}\frac{{\∂}^2{\mathrm{\ω}}_z}{\∂z^2}-\frac{\∂\mathrm{\φ}}{\∂z}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂y}^2}\\ \frac{\∂\mathrm{\φ}}{\∂z}\frac{{\∂}^2{\mathrm{\ω}}_z}{\∂x^2}-\frac{\∂\mathrm{\φ}}{\∂x}\frac{{\∂}^2{\mathrm{\ω}}_z}{\∂z^2}\\ 0\\ u(\frac{\∂\mathrm{\φ}}{\∂y}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂z}^2}-\frac{\∂\mathrm{\φ}}{\∂z}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂y}^2})+v(\frac{\∂\mathrm{\φ}}{\∂z}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂x}^2}-\frac{\∂\mathrm{\φ}}{\∂x}\frac{{\∂}^2{\mathrm{\ω}}_z}{{\∂z}^2})\end{array}\right]n_z$

其中，u，v，w是直角坐标系下的速度的三维分量。