2.3　CFD理论

最早的非定常不可压缩流体Navier-Stokes方程是以原始变量的形式来处理流动问题。一方面，原始变量公式提供了一种独特的视角认识质量守恒定律与边界条件。类似的算法可以处理二维和三维的流动问题。在必要情况下可以直接使用直观的湍流模型。另一方面，不可压缩流动的原始变量公式缺乏一个物理意义明确的压力方程。通常采用Poisson的压力方程［30］，对动量方程进行发散。也可以采用另外一种解决办法，添加一个与速度场直接相关的压力修正项［31］，这个压力修正项最后会出现在连续性方程中，左边乘以一个线性算子使压力方程能够顺利求解。算子的选择是根据不同的压力校正方法来确定的。只有当矩阵对角线数值被保留在算子上时才能得到一个人为的压缩性公式［32］。如果采用不包含矩阵的共轭梯度法来处理原始变量方程组，就无法顺利求解压力方程［33］。

原始变量公式在经过反复离散化后存在速度-压力解耦合的问题。交错网格法［34］虽然可以进行补救，但大大增加了程序的复杂性。如果对速度和压力进行不同程度的插值（尤其是对有限元分析［35］ ）或者在连续性方程中引入高阶压力耗散项［36］，非交错网格的压力振荡则会被抑制。

在本文中，非定常不可压缩Navier-Stokes方程的原始变量形式用压力修正方法进行了处理。其中，压力修正采用的是不包含矩阵的求解方法。处理层流问题的相应方程组已由Chaviaropoulos［37］提出，它适用于内部、外部二维流场的计算。曲线坐标系上，控制方程以非定常的形式在结构化网格中进行离散。速度分量在网格节点上计算，而压力在网格中心计算。动量方程用一个基于MSIP算法［38］的不完全LU因式分解近似方法（An Incomplete LU Approximate Factorization Technique）进行隐式处理。通过建立压力修正和速度修正的关联，压力修正项将会最终出现在连续性方程中，然后通过重启线性GMRES （m）算法［39］求解压力场。本文所采用方法的优势是压力场和速度场在交互运算的过程中具有良好的兼容性，因为压力修正公式有“精确”的处理形式。由于这个原因，可以省去压力修正算法中弛豫时间中的压力参数，唯一需要用户给定的参数就是局部（或全局）时间步长，而时间步长理论上可以取无限大的值。不包含矩阵的压力求解方法和交错网格离散化相结合，允许压力场的迭代计算在没有压力边界条件的情况下进行。值得注意的是，在过去的Navier-Stokes方程中GMRES算法已被用来加快标准SIMPLE压力修正算法的运算速度［40］。

采用速度一压力的交错耦合本身不足以抑制压力振荡。因此，一个四阶压力耗散项添加到连续性方程中同样被不包含矩阵的求解器进行隐式处理。这一项代表质量流量，不会改变连续性方程的二阶精度。为了防止在高雷诺数的速度场中的奇偶解耦合可以采用乘方格式对动量通量进行离散［41］。

对于湍流问题的处理，本方法采用了三种湍流模型：Baldwin-Lomax［42］的代数模型、Jones和Launder［43］的k-ε模型和Wilcox［44］的k-ω模型。当遇到粗糙网格时，这三种模型都可适应“壁面函数”边界条件。　目前还没有考虑层流到湍流的过渡问题和低雷诺数的影响。然而，层流到湍流的过渡流动问题可以通过在固体壁面上设置过渡位置（强制过渡或自由过渡）来进行处理。

本方法适用于二维定常和非定常翼型绕流问题的研究。在这个例子中，非定常可能是由随时间变化的来流条件或机翼变桨距运动引起的。当机翼变桨的时候，可以使用随翼型运动的相对参考系来研究流动问题，此时，需要将相应的惯性项添加入动量方程中，同时也要对远场的边界条件进行重新定义。但是，为了简化问题，本文将湍流视作准稳态问题来研究，没有对湍流模型进行特别处理。由于这个原因，本方法对湍流的处理仅仅适用于定常的、不存在大量流动分离的情况。

2.3.1　控制方程和湍流模型

变桨翼型的雷诺平均Navier-Stokes方程可通过无量纲的矢量形式表达。相应地，可以将此方程运用在不同位移（比如轴向位移）的情况下。控制方程建立在随翼型运动的相对坐标系上，翼型本身的形变不在考虑范围内。考虑到翼型在绝对坐标系上进行变桨距运动，不可压缩Navier-Stokes的连续性方程和动量方程为

绝对坐标系中的速度矢量V和相对坐标系中的速度矢量W相对应，公式为

式中：Ω是旋转速度矢里（变桨距速度）；r代表位置矢量；下标R代表相对坐标系。

速度被参考速度W∞归一化为尤量纲参数，压力（或者说压力除以密度项ρ）被归一化，VT是被V∞归一化的湍流运动黏度，Re是雷诺数。长度被翼型弦长归一化。

动量方程中包含有翼型变桨距所产生的三个惯性项，代表Coriolis力、角加速度、离心力。当Ω等于0时，W与V相等，公式将退化为标准Navier-Stokes方程组。

本算法中采用了三种湍流模型，分别是Baldwin-Lomax代数模型、k-ε湍流模型和k-ω湍流模型。Baldwin-Lomax代数模型严格遵循Baldwin和Lomax［42］介绍的使用规则，在此不再赘述。

Jones和Launder［43］提出的高雷诺数的标准k-ε湍流模型在此作简要介绍。k-ε方程以如下无量纲形式进行求解

式中：W代表相对速度；Re代表雷诺数。

Pk表达式为

湍流运动黏度vT为

对于标准k-ε模型，相关常数为Prk=1，Pre=1.3，Ce1=1.44，Ce2=1.92，Cm=0.09。在没有低雷诺数的影响下，对于k和ε方程可以通过经典的壁面函数（Wall-Function）确定。

除了k-ε模型，由Wilcox［44］提出来的k-ω模型也运用到程序中。k-ω的方程为

在k-ω模型中，相关常数为α=5/9，β=3/40，β*=9/100，σ=1/2，σ*=1/2。k-ω壁面条件通过壁面函数（wall function）来表示，当网格质量较好的时候也可直接给定壁面数值。当网格质量较好时，k和ω壁面条件为

式中：Δy是壁面到相邻点的距离。

三种湍流模型都被内置于程序中，用户可以在“wall function”边界条件的选项中进行选取。对于壁面函数的运算方式，有以下步骤。

第一步：取平行于壁面的速度分量和壁面相邻位置处的距离△y计算出翼型表面每个网格节点处的摩擦速度UT，摩擦速度UT的求解利用非线性系统的牛顿迭代法

如果则w/UT =κy+，否则w/UT=1/κln （Ey+ ）

式中：K=0.41，是冯卡门常数；的下标Is代表翼型表面的层流边界层的边缘。当使用Baldwin-Lomax模型时，摩擦速度是壁面涡量的函数。

第二步：当采用双方程的k-ε模型或k-ω模型时，壁面相邻处k，ε或ω的值通过下面的公式进行迭代计算。y+定义为

第三步：在动量方程中添加剪切应力。为简化计算，剪切应力是在壁面中段节点位置的湍流运动黏度中人为地添加一个数值实现的。这个人为添加进去的壁面数值随后通过线性外插值法进行迭代运算。

中段节点的数值如下

2.3.2　线性化与时间步长

对于绝对来流速度为V（假设Ω=0），下面将介绍层流非定常Navier-Stokes方程的线性化和时间步长算法。在此基础上可以很容易地将此线性化和时间步长的算法扩展到湍流和变桨距情况。

考虑一个迭代过程，令n代表当前时间点，通过一阶向后欧拉方法对时间导数项进行离散，同时将对流项在n+1附近线性化，由*表示，见图2-1。

由此可得一组隐性时间离散的非定常Navier-Stokes方程为

在此引入的速度和压力修正项δV和δp代表n+1和*时刻的差值。上式Navier-Stokes方程被写成准定常的形式，还需要进行额外的迭代使n+1时刻收敛（δV，δp →0）。将速度修正和压力修正引入动量方程，表达式为

式中：N代表*时刻的非定常（标量）Navier-Stokes算子；R代表*时刻的非定常残差。

将上式引入连续性方程，可以得到“准确的”压力修正公式

由于N-1的存在，压力修正算子▽·｛N-1▽（）｝的解析表达式可以消去。然而在某些简单情况下，例如线性Stokes问题，运用这个算子可以得到一些有用结论。例如，当公式中没有对流项时

前一种情况适用于Re＞＞1的情况，压力修正算子退化成Laplace算子，这也解释了为什么Poisson方程经常被用在压力修正算法中。这个结论只在忽略对流项影响的前提下有效。其他的算法中一般通过N （）矩阵的对角线数值进行近似

上述压力修正算子的任一近似方法都对压力修正项和速度修正项的兼容性有影响，同时也损害迭代过程的收敛性。本论文的想法是，没有必要对压力修正算子进行近似，也无需获得解析表达式，因为式（2-31）可以通过不包含矩阵的共轭梯度方法进行求解。由于压力修正算子是非对称的，在此处选用重启GMRES （m）共轭梯度法，因为这种方法具有较好的收敛性。值得注意的是，这种方法必须进行矩阵N求逆来保持速度场和压力场的交互迭代。为了简化数值运算，矩阵N由因式分解后的近似矩阵M代替。MSIP算法用来实现矩阵N的不完全LU因式分解。

下列过程实现不可压缩Navier-Stokes方程的数值运算过程。

第一步：通过N的值初始化*的值。

第二步：运算矩阵N且对其因式分解，得到矩阵M。

第三步：通过下式计算中间速度场。

第四步：通过重启型GMRES （m）算法对压力修正式（2-31）进行求解，得到δp，更新压力场。

第五步：通过下式更新速度场。

第六步：将*的值替代n+1的值，返回第二步迭代直到非定常动量和连续性残差满足收敛条件。

第七步：如果完成收敛，继续进入下一个时间点，从步骤一开始。

对于定常流动的运算，控制方程在每个时间点并不严格收敛。在每个时间步长内，一般需要GMRES （m）的7到20次内部迭代才能够实现最终收敛，其中m代表选定的Krylov子空间的维度。显然，对于非定常流动问题的运算需要使用内部迭代。根据选定时间步长的不同，内部迭代的次数在5到20次左右。需要注意的是，算法中耗时的部分是速度迭代以及矩阵N的因式分解，因为在每一个内部迭代中都要重复一次。当矩阵M保存之后，压力修正公式的计算都能快速完成。由于压力和速度的更新是兼容的，为了实现结果的收敛一般不需要外部迭代。

在双方程的湍流运算中，湍流模型k-ε或k-ω一般使用类似的方式进行线性化和离散化。

2.3.3　坐标变换

控制方程在计算域的空间中进行离散，计算域是由贴体空间坐标变换而来。（ξ，η）作为二维欧几里得空间（x，y）的贴体参数进行坐标变换x=x （ξ，η），y=y （ξ，η），而（x，y）表示位置矢量r的笛卡尔坐标。共变量为gi，逆变量为gj，i，j=1，2是（ξ，η）的基，其定义如下，下标代表偏导数。

这些基是正交的，满足方程，其中是Kronecker符号。j与关系为

由于正交性条件，有（gij）=1/（gij）。Jacobian行列式为

速度V的散度和对流—扩散算子N（）在（ξ，η）坐标系下表示为

式中：U和V是速度的逆变量。表示为

（u，v）是笛卡尔速度分量。压力梯度项可表述为

2.3.4　离散化

运用有限体积法，控制方程可以在（ξ，η）坐标系的均匀网格上进行离散，此时Δξ=Δη=1。令i和j分别作为ξ和η方向上网格节点的下标，笛卡尔速度分量被存储在网格节点（i，j）上，压力在网格中心位置（i+1/2，j+1/2）进行运算。动量方程的控制体积的中心与网格节点重合，而连续性方程的控制体积的中心与网格细胞中心位置重合，见图2-2。

由于通量（JU）和（JV）在网格细胞中线性分布，速度V的散度如下

动量方程中的对流—扩散算子可以通过Patankar提出的乘方格式进行离散。令

根据乘方格式，可得

其中，Peclet数可以定义为Pe=C/D。对流项和扩散项为

混合的导数系数如ai+1，j+1与1/Re （Jg12）i+1/2j相关。矩阵N （）的对角线值，可以由下式给出

值得注意的是，动量方程的残差也可以运用类似的方法进行偏导离散化。

2.3.5　速度-压力耦合

通过计算离散div.grad （）算子可以发现，速度—压力交错方法本身不足以抑制压力震荡。需要在连续性方程中增加一个四阶线性压力耗散项，式（2-31）由下式代替

式中：L4 （）代表ξ和η方向上四阶导数的和；ε是用户自定义的一阶参数。