MATLAB 飞行器等效拟合配平程序

MATLAB 飞行器等效拟合配平程序，适用于常规布局固定翼飞机，支持：

• 静稳定性导数建模
• 多操纵面等效舵面（δ_e, δ_a, δ_r）
• 非线性升力/阻力/力矩拟合
• 指定飞行状态（Ma, H, α）下的全机配平求解
• 配平曲线自动扫描（α–δ_e, Ma–C_L 等）

这套程序介于教材解析法和CFD/风洞数据驱动法之间，是很多高校/院所做飞行动力学、控制律设计前的标准前置工具。

---

一、配平问题数学描述

1、纵向配平方程（力矩平衡）

\[C_m(\alpha,\delta_e,Ma,H) = 0 \]

常用线性化形式（等效舵面）：

\[C_m = C_{m0} + C_{m\alpha}\alpha + C_{m\delta_e}\delta_e + C_{mq}\frac{\bar{c}}{2V}q \]

配平时 (q=0)，简化为：

\[\delta_e^{trim} = -\frac{C_{m0}+C_{m\alpha}\alpha}{C_{m\delta_e}} \]

---

2、横航向配平

\[C_l = C_n = 0 \Rightarrow \begin{cases} \delta_a^{trim} = -\dfrac{C_{l\beta}\beta + C_{lp}p + C_{lr}r}{C_{l\delta_a}}\\ \delta_r^{trim} = -\dfrac{C_{n\beta}\beta + C_{np}p + C_{nr}r}{C_{n\delta_r}} \end{cases}\]

---

二、核心程序结构

Trim_Aircraft/
├── main_trim.m              % 主程序
├── aircraft_data.m          % 飞机参数 & 气动导数
├── aerodyn_lookup.m         % 气动力/力矩查表/拟合
├── trim_equations.m         % 配平方程（fmincon）
├── plot_trim_curves.m       % 配平曲线绘制
└── utils/
    └── atmosISA.m           % 标准大气

---

三、飞机参数与气动导数（aircraft_data.m）

%% aircraft_data.m
function data = aircraft_data()
% 飞机基本参数
data.m     = 8500;        % 质量 kg
data.S     = 28;          % 参考面积 m^2
data.c_bar = 3.2;         % 平均气动弦长 m
data.b     = 10.5;        % 展长 m
data.Iyy   = 42000;       % 俯仰转动惯量 kg·m^2

% 重心位置（MAC 百分比）
data.xcg = 0.25;

% === 气动导数（示例值，可替换为风洞/CFD）===
% 纵向
data.CL0    = 0.25;
data.CLa    = 5.2;        % /rad
data.CLde   = 0.35;       % /rad
data.CD0    = 0.02;
data.K      = 0.045;      % 诱导阻力因子
data.Cm0    = 0.02;
data.Cma    = -0.8;       % /rad
data.Cmde   = -1.1;       % /rad
data.Cmq    = -12;        % /rad

% 横航向
data.Clb    = -0.12;
data.Clda   = 0.08;
data.Cldb   = 0.03;
data.Cnb    = 0.06;
data.Cnda   = 0.015;
data.Cndr   = -0.05;

% 发动机
data.Tmax  = 45000;       % N
data.xe    = 0.0;         % 推力线相对 c.g. 垂距 m
end

---

四、气动力/力矩拟合（aerodyn_lookup.m）

%% aerodyn_lookup.m
function [CL, CD, Cm] = aerodyn_lookup(alpha_deg, delta_e_deg, Ma, data)
% alpha_deg : 迎角 (deg)
% delta_e   : 升降舵偏角 (deg)

alpha = deg2rad(alpha_deg);
delta_e = deg2rad(delta_e_deg);

% 升力
CL = data.CL0 + data.CLa*alpha + data.CLde*delta_e;

% 阻力（抛物线）
CD = data.CD0 + data.K*(CL^2);

% 俯仰力矩
Cm = data.Cm0 + data.Cma*alpha + data.Cmde*delta_e;

% Ma 修正（简单示例）
if Ma > 0.8
    CL = CL * (1 - 0.1*(Ma-0.8));
    Cm = Cm * (1 + 0.05*(Ma-0.8));
end
end

---

五、配平方程（trim_equations.m）

%% trim_equations.m
function F = trim_equations(x, data, Ma, H, gamma_deg)
% x = [alpha; delta_e]
% 输出残差：俯仰力矩 + 法向力平衡

alpha_deg = x(1);
delta_e   = x(2);

[CL, CD, Cm] = aerodyn_lookup(alpha_deg, delta_e, Ma, data);

% 大气参数
[~, a, ~, rho] = atmosISA(H);
V = Ma * a;

% 气动力
L = 0.5*rho*V^2*data.S*CL;
D = 0.5*rho*V^2*data.S*CD;

% 重力分量
gamma = deg2rad(gamma_deg);
W = data.m * 9.81;

% 法向力平衡（配平爬升角）
Fn = L*cos(gamma) + D*sin(gamma) - W*cos(gamma);

% 俯仰力矩平衡
M = 0.5*rho*V^2*data.S*data.c_bar*Cm;

% 推力力矩（简化）
T = data.Tmax;
Mt = T * data.xe;

F = [M + Mt; Fn];   % 残差向量
end

---

六、主程序：配平求解（main_trim.m）

%% main_trim.m
clear; clc; close all;

data = aircraft_data();

% 飞行条件
Ma = 0.6;
H  = 5000;           % m
gamma_deg = 0;       % 水平飞行

% 初值
x0 = [3; -2];        % [alpha; delta_e]

options = optimoptions('fsolve',...
    'Display','iter',...
    'TolFun',1e-6,...
    'TolX',1e-6);

[x_trim, fval] = fsolve(@(x) trim_equations(x, data, Ma, H, gamma_deg),...
                         x0, options);

alpha_trim = x_trim(1);
delta_e_trim = x_trim(2);

fprintf('=== Trim Results ===\n');
fprintf('Alpha   = %.3f deg\n', alpha_trim);
fprintf('Delta_e = %.3f deg\n', delta_e_trim);
fprintf('Residual = [%.2e, %.2e]\n', fval(1), fval(2));

---

七、配平曲线扫描（plot_trim_curves.m）

%% plot_trim_curves.m
Ma_list = linspace(0.2, 0.9, 20);
alpha_list = zeros(size(Ma_list));
deltae_list = zeros(size(Ma_list));

for i = 1:length(Ma_list)
    x0 = [3; -2];
    [x_trim] = fsolve(@(x) trim_equations(x, data, Ma_list(i), H, 0), x0);
    alpha_list(i) = x_trim(1);
    deltae_list(i) = x_trim(2);
end

figure('Color','white')
subplot(1,2,1)
plot(Ma_list, alpha_list,'bo-','LineWidth',1.5)
xlabel('Mach Number'); ylabel('\alpha_{trim} (deg)')
grid on

subplot(1,2,2)
plot(Ma_list, deltae_list,'ro-','LineWidth',1.5)
xlabel('Mach Number'); ylabel('\delta_e_{trim} (deg)')
grid on
sgtitle('Longitudinal Trim Curves')

---

八、工程级增强

CFD/风洞数据驱动

• 把 aerodyn_lookup.m 换成 二维/三维查表

CL = interp2(alpha_grid, Ma_grid, CL_table, alpha, Ma);

多操纵面等效舵面

• 引入权重因子：

delta_e_eq = w1*delta_e + w2*canard + w3*flap;

非线性优化（fmincon）

• 加入舵偏限制、迎角限制

lb = [ -5; -25];
ub = [ 15;  25];

六自由度配平

• 同时配平：Fx=Fy=Fz=Mx=My=Mz=0
• 变量：α, β, p, q, r, δe, δa, δr, T

参考代码 matlab编写的飞行器等效拟合配平程序 www.youwenfan.com/contentcnw/82504.html

九、适用场景

用途	是否支持
静稳定性分析
控制律初始配平
飞行包线扫描
CFD 数据验证
非线性飞行动力学	（需扩展）