一种基于粗网格快速湍流壁面函数气动力预测方法与流程

本发明涉及计算流体力学雷诺平均湍流数值模拟技术领域，具体是一种基于粗网格快速湍流壁面函数气动力预测方法。

背景技术：

湍流数值模拟是高超声速飞行器设计和推进系统一体化设计过程中具有广泛的应用需求，湍流数值模拟效果关系到是否能为高超声速飞行器气动外形优化(包括飞行器头部、翼型、翼展、尾翼、内流道等)和飞行器弹道优化提供更好的研究支撑，在实际应用中起到非常重要的作用，具有很强的现实意义。目前，目前大多采用计算资源消耗较小的rans方法。对于工程中常用的sstk-ω模型，如需精确模拟壁面摩阻和热流，要求壁面网格足够密，通常y+≈1，在工程复杂湍流流动模拟时导致壁面附近网格骤增，计算收敛速度变慢。而工程外形大多非常复杂，局部网格加密不仅影响计算效应，更造成局部网格正交性变差，带来严重的数值刚性问题。

技术实现要素：

为克服现有技术的不足，本发明提供了一种基于粗网格快速湍流壁面函数气动力预测方法，解决现有技术存在的在工程复杂湍流流动模拟时壁面附近网格骤增导致的计算收敛速度变慢、占用计算机内部资源多、具有严重的数值刚性问题等不足。

本发明解决上述问题所采用的技术方案是：

一种基于粗网格快速湍流壁面函数气动力预测方法，包括以下步骤：

s1，建立飞行器的壁面网格模型，获取壁面网格内的飞行器的壁面剪切应力的初始值，其公式为：

其中，τ为剪切应力，τw为壁面剪切应力，μw为壁面分子粘性系数，μt为湍流粘性系数，u为壁面网格内壁面平行速度，y为壁面距离，为偏导数符号，为壁面网格内壁面平行速度在壁面法向的偏导数，up为壁面网格点的速度；

s2，获取壁面网格内的飞行器的摩擦速度uτ，其公式为：

其中，ρw为壁面密度；

s3，计算壁面网格内的无量纲壁面距离y+，其公式为：

比较y+的大小，若y+＞l，则采用壁面函数边界条件，执行uτ的迭代，执行公式为：

然后进入步骤s4；

若y+≤l，则跳出壁面函数，跳转至步骤s8；

其中，ν为运动粘性系数，l为粘性底层边缘对应的无量纲壁面距离，10≤l≤12；

s4，计算壁面网格内的摩擦速度uτ，其公式为：

其中，κ为冯卡门常数，e为积分常数；

s5，计算残差eps，其公式为：

其中，small为无穷小量，

10^-6≤small≤10^-12；

old表示上一循环状态，表示上一循环状态的摩擦速度；

s6，判断摩擦速度uτ是否收敛，若收敛则进入步骤s7，若不收敛则跳转至步骤s3，且设置ite加1，收敛判断公式为：

eps＜morite＞n，

其中，ite表示判断步数，ite初始值为0，m为残差控制量，n为循环控制量，

10^-2≤m≤10^-6，10≤n≤40；

s7，计算壁面网格内的湍流粘性系数，其公式为：

s8，更新壁面网格内的虚拟点湍流粘性系数，其公式为：

其中，为虚拟点湍流粘性系数，为壁面网格点湍流粘性系数。

本发明中，在步骤s3至s6中，应用壁面函数通过子迭代方式更新虚拟点湍流粘性系数，构造一种适用于粗网格的湍流快速气动力预测方法，突破了壁面网格距离限制，不仅在湍流边界层预测方面取得和密网格预测精度非常接近的预测效果，对于存在分离和再附的流动也表现良好，在应用到工程外形上时，在提高两个量阶的壁面距离后，仍可以得到和密网格精度相当的湍流气动力预测结果。而且，由于不用加密网格，避免了由于网格数量增加占用过多计算资源，避免了因占用计算机资源过多、计算量过大而导致计算机程序崩溃、系统过载死机等不足，从而提高了预测效率和可靠性。提供了一种高效率、高精度、高鲁棒性的湍流气动力预测方法。

本发明的目的在于为湍流数值仿真提供一种高效率、高精度、高鲁棒性的湍流气动力模拟方法，可满足航天工程中精细化气动特性要求。

为实现上述目的，本发明采用的总体技术方案是：通过子迭代得到壁面虚拟点湍流粘性系数，从而通过剪切应力的定义修正壁面气动力参数。本方法在程序实施过程中，对源代码改动很小，仅需要设计相应技术模块，更新虚拟点湍流粘性系数，提高了处理效率。

本发明通过在壁面粗网格内基于解析速度分布的壁面函数，通过构造子迭代的方式得到摩擦速度，可有效提高壁面粗网格预测精度，且随着壁面网格距离的增加，网格量整体减小，可显著节约计算时间和计算内存。

本发明通过子迭代方式运用于壁面函数构造一种适用于粗网格的湍流快速气动力预测技术，突破了壁面网格距离限制，不仅在湍流边界层预测方面取得和密网格非常接近的预测效果，对于存在分离和再附的流动也表现良好，在应用到工程外形上时，在提高两个量阶的壁面距离后，仍可以得到和密网格相同的湍流气动力预测结果。同时，随着壁面距离的增加，网格正交性变好，可有效解决密网格带来的数值刚性问题。

作为一种优选的技术方案，步骤s3中，选取粘性底层边缘对应的无量纲壁面距离l为11.06。

作为一种优选的技术方案，步骤s4中，选取冯卡门常数κ为0.4125。

作为一种优选的技术方案，步骤s4中，选取积分常数e为9.793。

作为一种优选的技术方案，步骤s5中，选取无穷小量small为10^-8。

作为一种优选的技术方案，步骤s6中，选取残差控制量m为0.001。

作为一种优选的技术方案，步骤s6中，选取循环控制量n为20。

本发明采用以上数值或参数的近似值，方便更好地计算和提升计算效率，在保证占用资源少、计算精度高、鲁棒性高的前提下，便于提高计算效率和减少工作量。

本发明相比于现有技术，具有以下有益效果：

(1)本发明本发明通过在壁面粗网格内基于解析速度分布的壁面函数，通过构造子迭代的方式得到摩擦速度，可有效提高壁面粗网格预测精度，且随着壁面网格距离的增加，网格量整体减小，可显著节约计算时间和计算内存；

(2)本发明通过子迭代方式运用于壁面函数构造一种适用于粗网格的湍流快速气动力预测技术，突破了壁面网格距离限制，不仅在湍流边界层预测方面取得和密网格非常接近的预测效果，对于存在分离和再附的流动也表现良好，在应用到工程外形上时，在提高两个量阶的壁面距离后，仍可以得到和密网格相同的湍流气动力预测结果；同时，随着壁面距离的增加，网格正交性变好，可有效解决密网格带来的数值刚性问题；

(3)本发明采用以上数值或参数的近似值，方便更好地计算和提升计算效率，在保证占用资源少、计算精度高、鲁棒性高的前提下，便于提高计算效率和减少工作量。

附图说明

图1为本发明所述方法的流程图；

图2、图3为ma＝2.85湍流边界层内平行速度/来流速度随壁面法向距离变化的分布图；

其中，图2为不同网格间距时不用本发明时ma＝2.85湍流边界层内平行速度/来流速度随壁面法向距离变化的分布图，图3为不同网格间距时使用本发明时ma＝2.85湍流边界层内平行速度/来流速度随壁面法向距离变化的分布图；

图4、图5为ma＝2.85压缩拐角不使用和使用本发明壁面摩擦系数对比图；

其中，图4为不使用本发明壁面摩擦系数图，图5为使用本发明壁面摩擦系数图；

图6为ma＝2.85压缩拐角数值平均残差收敛曲线，横坐标为所消耗的计算机时间；

图7为工程复杂高超声速飞行器在来流马赫数为4.0时，攻角和轴向力系数曲线；

图8为本发明所述实施例2的流程图。

附图中字母及相应名称：∞、来流值，u∞、来流速度，cf、摩擦力系数，cfx、轴向力系数，ma，马赫数，α、攻角，β、压缩拐角的倾斜角，wf、壁面函数，△、第一层网格距离，exp、实验值，δo、边界层厚度，resavg、数值平均残差，cputime、计算时间。

具体实施方式

下面结合实施例及附图，对本发明作进一步的详细说明，但本发明的实施方式不限于此。

实施例1

如图1至图8所示，一种基于粗网格快速湍流壁面函数气动力预测方法，包括以下步骤：

s1，建立飞行器的壁面网格模型，获取壁面网格内的飞行器的壁面剪切应力的初始值，其公式为：

其中，τ为剪切应力，τw为壁面剪切应力，μw为壁面分子粘性系数，μt为湍流粘性系数，u为壁面网格内壁面平行速度，y为壁面距离，为偏导数符号，为壁面网格内壁面平行速度在壁面法向的偏导数，up为壁面网格点的速度；

s2，获取壁面网格内的飞行器的摩擦速度uτ，其公式为：

其中，ρw为壁面密度；

s3，计算壁面网格内的无量纲壁面距离y+，其公式为：

比较y+的大小，若y+＞l，则采用壁面函数边界条件，执行uτ的迭代，执行公式为：

然后进入步骤s4；

若y+≤l，则跳出壁面函数，跳转至步骤s8；

其中，v为运动粘性系数，l为粘性底层边缘对应的无量纲壁面距离，10≤l≤12；

s4，计算壁面网格内的摩擦速度uτ，其公式为：

其中，κ为冯卡门常数，e为积分常数；

s5，计算残差eps，其公式为：

其中，small为无穷小量，

10^-6≤small≤10^-12；

old表示上一循环状态，表示上一循环状态的摩擦速度；

s6，判断摩擦速度uτ是否收敛，若收敛则进入步骤s7，若不收敛则跳转至步骤s3，且设置ite加1，收敛判断公式为：

eps＜morite＞n，

其中，ite表示判断步数，ite初始值为0，m为残差控制量，n为循环控制量，

10^-2≤m≤10^-6，10≤n≤40；

s7，计算壁面网格内的湍流粘性系数，其公式为：

s8，更新壁面网格内的虚拟点湍流粘性系数，其公式为：

其中，为虚拟点湍流粘性系数，为壁面网格点湍流粘性系数。

本发明中，在步骤s3至s6中，应用壁面函数通过子迭代方式更新虚拟点湍流粘性系数，构造一种适用于粗网格的湍流快速气动力预测方法，突破了壁面网格距离限制，不仅在湍流边界层预测方面取得和密网格预测精度非常接近的预测效果，对于存在分离和再附的流动也表现良好，在应用到工程外形上时，在提高两个量阶的壁面距离后，仍可以得到和密网格精度相当的湍流气动力预测结果。而且，由于不用加密网格，避免了由于网格数量增加占用过多计算资源，避免了因占用计算机资源过多、计算量过大而导致计算机程序崩溃、系统过载死机等不足，从而提高了预测效率和可靠性。提供了一种高效率、高精度、高鲁棒性的湍流气动力预测方法。

本发明的目的在于为湍流数值仿真提供一种高效率、高精度、高鲁棒性的湍流气动力模拟方法，可满足航天工程中精细化气动特性要求。

为实现上述目的，本发明采用的总体技术方案是：通过子迭代得到壁面虚拟点湍流粘性系数，从而通过剪切应力的定义修正壁面气动力参数。本方法在程序实施过程中，对源代码改动很小，仅需要设计相应技术模块，更新虚拟点湍流粘性系数，提高了处理效率。

本发明通过在壁面粗网格内基于解析速度分布的壁面函数，通过构造子迭代的方式得到摩擦速度，可有效提高壁面粗网格预测精度，且随着壁面网格距离的增加，网格量整体减小(壁面距离y增大，则y+增大，壁面网格变稀，网格量整体减小)，可显著节约计算时间和计算内存。

本发明通过子迭代方式运用于壁面函数构造一种适用于粗网格的湍流快速气动力预测技术，突破了壁面网格距离限制，不仅在湍流边界层预测方面取得和密网格非常接近的预测效果，对于存在分离和再附的流动也表现良好，在应用到工程外形上时，在提高两个量阶的壁面距离后，仍可以得到和密网格相同的湍流气动力预测结果。同时，随着壁面距离的增加，网格正交性变好，可有效解决密网格带来的数值刚性问题。

值得说明的是，在使用时，首先步骤s2根据摩擦速度定义预估一个摩擦速度uτ，为其初始值，在步骤s3中，利用壁面函数符合相应条件则进行迭代替换，更新原有的摩擦速度uτ。

值得说明的是，本发明采用了通过子迭代方式运用于壁面函数构造的技术手段，解决了现有技术存在的在工程复杂湍流流动模拟时计算收敛速度变慢、占用计算机内部资源多、具有严重的数值刚性问题等技术问题，达到了高效率、高精度、高鲁棒性、数值刚性问题得以解决的技术效果。

值得说明的是，一方面，本发明所述方法对计算机外部的处理对象为外部技术数据(如剪切应力，壁面数据，湍流，分子粘性数据，速度数据等具有特定物理含义的数据等，并非单纯的数学数值数据)，本发明通过对外部技术数据的处理，可满足航天工程中精细化气动特性要求；另一方面，本发明所述方法对计算机内部性能做了改善，通过在壁面粗网格内基于解析速度分布的壁面函数，通过构造子迭代的方式得到摩擦速度与计算机内部构建进行交互，减少了计算机内部资源的消耗，提高了计算速度。因此，本发明采用遵循自然规律的技术手段、解决了实际的技术问题、获得了符合自然规律的技术效果，本发明要保护的方案属于技术方案，本发明要保护的技术方案符合专利法第二条规定，符合专利法保护的客体。

作为一种优选的技术方案，步骤s4中，选取冯卡门常数κ为0.4125。

作为一种优选的技术方案，步骤s4中，选取积分常数e为9.793。

作为一种优选的技术方案，步骤s5中，选取无穷小量small为10^-8。

作为一种优选的技术方案，步骤s6中，选取残差控制量m为0.001。

作为一种优选的技术方案，步骤s6中，选取循环控制量n为20。

本发明采用以上数值或参数的近似值，方便更好地计算和提升计算效率，在保证占用资源少、计算精度高、鲁棒性高的前提下，便于提高计算效率和减少工作量。

实施例2

如图1至图8所示，本实施例基于实施例1，反映采用本发明的技术效果。

该技术通过在壁面粗网格内基于解析速度分布的壁面函数，通过构造子迭代的方式得到摩擦速度，可有效提高壁面粗网格预测精度，且随着壁面网格距离的增加，网格量整体减小，可显著节约计算时间和计算内存。

图2、图3是来流ma＝2.85的高超平板湍流边界层速度分布，当不使用壁面函数时(图2)，粗网格数值结果和密网格结果、试验结果偏差较大，且随着网格壁面距离的增加，精度越差，而图3中使用壁面函数后，粗网格结果得到极大改善，几乎和密网格、试验结果重合。

具体地：

a)图2给出了不同网格间距时不用本发明时ma＝2.85湍流边界层内平行速度/来流速度随壁面法向距离变化的分布图，从图中可知湍流模型不使用该发明时，随着第一层壁面网格距离的增大，网格逐渐变稀，边界层内速度分布和实验值偏差逐步扩大，对于壁面网格间距△＝1.0mm网格，和实验值差异最大可以达到30％。而速度分布变差意味着速度梯度的差异，即影响摩擦系数的预测，对于高超声速飞行器来说，飞行器大部分的表面都处于湍流边界层内，当摩擦系数预测变差，会影响飞行器阻力的预测精度，而阻力是飞行器设计时的一项重要指标。

b)图3给出了不同网格间距时使用本发明时ma＝2.85湍流边界层内平行速度/来流速度随壁面法向距离变化的分布图，从图中可知使用本发明后，即使网格从最密的△＝0.001mm逐渐变稀到△＝1mm边界层内速度分布依然和实验值高度吻合，表明使用本发明后可以极大的提高粗网格壁面摩擦系数，即摩阻的预测能力。

图4、图5给出了ma＝2.85压缩拐角壁面摩擦系数分布，流动中存在复杂的分离和再附现象(摩擦系数cf由正转负的点为分离点，右负转正的点为再附点，分离点在再附点之间的区域为分离区)，其中图4是不使用壁面函数时结果，壁面间距增大，预测结果大幅下降，即使对于分离区前(横坐标小于-3部分)，预测结果也非常糟糕。图5给出使用壁面函数后数值结果，粗网格预测精度相较与图4明显改善，且在分离区前(横坐标小于-3部分)，摩擦系数预测结果接近密网格和试验结果。

具体地：

c)图4给出了不同网格间距时不使用本发明时ma＝2.85压缩拐角壁面摩擦系数分布，从图中可知△＝0.001mm密网格同实验值吻合最好，随着网格间距增大，壁面摩擦系数同实验值差异越大，其中粗网格△＝1.0mm预测的摩擦系数仅为实验值20％，严重低估摩擦系数。压缩拐角的激波边界层干扰是高超声速飞行器中常见的现象，例如机翼、尾翼和机身的链接处，内流道中的突起物等。

图5给出了不同网格间距时使用本发明时ma＝2.85压缩拐角壁面摩擦系数分布，其中粗网格的预测能力相较与图4得到大幅提升，虽然较密网格在分离区附近任有一定差距，但考虑到使用粗网格后壁面网格可以大幅减少，且随着壁面距离的增加，在相同计算设置下，可大幅提高计算步长，因此计算收敛速度可以得到大幅提升，而工程中使用的阻力是这个飞行器表面的积分量，因此小范围的预测偏差并不影响整体的预测精度。

图6给出了ma＝2.85压缩拐角数值平均残差收敛曲线，横坐标为所消耗的计算机时间，从图中可知当残差下降到10^-4时，粗网格(△＝1.0mm)仅需密网格(△＝0.001mm)的25％，计算时间减少了约75％。图7给出了工程复杂高超声速飞行器在来流马赫数为4.0时，攻角和轴向力系数曲线，从图中可知，粗网格(△＝0.1mm)在使用该技术后计算精度显著提升，接近密网格(△＝0.0002mm)结果。

具体地：

图6比较ma＝2.85压缩拐角算例消耗的计算时间和整体残差之间的对比，从图中可知粗网格△＝1mm在下降到残差10^-4时，消耗的计算时间仅为密网格△＝0.001mm计算时间的四分之一，表面使用本发明后在提升壁面摩擦系数(如图5所示)的同时，可以进一步节约计算时间，同时由于网格量的减少还可以节约计算内存。

图7给出了典型高超声速飞行器攻角和轴向力系数曲线，飞行马赫数ma＝4.0其中横坐标表示飞行器飞行攻角，纵坐标为轴向力，其中飞行器包含机翼、v尾、进气道等。从图中可知对于粗网格△＝0.1mm不使用本发明时，轴向力系数同密网格△＝0.002mm的差距较大，大约为10％，而使用该发明后(图中wf＝1)轴向力系数预测精度显著提升，周向力差距在1％左右。轴向力系数是计算阻力的重要组成，是高超声速飞行器外形设计、弹道设计中非常重要的参数之一。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，依据本发明的技术实质，在本发明的精神和原则之内，对以上实施例所作的任何简单的修改、等同替换与改进等，均仍属于本发明技术方案的保护范围之内。