固体火箭发动机随动推力作用下细长体亚临界试车及动力学失稳临界推力预示

史晓鸣，侯凯宇，李海东，夏 鹏，刘陆广，高 阳，王志新，强科杰，赵志茹; SHI Xiao-ming; HOU Kai-yu; LI Hai-dong; XIA Peng; LIU Lu-guang; GAO Yang; WANG Zhi-xin; QIANG Ke-jie; ZHAO Zhi-ru

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

固体火箭发动机随动推力作用下细长体亚临界试车及动力学失稳临界推力预示 PDF

- ORCID：
史晓鸣 ¹
✉
- ORCID：
侯凯宇 ²
✉
- ORCID：
李海东 ¹
- ORCID：
夏鹏 ¹
- ORCID：
刘陆广 ¹
- ORCID：
高阳 ¹
- ORCID：
王志新 ³
- ORCID：
强科杰 ³
- ORCID：
赵志茹 ³

1. 上海机电工程研究所, 上海 201109； 2. 上海航天技术研究院, 上海 201109； 3. 上海航天动力技术研究所,上海 201109

中图分类号： V435； V414.3⁺4

最近更新：2023-09-04

DOI：10.16385/j.cnki.issn.1004-4523.2023.04.001

摘要

设计了随动推力作用下细长体模拟件试车试验系统，开展固体火箭发动机试车；理论推导了耦合模态频率差‑推力关系的多项式，结合亚临界试车数据外推预示失稳临界推力。研究结果表明：理论分析与地面试验结果吻合，只需开展少数几次安全可控的亚临界试车试验，获取系统的振动响应数据，即可准确预测临界推力；亚临界试验验证了随动推力横向分量引起模态间的刚度耦合，对于弯曲振动明显的细长体飞行器，推力对结构的影响应视为随动载荷。

关键词

固体火箭发动机; 随动推力; 细长体; 亚临界试车; 稳定边界预示

1 概述

固体火箭发动机被广泛用作飞行器的动力装置。随着飞行器高加速、远射程等技战指标的不断提升，飞行器广泛运用轻质材料并采用大型薄壁结构设计，气动布局一般设计为细长体升力体布局^［

1］，细长体结构横向弯曲特征愈发显著，此时固体火箭发动机推力作用方向将随着细长体弯曲振动而不断改变，形成随动力；另外，变推力、双脉冲等新型能量管理技术^{［参考文献 2

百度学术}2］在固体火箭发动机上的应用，使得飞行过程中发动机的输出推力幅值变化较大。上述两点使得随动推力作用下细长体振动及其气动弹性行为更为突出及复杂，国内外学者对此开展了大量的研究。Fazelzadeh等^{［参考文献 3

百度学术}3］研究了预扭角度和随动力加载角度对分布式随动力加载下预扭曲悬臂梁的动力学稳定性的影响。Park等^{［参考文献 4

百度学术}4］研究了随动力作用下两端自由圆柱壳的动力学稳定性问题，对于常数推力情况，圆柱壳与梁的临界推力保持一致；但对于脉动推力情况，圆柱壳高阶周向模态会使得耦合失稳形态发生改变。Pradhan等^{［参考文献 5

百度学术}5］研究了导弹弹体结构在受控随动推力作用下的动力学稳定性，并分析了弹上设备集中质量对稳定性的影响。Wu等^{［参考文献 6

百度学术}6］研究了推力对细长体导弹气动弹性的影响，指出推力的增大会导致颤振临界速度的下降。Ahmadian等^{［参考文献 7

百度学术}7］研究了非线性连接梁在分布式随动载荷作用下的稳定性。荣吉利等^{［参考文献 8‑9}8‑9］研究了随动推力作用下柔性旋转飞行器的振动稳定性，发现转速增加会诱发非均匀转子系统的弯曲模态与刚体模态间的耦合，降低系统的临界推力，致使系统动态失稳。全景阁等^{［参考文献 10

百度学术}10］针对轴向推力过载作用下大长细比导弹的研究发现，其刚体模态和弹性模态将出现耦合失稳，且随着轴向过载的增加，弹体的结构刚度特性下降，失稳速度降低。许赟等^{［参考文献 11

百度学术}11］研究发现随动推力会诱发弹箭飞行器的动力学失稳，推力同时影响结构横向振动的频率和振型特性，在推力作用下模态振型波谷发生移动，振幅也会发生改变。

值得指出的是，关于随动力的研究也一度存在争议^［

12］，是否将推力视为随动载荷会使结构的振动特性和动力学稳定性分析得出不同的结论^{［参考文献 1

百度学术}1，13］。因此，通过试验验证及研究工程应用中是否考虑随动力及其影响是十分必要的。目前常规的固体火箭发动机地面试车试验中将固体火箭发动机固定安装于试车台上，如图1所示。6个自由度完全被约束，此时的推力只能沿固定轴向作用于发动机尾端，不能模拟推力作为随动载荷的工况。而且受限于固体火箭发动机试车存在很大的危险性，具备固体火箭发动机试车资质的研究机构及学者较少。公开发表的文献中显示，只有Sugiyama等^{［参考文献 14‑15}14‑15］开展了固体火箭发动机推力作用于悬臂梁自由端切线的相关试验，通过调整发动机推力幅值大小以及悬臂梁的长度、宽度等尺寸，进行了多次试验，根据每次试验中悬臂梁振幅发散与否判断是否到达临界推力值。但多次试验导致的固体火箭发动机设计、工艺、加工、装配、试车的经济成本以及试验场地安全防护、风险管控代价均很高昂；而且在超临界推力的试验工况下，细长体试验对象进入发散的振动状态，一旦细长体结构折断，失去约束的固体火箭发动机高速飞出是极为危险的，试车台及周边试验设备也将受损。如何设计较为安全可控的亚临界试验，通过较少次数的试验，获得有效的系统亚临界振动响应数据，并从中外推预示出系统的失稳临界推力，是当前固体火箭发动机随动推力作用下细长体地面试车试验研究中亟需解决的问题。

图1 固体火箭发动机地面试车试验常规固定方式^［

16］

Fig. 1 Routine fixed method of ground trial test of solid rocket motor^［

16］

本文为解决上述问题，提出了固体火箭发动机随动推力作用下细长体亚临界试车试验及动力学失稳临界推力预示方法。设计了随动推力作用下细长体模拟件试车试验系统，通过自由端初始位移释放使得细长体悬臂梁产生自由衰减振动，然后固体火箭发动机点火施加随动推力，激光多普勒测振仪测试细长体悬臂梁振动速度信号，推力架测力装置记录推力；根据随动推力作用下细长体振动的数学模型，从理论上推导了系统耦合模态频率差‑推力关系的多项式；从几次不同亚临界推力试验下的振动速度信号中分析系统振动模态频率随推力的变化规律，用耦合模态频率差‑推力关系的多项式外推临界推力。试验及预示结果表明，本文方法只需要开展若干次亚临界试验，即可通过对亚临界状态下系统振动响应数据的处理，有效准确地预测出随动推力作用下细长体动力学失稳临界推力。

2 试车系统及试验方法

试车系统如图2所示，前挡承力墙上固定安装推力架测力装置，并与长×宽×高为1200 mm×9 mm×30 mm的铝合金（2A12）细长体模拟件一端连接，细长体模拟件另一端通过螺纹转接与试验固体火箭发动机相连接，形成悬臂结构。同时，在细长体模拟件两侧及上方分别安装限位防护立柱及横梁，防止试验进入发散状态下细长体因摆动幅度过大而折断。

图2 试车系统示意图

Fig. 2 Schematic diagram of trial system

试验固体火箭发动机质量标称值为20 kg，为防止细长体在固体火箭发动机重力作用下的弯曲变形造成自由端下垂，在试车台体底座上安装一支撑装置，同时在固体火箭发动机外壳上箍一聚醚醚酮（PEEK）滚珠与支撑装置顶面点接触，由此将固体火箭发动机撑起以保持细长杆的水平。进一步以橡皮绳向上拉伸吊起固体火箭发动机至滚珠将将刚好脱离支撑装置顶面，以最大限度减小滚珠与支撑装置顶面的摩擦。由此搭建一轴向自由度约束，水平横向及竖直侧向自由度部分释放的试验系统。

在水平面内将细长体自由端及固体火箭发动机拨离平衡位置形成细长体初始弯曲变形，并用聚乙烯线固定，如图3所示。在聚乙烯线上捆绑一小包黑火药用于固体火箭发动机点火前烧断聚乙烯线，从而使得细长体在水平面内形成初始位移扰动下的自由衰减振动。在自由衰减振动过程中，固体火箭发动机点火，其输出推力方向沿着尾端切线方向，由于细长体的弯曲振动导致尾端切线方向也随之交变，从而形成随动推力。

图3 试验件初始状态

Fig. 3 Initial condition of test piece

试验过程中以4台Polytec OFV‑505/5000激光多普勒测振仪测试细长体上的水平横向振动速度响应，ZonicBook 618E振动数据采集器以128 Hz的采样频率采集数字信号；压强传感器及推力架测力装置记录固体火箭发动机燃烧室内压强及推力值；Phantom Miro M320s高速相机以500 fps帧率高速摄像。

3 失稳临界推力外推预示方法

自由端作用随动推力的细长体悬臂结构如图4所示，推力方向始终沿着细长体自由端轴线的切线方向，其振动方程^［

13］为：

m \ddot{w} + c \dot{w} + [k - T (k_{c} + k_{n c})] w = 0

（1）

式中 $m$ ， $c$ 和 $k$ 分别为结构质量、阻尼和刚度矩阵； $w = {\{w_{1} w_{2} \dots w_{e}\}}^{T}$ 为弯曲变形横向位移，下标e表示推力作用的细长体尾端节点； $\dot{w}$ 和 $\ddot{w}$ 分别为弯曲变形横向速度和加速度； $T$ 为推力； $k_{c}$ 为推力沿着细长体结构轴向分量 $T \cdot \sqrt[]{1 - {(\frac{\partial w_{e}}{\partial x})}^{2}}$ 引起的保守力附加刚度矩阵，对于弹性小变形的情况， $T \cdot \sqrt[]{1 - {(\frac{\partial w_{e}}{\partial x})}^{2}} \approx T$ ； $k_{n c}$ 为推力沿着细长体结构横向分量 $T \cdot \frac{\partial w_{e}}{\partial x}$ 引起的非保守力附加刚度矩阵。

图4 自由端作用随动推力的细长悬臂结构示意图

Fig. 4 Structural schematic diagram of slender cantilever beam with follower thrust at free end

需要指出的是，地面试验与飞行状态仍存在天地不一致的情况：地面试验中细长体边界条件为细长体头部固支、尾端自由，而飞行过程中为两端自由；地面试验中，细长体任意截面上的轴向力大小均与推力大小相同，而飞行过程中，细长体轴向力分布需与细长体轴向惯性力分布平衡；地面试验中无气动载荷，而飞行过程中，受气动载荷的影响，会出现气动‑结构‑推力三者耦合的动力学失稳现象。本文主要针对结构‑推力耦合进行了较为深入的地面试验验证方法研究。当本文研究方法及模型需用于飞行状态的气动‑结构‑推力耦合动力学稳定性分析时，可进一步修改边界条件，并考虑轴向惯性力、气动力等载荷的影响。

取无推力状态下细长体悬臂结构质量归一化振型 $Φ$ ，对弯曲变形横向位移 $w$ 分离变量：

w (x, t) = Φ (x) q (t)

（2）

将式（2）代入式（1），并左乘 $Φ^{T}$ 得：

\ddot{q} + C \dot{q} + [Λ_{0} - T (K_{c} + K_{n c})] q = 0

（3）

式中 $q$ ， $\dot{q}$ 和 $\ddot{q}$ 分别为广义位移、速度和加速度； $C$ 为广义阻尼矩阵； $Λ_{0} = d i a g (ω_{0_1}^{2}, ω_{0_2}^{2}, \dots, ω_{0_n}^{2})$ 为无推力状态下细长体悬臂结构1~n阶固有圆频率平方组成的对角阵； $K_{c}$ 为推力沿着细长体结构轴向分量引起的保守力附加广义刚度矩阵， $K_{c}$ 为对角阵，其对角线上元素 $K_{c} (i, i)$ 表征推力对第i阶横向弯曲模态刚度的削弱程度； $K_{n c}$ 为推力沿着细长体结构横向分量引起的非保守力附加广义刚度矩阵， $K_{n c}$ 为非对角阵，其非对角线元素 $K_{n c} (i, j), i \neq j$ 表征推力引起的第i，j阶模态间的刚度耦合。

忽略系统阻尼，式（3）对应的无阻尼系统振动方程为：

\ddot{q} + [Λ_{0} - T (K_{c} + K_{n c})] q = 0

（4）

求解其广义特征值，可得随动推力作用下细长体悬臂结构的1~n阶固有圆频率 $ω_{1}, ω_{2}, \dots, ω_{n}$ 。

对本文试验结构开展数值仿真，系统圆频率与随动推力关系曲线如图5所示，根轨迹如图6所示。在随动推力作用下，随着推力增大，系统前两阶频率将逐渐接近（见图5），当推力增大到1620 N时，前两阶频率差将率先降为0，此时根轨迹图上第2阶模态的根轨迹进入复平面的右半边（见图6），表明发散形态为第2阶模态分支发散。

图5 系统圆频率与随动推力关系

Fig. 5 Relationship between system circle frequencies and follower thrust

图6 根轨迹图

Fig. 6 Diagram of root locus

根据上述规律，仅保留前两阶模态，取 $q = {\{q_{1} q_{2}\}}^{T}$ 进行模态截断，并将 ${\ddot{q}}_{1} = - ω^{2} q_{1}$ ， ${\ddot{q}}_{2} = - ω^{2} q_{2}$ 代入式（4）整理得：

\{\begin{array}{l} - ω^{2} q_{1} + (ω_{0_1}^{2} - T a_{1}) q_{1} - T b_{1} q_{2} = 0 \\ - ω^{2} q_{2} + (ω_{0_2}^{2} - T a_{2}) q_{2} - T b_{2} q_{1} = 0 \end{array}

（5）

式中 $ω$ 为频率； $a_{1}$ 和 $a_{2}$ 为 $K_{c}$ 对角线上元素， $a_{1} = K_{c} (1,1)$ ， $a_{2} = K_{c} (2,2)$ ； $b_{1}$ 和 $b_{2}$ 为 $K_{n c}$ 非对角线上元素， $b_{1} = K_{n c} (1,2)$ ， $b_{2} = K_{n c} (2,1)$ 。

式（5）改写成矩阵形式：

[\begin{matrix} ω^{2} - ω_{0_1}^{2} + T a_{1} & T b_{1} \\ T b_{2} & ω^{2} - ω_{0_2}^{2} + T a_{2} \end{matrix}] \{\begin{matrix} q_{1} \\ q_{2} \end{matrix}\} = 0

（6）

式（6）成立需满足：

|\begin{matrix} ω^{2} - ω_{0_1}^{2} + T a_{1} & T b_{1} \\ T b_{2} & ω^{2} - ω_{0_2}^{2} + T a_{2} \end{matrix}| = 0

（7）

从而得频率方程为：

ω^{4} - (ω_{0_1}^{2} + ω_{0_2}^{2} - T a_{1} - T a_{2}) ω^{2} +

(ω_{0_1}^{2} - T a_{1}) (ω_{0_2}^{2} - T a_{2}) - T^{2} b_{1} b_{2} = 0

（8）

根据二次代数方程根与系数关系可得前两阶圆频率平方差的平方为：

{(Δ ω_{}^{2})}^{2} = {(ω_{2}^{2} - ω_{1}^{2})}^{2} =

{[(ω_{0_1}^{2} - ω_{0_2}^{2}) - T (a_{1} - a_{2})]}^{2} + 4 T^{2} b_{1} b_{2}

（9）

将式（9）化简，系统前两阶圆频率平方差的平方与推力可表示为二次多项式的关系：

{(Δ ω_{}^{2})}^{2} = A_{0} T^{2} + A_{1} T + A_{2},

A_{0} = {(a_{1} - a_{2})}^{2} + 4 b_{1} b_{2},

A_{1} = - 2 (a_{1} - a_{2}) (ω_{0_1}^{2} - ω_{0_2}^{2}),

A_{2} = {(ω_{0_1}^{2} - ω_{0_2}^{2})}^{2}

（10）

因此，至少需要3组数据才可以确定A₀， A₁， A₂，然后外推 ${(Δ ω_{}^{2})}^{2} = 0$ 时刻的临界推力T_cr。

当T=0时，随动推力作用下系统前两阶圆频率与无推力状态下系统前两阶固有圆频率相等， $ω_{1}^{2} = ω_{0_1}^{2}$ ， $ω_{2}^{2} = ω_{0_2}^{2}$ 。进一步开展至少两次不同亚临界推力下的试车试验即可获取另外两组数据，从而求得A₀和A₁。

4 试验验证及数据分析

试验工况及系统振动频率试验和理论分析结果如表1所示。序号0是初始位移释放后发动机不点火状态下的试验，用于获取无推力作用下系统的振动信号，分析系统固有频率；序号1~4为初始位移释放后发动机点火状态下的试验。

表1 试验工况及系统振动频率

Tab. 1 Trial states and system vibration frequencies

序号	标称推力	试验		理论
序号	标称推力	ω₁/(rad∙s^-1)	ω₂/(rad∙s^-1)	ω₁/(rad∙s^-1)	ω₂/(rad∙s^-1)
0	0	3.27	49.0	3.47	52.1
1	1200	7.53	32.7	7.64	33.5
2	1300	7.85	29.2	8.47	31.2
3	1500	10.4	25.1	11.2	25.2
4	1700	-	-	16.3	16.3

各序号点火试验中，固体火箭发动机不同大小的推力通过改变喷管喉径来调节实现。序号1~3的标称推力值分别设计为1200，1300和1500 N，小于前文中临界推力的数值仿真结果1620 N，属于亚临界推力试车试验，用于获取亚临界推力状态下的振动数据；序号4的标称推力值设计为1700 N，属于超临界推力试车试验，用于对比验证根据亚临界试验数据外推获得临界推力结果的准确性。序号1~4的固体火箭发动机试车推力实测值如图7所示。

图7 试验固体火箭发动机推力测试值

Fig. 7 Thrust trial value of test solid rocket motor

典型的亚临界试验状态下系统振动响应如图8所示。亚临界状态下振动过程分为3个阶段，第1阶段为初始位移释放后系统自由衰减振动，如图9（a）~（c）所示；第2阶段为固体火箭发动机点火后，系统在随动推力作用下的振动，如图9（d）~（e）所示；第3阶段为固体推进剂燃尽，火箭发动机熄火后系统自由衰减振动，如图9（f）~（g）所示。对序号0振动响应数据、序号1~3亚临界试验第2阶段的振动响应数据分别做频谱分析，如图10所示，可得随动推力作用下系统的前两阶频率（见表1）。从表1及图10中可见，随着推力逐渐增大，系统1阶频率逐渐上升，2阶频率逐渐下降，两者呈现相互靠近的趋势，试验验证了随着细长体弯曲振动，推力作为随动载荷始终作用于细长体自由端切线方向，其横向分量引起了模态间的刚度耦合。这与一般的压杆稳定问题不同，在压杆稳定问题中载荷作用方向始终沿着未变形前细长体轴线方向，各阶模态刚度及频率只会依次单调下降至零，并不会出现系统频率相互接近、模态耦合的现象。

图8 亚临界试验振动响应（序号1）

Fig. 8 Vibration response of subcritical trial （Serial number 1）

图9 亚临界试验高速摄影

Fig. 9 High-speed photography of subcritical trial

图10 序号0~3系统振动响应频谱（幅值归一化）

Fig. 10 System vibration response frequency spectrum of serial number 0~3 （Amplitude normalization）

根据表1中序号0无推力以及序号1~3亚临界推力状态下的系统耦合模态频率差，以最小二乘法拟和求解式（10）中的系数，得A₀=1.01337，A₁=-5160.06，A₂=5714120。进一步外推得到（Δω²）²=0时刻的临界推力T_cr=1628.7 N，如图11所示。

图11 亚临界试验结果外推预示临界推力

Fig. 11 Extrapolated prediction of critical thrust via subcritical trial results

在此基础上，为验证上述外推临界推力的准确性，进一步开展了序号4标称推力值1700 N，大于临界推力试验外推值1628.7 N的超临界状态试验，其系统振动响应如图12所示。超临界状态下振动过程分为两个阶段，第1阶段为初始位移释放后系统自由衰减振动，如图13（a）所示；第2阶段为固体火箭发动机点火后，系统在随动推力作用下，振幅急剧扩大，动态发散，如图13（b）~（c）所示，证实了通过亚临界试验数据外推临界推力方法的准确性。由于数据采集系统的量程设置为±500 mm/s，故图12测试数据中可见的最大振动速度幅值只有500 mm/s，但实际速度已远大于该值，且还存在细长体与两侧防护立柱的碰撞限位，振动响应数据已被削峰限幅，导致未能有效分析超临界试验第2阶段发散状态下的振动频率。发散过程中细长体的发散变形为2阶弯曲振型形态（见图13（c）），也验证了图6根轨迹中第2阶模态分支发散的仿真结果。振幅过大后，细长体与两侧预先安装的限位防护立柱碰撞（见图13（c）），其位移受限，细长体卡在防护立柱及横梁上出现弯曲塑性变形，如图13（d）所示。

图12 超临界试验振动响应（序号4）

Fig.12 Vibration response of supercritical trial （Serial number 4）

图13 超临界试验高速摄影

Fig.13 High-speed photography of supercritical trial

5 结论

（1）本文提供了一整套随动推力作用下细长体模拟件的亚临界试车试验及临界推力外推预示方法；从理论上推导了系统耦合模态频率差‑推力关系的多项式；结合少数几次不同亚临界推力试验下系统的振动响应数据预测临界推力。该方法经试验验证准确、可行，为工程设计提供了一种可信、有效的验证手段。

（2）通过亚临界试验验证了随动推力的横向分量引起模态间的刚度耦合。因此，在工程设计中，对于弹性弯曲振动明显的细长体飞行器，分析其推力对结构的影响应将其作为随动载荷来考虑。

参考文献

孟光，周徐斌，苗军. 航天重大工程中的力学问题［J］. 力学进展， 2016， 46： 267-322. [百度学术]

Meng G， Zhou X B， Miao J. Mechanical problems in momentous projects of aerospace engineering［J］. Advances in Mechanics， 2016， 46： 267-322. [百度学术]

侯晓，付鹏，武渊. 固体火箭发动机能量管理技术及其新进展［J］. 固体火箭技术， 2017， 40（1）： 1-6. [百度学术]

Hou X， Fu P， Wu Y. Energy management technology of SRM and its development［J］. Journal of Solid Rocket Technology， 2017， 40（1）： 1-6. [百度学术]

Fazelzadeh S A， Karimi-Nobandegani A， Mardanpour P. Dynamic stability of pretwisted cantilever beams subjected to distributed follower force［J］. AIAA Journal， 2017， 55（3）： 955-964. [百度学术]

Park S H， Kim J H. Dynamic stability of a free-free cylindrical shell under a follower force［J］. AIAA Journal， 2000， 38（6）： 1070-1077. [百度学术]

Pradhan S， Datta P K. Dynamic instability characteristics of a free-free missile structure under a controlled follower force［J］. Aircraft Engineering and Aerospace Technology， 2006， 78（6）： 509-514. [百度学术]

Wu L， Xie C C， Yang C. Aeroelastic stability of a slender missile with constant thrust［J］. Procedia Engineering， 2012， 31： 128-135. [百度学术]

Ahmadian H， Azizi H. Stability analysis of a nonlinear jointed beam under distributed follower force［J］. Journal of Vibration and Control， 2011， 17（1）： 27-38. [百度学术]

荣吉利，徐天富，王玺，等. 随动推力作用下柔性旋转飞行器稳定性分析［J］. 宇航学报， 2015， 36（1）： 18-24. [百度学术]

Rong J L， Xu T F， Wang X， et al. Dynamic stability analysis of flexible spinning flight vehicles under follower thrust［J］. Journal of Astronautics， 2015， 36（1）： 18-24. [百度学术]

荣吉利，徐天富，王玺，等. 随动推力作用下柔性自旋飞行器横向振动响应及失稳分析［J］. 兵工学报， 2015， 36（3）： 516-522. [百度学术]

Rong J L， Xu T F， Wang X， et al. Analysis of transverse vibration response and instability of flexible spinning flight vehicles under follower thrust［J］. Acta Armamentarii， 2015， 36（3）： 516-522. [百度学术]

全景阁，叶正寅，张伟伟. 轴向载荷对大长细比导弹稳定性的影响研究［J］. 兵工学报， 2015， 36（1）： 94-102. [百度学术]

Quan J G， Ye Z Y， Zhang W W. Analysis on stability of a slender missile under axial loads［J］. Acta Armamentarii， 2015， 36（1）： 94-102. [百度学术]

许赟，谢长川，杨超. 推力作用下细长弹箭横向振动及稳定性分析［J］. 工程力学， 2009， 26（12）： 211-215. [百度学术]

Xu Y， Xie C C， Yang C. Transverse vibration and dynamic stability analysis of slender projects under thrust［J］. Engineering Mechanics， 2009， 26（12）： 211-215. [百度学术]

Elishakoff I. Controversy associated with the so-called “follower forces”： critical overview［J］. Applied Mechanics Reviews， 2005， 58（2）： 117-142. [百度学术]

Datta P K， Biswas S. Aeroelastic behaviour of aerospace structural elements with follower force： a review［J］. International Journal of Aeronautical and Space Science， 2011， 12（2）： 134-148. [百度学术]

Sugiyama Y， Langthjem M A， Iwama T，et al. Shape optimization of cantilevered columns subjected to a rocket-based follower force and its experimental verification［J］. Structural and Multidisciplinary Optimization， 2012， 46（6）： 829-838. [百度学术]

Sugiyama Y， Katayama K， Kiriyama K， et al. Experimental verification of dynamics stability of vertical cantilevered columns subjected to a sub-tangential force［J］. Journal of Sound and Vibration， 2000， 236（2）： 193-207. [百度学术]

Stadler L J， Huber J， Friedemann D， et al. The double pulse motor demonstrator MSA［C］. 40th International Conference on Environmental Systems，AIAA.Barcelona， Spain， 2010： 6755. [百度学术]