2. 飞行动力学建模的物理基础

飞行动力学实时模型之所以能“算得准、跑得快”，前提是把物理规律写成计算机可迭代求解的形式。这一节围绕四件事展开：坐标系怎么定、六自由度方程怎么来、牛顿-欧拉方程如何落地、气动力与力矩从哪里产生。四者是一条从几何描述到力学因果的链条，缺一环模型就会失真。

一、坐标系：先约定，再计算

建模中最常见的错误不是公式写错，而是坐标系约定不统一。工程上通常使用三套右手坐标系：

• 地面坐标系 O_eX_eY_eZ_e：原点取地面某固定点，X_e指北、Y_e指东、Z_e垂直向下。用于描述航迹、位置和重力方向。
• 机体系 O_bX_bY_bZ_b：原点取飞机重心，X_b指机头、Y_b指右翼、Z_b指机腹。用于描述角速度、气动力和力矩。
• 气流坐标系（速度坐标系）：X_a沿空速矢量，Z_a在对称面内垂直于X_a向下。用于定义升力、阻力和侧力。

三套坐标系之间通过欧拉角连接：偏航角ψ、俯仰角θ、滚转角φ描述机体系相对地面系的姿态；迎角α、侧滑角β描述气流系相对机体系的方向。转换矩阵按“先偏航、再俯仰、后滚转”的Z-Y-X顺序构造，机体系到地面系的旋转矩阵为：

R_b2e = Rz(ψ)·Ry(θ)·Rx(φ)

其中：
Rx(φ) = [1, 0, 0; 0, cosφ, -sinφ; 0, sinφ, cosφ]
Ry(θ) = [cosθ, 0, sinθ; 0, 1, 0; -sinθ, 0, cosθ]
Rz(ψ) = [cosψ, -sinψ, 0; sinψ, cosψ, 0; 0, 0, 1]

实时模型中，姿态更新通常不直接对欧拉角积分，而是用四元数q=[q0,q1,q2,q3]推进，再转回欧拉角输出，避免θ=±90°时的万向锁。四元数微分方程为 q̇ = 0.5·q⊗[0, p, q, r]，其中p、q、r是机体系角速度分量。这一步是姿态解算的核心，也是实时性最容易出问题的地方。

二、六自由度方程：从牛顿-欧拉落地

刚体六自由度运动方程由牛顿第二定律和角动量定理在机体系下展开得到。平动方程中，m·(V̇b + ω×Vb) = Fb，其中ω×Vb是机体系旋转带来的附加项，不能省略；转动方程中，I·ω̇ + ω×(I·ω) = Mb，ω×(I·ω)是陀螺耦合项，在大角速率机动时影响显著。

以某轻型飞机为例，设质量m=1200 kg，惯量矩阵I=diag(1800, 3500, 4200) kg·m²，巡航状态V=60 m/s、α=4°、β=0°。若在t=0时刻施加δe=-5°的升降舵阶跃，俯仰力矩系数Cm≈-0.35，动压q̄=0.5×1.225×60²≈2205 Pa，S=16 m²，c̄=1.5 m，则俯仰力矩Mm≈2205×16×1.5×(-0.35)≈-18522 N·m。由Iyy=3500 kg·m²可得初始俯仰角加速度q̇≈-5.29 rad/s²。这个量级说明：实时模型必须以毫秒级步长积分，否则角速率会迅速发散。

实时实现时，通常把方程写成状态空间形式：ẋ = f(x, u)，x为12维状态向量，u为舵面偏度和油门。积分器可选四阶龙格-库塔或半隐式欧拉。以1 kHz步长为例，每步计算量约几百次浮点运算，在嵌入式飞控计算机上完全可实时。关键是把气动系数表做成查表加线性插值，避免在线求解复杂气动方程。

三、气动力与力矩：从迎角到舵面

气动力和力矩是飞行动力学模型中最“非线性”的部分。升力系数CL主要随迎角α线性增长，斜率约0.1/deg，直到失速迎角（约15°~18°）后急剧下降；阻力系数CD由零升阻力和诱导阻力组成，CD=CD0+k·CL²；俯仰力矩系数Cm随α增大而减小，体现静稳定性。舵面偏转改变有效迎角和弯度，从而改变CL、CD和Cm。

实时模型中，气动数据通常以三维查表形式存储：CL=f(α, Ma, δe)，CD=f(α, Ma, δe)，Cm=f(α, Ma, q, δe)。其中q是俯仰角速率，用于阻尼项。查表时先按Ma和δe定位，再在α方向线性插值。若Ma跨越0.7以上，还需考虑压缩性修正。侧滑角β产生侧力和偏航力矩，方向舵δr用于配平。

四、建模中的常见陷阱

• 坐标系混用：把机体系角速度直接当欧拉角速率积分，导致姿态错误。正确关系是φ̇、θ̇、ψ̇与p、q、r之间通过转换矩阵耦合。
• 忽略陀螺项：大角速率下ω×(I·ω)不可忽略，否则滚转-偏航耦合会失真。
• 气动系数符号：不同文献对Cm正负号约定不同，建模前必须统一，否则升降舵极性反了。
• 积分步长过大：角速率动态时间常数约0.1~0.5 s，步长超过10 ms就可能不稳定。