四旋翼无人机动力学数学模型与Matlab仿真
四旋翼无人机动力学数学模型与Matlab仿真
四旋翼无人机作为一种常见的无人机类型，其动力学特性是理解和控制无人机飞行的基础。本文将从力的来源、数学模型、模型总结和Matlab仿真四个方面，详细介绍四旋翼无人机的动力学原理。
力的来源
无人机的动力系统主要由电调-电机-螺旋桨组成。直观来看，电机带动螺旋桨旋转产生升力。螺旋桨桨面设计特殊，旋转时桨面上下空气流速不同，根据伯努利定理（流速大压强小，流速小压强大），会产生向上的推力。
数学模型
对于四旋翼无人机，可以将输入简化为四个电机的油门，每个电机的油门归一化到0-1区间。四旋翼动力学数学模型的目标是根据电机的油门计算出升力和各轴的力矩。
电机转速与油门关系
电机转速与油门的关系可以近似为一阶系统，公式如下：
$$
\dot{\bar{\omega}}(t) = \frac{C_m \sigma - \bar{\omega}(t) + \bar{\omega}_m}{T_m}
$$
其中：
- $C_m$：电机转速斜率
- $\bar{\omega}_m$：电机转速截距
- $T_m$：电机时间常数
升力与转速关系
升力与转速平方成正比，公式如下：
$$
T = C_T \sum_{i=1}^{4} \omega_i^2
$$
其中：
- $C_T$：升力系数
- $T$：无人机的合升力
力矩计算
力矩是描述力对物体产生旋转效果的物理量，公式为 $\tau = r \times F$。在四旋翼中，力、力臂、力矩的关系如下图所示：
其中：
- $F_1, F_2, F_3, F_4$：各电机产生的力
- $d$：力臂矢量
- $M_1, M_2, M_3, M_4$：各电机所产生的力矩
在z轴方向上，螺旋桨旋转会产生反扭矩，其大小与转速平方成正比：
$$
M_z = C_M (\omega_1^2 - \omega_2^2 + \omega_3^2 - \omega_4^2)
$$
其中：
- $C_M$：反扭力矩系数
数学模型总结
油门和电机转速的计算公式：
$$
\dot{\bar{\omega}}(t) = \frac{C_m \sigma - \bar{\omega}(t) + \bar{\omega}_m}{T_m}
$$转速和升力的计算公式：
$$
T = C_T \sum_{i=1}^{4} \omega_i^2
$$转速和力矩的计算公式：
$$
\tau_x = \frac{\sqrt{2}}{2} d C_T (-\omega_1^2 + \omega_2^2 + \omega_3^2 - \omega_4^2) \
\tau_y = \frac{\sqrt{2}}{2} d C_T (\omega_1^2 + \omega_2^2 - \omega_3^2 - \omega_4^2) \
\tau_z = C_M (\omega_1^2 - \omega_2^2 + \omega_3^2 - \omega_4^2)
$$
Matlab 仿真
油门与电机转速模型
%% 油门与电机转速模型测试
global dt Tm Cm varpim
dt = 1e-3; % 仿真时间步长
Cm = 706.01; % 油门增大1，电机转速变化(RPM)
varpim = 170.47; % 零占空比时电机转速(RPM)
Tm = 0.260; % 电机时间常数
N = 2000;
t = 0:dt:dt*(N-1);
sigma = [0.7; 0.6; 0.5; 0.4];
varpi = zeros(N, 4);
k=1;
for tt=0:dt:(N-2)*dt
k = k+1;
% 动力单元模型
varpi(k, 1) = motor(sigma(1), varpi(k-1, 1)); % 电机1转速
varpi(k, 2) = motor(sigma(2), varpi(k-1, 2)); % 电机2转速
varpi(k, 3) = motor(sigma(3), varpi(k-1, 3)); % 电机3转速
varpi(k, 4) = motor(sigma(4), varpi(k-1, 4)); % 电机4转速
end
figure(1);plot(t, varpi(:,1), 'LineWidth', 1.5); hold on
plot(t, varpi(:,2), 'LineWidth', 1.5);
plot(t, varpi(:,3), 'LineWidth', 1.5);
plot(t, varpi(:,4), 'LineWidth', 1.5); hold off
legend(['\sigma_1=' num2str(sigma(1))], ['\sigma_2=' num2str(sigma(2))],['\sigma_3=' num2str(sigma(3))],['\sigma_4=' num2str(sigma(4))]);
xlabel('时间 t (s)');ylabel('转速 \varpi (rad/s)');title('电机模型测试'); grid on; grid minor
%% 电机模型
% 输入：油门大小 sigma（0-1）
% 电机上一时刻的转速（rad/s）
% 输出：此时刻电机转速（rad/s）
function varpi = motor(sigma, varpi_)
global dt Tm Cm varpim;
dvarpi = (Cm * sigma + varpim - varpi_) / Tm * dt;
varpi = varpi_ + dvarpi;
end
油门与升力、力矩的关系仿真
%% 油门与升力、力矩模型测试
global dt Tm Cm varpim d cT cM
dt = 1e-3; % 仿真时间步长
Cm = 706.01; % 油门增大1，电机转速变化(RPM)
varpim = 170.47; % 零占空比时电机转速(RPM)
Tm = 0.260; % 电机时间常数
d = 0.225; % 450mm/2
cT = 1.201e-5; % 升力系数
cM = 1.574e-7; % 反扭力系数
N = 2000;
t = 0:dt:dt*(N-1);
sigma = [0.7; 0.6; 0.5; 0.4];
varpi = zeros(N, 4);
T = zeros(N, 1);
tau = zeros(N, 3);
k=1;
for tt=0:dt:(N-2)*dt
k = k+1;
% 电机模型
varpi(k, 1) = motor(sigma(1), varpi(k-1, 1)); % 电机1转速
varpi(k, 2) = motor(sigma(2), varpi(k-1, 2)); % 电机2转速
varpi(k, 3) = motor(sigma(3), varpi(k-1, 3)); % 电机3转速
varpi(k, 4) = motor(sigma(4), varpi(k-1, 4)); % 电机4转速
[T(k), tau(k,:)] = power_mix(varpi(k, :));
end
figure(1);subplot(211); plot(t, T, 'linewidth', 1.5); title('动力合成模型');ylabel('升力 (N)');
subplot(212);plot(t, tau(:,1), 'linewidth', 1.5);hold on
plot(t, tau(:,2),'linewidth', 1.5);plot(t, tau(:,3),'linewidth', 1.5);hold off
ylabel('力矩 (N\cdotm)');xlabel('时间 (t)'); legend('\tau_x', '\tau_y', '\tau_z');
%% 电机模型
% 输入：油门大小 sigma（0-1）
% 电机上一时刻的转速（rad/s）
% 输出：此时刻电机转速（rad/s）
function varpi = motor(sigma, varpi_)
global dt Tm Cm varpim;
dvarpi = (Cm * sigma + varpim - varpi_) / Tm * dt;
varpi = varpi_ + dvarpi;
end
%% 动力合成模型
% 输入：四个电机转速
% 输出：合升力与三轴力矩
function [T, tau] = power_mix(varpi)
global cT cM d;
T = cT * sum(varpi.^2);
tau(1) = sqrt(2)/2 * d * cT * (-varpi(1)^2 + varpi(2)^2 + varpi(3)^2 - varpi(4)^2);
tau(2) = sqrt(2)/2 * d * cT * ( varpi(1)^2 + varpi(2)^2 - varpi(3)^2 - varpi(4)^2);
tau(3) = cM * (varpi(1)^2 - varpi(2)^2 + varpi(3)^2 - varpi(4)^2);
end