ABSTRACT

A numerical simulation investigation based on the Navier Stokes equations of dynamic characteristics of missile geometries with fins in Supersonic Environment was presented. The 2-order Roe scheme and the dual-time-step method are respectively applied to discrete of the spatial and time derivative of the unsteady fluency. Based on the unsteady Etkin pneumatics model, dynamic stability derivatives were calculated from load history of the unsteady flow. The influence of frequency of force-oscillation, inner-iteration steps and mesh dimension to aerodynamic force/moment was investigated. The dynamic stability derivatives predicted by the present methods quite agree with the experiment results.

Keywords:Numerical Simulation, N-S Equations, Supersonic, Dynamic Derivative, Force-Harmonic Analysis Method

超音速条件下复杂带翼外形的强迫振动分析

刘绪，王光辉，范化喜，张金春

海军潜艇学院，山东 青岛

收稿日期：2019年2月9日；录用日期：2019年2月22日；发布日期：2019年3月1日

摘 要

基于三维非定常N-S (Navier Stokes)方程开展超音速条件下“十”字翼导弹动态特性数值研究。空间方向的离散采用二阶精度的Roe格式，非定常流动求解采用“双时间步”(Dual-time-step)方法的隐式LU-SGS格式。在Etkin非定常气动力模型基础上，通过强迫简谐运动过程辨识动态稳定性参数(动导数)，分析了减缩频率、子迭代步数和网格规模对动导数辨识结果的影响，并通过与实验结果的比较，验证了计算结果的可靠性。

关键词 :数值模拟，N-S方程，超音速，动导数，强迫简谐分析法

Copyright © 2019 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

动态稳定性参数(动导数)是飞行器控制系统设计，飞行器轨道设计和飞行器纵、横向动态稳定性分析的重要参数 [1] 。目前，评估飞行器飞行品质和动态特性的方法主要有飞行实验、风洞实验和理论计算，其重要内容之一是获取动导数。飞行实验和风洞实验确定动态参数难度大、周期长，且实验状态受实验条件的诸多限制 [2] 。从国内外发展来看，通过数值计算获得动导数是动态参数研究发展的重要趋势 [3] 。

动导数计算方法分为：1) 小扰动线化理论、牛顿理论与气动力工程模型相结合的动导数工程近似方法；2) 模拟动态实验的数值自由振荡法和数值强迫振荡法。动导数工程近似方法在90年代以前被国内外普遍采用，这是一种经验和半经验的方法，如内伏牛顿流理论等。工程近似方法的最大特点在于快捷，但依赖于经验性，考虑气流分离、再附和尾流等效应较为困难，对于超音速复杂外形流动难以适应，计算精度难以保障。因此考虑了流动非线性的模拟动态实验的全数值计算是动导数研究发展的重要方向 [4] 。

本文采用数值模拟强迫振荡法辨识动态稳定性参数，在Etkin非定常气动力模型的基础上，建立了单自由度俯仰及滚转运动时的非定常气动力/力矩系数和状态变量之间的数学模型，给出了强迫简谐分析方法求解俯仰阻尼导数及滚转阻尼导数的数值计算方法。此外，为了考虑俯仰引起的纵横向交叉耦合效应，本文还开展了“十”字翼导弹强迫俯仰振动时的交叉耦合导数计算。

2. 数值方法

在笛卡尔惯性坐标系下，对完全气体、忽略质量力下的三维无量纲Navier-Stokes方程形式如下：

$\frac{\partial U}{\partial t}+\frac{\partial E}{\partial x}+\frac{\partial F}{\partial y}+\frac{\partial G}{\partial z}=\frac{\partial E_v}{\partial x}+\frac{\partial F_v}{\partial y}+\frac{\partial G_v}{\partial z}$ (1)

式中：U为守恒变量，E、F、G为无粘通量，E_v、F_v、G_v为粘性通量。本文在空间上采用二阶精度的Roe格式，Roe格式是基于Godunov方法的基本思路发展起来的，是通量差分分裂(FDS)格式的一种，由于其优秀的间断分辨率和较小的数值耗散性，目前受到广泛的应用。

采用引入“双时间步”(Dual-time-step)方法的LU-SGS隐式格式离散流体运动方程时间项。动网格生成采用刚性网格生成方法。远场入流边界采用适用于动态边界条件下的一维Riemann不变量的无反射边界条件。对于定常/非定常超音速出流，边界点上的值由内流场计算结果外推得到。对于壁面边界，采用速度无滑移、绝热壁条件。奇性轴采用外插后周向平均处理。

3. 动导数计算方法

对于飞行器做单自由度的俯仰运动，如果其质心速度不变，则确定运动的独立状态变量只有攻角 $\alpha$ 和俯仰角速度q。对俯仰力矩系数C_m，由Etkin假定，俯仰力矩系数可写成 [5] ：

$C_m=C_m\left(\alpha \left(t\right),\stackrel{\dot{}}{\alpha }\left(t\right),\ddot{\alpha }\left(t\right),\cdots ,q\left(t\right),\stackrel{\dot{}}{q}\left(t\right),\ddot{q}\left(t\right),\cdots \right)$ (2)

给定强迫振动：

$\alpha \left(t\right)=\theta ={\alpha }_0+{\alpha }_m\sin kt$ (3)

式中 $\alpha$ 为瞬时攻角， $\theta$ 为俯仰角， ${\alpha }_0$ 为基准状态的攻角， ${\alpha }_0$ 为振动幅度，k为减缩频率。

在k不很大，且忽略高阶导数的影响时，当非定常振动过程达到谐振解，通过数值积分可以求出俯仰刚度 $C_{m_{\theta }}$ 和俯仰阻尼导数 $C_{m_{\stackrel{\dot{}}{\theta }}}$ ：

$C_{m_{\theta }}=\frac{k}{{\theta }_m\text{π}}{\displaystyle {\int }_{t_s}^{t_s+T_c}{C}_m\left(t\right)\sin kt\text{d}t}$ (4)

$C_{m_{\stackrel{\dot{}}{\theta }}}=\frac{1}{{\theta }_m\text{π}}{\displaystyle {\int }_{t_s}^{t_s+T_c}{C}_m\left(t\right)\cos kt\text{d}t}$ (5)

$t_s$ 为积分起始时刻， $T_c$ 为振动周期。同样的，强迫俯仰振动时引起的俯仰-滚转力矩交叉耦合导数为：

$C_{m_{\stackrel{\dot{}}{\alpha }}}+C_{m_q}=\frac{1}{{\theta }_m\text{π}}{\displaystyle {\int }_{t_s}^{t_s+T_c}{C}_l\left(t\right)\cos kt\text{d}t}$ (6)

与之相仿，给定强迫滚转振动：

$\psi ={\psi }_m\sin kt$ (7)

其中 $\psi$ 为滚转角， ${\psi }_m$ 为滚转角振幅。由此辨识出的滚转阻尼导数为：

$C_{l\stackrel{\dot{}}{\psi }}=C_{lp}+C_{l\stackrel{\dot{}}{\beta }}\sin \left({\alpha }_0\right)=\frac{1}{\text{π}{\psi }_m}{\displaystyle {\int }_{t_s}^{t_s+T_c}{C}_l\left(t\right)\cos kt\text{d}t}$ (8)

4. 验证算例

本文采用HBS(Hyper Ballistic Shape)外形作为验证算例 [6] ，HBS外形是高超声速导弹外形的标模。HBS外形参数见图1，其中： $$ 。图2为HBS模型对称面网格，网格数： $131\times 73\times 51$ (流向 × 周向 × 径向)。

图1. HBS外形参数

图2. HBS对称面网格

图3. 俯仰力矩系数迟滞曲线

计算条件为：马赫数 $M{a}_{\infty }=6.85$ ，以头部直径为参考长度的 ${\mathrm{Re}}_d=0.72\times {10}^6$ ，质心位置 $X_{cg}/L=0.72$ ，强迫俯仰振荡的振幅为1˚，减缩频率分别取0.2和0.05。图4是不同起始攻角下对俯仰阻尼导数的计算结果与文献 [6] 中的实验结果的比较。实验结果的误差带为±15%。从计算的总体情况来看：数值结果和实验结果在趋势上是一致的，都处于俯仰动态稳定状态。从数据的总体分布范围来说，数值计算结果与实验一致，最大差距为14%。HBS外形在攻角小于14˚范围内时的俯仰阻尼导数变化不大，而大于14˚后俯仰阻尼导数迅速增加，这和图3俯仰力矩系数迟滞曲线所反映的内容是一致的。

图5给出了不同起始攻角下的俯仰刚度与实验结果的比较。数值计算结果显示，随着攻角的增大，静稳定性参数由正变负，在攻角为8度附近进入静稳定状态。随着攻角进一步增大，静导数在攻角为17度附近再次改变符号，静稳定性第二次发生改变。对比减缩频率分别为0.2和0.05的两组数据，较小的减缩频率与实验结果更为接近。当攻角为20˚时，缺乏实验数据来评估数值模拟的结果。

综合图4和图5中的两组计算结果，可以看到本文基于Etkin非定常气动力模型所发展的强迫简谐分

图4. 俯仰阻尼导数计算结果与实验比较

图5. 俯仰刚度计算结果与实验比较

析法辨识动态稳定性参数是可靠的，能够满足强迫振动非定常流场的计算要求。

5. “十”字翼导弹俯仰阻尼导数计算

图6给出了“十”字翼导弹的外形尺寸。导弹全长L = 10 m，直径d = 1 m。对这类外形复杂的带控制舵导弹外形，结构化的单块网格已经不能满足其网格生成要求。本文采用块间完全对接的多块结构化网格，对称面网格与分区情况见图7。网格总量为200万。近壁第一层网格到壁面的距离为 $1\times {10}^{-4}L$ 。给定强迫俯仰振动形式 $\alpha \left(t\right)=\theta ={\alpha }_0+{\alpha }_m\sin kt$ ，其中初始攻角 ${\alpha }_0={1.5}^{\circ }$ ，振幅 ${\alpha }_m={1.5}^{\circ }$ ，减缩频率k取0.05，质心位置 $X_{cg}/L=0.5$ 。

图6. “十”字翼导弹外形尺寸

图7. 对称面网格与分区图

从图8俯仰力矩系数的时间历程曲线可以看出，俯仰力矩系数收敛情况很好，从第二个周期开始已经完全达到谐振，随时间按正弦规律变化，非定常滞后效应比较明显。随着马赫数增加，同一时刻下的俯仰力矩系数减小。在同一周期内，不同马赫数的俯仰力矩系数均有相同的正弦振动形式，其最小值差别不大，最大值之间有较大差距。图9和图10给出了滚转力矩系数和偏航力矩系数的时间历程曲线。我们看到滚转力矩基本按正弦规律呈周期性变化，而偏航力矩呈不规则的脉动状态变化，无周期性规律，说明强迫俯仰振动时俯仰方向和滚转方向的交叉耦合性相对较强，本文对俯仰-滚转交叉耦合导数开展了数值辨识，其结果列于表1。

图8. 俯仰力矩系数时间历程曲线

图9. 滚转力矩系数时间历程曲线

图10. 偏航力矩系数时间历程曲线

“十”字翼导弹的俯仰力矩静、动导数均为负值，说明导弹是俯仰静、动稳定的。随着马赫数的增加，静、动导数的绝对值减小，静、动不稳定性减弱。同时，俯仰–滚转交叉耦合导数为正值，表明在俯仰过程中可能出现滚转方向的动不稳定问题。从表1还可以看出，俯仰–滚转阻尼导数量级很小，比俯仰力矩动导数小4个量级，计算时涉及表面摩阻等复杂的计算问题，增加了计算难度，因此交叉耦合导数计算的正确性和计算精度还需相关实验数据进一步加以验证。

图11进行了俯仰阻尼导数的计算结果与文献 [7] 中实验数据和文献 [8] 中基于准定常欧拉方程的求解结果的相互比较。从中可以看出，本文预测的俯仰阻尼导数比文献中采用准定常欧拉方程的求解结果更加接近风洞实验的数据，最大误差不超过6%。这里还需要说明的是，文献中采用定速率运动形式得到气动力时间历程，本文采用强迫正弦运动得到气动力时间历程，两种方法得到的结果基本相同。

图11. 俯仰阻尼导数计算结果与实验比较

表1. “十”字翼导弹动导数辨识结果

6. “十”字翼导弹滚转阻尼导数辨识

“十”字翼导弹强迫滚转振动的计算条件为：攻角0˚，侧滑角0˚，飞行高度1 km。马赫数分别为1.53，2.03，2.27，2.54。给定强迫滚转振动形式 $\psi ={\psi }_m\sin kt$ ，其中振幅 ${\psi }_m$ 取2.5˚，减缩频率取0.05，质心位置 $X_{cg}/L=0.5$ 。

图12为不同Ma下滚转力矩的迟滞环和时间历程曲线。从第二个周期开始，强迫滚转振动达到谐振，滚转力矩的迟滞环收敛成为标准的椭圆，其长短轴均与坐标轴平行。随马赫数的增大，同一周期内的相同时刻滚转力矩越小。

通过积分强迫滚转振动滚转力矩的迟滞环曲线，图13给出了滚转阻尼导数的辨识结果。随着马赫数的增加，滚转阻尼动导数的绝对值减小，动稳定性减弱，但没有符号上的变化，没有改变“十”字翼导弹的滚转动稳定性。本文将强迫滚转辨识出的动导数与文献 [9] 与 [10] 中的实验数据进行了比较。滚转阻尼导数的实验精度一般不高，主要由于实验中存在支架、洞壁干扰等不确定因素，这使对同一模型在不同风洞中的实验结果也存在一定的差异，本文计算得到的滚转阻尼导数介于两个实验结果之间，验证了方法和程序的可靠性。

非定常流动计算的精度依赖于三个重要的计算参数：时间步长、子迭代步数、网格规模。在强迫振动辨识动态稳定性参数的过程中，减缩频率也是影响辨识精度的重要参数，必须像马赫数、雷诺数一样考虑到与实验的相似性。本文选取了几组不同的参数值，分别考察减缩频率、子迭代步数、网格规模对“十”字翼导弹滚转阻尼导数辨识的影响。

图12. 不同Ma下滚转力矩的迟滞环和时间历程曲线

图13. 滚转阻尼导数计算结果与实验比较

表2. 不同计算参数对滚转阻尼导数的影响

减缩频率实际上是强迫振动频率快慢的反映，其大小会影响动导数的量值乃至符号。计算中减缩频率的大小应和实验一致才能得到可以相互对比的结果。但在实验中减缩频率一般很小，甚至达到10⁻⁴的量级，这在计算中产生的工作量是让人无法接受的。为考察减缩频率对动导数辨识结果的影响，本文在Ma = 1.58时选取了减缩频率为0.01、0.05、0.2、1.0四个数值，在其他条件相同的情况下分别进行计算，辨识结果列于表2。随着减缩频率的增大，强迫滚转振动的频率增加，滚转阻尼导数绝对值略有减小，滚转动稳定性减弱。但对于“十”字翼导弹这类细长体来说，滚转阻尼导数变化并不显著。减缩频率从0.01变化到1.00，滚转阻尼导数之间差别约为10%。

时间步长和亚迭代步数这两个参数的组合实际上决定了非定常流动的收敛程度。减小时间步长、加大亚迭代步数都可以提高非定常流动的收敛程度，从而提高动态稳定性参数的辨识精度，但同时计算量也会大幅上升。本文无量纲时间步长统一取0.01。子迭代步数分别取1，5，30。在Ma = 1.58条件下计算得到的动导数结果列于表2。可以看出，程序对非定常流场的计算精度较高，子迭代步数变化对滚转阻尼导数的影响不大。

对于速度导数(摩阻)、温度导数(热流)等导数量，网格的影响十分巨大，动导数也是一种物理量的导数，特别是与摩阻直接相关的滚转阻尼导数，根据以往钝锥的研究经验，网格规模的影响肯定不可忽略。本文计算“十”字翼导弹采用的是块间完全对接的多块结构化网格，网格总量为200万，近壁第一层网格到壁面的距离为1 × 10⁻⁴ L，网格规模属于中小型。为了研究网格规模特别是近壁第一层网格到壁面的距离对滚转阻尼导数辨识的敏感程度，本文将近壁第一层网格到壁面的距离稍作加密，达到6 × 10⁻⁵ L，同时网格总量增加到260万，辨识了Ma = 1.58条件下的滚转阻尼导数是−13.439，与未加密前的结果差别不大，这说明对“十”字翼导弹，之前的网格密度对非定常气动力和气动阻尼导数的计算来说已经足够，网格规模对滚转阻尼导数辨识不敏感，因此对于“十”字翼导弹，本文没有进一步开展网格规模的影响研究。

7. 结论

本文基于三维非定常NS方程，采用数值模拟的方法开展了超音速条件下“十”字翼导弹动态特性研究。主要结论如下：

1) 本文数值模拟了HBS模型和“十”字翼导弹小振幅强迫振动下的非定常气动力/力矩特性，对强迫运动过程得到的气动力迟滞曲线辨识了俯仰和滚转阻尼导数。通过与实验结果的比较验证了动导数辨识方法的可靠性。此外还开展了减缩频率、子迭代步数、网格规模对动导数辨识结果的影响分析。

2) 超音速条件下“十”字翼导弹俯仰运动时，俯仰–滚转交叉耦合导数有较为明显的周期性规律(但量级很小，与直接阻尼导数相比差4个量级)，而对于俯仰–偏航交叉耦合导数来说，其量值呈不规则的脉动状态变化，与俯仰振动的交叉耦合性很弱。交叉耦合导数计算的正确性还需做进一步的深入研究。

文章引用

刘 绪,王光辉,范化喜,张金春. 超音速条件下复杂带翼外形的强迫振动分析
Forced Vibration Analysis of Complex Configurations with Fins in Supersonic Environment[J]. 声学与振动, 2019, 07(01): 1-11. https://doi.org/10.12677/OJAV.2019.71001

参考文献