柔性弹箭动力学建模通常采用基于能量原理的拉格朗日法建立封闭的柔性弹箭气动弹性和飞行力学方程组，在小攻角和小变形条件下，对弹体弹性变形和大运动进行耦合^[1-6].其中，对弹箭振动特性和气动力参数的准确获取是柔性弹道计算和稳定性分析的关键.线性多体系统传递矩阵法特别适合于解决线性多体系统振动特性、正交性、动力响应等问题^[7]，兼具方便建模，计算速度快和精度高的优点.计算流体力学(CFD)方法可以很好的预测复杂几何体的气动特性参数和流体流动现象^[8-11]，但是如果将其应用到整个的封闭方程里则计算量太大，计算周期长，不利于快速计算和通过系统方程进行理论分析.而基于气动力系数导数的气动力计算方法则在小变形和小攻角范围内可以兼顾精度高和计算简单的特点.

本文采用多体系统传递矩阵法和CFD与经验公式相结合的方法，分别计算了火箭弹的振动特性和气动力系数导数，在此基础上计算弯曲火箭弹沿轴向单位长度气动力系数分布，并通过CFD方法对该方法的精度进行了验证.

1 模型建立 1.1 火箭弹结构建模

随着长细比的增加，弹体的抗弯刚度和横向振动的固有频率降低，弹体的横向弯曲为影响飞行力学的主要因素，因此大长细比火箭弹可看作自由-自由非均质Euler-Bernoulli梁^{[3, 12-13]}.如图 1所示，非均质梁进一步简化为N段均质梁依次连接而成^[7]，第i段梁的输入端和输出端在物理坐标下的状态矢量为Z_i-1,i和Z_i,i+1，其中Z_i-1,i=[y_i-1,iθ_zi-1,im_zi-1,iq_yi-1,i]^T.如图 2所示，y_i-1,i、θ_zi-1,i、m_zi-1,i和q_yi-1,i分别为物理坐标系下梁输入端的位移、转动角度、所受力矩和剪切应力，m-_i、l_i、EI_i和S_i分别为第i段梁的线密度、长度、抗弯刚度和所受轴向力(如果为拉力则为负值).

火箭弹轴上的每一点所受的轴向力大小不同，等于该点至火箭弹顶端部分受到的惯性力.底部所受轴向力大小等于推力，顶端所受轴向力为零.S_i取第i段质心处所受的轴向力^[3]为

${{S}_{i}}=\frac{{{F}_{T}}}{m}\left( m-\sum\limits_{j=1}^{i-1}{{{{\bar{m}}}_{j}}{{l}_{j}}-\frac{1}{2}{{{\bar{m}}}_{j}}{{l}_{j}}} \right).$

式中:F_T为末端推力;m为整个火箭弹的质量;l_j为第j段的长度;m_j为第j段的线密度.

受均匀轴向压力作用的均质Euler-Bernoulli梁横向自由振动微分方程为^[14]

$EI\frac{{{\partial }^{4}}y}{\partial {{x}^{4}}}+S\frac{{{\partial }^{2}}y}{\partial {{x}^{2}}}+\bar{m}\frac{{{\partial }^{2}}y}{\partial {{t}^{2}}}=0.$

推导可得长度为l的Euler-Bernoulli梁的传递矩阵为

$B\left( x \right)=\left[ \begin{matrix} cosh\left( Mx \right) & sinh\left( Mx \right) & cos\left( Nx \right) & sin\left( Nx \right) \\ Msinh\left( Mx \right) & Mcosh\left( Mx \right) & -Nsin\left( Nx \right) & Ncos\left( Nx \right) \\ EI{{M}^{2}}cosh\left( Mx \right) & EI{{M}^{2}}sinh\left( Mx \right) & -EI{{N}^{2}}cos\left( Nx \right) & -EI{{N}^{2}}sin\left( Nx \right) \\ \left( EI{{M}^{3}}-SM \right)sinh\left( Mx \right) & \left( EI{{M}^{3}}-SM \right)cosh\left( Mx \right) & \left( EI{{N}^{3}}+SN \right)sin\left( Nx \right) & -\left( EI{{N}^{3}}+SN \right)cos\left( Nx \right) \\ \end{matrix} \right],$

其中：

$M=\sqrt{\frac{\frac{{{S}^{2}}}{{{\left( EI \right)}^{2}}}+\frac{4\bar{m}{{\omega }^{2}}}{EI}-\frac{S}{EI}}{2}},N=\sqrt{\frac{\frac{{{S}^{2}}}{{{\left( EI \right)}^{2}}}+\frac{4\bar{m}{{\omega }^{2}}}{EI}-\frac{S}{EI}}{2}}.$

对于非均质梁，连接点处在模态坐标下的状态矢量为Z_0,1，Z_1,2，Z_2,3，…，Z_i-1,i，…，Z_n,n+1，每一段的传递方程为

${{Z}_{i,i+1}}={{U}_{i}}{{Z}_{i-1,i}},\left( i=1,2,\cdots ,n \right).$

(1)

因此总传递方程为

${{Z}_{n,n+1}}={{U}_{all}}{{Z}_{0,1}},$

其中，U_all为总传递矩阵U_all=U_nU_n-1…U_i…U₂U₁.

$\Delta =\left| \begin{matrix} {{U}_{all}}_{3,1} & {{U}_{all}}_{3,2} \\ {{U}_{all}}_{4,1} & {{U}_{all}}_{4,2} \\ \end{matrix} \right|=0,$

(2)

式(2)即为特征方程，解此方程可得频率ω.当频率ω_k(k为模态阶数)得到后，M_k、N_k、U_i-1,i和U_all都可以确定.通过边界条件可得初始端的状态矢量Z_0,1=B0a，同时非均质梁任何位置的状态矢量可由式(1)得到.

1.2 气动力模型

火箭弹弹身细长，翼面较小，并且翼面在弦长方向的尺寸占弹身总长很少，对于这样的柔性弹，在线性理论范围内可以将弹体沿轴向划分为若干微段，各段所受到的气动力为作用于该段压力中心的合力.在准定常理论和线性理论的假设下，不计气流的尾流影响，各段的气动力只和该段的当地攻角有关，认为变形后的弹体各段当地柔性攻角等于该段的合成速度向量与该段弦线之间的瞬时夹角，采用数值模拟方法计算出不同Ma数下的气动力系数导数的分布，可直接将气动力表述成当地柔性攻角的函数^{[5-6, 15-16]}.

弹体微段的当地柔性攻角侧滑角α(x,t)和β(x,t)分别为^[15]

$\begin{align} & \alpha x,t=\frac{w}{u}-\sum i\varphi {{\prime }_{i}}\left( x \right){{\zeta }_{i}}t+\frac{1}{u}\sum\limits_{i}{{{\varphi }_{i}}\left( x \right){{{\dot{\zeta }}}_{i}}\left( t \right)-\frac{qx}{u}}, \\ & \beta \left( x,t \right)=\frac{v}{u}-\sum\limits_{i}{\varphi {{\prime }_{i}}\left( x \right){{\eta }_{i}}t+\frac{1}{u}\sum\limits_{i}{{{\varphi }_{i}}\left( x \right){{{\dot{\eta }}}_{i}}\left( t \right)+\frac{rx}{u}}}. \\ \end{align}$

式中:u、v、w分别为速度在体坐标系中的分量；p、q、r分别为角速度在体坐标系中的分量；η_it、ζ_i(t)分别为广义坐标；φ_i(x)为振型函数.

基于线性理论和准定常理论，弹体所受整体的阻力、侧向力和升力分别为

$\begin{align} & {{F}_{D}}=-{{\int }_{L}}\frac{1}{2}\rho {{V}^{2}}S{{C}^{0}}_{D}\left( x \right)+\frac{1}{2}\rho {{V}^{2}}SC{{'}_{D}}\left( x \right)({{\alpha }^{2}}\left( x,t \right)+ \\ & {{\beta }^{2}}\left( x,t \right)dx, \\ & {{F}_{S}}=-{{\int }_{L}}\frac{1}{2}\rho {{V}^{2}}SC{{'}_{S}}\left( x \right)\beta x,tdx, \\ & {{F}_{L}}=-{{\int }_{L}}12\rho {{V}^{2}}SC{{'}_{L}}\left( x \right)\alpha x,tdx. \\ \end{align}$

式中：C⁰_D(x)、C^′_D(x)、C^′_S(x)、C^′_L(x)分别为火箭弹沿弹身轴向的单位长度零升阻力系数、阻力系数导数、侧向力系数导数、升力系数导数的分布.对于轴对称火箭弹，C^′_S(x)和C^′_L(x)相等.上述的气动力系数导数可通过气动力系数计算得到：

${{C}^{\prime }}_{L}={{C}^{\alpha }}_{L}\left( x \right)/\alpha ,$

(3)

${{C}^{\prime }}_{D}={{C}^{\alpha }}_{D}\left( x \right)-{{C}^{0}}_{D}\left( x \right)/{{\alpha }^{2}}.$

(4)

式中：C^α_L(x)、C^α_D(x)分别为刚体火箭弹在小攻角α下沿弹身的单位长度升力系数、阻力系数；C⁰_D(x)为刚体火箭弹沿弹身轴向的单位长度零升阻力系数，通过CFD对火箭弹进行数值模拟即可计算出需要的参数值.

2 计算结果 2.1 火箭弹振动特性

利用已推导的多体系统传递矩阵法公式计算某火箭弹被动段时的频率和振型，被动段时质量不变，推力为零.如图 3所示，当火箭弹所分段数超过50时1阶频率的计算结果没有明显的变化，说明将火箭弹简化为50段均质梁即可以保证计算的精度.图 4为该火箭弹前4阶振型，x/L=0为火箭弹尾部.

2.2 火箭弹气动力系数计算

由于火箭弹为轴对称，其法向和侧向具有相同的气动特性，因此本文中对侧向气动力不予考虑.如图 5所示，火箭弹为刚体模型，沿着轴向将火箭弹切成23段，每一个小段的中间点作为气动力的作用点，采用CFD方法计算沿轴向的气动力系数分布^[17-18].计算域来流条件为静压101 325 Pa，静温288 K，计算域外边界设为压力远场，即基于Rieman不变量的无反射边界条件，火箭弹物面设为无滑移壁面.采用SST k-ω湍流模型，基于密度求解器进行求解，离散格式采用AUSM格式和隐式求解格式，对流通量采用二阶迎风格式.本文采用多Block生成流场拓扑结构，外O Block生成弹体边界层的方法，生成了高质量的结构化网格，在边界层内对网格进行加密，保证y⁺<1，网格数为310万.如图 7所示，获得模型在α=0°和α=2°情况下的气动参数并与以往的实验数据比较，C_L和C_D的误差均在10%以内，验证了该数值计算方法的可行性与准确性.图 6和图 8分别是Ma=2，α=2°工况下的速度云图和火箭弹的升力、阻力系数分布.通过式(3)~(4)可得该火箭弹的气动力系数导数分布.

2.3 弯曲模型气动力计算方法验证

由于弹体变形基本基于1阶振型，因此本文基于1阶振型建立了弯曲的弹体模型.考虑到建模时可能存在尺寸误差及CFD计算误差，同时为了分析本文中气动力计算方法的适用范围，分别选择广义坐标ζ₁=0.025和ζ₁=0.100，对应的模型示意，如图 9所示.图 10为两个模型的弹体斜率，ζ₁=0.025和ζ₁=0.100远大于火箭弹实际正常飞行时的变形量，但是当弹体发生颤振或者弯曲发散时变形剧烈，很有可能达到这样的变形程度.

分别基于气动力系数导数和采用CFD方法计算了火箭弹在这两种弯曲变形下，Ma=2、不同攻角的气动力系数分布，并对计算结果进行对比.为了捕捉到微小变形对气流产生的影响，第1层网格高度需要足够小，此处设置y⁺<1，网格数为380万.对于ζ₁=0.025，即变形相对较小的情况下，攻角取0°、2°、4°时的单位长度阻力系数、单位长度升力系数分布的对比结果如图 11所示，可以发现两种方法结果吻合较好，阻力系数偏差在5%以内，升力系数偏差在10%以内.对于ζ₁=0.100，即变形相对较大的的情况下，攻角取0°、2°、4°时的单位长度阻力系数、单位长度升力系数分布的对比结果如图 12所示，在0°攻角时，两种方法的结果吻合较好；在2°攻角时弹的头部和尾部等局部阻力系数有一定的误差，但是误差较小；在4°攻角时，头部和尾部等局部误差增大.ζ₁=0.100的变形下在4°攻角对称面速度云图如图 13所示，由于变形增大，攻角增大，头部及弹身的尾流对火箭弹尾部的流动产生了复杂影响，此时，线性理论已经不再适用.

3 结 论

1) 推导了考虑轴向力作用的Euler-Bernoulli梁的传递矩阵，将某真实火箭弹简化成自由-自由的非均质Euler-Bernoulli梁，应用多体系统传递矩阵法计算了其频率和振型.

2) 采用基于SST k-ω湍流模型的CFD方法计算了该火箭弹在不同Ma数、不同攻角下的气动力系数，与实验数据相比最大误差不超过10%，证明了采用数值方法计算火箭弹气动力系数分布可以保证足够的精度以及采用数值方法验证弯曲火箭弹气动力建模的可行性.

3) 基于该火箭弹1阶振型建立了两种不同变形程度的弯曲弹体模型，在准定常理论和线性理论的假设条件下，结合气动力系数导数进行两个弯曲模型的气动力计算，得到了Ma=2，不同攻角两个模型的单位长度升力系数分布和阻力系数分布，并采用CFD方法对计算结果进行验证.通过验证发现，该火箭弹在实际飞行可能的变形范围内，基于准定常理论和线性理论的柔性弹箭气动力计算方法结果精度较好，在远超于实际的大变形时，该方法在弹头和弹身处仍具有一定精度，但尾部误差较大，计算精度降低.