4、无人机动力学建模：刚体运动方程、空气动力学简化模型、推力与力矩计算
好，咱们进入正题。无人机动力学建模，说白了就是回答一个问题：这玩意儿到底是怎么飞起来的？
我刚开始做飞控那会儿，总觉得动力学是理论派的事，代码能跑就行。直到有一次，我调了一周的PID参数，飞机就是稳不住——后来才发现，是模型里的转动惯量给错了。嗯，从那以后，我再也不敢轻视建模这一步。
这一节，我会带你从三个层面拆解无人机动力学：
- 刚体运动方程——描述飞机在空中的位置和姿态变化
- 空气动力学简化模型——升力、阻力怎么来的
- 推力与力矩计算——电机和螺旋桨到底给了飞机什么力
咱们先看一张整体框架图，帮你建立全局认知。
4.1 刚体运动方程：位置与姿态的数学描述
无人机在空中，本质上就是一个六自由度的刚体。三个平动自由度（x, y, z），三个转动自由度（roll, pitch, yaw）。
我个人习惯把运动方程拆成两部分：
- 平动动力学：由牛顿第二定律描述，F = m * a
- 转动动力学：由欧拉方程描述，M = I * α + ω × (I * ω)
你想想看，飞机在空中的加速度，其实来自两个部分：一个是发动机推力，一个是重力。阻力嘛，我们后面再聊。
核心公式（机体坐标系下）：
m * dV/dt = F_thrust + F_gravity + F_aero
I * dω/dt = M_thrust + M_gyro - ω × (I * ω)其中：
- V = [u, v, w]^T —— 机体坐标系下的速度
- ω = [p, q, r]^T —— 角速度
- I —— 转动惯量矩阵（3x3）
我在项目中遇到过一个问题：很多人直接把地面坐标系下的加速度方程拿来用，结果仿真出来的轨迹完全不对。为什么？因为推力方向是随姿态变化的，必须做坐标系旋转。
所以，Simulink里建模时，我建议你这样做：
- 用 旋转矩阵 R 把机体坐标系下的推力转换到地面坐标系
- 在地面坐标系下写平动方程
- 用 四元数 或 欧拉角 更新姿态
小技巧： 我个人偏爱用四元数，因为它没有万向锁问题。欧拉角虽然直观，但俯仰角接近±90°时，仿真会崩掉。你可以在Simulink里用“Quaternion”模块，省去手动推导的麻烦。
4.2 空气动力学简化模型：别把简单问题复杂化
说到空气动力学，很多初学者就头大。其实对于多旋翼无人机，我们不需要搞CFD（计算流体力学）。
说白了，我们只需要两个力：升力和阻力。
升力公式：
L = 0.5 * ρ * v² * S * CL阻力公式：
D = 0.5 * ρ * v² * S * CD其中：
- ρ —— 空气密度（标准大气压下约1.225 kg/m³）
- v —— 空速（相对空气的速度）
- S —— 参考面积（机翼面积或桨盘面积）
- CL, CD —— 升力系数和阻力系数
嗯，这里要注意：对于多旋翼，升力主要来自螺旋桨，而不是机翼。所以S通常取桨盘面积，CL和CD需要查桨叶数据或做实验标定。
避坑指南： 我曾经在仿真里直接用常数CL和CD，结果发现飞机在高速飞行时阻力明显偏小。后来才意识到，CL和CD其实是攻角的函数。如果你做的是高速飞行仿真（比如固定翼），一定要用查表法或多项式拟合。多旋翼低速悬停时，用常数问题不大。
在Simulink里，我习惯把气动模型封装成一个子系统：
- 输入：空速、攻角、侧滑角
- 输出：气动力（Fx, Fy, Fz）和气动力矩（Mx, My, Mz）
- 内部：用MATLAB Function或Lookup Table实现CL、CD曲线
4.3 推力与力矩计算：电机和螺旋桨的配合
这部分是实际工程中最容易出问题的地方。你想想看，四个电机转起来，到底给了飞机多大的推力和力矩？
先看单个电机-螺旋桨模型：
推力: T = kT * ω²
反扭矩: Q = kQ * ω²其中：
- ω —— 电机转速（rad/s）
- kT —— 推力系数
- kQ —— 扭矩系数
这两个系数通常由桨叶的几何参数决定。我在项目中一般通过静拉力测试来标定，简单说就是：把电机固定在测力台上，测不同油门下的推力和扭矩，然后拟合出kT和kQ。
接下来是力矩分配。对于四旋翼（X型布局），总推力和三个轴的力矩可以写成：
| 控制量 | 物理意义 | 与电机转速的关系 |
|---|---|---|
| F_total | 总推力（z轴） | kT * (ω₁² + ω₂² + ω₃² + ω₄²) |
| M_roll | 滚转力矩（x轴） | L * kT * (ω₂² - ω₄²) |
| M_pitch | 俯仰力矩（y轴） | L * kT * (ω₁² - ω₃²) |
| M_yaw | 偏航力矩（z轴） | kQ * (ω₁² - ω₂² + ω₃² - ω₄²) |
其中L是电机到重心的距离。你看，偏航力矩靠的是反扭矩差，而不是推力差——这一点很多新手会搞混。
Simulink实现建议：
我一般会建一个“Mixer”模块，输入是期望的[F_total, M_roll, M_pitch, M_yaw]，输出是四个电机的转速指令。这个模块本质上就是一个4x4矩阵的逆运算。
% 伪代码示例
% 分配矩阵 (X型四旋翼)
mix_matrix = [1, 1, 1, 1;
0, 1, 0, -1;
-1, 0, 1, 0;
1, -1, 1, -1];
% 注意：这里需要根据实际布局调整符号个人经验： 电机响应是有延迟的，一般一阶惯性环节就能模拟。时间常数τ大概在0.01~0.05秒之间，具体看电机和电调。别忘了在Simulink里加一个“Transport Delay”或“Transfer Fcn”模块，否则仿真出来的响应会过于理想。
4.4 把模型组装到Simulink里
好了，三个部分都讲完了。怎么把它们拼起来？
我建议的模块结构是这样的：
- 输入层：四个电机的PWM指令（0~1归一化）
- 电机模型：PWM → 转速（含一阶延迟）
- 推力/力矩计算：转速 → 推力和力矩（用上面的公式）
- 刚体动力学：推力和力矩 → 加速度、角加速度 → 积分得到速度和位置
- 气动模型：速度 → 气动力/力矩（反馈到刚体动力学）
- 输出层：位置、姿态、速度、角速度
你想想看，这个结构是不是很清晰？每个模块都可以单独测试，出了问题也好定位。
注意： 仿真步长不要太大。我一般用固定步长1ms，四阶龙格-库塔法。步长太大，刚体动力学积分会发散；步长太小，仿真跑不动。1ms是个不错的折中。
最后，送你一句话：模型是仿真的灵魂，参数是模型的生命。 花时间标定好kT、kQ、I这些参数，比调PID参数重要得多。
公众号:蓝海资料掘金营，微信deep3321