飞行器空间模拟-F-16飞行动力学模型代码分析

https://github.com/EthanJamesLew/f16-flight-dynamics

一. 函数模块

1. 工具函数模块 (Utility Functions)

1.1 基础数学函数

inline int fix(double ele)
inline int sign(double ele)

fix函数: 根据数值符号进行向下取整(正数)或向上取整(负数)，返回整数
sign函数: 返回数值的符号(-1, 0, 1)

2. 查找表模块 (Lookup Tables)

2.1 气动系数查找表

cl_a[12][7]: 滚转力矩系数查找表 (Rolling moment coefficient)
cm_a[12][5]: 俯仰力矩系数查找表 (Pitching moment coefficient)
cn_a[12][7]: 偏航力矩系数查找表 (Yawing moment coefficient)
cx_a[12][5]: X轴气动力系数查找表 (X-axis aerodynamic force coefficient)
cz_a[12]: Z轴气动力系数查找表 (Z-axis aerodynamic force coefficient)

2.2 阻尼系数查找表

dampp_a[12][9]: 各种阻尼系数查找表

2.3 控制面效应查找表

dlda_a[12][7]: 副翼对滚转力矩的影响
dldr_a[12][7]: 方向舵对滚转力矩的影响
dnda_a[12][7]: 副翼对偏航力矩的影响
dndr_a[12][7]: 方向舵对偏航力矩的影响

2.4 推力查找表

thrust_a[6][6]: 怠速推力表
thrust_b[6][6]: 军用推力表
thrust_c[6][6]: 最大推力表

3. 大气数据计算模块 (Atmospheric Data Computer)

3.1 adc函数

adc_return adc(double vt, double alt)

输入: 真空速(vt)、高度(alt)
输出: 马赫数(amach)、动压(qbar)
功能: 根据标准大气模型计算空气密度、声速、马赫数和动压

计算过程:

温度修正因子计算
根据高度计算温度和密度
计算当地声速
计算马赫数和动压

4. 气动系数计算模块 (Aerodynamic Coefficients)

4.1 力矩系数函数

cl函数 – 滚转力矩系数

double cl(double alpha, double beta)

输入: 攻角(alpha)、侧滑角(beta)
功能: 通过双线性插值计算滚转力矩系数

cm函数 – 俯仰力矩系数

double cm(double alpha, double el)

输入: 攻角(alpha)、升降舵偏角(el)
功能: 计算俯仰力矩系数

cn函数 – 偏航力矩系数

double cn(double alpha, double beta)

输入: 攻角(alpha)、侧滑角(beta)
功能: 计算偏航力矩系数

4.2 力系数函数

cx函数 – X轴力系数

double cx(double alpha, double el)

输入: 攻角(alpha)、升降舵偏角(el)
功能: 计算X轴气动力系数(阻力相关)

cy函数 – Y轴力系数(侧力)

double cy(double beta, double ail, double rdr)

输入: 侧滑角(beta)、副翼偏角(ail)、方向舵偏角(rdr)
功能: 计算侧力系数，使用线性组合

cz函数 – Z轴力系数

double cz(double alpha, double beta, double el)

输入: 攻角(alpha)、侧滑角(beta)、升降舵偏角(el)
功能: 计算Z轴气动力系数(升力相关)

5. 阻尼系数模块 (Damping Module)

5.1 dampp函数

boost::array<double, 9> dampp(double alpha)

输入: 攻角(alpha)
输出: 9个阻尼系数的数组
功能: 根据攻角插值计算各种阻尼系数

6. 控制导数模块 (Control Derivatives)

6.1 副翼效应函数

dlda: 副翼对滚转力矩的影响
dnda: 副翼对偏航力矩的影响

6.2 方向舵效应函数

dldr: 方向舵对滚转力矩的影响
dndr: 方向舵对偏航力矩的影响

特点: 所有控制导数函数都使用相同的双线性插值算法，基于攻角和侧滑角进行插值。

7. 推力模型模块 (Thrust Model)

7.1 推力相关函数

tgear函数 – 油门-功率转换

double tgear(double thtl)

功能: 将油门位置转换为发动机功率百分比
特点: 分段线性函数，0.77为分界点

rtau函数 – 时间常数计算

inline double rtau(double dp)

功能: 计算功率变化的时间常数

pdot函数 – 功率变化率

double pdot(double p3, double p1)

功能: 计算发动机功率的变化率，模拟发动机动态响应

thrust函数 – 推力计算

double thrust(double power, double alt, double rmach)

输入: 功率、高度、马赫数
功能: 通过三维插值计算发动机推力
特点: 根据功率范围选择不同的推力表进行插值

8. Morelli模型模块 (Advanced Aerodynamic Model)

8.1 morelli函数

boost::array<double, 6> morelli(各种参数...)

功能: 高级气动模型，计算6个气动系数
输入: 攻角、侧滑角、控制面偏角、角速度等
输出: [Cx, Cy, Cz, Cl, Cm, Cn]
特点: 使用多项式拟合，考虑重心位置影响

9. 插值算法分析

9.1 双线性插值流程

大部分气动系数计算函数使用相同的插值算法：

索引计算: 将输入参数映射到查找表索引
边界处理: 确保索引在有效范围内
权重计算: 计算插值权重
双线性插值: 在两个维度上进行线性插值
符号处理: 根据侧滑角符号调整结果

9.2 常数定义

SIGN_FACTOR: 1.1f，用于插值计算中的符号处理

10. 技术特点总结

数据驱动: 大量使用查找表存储实验或CFD数据
高效插值: 统一的插值算法保证计算效率
物理建模: 考虑了重心位置、阻尼效应等物理因素
模块化设计: 不同功能模块清晰分离
工程实用: 基于Stevens & Lewis经典教材的成熟模型

二. F-16飞行动力学仿真模型分析

概述

这是一个F-16战斗机的六自由度飞行动力学仿真模型，基于Stevens & Lewis的飞行动力学理论实现。F-16是一款不稳定的”电传飞控”飞机，在无控制状态下会表现出不稳定行为。

核心结构体参数分析

F16PlantParameters 飞机物理参数

参数	数值	单位	物理意义
`xcg`	0.35	–	重心位置（相对于平均气动弦长的比例）
`s`	300.0	ft²	机翼参考面积
`b`	30.0	ft	机翼翼展
`cbar`	11.32	ft	平均气动弦长
`rm`	1.57e-3	slug/ft³	空气密度相关常数
`xcgr`	0.35	–	参考重心位置
`he`	160.0	ft·lb·s	发动机陀螺效应参数
`c1~c9`	各异	–	惯性矩和产品惯性矩的组合系数
`rtod`	57.29578	–	弧度转角度系数（180/π）
`g`	32.17	ft/s²	重力加速度

状态变量定义

模型使用13个状态变量描述飞机运动：

基本运动状态（0-8）

x[0] – vt: 真空速（ft/s）
x[1] – alpha: 攻角（弧度）
x[2] – beta: 侧滑角（弧度）
x[3] – phi: 滚转角（弧度）
x[4] – theta: 俯仰角（弧度）
x[5] – psi: 偏航角（弧度）
x[6] – p: 滚转角速度（弧度/s）
x[7] – q: 俯仰角速度（弧度/s）
x[8] – r: 偏航角速度（弧度/s）

位置与推力状态（9-12）

x[9-10] – 水平位置坐标
x[11] – alt: 高度（ft）
x[12] – power: 发动机功率状态

控制输入变量

uinput[0] – thtlc: 油门杆位置（0-1）
uinput[1] – el: 升降舵偏转角（度）
uinput[2] – ail: 副翼偏转角（度）
uinput[3] – rdr: 方向舵偏转角（度）

核心功能模块

1. 气动力计算

// 两种气动力模型可选
if (model_type == STEVENS) {
    // Stevens & Lewis 查表法
    cxt = cx(alpha, el);  // 轴向力系数
    cyt = cy(beta, ail, rdr);  // 侧向力系数
    czt = cz(alpha, beta, el);  // 法向力系数
} else if (model_type == MORELLI) {
    // Morelli 模型（更精确）
    // 计算六个气动力和力矩系数
}

2. 力方程求解

// 机体坐标系下的力方程
udot = r*v - q*w - g*sth + rm*(qs*cxt + t);
vdot = p*w - r*u + gcth*sph + ay;
wdot = q*u - p*v + gcth*cph + az;

3. 运动学方程

// 欧拉角变化率
dxdt[3] = p + (sth/cth) * (qsph + r*cph);  // 滚转角速度
dxdt[4] = q*cph - r*sph;                   // 俯仰角速度
dxdt[5] = (qsph + r*cph) / cth;            // 偏航角速度

4. 力矩方程

// 三个力矩方程，描述角加速度
dxdt[6] = (c2*p + c1*r + c4*he)*q + qsb*(c3*clt + c4*cnt);
dxdt[7] = (c5*p - c7*he)*r + c6*(r*r - p*p) + qs*cbar*c7*cmt;
dxdt[8] = (c8*p - c2*r + c9*he)*q + qsb*(c4*clt + c9*cnt);

输出计算

加速度输出

Nz: 法向过载（g）
Ny: 侧向过载（g）
az, ay: 法向和侧向加速度（ft/s²）

主要作用和应用

1. 飞行仿真

提供F-16飞机的高保真度飞行动力学模型
支持实时飞行仿真和训练

技术特点

高精度建模: 采用非线性六自由度模型，考虑了气动力、推力、重力等所有主要作用力
多模型支持: 支持Stevens和Morelli两种气动力模型，可根据精度需求选择
实时仿真: 结构化设计支持实时仿真应用
模块化设计: 继承自ComponentBase，便于集成到更大的仿真系统中

三. 飞行控制系统

主要公共函数

1. 构造函数

LowLevelController() – 默认构造函数，使用预设的平衡状态和控制输入
LowLevelController(f16_state_type state_trim, llc_output_type control_trim) – 自定义构造函数，允许指定平衡状态和控制配平

2. 核心控制函数

void output(double time, const llc_state_type &state, llc_output_type &output, const llc_input_type &input)
- 功能: 主要的控制输出函数，基于LQR反馈控制计算控制面偏转角度
- 输入: 时间、控制器状态、输入（包含F-16状态和参考控制指令）
- 输出: 控制输出（油门、升降舵、副翼、方向舵）
- 特点: 包含控制限幅保护
void update(double time, const llc_state_type &state, llc_state_type &state_up, const llc_input_type &input)
- 功能: 更新控制器内部状态
- 计算: 法向载荷因子误差、滚转速率、偏航速率误差

3. 辅助函数

void printMatrix(const matrix<double> &m) – 矩阵打印函数，用于调试

主要私有/内联函数

4. 内部计算函数

inline matrixd reorder_xctrl(const f16_state_type &state, const llc_state_type &cstate)
- 功能: 重新排列状态向量以匹配控制器格式
- 返回: 8×1的状态矩阵
inline matrixd build_lqr()
- 功能: 构建LQR控制增益矩阵
- 返回: 3×8的增益矩阵
- 特点: 包含CSAF配置的增益参数

关键数据结构

5. 控制限制结构体

CtrlLimits – 定义各控制面的最大最小限制值
- 油门: 0.0 – 1.0
- 升降舵: ±25.0度
- 副翼: ±21.5度
- 方向舵: ±30.0度
- 滚转角限制: ±60.0度
- 法向载荷因子: -1.0 到 6.0

6. 状态索引枚举

F16StateIdxs – F-16状态变量的索引定义，包含速度、姿态角、角速度、位置、功率等

使用场景

这个低级控制器主要用于：

F-16飞机的姿态控制
基于LQR算法的反馈控制
控制面偏转的计算和限幅
飞行状态的实时更新

核心的使用流程是通过 output() 函数计算控制输出，通过 update() 函数更新控制器状态。