一种采用RANS/LES混合技术模拟飞行器摇滚运动的方法

摘要

本发明公开了一种采用RANS/LES混合技术模拟飞行器摇滚运动的方法，目的在于解决现有技术无法准确模拟飞行器大范围分离流动条件下的自激摇滚运动的问题，主要包括：在传统的非定常RANS方法基础上对湍流模型长度尺度进行适当修正，构造出DES类方法；采用预估-校正算法，在双时间步的伪时间迭代过程中实现流场控制方程和刚体运动方程的紧耦合求解，得到飞行器气动参数和运动参数随时间的变化历程。本发明能够在保证一定计算效率的前提下，较准确地模拟复杂流动条件下飞行器自激摇滚运动，为飞行器气动/运动的非线性耦合特性评估提供有效手段。

说明书

技术领域

本发明涉及的是飞行器气动模拟方法，尤其是一种采用RANS/LES 混合技术模拟飞行器摇滚运动的方法。

背景技术

摇滚(Rock)是现代战斗机及导弹武器经常会遇到的气动/运动 耦合现象，它是飞行器在滚转方向的振荡运动。某些情况下伴随滚转 振荡，偏航方向也会发生振荡，称为荷兰滚(Dutch-roll)。随着迎 角增大，荷兰滚逐渐演变成以滚转为主的摇滚运动。摇滚运动通常是 以极限环振荡(Limit Cycle Oscillation，LCO)形式出现。极限环 振荡意味着在一个周期内运动的能量变化为零，即吸收的能量与耗散 的能量相等，振荡既不发散也不收敛，形成等幅等周期振荡。理论上， 在扰动作用下从任意初始状态都会进入摇滚，最后稳定在该极限环状 态，然而一些极限环振荡是不稳定的，此时相平面可能存在其它极限 环，系统在扰动作用下会在极限环间跳变，进入混沌运动。摇滚运动 是非常危险的状态，如果不了解所发生的摇滚运动的特性和流动机理， 控制系统没有合适的控制方法去修正，那么从理论上讲就已经进入失 控飞行。

目前，针对飞行器摇滚运动的预测和模拟，主要通过风洞试验和 数值计算两种方法来开展。

风洞试验能够较真实地模拟飞行器的绕流流态，在模拟摇滚运动 上具有一定的准确性。但风洞试验存在不均匀来流、洞壁干扰、支撑 干扰、轴承摩擦等多方面干扰，而飞行器的自激摇滚运动对来流条件 较为敏感，使风洞模拟结果与真实飞行条件存在一定差别。

数值计算主要是划分飞行器外流场的空间离散网格，并耦合流场 控制方程和刚体运动方程，进行非定常时间推进求解。这种方法不存 在流场干扰，接近真实飞行条件。

目前在流场迭代计算中普遍采用非定常RANS方法。其思想是基 于雷诺平均，将流场信息(速度、压力等)分解成两部分，第一部分 为与时间无关的定常部分，另一部分为脉动部分。因此而产生的雷诺 应力项，通常采用基于Boussinesq涡粘性假设，转化为涡粘系数的 求解问题，从而借助湍流模型进行求解。工程中常用的湍流模型有一 方程S-A模型、双方程、模型等。

RANS方法在附着流动在附着流动或小分离流动的预测中可以得 到较准确得结果，然而当流动存在大面积分离时，RANS方法由于引 入了较大耗散，抑制了小尺度涡运动的解析，因此随着迎角的增大， 模拟精度急剧下降，对准确模拟摇滚运动提出了一定的挑战。

为准确模拟大迎角、复杂流场条件下的飞行器摇滚运动，需要对 每一时刻的流场涡结构进行较为准确的模拟。而RANS/LES混合技术， 结合了RANS方法和LES方法的优点，在高雷诺数、大分离流动的计 算中已经得到了广泛的应用。RANS/LES混合技术的基本思想是对流 场采用分区处理，在近壁面流动中采用RANS方法使得计算量可以承 受，在分离流中采用LES方法使得流场涡结构更加精细。

目前RANS/LES混合技术主要被应用在模型静止或指定运动方式 的非定常分离流动计算中，还没有在气动/运动耦合计算中得到运用， 前者中前一时刻的计算结果只影响后面时刻的流场演化，而在后者中 同时还会对后面时刻的模型运动造成影响。而RANS/LES混合方法将 引入更多的气流脉动信息，与传统的非定常RANS方法相比，收敛速 度变慢，计算误差增大，而计算误差随时间积累对后续流场结构的演 化及飞行器的姿态运动的影响很容易导致模拟出错误的结果。因此， 一方面需要采用更小的时间步长进行时间精确求解，并进行不同时间 步长的比较验证，另一方面需要采用时间精度更高的紧耦合方法开展 气动/运动的耦合计算。

发明内容

本发明的目的，就是针对现有技术所存在的不足，而提供一种采 用RANS/LES混合技术模拟飞行器摇滚运动的方法的技术方案，该方 案采用RANS/LES混合技术模拟飞行器摇滚运动，在保证求解效率 的前提下，更精细地解析复杂分离流动中的涡结构，从而对涡运动主 导的飞行器摇滚运动作出较准确的预测和模拟。

本方案是通过如下技术措施来实现的：一种采用RANS/LES混合 技术模拟飞行器摇滚运动的方法，其特征在于：包括以下步骤：

步骤一、在传统的非定常RANS方法基础上，将湍流模型中长度 尺度替换为DES类方法的长度尺度，使流场解算过程中在壁面附近继 续用RANS方法模拟小尺度的湍流，而在分离区采用类似Smagorinsky 亚格子应力模型来模拟涡结构，从而得到模拟飞行器的流场解算方法；

步骤二、建立滚转单自由度刚体运动方程，并结合步骤一得到的 流场解算方法，采用预估-校正算法，构造气动/运动紧耦合求解方法；

滚转单自由度刚体运动方程如下：

$\left(\begin{array}{c}\frac{{\mathrm{d\omega }}_x}{\mathrm{dt}}=\frac{M_x}{I_{\mathrm{xx}}}\\ \mathrm{d\gamma }/\mathrm{dt}={\omega }_x\end{array}\right)$

式中，ω_x为模型滚转角速度，M_x为模型所受滚转力矩，I_xx为 模型滚转方向的转动惯量，γ为模型滚转角，t为有量纲时间；

在第n个真实时间步计算开始时，使用之前时刻计算得到的气动 力对当前时刻的滚转角速度和滚转角进行预测，预测算法为：

$\left(\begin{array}{c}{\omega }_x^n={\omega }_x^{n-1}+\frac{\mathrm{dt}}{12}·[23{\left(\frac{M_x}{I_{\mathrm{xx}}}\right)}^{n-1}-16{\left(\frac{M_x}{I_{\mathrm{xx}}}\right)}^{n-2}+5{\left(\frac{M_x}{I_{\mathrm{xx}}}\right)}^{n-3}\\ {\gamma }^n={\gamma }^{n-1}+\frac{\mathrm{dt}}{12}·[{23\omega }_x^{n-1}-{16\omega }_x^{n-2}+{5\omega }_x^{n-3}]\end{array}\right)$

式中，上标n表示第n个真实时间步；

随后在冻结当前真实时刻，进行伪时间迭代求解流场的过程中， 使用校正算法对当前运动参数进行校正，校正算法为：

$\left(\begin{array}{c}{{\omega }_x}^n={{\omega }_x}^{n-1}+\frac{\mathrm{dt}}{12}·[5{\left(\frac{M_x}{I_{\mathrm{xx}}}\right)}^{n,\mathrm{new}}+8{\left(\frac{M_x}{I_{\mathrm{xx}}}\right)}^{n-1}-{\left(\frac{M_x}{I_{\mathrm{xx}}}\right)}^{n-2}\\ {\gamma }^n={\gamma }^{n-1}+\frac{\mathrm{dt}}{12}·[{5\omega }_x^{n,\mathrm{new}}+{8\omega }_x^{n-1}-{\omega }_x^{n-2}]\end{array}\right)$

式中，上标“new”表示通过最近一步伪时间迭代得到的值；

步骤三、针对特定飞行器外形划分网格；

步骤四、采用步骤二中的计算方法，针对步骤三得到的空间网格， 开展计算，得到飞行器气动参数和运动参数随时间的变化历程。

作为本方案的优选：步骤三中针对特定飞行器外形划分的网格的 规模稍大于RANS计算所要求网格规模。

本方案的有益效果可根据对上述方案的叙述得知，由于在该方案 中在传统的非定常RANS方法基础上，将湍流模型中长度尺度替换为 DES类方法的长度尺度，使流场解算过程中在壁面附近继续用RANS 方法模拟小尺度的湍流，能够更加准确地模拟出湍流模型的摇滚运动， 能够在保证求解效率的前提下，更精细地解析复杂分离流动中的涡结 构，从而对涡运动主导的飞行器摇滚运动作出较准确的预测和模拟。

由此可见，本发明与现有技术相比，具有实质性特点和进步，其 实施的有益效果也是显而易见的。

附图说明

图1是翼身组合体模型外形示意图。

图2是翼身组合体模型的空间网格划分示意图。

图3是采用RANS方法得到的模型滚转角随时间的变化曲线图。

图4是采用DDES方法得到的模型滚转角随时间的变化曲线图。

图5是以滚转角速度-滚转角相图的形式给出的风洞自由摇滚试 验与采用DDES方法进行自由摇滚数值模拟的结果对比图。

具体实施方式

本说明书中公开的所有特征，或公开的所有方法或过程中的步骤， 除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。

本说明书(包括任何附加权利要求、摘要和附图)中公开的任一 特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加 以替换。即，除非特别叙述，每个特征只是一系列等效或类似特征中 的一个例子而已。

本发明包括以下步骤：

步骤一、在传统的非定常RANS方法基础上，将湍流模型中长度 尺度替换为DES类方法的长度尺度，使流场解算过程中在壁面附近继 续用RANS方法模拟小尺度的湍流，而在分离区采用类似Smagorinsky 亚格子应力模型来模拟涡结构，从而得到模拟飞行器的流场解算方法；

步骤二、建立滚转单自由度刚体运动方程，并结合步骤一得到的 流场解算方法，采用预估-校正算法，构造气动/运动紧耦合求解方法；

滚转单自由度刚体运动方程如下：

$\left(\begin{array}{c}\frac{{\mathrm{d\omega }}_x}{\mathrm{dt}}=\frac{M_x}{I_{\mathrm{xx}}}\\ \mathrm{d\gamma }/\mathrm{dt}={\omega }_x\end{array}\right)$

式中，ω_x为模型滚转角速度，M_x为模型所受滚转力矩，I_xx为 模型滚转方向的转动惯量，γ为模型滚转角，t为有量纲时间；

在第n个真实时间步计算开始时，使用之前时刻计算得到的气动 力对当前时刻的滚转角速度和滚转角进行预测，预测算法为：

$\left(\begin{array}{c}{\omega }_x^n={\omega }_x^{n-1}+\frac{\mathrm{dt}}{12}·[23{\left(\frac{M_x}{I_{\mathrm{xx}}}\right)}^{n-1}-16{\left(\frac{M_x}{I_{\mathrm{xx}}}\right)}^{n-2}+5{\left(\frac{M_x}{I_{\mathrm{xx}}}\right)}^{n-3}\\ {\gamma }^n={\gamma }^{n-1}+\frac{\mathrm{dt}}{12}·[{23\omega }_x^{n-1}-{16\omega }_x^{n-2}+{5\omega }_x^{n-3}]\end{array}\right)$

式中，上标n表示第n个真实时间步；

随后在冻结当前真实时刻，进行伪时间迭代求解流场的过程中， 使用校正算法对当前运动参数进行校正，校正算法为：

$\left(\begin{array}{c}{{\omega }_x}^n={{\omega }_x}^{n-1}+\frac{\mathrm{dt}}{12}·[5{\left(\frac{M_x}{I_{\mathrm{xx}}}\right)}^{n,\mathrm{new}}+8{\left(\frac{M_x}{I_{\mathrm{xx}}}\right)}^{n-1}-{\left(\frac{M_x}{I_{\mathrm{xx}}}\right)}^{n-2}\\ {\gamma }^n={\gamma }^{n-1}+\frac{\mathrm{dt}}{12}·[{5\omega }_x^{n,\mathrm{new}}+{8\omega }_x^{n-1}-{\omega }_x^{n-2}]\end{array}\right)$

式中，上标“new”表示通过最近一步伪时间迭代得到的值；

步骤三、针对特定飞行器外形划分网格；

步骤四、采用步骤二中的计算方法，针对步骤三得到的空间网格， 开展计算，得到飞行器气动参数和运动参数随时间的变化历程。

步骤三中针对特定飞行器外形划分的网格的规模稍大于RANS计 算所要求网格规模。

将本发明应用于翼身组合体在M＝0.6，α＝35°条件下的自由 摇滚运动模拟。

计算采用外形和转动惯量参考某翼身组合体自由摇滚风洞试验 模型，模型全长为504mm，展长为281.4mm，机身直径为63mm，绕模型 x轴(机身中轴)的转动惯量为0.003675kg·m2。

本实施例以Spalart-Allmaras方程湍流模型为基础，修改其中 长度尺度d_w为DDES的长度尺度，构造得到DDES方法。

本实例以翼身组合体模型为基础，划分空间网格，网格单元总量 约592万。如图1、图2所示。

计算条件为M＝0.6，α＝35°，模型参考长度L_ref取机身最大直 径63mm，参考面积S_ref取机身最大横截面积3117.245mm2，基于模型 参考长度的雷诺数Re_L＝1×10⁶。

为了使模型更快进入极限环摇滚运动状态，数值模拟中将模型初 始滚转角设定为5°，并释放模型滚转方向的自由度。

选用不同时间步长，并分别采用RANS方法和DDES方法，结合刚 体运动方程，开展紧耦合计算，得到模型滚转角、滚转角速度、滚转 力矩系数等值随时间的变化历程。

根据不同计算方法的数值模拟结果，分别做出模型滚转角随时间 的变化曲线，如图3、图4所示。采用RANS方法，模型滚转角发生振 荡并很快收敛到平衡位置；采用DDES方法，模型快速进入极限环摇 滚运动。对两种计算方法来说，不同时间步长对模拟结果没有本质区 别。

根据数值模拟结果和风洞试验结果，做出模型滚转角速度-滚转 角相位图，如图5所示。采用DDES方法所得结果与风洞试验对应结果 吻合较好，由于风洞试验存在较多的干扰因素，因此结果中不同周期 的随机性要大一些，而计算结果的随机性相对较小。

风洞试验和不同计算方法对翼身组合体摇滚运动模拟结果如表1 所示：

表1 风洞试验和数值模拟对翼身组合体摇滚运动模拟结果

模拟方法 振幅(°) 频率(Hz) 风洞自由摇滚试验 35.271 7.418 数值模拟(RANS) 运动收敛 7.993 数值模拟(DDES) 35.732 7.621

由表1可知，采用传统的RANS方法对翼体摇滚运动的模拟效果较 差，而采用DDES方法可以得到与风洞试验更接近的结果。

本发明采用了DES类方法对飞行器摇滚运动进行数值模拟，能够 对大面积分离流动的复杂涡运动做出精细解析，从而模拟得到与风洞 试验结果较为一致的摇滚运动。

本发明并不局限于前述的具体实施方式。本发明扩展到任何在本 说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或 过程的步骤或任何新的组合。