1. 引言
动态稳定性参数(动导数)是飞行器控制系统设计,飞行器轨道设计和飞行器纵、横向动态稳定性分析的重要参数 [1] 。目前,评估飞行器飞行品质和动态特性的方法主要有飞行实验、风洞实验和理论计算,其重要内容之一是获取动导数。飞行实验和风洞实验确定动态参数难度大、周期长,且实验状态受实验条件的诸多限制 [2] 。从国内外发展来看,通过数值计算获得动导数是动态参数研究发展的重要趋势 [3] 。
动导数计算方法分为:1) 小扰动线化理论、牛顿理论与气动力工程模型相结合的动导数工程近似方法;2) 模拟动态实验的数值自由振荡法和数值强迫振荡法。动导数工程近似方法在90年代以前被国内外普遍采用,这是一种经验和半经验的方法,如内伏牛顿流理论等。工程近似方法的最大特点在于快捷,但依赖于经验性,考虑气流分离、再附和尾流等效应较为困难,对于超音速复杂外形流动难以适应,计算精度难以保障。因此考虑了流动非线性的模拟动态实验的全数值计算是动导数研究发展的重要方向 [4] 。
本文采用数值模拟强迫振荡法辨识动态稳定性参数,在Etkin非定常气动力模型的基础上,建立了单自由度俯仰及滚转运动时的非定常气动力/力矩系数和状态变量之间的数学模型,给出了强迫简谐分析方法求解俯仰阻尼导数及滚转阻尼导数的数值计算方法。此外,为了考虑俯仰引起的纵横向交叉耦合效应,本文还开展了“十”字翼导弹强迫俯仰振动时的交叉耦合导数计算。
2. 数值方法
在笛卡尔惯性坐标系下,对完全气体、忽略质量力下的三维无量纲Navier-Stokes方程形式如下:
(1)
式中:U为守恒变量,E、F、G为无粘通量,Ev、Fv、Gv为粘性通量。本文在空间上采用二阶精度的Roe格式,Roe格式是基于Godunov方法的基本思路发展起来的,是通量差分分裂(FDS)格式的一种,由于其优秀的间断分辨率和较小的数值耗散性,目前受到广泛的应用。
采用引入“双时间步”(Dual-time-step)方法的LU-SGS隐式格式离散流体运动方程时间项。动网格生成采用刚性网格生成方法。远场入流边界采用适用于动态边界条件下的一维Riemann不变量的无反射边界条件。对于定常/非定常超音速出流,边界点上的值由内流场计算结果外推得到。对于壁面边界,采用速度无滑移、绝热壁条件。奇性轴采用外插后周向平均处理。
3. 动导数计算方法
对于飞行器做单自由度的俯仰运动,如果其质心速度不变,则确定运动的独立状态变量只有攻角
和俯仰角速度q。对俯仰力矩系数Cm,由Etkin假定,俯仰力矩系数可写成 [5] :
(2)
给定强迫振动:
(3)
式中
为瞬时攻角,
为俯仰角,
为基准状态的攻角,
为振动幅度,k为减缩频率。
在k不很大,且忽略高阶导数的影响时,当非定常振动过程达到谐振解,通过数值积分可以求出俯仰刚度
和俯仰阻尼导数
:
(4)
(5)
为积分起始时刻,
为振动周期。同样的,强迫俯仰振动时引起的俯仰-滚转力矩交叉耦合导数为:
(6)
与之相仿,给定强迫滚转振动:
(7)
其中
为滚转角,
为滚转角振幅。由此辨识出的滚转阻尼导数为:
(8)
4. 验证算例
本文采用HBS(Hyper Ballistic Shape)外形作为验证算例 [6] ,HBS外形是高超声速导弹外形的标模。HBS外形参数见图1,其中:
。图2为HBS模型对称面网格,网格数:
(流向 × 周向 × 径向)。
![](//html.hanspub.org/file/1-2820067x31_hanspub.png)
Figure 2. The symmetry plane grid of HBS
图2. HBS对称面网格
![](//html.hanspub.org/file/1-2820067x32_hanspub.png)
Figure 3. Lagging curves of pitch moment vs attack angle
图3. 俯仰力矩系数迟滞曲线
计算条件为:马赫数
,以头部直径为参考长度的
,质心位置
,强迫俯仰振荡的振幅为1˚,减缩频率分别取0.2和0.05。图4是不同起始攻角下对俯仰阻尼导数的计算结果与文献 [6] 中的实验结果的比较。实验结果的误差带为±15%。从计算的总体情况来看:数值结果和实验结果在趋势上是一致的,都处于俯仰动态稳定状态。从数据的总体分布范围来说,数值计算结果与实验一致,最大差距为14%。HBS外形在攻角小于14˚范围内时的俯仰阻尼导数变化不大,而大于14˚后俯仰阻尼导数迅速增加,这和图3俯仰力矩系数迟滞曲线所反映的内容是一致的。
图5给出了不同起始攻角下的俯仰刚度与实验结果的比较。数值计算结果显示,随着攻角的增大,静稳定性参数由正变负,在攻角为8度附近进入静稳定状态。随着攻角进一步增大,静导数在攻角为17度附近再次改变符号,静稳定性第二次发生改变。对比减缩频率分别为0.2和0.05的两组数据,较小的减缩频率与实验结果更为接近。当攻角为20˚时,缺乏实验数据来评估数值模拟的结果。
综合图4和图5中的两组计算结果,可以看到本文基于Etkin非定常气动力模型所发展的强迫简谐分
![](//html.hanspub.org/file/1-2820067x36_hanspub.png)
Figure 4. Comparison of computational and experimental pitch damping derivative
图4. 俯仰阻尼导数计算结果与实验比较
![](//html.hanspub.org/file/1-2820067x37_hanspub.png)
Figure 5. Comparison of computational and experimental pitching moment slope
图5. 俯仰刚度计算结果与实验比较
析法辨识动态稳定性参数是可靠的,能够满足强迫振动非定常流场的计算要求。
5. “十”字翼导弹俯仰阻尼导数计算
图6给出了“十”字翼导弹的外形尺寸。导弹全长L = 10 m,直径d = 1 m。对这类外形复杂的带控制舵导弹外形,结构化的单块网格已经不能满足其网格生成要求。本文采用块间完全对接的多块结构化网格,对称面网格与分区情况见图7。网格总量为200万。近壁第一层网格到壁面的距离为
。给定强迫俯仰振动形式
,其中初始攻角
,振幅
,减缩频率k取0.05,质心位置
。
![](//html.hanspub.org/file/1-2820067x43_hanspub.png)
Figure 6. Shape of missile geometries with fins
图6. “十”字翼导弹外形尺寸
![](//html.hanspub.org/file/1-2820067x44_hanspub.png)
Figure 7. The symmetry plane grid and Domain decomposition
图7. 对称面网格与分区图
从图8俯仰力矩系数的时间历程曲线可以看出,俯仰力矩系数收敛情况很好,从第二个周期开始已经完全达到谐振,随时间按正弦规律变化,非定常滞后效应比较明显。随着马赫数增加,同一时刻下的俯仰力矩系数减小。在同一周期内,不同马赫数的俯仰力矩系数均有相同的正弦振动形式,其最小值差别不大,最大值之间有较大差距。图9和图10给出了滚转力矩系数和偏航力矩系数的时间历程曲线。我们看到滚转力矩基本按正弦规律呈周期性变化,而偏航力矩呈不规则的脉动状态变化,无周期性规律,说明强迫俯仰振动时俯仰方向和滚转方向的交叉耦合性相对较强,本文对俯仰-滚转交叉耦合导数开展了数值辨识,其结果列于表1。
![](//html.hanspub.org/file/1-2820067x45_hanspub.png)
Figure 8. Curves of pitch moment vs time
图8. 俯仰力矩系数时间历程曲线
![](//html.hanspub.org/file/1-2820067x46_hanspub.png)
Figure 9. Curves of roll moment vs time
图9. 滚转力矩系数时间历程曲线
“十”字翼导弹的俯仰力矩静、动导数均为负值,说明导弹是俯仰静、动稳定的。随着马赫数的增加,静、动导数的绝对值减小,静、动不稳定性减弱。同时,俯仰–滚转交叉耦合导数为正值,表明在俯仰过程中可能出现滚转方向的动不稳定问题。从表1还可以看出,俯仰–滚转阻尼导数量级很小,比俯仰力矩动导数小4个量级,计算时涉及表面摩阻等复杂的计算问题,增加了计算难度,因此交叉耦合导数计算的正确性和计算精度还需相关实验数据进一步加以验证。
图11进行了俯仰阻尼导数的计算结果与文献 [7] 中实验数据和文献 [8] 中基于准定常欧拉方程的求解结果的相互比较。从中可以看出,本文预测的俯仰阻尼导数比文献中采用准定常欧拉方程的求解结果更加接近风洞实验的数据,最大误差不超过6%。这里还需要说明的是,文献中采用定速率运动形式得到气动力时间历程,本文采用强迫正弦运动得到气动力时间历程,两种方法得到的结果基本相同。
![](//html.hanspub.org/file/1-2820067x48_hanspub.png)
Figure 11. Comparison of computational and experimental pitch damping derivative
图11. 俯仰阻尼导数计算结果与实验比较
![](Images/Table_Tmp.jpg)
Table 1. Computation results of dynamic derivatives for missile geometries with fins
表1. “十”字翼导弹动导数辨识结果
6. “十”字翼导弹滚转阻尼导数辨识
“十”字翼导弹强迫滚转振动的计算条件为:攻角0˚,侧滑角0˚,飞行高度1 km。马赫数分别为1.53,2.03,2.27,2.54。给定强迫滚转振动形式
,其中振幅
取2.5˚,减缩频率取0.05,质心位置
。
图12为不同Ma下滚转力矩的迟滞环和时间历程曲线。从第二个周期开始,强迫滚转振动达到谐振,滚转力矩的迟滞环收敛成为标准的椭圆,其长短轴均与坐标轴平行。随马赫数的增大,同一周期内的相同时刻滚转力矩越小。
通过积分强迫滚转振动滚转力矩的迟滞环曲线,图13给出了滚转阻尼导数的辨识结果。随着马赫数的增加,滚转阻尼动导数的绝对值减小,动稳定性减弱,但没有符号上的变化,没有改变“十”字翼导弹的滚转动稳定性。本文将强迫滚转辨识出的动导数与文献 [9] 与 [10] 中的实验数据进行了比较。滚转阻尼导数的实验精度一般不高,主要由于实验中存在支架、洞壁干扰等不确定因素,这使对同一模型在不同风洞中的实验结果也存在一定的差异,本文计算得到的滚转阻尼导数介于两个实验结果之间,验证了方法和程序的可靠性。
非定常流动计算的精度依赖于三个重要的计算参数:时间步长、子迭代步数、网格规模。在强迫振动辨识动态稳定性参数的过程中,减缩频率也是影响辨识精度的重要参数,必须像马赫数、雷诺数一样考虑到与实验的相似性。本文选取了几组不同的参数值,分别考察减缩频率、子迭代步数、网格规模对“十”字翼导弹滚转阻尼导数辨识的影响。
(a) 时间历程曲线(b) 迟滞环
Figure 12. Roll moment vs attack angle and time among different Mach numbers
图12. 不同Ma下滚转力矩的迟滞环和时间历程曲线
![](//html.hanspub.org/file/1-2820067x56_hanspub.png)
Figure 13. Comparison of computational and experimental roll damping derivative
图13. 滚转阻尼导数计算结果与实验比较
![](Images/Table_Tmp.jpg)
Table 2. Influence of different parameters on the roll damping derivative
表2. 不同计算参数对滚转阻尼导数的影响
减缩频率实际上是强迫振动频率快慢的反映,其大小会影响动导数的量值乃至符号。计算中减缩频率的大小应和实验一致才能得到可以相互对比的结果。但在实验中减缩频率一般很小,甚至达到10−4的量级,这在计算中产生的工作量是让人无法接受的。为考察减缩频率对动导数辨识结果的影响,本文在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−4 L,网格规模属于中小型。为了研究网格规模特别是近壁第一层网格到壁面的距离对滚转阻尼导数辨识的敏感程度,本文将近壁第一层网格到壁面的距离稍作加密,达到6 × 10−5 L,同时网格总量增加到260万,辨识了Ma = 1.58条件下的滚转阻尼导数是−13.439,与未加密前的结果差别不大,这说明对“十”字翼导弹,之前的网格密度对非定常气动力和气动阻尼导数的计算来说已经足够,网格规模对滚转阻尼导数辨识不敏感,因此对于“十”字翼导弹,本文没有进一步开展网格规模的影响研究。
7. 结论
本文基于三维非定常NS方程,采用数值模拟的方法开展了超音速条件下“十”字翼导弹动态特性研究。主要结论如下:
1) 本文数值模拟了HBS模型和“十”字翼导弹小振幅强迫振动下的非定常气动力/力矩特性,对强迫运动过程得到的气动力迟滞曲线辨识了俯仰和滚转阻尼导数。通过与实验结果的比较验证了动导数辨识方法的可靠性。此外还开展了减缩频率、子迭代步数、网格规模对动导数辨识结果的影响分析。
2) 超音速条件下“十”字翼导弹俯仰运动时,俯仰–滚转交叉耦合导数有较为明显的周期性规律(但量级很小,与直接阻尼导数相比差4个量级),而对于俯仰–偏航交叉耦合导数来说,其量值呈不规则的脉动状态变化,与俯仰振动的交叉耦合性很弱。交叉耦合导数计算的正确性还需做进一步的深入研究。