3、飞行动力学基础：刚体运动方程、牛顿-欧拉方程、空气动力学简化模型

好，咱们进入第三讲。说实话，飞行动力学是很多做控制的人觉得「劝退」的一章。为什么？因为公式多、概念杂，而且看起来跟控制算法没啥直接关系。

但我得说句实话：你如果不懂飞机怎么动，你写的控制律就是空中楼阁。我在项目里见过太多人，PID参数调得飞起，结果真机一飞就炸——原因很简单，他压根没搞懂飞机在空中的运动规律。

这一章，咱们就把飞行动力学的核心骨架拆出来。不搞复杂的推导，只讲你写代码、调参数真正需要的东西。

3.1 刚体运动方程：飞机到底怎么动？

先问个问题：一架四旋翼在空中，它有几个自由度？

答案是6个。3个平动（前后、左右、上下），3个转动（俯仰、滚转、偏航）。这就是所谓的「六自由度刚体运动」。

我个人习惯把运动方程拆成两部分来看：

• 运动学：描述位置和姿态怎么随时间变化，跟力没关系
• 动力学：描述力和力矩怎么产生加速度和角加速度

先看运动学。位置变化很简单：

// 位置变化率 = 速度
p_dot = v

姿态变化稍微麻烦点。我建议用四元数，别用欧拉角。为什么？欧拉角有万向锁问题，而且插值不线性。我在做倾转旋翼机项目时，就吃过欧拉角的亏——飞到90度俯仰时，姿态解算直接炸了。

四元数运动学方程：

// q = [qw, qx, qy, qz] 是单位四元数
// ω = [p, q, r] 是机体角速度
q_dot = 0.5 * q ⊗ ω

这里 ⊗ 表示四元数乘法。代码实现时，我习惯写成矩阵形式：

void quaternion_kinematics(float q[4], float omega[3], float q_dot[4]) {
    q_dot[0] = -0.5 * (q[1]*omega[0] + q[2]*omega[1] + q[3]*omega[2]);
    q_dot[1] =  0.5 * (q[0]*omega[0] - q[3]*omega[1] + q[2]*omega[2]);
    q_dot[2] =  0.5 * (q[3]*omega[0] + q[0]*omega[1] - q[1]*omega[2]);
    q_dot[3] = -0.5 * (q[2]*omega[0] - q[1]*omega[1] - q[0]*omega[2]);
}

3.2 牛顿-欧拉方程：力和力矩怎么算？

好，运动学搞定了。接下来是动力学——说白了就是牛顿第二定律的空中版。

平动部分（牛顿方程）：

m * a = F_gravity + F_thrust + F_aero

这里要注意坐标系。重力在地面系，推力和气动力在机体系。我建议统一转到地面系再算加速度，这样位置积分方便。

转动部分（欧拉方程）：

J * ω_dot = M - ω × (J * ω)

其中 J 是转动惯量矩阵，M 是合力矩。这个 × 是叉乘，在代码里要小心实现：

void euler_dynamics(float J[3][3], float omega[3], float M[3], float omega_dot[3]) {
    float J_omega[3];
    mat_vec_mul(J, omega, J_omega);  // J * ω
    
    float cross[3];
    cross_product(omega, J_omega, cross);  // ω × (J*ω)
    
    float temp[3];
    for(int i=0; i<3; i++) temp[i] = M[i] - cross[i];
    
    // 解 J * ω_dot = temp
    solve_linear_system(J, temp, omega_dot);
}

3.3 空气动力学简化模型：别被复杂公式吓到

说到空气动力学，很多做控制的同学就头疼。其实对于飞行控制来说，你不需要搞CFD仿真，只需要一个能用的简化模型。

我个人常用的简化思路：

• 升力：跟迎角成正比，跟速度平方成正比
• 阻力：跟速度平方成正比，方向与速度相反
• 力矩：跟角速度成正比（阻尼项）

对于固定翼，升力系数可以近似为：

CL = CL0 + CL_alpha * alpha

其中 alpha 是迎角。对于多旋翼，更简单——升力直接跟螺旋桨转速平方成正比：

F_thrust = k_f * omega_motor^2
M_torque = k_m * omega_motor^2

这里 k_f 和 k_m 是螺旋桨的力系数和力矩系数，一般通过实验标定得到。

参数	物理意义	典型值（四旋翼）
k_f	力系数 (N/(rad/s)^2)	1.0e-5 ~ 1.0e-4
k_m	力矩系数 (Nm/(rad/s)^2)	1.0e-7 ~ 1.0e-6
CL_alpha	升力线斜率 (1/rad)	4.0 ~ 6.0

3.4 把这些方程串起来：仿真框架

好了，现在我们把所有东西拼在一起。一个完整的飞行动力学仿真，大概长这样：

// 伪代码：一个控制周期
void dynamics_step(float dt) {
    // 1. 计算力和力矩
    float F[3], M[3];
    compute_forces_and_torques(state, controls, &F, &M);
    
    // 2. 平动动力学（地面系）
    float acc[3];
    acc[0] = F[0] / mass + gravity[0];  // 包含重力
    acc[1] = F[1] / mass + gravity[1];
    acc[2] = F[2] / mass + gravity[2];
    
    // 3. 转动动力学（机体系）
    float omega_dot[3];
    euler_dynamics(J, state.omega, M, omega_dot);
    
    // 4. 积分更新
    state.vel += acc * dt;
    state.pos += state.vel * dt;
    state.omega += omega_dot * dt;
    
    // 5. 姿态更新（四元数）
    quaternion_kinematics(state.q, state.omega, q_dot);
    state.q += q_dot * dt;
    normalize_quaternion(&state.q);
}

这个框架看起来简单，但实际用的时候有几个坑：

• 积分方法：我用的是最简单的欧拉法，但如果你要精度高，建议用四阶龙格-库塔
• 坐标系转换：力和力矩一定要搞清楚是在哪个坐标系下算的
• 数值稳定性：dt 太大时，欧拉法会发散。我一般控制在 1ms 以内

嗯，这一章内容不少。但说白了，飞行动力学就是「力产生加速度，加速度积分成速度，速度积分成位置」——只不过在三维空间里，多了旋转和坐标系变换。

下一章我们会讲控制算法本身。但在这之前，我建议你把这一章的仿真框架自己写一遍。代码跑通了，你对飞机怎么动就有了直觉。这个直觉，比背一百个公式都管用。