MATLAB】【Simulink】【固定翼无人机】【六自由度方程】【飞行动力学】【控制仿真】

这是一个基于 Simulink 的 固定翼无人机飞行动力学仿真程序，核心采用 六自由度非线性运动方程。已实现完整的 气动力计算、力矩计算、状态方程，支持固定翼无人机在不同工况下的飞行模拟（滑翔 / 转弯等）。
功能特色

▷ 非线性 6-DOF 方程：涵盖纵向、横向、方向、运动学方程
▷ 气动力建模：升力、阻力、侧力、俯仰/偏航/滚转力矩计算
▷ 状态变量齐全：速度、航迹角、攻角、滚转角、航向角、位置/高度等
▷ 输入可控：油门、升降舵、方向舵、副翼偏角均可设定
▷ Simulink 接口：可嵌入 Simulink 仿真框架，支持进一步控制律设计
▷ 注释清晰：物理量单位与来源说明详细，便于二次开发和论文使用

仿真亮点

→ 支持 不同机动模式（滑翔、转弯）
→ 参数包含 机翼/尾翼/机身气动特性，参考实际气动数据库
→ 易于扩展：可加入扰动、大气模型、控制器模块等
→ 科研/课程友好：适合仿真、飞行力学课程

内容完整性：一个完整的六自由度飞行动力学仿真系统包含多个模块（如气动力查表、质量特性、坐标系转换、推进系统、控制律等），
以下示例

1. 六自由度方程核心框架（Matlab伪代码）

matlab
% 六自由度非线性运动方程 - 简化版示意
function dxdt = rigidBodyEOM(t, x, u, params)
% 状态变量 x: [pn, pe, h, u, v, w, phi, theta, psi, p, q, r]
% 控制输入 u: [delta_e, delta_a, delta_r, delta_t] (升降舵、副翼、方向舵、油门)

% 解包状态
pn = x(1); pe = x(2); h = x(3);
u_body = x(4); v_body = x(5); w_body = x(6);
phi = x(7); theta = x(8); psi = x(9);
p = x(10); q = x(11); r = x(12);

% 计算空速、攻角、侧滑角
Vt = sqrt(u_body^2 + v_body^2 + w_body^2);
alpha = atan2(w_body, u_body);
beta = asin(v_body / Vt);

% 气动力/力矩系数计算（可来自查表或经验公式）
[Cx, Cy, Cz, Cl, Cm, Cn] = aerodynamicCoefficients(alpha, beta, p, q, r, u);

% 计算力和力矩
Q = 0.5 params.rho Vt^2;
L = -Q params.S Cz; % 升力
D = -Q params.S Cx; % 阻力
Y = -Q params.S Cy; % 侧向力

Ml = Q params.S params.b Cl; % 滚转力矩
Mm = Q params.S params.c Cm; % 俯仰力矩
Mn = Q params.S params.b Cn; % 偏航力矩

% 牛顿第二定律：F = ma 和 M = I alpha
% 平动方程（地面坐标系）
pn_dot = (cos(theta)cos(psi)u_body + …
(sin(phi)sin(theta)cos(psi) - cos(phi)sin(psi))v_body + …
(cos(phi)sin(theta)cos(psi) + sin(phi)sin(psi))w_body);

pe_dot = (cos(theta)sin(psi)u_body + …
(sin(phi)sin(theta)sin(psi) + cos(phi)cos(psi))v_body + …
(cos(phi)sin(theta)sin(psi) - sin(phi)cos(psi))w_body);

h_dot = (-sin(theta)u_body + sin(phi)cos(theta)v_body + cos(phi)cos(theta)w_body);

% 机体加速度（考虑重力和外力）
u_dot = rv_body - qw_body + Lsin(alpha) + Dcos(alpha) + params.gsin(theta) + …;
v_dot = pw_body - ru_body + Y + params.g(sin(phi)cos(theta)) + …;
w_dot = qu_body - pv_body - Lcos(alpha) + Dsin(alpha) + params.g(cos(phi)cos(theta)) + …;

% 转动方程
omega = [p; q; r];
I_body = params.I; % 惯性张量
M = [Ml; Mm; Mn];
omega_dot = I_body \ (M - cross(omega, I_body omega));

p_dot = omega_dot(1);
q_dot = omega_dot(2);
r_dot = omega_dot(3);

% 姿态更新（欧拉角微分方程）
phi_dot = p + qsin(phi)tan(theta) + rcos(phi)tan(theta);
theta_dot = qcos(phi) - rsin(phi);
psi_dot = (qsin(phi) + r*cos(phi)) / cos(theta);

dxdt = [pn_dot; pe_dot; h_dot; u_dot; v_dot; w_dot; …
phi_dot; theta_dot; psi_dot; p_dot; q_dot; r_dot];
end

2. Simulink 模型建议结构
   您可以构建如下模块：
   State Integrator: 使用 Integrator 模块积分状态变量。
   Kinematics & Dynamics: 分别实现运动学和动力学方程。
   Aerodynamics Subsystem: 查表 lookup tables 实现 (C_L(\alpha)), (C_D(\alpha)), (C_m(\alpha)) 等。
   Atmosphere Model: 使用标准大气模型（如 COESA）。
   Propulsion: 推力模型（与油门、空速相关）。
   Actuator Models: 一阶惯性环节模拟舵机响应。

您可以按照以下结构在 Simulink 中搭建模型：

Text
编辑
[Input: 控制输入 (δe, δa, δr, throttle)]
↓
[Propulsion Model] → [Thrust Force]
↓
[Aerodynamics Subsystem]
├─ α, β, p, q, r → [Lookup Tables] → CL, CD, CY, Cl, Cm, Cn
└─ Q = 0.5 * ρ * V² → 力和力矩计算
↓
[Forces & Moments] → [Newton-Euler Equations]
↓
[Integrator Block] → [State Vector: pn, pe, h, u,v,w, φ,θ,ψ, p,q,r]
↓
[Coordinate Transformations]
├─ Body to Earth: Velocity → Position (pn, pe, h)
└─ Euler Angles → Direction Cosine Matrix
↓
[Output: State Plots / Scope]

Matlab
编辑
function [CL, CD, CM] = aerodynamicCoefficients(alpha, beta, p, q, r, u)
% 示例：基于攻角α查表计算升力、阻力、俯仰力矩系数
% 假设数据来自实际机型

alpha_deg = alpha * 180 / pi;

% 查表数据（示例）
alpha_tab = [-10:1:10]; % 攻角范围（度）
CL_tab = 0.1 * alpha_tab + 0.2; % 简化线性关系
CD_tab = 0.01 + 0.005 * alpha_tab.^2;
CM_tab = -0.04 * alpha_tab;

% 插值
CL = interp1(alpha_tab, CL_tab, alpha_deg, ‘linear’, ‘extrap’);
CD = interp1(alpha_tab, CD_tab, alpha_deg, ‘linear’, ‘extrap’);
CM = interp1(alpha_tab, CM_tab, alpha_deg, ‘linear’, ‘extrap’);

% 侧滑角影响（简化）
CY = 0.1 * beta;
Cl = -0.05 * beta;
Cn = 0.03 * beta;

end
📌 六自由度运动方程（S-Function 或 MATLAB Function 模块中使用）
Matlab
编辑
function dxdt = sixDOF(t, x, u, params)
% x: [pn, pe, h, u, v, w, phi, theta, psi, p, q, r]
% u: [delta_e, delta_a, delta_r, throttle]

% 解包
pn = x(1); pe = x(2); h = x(3);
u_b = x(4); v_b = x(5); w_b = x(6);
phi = x(7); theta = x(8); psi = x(9);
p = x(10); q = x(11); r = x(12);

% 计算空速
Vt = sqrt(u_b^2 + v_b^2 + w_b^2);

% 攻角与侧滑角
alpha = atan2(w_b, u_b);
beta = asin(v_b / Vt);

% 气动力计算
[Cx, Cy, Cz, Cl, Cm, Cn] = aerodynamicCoefficients(alpha, beta, p, q, r, u_b);

% 密度（标准大气）
rho = 1.225; % kg/m³（海平面）

% 力计算
Q = 0.5 * rho * Vt^2;
L = -Q * params.S * Cz; % 升力
D = -Q * params.S * Cx; % 阻力
Y = -Q * params.S * Cy; % 侧向力

% 力矩计算
Ml = Q * params.S * params.b * Cl;
Mm = Q * params.S * params.c * Cm;
Mn = Q * params.S * params.b * Cn;

% 推力（简化）
T = params.T_max * u(4); % throttle

% 质量与惯性
m = params.mass;
Ixx = params.I(1); Iyy = params.I(2); Izz = params.I(3);

% 平动方程（地面坐标系）
% 速度分量转换到地面系
R = rotationMatrix(phi, theta, psi);
v_earth = R * [u_b; v_b; w_b];

% 地面坐标系加速度
acc_north = R(1,1)*dudt + R(1,2)*dvdt + R(1,3)*dwdt;
acc_east = R(2,1)*dudt + R(2,2)*dvdt + R(2,3)*dwdt;
acc_down = R(3,1)*dudt + R(3,2)*dvdt + R(3,3)*dwdt;

% 姿态微分方程（欧拉角）
phi_dot = p + qsin(phi)tan(theta) + rcos(phi)tan(theta);
theta_dot = qcos(phi) - rsin(phi);
psi_dot = (qsin(phi) + rcos(phi)) / cos(theta);

% 角加速度（转动方程）
omega = [p; q; r];
M_vec = [Ml; Mm; Mn];
omega_dot = I_body \ (M_vec - cross(omega, I_body * omega));

% 输出导数
dxdt = [acc_north; acc_east; acc_down; …
dudt; dvdt; dwdt; …
phi_dot; theta_dot; psi_dot; …
omega_dot(1); omega_dot(2); omega_dot(3)];
end