15、自由杆冲击动力学的MATLAB模拟与分析

自由杆冲击动力学的MATLAB模拟与分析

1. 引言

在动力学研究中，自由杆的冲击问题是一个经典且重要的研究方向。通过MATLAB进行数值模拟，可以深入了解自由杆在冲击过程中的运动特性和力学行为。本文将详细介绍如何使用MATLAB对自由杆的自由落体和冲击过程进行模拟，并分析其动力学特性。

2. 自由杆的基本参数与初始条件

首先，我们需要定义自由杆的基本参数，包括长度、直径、密度、质量等。以下是相关的MATLAB代码：

L = 0.250; % 杆的长度 (m)
d = 0.012; % 杆的直径 (m)
R = d/2; % 杆的半径 (m)
rho = 7800; % 材料的密度 (kg/m^3)
m = pi*R^2*L*rho; % 杆的质量 (kg)
g = 9.81; % 重力加速度 (m/s^2)
% 杆绕质心C的转动惯量
IC = m*(L^2+3*R^2)/12; % (kg m^2)
H = 0.18; % 下落高度 (m)

这些参数将用于后续的动力学计算。

3. 自由杆的位置与运动描述

自由杆的位置由三个广义坐标(q_1(t))、(q_2(t))和(q_3(t))定义。其中，(q_1)是杆端点(T)的(x)坐标，(q_2)是端点(T)的(y)坐标，(q_3)是杆与(x)轴的夹角。杆的角速度(\omega)和角加速度(\alpha)可以通过对(q_3)求导得到。以下是相关的MATLAB代码：

syms q1(t) q2(t) q3(t)
% q1 x-coordinate of the end T
% q2 y-coordinate of the end T
% q3 angle of the bar with x-axis
omega_ = diff([0 0 q3(t)],t); % 杆的角速度
alpha_ = diff(omega_,t); % 杆的角加速度

杆的接触点(T)、质心(C)和非冲击端(A)的位置向量可以通过广义坐标表示。以下是相关的MATLAB代码：

% T contact point
rT_ = [q1(t) q2(t) 0];
vT_ = diff(rT_,t);
aT_ = diff(vT_,t);
% C mass center
rC_ = rT_+0.5*[L*cos(q3(t)),L*sin(q3(t)),0];
vC_ = simplify(diff(rC_,t));
aC_ = simplify(diff(vC_,t));

4. 自由落体阶段的动力学方程

在自由落体阶段，杆的运动可以通过牛顿 - 欧拉方程描述。以下是相关的MATLAB代码：

G_ = [0 -m*g 0]; % 重力
% Newton Euler equations of motion
Neom =-m*aC_+G_;
Eeom =-IC*alpha_;
eqNE1 = Neom(1);
eqNE2 = Neom(2);
eqNE3 = Eeom(3);

为了求解这些方程，需要将其转化为一阶系统。以下是相关的MATLAB代码：

[V,S] = odeToVectorField(eqNE1, eqNE2, eqNE3);
syms Y
eomNE_ = matlabFunction(V,'vars', {t,Y});

自由落体的初始条件如下：

time0 = [0 1];
% initial conditions
% starts from rest
q10=0; % q1(0)
q20=H; % q2(0)
q30=pi/4; % q3(0) dropping angle
dq10=0; % dq1(0)/dt
dq20=0; % dq2(0)/dt
dq30=0; % dq3(0)/dt
x0 = [q30 dq30 q10 dq10 q20 dq20];

使用 ode45 函数求解非线性系统的数值解，并使用事件函数检测冲击时刻：

option0 = odeset('RelTol',1e-3,'MaxStep',1e-3,'Events',@eventFL0);
[tt0, xs0, t0, y0] = ode45(@eomNE_,time0,x0,option0);
function [value,isterminal,direction] = eventFL0(~,x)
    value = x(5);
    isterminal = 1;
    direction = 0;
end

5. 弹性冲击阶段的动力学方程

当杆的端点(T)接触到水平表面时，进入弹性冲击阶段。此时，法向接触力(P_e)和摩擦力(F_f)可以通过以下公式计算：

E1 = 200e9; % 杨氏模量
nu = 0.29; % 泊松比
k1 = (2/(3*(1-nu^2))*E1*sqrt(R)); % 弹性常数
delta = abs(q2(t));
% 弹性法向接触力
Pe = k1*delta^1.5;
mu = 0.1; % 摩擦系数
P_ = [mu*Pe, Pe, 0]; % 总弹性接触力

弹性冲击阶段的牛顿 - 欧拉方程如下：

% Newton-Euler eom
NeomE =-m*aC_ + G_ + P_;
EeomE =-IC*alpha_ + cross(rT_-rC_, P_);
eqNE1E = NeomE(1);
eqNE2E = NeomE(2);
eqNE3E = EeomE(3);
[V_,S] = odeToVectorField(eqNE1E, eqNE2E, eqNE3E);
eomI_ = matlabFunction(V_,'vars', {t,Y});

弹性冲击的初始条件如下：

q3i = y0(1);
dq3i = y0(2);
q1i = y0(3);
dq1i = y0(4);
q2i = y0(5);
dq2i = y0(6);
xi = [q3i dq3i q1i dq1i q2i dq2i];

使用 ode45 函数求解弹性冲击阶段的数值解：

time0 = [0 1];
[tte, xsE, tE, yE] = ode45(eomI_,time0,xi,option0);

6. 结果分析

在冲击结束时，可以得到以下结果：
|参数|数值|
|----|----|
|(q_1) (m)|0.00|
|(dq_1) (m/s)|-1.99|
|(q_2) (m)|0.00|
|(dq_2) (m/s)|1.88|
|(q_3) (rad)|0.78|
|(dq_3) (rad/s)|-24.38|

可以绘制垂直弹性位移(q_2)、垂直速度(\frac{dq_2}{dt})和弹性力(F_e)的曲线，以及动能(T)的曲线。以下是相关的MATLAB代码：

q1E = yE(3);
dq1E = yE(4);
q2E = yE(5);
dq2E = yE(6);
q3E = yE(1);
dq3E = yE(2);
xE = [q1E dq1E q2E dq2E q3E dq3E];
pe = double(subs(Pe, q2(t), q2E));
q1se = real(xsE(:,3));
dq1se = real(xsE(:,4));
q2se = real(xsE(:,5));
dq2se = real(xsE(:,6));
q3se = real(xsE(:,1));
dq3se = real(xsE(:,2));
q2max=max(abs(q2se));
Fen=subs(Pe,q2(t),q2se);
ql = {diff(q1(t),t), diff(q2(t),t), diff(q3(t),t),...
    q1(t), q2(t), q3(t)};
qne = {dq1se, dq2se, dq3se,...
    q1se, q2se, q3se};
Te = subs(T,ql,qne);
figure(1)
subplot(3,1,1)
plot(tte,q2se,'b','LineWidth',1.5),...
    ylabel('q2 (m)'),grid;
subplot(3,1,2)
plot(tte,dq2se,'b','LineWidth',1.5),...
    ylabel('dq2/dt (m/s)'),grid;
subplot(3,1,3)
plot(tte,Fen,'b','LineWidth',1.5),...
    xlabel('t (s)'),ylabel('Fe (N)'),grid;
figure(2)
plot(tte,Te,'b','LineWidth',1.5),...
    xlabel('t (s)'),ylabel('T (J)')
grid on

7. 拉格朗日方法

除了牛顿 - 欧拉方法，还可以使用拉格朗日方法求解自由杆的动力学方程。拉格朗日方程如下：
(\frac{d}{dt}(\frac{\partial T}{\partial \dot{q}_i}) - \frac{\partial T}{\partial q_i} = Q_i)，(i = 1, 2, 3)
其中，(T)是系统的动能，(Q_i)是广义力。以下是使用拉格朗日方法的MATLAB代码：

% 杆在一般平面运动中的动能
T = 0.5*( m*(vC_*vC_.')+IC*(omega_*omega_.') );
% dT/d(dq/dt)
Tdq1 = diff(T, diff(q1(t),t)); % dT/dq1'
Tdq2 = diff(T, diff(q2(t),t)); % dT/dq2'
Tdq3 = diff(T, diff(q3(t),t)); % dT/dq3'
% d(dT/d(dq/dt))/dt
Tt1 = diff(Tdq1, t); % d(dT/dq1')/dt
Tt2 = diff(Tdq2, t); % d(dT/dq2')/dt
Tt3 = diff(Tdq3, t); % d(dT/dq3')/dt
% dT/dq
Tq1 = diff(T, q1(t)); % dT/dq1
Tq2 = diff(T, q2(t)); % dT/dq2
Tq3 = diff(T, q3(t)); % dT/dq3
% 拉格朗日方程的左边
LHS1 = Tt1 - Tq1; % d(dT/dq1')/dt - dT/dq1
LHS2 = Tt2 - Tq2; % d(dT/dq2')/dt - dT/dq2
LHS3 = Tt3 - Tq3; % d(dT/dq3')/dt - dT/dq3
% 偏导数
rC_1_ = diff(rC_, q1(t)); % drC_/dq1
rC_2_ = diff(rC_, q2(t)); % drC_/dq2
rC_3_ = diff(rC_, q3(t)); % drC_/dq3
rT_1_ = diff(rT_, q1(t)); % drT_/dq1
rT_2_ = diff(rT_, q2(t)); % drT_/dq2
rT_3_ = diff(rT_, q3(t)); % drT_/dq3
% 广义主动力Q
G_ = [0 -m*g 0]; % 重力
% 对应于q1的广义力
Q1 = rC_1_*G_.'; % Q1 = G_.drC_/dq1
% 对应于q2的广义力
Q2 = rC_2_*G_.'; % Q2 = G_.drC_/dq2
% 对应于q3的广义力
Q3 = rC_3_*G_.'; % Q3 = G_.drC_/dq3
% 拉格朗日方程
Lagrange1 = LHS1 - Q1;
Lagrange2 = LHS2 - Q2;
Lagrange3 = LHS3 - Q3;
[V,S] = odeToVectorField(Lagrange1,Lagrange2,Lagrange3);
syms Y
eomLagra_ = matlabFunction(V,'vars', {t,Y});

求解自由落体和弹性冲击阶段的数值解的步骤与牛顿 - 欧拉方法类似。

8. 问题求解

以下是几个相关的问题及其求解思路：
- 问题9.1 ：一个质量为(m)、长度为(L)的细长旋转杆从(\theta(0) = \frac{\pi}{4})的位置静止释放。当角度(\theta)为(\theta_0)时，杆的端点(A)接触到弹簧和阻尼器系统。求杆的最大力以及从杆释放到与弹簧和阻尼器系统接触结束的时间。
- 已知参数：

L = 0.5; % 杆的长度 (m)
m = 2; % 质量 (kg)
g = 9.81; % 重力加速度 (m/s^2)
k = 1500; % 弹簧常数 (N/m)
c = 20; % 粘性阻尼系数 (N s/m)
theta_0 = -pi/4; % 接触角度

- 求解思路：建立杆的动力学方程，考虑弹簧和阻尼器的作用，使用数值方法求解。

• 问题9.2 ：一个质量为(m)、长度为(L)的细长杆，在销钉支撑(O)处有一个与角速度(\frac{d\theta}{dt})相反的阻力矩(b\frac{d\theta}{dt})。杆从初始角度(\theta_0)静止开始摆动。有两组初始条件（小角度和大角度）和两个特定时刻(t_1)和(t_2)。
  • 已知参数：

m = 4; % 质量 (kg)
L = 2; % 长度 (m)
b = 1.4;
g = 9.81; % 重力加速度 (m/s^2)
theta01 = 2*pi/180; % 小角度 (rad)
theta02 = 75*pi/180; % 大角度 (rad)
t1 = 2; % 特定时刻 (s)
t2 = 10; % 特定时刻 (s)

- 求解思路：分别建立线性化和非线性的动力学方程，使用数值方法求解不同初始条件下的响应\(\theta(t)\)，并绘制相应的曲线。

9. 总结

本文详细介绍了如何使用MATLAB对自由杆的自由落体和弹性冲击过程进行模拟和分析。通过牛顿 - 欧拉方法和拉格朗日方法建立了动力学方程，并使用数值方法求解。同时，还给出了几个相关问题的求解思路。这些方法和代码可以帮助读者更好地理解自由杆的动力学特性，为实际工程应用提供参考。

以下是整个过程的mermaid流程图：

graph TD;
    A[定义基本参数] --> B[描述自由杆位置与运动];
    B --> C[求解自由落体阶段动力学方程];
    C --> D[检测冲击时刻];
    D --> E[求解弹性冲击阶段动力学方程];
    E --> F[结果分析];
    F --> G[拉格朗日方法求解];
    G --> H[问题求解];

自由杆冲击动力学的MATLAB模拟与分析

10. 详细操作步骤总结

为了更清晰地展示整个模拟和分析过程，下面将上述内容总结为详细的操作步骤：
1. 定义基本参数
- 输入杆的长度、直径、密度、重力加速度、下落高度等参数。
- 计算杆的质量、转动惯量等。
matlab L = 0.250; % 杆的长度 (m) d = 0.012; % 杆的直径 (m) R = d/2; % 杆的半径 (m) rho = 7800; % 材料的密度 (kg/m^3) m = pi*R^2*L*rho; % 杆的质量 (kg) g = 9.81; % 重力加速度 (m/s^2) IC = m*(L^2+3*R^2)/12; % 杆绕质心C的转动惯量 (kg m^2) H = 0.18; % 下落高度 (m)
2. 描述自由杆位置与运动
- 定义广义坐标(q_1(t))、(q_2(t))和(q_3(t))。
- 计算杆的角速度、角加速度。
- 计算接触点(T)、质心(C)和非冲击端(A)的位置向量、速度向量和加速度向量。
matlab syms q1(t) q2(t) q3(t) omega_ = diff([0 0 q3(t)],t); % 杆的角速度 alpha_ = diff(omega_,t); % 杆的角加速度 rT_ = [q1(t) q2(t) 0]; vT_ = diff(rT_,t); aT_ = diff(vT_,t); rC_ = rT_+0.5*[L*cos(q3(t)),L*sin(q3(t)),0]; vC_ = simplify(diff(rC_,t)); aC_ = simplify(diff(vC_,t));
3. 求解自由落体阶段动力学方程
- 建立牛顿 - 欧拉方程。
- 将方程转化为一阶系统。
- 定义初始条件。
- 使用 ode45 函数求解数值解。
matlab G_ = [0 -m*g 0]; % 重力 Neom =-m*aC_+G_; Eeom =-IC*alpha_; eqNE1 = Neom(1); eqNE2 = Neom(2); eqNE3 = Eeom(3); [V,S] = odeToVectorField(eqNE1, eqNE2, eqNE3); syms Y eomNE_ = matlabFunction(V,'vars', {t,Y}); time0 = [0 1]; q10=0; % q1(0) q20=H; % q2(0) q30=pi/4; % q3(0) dropping angle dq10=0; % dq1(0)/dt dq20=0; % dq2(0)/dt dq30=0; % dq3(0)/dt x0 = [q30 dq30 q10 dq10 q20 dq20]; option0 = odeset('RelTol',1e-3,'MaxStep',1e-3,'Events',@eventFL0); [tt0, xs0, t0, y0] = ode45(@eomNE_,time0,x0,option0); function [value,isterminal,direction] = eventFL0(~,x) value = x(5); isterminal = 1; direction = 0; end
4. 检测冲击时刻
- 使用事件函数检测杆端点(T)接触到水平表面的时刻。
5. 求解弹性冲击阶段动力学方程
- 计算法向接触力和摩擦力。
- 建立弹性冲击阶段的牛顿 - 欧拉方程。
- 将方程转化为一阶系统。
- 定义初始条件。
- 使用 ode45 函数求解数值解。
matlab E1 = 200e9; % 杨氏模量 nu = 0.29; % 泊松比 k1 = (2/(3*(1-nu^2))*E1*sqrt(R)); % 弹性常数 delta = abs(q2(t)); Pe = k1*delta^1.5; mu = 0.1; % 摩擦系数 P_ = [mu*Pe, Pe, 0]; % 总弹性接触力 NeomE =-m*aC_ + G_ + P_; EeomE =-IC*alpha_ + cross(rT_-rC_, P_); eqNE1E = NeomE(1); eqNE2E = NeomE(2); eqNE3E = EeomE(3); [V_,S] = odeToVectorField(eqNE1E, eqNE2E, eqNE3E); eomI_ = matlabFunction(V_,'vars', {t,Y}); q3i = y0(1); dq3i = y0(2); q1i = y0(3); dq1i = y0(4); q2i = y0(5); dq2i = y0(6); xi = [q3i dq3i q1i dq1i q2i dq2i]; time0 = [0 1]; [tte, xsE, tE, yE] = ode45(eomI_,time0,xi,option0);
6. 结果分析
- 提取冲击结束时的参数。
- 绘制垂直弹性位移、垂直速度、弹性力和动能的曲线。
matlab q1E = yE(3); dq1E = yE(4); q2E = yE(5); dq2E = yE(6); q3E = yE(1); dq3E = yE(2); xE = [q1E dq1E q2E dq2E q3E dq3E]; pe = double(subs(Pe, q2(t), q2E)); q1se = real(xsE(:,3)); dq1se = real(xsE(:,4)); q2se = real(xsE(:,5)); dq2se = real(xsE(:,6)); q3se = real(xsE(:,1)); dq3se = real(xsE(:,2)); q2max=max(abs(q2se)); Fen=subs(Pe,q2(t),q2se); ql = {diff(q1(t),t), diff(q2(t),t), diff(q3(t),t),... q1(t), q2(t), q3(t)}; qne = {dq1se, dq2se, dq3se,... q1se, q2se, q3se}; Te = subs(T,ql,qne); figure(1) subplot(3,1,1) plot(tte,q2se,'b','LineWidth',1.5),... ylabel('q2 (m)'),grid; subplot(3,1,2) plot(tte,dq2se,'b','LineWidth',1.5),... ylabel('dq2/dt (m/s)'),grid; subplot(3,1,3) plot(tte,Fen,'b','LineWidth',1.5),... xlabel('t (s)'),ylabel('Fe (N)'),grid; figure(2) plot(tte,Te,'b','LineWidth',1.5),... xlabel('t (s)'),ylabel('T (J)') grid on
7. 拉格朗日方法求解
- 计算系统的动能。
- 计算拉格朗日方程的各项。
- 建立拉格朗日方程。
- 将方程转化为一阶系统。
- 求解自由落体和弹性冲击阶段的数值解。
matlab T = 0.5*( m*(vC_*vC_.')+IC*(omega_*omega_.') ); Tdq1 = diff(T, diff(q1(t),t)); % dT/dq1' Tdq2 = diff(T, diff(q2(t),t)); % dT/dq2' Tdq3 = diff(T, diff(q3(t),t)); % dT/dq3' Tt1 = diff(Tdq1, t); % d(dT/dq1')/dt Tt2 = diff(Tdq2, t); % d(dT/dq2')/dt Tt3 = diff(Tdq3, t); % d(dT/dq3')/dt Tq1 = diff(T, q1(t)); % dT/dq1 Tq2 = diff(T, q2(t)); % dT/dq2 Tq3 = diff(T, q3(t)); % dT/dq3 LHS1 = Tt1 - Tq1; % d(dT/dq1')/dt - dT/dq1 LHS2 = Tt2 - Tq2; % d(dT/dq2')/dt - dT/dq2 LHS3 = Tt3 - Tq3; % d(dT/dq3')/dt - dT/dq3 rC_1_ = diff(rC_, q1(t)); % drC_/dq1 rC_2_ = diff(rC_, q2(t)); % drC_/dq2 rC_3_ = diff(rC_, q3(t)); % drC_/dq3 rT_1_ = diff(rT_, q1(t)); % drT_/dq1 rT_2_ = diff(rT_, q2(t)); % drT_/dq2 rT_3_ = diff(rT_, q3(t)); % drT_/dq3 G_ = [0 -m*g 0]; % 重力 Q1 = rC_1_*G_.'; % Q1 = G_.drC_/dq1 Q2 = rC_2_*G_.'; % Q2 = G_.drC_/dq2 Q3 = rC_3_*G_.'; % Q3 = G_.drC_/dq3 Lagrange1 = LHS1 - Q1; Lagrange2 = LHS2 - Q2; Lagrange3 = LHS3 - Q3; [V,S] = odeToVectorField(Lagrange1,Lagrange2,Lagrange3); syms Y eomLagra_ = matlabFunction(V,'vars', {t,Y});
8. 问题求解
- 针对不同问题，建立相应的动力学方程，使用数值方法求解。

11. 方法对比

牛顿 - 欧拉方法和拉格朗日方法是求解动力学问题的两种常用方法，它们各有优缺点，下面通过表格进行对比：
| 方法 | 优点 | 缺点 |
| ---- | ---- | ---- |
| 牛顿 - 欧拉方法 | - 物理意义明确，直观地反映了力和加速度的关系。
- 对于简单系统，方程建立相对容易。 | - 对于复杂系统，需要考虑较多的力和力矩，方程推导过程繁琐。
- 约束处理相对复杂。 |
| 拉格朗日方法 | - 只需要考虑系统的动能和广义力，不需要直接分析力和力矩，方程推导相对简单。
- 对于有约束的系统，处理起来较为方便。 | - 物理意义不如牛顿 - 欧拉方法直观。
- 计算过程中可能涉及较多的符号运算，计算量较大。 |

在实际应用中，可以根据系统的复杂程度和具体需求选择合适的方法。如果系统结构简单，力和力矩的分析比较直观，牛顿 - 欧拉方法可能更合适；如果系统存在较多的约束，拉格朗日方法可能更具优势。

12. 注意事项

在使用MATLAB进行模拟和分析时，需要注意以下几点：
1. 符号运算 ：Lagrange方法中涉及大量的符号运算，可能会导致计算速度较慢。可以在适当的时候进行化简操作，提高计算效率。
2. 初始条件 ：初始条件的设置对模拟结果有很大影响，需要根据实际情况准确设置。
3. 数值求解 ：使用 ode45 等数值求解函数时，需要合理设置求解选项，如相对误差容限 RelTol 和最大步长 MaxStep ，以保证求解的准确性和稳定性。
4. 事件函数 ：事件函数用于检测特定事件的发生，如冲击时刻的检测。需要确保事件函数的定义正确，以准确捕捉事件。

13. 拓展应用

自由杆冲击动力学的研究在很多领域都有重要的应用，例如：
1. 机械工程 ：在机械设计中，了解杆类零件的冲击特性可以帮助优化设计，提高机械系统的可靠性和安全性。
2. 航空航天 ：飞行器在起飞、着陆等过程中可能会受到冲击，研究自由杆的冲击动力学可以为飞行器的结构设计提供参考。
3. 材料科学 ：通过研究杆与材料的冲击相互作用，可以评估材料的抗冲击性能，为材料的选择和开发提供依据。

14. 总结与展望

本文详细介绍了使用MATLAB对自由杆的自由落体和弹性冲击过程进行模拟和分析的方法，包括牛顿 - 欧拉方法和拉格朗日方法。通过具体的代码示例和操作步骤，展示了如何建立动力学方程、求解数值解和进行结果分析。同时，还对两种方法进行了对比，指出了各自的优缺点和适用场景。

在未来的研究中，可以进一步拓展该模型，考虑更多的因素，如杆的非线性特性、材料的阻尼特性等。此外，还可以将该模型应用到更复杂的系统中，如多杆机构的冲击问题，以更好地满足实际工程的需求。

以下是整个过程的另一个mermaid流程图，重点突出不同阶段的输入输出：

graph LR;
    classDef startend fill:#F5EBFF,stroke:#BE8FED,stroke-width:2px;
    classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;
    classDef decision fill:#FFF6CC,stroke:#FFBC52,stroke-width:2px;

    A([定义基本参数]):::startend --> B(描述自由杆位置与运动):::process;
    B --> C(求解自由落体阶段动力学方程):::process;
    C --> D{是否检测到冲击?}:::decision;
    D -->|否| C;
    D -->|是| E(求解弹性冲击阶段动力学方程):::process;
    E --> F(结果分析):::process;
    F --> G(拉格朗日方法求解):::process;
    G --> H(问题求解):::process;
    H --> I([结束]):::startend;

通过以上内容，希望读者能够掌握使用MATLAB进行自由杆冲击动力学模拟和分析的方法，为相关领域的研究和应用提供帮助。