结构随机细长体飞行器动力学建模与分析

齐 麟，杨 亚，陈意芬，殷 玮，赵宏宇

（1.上海交通大学航空航天学院，上海 200240；2.上海机电工程研究所，上海 201109）

细长体构型在飞行器外形设计中有着广泛的应用，如运载火箭、导弹及飞机机身。细长体外形的自身性质导致此类飞行器在飞行过程中容易产生弹性形变，并且表现出复杂的振动特性。振动模态呈现出低频、密频的特点［1］。这些复杂密集的低频振动模态在控制系统的设计过程中起着重要的作用。由于模态频率降低至控制系统带宽内，从而导致控制系统的输出会受到弹性振动的影响，伺服机构输出的控制力矩又导致气动力重新分布，进一步影响细长体变形。与此同时，随着大量先进新型材料在航空航天领域的广泛应用，在一定程度上增强了飞行器结构参数的随机性（如弹性模量、泊松比等）［2］。结构随机因素使得结构动力学稳定性问题呈现出新的特点［3］。

导致结构失稳的原因有很多［4］。其中一种为推力及气动力相互作用下导致的结构失稳［5］。Beal［6］开展随动推力影响下的稳定性问题研究较早，对均匀梁模型展开了常值及脉冲推力影响的稳定性分析。Wu［7-8］将细长体飞行器视为两端自由的欧拉-伯努利梁，采用有限元方法研究了细长飞行器在矢量控制推力作用下的失稳机理，分析了结构稳定性与弹性模态的关系。Park［9］用均匀、两端自由的Timoshenko梁代替细长体，利用扩展的Hamilton 原理，建立了描述横向平面运动的动力学方程，研究了转动惯量及剪切变形对稳定性的影响。Trikha［10-11］考虑了刚性、纵向和横向振动模态的耦合，建立了分析细长体飞行器在推力和空气动力作用下稳定性的动力学模型。采用时频Fourier 谱-有限元法和H-P 有限元法对动力学模型进行离散，并对其发散和颤振现象进行了分析。宋健［12］推导了推力及空气动力作用下细长体飞行器横向运动方程及边界条件，分析了横向振动模态出现低频特性的原因。许赟等［13］将细长体构型简化为变质量非均匀Euler-Bernoulli梁模型，将推力定义为随动力，依据Hamilton原理，建立了描述结构变形的动力学方程。通过有限元方法对方程进行解算，研究了不同推力作用下结构失稳的特性。荣吉利等［14］将旋转细长体飞行器简化为带有附加质量的Timoshenko 梁模型，同时考虑旋转导致的陀螺力矩及随动推力的影响，建立了结构动力学方程。通过有限元方法进行了仿真分析，研究了旋转速度、剪切变形、随动推力及附加质量等因素对结构稳定性的影响。鄢雄伟等［15］建立了旋转导弹在变推力影响下的耦合动力学模型，研究了推力快速变化对弹体动力稳定性的影响，并给出了推力变化率设计边界。

上述在细长体构型飞行器结构稳定性方面的研究，均建立在结构参数为确定值的基础上。然而，随机性是材料的固有属性，尤其随着先进材料的应用及飞行任务需求的发展，随机性的影响日益增强。因此，有必要对结构随机因素影响下的细长体稳定性问题给予更多的关注。

本文将对结构随机因素及推力作用下的结构稳定性进行研究。首先基于矩阵摄动理论，对结构随机因素导致相应参数矩阵发生改变时，细长体动力学模型进行推导。然后，求解反映特征值变化的显式表达式，并得到其期望与方差，形成一套随机因素在推力作用下细长体飞行器动力学稳定性问题中影响的完整研究方法。同时，为验证方法有效性，对本文所得到的结果进行了数值仿真验证。

1 推力作用下的细长体飞行器动力学建模

随着对细长体飞行器有效载荷及射程需求的增加，设计上往往要求飞行器具有较高的推重比。此时，推力对于细长体弹箭动力学稳定性的影响需要给予考虑。首先，对推力作用下的动力学稳定性分析模型进行推导。飞行器受力如图1所示，其中fp为发动机推力，fadw为气动阻力。

图1 飞行器受力Fig.1 Definition of force of slender vehicle

将弹体横向位移w(x，t)按振动模态进行展开，有

式（1）中，ϑi(x)为振动模态，qi(t)为广义坐标。不考虑结构阻尼，有

式（2）中，及分别为质量矩阵、结构刚度矩阵及推力引起的附加刚度矩阵，具体形式如下：

其中，m(x)为质量密度；E为弹性模量；I为惯性矩；N(x)为轴向力，具体形式如下：

对于动力学稳定性问题，可将广义坐标qi(t)表示为

式（4）中，为特征值。将式（4）带入式（3）可得

式（5）中，

由式（5）可知，A，B矩阵均为非对称矩阵，所以特征值为复数，形式如实部表示振动频率，虚部表示振动阻尼。根据稳定性概念，当实部＜0 时，结构振动会收敛，系统为稳定状态。当实部＞0时，结构振动会发散，则系统不稳定。当实部=0 时，结构振动位移为等幅振荡，系统处于临界状态。此时的推力值即为临界推力。临界推力的计算对于结构稳定来说具有重要的意义，它决定了飞行包线的形状以及任务剖面的设计。

结构随机因素的随机性会使得振动模态ϑi(x)具有一定的随机性。而振动模态的变化会让及发生相应的改变，从而使得A，B矩阵具有相应的随机性，由此计算得到的推力临界值也为随机参量。

结构随机因素的改变会引起动力学稳定边界的变化，若采用传统的分析方式研究稳定边界的变化规律，则结构随机因素每变化一次就要对细长体飞行器的广义特征值问题求解一次，需进行多次迭代计算，分析过程非常耗时和麻烦。因此，本文采用矩阵摄动理论建立一种能够快速分析细长体飞行器动力学稳定边界变化规律的分析方法，以提高计算效率。

下面将基于矩阵摄动理论，对考虑结构随机因素影响时的细长体飞行器动力学稳定性问题进行建模，并求解特征值的解及其期望与方差。

2 结构随机因素影响下的稳定性问题建模

2.1 二阶递推模型建立

对于描述细长体飞行器动力学稳定性的特征值问题，如式（5），与其相对应的伴随特征值问题可表示为

易证特征值问题式（5）及伴随特征值问题式（6）具有相同的特征值。并且模态向量和满足以下条件：

由于弹性模量、泊松比、几何尺寸及质量密度的结构参数具有随机性，因此矩阵A，B为随机矩阵。若采用上一节的方式对A，B矩阵进行泰勒展开，则由于其A，B矩阵的特点，很难得到特征值关于随机因素的解析解。因此，本小节采用矩阵摄动方法研究随机因素对特征值的影响，即研究随机因素对细长体飞行器在推力作用下动力学稳定性的影响。

当结构随机因素产生有限变化时，质量矩阵、结构刚度矩阵及附加矩阵也会有所改变，其变化可表示为

相应的有

式（8）-（9）中，ε表示摄动系数，为小参数；及分别为质量矩阵、结构刚度矩阵及附加矩阵在随机因素均值处的取值。分别表示随机因素引起的变化。当ε→0 时，即为确定性的复特征值问题。

由于上述矩阵发生变化，因此，特征值及模态向量也会发生较小变化。将特征值及模态向量按小参数ε展开成幂级数，有

模态向量满足如下正则化条件：

将式（8）-（10）带入式（5）、（6）及式（11），可建立求解一阶及二阶摄动的递推模型。将递推公式按小参数ε的阶次进行整理。

零阶方程为

一阶方程为

二阶方程为

通过求解一阶摄动方程可以得到特征值及模态向量的一阶摄动值。利用二阶摄动方程可以求出特征值及特征向量的二阶摄动。

2.2 一阶摄动值求解

将模态向量的一阶摄动按展开，即

式（15）中，G0矩阵为待定系数矩阵。将式（15）带入式（13）的第一式有

将左乘式（16），得

由模态向量的正交性可得

式（18）中

同时，待定系数矩阵G0的非对角元素可表示为

与式（15）相类似，将模态向量按展开，即

式（21）中，H0为待定系数矩阵。将式（21）带入到式（13）的第二式可得

左乘上式后得

根据模态向量的正交化性质，上式可整理为

式（24）中

对于待定系数矩阵的对角元素，可由模态向量的正交条件求得。将式（15）及式（21）带入到式（13）的第四式，可得

整理后可得

由于=1，因此上式可改写为

设Gi0i=Hi0i，因此有

通过求解式（19）、（20）、（25）、（29）后可以得到特征值及模态向量的一阶摄动值。当结构随机因素的变异系数较小时，即结构参数发生较小变化，一阶摄动可以达到较理想的精度。但当结构随机因素具有较大的变异系数时，需要对二阶摄动进行求解。

2.3 二阶摄动值求解

对模态向量的二阶摄动按展开，即

将式（30）带入到式（14）中的第一式可得

左乘后得

根据模态向量的正交化条件，式（32）可整理为

式（33）中

从而待定系数矩阵G1的非对角元素可表示为

与模态向量的二阶摄动推导过程相类似，下面推导模态向量的二阶摄动值。将按展开，即

将式（36）带入式（14）的第二式，可得

左乘式（37）后可得

利用模态向量的正交性，式（38）可整理为

式（39）中

由于为对角阵，所以待定系数矩阵H1的非对角元素可表示为

对于矩阵G1及H1的对角元素可通过模态向量的正交特性求得，将式（30）及式（36）带入式（14）中的第四式可得

整理后可得

根据模态向量的正交性，并设Gi1i=Hi1i，可求得矩阵G1及H1的对角元素为

通过求解式（34）、（35）、（40）、（43）可得到特征值及模态向量的二阶摄动。将特征值及模态向量的一阶及二阶摄动值带入式（10）即可得到质量及刚度矩阵发生变化时的特征值及模态向量。

前述推导得到了当矩阵因结构随机因素而产生的一定的摄动时，特征值及模态向量的变化结果。接下来将基于特征值的二阶摄动展开，对特征值的数学期望及方差的显式表达式进行推导。

2.4 特征值的期望与方差

由于结构随机因素的影响，摄动量及均可视为随机矩阵，因此特征值也具有随机性。为研究其统计特性，求取其数学期望与方差，本部分假设矩阵A0、B0、A1、B1中各元素统计特性已知，基于前述得到特征值的二阶摄动展开式，对其期望及方差进行推导。

依据式（10），第i个特征值的数学期望为

由于一阶摄动量及二阶摄动量的期望为0，因此，特征值的期望为

考虑特征值的二阶摄动结果，假设特征值λˉ、一阶摄动量及二阶摄动量相互独立，则方差可表示为

假设A0、B0相互独立，且E(A0)=0、E(B0)=0，则式（47）中展开后的交叉项的期望为0，因此式（48）可写为

将式（49）的平方项展开为标量形式，有

将式（50）-（51）带入式（49）后有

假设矩阵A0、B0中各元素相互独立，则式（52）可整理为

综合式（46）、（53），即可得到特征值的方差。

至此，二阶精度的特征值期望及具有一阶精度的特征值方差的计算模型推导完成。为对本小节推导的当质量矩阵、刚度矩阵变化时的特征值的解以及方差有效性进行验证，下面将对上述推导进行数值仿真分析。

3 数值仿真验证及分析

在本文基础上可以进一步分析推力临界值在结构随机因素作用下的变化。本小节将计算当结构随机因素导致质量矩阵及刚度矩阵发生变化时，特征值λ的解，并与矩阵摄动变化后的理论计算结果进行对比，仿真条件如表1所示。

表1 仿真条件Tab.1 Simulation data of vehicle

为验证本文方法的有效性，当摄动系数分别为0.01、0.02、0.03、0.04及0.05时，进行数值仿真计算，计算结果如表2及表3所示。表2及表3分别为一阶摄动方程及二阶摄动方程的计算结果与理论值的对比。从表中可以看出，本文所建立的方法能够较好地计算质量矩阵及刚度矩阵改变时，特征值的变化。并且二阶摄动得到的结果整体优于一阶摄动。在摄动系数较小时，即ε≤0.03时，一阶摄动所得到的结果也具有一定的精度。表4分别给出了摄动值为0.01 时的理论值、一阶摄动值及二阶摄动值的计算时间，计算条件为CPU：Intel酷睿i5 7200U、2.5 GHz主频，内存：8 G。一阶摄动及二阶摄动计算时间均小于理论值，并且一阶计算时间约为二阶计算的1/2，约为理论值计算的1/4。由于一阶摄动方程的计算量明显小于二阶摄动。因此，在质量矩阵及刚度矩阵变化较小时，可采用一阶摄动方程进行相关计算，从而节省计算时间。但当ε＞0.03 时，一阶摄动的计算结果误差较大，而二阶摄动依然保持较好的误差精度。

表2 一阶摄动计算结果对比Tab.2 The calculation result comparison of the first order perturbation

表3 二阶摄动计算结果对比Tab.3 The calculation result comparison of the second order perturbation

表4 计算时间对比Tab.4 The comparison of computing time

为研究二阶摄动方程解的精度及其所能适应的摄动系数变化范围，下面将进一步改变摄动系数，研究本方法的适用性。

当摄动系数ε从0 变化到0.1 时，本文计算所得的二阶特征值与理论值的比较如图2所示。图2（a）及（b）分别为特征值实部及虚部随变异系数的变化曲线。图2（c）描述的是变异系数变化引起的特征值变化，横坐标为特征值实部，纵坐标为虚部。从图2可以看出，随着摄动系数的增大，本文计算值与理论值的误差逐渐增大。当ε＜0.05 时，即图2中第一个方框中的部分，表明本文计算结果具有较好的准确性。然而当ε＞0.05 时，随着ε的增大，本文计算结果逐渐偏离理论值。这说明本文计算结果的准确性建立在摄动系数较小的前提上，当ε增大到一定值时，本文方法便不再适用。通过图2（a）、（b）可以发现，当不考虑结构阻尼时，质量矩阵与刚度矩阵的变化对特征值虚部的影响较小。

图2 不同ε时，本文计算结果与理论值的对比Fig.2 Comparison between the calculated results and theoretical values for different ε

为了验证本文推力临界值结果的有效性，与参考文献［7］中的方法所得到的结果进行了对比，计算结果如表5及图3。结果表明本文所得的推力临界值具有较好的准确性。至此已对结构参数发生变化时，本文特征值计算结果及推力临界值的正确性进行了仿真验证。然而，推力影响下的细长体飞行器动力学稳定性是最为关注的问题。下面将通过本文方法研究质量矩阵、刚度矩阵的变化对推力临界值的影响。

表5 推力临界值对比Tab.5 Comparison of critical thrust values

图3 特征值随推力的变化Fig.3 Variation of eigenvalue with thrust

根据前述仿真分析，当摄动系数ε的变化范围为0 →0.05 时，本文特征值计算结果较为准确。因此，选择ε的摄动范围为0 →0.05，研究质量矩阵、刚度矩阵发生小幅度改变时，推力临界值的变化。

图4表示的是摄动系数ε的变化对推力临界值的影响。图4（a）为ε从0 递增到0.05 时，特征值实部随推力的变化。可以看出当质量与刚度矩阵发生变化时，推力临界值也会产生相应的改变。摄动系数ε从0递增到0.05，推力临界值从23.75 kN 变化到26.13 kN，变化幅度为10%，并且两者变化呈线性关系，如图5所示。图4（b）为特征值虚部的变化。可以看出，在推力临界值附近，虚部的变化趋势发生转折。对比图4（a）和（b）发现，质量矩阵及刚度矩阵的变化对特征值虚部的影响小于特征值实部。这是因为在计算中忽略了结构阻尼矩阵，而特征值的虚部代表了振动特性的阻尼，因此虚部的变化幅度小于实部。

图4 ε变化对推力临界值的影响Fig.4 Influence of ε on the thrust critical value

图5 ε与推力临界值的关系Fig.5 Relationship between ε and thrust critical value

弹性模量及质量密度等结构参数所具有的随机特性会带来细长体飞行器固有频率及位移模态的变化，这些变化会导致动力学稳定性分析中的质量矩阵及刚度矩阵产生相应的改变，进而影响动力学发散的推力临界值。综合以上分析发现，即使当质量矩阵及发生小幅度改变时，推力临界值也会产生较大幅度的变化（约为前者改变量的两倍）。

当质量矩阵及发生随机小幅度改变时，特征值及模态向量的值也就具有了随机性。本文在假设矩阵A0、B0中各元素统计特性已知的情况下，基于矩阵摄动理论推导了特征值的期望及方差。下面采用Monte-Carlo 方法对本文推导的特征值期望及方差进行验证，样本数量为2000，验证结果见表6。

从表6可以看出，本文方法计算得到的数学期望及方差具有一定的准确性。然而方差的误差整体上比期望的大，这是因为在求导过程中，认为矩阵A0、B0中的各元素相互独立，忽略了其相关性。

表6 特征值期望及方差验证Tab.6 Expectation and variance verification of Eigenvalue

表6 特征值期望及方差验证Tab.6 Expectation and variance verification of Eigenvalue

方差/rad2·s-2推力期望/rad·s-1本文结果MCS 本文结果MCS-0.1201+1.2338i-0.2393+1.2802i-0.3055+1.3218i-0.3526+1.3596i 5 kN 10 kN 15 kN 20 kN-0.1250+1.2289i-0.2412+1.2574i-0.3066+1.3171i-0.3532+1.3549i 1.4785×10-5 9.2746×10-5 8.2458×10-5 8.0021×10-5 1.6109×10-4 9.3801×10-5 8.4660×10-5 8.1478×10-5

4 总 结

本文探究了结构随机因素对推力及气动力作用下的细长体飞行器动力学稳定性问题的影响。基于矩阵摄动理论对结构随机因素导致矩阵及发生相应的改变时，细长体飞行器动力学模型进行了推导，并给出了广义特征值的解及其期望与方差的显式表达式。形成了一套结构随机因素在推力作用下对细长体飞行器动力学稳定性问题中影响的完整研究方法，并对结果有效性进行了仿真验证。同时，当矩阵及发生相应的改变时，对导致动力学失稳的推力临界值的变化进行了仿真分析。结果表明，广义特征值的期望值及方差计算具有较好的精度，并且即使当矩阵及发生小幅度改变时，推力临界值也会产生较大幅度的变化。

相比于传统的蒙特卡洛方法，本文所提出的方法具有不需要大量计算、计算成本低的优势，同时能够对固有频率及模态向量的数字特征进行直接计算。在结构随机因素导致质量矩阵、刚度矩阵发生变化的情况下，不需要反复进行特征值解算即可对临界推力进行分析，提高了计算效率。