14、平面运动物体的直接动力学分析与MATLAB实现

平面运动物体的直接动力学分析与MATLAB实现

在机械系统的动力学分析中，我们常常需要研究物体在不同力作用下的运动状态。本文将借助MATLAB工具，深入探讨两个典型的动力学问题：球体与弹簧的相互作用，以及旋转连杆在弹性力作用下的运动。

1. 球体与弹簧的动力学分析

当一个球体落在垂直弹簧上时，会经历压缩和恢复两个阶段。下面我们将详细分析这个过程，并通过MATLAB进行数值模拟。

1.1 运动方程的建立

假设在 $t = 0$ 时刻，球体与弹簧接触，此时球体的速度为 $v_0$。球体与弹簧接触时的运动方程为：
[m\ddot{x} = mg - kx]
其中，$m$ 是球体的质量，$g$ 是重力加速度，$k$ 是弹簧的劲度系数，$x$ 是球体的位移。

接触力由弹性力产生，表达式为：
[P = -kx\mathbf{i}]
初始条件为：
[x(0) = 0, \quad \dot{x}(0) = v_0]

1.2 MATLAB求解符号解

使用MATLAB的 dsolve 函数求解符号解，代码如下：

syms x(t)
syms m k g h xo vo
assume(m,'positive') % 小球体的质量
assume(k,'positive') % 弹簧常数
assume(g,'positive') % 重力加速度
assume(xo,'positive') % 位移的初始条件
assume(vo,'positive') % 速度的初始条件
% 运动方程
eqn = m*diff(x,2) + k*x - m*g;
v = diff(x);
xSol = dsolve(eqn, x(0)==xo, v(0)==vo);
vSol = diff(xSol,t);
fprintf('x(t)=\n')
pretty(xSol)
fprintf('v(t)=\n')
pretty(vSol)
fprintf('\n')

1.3 数值输入与求解

给定数值输入：

lists = { m, g, k, h}; % 符号列表
listn = {10, 9.81, 294e3, 1}; % 数值列表
% m = 10 (kg) 质量
% g = 9.81 (m/s^2) 重力加速度
% k = 294e3 (N/m) 弹簧常数
% h = 1 (m) 下落高度
% v0 = sqrt(2*g*h) 初始速度 (m/s)
x0 = 0; % (m) x(0) = 0
v0 = subs(sqrt(2*g*h), lists, listn); % (m/s)
xn = subs(xSol, lists, listn);
xn = subs(xn, {xo,vo}, {x0, v0});
xn = vpa(xn,3);
vn = subs(vSol, lists, listn);
vn = subs(vn, {xo,vo}, {x0, v0});
vn = vpa(vn,3);
fprintf('x(t)=\n')
pretty(xn)
fprintf('v(t)=\n')
pretty(vn)
fprintf('\n')

1.4 求解关键时间点

求解球体达到最大位移的时间 $t_v$ 和分离时间 $t_f$：

t0 = 0; % 初始时间
t2 = 0.02; % 任意选择的最终时间
tv = vpasolve(vn, t, [t0 t2]);
tv = double(tv);
fprintf('tv = %6.4f (s)\n',tv)
tf = vpasolve(xn, t, [tv t2]);
tf = double(tf);
fprintf('tf = %6.4f (s)\n',tf)

1.5 计算接触弹性力

接触弹性力函数为：

k = 294e3; % 弹簧常数 (N/m)
Pen = k*xn; % (N) 接触弹性力
xmax = subs(xn, t, tv); % tv 时刻的最大位移
Pmax = subs(Pen,t, tv); % tv 时刻的最大接触弹性力

1.6 绘制结果

使用以下代码绘制位移、速度和弹性力随时间的变化曲线：

figure
subplot(3,1,1)
fplot(xn, [0 tf],'-k','Linewidth',2)
ylabel('x (m)')
xlabel('time t (s)')
grid on
subplot(3,1,2)
fplot(vn, [0 tf],'-k','Linewidth',2)
ylabel('dx/dt (m/s)')
xlabel('time t (s)')
grid on
subplot(3,1,3)
fplot(Pen, [0 tf],'-k','Linewidth',2)
ylabel('Pe (N)'), xlabel('t (s)')
grid on

1.7 数值求解运动方程

将二阶微分方程改写为一阶系统：
[\dot{x}_1 = x_2]
[\dot{x}_2 = g - \frac{k}{m}x_1]
MATLAB代码如下：

syms t x1 x2
dx1 = x2;
dx2 = g-k*x1/m;
dx = [dx1; dx2];
g = matlabFunction(dx,'vars',{t, [x1; x2]});

创建 eventx1.m 函数来停止积分：

function [value,isterminal,direction] = eventx1(t,x)
value = x(1);
isterminal = 1;
direction = 0;
end

设置ODE求解器的选项：

optionx1 = odeset('RelTol',1e-6,'MaxStep',1e-6,...
    'Events',@eventx1);
ic0 = [0 v0];
t0 = 0; tf = 1; ti = [t0 tf];
[tte, xse, te, ye] = ode45(g, ti, ic0, optionx1);

获取时间向量和解决方案：

time = tte;
dis = xse(:,1);
vel = xse(:,2);
Pen = k*dis;
xmax = max(dis);
Pmax = max(Pen);

绘制数值解的曲线：

figure(1)
grid on
subplot(3,1,1)
plot(time,dis,'LineWidth',1.5)
ylabel('x (m)')
grid on
subplot(3,1,2)
xlabel('t (s)')
plot(time,vel,'LineWidth',1.5)
ylabel('dx/dt (m/s)')
grid on
subplot(3,1,3)
plot(time,Pen,'LineWidth',1.5)
ylabel('Pe (N)'), xlabel('t (s)')
grid on

2. 旋转连杆在弹性力作用下的动力学分析

考虑一个均匀圆柱形复合摆，在特定角度时受到垂直弹性力的作用。下面我们将分阶段分析其运动，并使用MATLAB进行模拟。

2.1 问题描述与参数设置

摆杆长度 $L = 0.50$ m，直径 $d = 0.01$ m，密度 $\rho = 7800$ kg/m³。弹性力 $F_e$ 在 $\theta(t_0) = \theta_0 = -\frac{\pi}{6}$ 时施加，弹性常数 $k = 2000$，$\alpha = 1.5$。

2.2 运动学分析

计算质心和端点的位置向量：

L = 0.50; % 连杆长度 (m)
d = 0.01; % 直径 (m)
R = d/2; % 半径 (m)
rho = 7800; % 密度 (kg/m^3)
m = pi*R^2*L*rho; % 质量 (kg)
g = 9.81; % 重力加速度 (m/s^2)
alpha = 1.5;
k = 2000; % 弹簧常数
syms theta(t)
omega_ = [0 0 diff(theta(t),t)];
alpha_ = diff(omega_,t);
c = cos(theta);
s = sin(theta);
xC = (L/2)*c;
yC = (L/2)*s;
rC_ = [xC yC 0];
xA = L*c;
yA = L*s;
rA_ = [xA yA 0];
vA_ = diff(rA_,t);

计算绕固定点 $O$ 的转动惯量：

IC = m*(L^2+3*R^2)/12; % 绕质心C的转动惯量
IO = IC + m*(L/2)^2; % 绕固定点O的转动惯量

2.3 分阶段分析

• 阶段一：弹性力施加前
  重力驱动连杆运动，运动方程为：
  [I_O\alpha = \sum M_O = \mathbf{r}_C \times \mathbf{G}]
  其中，$\mathbf{G} = -mg\mathbf{j}$。
  MATLAB代码如下：

G_ = [0 -m*g 0]; % 质心处的重力
MO_ = cross(rC_,G_); % 绕O点的力矩
eq_ = -IO*alpha_ + MO_;
eqz = eq_(3);
[V,S] = odeToVectorField(eqz);
syms Y
eom0_ = matlabFunction(V,'vars', {t,Y});

创建 event1.m 函数来检测特定角度：

function[value,ist,dir]=event1(~,x)
value = x(1)+pi/6;
ist = 1;
dir = 0;
end

求解运动方程：

time0_ = [0 1]; % 时间域
x0_ = [0 0]; % 初始条件
option0 = ...
    odeset('RelTol',1e-3,'MaxStep',1e-3,'Events',@event1);
[tt0, xs0, t0, ye0] = ...
    ode45(eom0_,time0_,x0_,option0);
theta0 = ye0(1);
omega0 = ye0(2);
fprintf('initial conditions for contact\n')
fprintf('theta0 = %6.3g (rad) = %g (degrees)\n',...
    theta0, theta0*180/pi)
fprintf('omega0 = %6.3g (rad/s)\n', omega0)
qt = {diff(theta,t),theta};
list0 = {omega0, theta0};
vAn_ = subs(vA_, qt, list0);
fprintf('vA_ = [%6.3g, %6.3g, %g] (m/s)\n', vAn_)

• 阶段二：弹性力施加后
  弹性力为：
  [F_e = k|\delta|^\alpha\mathbf{j}]
  其中，$\delta = L\sin\theta - L\sin\theta_0$。
  运动方程为：
  [I_O\alpha = \sum M_O = \mathbf{r}_C \times \mathbf{G} + \mathbf{r}_A \times \mathbf{F}_e]
  MATLAB代码如下：

delta = L*sin(theta(t))-L*sin(theta0);
Fe = k*abs(delta)^alpha; % 弹性力
Fe_ = [0, Fe, 0];
MOe_ = cross(rC_,G_) + cross(rA_,Fe_);
eqee_ = -IO*alpha_ + MOe_;
eqze = eqee_(3);
[V,S] = odeToVectorField(eqze);
eome_ = matlabFunction(V,'vars', {t,Y});

求解运动方程：

timee_ = [0 1];
[tte, xse, te, yee] = ...
    ode45(eome_,timee_, [theta0 omega0], option0);
thetaf = yee(1);
omegaf = yee(2);
fprintf('final conditions for contact\n')
fprintf('tf = %6.3g (s) \n',te);
fprintf('thetaf = %6.3g (rad) = %g (degrees)\n',...
    thetaf, thetaf*180/pi)
fprintf('omegaf = %6.3g (rad/s)\n', omegaf)
thetan = xse(:,1);
omegan = xse(:,2);
deltan = L*(sin(thetan) - sin(theta0)*ones(length(thetan),1));
Fen = k*abs(deltan).^alpha;
thetam = min(thetan);
Fm = max(Fen);
fprintf('thetam = %6.3g (rad) = %6.3g (degrees)\n',...
    thetam, thetam*180/pi)
fprintf('Fm = %6.3g (N)\n', Fm)

2.4 动画展示

使用以下代码展示连杆在弹性力作用下的运动：

for i = 1:10:length(thetan)
    clf
    xAn = L*cos(thetan(i));
    yAn = L*sin(thetan(i));
    Ov = zeros(length(thetan(i)),1);
    xO = Ov; yO = Ov;
    xlabel('x(m)'), ylabel('y(m)')
    sf=0.5;
    axis([-sf sf -sf sf])
    grid on
    hold on
    % 绘制关节O
    plot_pin0(0,0,pin_scale_factor,0)
    plot([xO xAn],[yO yAn],'k','LineWidth',1.5)
    % 标记A点
    textAx = xAn+0.01;
    textAy = yAn;
    text(textAx, textAy,'A','FontSize',text_size)
    % 绘制弹性力向量
    quiver(xAn,yAn,0,.15,0,'Color','r','LineWidth',2.5)
    text(xAn+0.01,yAn+.07,'F_e',...
        'fontsize',14,'fontweight','b')
    pause(0.2);
end

3. 拉格朗日方程求解

使用拉格朗日方程也可以得到相同的结果。拉格朗日方程为：
[\frac{d}{dt}\left(\frac{\partial T}{\partial \dot{\theta}}\right) - \frac{\partial T}{\partial \theta} = Q]
其中，$T$ 是系统的动能，$Q$ 是广义主动力。

3.1 计算动能和广义主动力

T = 0.5*IO*(omega_*omega_.');
Tdq = diff(T, diff(theta(t),t));
Tt = diff(Tdq, t);
Tq = diff(T, theta(t));
LHS = Tt - Tq;
drC_ = diff(rC_,theta(t));
G_ = [0 -m*g 0];
Q = drC_*G_.';
Lagrange = LHS - Q;
[V,S] = odeToVectorField(Lagrange);
syms Y
eom0_ = matlabFunction(V,'vars', {t,Y});

3.2 求解运动方程

time0 = [0 1];
x0 = [0 0]; % 初始条件
option0 = ...
    odeset('RelTol',1e-3,'MaxStep',1e-3,'Events',@event1);
[tt0, xs0, t0, ye0] = ...
    ode45(eom0_,time0,x0,option0);
theta0 = ye0(1);
omega0 = ye0(2);
fprintf('initial conditions for contact\n')
fprintf('theta0 = %6.3g (rad) = %g (degrees)\n',...
    theta0, theta0*180/pi)
fprintf('omega0 = %6.3g (rad/s)\n', omega0)
qt = {diff(theta,t),theta};
list0 = {omega0, theta0} ;
vAn_ = subs(vA_, qt, list0);
vA_ = double(vAn_);
fprintf('vA_ = [%6.3g, %6.3g, %g] (m/s)\n', vA_)

3.3 弹性接触阶段

delta = L*sin(theta(t))-L*sin(theta0);
Fe = k*abs(delta)^alpha;
Fe_ = [0, Fe, 0];
drA_ = diff(rA_,theta(t));
Q = drC_*G_.'+drA_*Fe_.';
Lagrange = LHS - Q;
[V,S] = odeToVectorField(Lagrange);
syms Y
eome_ = matlabFunction(V,'vars', {t,Y});
[tte, xse, te, yee] = ...
    ode45(eome_,time0, [theta0 omega0],option0);
tf = te(2);
thetaf = yee(2);
omegaf = yee(4);
fprintf('final conditions for contact\n')
fprintf('tf = %6.3g (s) \n',tf);
fprintf('thetaf = %6.3g (rad) = %g (degrees)\n',...
    thetaf, thetaf*180/pi)
fprintf('omegaf = %6.3g (rad/s)\n', omegaf)
thetan = real(xse(:,1));
omegan = real(xse(:,2));
deltan = L*(sin(thetan) - sin(theta0)*ones(length(thetan),1));

4. 总结

通过以上分析，我们可以看到：
- 对于球体与弹簧的问题，我们可以通过建立运动方程并使用MATLAB求解，得到球体的位移、速度和弹性力随时间的变化曲线。
- 对于旋转连杆在弹性力作用下的问题，我们分阶段分析了其运动，并使用牛顿 - 欧拉方法和拉格朗日方程得到了相同的结果。
- MATLAB提供了强大的工具来解决复杂的动力学问题，包括符号求解、数值求解和动画展示。

通过这些方法，我们可以更深入地理解机械系统的动力学特性，为工程设计和分析提供有力的支持。

关键步骤总结

球体与弹簧问题

步骤	操作
1	建立运动方程
2	使用 dsolve 求解符号解
3	输入数值参数并求解数值解
4	求解关键时间点和最大位移、最大力
5	绘制结果曲线
6	将二阶微分方程改写为一阶系统并使用 ode45 求解

旋转连杆问题

步骤	操作
1	设置连杆参数和符号变量
2	计算质心和端点的位置向量、转动惯量
3	分阶段分析运动，建立运动方程
4	使用 ode45 求解运动方程
5	使用拉格朗日方程求解并验证结果
6	展示连杆运动的动画

流程图

graph TD;
    A[开始] --> B[球体与弹簧问题];
    B --> B1[建立运动方程];
    B1 --> B2[求解符号解];
    B2 --> B3[输入数值参数];
    B3 --> B4[求解数值解];
    B4 --> B5[求解关键时间点和最大量];
    B5 --> B6[绘制结果曲线];
    B6 --> B7[改写为一阶系统并求解];
    A --> C[旋转连杆问题];
    C --> C1[设置参数和符号变量];
    C1 --> C2[计算位置向量和转动惯量];
    C2 --> C3[分阶段分析运动];
    C3 --> C4[建立运动方程];
    C4 --> C5[求解运动方程];
    C5 --> C6[使用拉格朗日方程验证];
    C6 --> C7[展示动画];
    B7 --> D[结束];
    C7 --> D;

通过以上的分析和操作，我们可以全面地研究平面运动物体在不同力作用下的动力学特性。希望这些内容对你有所帮助！

平面运动物体的直接动力学分析与MATLAB实现

5. 技术点分析

在上述的动力学分析过程中，涉及到了多个关键的技术点，下面我们对这些技术点进行详细分析。

5.1 符号计算与数值计算

在MATLAB中，符号计算和数值计算是解决动力学问题的重要手段。
- 符号计算 ：使用 syms 定义符号变量， dsolve 求解符号微分方程， subs 进行符号替换等。例如在球体与弹簧问题中，通过符号计算得到了位移和速度的符号解：

syms x(t)
syms m k g h xo vo
assume(m,'positive') 
assume(k,'positive') 
assume(g,'positive') 
assume(xo,'positive') 
assume(vo,'positive') 
eqn = m*diff(x,2) + k*x - m*g;
v = diff(x);
xSol = dsolve(eqn, x(0)==xo, v(0)==vo);
vSol = diff(xSol,t);

• 数值计算 ：使用 vpasolve 求解数值方程， ode45 求解常微分方程的数值解等。如在求解球体达到最大位移的时间 $t_v$ 和分离时间 $t_f$ 时：

tv = vpasolve(vn, t, [t0 t2]);
tf = vpasolve(xn, t, [tv t2]);

5.2 微分方程的处理

对于二阶微分方程，通常需要将其转化为一阶系统进行求解。以球体与弹簧问题为例，二阶微分方程 $m\ddot{x} = mg - kx$ 转化为一阶系统：
[\dot{x}_1 = x_2]
[\dot{x}_2 = g - \frac{k}{m}x_1]
对应的MATLAB代码为：

syms t x1 x2
dx1 = x2;
dx2 = g-k*x1/m;
dx = [dx1; dx2];
g = matlabFunction(dx,'vars',{t, [x1; x2]});

5.3 事件检测

在求解微分方程时，常常需要检测特定的事件，如球体与弹簧分离的时刻。通过定义事件函数并使用 odeset 设置事件选项来实现。例如在球体与弹簧问题中，定义 eventx1.m 函数：

function [value,isterminal,direction] = eventx1(t,x)
value = x(1);
isterminal = 1;
direction = 0;
end

然后设置ODE求解器的选项：

optionx1 = odeset('RelTol',1e-6,'MaxStep',1e-6,...
    'Events',@eventx1);

6. 应用拓展

上述的动力学分析方法可以应用到更复杂的机械系统中，以下是一些可能的拓展应用：

6.1 多体系统动力学

对于多个物体组成的系统，可以将每个物体的运动方程联立起来，形成一个更大的方程组进行求解。例如，多个球体与弹簧的组合系统，或者多个连杆组成的机械臂系统。

6.2 非线性动力学问题

在实际应用中，很多系统存在非线性因素，如非线性弹簧、摩擦力等。可以通过修改运动方程，考虑这些非线性因素，然后使用数值方法求解。

6.3 优化设计

通过对动力学分析结果的研究，可以对机械系统进行优化设计。例如，调整弹簧的刚度、连杆的长度等参数，以达到更好的性能指标。

7. 注意事项

在使用MATLAB进行动力学分析时，需要注意以下几点：
- 符号变量的定义 ：在使用符号计算时，要确保符号变量的定义正确，并且在进行符号替换时，变量的对应关系要准确。
- 数值求解的精度 ：在使用数值方法求解时，要合理设置相对误差容限 RelTol 和最大步长 MaxStep 等参数，以保证求解的精度。
- 事件检测的准确性 ：事件函数的定义要准确，确保能够正确检测到所需的事件。

8. 总结与展望

本文详细介绍了平面运动物体的直接动力学分析方法，包括球体与弹簧问题和旋转连杆在弹性力作用下的问题。通过建立运动方程、使用MATLAB进行符号求解和数值求解，得到了物体的运动状态和相关物理量的变化曲线。同时，还介绍了拉格朗日方程的应用，并验证了其与牛顿 - 欧拉方法的一致性。

未来，随着机械系统的不断复杂化和智能化，动力学分析的需求也会越来越高。我们可以进一步研究更复杂的系统，考虑更多的因素，如材料的非线性、流体的影响等。同时，结合机器学习和人工智能技术，实现对机械系统的智能优化和故障诊断。

关键技术点总结

技术点	描述	应用场景
符号计算	使用 syms 、 dsolve 等函数进行符号变量定义和符号方程求解	求解运动方程的符号解
数值计算	使用 vpasolve 、 ode45 等函数进行数值方程求解和常微分方程数值解	求解关键时间点和运动方程的数值解
微分方程转化	将二阶微分方程转化为一阶系统	便于使用ODE求解器求解
事件检测	定义事件函数并使用 odeset 设置事件选项	检测特定事件，如分离时刻

拓展应用流程图

graph TD;
    A[动力学分析方法] --> B[多体系统动力学];
    A --> C[非线性动力学问题];
    A --> D[优化设计];
    B --> B1[联立运动方程];
    B1 --> B2[求解方程组];
    C --> C1[考虑非线性因素];
    C1 --> C2[修改运动方程];
    C2 --> C3[数值求解];
    D --> D1[分析结果];
    D1 --> D2[调整参数];
    D2 --> D3[验证性能];

通过对这些技术点的掌握和应用，我们可以更好地解决各种机械系统的动力学问题，为工程实践提供有力的支持。希望本文能够对你有所帮助，激发你对动力学分析的兴趣和探索。