CN107742050A - 一种面向混合网格高精度阻力预测的熵修正方法 - Google Patents

Abstract

本发明公开了一种面向混合网格高精度阻力预测的熵修正方法，主控方程对流项采用二阶迎风Roe通量差分分裂格式进行离散，粘性项采用中心差分格式离散；在传统Harten‑Yee熵修正基础上，以界面上的法向速度与界面上速度的比值为判断依据，对可能是驻点或者激波所在区域保留熵修正的效果，对平行于翼面流向的界面不使用熵修正或者熵修正尽量小；湍流模型采用S‑A一方程湍流模型，湍流控制方程空间离散采用一阶迎风格式。本发明在不改变其它区域熵修正的同时，在平行于翼面流向的界面上不做熵修正或者减小熵修正，尽量不改变此处的真实耗散，从而把熵修正对阻力预测精度的影响降到了最低，同时保留使用熵修正带来的程序鲁棒性等优点。

Description

一种面向混合网格高精度阻力预测的熵修正方法

技术领域

本发明涉及空气动力学领域，具体涉及一种面向混合网格高精度阻力预测的熵修正方法。

背景技术

Roe格式通量差分分裂格式是以近似Riemann分解为基础的Godunov类求解器，具有耗散小、接触间断分辨率高和激波捕获性能强、可以较好模拟边界层流动等优点，在亚跨超声速流场模拟中得到广泛应用。但对于特定问题Roe格式有时不满足熵条件，会产生非真实的膨胀波，导致非物理解，在强激波和音速点附近，需要对原始的Roe格式进行熵修正。

通常的Harten-Yee熵修正方法，由于翼面法向的速度是小量，导致平行于流向的界面处出现很小的特征值，Harten-Yee熵修正可能会显著提高这一特征值，以避免出现非物理解。但是这样就加大了该界面上物理量的耗散，特别是在附面层内对物面法向速度的求解精度影响较大，导致不能准确地再现附面层的速度剖面和壁面附近的切应力，使摩擦阻力的预测出现偏差。已有的有关Roe格式熵修正方法研究的文章，要么是对已有熵修正的性能做评估，要么是通过改进声速膨胀波附近的熵修正以消除非物理膨胀激波，而有关通过改进熵修正方法提高附面层内阻力预测精度，尤其是提高三维非结构网格阻力预测精度的文章还未见到。

本发明针对非结构混合网格中Roe格式熵修正的特点，通过改进传统的Roe格式Harten-Yee熵修正方法，提出了一种可提高非结构混合网格粘性计算精度的Harten-Yee熵修正改进方法，改进后的熵修正方法残差收敛特性与原始Harten-Yee熵修正一致，计算结果和无熵修正时的计算结果基本一致，改进后的熵修正方法既保留了使用熵修正带来的程序鲁棒性等优点，同时把熵修正对阻力预测精度的影响降到最低，提高了阻力预测精度。

发明内容

为了克服现有技术的上述缺点，本发明提出了一种面向混合网格高精度阻力预测的熵修正方法。

本发明解决其技术问题所采用的技术方案是：一种面向混合网格高精度阻力预测的熵修正方法，包括如下步骤：

步骤一、根据计算模型的数字外形，进行空间离散，生成非结构混合网格；

步骤二、主控方程对流项采用二阶迎风Roe通量差分分裂格式进行离散，粘性项采用中心差分格式离散；

步骤三、对Roe格式进行熵修正：在传统Harten-Yee熵修正基础上，以界面上的法向速度与界面上速度的比值为判断依据，对可能是驻点或者激波所在区域保留熵修正的效果，对平行于翼面流向的界面不使用熵修正或者熵修正尽量小；

步骤四、湍流模型采用S-A一方程湍流模型，湍流控制方程空间离散采用一阶迎风格式；

步骤五、通过迭代运算，得到计算模型的气动特性数据和流场。

与现有技术相比，本发明的积极效果是：本发明提出一种面向非结构混合网格高精度阻力预测的Roe格式熵修正方法，该方法通过改进传统的Roe格式Harten-Yee熵修正方法，在不改变其它区域熵修正的同时，在平行于翼面流向的界面上不做熵修正或者减小熵修正，尽量不改变此处的真实耗散，从而把熵修正对阻力预测精度的影响降到了最低，同时保留使用熵修正带来的程序鲁棒性等优点。

附图说明

本发明将通过例子并参照附图的方式说明，其中：

图1是DLR-F4翼身组合体空间离散网格；

图2是网格不同位置示意图；

图3是无熵修正、Harten-Yee熵修正和本文熵修正三种方法对主控方程变量ρ的残差收敛曲线影响；

图4是无熵修正、Harten-Yee熵修正和本文熵修正三种方法计算的气动力和力矩的值，以及与CFD++、CFL3D、NSU3D等知名CFD软件和NLR、ONERA、DRA等风洞得到的结果比较；

图5是机翼剖面弦向中段速度型比较及放大图；

图6是机翼前缘附近速度型比较及放大图；

其中：no表示无熵修正，modified表示本发明改进的传统Harten-Yee熵修正，original表示传统的Harten-Yee熵修正。

具体实施方式

一种面向混合网格高精度阻力预测的熵修正方法，控制方程采用雷诺平均Navier-Stokes方程，空间离散采用有限体积法，未知变量位于网格单元的体心，计算网格为非结构混合网格，在物面附近附面层内采用的是扁长型的三棱柱网格，主控方程对流项采用二阶迎风Roe通量差分分裂格式进行离散，粘性项采用中心差分格式离散，并采用多重网格技术进行收敛加速。Roe通量差分分裂格式的熵修正是在传统Harten-Yee熵修正基础上，根据网格不同位置，以界面上的法向速度与界面上速度的比值为判断依据，实现不同区域熵修正的不同，降低熵修正对附面层速度剖面和壁面附近切应力的影响，提高阻力的预测精度。湍流模型采用S-A一方程湍流模型，湍流控制方程空间离散采用一阶迎风格式。主控方程和湍流方程时间迭代均采用LU-SGS方法。

本发明方法具体包括以下步骤：

步骤一、控制方程采用非定常可压缩雷诺平均Navier-Stokes方程：

守恒形式的非定常可压缩N-S方程可写成如下积分形式：

其中，Ω为控制体体积，为控制体封闭面面积，表示守恒变量，表示无粘通量，表示粘性通量；

步骤二、根据计算模型的数字外形，进行空间离散，生成非结构混合网格：

根据DLR-F4翼身组合体的数模，进行空间离散，生成非结构混合网格，计算时采用半模进行计算，远场边界取约50倍的机翼平均气动弦长，半模网格单元总数为2164万，其中三棱柱795万,四面体1368万。物面法向三棱柱层数为27层，第一层间距约为1.0×10^-6m(y⁺≈1)，物面单元数为29.5万。机翼后缘采用各向异性三角形网格，单元数为32个，生成的网格见图1。

步骤三、控制方程对流项采用二阶迎风Roe通量差分分裂格式进行离散，粘性项采用中心差分格式离散：

Roe格式在控制体单元面上的通量表达式为：

Roe平均矩阵和左右状态差的乘积的求法如下：

其中：

当特征值很小时，Roe格式会违反熵条件，产生非物理解，如膨胀激波、Carbuncle现象等。为避免出现非物理解，需要对Roe平均矩阵的特征值进行熵修正。

步骤四、对Roe格式进行熵修正，在传统的Harten-Yee熵修正基础上，以界面上的法向速度与界面上速度的比值为判断依据，对可能是驻点，或者激波所在区域保留熵修正的效果，对平行于翼面流向的界面不使用熵修正或者熵修正尽量小：

传统的Harten-Yee熵修正方法，公式如下：

其中δ定义为：

δ^*是一个小值，一般取0.0～0.4。

传统的Harten-Yee熵修正方法，由于翼面法向的速度是小量，如图2的A处，导致平行于流向的界面处出现很小的特征值，Harten-Yee熵修正可能会显著提高这一特征值，以避免出现非物理解。但是这样就加大了该界面上物理量的耗散，导致不能准确地再现附面层的速度剖面和壁面附近的切应力，使摩擦阻力的预测出现偏差。

对于每个网格面我们定义一个参数c，将每个面的c值乘以传统Harten-Yee熵修正公式中的参数δ^*，得到新的δ^*值：

δ^*＝cδ^*

参数c的值根据网格位置不同取值为0～1，当c的值为0时，表示此位置不使用熵修正，当c的值为1时，表示此位置采用传统的Harten-Yee熵修正，不影响熵修正的效果。

我们希望在物面附近类似图2中A位置的网格区域不使用熵修正或者熵修正尽量小；但对物面附近类似图2中B位置的网格区域，虽然同样有扁长型网格，因为这个位置可能会是驻点，或者激波所在区域，我们不希望在这些网格的任何界面上减小Harten-Yee熵修正，因此需要将这些区域剔除，在此使用界面上的法向速度与界面上速度的比值来剔除这些界面。为此，我们定义的参数c表达式如下：

其中d表示所求面相邻的两个单元的体心之间的距离，S是该面的面积。对于平行于物面流向的界面来说，d²<<S，因此c≈0；而对于垂直于翼面的界面，d²>>S，因此c≈1；而对于附面层外的网格，一般来说c≈1，不影响熵修正的效果；对于激波区域或驻点附近的网格，公式中的第一项虽然很小，但是由于V_n≈V，第二项的值接近1，因此c≈1。

本步骤中对Roe格式采用的熵修正方法，是对传统Harten-Yee熵修正的改进，可以实现基本不改变其它区域熵修正的同时，在平行于翼面流向的界面上不做熵修正或减小熵修正，尽量不改变此处的真实耗散，既保证了计算程序的鲁棒性，又把熵修正对阻力预测精度的影响降到了最低，提高了阻力的预测精度。

步骤五、湍流模型采用S-A一方程湍流模型，湍流控制方程空间离散采用一阶迎风格式，时间迭代均采用LU-SGS方法：

在S-A湍流模型中，湍流粘性系数μ_t通过一个中间变量来计算：

由以下控制方程计算得到：

其中：G_ν是生成的湍流粘度，Y_ν是壁面附近区域由于壁面阻尼效应消耗的湍流粘度，C_b2是常量，ν是分子动粘度，是自定义的源项，μ_t是湍流粘度。

步骤六、通过迭代运算，得到计算模型的气动特性数据和流场。

图3给出了无熵修正、Harten-Yee熵修正和改进后的Harten-Yee熵修正三种方法对主控方程变量ρ的残差收敛曲线影响。可以看出，无熵修正方法的残差只下降了不到4个量级，而使用了熵修正方法后下降了7个量级，在6000步后还将继续下降，并且Harten-Yee熵修正和改进后的Harten-Yee熵修正残差收敛曲线基本一致，收敛性都很好。图4的表中列出了无熵修正、Harten-Yee熵修正和改进后的Harten-Yee熵修正三种方法计算的气动力和力矩的值，以及CFD++、CFL3D、NSU3D等知名CFD软件和NLR、ONERA、DRA等风洞得到的结果，可以看出，本文三种不同熵修正方法计算得到的结果与知名CFD软件计算得到的结果基本一致，与风洞试验结果也吻合较好。总的来说，两种熵修正方法对气动力系数和力矩系数影响较小，相较而言，改进后的熵修正对气动力和力矩影响更小，主控方程残差收敛特性更好，下降量级更多，计算得到的气动力系数和力矩系数与试验结果更接近。改进后的Harten-Yee熵修正方法计算结果与原始Harten-Yee熵修正方法计算结果相比，升力减小约0.7％，阻力减小3.6个阻力单位，低头俯仰力矩减小约0.6％。改进后的Harten-Yee熵修正方法计算结果与无熵修正计算结果相比，在取四位有效数字情况下气动力和力矩值完全相同，升力只有在第5位有效数字，阻力在第6位有效数字才不相同，说明本文的修改将Harten熵修正对气动力的影响降到最低限度，基本上可以忽略不计。可以看出，改进后的熵修正，软件鲁棒性和原始的Harten-Yee熵修正方法基本一样，但计算精度提高，计算结果和无熵修正时的计算结果基本一致，达到了增加软件鲁棒性，同时提高软件阻力预测准度的目的。

我们改进熵修正是希望在类似图2中A处的翼型中部区域附面层内尽量减小或者不加熵修正，而在类似图2中B处的翼型前缘区域附面层内熵修正尽量与原始的Harten-Yee熵修正保持一致。图5是机翼剖面弦向中段速度型的比较，其中右图是左图箭头部分的放大。熵修正对速度型影响比较小，但是从放大图可见，原始的熵修正不仅影响速度矢量的大小，而且影响方向，而本文改进的熵修正，速度型与没用熵修正的结果基本上重合在一起。图6是前缘附近的速度型比较，由图可见，改进后的熵修正与原始的熵修正基本重合，而与无熵修正的结果不同。因此，本文对熵修正的改进达到了我们所希望的效果。

Claims (6)

1.一种面向混合网格高精度阻力预测的熵修正方法，其特征在于：包括如下步骤：

步骤一、根据计算模型的数字外形，进行空间离散，生成非结构混合网格；

步骤二、主控方程对流项采用二阶迎风Roe通量差分分裂格式进行离散，粘性项采用中心差分格式离散；

步骤三、对Roe格式进行熵修正：在传统Harten-Yee熵修正基础上，以界面上的法向速度与界面上速度的比值为判断依据，对可能是驻点或者激波所在区域保留熵修正的效果，对平行于翼面流向的界面不使用熵修正或者熵修正尽量小；

步骤四、湍流模型采用S-A一方程湍流模型，湍流控制方程空间离散采用一阶迎风格式；

步骤五、通过迭代运算，得到计算模型的气动特性数据和流场。

2.根据权利要求1所述的一种面向混合网格高精度阻力预测的熵修正方法，其特征在于：步骤一所述的空间离散采用有限体积法，未知变量位于网格单元的体心，在物面附近附面层内采用的是扁长型的三棱柱网格。

3.根据权利要求1所述的一种面向混合网格高精度阻力预测的熵修正方法，其特征在于：步骤二所述主控方程采用非定常可压缩雷诺平均Navier-Stokes方程。

4.根据权利要求1所述的一种面向混合网格高精度阻力预测的熵修正方法，其特征在于：所述主控方程和湍流控制方程时间迭代均采用LU-SGS方法。

5.根据权利要求1所述的一种面向混合网格高精度阻力预测的熵修正方法，其特征在于：步骤三所述对Roe格式进行熵修正时，对于每个网格面定义一个参数c，将每个网格面的参数c值乘以传统Harten-Yee熵修正公式中的参数δ^*，得到新的δ^*值：

δ^*＝cδ^*

其中，参数c的值根据网格位置不同取值为0～1。

6.根据权利要求5所述的一种面向混合网格高精度阻力预测的熵修正方法，其特征在于：参数c采用如下公式确定：

<mrow> <mi>c</mi> <mo>=</mo> <mi>M</mi> <mi>A</mi> <mi>X</mi> <mrow> <mo>(</mo> <mfrac> <mfrac> <mi>d</mi> <msqrt> <mi>S</mi> </msqrt> </mfrac> <msqrt> <mrow> <mn>1</mn> <mo>+</mo> <msup> <mfrac> <mi>d</mi> <mi>S</mi> </mfrac> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> <mo>,</mo> <mfrac> <msub> <mi>V</mi> <mi>n</mi> </msub> <mi>V</mi> </mfrac> <mo>)</mo> </mrow> </mrow>

其中，V_n表示界面上的法向速度，V表示界面上速度，d表示所求面相邻的两个单元的体心之间的距离，S表示该面的面积。