CN103914608A - 不可压缩旋流场的数值模拟中使用的涡量保持技术 - Google Patents

Abstract

本发明是一种不可压缩流的旋涡运动的数值模拟方法，被称作涡量保持技术。根据不可压缩流的特点，在动量方程中通过加入两种不同形式的力，以提高一类以旋涡运动为主的流场的数值模拟精度。这两种形式的力分别是涡量在变化梯度方向的螺旋力和涡量在变化梯度方向的粘性耗散力。该方法使计算网格内的涡量在变化梯度方向的螺旋力的积分计算转化为计算网格边界上的上的力的通量计算，可以使其空间离散具有高阶精度的格式；同时动量方程的源项保留涡量在变化梯度方向的粘性耗散力，用来提高数值解的收敛性和稳定性。这两个力采用不同的放大系数，可以进一步保持涡量的精度，更精确地用模拟流场中的旋涡运动。

Description

不可压缩旋流场的数值模拟中使用的涡量保持技术

技术领域

本发明涉及计算流体力学(CFD：Computational Fluid Dynamics)领域中的一种数值方法，具体是一种在不可压缩旋流场的数值模拟中使用的涡量保持技术。

背景技术

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

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

数值耗散对流场的数值模拟结果的最明显的影响体现在对流动变量的间断界面的捕捉。如前所述，激波是强间断界面，对其捕捉必须加入一定的人工数值耗散。因为激波前后存在熵增，即能量的损失，所以，通过加入人工数值耗散捕捉激波具有合理的物理意义。但是，过大的数值耗散会使数值解获得的激波界面变得模糊，降低了数值解对激波强度和空间位置的预测精度。流场中存在另一类流动不连续现象，即接触不连续(contact discontinuity)。相对激波而言，接触不连续的特点是弱不连续，跨过不连续界面，压力和法向速度是连续的。工程实践中，这种类型的流动现象大量存在。例如，流场中钝体后方的涡脱落、自由剪切流的流动等等。数值模拟中对于这种接触不连续的捕捉更加困难，因为数值方法中的数值耗散即使很小也会使弱不连续界面变得模糊，降低数值解对流场的预测精度，这也是旋流场(Vortex-dominated Flows)的数值模拟技术成为CFD领域的重大挑战的原因。

为了提高旋流场的数值模拟的精度，一种方法是加密计算网格，在更加细小的空间尺度内求解流体控制方程。加密计算网格首先会使计算量加大，增加计算成本。此外，数值计算的误差随着计算网格的增加会不断积累，在一定程度上造成相反的效果。另一种方法是在流场中采用物理模型来增加流场中描述旋流流动的变量-涡量(vorticity)的强度。例如在流场中加入点涡模型，可以人为地增加涡量；或者在流场局部直接求解涡量方程，以减小涡量的输运过程中的耗散。但是，这些方法在应用上仍受到一定限制。点涡模型是在预先明确旋涡发生位置的前提下才能使用，仅适合一些简单的流动现象。除了二维不可压缩正压流场，涡量方程比与欲求解的Euler，Navier-Stokes方程更为复杂。

二十世纪初期，提出了一种提高不可压缩(incompressible)旋流场的求解精度的数值方法，涡量限制法(Vorticity Confinement)。该方法的原理是在流体控制方程的动量方程中的源项位置，加入一个涡量形式的体积力项，从数值耗散中将涡量减去，克服数值耗散造成的旋涡场的接触不连续界面的模糊，从而更精确地捕捉旋涡结构，实现提高旋涡场的计算精度的目的。涡量限制法的具体内容如下：首先写出不可压缩、无粘流的控制方程，包括连续方程和动量方程，分别为

$\▿\·\stackrel{\→}{V}=0,---\left(1\right)$

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

式中，V是速度矢量，包含直角坐标系下的三维i，j，k分量u，v，w；ρ、p、t分别是密度、压力和时间。式中算子表示为符号●代表内积计算。动量方程(2)的右边是体积力项形式的源项按照涡量限制法中的定义

$\stackrel{\→}{B}={\stackrel{\→}{n}}_c\×\stackrel{\→}{\mathrm{\ω}},---\left(3\right)$

式中，代表涡量，在直角坐标系下有按照涡量的定义，

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

式中符号×代表叉乘运算，代表涡量梯度变化的方向，即

${\stackrel{\→}{n}}_c=\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|},---\left(4\right)$

其中，φ是涡量的模；是涡量的模的梯度的模，即

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

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

动量方程(2)中的源项需要公式(3)中的乘以一个比例因子ε，则

$\stackrel{\→}{B}=\mathrm{\ϵ}({\stackrel{\→}{n}}_c\×\stackrel{\→}{\mathrm{\ω}}),---\left(7\right)$

其物理意义是指向涡量中心的力。Steinhoff的原始的涡量限制法中给出ε为常数，为0.01-0.05，其作用是用来调整涡量限制的大小。该方法形式简单，不需要加密计算网格，在旋流场的数值模拟中得到广泛应用。对该方法在随后的改进主要集中在参数ε的表达。

涡量限制法是依靠在流场中加入限定的涡量以抵消数值耗散的原理来模拟旋流场的流动状态。尽管该方法在捕捉不可压缩流场中的接触不连续的方面已经取得了明显的改进效果，但存在以下缺陷：

1.对加入的涡量调整完全靠系数ε；

2.加入的涡量的空间离散精度是被限定的，无法进一步提高；

3.源项对动量方程的数值解的收敛的稳定性影响是不确定的。

为了更加精确地模拟不可压缩流的旋涡运动，需要进一步改进对于流场中接触不连续的捕捉机制。一种途径是根据不可压缩流的特点，改进在流场中通过加入涡量来抵消数值耗散的内在工作机理，通过保持涡量的精度，更精确地用模拟流场中的旋涡运动，形成一种新的不可压缩流的旋涡运动的数值模拟技术，涡量保持技术(Vorticity Refinement)。该技术可以使涡量的空间离散具有高阶精度的格式，同时该可以利用源项增进收敛进程的稳定性。

发明内容

本发明涉及计算流体力学领域中一种模拟不可压缩旋流场的数值方法，具体是根据不可压缩流的特点，在流场中通过加入两种不同形式的涡量力来抵消数值耗散的数值方法。该方法可以使加入的涡量的空间离散具有高阶精度的格式，同时还可以利用源项增进数值解的收敛进程的稳定性。通过保持涡量的精度，更精确地用模拟流场中的旋涡运动，是一种新的不可压缩流的旋涡运动的数值模拟技术-涡量保持技术(Vorticity Refinement)。

本发明采用的技术方案：

如公式(2)表示的，涡量限制法中在不可压缩流动的动量方程的等号右边加入了一个体积力项，它代表涡量在其变化梯度相反的方向的受力，如公式(7)所示。实际上，如果将公式(4)带入公式示(7)并运用算子的运算规律可得

$\stackrel{\→}{B}=\mathrm{\ϵ}({\stackrel{\→}{n}}_c\×\stackrel{\→}{\mathrm{\ω}})=\mathrm{\ϵ}\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\×\stackrel{\→}{\mathrm{\ω}}=\mathrm{\ϵ}\▿\×\left(\frac{\mathrm{\φ}}{|\▿\mathrm{\φ}|}\stackrel{\→}{\mathrm{\ω}}\right)-\mathrm{\ϵ}\frac{\mathrm{\φ}}{|\▿\mathrm{\φ}|}(\▿\×\stackrel{\→}{\mathrm{\ω}}).---\left(8\right)$

在上式中再次运用算子的运算规律可得

$\▿\×\stackrel{\→}{\mathrm{\ω}}=\▿\×(\▿\×\stackrel{\→}{V})=\▿(\▿\·\stackrel{\→}{V})-{\▿}^2\stackrel{\→}{V},---\left(9\right)$

式中是拉普拉斯算子，定义为又根据公式(1)不可压缩流的连续方程，将公式(9)代入公式(8)，可得

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

对于不可压缩流，公式(10)与公式(7)完全等价。现将公式(10)中的重新写成如下形式，即

$\stackrel{\→}{B}={\stackrel{\→}{B}}_1+{\stackrel{\→}{B}}_2,---\left(11\right)$

其中，

${\stackrel{\→}{B}}_1={\mathrm{\ϵ}}_1\▿\×\left(\frac{\mathrm{\φ}}{|\▿\mathrm{\φ}|}\stackrel{\→}{\mathrm{\ω}}\right),---\left(12\right)$

因为是螺旋度的定以，所以被称作涡量在变化梯度方向的螺旋力；

${\stackrel{\→}{B}}_2={\mathrm{\ϵ}}_2\frac{\mathrm{\φ}}{|\▿\mathrm{\φ}|}\left({\▿}^2\stackrel{\→}{V}\right),---\left(13\right)$

因为拉普拉斯算子的耗散特性，被称作在涡量在变化梯度方向的粘性耗散力。

至此公式(2)中以动量方程的源项形式表示的涡量在其变化梯度相反的方向的受力(涡量限制)限制被分解为两部分，形成两种力的形式。其中涡量在变化梯度方向的螺旋力 与涡量的变化方向和涡量大小有关；而涡量在变化梯度方向的粘性耗散力 仅与涡量的变化梯度方向有关。图1给出两种形式的力的形成原理和关系图。图中表明，由涡量和涡量梯度变化的方向向量构成平面1，被称作涡量作用平面。原始的涡量限制与该平面垂直，并与它的两个分量，涡量在变化梯度方向的螺旋力和涡量变化梯度方向的粘性耗散力共同构成平面2，被称作涡量限制平面。平面1垂直与平面2。和的夹角被称作螺旋角，该夹角的余弦即为涡量限制被分解后形成的两种形式的力，和将在数值方法中起到不同的作用。

公式(12)和(13)中是用了两个参数：ε₁和ε₂，分别作为不同的两个力的放大系数。如果ε₁和ε₂相等，且等于公式(7)中的ε，则成为同向涡量限制(isotropic vorticity  confinement)，其效果等同于公式(7)表示的原始的涡量限制；而ε₁和ε₂不相等时，成为异向涡量限制(anisotropic vorticity confinement)。而不可压缩旋流场的数值模拟精度的提高可以通过异性涡量限制获得。

作为源项形式的体积力对偏微分方程(2)的数值解的收敛过程有重要影响。其作用体现在加入源项后，方程(2)有可能变成具有刚性的(stiff)的性质，致使一类依靠显式时间迭代获得数值解的收敛性降低，同时也降低数值解的稳定性。因而，此类方程的求解为了保证时间精度，必须采用较小的时间步长，结果造成计算速度变慢。衡量方程的刚性的指标是源项的雅各比矩阵的特征值。如果雅各比矩阵的最大特征值和最小特征值的实部的比例小，则说明偏微分方程的刚性小，数值解容易收敛，而且稳定性强。

考虑两种情况：第一，源项分别为即原始的涡量限制法中的体积力，源项的雅各比矩阵为特征值为i等于2或3；第二，源项为即涡量变化梯度方向的粘性耗散力，源项的雅各比矩阵为特征值为在两种形式的体积力的情况下，源项的雅各比矩阵分别表示为

$\stackrel{\‾}{\stackrel{\‾}{S}}=\frac{\∂\stackrel{\→}{B}}{\∂\stackrel{\→}{V}},---\left(14\right)$

$\stackrel{\‾}{\stackrel{\‾}{S}}=\frac{\∂{\stackrel{\→}{B}}_2}{\∂\stackrel{\→}{V}}.---\left(15\right)$

按照雅各比矩阵的特征值的求解方法，求解下式

$|{\mathrm{\λ}}_i^{B}I-\stackrel{\‾}{\stackrel{\‾}{S}}|=0$ 和 $|{\mathrm{\λ}}_i^{B_2}I-{\stackrel{\‾}{\stackrel{\‾}{S}}}_2|=0,---\left(16\right)$

其中I是单位矩阵，即可获得两种源项的特征值，而且有

$\frac{\min \left[\mathrm{Re}\left({\mathrm{\λ}}_i^B\right)\right]}{\max \left[\mathrm{Re}\left({\mathrm{\λ}}_i^B\right)\right]}\≤\frac{\min \left[\mathrm{Re}\left({\mathrm{\λ}}_i^{B_2}\right)\right]}{\max \left[\mathrm{Re}\left({\mathrm{\λ}}_i^{B_2}\right)\right]}.---\left(17\right)$

上式中Re表示取实部运算。表明了以为源项时，偏微分方程的刚性比以为源项时小。因而，在使用一类显式时间积分方法求解不可压缩流的动量方程(2)时会增强收敛性和数值稳定性。所以，利用这个特性，可以将(涡量在变化梯度方向的螺旋力)从等号右边移动到等号左边，等号右边的源项仅保留(涡量变化梯度方向的粘性耗散力)。

如公式(12)所示，涡量在变化梯度方向的螺旋力与涡量的变化方向和涡量大小有关，数值解中需要高精度的空间离散，以提高模拟精度。源项中的被移动等号左边后，可以利用高斯定理，将控制体单元内的体积力形式转变为表面通量的积分形式。有限体积法(FiniteVolume Method)中，某一控制体单元，即某一计算网格的体积为Ω，而其表面积为面积元为dS。根据高斯定理，可将体积积分变换为面积积分，则涡量在变化梯度方向的螺旋力从体积积分变为面积积分的过程可表示为

$\underset{\mathrm{\Ω}}{\∫}{\stackrel{\→}{B}}_1\mathrm{d\Ω}=\underset{\mathrm{\Ω}}{\∫}{\mathrm{\ϵ}}_1\▿\×\left(\frac{\mathrm{\φ}}{|\▿\mathrm{\φ}|}\stackrel{\→}{\mathrm{\ω}}\right)\mathrm{d\Ω}$ $\left(18\right)$

$=\underset{\∂\mathrm{\Ω}}{\∫}{\mathrm{\ϵ}}_1\stackrel{\→}{n}\×\left(\frac{\mathrm{\φ}}{|\▿\mathrm{\φ}|}\stackrel{\→}{\mathrm{\ω}}\right)\mathrm{dS},$

其中计算网格边界的单位向量。于是，在计算网格边界上产生了一个由于涡量在变化梯度方向的螺旋力产生的向量，被称作涡量的螺旋力向量

${\stackrel{\→}{F}}_{\mathrm{\ω}}={\mathrm{\ϵ}}_1\stackrel{\→}{n}\×\left(\frac{\mathrm{\φ}}{|\▿\mathrm{\φ}|}\stackrel{\→}{\mathrm{\ω}}\right).---\left(19\right)$

在数值解的求解过程中，可以对该向量在计算网格边界上进行高阶精度的空间离散。由于该向量与涡量的变化方向和涡量大小有关，为进一步提高了旋涡流场的空间模拟精度提供了可能性。例如可以同过流动变量的高阶插值获得计算网格界面上用来计算向量的值。

此外，如果适当增大ε₁、减小ε₂，意味着增大涡量在变化梯度方向的螺旋力、降低涡量变化梯度方向的粘性耗散力，即采用异向涡量限制，利用调节这两个力的内部关系的机制，可以进一步调整旋涡运动中的接触不连续界面的捕捉效果。

总之，本发明提出的提高可压缩旋流场的数值模拟精度的数值方法包括四个技术元素，分别为：

1.将涡量限制法中的源项分解称为两个力，涡量在变化梯度方向的螺旋力和涡量变化梯度方向的粘性耗散力，并将涡量在变化梯度方向的螺旋力移动到不可压缩流的动量方程的等号左变；

2.通过高斯定理将计算网格内涡量在变化梯度方向的螺旋力转化为计算网格边界上的上的力，进行高精度的空间离散，按照通量进行计算；

3.源项保留涡量变化梯度方向的粘性耗散力用来提高数值解的收敛性和稳定性；

4.涡量在变化梯度方向的螺旋力和涡量变化梯度方向的粘性耗散力采用不同的放大系数，即采用异向的涡量限制。

该数值方法被称为可压缩旋流场的涡量保持技术(Vorticity Refinement)。

本发明的优点：

本发明提出的提高不可压缩旋流场的数值模拟精度的数值方法是根据不可压缩流的特点，通过加入两种不同形式的力来改进在流场中抵消数值耗散的内在工作机制，通过保持涡量的精度，更精确地模拟流场中的旋涡运动。该方法使计算网格内的涡量在变化梯度方向的螺旋力转化为计算网格边界上的上的力，可以使其空间离散具有高阶精度的格式；同时源项保留涡量变化梯度方向的粘性耗散力，用来提高数值解的收敛性和稳定性。这两个力采用不同的放大系数，可以进一步保持涡量的精度，更精确地用模拟流场中的旋涡运动。

附图说明

图1涡量在变化梯度方向的螺旋力和涡量变化梯度方向的粘性耗散力形成原理和关系图

图中，1：涡量作用平面、2：螺旋角、3涡量限制平面。

图2三角机翼的外形和计算网格。

图3数值方法的验证，其中包括：

图3(a)源项对收敛性作用的比较；

图3(b)异向涡量限制效果的比较。

具体实施方式

以下以一个具体实施方案进一步说明本发明提出的，，这种数值模拟技术可以用Fortran90计算机高级程序语言实现，并通过计算机运行来模拟三维三角机翼在高攻角飞行状态下的周围流场。

具体实施方案中，采用大展翼比的三角机翼在高攻角下的不可压缩流场。大量实验表明三角机翼在高攻角下飞行，其周围的流场是以旋涡运动为主的流场。三角翼弦长为c、翼展为b，展翼比为7。计算网格采用多块结构化适体六面体网格，机翼周围采用O-型网格，围绕机翼剖面方向192个网格；法线方向68个网格；翼展方向96个网格。图2表示了三角机翼的几何形状和围绕三角机翼生成的计算网格。数值模拟的空间离散采用有限体积法(FiniteVolume Method)，每个计算网格成为控制单元，流动变量在控制单元的中心。

数值模拟的飞行条件是：

来流速度：57m/s

大气压力：100kPa

大气密度：1.2kg/m³

攻角：20°

不可压缩无粘流的动量方程，即方程(2)中不计粘性项本发明中要求加上在计算网格边界上产生了一个由于涡量在变化梯度方向的螺旋力产生的向量意味着在动量方程左边的对流项后面加上这个向量，在等号右边，即源项的位置加上涡量变化梯度方向的粘性耗散力按照三维守恒的积分形式重新写出动量方程，

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

其中，守恒变量对流项向量和涡量的螺旋力向量(见公式(19))分别为

$\stackrel{\→}{V}=\left[\begin{array}{c}u\\ v\\ w\end{array}\right];$ ${\stackrel{\→}{F}}_c=\left[\begin{array}{c}\mathrm{uV}+n_{x}p/\mathrm{\ρ}\\ \mathrm{vV}+n_{y}p/\mathrm{\ρ}\\ \mathrm{wV}+n_{z}p/\mathrm{\ρ}\end{array}\right];$ ${\stackrel{\→}{F}}_{\mathrm{\ω}}={\mathrm{\ϵ}}_1\frac{\mathrm{\φ}}{|\▿\mathrm{\φ}|}\left[\begin{array}{c}{\mathrm{\ω}}_z{n}_y-{\mathrm{\ω}}_y{n}_z\\ {\mathrm{\ω}}_x{n}_z-{\mathrm{\ω}}_z{n}_x\\ {\mathrm{\ω}}_y{n}_x-{\mathrm{\ω}}_x{n}_y\end{array}\right].---\left(21\right)$

上式中，密度ρ为一固定值；V代表速度不变量，是网格边界的外法线向量与速度向量的点乘，即

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

可以将对流项向量和涡量的螺旋力向量写成展开形式，

${\stackrel{\→}{F}}_c=\left[\begin{array}{c}{u}^2+p/\mathrm{\ρ}\\ \mathrm{uv}\\ \mathrm{uw}\end{array}\right]n_x+\left[\begin{array}{c}\mathrm{uv}\\ v^2+p/\mathrm{\ρ}\\ \mathrm{vw}\end{array}\right]n_y+\left[\begin{array}{c}\mathrm{uw}\\ \mathrm{vw}\\ w^2+p/\mathrm{\ρ}\end{array}\right]n_z.---\left(24\right)$

对流项经压力分离后，可以写成不计压力的对流项向量和压力向量之和，显然有，

${\stackrel{\→}{F}}_c={\stackrel{\→}{F}}_{\mathrm{cV}}+{\stackrel{\→}{F}}_{\mathrm{cp}},---\left(25\right)$

其中，

${\stackrel{\→}{F}}_{\mathrm{cV}}=\left[\begin{array}{c}{u}^2\\ \mathrm{uv}\\ \mathrm{uw}\end{array}\right]n_x+\left[\begin{array}{c}\mathrm{uv}\\ v^2\\ \mathrm{vw}\end{array}\right]n_y+\left[\begin{array}{c}\mathrm{uw}\\ \mathrm{vw}\\ w^2\end{array}\right]n_z;---\left(26\right)$

${\stackrel{\→}{F}}_{\mathrm{cp}}=\left[\begin{array}{c}p/\mathrm{\ρ}\\ 0\\ 0\end{array}\right]n_x+\left[\begin{array}{c}0\\ p/\mathrm{\ρ}\\ 0\end{array}\right]n_y+\left[\begin{array}{c}0\\ 0\\ p/\mathrm{\ρ}\end{array}\right]n_z.---\left(27\right)$

进一步可以写出本发明提出的涡量在变化梯度方向的螺旋力向量 和源项形式的涡量 变化梯度方向的粘性耗散力 的具体形式分别为，

${\stackrel{\→}{F}}_{\mathrm{\ω}}={\mathrm{\ϵ}}_1\frac{\mathrm{\φ}}{|\▿\mathrm{\φ}|}\left[\begin{array}{c}0\\ \frac{\∂u}{\∂y}-\frac{\∂v}{\∂x}\\ \frac{\∂u}{\∂z}-\frac{\∂w}{\∂x}\end{array}\right]n_x+{\mathrm{\ϵ}}_1\frac{\mathrm{\φ}}{|\▿\mathrm{\φ}|}\left[\begin{array}{c}\frac{\∂v}{\∂x}-\frac{\∂u}{\∂y}\\ 0\\ \frac{\∂v}{\∂z}-\frac{\∂w}{\∂v}\end{array}\right]n_y+{\mathrm{\ϵ}}_1\frac{\mathrm{\φ}}{|\▿\mathrm{\φ}|}\left[\begin{array}{c}\frac{\∂w}{\∂x}-\frac{\∂u}{z}\\ \frac{\∂w}{\∂y}-\frac{\∂v}{\∂z}\\ 0\end{array}\right]n_z;---\left(28\right)$

${\stackrel{\→}{B}}_2={\mathrm{\ϵ}}_2\frac{\mathrm{\φ}}{|\▿\mathrm{\φ}|}\left[\begin{array}{c}\frac{{\∂}^{2}u}{\∂x^2}+\frac{{\∂}^{2}u}{\∂y^2}+\frac{{\∂}^{2}u}{z^2}\\ \frac{{\∂}^{2}u}{\∂x^2}+\frac{{\∂}^{2}v}{\∂y^2}+\frac{{\∂}^{2}v}{\∂z^2}\\ \frac{{\∂}^{2}w}{\∂x^2}+\frac{{\∂}^{2}w}{\∂y^2}+\frac{{\∂}^{2}w}{\∂z^2}\end{array}\right],---\left(29\right)$

其中，涡量的模φ和涡量的模的梯度的模由公式(5)和(6)给出。公式(28)和(29)中涉及的一阶导数和二阶导数项均可利用格林公式求取或者直接按照泰勒展开求取，其过程是公知的，不再叙述。

以下是在同位网格(Collocated Grid)中运用基于压力的方法(Pressure-Based Method)求解包含涡量在变化梯度方向的螺旋力向量和涡量变化梯度方向的粘性耗散力的不可压缩流控制方程的过程。

首先将公式(20)写成离散形式，

$\left(\frac{{\stackrel{\→}{V}}^{t+\mathrm{\Δt}}-\stackrel{\→}{V}}{\mathrm{\Δt}}\right){\mathrm{\Ω}}_{i,j,k}+\underset{m}{\overset{N_F}{\mathrm{\Σ}}}{({\stackrel{\→}{F}}_{\mathrm{cV}}+{\stackrel{\→}{F}}_{\mathrm{cp}}-{\stackrel{\→}{F}}_{\mathrm{\ω}})}_m\mathrm{\Δ}{S}_m={\left({\stackrel{\→}{B}}_2\mathrm{\Ω}\right)}_{i,j,k},---\left(31\right)$

其中，Ω_i，j，k代表空间任一离散控制体体积；Δt代表积分时间步长；代表流场更新后速度；ΔS_m代表控制体任一表面面积；N_F代表控制体表面个数，三维情况下为6。重新整理右手项，有

$\stackrel{\→}{R}\left(\stackrel{\→}{V}\right)=\underset{m}{\overset{N_F}{\mathrm{\Σ}}}{({\stackrel{\→}{F}}_{\mathrm{cV}}-{\stackrel{\→}{F}}_{\mathrm{\ω}})}_m{\mathrm{\ΔS}}_m-{\left({\stackrel{\→}{B}}_2\mathrm{\Ω}\right)}_{i,j,k}+\underset{m}{\overset{N_F}{\mathrm{\Σ}}}{\left({\stackrel{\→}{F}}_{\mathrm{cp}}\right)}_m{\mathrm{\ΔS}}_m;---\left(32\right)$

而使用一个假设的初始速度计算一个不计压力项的右手项即

${\stackrel{\→}{R}}_c\left({\stackrel{\→}{V}}^{*}\right)=\underset{m}{\overset{N_F}{\mathrm{\Σ}}}{[{\stackrel{\→}{F}}_{\mathrm{cV}}\left({\stackrel{\→}{V}}^{*}\right)-{\stackrel{\→}{F}}_{\mathrm{\ω}}{\left({\stackrel{\→}{V}}^{*}\right)}_m]\mathrm{\ΔS}}_m-{\stackrel{\→}{B}}_2\left({\stackrel{\→}{V}}^{*}\right){\mathrm{\Ω}}_{i,j,k},---\left(33\right)$

显然有，

$\stackrel{\→}{R}\left(\stackrel{\→}{V}\right)={\stackrel{\→}{R}}_c\left({\stackrel{\→}{V}}^{*}\right)+\underset{m}{\overset{N_F}{\mathrm{\Σ}}}{\left({\stackrel{\→}{F}}_{\mathrm{cp}}\right)}_m\mathrm{\Δ}{S}_m.---\left(34\right)$

不计压力项时候，可以从公式(31)获得这个个初始速度场

${\stackrel{\→}{V}}^{*}=\stackrel{\→}{V}-\frac{\mathrm{\Δt}}{{\mathrm{\Ω}}_{i,j,k}}{\stackrel{\→}{R}}_c\left(\stackrel{\→}{V}\right),---\left(35\right)$

其中需要在计算网格边界m上的不计压力的对流项向量涡量在变化梯度方向的螺旋力向量以及在计算网格内的涡量变化梯度方向的粘性耗散力再带入公式(33)以获得不计压力的右手项

在按照公式(28)计算网格边界m处的涡量在变化梯度方向的螺旋力向量的时候，对于其中的一阶导数项采用高精度离散方法，例如在i-方向的三阶空间离散表示为，

$\frac{\∂u}{\∂x}=\frac{2u_{i+1}+3u_i-6u_{i-1}+u_{i-2}}{6\mathrm{\Δx}}.---\left(36\right)$

网格中心点的初始速度在计算网格边界上利用公知的MUSCL方法进行二阶精度插值，可获得在网格边界m上的速度因而，根据公式(35)不计压力项的情况下，在边界m处的速度表示为

${\stackrel{\→}{\stackrel{\‾}{V}}}_m^{*}={\stackrel{\→}{V}}_m^{*}-\frac{\mathrm{\Δt}}{{{\mathrm{\Ω}}^{*}}_m}{\stackrel{\→}{R}}_c\left({\stackrel{\→}{V}}_m^{*}\right),---\left(37\right)$

其中，是包含网格边界m的辅助网格的体积。在边界m处，考虑压力项的速度为

${\stackrel{\→}{V}}_m={\stackrel{\→}{V}}_m^{*}-\frac{\mathrm{\Δt}}{{\mathrm{\Ω}}_{i,j,k}}{\stackrel{\→}{R}}_c\left({\stackrel{\→}{V}}_m^{*}\right)-\frac{\mathrm{\Δt}}{{\mathrm{\Ω}}_{i,j,k}}{\stackrel{\→}{R}}_{\mathrm{pm}},---\left(38\right)$

其中，是公式(32)在网格边界m处的值，任然用MUSCL方法进行二阶插值获得。整理公式(37)和(38)，有

${\stackrel{\→}{V}}_m={\stackrel{\→}{\stackrel{\‾}{V}}}_m^{*}-\frac{\mathrm{\Δt}}{{\mathrm{\Ω}}_m}{\stackrel{\→}{R}}_{\mathrm{pm}}.---\left(39\right)$

考虑压力项的速度应该满足不可压缩流的连续方程(1)，即将连续方程(1)写成离散形式，

$\underset{m}{\overset{N_F}{\mathrm{\Σ}}}{\stackrel{\→}{V}}_m{\mathrm{\ΔS}}_m=0.---\left(40\right)$

将公式(39)带入(40)，有

$\underset{m}{\overset{N_F}{\mathrm{\Σ}}}({\stackrel{\→}{\stackrel{\‾}{V}}}_m^{*}-\frac{\mathrm{\Δt}}{{\mathrm{\Ω}}_{i,j,k}}{\stackrel{\→}{R}}_{\mathrm{pm}}){\mathrm{\ΔS}}_m=0.---\left(41\right)$

求解方程(41)，即可获得t+Δt时刻，流场的计算网格边界m处的速度和压力值p_m，这两个值可以通过集合平均处理，转化为计算网格中心处的值。

图3给出数值方法的验证。其中图3(a)表示源项对收敛性作用的比较，图中表明了用两种方法获得的残差收敛曲线。采用本发明提出的异向涡量限制，ε₁＝0.07、ε₂＝0.03比采用原始的涡量限制法，ε＝0.05时，收敛速度加快，原因在与源项的形式加强了数值解的收敛性和稳定性。图3(b)表明异向涡量限制效果的比较，图中表明了两种方法计算的三角机翼上部涡量云图。采用本发明提出的异向涡量限制，ε₁＝0.07、ε₂＝0.03比采用同向涡量限制，ε＝0.05时可以捕捉到更强的涡量，原因在于采用了高阶精度离散和较大的放大系数。

Claims (7)

1.一种在不可压缩旋流场的数值模拟中使用的涡量保持技术，其特征在于在不可压缩流的动   量方程中加入两种不同形式的力，分别是涡量在变化梯度方向的螺旋力和涡量在变化梯度方向的粘性耗散力。

2.根据权利要求1所述的一种在不可压缩旋流场的数值模拟中使用的涡量保持技术，其特征 在于所述的涡量在变化梯度方向的螺旋力的形式是

${\stackrel{\→}{B}}_1={\mathrm{\ϵ}}_1\▿\×\left(\frac{\mathrm{\φ}}{|\▿\mathrm{\φ}|}\stackrel{\→}{\mathrm{\ω}}\right)$ ，

其中，代表涡量，在直角坐标系下有；φ是涡量的模；是涡量的模的梯度的模，即；；ε₁是的放大系数。

3.根据权利要求1所述的一种在不可压缩旋流场的数值模拟中使用的涡量保持技术，其特征 在于所述的涡量变化梯度方向的粘性耗散力的形式是。

${\stackrel{\→}{B}}_2={\mathrm{\ϵ}}_2\frac{\mathrm{\φ}}{|\▿\mathrm{\φ}|}\left({\▿}^2\stackrel{\→}{V}\right)$ ，

其中，V是速度矢量，；是拉普拉斯算子；ε₂是的放大系数。

4.根据权利要求2所述涡量在变化梯度方向的螺旋力的形式通过高斯定理，将计算网格单元内的力的积分形式转变为计算网格表面通量的积分形式，并且保留在动量方程等号的左边，采用二阶以上空间离散精度进行空间离散。

5.根据权利要求2所述的涡量在变化梯度方向的粘性耗散力的形式保留在动量方程右边的源项位置。

6.根据权利要求2所述的涡量在变化梯度方向的螺旋力的放大系数ε₁，其特征在于ε₁的范围为0.05<ε₁<0.07。

7.根据权利要求2所述的涡量在变化梯度方向的粘性耗散力的放大系数ε₂，其特征在于ε₂的范围为0.03<ε₂<0.05。