3. 六自由度动力学方程：牛顿-欧拉方程推导
好，咱们进入正题。六自由度动力学方程，说白了就是回答一个问题：飞行器为什么会动？ 以及 它怎么动？
我个人习惯把这个问题拆成两半来看：平移 和 旋转。平移靠牛顿第二定律，旋转靠欧拉方程。合起来，就是经典的牛顿-欧拉法。
3.1 牛顿-欧拉方程推导
先看平移。牛顿第二定律大家都熟：
F = m * a但在飞行器上，这个 a 不是随便取的。你得注意坐标系。
我一般把力 F 分解到机体坐标系（Body Frame）上。但加速度 a 必须是在惯性系（Inertial Frame）下测量的。为什么？因为牛顿定律只对惯性系成立。
所以，实际写方程时，我们得把机体坐标系下的力，转换到惯性系下。或者反过来，把惯性系下的加速度投影到机体轴上。
嗯，这里有个坑。我曾经在项目里直接拿机体坐标系下的加速度去套牛顿第二定律，结果算出来的轨迹全是歪的。后来才发现，忘了考虑哥氏加速度和向心加速度。
完整的平移动力学方程，在机体坐标系下写成这样：
m * (dV/dt + ω × V) = F_ext其中：
- V 是速度向量（在机体坐标系下表示）
- ω 是角速度向量
- F_ext 是合外力（重力、推力、气动力等）
- × 表示叉乘
这个 ω × V 项，就是哥氏加速度。你想想看，飞行器一边转一边往前飞，这两个运动耦合在一起，就会产生额外的惯性力。
3.2 平移与旋转动力学
说完平移，再看旋转。旋转动力学用的是欧拉方程：
I * dω/dt + ω × (I * ω) = M_ext这里 I 是惯性张量，M_ext 是合外力矩。
这个方程看着简单，但实际用起来，I 矩阵的处理是个大麻烦。我刚开始做飞控时，以为惯性张量就是个常数，直接写死。结果飞行器一做大机动，模型就崩了。
为什么会这样？因为惯性张量是相对于质心定义的。如果飞行器的质量分布发生变化（比如燃料消耗、载荷移动），I 就会变。
所以，我建议你在做参数辨识时，一定要把惯性张量作为待辨识参数，而不是直接查手册。
核心要点：
- 平移方程：m * (dV/dt + ω × V) = F
- 旋转方程：I * dω/dt + ω × (I * ω) = M
- 两个方程通过 ω 和 V 耦合在一起
3.3 惯性张量与质量特性
惯性张量 I 是个 3×3 的对称矩阵：
I = | Ixx -Ixy -Ixz |
| -Ixy Iyy -Iyz |
| -Ixz -Iyz Izz |对角线元素是绕各轴的转动惯量，非对角线元素是惯性积。
我见过不少同学，直接假设惯性积为零。这在理想对称的飞行器上勉强成立，但实际工程中，几乎没有完全对称的飞行器。
举个例子。我曾经调试一架固定翼无人机，发现它滚转时总会带出偏航。查了半天，发现是机翼两侧的油箱不对称，导致 Ixz 不为零。滚转力矩通过惯性积耦合到了偏航轴上。
所以，惯性积不是小量，不能随便忽略。
| 参数 | 物理意义 | 常见取值（小型无人机） |
|---|---|---|
| Ixx | 绕X轴转动惯量 | 0.01 ~ 0.1 kg·m² |
| Iyy | 绕Y轴转动惯量 | 0.01 ~ 0.1 kg·m² |
| Izz | 绕Z轴转动惯量 | 0.02 ~ 0.2 kg·m² |
| Ixy, Ixz, Iyz | 惯性积 | 通常为 Ixx 的 1%~10% |
我的经验： 惯性张量可以通过 CAD 模型估算，但最好用实验辨识。常用的方法有：
- 扭摆法：测量绕单轴的转动惯量
- 三线摆法：测量绕垂直轴的转动惯量
- 参数辨识法：通过飞行数据反推惯性张量
3.4 完整六自由度EOM
把平移和旋转方程合在一起，再加上运动学关系，就得到了完整的六自由度方程（EOM）。
我习惯把 EOM 写成状态空间形式：
状态向量：x = [x, y, z, φ, θ, ψ, u, v, w, p, q, r]^T
其中：
- x, y, z：位置（惯性系）
- φ, θ, ψ：姿态角（欧拉角）
- u, v, w：速度（机体坐标系）
- p, q, r：角速度（机体坐标系）
动力学方程：
- 平移：m * (du/dt + q*w - r*v) = Fx
m * (dv/dt + r*u - p*w) = Fy
m * (dw/dt + p*v - q*u) = Fz
- 旋转：Ixx*dp/dt - Ixy*dq/dt - Ixz*dr/dt + (Izz-Iyy)*q*r + Ixz*p*q - Iyz*q*q + Ixy*p*r - Ixz*r*r = Mx
...（类似，共三个方程）
运动学方程：
- 位置：dx/dt = u*cosθ*cosψ + v*(sinφ*sinθ*cosψ - cosφ*sinψ) + w*(cosφ*sinθ*cosψ + sinφ*sinψ)
...（类似，共三个方程）
- 姿态：dφ/dt = p + q*sinφ*tanθ + r*cosφ*tanθ
dθ/dt = q*cosφ - r*sinφ
dψ/dt = (q*sinφ + r*cosφ) / cosθ看着很复杂对吧？别怕。实际做参数辨识时，我们通常不会直接用这么完整的方程。而是根据飞行器的特点做简化。
注意： 欧拉角在俯仰角 θ = ±90° 时会出现奇异性（万向锁）。如果你的飞行器会做大角度机动，建议改用四元数表示姿态。
我曾经在一个垂直起降项目中吃过这个亏。飞行器从悬停转平飞时，俯仰角接近90度，欧拉角直接发散，仿真崩了。后来换成四元数，问题才解决。
3.5 知识体系总览
下面这张图，是我自己总结的六自由度动力学知识框架。你可以把它当作一个「地图」，随时回来对照。
这张图把整个章节的逻辑串起来了。从上到下，从牛顿-欧拉法出发，分支出平移和旋转两条线，最后汇合到完整的 EOM。你写代码做仿真时，也可以按这个结构来组织代码。
本章小结：
- 六自由度动力学 = 平移动力学 + 旋转动力学
- 平移方程别忘了哥氏加速度项
- 旋转方程的核心是惯性张量，惯性积不能随便忽略
- 完整 EOM 是 12 阶状态空间，但实际应用时可简化
- 大角度机动时，用四元数代替欧拉角
好了，这一章就到这里。下一章我们会聊怎么用实验数据来辨识这些参数。到时候，我会分享一些我在实际项目中踩过的坑，以及怎么用最小二乘法把参数「算」出来。