简介：面向吸气式高超声速飞行器纵向动态建模与仿真需求，压缩包内提供刚体与弹性体两套Simulink模型及对应绘图脚本，可直接运行并输出可对比的速度、加速度、姿态角等飞行状态曲线。刚体模型从整体运动学与动力学出发，忽略结构变形，计算简单、运行快，适合概念设计阶段的快速摸底；弹性体模型引入结构弹性特性，能反映高超声速气动力与热效应引起的变形影响，更接近真实动态响应。纵向刚体模型主要描述俯仰方向的力和力矩平衡，弹性体模型则计入结构弹性模态对纵向运动的影响，两者联合有助于理解刚弹耦合问题。资源共5个文件，包括2个.slx模型和3个.m脚本，整体仅32KB，结构紧凑，适合航空航天、控制工程方向的学生与工程师学习使用。目前已有389人学习，模型注释清晰、层级简单，便于二次开发和参数调整。通过两模型的对比结果，可清晰评估弹性效应对飞行稳定性和控制律设计的影响，为后续控制器设计与结构优化提供直观依据。
1. 为什么吸气式高超声速飞行器要把刚体和弹性体分开建模
吸气式高超声速飞行器飞行速度一般在5马赫以上，机身细长，配平攻角不大，机体的结构模态频率往往不会很高。尤其是吸气式构型，进气道和尾部喷管向外延伸，前体长度占全机比例大，第一阶弯曲频率可能只有几十赫兹，这与刚体短周期频率只差一个数量级。若纵向仿真模型仍把机体当成无限刚性，控制器带宽稍微拉高，就会在仿真中看不到弹性回路的真实作用，等到实际飞行时控制面激励起机身抖动，才意识到问题。因此在做控制律验证时，通常要同时准备两套模型：一套是纵向刚体仿真模型，用于传统飞控设计；另一套是弹性体仿真模型，用于检查结构模态是否被激励。把它们放在 Simulink 里用同一激励跑一遍，再用 MATLAB 把响应曲线画在一起，可以直观看出哪些频段、哪些操纵量影响弹性模态。这个做法对做吸气式高超声速飞行器总体、飞控和气动弹性耦合的工程师都适用。
2. 纵向刚体模型与弹性体模型的动力学差异
2.1 刚体纵向动力学方程是怎么建立的
吸气式高超声速飞行器纵向刚体模型通常取速度、攻角、俯仰角速率和俯仰角作为核心状态。以巡航状态平飞配平点为参考，使用小扰动线性模型，常见形式是：
$\dot{x}_r = A_r x_r + B_r u$
其中 $x_r = [V, \alpha, q, \theta]^T$，$u = [\delta_e; \phi]$ 包含升降舵偏角和燃油当量比。刚体模型把飞行器看作质量集中刚体，忽略结构变形。对吸气式飞行器来说，推进系统与机体耦合严重，推力线和气动中心之间距离大，所以在 $B_r$ 里往往会加入推力对攻角变化率的扰动项。
以短周期模态为例，攻角和俯仰速率的子系统可以单独写出来：
$\begin{bmatrix} \dot{\alpha} \ \dot{q} \end{bmatrix} = \begin{bmatrix} Z_\alpha & 1 \ M_\alpha & M_q \end{bmatrix} \begin{bmatrix} \alpha \ q \end{bmatrix} + \begin{bmatrix} Z_{\delta_e} \ M_{\delta_e} \end{bmatrix} \delta_e$
这里的 $Z_\alpha$、$M_\alpha$ 等是纵向稳定性导数，来自配平点的气动系数插值。对于吸气式高超声速飞行器，$M_\alpha$ 往往比普通飞机更大，因为气动中心与重心的距离经过精细匹配后很短，静稳定裕度不大；而推进系统引起的力矩变化还会叠加到 $M_\phi$ 上。建模时要把燃油当量比与推力系数之间的延迟也考虑进去，通常用一个一阶惯性环节表示。
在得到 A、B 矩阵之后，用 MATLAB 的ss函数可以快速生成连续状态空间模型，然后检查特征值：
% 刚体模型状态矩阵示例（巡航配平点小扰动） Ar = [ -1.4e-3 -200 0 -9.8; -2.5e-4 -0.02 1 0; 8.0e-5 -3.2 -0.35 0; 0 0 1 0 ]; Br = [ 8.2e-3 0; -1.4e-4 2.1e-4; -7.9e-2 0.8; 0 0 ]; sys_r = ss(Ar, Br, eye(4), zeros(4,2)); eig(sys_r)代码里的第一行是速度方程对速度的阻尼导数，数值由轴向气动阻力与质量比决定；第二行攻角项来自升力线斜率；第三行俯仰力矩项是短周期核心。eye(4)作为输出矩阵是为了把状态直接导出来供 MATLAB 后处理。这里的数值不是某架真实飞行器的定版参数，只用于说明矩阵结构，正式项目中通常从气动文件读表，再用脚本自动装配。
2.2 弹性体模型在哪些地方比刚体模型多东西
弹性体仿真模型在刚体的运动自由度上，再叠加若干结构模态坐标。纵向弹性模型常用模态叠加法，把结构变形表示为广义坐标 $\eta_i$ 乘模态振型。
每个弹性模态都用一个等效二阶系统表示：
$\ddot{\eta}_i + 2\zeta_i \omega_i \dot{\eta}_i + \omega_i^2 \eta_i = N_i$
右边 $N_i$ 是模态广义力。与刚体模型的差异在于，这个广义力不仅包含气动力，还包含刚体运动加速度引起的惯性耦合项。在吸气式高超声速飞行器中，柔性前体和尾部喷管壁面会随攻角变形，改变进气道捕获面积和喷管膨胀比，因此弹性模态对推进系统也有反作用。
在 Simulink 中搭建弹性体模型时，我一般不用有限元程序直接算全机瞬态响应，因为时步太长。常见做法是先把机体的有限元模型做模态分析，提取前几阶纵向弯曲模态的频率、阻尼和振型斜率，然后作为参数放进 Simulink。这样既能反映弹性对传感器信号的贡献，又不必在每个积分步求解大规模矩阵。
2.3 气动弹性耦合项的处理方式
在简化工程模型里，常把刚体状态和弹性状态叠加到一个状态空间中：
$\begin{bmatrix}\dot{x}r\\ddot{\eta}\end{bmatrix} = \begin{bmatrix}A{rr} & A_{re}\A_{er} & A_{ee}\end{bmatrix}\begin{bmatrix}x_r\\dot{\eta}\end{bmatrix} + \begin{bmatrix}B_r\B_e\end{bmatrix}u$
下表是二者在处理耦合时的典型差异：
| 对比项 | 刚体模型 | 弹性体模型 |
|---|---|---|
| 状态量 | 速度、攻角、俯仰率、俯仰角 | 刚体状态加弹性模态位移、速度 |
| 运动方程 | 质量、力与力矩平衡 | 力平衡加模态二阶方程 |
| 弹性模态频率 | 假设无穷大 | 用几十到几百 rad/s 的模态表示 |
| 对传感器位置 | 不关心 | 速率陀螺测到的角速度包含模态振型斜率贡献 |
| 对控制器频率上限 | 可按刚体带宽设计 | 需避免激励低阶模态 |
表达式中的 $A_{re}$ 表示结构变形产生的等效攻角和气动力矩；$A_{er}$ 表示刚体加速度与弹性广义力的耦合；$A_{ee}$ 是一个对角阵，内部是 $-\omega_i^2$ 和 $-2\zeta_i\omega_i$ 的组合。
处理耦合时有一个容易忽略却很重要的点：弹性体模型里的攻角并不是刚体模型的攻角，而是各弹性振型在传感器或气动参考点处的弹性偏转角度叠加出来的“当地攻角”。因此，如果你用刚体模型设计控制律，再把控制律接到弹性体模型上，两个模型内部的状态量定义必须一致。Simulink 模型里推荐用总线对象把这些不同含义的信号打包，减少接线配错。
2.4 弹性模态参数从哪里来
弹性模态的频率和振型依赖结构质量和刚度分布。高超声速飞行器在飞行中受气动加热影响，结构材料刚度会随温度下降，弹性频率也可能降低 5%～15%。所以常温地面试验得到的频率不能直接用于热状态仿真。常见做法是取冷态和热态两组参数，分别跑弹性体仿真模型，对比控制器稳定性边界。
参数数量太多时可以先回归刚体模型，通过落振实验辨识前两阶频率和阻尼比，然后在 MATLAB 里用ssest或简单的最小二乘拟合得到模态参数。需要说明的是，阻尼比通常很小，在 0.01 到 0.04 之间，测量误差对幅频响应峰值影响很大，仿真时要带范围再跑一遍。
3. 在 Simulink 中搭建刚体与弹性体两个纵向仿真模型
3.1 顶层模块的划分
常见做法是建立一个顶层模型，内部按刚体、弹性、气动、推进分成多个 subsystem。刚体模型用连续时间状态方程；弹性体在刚体子系统外面并联模态子系统。这样同一个激励信号可以同时喂给两个模型，方便后续对比。
顶层中，先由输入生成模块产生升降舵和燃油当量比的控制指令，再分成两路进入刚体模型和弹性体模型。为了让两者的输出能在同一个时间基准上比较，两个子系统的积分器必须在同一个求解器任务中运行，不能给其中一个单独设固定步长而另一个用变步长。虽然 Simulink 允许每个子系统有自己的采样时间，但整个模型只有一个主求解器。
我一般会把刚性模型的仿真配置为 ode45，相对误差设为 1e-6；弹性模型则容易因为高频模态出现数值刚性，改用 ode15s 更稳。为了对比公平，两边的求解器步长上限都设成 1e-4 秒。实际上，如果你在一个顶层模型里同时包含两个模型，Simulink 只能选一个求解器，所以更常用的做法是分成两个独立模型，用sim命令分别运行，再把结果取回 MATLAB 统一画图。这样既能在刚体模型上保留 ode45，又能在弹性体模型上用 ode15s。
3.2 刚体模型的最小 Simulink 实现
刚体模型建议用 MATLAB 初始化脚本设置变量，再由 Simulink 模块引用参数。下面的脚本放在模型的InitFcn回调中执行：
% 吸气式高超声速飞行器纵向刚体模型初始化 V0 = 2200; % 巡航速度 m/s alpha0 = 3.2; % 配平攻角 deg q0 = 0; theta0 = alpha0 * pi/180; Ar = [ -1.2e-3 -180 0 -9.81; -2.1e-4 -0.015 1 0; 7.2e-5 -2.9 -0.42 0; 0 0 1 0 ]; Br = [ 8.2e-3 0; -1.3e-4 2.0e-4; -8.3e-2 0.75; 0 0 ]; x0 = [V0; 0; q0; theta0];各参数含义：Ar(1,3)为 0，说明速度方程里俯仰速率不直接产生轴向加速度；Ar(2,3)=1表示攻角与俯仰速率之间的运动学关系；Ar(3,2)是俯仰力矩对攻角的导数，负值代表静稳定。Br(3,1)是升降舵操纵力矩效率，数值过大容易造成操纵过强，仿真中如果看到高频振荡可以先检查这里。
在 Simulink 中，用 4 个积分器、增益矩阵Ar和Br按状态方程搭接即可。具体连线步骤是：把积分器输出组合成列向量，复制一份乘以Ar，把控制输入乘以Br，两者相加后送入积分器输入。矩阵乘法用Matrix Multiply模块，但要注意维度。Simulink 里推荐使用Integrator模块的状态端口输出来避免代数环。
3.3 弹性体模态子系统的搭法
弹性体模型在刚体模型基础上加一个“广义力计算”子系统和几个模态二阶子系统。每个模态用两个积分器搭，输入为广义力，输出为模态位移和速度。频率和阻尼参数用常数模块。
初始化脚本添加：
% 弹性模态参数（前两阶） omega = [60; 125] * 2 * pi; % 模态频率 rad/s zeta = [0.02; 0.015]; % 阻尼比 % 模态广义力耦合系数 Q_alpha = [-0.26; -0.11]; % 攻角对广义力的影响 Q_ele = [1.8; 1.15]; % 升降舵对广义力的影响每个模态子系统的内部结构是：广义力N = Q_alpha*alpha + Q_ele*delta_e + N_coupling，其中N_coupling来自刚体加速度。二阶系统可用两个积分器串联：前一个积分器输出模态速度，后一个输出模态位移，再用omega^2和2*zeta*omega反馈到输入端。给广义力加限幅，防止仿真发散。
耦合项的接入需要从刚体子系统引出加速度信号。加速度信号在刚体模型中可以通过Derivative模块从速度求导得到，但噪声大；更好的做法是在刚体模型内部直接用动力学方程的右端项作为加速度输出，用 Goto 传给弹性子系统。
3.4 关键参数设置与常见坑
下表列出两个模型的关键参数：
| 参数 | 刚体模型 | 弹性体模型 | 影响 |
|---|---|---|---|
| 求解器 | ode45 | ode15s | 弹性模型高频模态易刚性 |
| 最大步长 | 1e-3s | 1e-4s | 防止漏掉模态峰值 |
| 弹性模态阶数 | 0 | 2~5 | 太少丢频率特性 |
| 传感器信号 | 角度、加速度 | 需加附加振型斜率项 | 否则姿态响应偏乐观 |
| 推进延迟 | 一阶环节 | 与模态同入推进系统 | 延迟影响相位 |
发现模型不稳定时先不要加控制律，先看开环阶跃响应。刚体模型的开环响应应该是发散的但频率有规律；弹性体模型则必须看到叠加的高频波纹，看不到说明耦合系数太大或模态阻尼被设置错误。另外，MATLAB 2023b 之后的 Simulink 对矩阵参数编码更严格，parameter对象必须是 double 类型，不要在初始化脚本里用single或int8。
3.5 从工作区加载输入信号做对比驱动
在对比实验中，控制输入信号通常从 MATLAB 工作区读入，比如阶跃、倍频扫频或真实轨迹。使用 From Workspace 模块，将输入定义为结构体：
u.time = (0:0.0001:5)'; u.signals.values = [delta_e_cmd, phi_cmd]; u.signals.dimensions = 2;在模型参数设置里，把输入端口接到 From Workspace。这样修改输入波形时不需要改模型结构，只需重新生成这个变量。注意别把u.time写成非等间距，求解器变步长会插值，但过大的时间间隙会被忽略，导致激励信号台阶化。
4. 用 MATLAB 跑完仿真并画出刚体/弹性体对比曲线
4.1 在 Simulink 里配置输出端口
要对比，先保证两个模型接收到相同控制输入。常见做法是从工作区读入时间序列信号u.time、u.signals.values。用sim命令运行两个不同模型或同一个模型的不同模式。如果是同一个顶层模型，用set_param切换弹性模态开关。
示例：
% 运行刚体模型 set_param('hypersonic_longi', 'SimulationCommand', 'stop'); set_param('hypersonic_longi', 'StopTime', '5'); simOut_r = sim('hypersonic_longi', 'StopTime', '5', ... 'SolverType', 'Variable-step', 'Solver', 'ode45'); % 运行弹性体模型 simOut_e = sim('hypersonic_longi', 'StopTime', '5', ... 'SolverType', 'Variable-step', 'Solver', 'ode15s', ... 'MaxStep', '1e-4');两个simOut中分别包含仿真时间、状态和输出信号。由于两个模型内部状态数量不同，输出端口要固定为同样的维度，一般只输出需要对比的物理量。
4.2 写对比脚本来画响应曲线
% 提取结果 t_r = simOut_r.tout; t_e = simOut_e.tout; alpha_r = simOut_r.yout{1}.Values.Data; % 攻角 rad alpha_e = simOut_e.yout{1}.Values.Data; figure('Color','w'); subplot(2,1,1); plot(t_r, rad2deg(alpha_r), 'b-', 'LineWidth',1.5); hold on; plot(t_e, rad2deg(alpha_e), 'r--', 'LineWidth',1.5); xlabel('时间 (s)'); ylabel('攻角 (deg)'); legend('刚体模型','弹性体模型','Location','best'); title('纵向短周期响应对比'); grid on; axis tight;解释：时间向量长度不同，plot 会自己匹配 x 值。如果 y 维度不对，检查To Workspace是否选了Timeseries。rad2deg是 MATLAB 自带函数，不需要自己乘以 180/pi。
为了得到更完整的对比，还可以增加一个子图绘制俯仰速率：
subplot(2,1,2); plot(t_r, simOut_r.yout{2}.Values.Data, 'b-'); hold on; plot(t_e, simOut_e.yout{2}.Values.Data, 'r--'); xlabel('时间 (s)'); ylabel('俯仰角速率 (deg/s)'); grid on;提示：
To Workspace模块最好选Timeseries输出格式，这样提取数据时不会出现维度被打平或变量名被自动改的情况。
4.3 从对比图里识别弹性耦合特征
图上通常看到三类现象：响应初始段出现刚体模型没有的高频波纹；转折时间点提前或滞后；稳态值分离。第一类说明模态被控制输入激励，第二类说明刚体-弹性耦合改变了等效阻尼，第三类说明推进或弹性形变改变了配平。其中第二类对控制增益设计影响最大。
用findpeaks自动读取峰值并比较超调量：
[pks_r, loc_r] = findpeaks(alpha_r); [pks_e, loc_e] = findpeaks(alpha_e); peak_diff = pks_e - pks_r;如果alpha_e的高频分量太多，直接findpeaks会找到很多毛刺，先做低通滤波或先对弹性响应提取包络。包络可以用envelope函数。
如果要观察频率成分，用fft画频谱：
fs = 1000; N = length(alpha_e); freq = (0:N-1)*fs/N; alpha_fft = abs(fft(alpha_e-mean(alpha_e))); plot(freq(1:N/2), alpha_fft(1:N/2));看频谱峰值是否出现在设定的模态频率周围。适当增加仿真时长，频率分辨率会更好。
4.4 配平点的差别也要一起记录
刚体模型和弹性体模型的配平攻角通常不同。用 Trim Analysis 或者手动仿真到稳态后读取。在 MATLAB 中可以用findop和linearize在指定状态点线性化弹性模型，得出新的短周期频率和阻尼。把这些值列成表，作为控制器增益重新整定的参考。配平差异要单独记录，不要混在动态响应对比图里，因为配平和动态是不同时间尺度的问题。
5. 比较模型时的三个工程技巧
5.1 模态取舍：不要一次把弹性模态全加上
实际飞行器结构模态很多，但纵向控制器带宽通常在 10 rad/s 以下，只保留最低一两阶弹性模态即可。加模态时，先在频率响应里确保模态频率离刚体频率有足够差异，否则把它视为刚体的一部分，数值上容易产生病态矩阵。
在 Simulink 中可以通过 enable 子系统切换模态数量，或者用变量控制模态子系统的通断。先用一阶模态粗略看趋势，再加上第二阶验证耦合。不要在一开始就把五阶模态都加进去，否则调试每一步都要处理大量波形，很难定位问题。
5.2 用滤波器去掉仿真噪声，但不滤掉物理模态
比较两者曲线时，如果弹性模型输出直接用，高频模态可能与求解器步长共振，看起来噪声很大。此时可以在 MATLAB 绘图中用smoothdata或butter低通滤波，但截止频率必须低于第一阶弹性模态频率而高于刚体短周期频率。这个窗口需要根据仿真采样频率调整。
例如第一阶模态频率是 60 rad/s 约 9.5Hz，刚体短周期约 2Hz，截止可取 20Hz：
fs = 1000; fc = 20; % Hz，低于第一阶模态95Hz但要高于短周期 [b,a] = butter(4, fc/(fs/2), 'low'); alpha_e_filtered = filtfilt(b,a,alpha_e);注意，filtfilt是零相位滤波，适合画图对比；但用于控制律反馈时，零相位滤波会引入未来信息，实际实现不建议用。只在离线和对比数据分析时使用。
5.3 验证弹性模态是否正确激发
可以用正弦扫频输入或快速升降舵阶跃。给 30% 杆力阶跃，快速傅里叶变变换频谱中看到弹性频率峰值就说明激发成功。然后用modalfit或手动测量峰值间距与阻尼比，与理论值比较。如果激励不到，检查模态广义力中的耦合系数符号和大小，符号错误会让响应反向叠加，模态被抵消。
具体操作：
% 正弦扫频输入 0.5~40Hz，幅值 5deg f0 = 0.5; f1 = 40; T = 30; t_sweep = (0:1/fs:T)'; u_sweep = 5 * sin(2*pi*(f0*t_sweep + (f1-f0)*t_sweep.^2/(2*T)));把u_sweep写入工作区后从 From Workspace 输入，运行弹性体模型，对攻角响应做短时傅里叶变换spectrogram，横轴时间、纵轴频率，能看到刚体模态和弹性模态的贯穿线。如果弹性模态线不连续，说明激励能量不够或模态阻尼太大。调整扫频幅值时不要超过升降舵偏转角限，否则光看到饱和段。
本文还有配套的精品资源，点击获取
