四旋翼无人机动力学数学模型与Matlab仿真
四旋翼无人机动力学数学模型与Matlab仿真
四旋翼无人机是一种常见的无人机类型，其动力学特性对于设计和控制至关重要。本文将从力的来源、数学模型、模型总结和Matlab仿真四个方面详细介绍四旋翼无人机的动力学理论，并通过具体的仿真代码帮助读者理解理论模型的实际应用。
力的来源
无人机的动力系统主要由电调-电机-螺旋桨组成。电机带动螺旋桨旋转产生升力，其原理可以用伯努利定理解释：流速大，压强小；流速小，压强大。螺旋桨的桨面设计使得旋转时桨面上下的空气流速不同，从而产生向上的推力。
数学模型
对于四旋翼无人机，可以将输入简化为四个电机的油门，每个电机的油门归一化到0-1区间。四旋翼动力学数学模型的目标是根据电机的油门计算出升力和各轴的力矩。以下是三个关键模型近似：
- 电机转速与油门的关系：
电机稳态转速与油门成线性关系，实际中存在截距。
电机近似为一阶系统逐渐达到稳态转速。
公式如下：
$$\dot{\bar{\omega}}(t) = \frac{C_m \sigma + \bar{\omega}_m - \bar{\omega}(t)}{T_m}$$
其中，$C_m$为电机转速斜率，$\bar{\omega}_m$为电机转速截距，$T_m$为电机时间常数。
- 螺旋桨升力与转速的关系：
产生的升力与转速平方成正比。
公式如下：
$$T = C_T \sum_{i=1}^{4} \omega_i^2$$
其中，$C_T$为升力系数。
- 合力与力矩的产生：
四个螺旋桨的合力共同作用在机体系z轴。
四个螺旋桨力的差异在机体系三个轴产生力矩。
x和y轴靠升力的不平衡来产生力矩；z轴力矩的产生靠反扭矩。
力矩的计算公式为：
$$\tau = r \times F$$
其中，$\tau$是力矩（单位通常是牛顿米，N·m）；$F$是作用在物体上的力（单位通常是牛顿，N）；$r$是力臂的矢量（单位通常是米，m）。在四旋翼中，力、力臂、力矩的关系可以用下图表示：
其中绿色F1、F2、F3、F4为各电机产生的力，方向为垂直xy平面向上，黄色d为力臂矢量，则力矩Mi的方向通过右手定则可以得到方向，橙色M1、M2、M3、M4则为各电机所产生的力矩。
以电机1为例，其产生的力矩M1在机体系x,y轴的分量为：
$$M_{1x} = d \cdot F_1 \cdot \sin(\theta)$$
$$M_{1y} = d \cdot F_1 \cdot \cos(\theta)$$同理可得到四个电机产生的合力矩，在x和y轴为：
$$\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)$$在z轴方向上，螺旋桨旋转产生的反扭力矩也和螺旋桨的转速平方成正比：
$$\tau_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}_m - \bar{\omega}(t)}{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