大长细比火箭气动弹性分析①

2020-09-05龚学兵胡大庆王永平

固体火箭技术 2020年4期

龚学兵，胡大庆，王永平，郑佩

(中国航天科技集团有限公司四院四十一所，西安 710025)

0 引言

随着火箭总体设计技术的快速发展，火箭的长细比越来越大，相对弯曲刚度越来越小，箭体的气动载荷分布更加复杂，控制模型的非线性特征更强，火箭气动弹性控制问题越来越突出。气动弹性控制技术是未来火箭总体设计的关键技术之一，该项技术能有效提升火箭控制精度、降低结构的消极质量、提升弹道规划的能力，有利于分析总体参数拉偏范围对产品性能的影响以及飞行试验的危险因素来源。

针对导弹、飞机的气动伺服弹性问题，美国航空航天学会作了一系列有关气动弹性现象、气动弹性分析和气动弹性试验等方面的专题综述，并对气动伺服弹性研究的未来发展趋势提供指导[1]。目前，针对气动弹性问题，主要从理论分析和数值仿真计算、试验分析等三方面开展相关研究[2]。

理论分析方面，通过在飞行器动力学方程的基础上引入结构动力学模型，并将气动力通过有理函数(振型函数)进行广义非定常气动力进行拟合[3-7]，从而实现了通过理论的数学模型研究柔性变形、横向弯曲振动、旋转锁定等气动弹性问题的机理，并研究弹性变形对箭体姿态、控制系统设计的影响；而针对箭体姿态运动与弹性变形的耦合作用，部分作者[8-10]通过引入模态化方程组、大运动弹性体线性化处理、非均匀梁弯曲运动方程组等数学模型实现摆动与弹体振动耦合分析、惯量矩阵的时变特征分析、以及弹性导弹的侧向运动分析。虽然上述方法能够从理论角度给出弹性运动的精确结果，但非线性特征更加突出，导致控制系统设计存在困难。

数值仿真计算方面，采用计算流体力学程序计算非线性非定常气动力是当前跨声速气动伺服弹性系统建模研究的主要手段。气动伺服弹性的非线性特征导致流固耦合计算陷入开环系统模型阶数过高的问题，不利于控制率的设计。目前,相关学者[11-13]开展了相关模型降阶、提升计算效率、增强算法的鲁棒性和提供求解精度的相关学术研究。气动仿真[14-17]通过数字求解，能将复杂的模型通过分段线性化插值的手段，将非线性问题转为线性化问题，有利于准确的分析弹性变形对箭体的升力系数、阻力系数和俯仰力矩系数的影响，以及在不同马赫数和攻角下的气动力变化规律，同时进行结构强度校核与系统颤振速度分析，但上述方法的定量结果难以通过线性化的数学模型进行表征。

综上所述，虽然国、内外研究学者在气动弹性控制问题领域开展了大量的理论与仿真分析工作，但上述研究未从总体设计角度考虑，难以在线性化处理后的气动伺服弹性动力学模型的基础上，为控制系统设计提供气动参数、气动弹性系数的准确拉偏范围。工程上常用的气动伺服弹性分析方法是利用线性化的传递函数模型，分析弹箭开环系统的幅频、相频特性，开展控制系统的频域PID控制系统设计，通过添加陷幅滤波器抑制弹箭开环系统与飞行控制系统的不利耦合作用。而弹性箭体的气动参数具有非线性特征，不利于控制系统设计。本文通过理论分析部分提供控制系统所需的线性化弹性动力系数，并通过数值仿真技术为控制系统参数拉偏提供准确的变化范围，结合伯德图的稳定性分析结果，为气动弹性控制的总体设计提供指导。

1 未变形箭体的气动计算

某火箭的箭体结构与气动外形见图1，长细比达到16，气动弹性问题较为突出。在飞行过程中，各点的弹道和气动力参数均处于不断的变化之中，在气动弹性参数计算时不可能对全弹道各点的参数都进行考察。因此，本文给出了火箭在飞行过程中攻角取0°、2°、4°和马赫数分别为2、3、4、5状态组合下的气动参数。

利用ICEM CFD对上述理论化模型进行网格划分。其中，尾部流场区域选取箭体参考长度的7.5倍左右，头部区域选取箭体参考长度3.75倍左右，径向外流场取箭体参考长度的3倍左右(图2)。在Fluent中设置了壁面函数，在划分网格之前根据相似型号的工程经验确定第一层网格的厚度值(商业上可以用Pointwise与CFD-online计算yplus值，从而确定第一层网格厚度值。k-ε模型的yplus值为30～300之间)。本文的首层网格的厚度值设置为1 mm，附面层的厚度为15 mm，网格的厚度增长率取1.1(图3)。网格的无关性是通过对比多个设计师的仿真计算结果，确认网格数量和网格质量不影响气动参数的计算精度。

图2 未变形箭体的ICEMCF气动网格划分Fig.2 Aerodynamic grid generation of un- deformed rocket by ICEMCFD

图3 尾翼位置的局部网格及边界层网格示意图Fig.3 Aerodynamic grid generation of un- deformed rocket by ICEMCFD

气动计算的控制方程为雷诺平均Navier-Stokes(RANS)方程，空间离散格式为一阶迎风格式，采用双时间步推进进行求解。通过k-ε湍流模型和壁面函数法求解近壁面处边界层，远场边界条件为压力远场，箭体壁面边界条件为无滑移壁面。在开展箭体的气动参数仿真计算时，单个状态的计算需要迭代到5000步。针对高超声速状态，需要选取不同的Courant数，其默认值是5。在实际计算中，初始Courant设置为1，逐步增加为3、5、10，每次增加Courant参数后，在上一步计算结果的基础上重新进行气动仿真计算。每个Case的工况由马赫数、大气压共同决定。

气动参数计算结果如表1所示。

表1 未变形的轴向力系数Ca1、向力系数Cn1及俯仰力矩系数Mz1Table 1 Axial force coefficient(Ca1)，normal force coefficient(Cn1)and pitching moment coefficient(Mz1) of un-deformed rocket

2 变形箭体的理论气动参数计算

变形箭体的理论计算为控制系统设计提供线性化的理论模型，便于控制系统通过经典PID控制系统设计方法进行气动弹性控制系统设计。

2.1 弹性动力系数的理论计算公式

假设火箭弹性振动用广义坐标qi(t)表示，它是随时间变化的量，由下列二阶常微分方程确定：

(1)

式中ωi为火箭第i次振型的固有频率；ξi为火箭第i次振型的阻尼系数；Qi为火箭对应的第i次振型的广义力；Mi为火箭对应的第i次振型的广义质量。

(2)

(3)

式中fy1(x)为作用于火箭的外力在箭体轴y1上的投影；m(x)为沿箭体纵轴的弹体质量分布。

频率ωi和固有振型函数Wi(x)为壳体结构特性和沿纵轴质量分布的函数；Qi为影响弹体振动的广义力。由式(1)可知，只要知道fy1(x,t)即可求出其表达式，fy1(x,t)包括控制力、气动力、推力以及舵的惯性力等诸力沿箭体轴oy1的投影。

在未变形的气动参数基础上，进行变形箭体的气动参数计算。通过振型函数、模态质量进行气动参数的重新计算，主要依据如下：

(1)对于未变形的箭体

(2)对于变形后的箭体

除气动力外，上面提到的各力都作用在火箭箭体的某一位置，所以是集中力，由(1)和(2)可直接得到相应的广义力。

将广义力的具体表达式代入式(4)，略去各次振型之间的相互影响，即得到变形火箭在某次振型影响下的振动方程：

(4)

其中，

通过式(4)结合弹性振动的弹体小扰动方程，获取弹性振动下舵偏角到弹体姿态角的传递函数，结合伯德图的幅值裕度和相角裕度，即可完成弹性变形对火箭飞行稳定性的影响。

2.2 模态分析

模态分析主要用于提取振型，对实体几何模型优化，要求对大的圆孔进行保留，确保变形斜率的准确性。图4给出了实体结构模型的大尺寸的圆孔等特征保留的示意图。

图4 结构圆孔尺寸保留示意图Fig.4 Schematic diagram of structural with round hole size reservation

模型简化后，为了保证质量分布符合实体结构的质量分布特征，必须在模型中加入实体模型的其他附属件的质量，质量分布如图5所示。在质量分布设置过程中，应该根据质量分布以及刚度影响等因素，设置相应的点质量和面质量。

图5 箭体质量分布示意图Fig.5 Diagram of mass distribution of rocket

全箭模态仿真结果如表2所示，同一阶弯曲模态包含两个正交方向上的模态振型。由于Y、Z方向的质量分布存在差异，导致不同象限(I、III与II、IV)的频率分布存在一定的差异性，该频率计算结果与相同型号产品地面试验结果修正后的频率偏差在2 Hz范围以内。通常振型斜率应由速率陀螺进行测量，但由于试验条件限制，当前只能以理论分析为主。

表2 箭体模态仿真结果Table 2 Modal analysis results of rocket body Hz

2.3 振型函数的提取

根据上述结构模态的分析结果，采用APDL语言，提取弹体坐标位置以及相应的振型幅值，图6为箭体振型随着弹体位置变化的振型幅值函数图。坐标原点选取弹体顶点作为参考点。

(a) Ma=2 (b) Ma=3

2.4 基于理论分析的气动升力对比分析

表3 一阶振型修正的法向力系数斜率Table 3 Slope of normal force revised by first bend mode

对比分析计算结果可知，由于马赫数的增加，导致理论计算的弹体气动升力斜率逐渐加大。这是由于马赫数的增加，导致一阶振型的幅值变大，振型引起的附加攻角增大，变形箭体的法向力系数斜率逐渐增大。

3 变形箭体的气动对比分析

虽然上述理论分析可以用于气动弹性控制的动力系数计算，但变形后的箭体气动特性与未变形箭体的气动特性存在一定的差异性，例如：弹体压心前移、气动系数不对称等，而且线性的方程与非线性气动弹性方程求解也存在一定的误差，需要开展风洞试验进行修正，而风洞试验费用昂贵。本文通过Fluent仿真模拟试验为理论计算结果提供拉偏范围，用于降低气动伺服弹性控制系统的设计风险，并给控制系统的弹道拉偏提供合理的变化范围。

根据弹体变形的结果，获得弹体一阶变形的振型函数。忽略弹翼、舵面的变形。变形后的几何模型如图7所示。

图7 变形后的箭体模型Fig.7 Geometric model of deformed rocket

根据变形后的弹体通过ICEM CFD软件划分变形火箭流场网格，如图8所示。

图8 变形箭体的气动网格划分结果Fig.8 Aerodynamic meshing result of the deformed rocket

计算获得火箭在一阶变形后的气动阻力、升力、俯仰力矩系数(表4)，并与理论计算结果进行对比。

对比表1和表4可以看出：

表4 变形的轴向力系数Ca1、向力系数Cn1及俯仰力矩系数Mz1Table 4 Axial force coefficient(Ca1)，normal force coefficient(Cn1)and pitching moment coefficient(Mz1) of deformed rocket

(1)弹体在同一马赫数条件下的阻力会逐渐减小；

(2)零度攻角下弹体的升力不为零，这是由于附加攻角导致弹体在理论0°攻角下仍然存在升力；

(3)零度攻角下弹体的俯仰力矩不为零，这是由于附加攻角导致弹体在理论0°攻角下仍然存在升力，进而产生俯仰力矩。

箭体载荷设计一般考虑极限条件下，本文选取4°攻角作为载荷设计依据，如表5所示。由表5可以看出，由于弹性变形，箭体的轴向力矩系数变小，而法相力矩系数和俯仰力矩系数变大。最大偏差分别为-15.59%、10.34%和16.91%。在箭体载荷设计时应考虑拉弯组合效应。同时在弹道拉偏过程中，应考虑箭体弹性变形引起的气动参数变化。

此外，表6为变形前后的压心系数，通过对比可以看出弹体压心向前移动。随着马赫数增加，在4°攻角条件下，弹体压心前移最大为13.4%，压心的变化对控制系统的稳定性设计应着重考虑。

表6 变形前、后的压心系数Xp(α)Table 6 Pressure center coefficient(Xp(α))of undeformed rocket and deformed rocket

计算获得弹体在一阶变形后的气动参数，并与未变形的气动参数、气弹理论计算的结果对比分析，如表7所示。可以看出，由模态叠加法获得的弹体升力斜率随着马赫数增加而变小，理论计算结果与Fluent仿真计算结果在Ma=5条件下出现较大的偏差，达到14%。因此，在弹性动力系数计算时，应提供14%左右的拉偏量。

表7 变形火箭一阶振型修正的法向力系数斜率的理论计算结果与Fluent计算结果对比Table 7 Slope of normal force coefficient revised by first bend mode of the deformed rocket with the theoretical and Fluent results

4 静气弹稳定性分析

通过式(2)～式(4)，结合表7给出的参数，即可完成弹性动力学系数的计算。

表8 一级飞行横向运动方程式系数Table 8 Coefficient of lateral motion equation of first flight

在弹性动力系数计算的基础上，结合弹体的传递函数，即可通过伯德图开展气弹性稳定裕度分析(图9和表9)。由图9可见，一阶振型频率对应的幅值曲线处于0 dB以下；由表9可见，一阶弹性幅值裕度大于0 dB以上，满足弹性运动幅值稳定，系统相位滞后也在可以接受的范围内。通过引入陷波滤波器，实现对箭体弹性振动的主动抑制。在拉偏下限状态下，系统依然能够满足弹性幅值运动稳定。