0 引言
国内外学者在旋转弹动稳定性方面进行了大量研究,Murphy[10-11]在弱非线性假设的前提下获得考虑非线性气动力矩的锥形运动解析解;Schiff[12]提出了一种锥形运动稳态、无粘流场的数值求解方法;Platus[13-14]阐述了马格努斯力矩产生的机理及气动偏心和质量偏心、内部质量特性变化等是导致旋转弹动态失稳的因素;闫晓勇[15]通过引入广义章动角和进动角速率将非线性锥形运动数学模型进行线性化,获得了稳定锥形运动时临界转速的解析解;王华毕[16]通过将非线性动力学方程线性化,获得侧向力作用下的稳定性判据,并给出锥形运动稳定的临界转速和临界锥角。陈尔康等[17,18]将喷气阻尼力矩和推力偏斜引入锥形运动方程,并推导得出了喷气阻尼和推力偏斜对角运动的影响,给出了稳定性判据。
本文通过建立旋转弹短周期动力学模型,得到由欧拉角表示的锥形运动方程,通过引入复数攻角获得锥形运动方程的解析解,给出了陀螺效应、气动静力矩作用下弹体的锥形运动形态,并推导了加入马格努斯效应、阻尼效应作用下的锥形运动方程解析解,同时对稳定性判据进行了推导,最后分析了自旋转速、静稳定性对锥形运动稳定性的影响。
1 建模
1.1 动力学方程
在准弹体坐标系下弹体绕质心转动的动力学方程为[19]:
$\left\{\begin{array}{l}{J}_{x4}\frac{d{\omega }_{x4}}{dt}+({J}_{z4}-{J}_{y4}){\omega }_{z4}{\omega }_{y4}={M}_{sx}\\ {J}_{y4}\frac{d{\omega }_{y4}}{dt}+({J}_{x4}-{J}_{z4}){\omega }_{x4}{\omega }_{z4}+{J}_{z4}{\omega }_{z4}\stackrel{ ·}{\gamma }={M}_{sy}\\ {J}_{z4}\frac{d{\omega }_{z4}}{dt}+({J}_{y4}-{J}_{x4}){\omega }_{y4}{\omega }_{x4}-{J}_{y4}{\omega }_{y4}\stackrel{ ·}{\gamma }={M}_{sz}\end{array}\right.$
式中Jx4,Jy4,Jz4分别为转动惯量在准弹体坐标系下的分量,Msx,Msy,Msz为气动力矩相对准弹体坐标系的分量,ωx4,ωy4,ωz4为旋转角速度在准弹体坐标系的分量。由于弹体姿态角的变化可认为是短周期的动态变化,因此建立短周期动力学模型,即在一小段飞行时间内,认为速度矢量不变,有
$\left\{\begin{array}{l}\stackrel{ ·}{\gamma }={\omega }_{x4}\\ \stackrel{ ·}{\psi }=\stackrel{ ·}{\beta }{ }^{*}={\omega }_{y4}\\ \stackrel{ ·}{\vartheta }=\stackrel{ ·}{\alpha }{ }^{*}={\omega }_{z4}\end{array}\right.$
式中$\stackrel{ ·}{\gamma }$、$\stackrel{ ·}{\psi }$、$\stackrel{ ·}{\vartheta }$分别为滚转角速度、偏航角速度和俯仰角速度,α*和β*为准攻角和准侧滑角。气动力矩Msx,Msy,Msz可以表示为:
$\left\{\begin{array}{l}{M}_{sx}={M}_{x}^{{\sigma }_{x}}{\sigma }_{x}+{M}_{x}^{{\omega }_{x4}}{\omega }_{x4}\\ {M}_{sy}={M}_{y}^{{\beta }^{*}}{\beta }^{*}+{M}_{y}^{{\omega }_{y4}}{\omega }_{y4}+M″{ }_{y}\left(\frac{\stackrel{ ·}{\gamma }D}{V}\right){\alpha }^{*}\\ {M}_{sz}={M}_{z}^{{\alpha }^{*}}{\alpha }^{*}+{M}_{z}^{{\omega }_{z4}}{\omega }_{z4}+M″{ }_{z}\left(\frac{\stackrel{ ·}{\gamma }D}{V}\right){\beta }^{*}\end{array}\right.$
式中${M}_{y}^{{\beta }^{*}}$、${M}_{z}^{{\alpha }^{*}}$为气动稳定力矩系数,${M}_{x}^{{\omega }_{x4}}$、${M}_{y}^{{\omega }_{y4}}$、${M}_{z}^{{\omega }_{z4}}$为气动阻尼力矩系数,M″y、 M″z为马格努斯力矩系数。
对于轴对称导弹,令Jx4=Jx为极转动惯量,Jy4=Jz4=Je为赤道转动惯量。将式(2)、(3)代入(1)中,可得旋转导弹横向转动的动力学方程为:
$\left\{\begin{array}{l}{J}_{e}\stackrel{ ··}{\beta }{ }^{*}+{J}_{x}\stackrel{ ·}{\gamma }\stackrel{ ·}{\alpha }{ }^{*}={m}_{y}^{{\beta }^{*}}qSL{\beta }^{*}+{{m}_{y4}^{\omega }}_{y}qSL{\omega }_{y4}+m″{ }_{y}qSL\left(\frac{\stackrel{ ·}{\gamma }D}{V}\right){\alpha }^{*}\\ {J}_{e}\stackrel{ ··}{\alpha }{ }^{*}-{J}_{x}\stackrel{ ·}{\gamma }\stackrel{ ·}{\beta }{ }^{*}={m}_{z}^{{\alpha }^{*}}qSL{\alpha }^{*}+{m}_{z}^{{\omega }_{z4}}qSL{\omega }_{z4}+m″{ }_{z}qSL\left(\frac{\stackrel{ ·}{\gamma }D}{V}\right){\beta }^{*}\end{array}\right.$
由于阻尼力矩${M}_{y}^{{\omega }_{y4}}$ωy4、${M}_{z}^{{\omega }_{z4}}$ωz4分别与ωy4、ωz4符号相反,若令${m}_{y}^{{\omega }_{y4}}$、${m}_{z}^{{\omega }_{z4}}$为正值,则式(4)右端阻尼力矩可用负号表示。对于静稳定弹体,正攻角产生负偏航方向的马格努斯力矩、正侧滑产生正俯仰方向的马格努斯力矩,若令m″y、m″z为正,则式(4)右端马格努斯力矩项同样可引入负号进行表示。
对于无控轴对称旋转弹,有:
$\left\{\begin{array}{l}{m}_{y}^{{\beta }^{*}}={m}_{z}^{{\alpha }^{*}}=m{\delta }_{s}\\ {m}_{y}^{{\omega }_{y4}}={m}_{z}^{{\omega }_{z4}}={m}_{d}^{{\omega }_{b}}\\ m″{ }_{y}=m″{ }_{z}=m″{ }_{m}\end{array}\right.$
式中${m}_{s}^{\delta }$、${m}_{d}^{{\omega }_{b}}$、m″m分别为章动静稳定力矩系数对章动角δ的导数、章进动阻尼力矩系数对ωb的导数、马格努斯力矩系数对弹体无量纲旋转角速度$\stackrel{ ·}{\gamma }$D/V和章动角δ的二阶偏导数。将式(5)代入(4),同时根据前述整理各力矩项的符号,可得由欧拉角表示的弹体锥形运动方程:
$\left\{\begin{array}{l}\stackrel{ ··}{\beta }{ }^{*}=\frac{qSL}{{J}_{e}}{m}_{s}^{\delta }{\beta }^{*}-\frac{qSL}{{J}_{e}}{m}_{d}^{{\omega }_{d}}{\stackrel{ ·}{\beta }}^{*}-\frac{{J}_{x}}{{J}_{e}}\stackrel{ ·}{\gamma }{\stackrel{ ·}{\alpha }}^{*}-\frac{qSL}{{J}_{e}}m″{ }_{m}\left(\frac{\stackrel{ ·}{\gamma }D}{V}\right){\alpha }^{*}\\ \stackrel{ ··}{\alpha }{ }^{*}=\frac{qSL}{{J}_{e}}{m}_{s}^{\delta }{\alpha }^{*}-\frac{qSL}{{J}_{e}}{m}_{d}^{{\omega }_{d}}{\stackrel{ ·}{\alpha }}^{*}+\frac{{J}_{x}}{{J}_{e}}\stackrel{ ·}{\gamma }{\stackrel{ ·}{\beta }}^{*}+\frac{qSL}{{J}_{e}}m″{ }_{m}\left(\frac{\stackrel{ ·}{\gamma }D}{V}\right){\beta }^{*}\end{array}\right.$
1.2 初始扰动对锥形运动稳定性的影响分析
将复数变量Δ=β*+iα*代入式(6),不考虑马格努斯项和阻尼项的影响,则式(6)可以表示为:
$\stackrel{ ··}{\Delta }$-i$\frac{{J}_{x}\stackrel{ ·}{\gamma }}{{J}_{e}}\stackrel{ ·}{\Delta }$-$\left(\frac{{J}_{x}\stackrel{ ·}{\gamma }}{2{J}_{e}}\right)$2$\frac{1}{{S}_{g}}$Δ=0
式中Sg为陀螺稳定因子,其表达式为:
Sg=$\frac{\left({J}_{x}\stackrel{ ·}{\gamma }\right)2}{4qSL{J}_{e}m{\delta }_{s}}$
式(7)所描述的系统的特征方程为:
$\stackrel{ ··}{\lambda }$-i$\frac{{J}_{x}\stackrel{ ·}{\gamma }}{{J}_{e}}\stackrel{ ·}{\lambda }$-$\left(\frac{{J}_{x}\stackrel{ ·}{\gamma }}{2{J}_{e}}\right)$2$\frac{1}{{S}_{g}}$λ=0
特征方程(9)的根为:
λ1,2=i$\frac{{J}_{x}\stackrel{ ·}{\gamma }}{2{J}_{e}}\left(1\pm \sqrt[ ]{1-\frac{1}{{S}_{g}}}\right)$
令ω1、ω2分别表示运动模态中的快、慢频率:
$\left\{\begin{array}{l}{\omega }_{1}=\frac{{J}_{x}\stackrel{ ·}{\gamma }}{2{J}_{e}}\left(1+\sqrt[ ]{1-\frac{1}{{S}_{g}}}\right)\\ {\omega }_{2}=\frac{{J}_{x}\stackrel{ ·}{\gamma }}{2{J}_{e}}\left(1-\sqrt[ ]{1-\frac{1}{{S}_{g}}}\right)\end{array}\right.$
则式(7)的通解可表示为:
Δ=C1${e}^{i{\omega }_{1}t}$+C2${e}^{i{\omega }_{2}t}$=ρ1cos(ω1t+ε1)+ρ2cos(ω2t+ε2)+i[ρ1sin(ω1t+ε1)+ρ2sin(ω2t+ε2)]
根据复数变量的定义可得:
$\left\{\begin{array}{l}{\beta }^{*}={\rho }_{1}cos({\omega }_{1}t+{\epsilon }_{1})+{\rho }_{2}cos({\omega }_{2}t+{\epsilon }_{2})\\ {\alpha }^{*}={\rho }_{1}sin({\omega }_{1}t+{\epsilon }_{1})+{\rho }_{2}sin({\omega }_{2}t+{\epsilon }_{2})\end{array}\right.$
式(13)包含四个待定系数ρ1、ρ2、ε1、ε2,假设弹体在t=0初始时刻存在扰动:α*=${\alpha }_{0}^{*}$、β*=${\beta }_{0}^{*}$、${\stackrel{ ·}{\alpha }}^{*}$=${\stackrel{ ·}{\alpha }}_{0}^{*}$、${\stackrel{ ·}{\beta }}^{*}$=${\stackrel{ ·}{\beta }}_{0}^{*}$,代入式(13)可得系数ρ1、ρ2、ε1、ε2,进而得到弹体的锥形运动轨迹,直接反映旋转弹章、进动耦合的运动形态。
2 马格努斯效应和阻尼效应对锥形运动稳定性的影响分析
考虑陀螺效应、气动静力矩、马格努斯效应和阻尼效应等全部气动力项时,同样将复数变量Δ=β*+iα*代入式(6),则式(6)可表示为:
$\stackrel{ ··}{\Delta }$-$\left(-\frac{qSL{m}_{d}^{{\omega }_{b}}}{{J}_{e}}+i\frac{{J}_{x}\stackrel{ ·}{\gamma }}{{J}_{e}}\right)\stackrel{ ·}{\Delta }$-$\left(\frac{qSL{m}_{s}^{\delta }}{{J}_{e}}+i\frac{qSLm″{ }_{m}\stackrel{ ·}{\gamma }}{{J}_{e}}\right)$Δ=0
式(14)所描述的系统特征方程的根为:
λ1,2=$\frac{1}{2}\left(-\frac{qSL{m}_{d}^{{\omega }_{b}}}{{J}_{e}}+i\frac{{J}_{x}\stackrel{ ·}{\gamma }}{{J}_{e}}\right)$±$\frac{1}{2}\sqrt[ ]{{\left(-\frac{qSL{m}_{d}^{{\omega }_{b}}}{{J}_{e}}+i\frac{{J}_{x}\stackrel{ ·}{\gamma }}{{J}_{e}}\right)}^{2}+4\left(\frac{qSL{m}_{s}^{\delta }}{{J}_{e}}+i\frac{qSLm″{ }_{m}\stackrel{ ·}{\gamma }}{{J}_{e}}\right)}$
令RRe=${\left(\frac{qSL{m}_{d}^{{\omega }_{b}}}{{J}_{e}}\right)}^{2}$-${\left(\frac{{J}_{x}\stackrel{ ·}{\gamma }}{{J}_{e}}\right)}^{2}\left(1-\frac{1}{{S}_{g}}\right)$,RIm=2$\frac{{J}_{x}\stackrel{ ·}{\gamma }qSL{m}_{d}^{{\omega }_{b}}}{{J}_{e}^{2}}$Sd,R=$\sqrt[ ]{{R}_{Re}^{2}+{R}_{Im}^{2}}$,则式(15)可表示为:
λ1,2=$\frac{1}{2}\left(-\frac{qSL{m}_{d}^{{\omega }_{b}}}{{J}_{e}}+i\frac{{J}_{x}\stackrel{ ·}{\gamma }}{{J}_{e}}\right)$±$\frac{1}{2}\left(\pm \sqrt[ ]{\frac{R+{R}_{Re}}{2}}\pm i\sqrt[ ]{\frac{R-{R}_{Re}}{2}}\right)$
式中Sd为动稳定因子,表达式为:
Sd=2$\frac{{J}_{e}m″{ }_{m}}{{J}_{x}{m}_{d}^{{\omega }_{b}}}$-1
式(16)中,特征根正负号搭配取决于RIm的正负,由RIm的表达式可以看出RIm与Sd同号,下面分情况讨论特征根的表达式:
(1)当动稳定因子Sd>0时,特征根为:
λ1,2=$\frac{1}{2}\left(-\frac{qSL{m}_{d}^{{\omega }_{b}}}{{J}_{e}}\pm \sqrt[ ]{\frac{R+{R}_{Re}}{2}}\right)$+$\frac{1}{2}$i$\left(\frac{{J}_{x}\stackrel{ ·}{\gamma }}{{J}_{e}}\pm \sqrt[ ]{\frac{R-{R}_{Re}}{2}}\right)$
令
$\left\{\begin{array}{l}{\sigma }_{1}=\frac{1}{2}\left(-\frac{qSL{m}_{d}^{{\omega }_{b}}}{{J}_{e}}+\sqrt[ ]{\frac{R+{R}_{Re}}{2}}\right)\\ {\sigma }_{2}=\frac{1}{2}\left(-\frac{qSL{m}_{d}^{{\omega }_{b}}}{{J}_{e}}-\sqrt[ ]{\frac{R+{R}_{Re}}{2}}\right)\end{array}\right.$
$\left\{\begin{array}{l}{\omega }_{1}=\frac{1}{2}\left(\frac{{J}_{x}}{{J}_{e}}\stackrel{ ·}{\gamma }+\sqrt[ ]{\frac{R-{R}_{Re}}{2}}\right)\\ {\omega }_{2}=\frac{1}{2}\left(\frac{{J}_{x}}{{J}_{e}}\stackrel{ ·}{\gamma }-\sqrt[ ]{\frac{R-{R}_{Re}}{2}}\right)\end{array}\right.$
则式(14)的通解可表示为:
Δ=C1eσ1t${e}^{i{\omega }_{1}t}$+C2${e}^{{\sigma }_{2}t}{e}^{i{\omega }_{2}t}$
其中C1=ρ1${e}^{i{\epsilon }_{1}}$,C2=ρ2${e}^{i{\epsilon }_{2}}$。根据欧拉公式eix=cosx+isinx,式(21)可写为:
Δ=[ρ1${e}^{{\sigma }_{1}t}$cos(ω1t+ε1)+ρ2${e}^{{\sigma }_{2}t}$cos(ω2t+ε2)]+i[ρ1${e}^{{\sigma }_{1}t}$sin(ω1t+ε1)+ρ2${e}^{{\sigma }_{2}t}$sin(ω2t+ε2)]
根据复数角的定义可得β*和α*的表达式为:
$\left\{\begin{array}{l}{\beta }^{*}={\rho }_{1}{e}^{{\sigma }_{1}t}cos({\omega }_{1}t+{\epsilon }_{1})+{\rho }_{2}{e}^{{\sigma }_{2}t}cos({\omega }_{2}t+{\epsilon }_{2})\\ {\alpha }^{*}={\rho }_{1}{e}^{{\sigma }_{1}t}sin({\omega }_{1}t+{\epsilon }_{1})+{\rho }_{2}{e}^{{\sigma }_{2}t}sin({\omega }_{2}t+{\epsilon }_{2})\end{array}\right.$
(2)当动稳定因子Sd<0时,特征根为:
λ1,2=$\frac{1}{2}\left(-\frac{qSL{m}_{d}^{{\omega }_{b}}}{{J}_{e}}\pm \sqrt[ ]{\frac{R+{R}_{Re}}{2}}\right)$+$\frac{1}{2}$i$\left(\frac{{J}_{x}\stackrel{ ·}{\gamma }}{{J}_{e}}\mp \sqrt[ ]{\frac{R-{R}_{Re}}{2}}\right)$
此时式(14)的通解为:
$\begin{aligned} \Delta= & C_{1} e^{\sigma_{1} t} e^{\mathrm{i} \omega_{2} t}+C_{2} e^{\sigma_{2} t} e^{\mathrm{i} \omega_{1} t} \\ = & {\left[\rho_{1} e^{\sigma_{1} t} \cos \left(\omega_{2} t+\varepsilon_{1}\right)+\rho_{2} e^{\sigma_{2} t} \cos \left(\omega_{1} t+\varepsilon_{2}\right)\right]+} \\ & i\left[\rho_{1} e^{\sigma_{1} t} \sin \left(\omega_{2} t+\varepsilon_{1}\right)+\rho_{2} e^{\sigma_{2} t} \sin \left(\omega_{1} t+\varepsilon_{2}\right)\right] \end{aligned}$
同样根据复数角定义可得β*和α*的表达式为:
$\left\{\begin{array}{l}{\beta }^{*}={\rho }_{1}{e}^{{\sigma }_{1}t}cos({\omega }_{2}t+{\epsilon }_{1})+{\rho }_{2}{e}^{{\sigma }_{2}t}cos({\omega }_{1}t+{\epsilon }_{2})\\ {\alpha }^{*}={\rho }_{1}{e}^{{\sigma }_{1}t}sin({\omega }_{2}t+{\epsilon }_{1})+{\rho }_{2}{e}^{{\sigma }_{2}t}sin({\omega }_{1}t+{\epsilon }_{2})\end{array}\right.$
由式(18)、(24),系统稳定的条件为-$\frac{qSL{m}_{d}^{{\omega }_{b}}}{{J}_{e}}$+$\sqrt[ ]{\frac{R+{R}_{Re}}{2}}$<0,整理可得锥形运动稳定性的条件为:
${S}_{d}^{2}$<1-$\frac{1}{{S}_{g}}$
3 仿真与分析
根据短周期动力学模型假设,认为弹体速度矢量保持不变,假设速度矢量沿水平方向,弹道倾角θ=0°,弹道偏角ψv=0°,弹体自旋转速为$\stackrel{ ·}{\gamma }$=8r/s,飞行速度为v=270m/s,极转动惯量Jx=0.2,赤道转动惯量Je=50。同时假设初始状态俯仰角ϑ=3°,偏航角ψ=3°,绕Oy4轴角速度ωy4=2 (°)/s,绕Oz4轴角速度ωz4=2 (°)/s。采用系数“冻结”法,令弹体为静稳定,力矩系数导数${m}_{s}^{\delta }$=-0.03,不考虑马格努斯力矩项和阻尼力矩项时,存在初始扰动状态的弹体在陀螺效应和气动静力矩作用下的锥形运动曲线如图2(a)所示,角速度ωy4和ωz4的变化如图2(b)和图2(c)所示。
图2(a)可以看出在陀螺效应和气动静力矩作用下,锥形运动为二圆运动的合成运动,由式(8)计算得此时陀螺稳定因子Sg<0,因此对于静稳定弹体,锥形运动为稳定的二圆运动。当弹体为静不稳定时,即mδs>0,对于低转速、长细比较大的弹体,Jx≪Je,此时Sg为正且0<Sg<1,因此1-1/Sg<0,故此时弹体运动模态的频率式(11)应表示为ω1,2=$\frac{{J}_{x}\stackrel{ ·}{\gamma }}{2{J}_{e}}$(1±i$\sqrt[ ]{\frac{1}{{S}_{g}}-1}$),当t→∞时,式(12)中第二项C2${e}^{i{\omega }_{2}t}$→∞,即二圆运动模态中的慢圆运动将会发散[20],从而引起锥形运动发散。也就是说,对于静不稳定弹体,陀螺稳定[19](Sg>1)是锥形运动稳定的前提。由式(8)的表达式可以看出,对于低转速、静不稳定的弹体,通过匹配极转动惯量与赤道转动惯量来实现陀螺稳定比较困难,通过提高弹体的自旋转速$\stackrel{ ·}{\gamma }$来实现更为可行。图3为提高转速后,静稳定和静不稳定条件下弹体锥形运动曲线,此时锥形运动均为稳定的二圆运动合成运动,静不稳定情况下的角运动幅值大于静稳定情况。
考虑马格努斯和阻尼效应时,保持其它仿真条件不变,对于静稳定弹体,令${m}_{s}^{\delta }$=-0.03,${m}_{d}^{{\omega }_{b}}$=0.5,m″m=0.05,由式(17)计算可得,动稳定因子Sd>0,符合稳定条件式(27),此时锥形运动是稳定的,弹体静稳定时的收敛的锥形运动曲线如图4所示。由稳定条件式(27)可以看出,对于静稳定的弹体,随着自旋转速$\stackrel{ ·}{\gamma }$的提高,稳定范围将会缩小,甚至造成锥形运动的发散。
在锥形运动收敛的前提下,图5对比了不同转速的锥形运动形态,可以看出转速的降低可以减小锥形运动幅值,加速其收敛,有效抑制了角运动。
当弹体静不稳定时,Sg> 0,同样当弹体的极转动惯量Jx远小于赤道转动惯量Je时,转速的高低均很难使锥形运动收敛,尤其对于低速旋转弹体,0<Sg<1,此时${S}_{d}^{2}$>1-$\frac{1}{{S}_{g}}$,则锥形运动必然是动不稳定的,图6为静不稳定情况下发散的锥形运动。对于某些特定情况,通过匹配一定的赤道转动惯量和极转动惯量,同时提高转速则可以实现静不稳定弹体的动稳定。
4 结论
本文分别推导了陀螺效应、气动静力矩作用下的锥形运动形态和同时考虑马格努斯、阻尼效应作用下的锥形运动方程的解析解,并给出了稳定性条件,所得结论如下:
对于存在初始扰动的无控旋转弹,仅在陀螺效应和静力矩的影响下,低速旋转的静稳定弹体的锥形运动是稳定的二圆运动。对于静不稳定的弹体可以通过提高自旋转速来实现动稳定,在转速相同的条件下,静不稳定弹体的角运动幅值要大于静稳定的弹体。
在考虑陀螺效应和静力矩作用的同时加入马格努斯和阻尼效应的影响,此时旋转弹的锥形运动曲线为收敛或发散的圆运动曲线。对于静稳定的弹体,自旋转速越高锥形运动的幅值越大,收敛速度越慢,同时随着转速的提高,锥形运动的稳定边界将会缩小甚至造成发散。对于具有一定转动惯量的静不稳定弹体,通过提高转速可以实现锥形运动的收敛。