simulink 蹦极跳动力学仿真全过程：单自由度时不变连续系统

1 动力学方程的建立
2 微分方程到框图
3 第一次仿真
4 框图与m文件的互动

1 动力学方程的建立

蹦极跳的动力学方程，是按照牛顿运动定律建立的，分析蹦极人的受力，包括重力$mg$、阻力$f$、绳子拉力$F$，以上各力的矢量和等于$m{x}^{\prime \prime }$。
其中，阻力采用一次阻尼和二次阻尼联合模型，即$$f=a_1{x}^{\prime }+a_2|x^{\prime }|x^{\prime }$$绳子的拉力采用分段函数模型，蹦极绳长度为$L$；蹦极绳仅在蹦极人位移距离超过绳长的时候起拉力作用，拉力系数为$k$
$$F=\{\begin{array}{ll}k(x-L)& x>L\\ 0& x\le L\end{array}$$以上函数中，$x$为蹦极人坐标，取向下为正，并取原点为蹦极人起跳位置，即$x(0)=0$；本项目暂时不考虑蹦极人的起跳速度，即认为$x^{\prime }(0)=0$；。
最终建立的蹦极跳的动力学微分方程为$$m{x}^{\prime \prime }=mg-a_1{x}^{\prime }-a_2{x}^{\prime }|x^{\prime }|-F\text{\hspace{1em}(1)}$$注意：建立微分方程一定有边界条件或初值条件，否则仿真是无法进行的
上述方程中，按照目前网上较为流行的参数，取$a_1=1,a_2=1,k=20N/m,m=70kg，L=30，桥高H=80$
这是一个单自由度时不变连续系统微分方程

2 微分方程到框图

本节介绍如何从一个单自由度时不变连续系统微分方程建立相应的simulink框图。

2.1 框图主线的确定

首先确定微分方程的最高阶导数，本项目为2阶，因此需要2个积分模块（Integrator）才能回到$x$,即$x^{\prime \prime }\mathrm{to}{x}^{\prime }$和$x^{\prime }\mathrm{to}x$两个模块（Block）

其次是将微分方程改写为如下形式$$C_n{x}^{(n)}=\sum _{i=0}^n{C}_{i-1}{x}^{(i-1)}\text{\hspace{1em}(2)}$$本项目的方程（1）恰好是这种形式，就不加以改写了，其中$n=4$，因此加入ADD模块，形成一个等号态势，其中1个加法，三个减法，即

2.2 方程（2）右侧各项的添加

逐步添加方程中（2）各项的系数，即增益。
从布局的角度，建议除常数项外，从低导数阶向高导数阶逐个开展

2.2.1 弹性绳项

弹性绳项是导数0阶的，采用Switch模块实现分段函数，首先引入绳长L，做减法后输入Switch模块，输出再增加增益$k=20$

2.2.2 一次阻尼项

一次阻尼项是$x^{\prime }$相关项$a_1{x}^{\prime }$，从$x^{\prime }$线引出，并增加增益Gain模块，即

2.2.3 二次阻尼项

二次阻尼项$a_2{x}^{\prime }|x‘’|$仍与$x^{\prime }$相关，仍从$x^{\prime }$线引出，增加绝对值模块Abs和乘法模块Product

2.2.4 质量相关$mg$

采用增益的方式处理重力，将质量单独作为常数输入

2.3 方程的封闭

增加方程的左端项，实现方程的封闭，上面ADD模块的输出等于$m{x}^{\prime \prime }$，除以$m$后就和主线空悬的$x^{\prime \prime }$对应上了，此处引入除法模块Divide，实现对质量的同一输入，便于后续m文件调佣

2.4 结果输出

输出位移$x$，速度$x^{\prime }$和加速度$x^{\prime \prime }$输出到workspace中
还计算蹦极人到地面距离，即$H-x$，采用Scope展示

3 第一次仿真

3.1 参数与初值

第一次仿真参数按上文执行，$a_1=1,a_2=1,k=20N/m,m=70kg，L=30，桥高H=80$
初值按$x(0)=0，x^{\prime }(0)=0$；
仿真时间300s

3.2 第一次仿真结果

第一次仿真结果如下图所示，可见第一次下落是蹦极人到地面的距离已经是负值，安全性评估不合格。

3.3 第一次仿真结果的评估

这里的评估不是安全性的评估，而是仿真正确性的基本评估。

3.3.1 无阻尼极值的评估

无阻尼的情况下，以0m处为零重力势能面，系统总能量$E=0J$
在$x<30m$阶段，蹦极人为自由落体状态，耗时$t_1=\sqrt{2L/g}=2.47s$
在$x>30m$阶段终末，系统总能量保持不变，重力势能$-mgx$，弹性势能$1/2k(x-L{)}^2$，两者之和为零，即$x^2-128.6x+900=0$,解得$x=121.2$
仿真结果如下，可见$x$极值，及位移速度加速度的周期性表现良好，符合物理预期

3.3.2 平衡位置评估

采用较长的仿真时间，如500s，结果如下，其平衡位置值约为15.6~15.9m

从方程（1）出发，考虑平衡状态时，$x^{\prime \prime }=0,x^{\prime }=0$，则方程（1）变为
$$F=k(x-L)=mg$$解算x后计算$H-x$即
$$H-mg/k-L=80-70\ast 9.8/20-30=15.7m$$
与仿真值吻合良好

4 框图与m文件的互动

本项目不涉及从Matlab的Workspace读取文件，只涉及采用常量模块的变量名输入；
本项目设计从SIMULINK的结果输出到Workspace后读取进行后续操作

4.1 常量模块的变量名输入

常量模块中，调整质量$m$和弹性系数$k$，注意，为了方便调试，在写m
文件中，先定义m和k。注意框图中两个常量模块的变化。

clc;clear

m=70;k=20；

4.2 仿真结果的输出

采用整体结果全部输入到一个结构体simOut的方式输出，便于后续变更，绘图展示结果一致性

clc;clear

m=70;k=20;
simOut = sim('csdn01');
t=simOut.tout;
x=simOut.x.Data;
dx=simOut.dx.Data;
d2x=simOut.d2x.Data;
plot(t,x,t,dx,t,d2x)

SIMULINK结果

M文件结果