刚体动力学方程：牛顿-欧拉方程
做飞控仿真，绕不开的就是刚体动力学。说白了，就是研究飞行器这个「铁疙瘩」在力和力矩作用下怎么动。我个人习惯把这个问题拆成两半来看：平动和转动。牛顿管平动，欧拉管转动，合起来就是牛顿-欧拉方程。
核心思想：飞行器视为刚体，质量分布不变，形变忽略不计。六自由度 = 三个平动自由度 + 三个转动自由度。
平动动力学
先看平动。牛顿第二定律，F = ma，这个大家都熟。但在飞行器上，力是作用在机体上的，而我们要看的是在地面坐标系下的运动。
我遇到过不少新手，直接把机体坐标系下的力拿去算位置，结果仿真跑起来飞机横着飞。嗯，这里要注意坐标系变换。
平动方程的标准形式：
m * (dV/dt + ω × V) = F_gravity + F_aero + F_thrust其中：
- m — 飞行器质量（常数，不考虑燃料消耗时）
- V — 速度向量（机体坐标系下）
- ω — 角速度向量
- F_gravity — 重力（地面坐标系下，需转换到机体）
- F_aero — 气动力（升力、阻力、侧力）
- F_thrust — 发动机推力
你想想看，为什么会有 ω × V 这一项？因为机体坐标系在旋转。我在做第一个六自由度模型时，忘了加这一项，结果仿真出来的爬升率怎么都对不上。后来查了半天，才发现是哥氏加速度没算。
实战技巧：建议先把重力单独拎出来算，再叠加气动力和推力。这样调试时容易定位问题。我曾经用这个办法，半小时就找到了一个气动系数符号错误。
转动动力学
转动部分用的是欧拉方程。说白了就是角动量定理在旋转坐标系下的表达。
欧拉方程长这样：
I * dω/dt + ω × (I * ω) = M这里：
- I — 惯性张量（3×3矩阵）
- ω — 角速度向量 [p, q, r]ᵀ
- M — 合外力矩
为什么会有 ω × (Iω) 这一项？因为飞行器在转，惯性主轴也在跟着转。这个耦合项是造成飞机横滚-偏航耦合的根本原因。我记得有一次做横滚机动仿真，不加这一项，滚转响应快得像战斗机；加上之后，才看到真实的荷兰滚模态。
注意：欧拉方程在机体坐标系下写最方便。但要注意，惯性张量 I 是常数——前提是你把坐标系固连在机体上，且质量分布不变。如果挂载变化了，I 要重新算。
惯性张量与力矩
惯性张量，说白了就是描述「质量怎么分布」的。对角线是绕各轴的转动惯量，非对角线是惯性积。
标准形式：
I = | Ixx -Ixy -Ixz |
| -Ixy Iyy -Iyz |
| -Ixz -Iyz Izz |对于大多数固定翼无人机，对称面是 XZ 平面，所以 Ixy = Iyz = 0。但 Ixz 不一定为零，这会导致横滚和偏航耦合。
我在做一款飞翼布局时，Ixz 特别大，结果横滚指令下去，飞机自己就开始偏航。后来在控制律里加了解耦补偿才搞定。
力矩来源主要有三个：
- 气动力矩 — 由迎角、侧滑角、舵面偏转产生
- 推力力矩 — 发动机安装位置偏离重心时产生
- 陀螺力矩 — 螺旋桨或涡轮旋转产生的进动效应
避坑指南：我曾经忽略过陀螺力矩，结果仿真大迎角机动时，横滚响应总是慢半拍。加上之后才匹配上试飞数据。对于螺旋桨飞机，这个力矩不能省。
知识体系总览
下面这张图，是我做仿真时经常贴在墙上的。它把刚体动力学的核心逻辑串起来了：
这张图你看懂了吗？从上往下，牛顿-欧拉方程分出平动和转动两条线。平动关心力，转动关心力矩。而惯性张量是连接两者的桥梁——它决定了同样大小的力矩，会产生多大的角加速度。
仿真实现要点
实际写代码时，我建议按这个顺序来：
- 初始化 — 设置质量 m、惯性张量 I、初始位置和姿态
- 力与力矩计算 — 根据当前状态算气动力/力矩、推力/力矩
- 加速度求解 — 用牛顿-欧拉方程解出线加速度和角加速度
- 积分更新 — 用数值积分（四阶龙格-库塔最稳）更新速度和位置
我的习惯：积分步长取 0.001 秒，对于大多数固定翼无人机足够了。如果算力紧张，0.005 秒也能接受，但要注意高频动态可能会失真。
嗯，刚体动力学这块就聊到这儿。记住一句话：力决定怎么飞，力矩决定怎么转，惯性张量决定转得有多「费劲」。搞懂了这三者的关系，六自由度仿真就拿下了一半。
公众号:蓝海资料掘金营，微信deep3321