Matlab求解理论力学问题系列(二)典型机构的运动分析

1、Matlab求解理论力学问题系列（二）典型机构的运动分析摘要运动学计算机辅助分析可以帮助人们方便了解机械系统的整个运动过程，获得系统中任意点的运 动轨迹、速度和加速度的变化规律，直观明了。Matlab软件是实现计算机辅助分析的工具之一。本文着重介绍 了 2个典型机构的运动分析，一个案例侧重从数值计算的角度介绍Matlab如何求解机构运动涉及的非线性方 程组和如何进行动画演示；另一个案例着重介绍如何利用Matlab中的符号推导功能，求解全参数化的运动学 问题。关键词 非线性代数方程组，牛顿-拉普森方法，运动轨迹，动画显示，符号推导KINEMATICS ANALYSIS OF TYPICAL ME

2、CHANISMSAbstract Computer-aided analysis of kinematics can help understand the whole movement process of a mechanical system, along with the motion trajectory, the velocity and the acceleration of the system at any time or position. Matlab can be used to realize the computer-aided analysis. This pap

3、er focuses on the kinematics analysis of two typical mechanisms. In one case, the focus is on how to solve the nonlinear algebraic equations by Matlab and how to perform the animation demonstration. In another case, the focus is on how to solve the fully parametric kinematics problem by using the sy

4、mbolic derivation function in Matlab.Key words nonlinear algebraic equations, Newton-Raphson method, motion trajectory, animation demonstration, symbolic derivation在传统的运动学教学中，为了避免复杂的求导 而陷入纯数学演算，通常强调基点法和点的复合运 动等具有物理含义的方法，求解系统在特定时刻或 位置的运动信息1。这样处理的优点是强调物理概 念，缺点是对系统的整体运动不容易了解，某些点的 运动轨迹有时难以想象。这样训练出来的学生，对

5、于 机械运动缺乏了解，也难以胜任未来的设计工作。利用运动学计算机辅助分析方法和相关软件，可 以很容易演示系统的整个运动过程、求出任意点或刚 体在任意时刻的速度、加速度、角速度、角加速度等 运动量。Matlab是实现计算机辅助分析的工具之一, 可以帮助学生加深理解机械运动的特点，为将来学 生自己设计装置提供了方便。本文着重介绍2个典型机构的运动轨迹、速度 和加速度分析，以及利用Matlab中的符号推导功 能，从位移函数求导得到速度和加速度的表达式，并 进行画图和动画显示。在传统点的复合运动中，要特别注意动点、动系 的选择，如果选择不当，就分析不下去，因此对很多 同学而言有一定的难度。而采用Mat

6、lab，可以采用统一的方法，只要列出一般位置的表达式就可以利 用符号推导功能，得到速度、加速度等运动量的表 达式。1 Matlab中非线性方程的求解及动画演示案例1：如图1,已知四连杆机构ABCD, AB 杆长为a1, BC杆长为a, CD杆长为a3, AD距离 为a4。若AB杆以匀角速度钏转动，初始Oo = 0。 求BC和CD杆的角度、角速度变化规律。图1四连杆机构这个问题中尺寸不是特别给定的值，一般学生 分析起来会有困难。而用Matlab可以进行规范化处 理：建立坐标Axy, O = t已知，广义坐标为初和 02。系统的约束方程为2a1 cos O + a2 cos 01 + a3 cos

7、 02 a4 = 0a1 sin O + a2 sin 01 + a3 sin 02 = 0 J方程(1)是关于转角01和公的非线性方程组，通 常没有解析解，下面给出一般的处理方法。设有非线性方程组fi (xi x, x), i = 1, 2, . ,n,没有解析解，可用多种数值计算的方法求解, 这里以牛顿-拉普森法为例，其主要步骤如下：步骤(1)：给定计算中的允许误差8,给出一组 粗略的估计值x* = x x, ,xn T。步骤(2)：计算fi(x*),判断jfi(x*)1 8是否 成立，若成立，当前估计值即为一组数值解，求解结 束。若不成立，到步骤(3)。步骤(3)：计算非线性方程组的雅可

8、比矩阵 J(x),代入当前估计值得到J(x*)。多元函数fi(xi, x2，,xn)的雅可比矩阵定义为J (x) =dfi/dxi df /dxiJ (x) =dfi/dxi df /dxidfn/dx1dfn/dX2 . dfdfn/dx1dfn/dX2 . dfn/dXn步骤(4)：类似一元函数的泰勒展开式，f(x)= f (xo) + r(xo)(x - xo) + o(x xo),多元函数为fi (x) = fi (x*) + J (x* )dx + o(dx)略去高阶小量，从而得到关于修正值dx =dxi, dx2, , dxnT的一组线性方程组J (x* )dx + fi (x*)

9、 = 0步骤(5)：解线性方程组，得到修正值dx。步骤(6)：得到新的估计值x* = x* +dx,返回 到步骤(2)。具体回到案例1,把方程(1)改写为fi(01,02) = a2 cos 01 + a3 cos 02 + ai cos O a4 f2 (01,02) = a2 sin 0i + a3 sin 02 + ai sin O J 已知ai,a2,a3a和O,估计值是(0*,0*)。计算雅可 比矩阵，并代入当前估计值J(0*,0*)= F s、a3 的0a2 cos i* a3 cos 2*(3)得到一组关于修正值d01, d02T的线性方程组(3)a2 sin 0* a3 sin

10、 0*a2 cos 1* a3 cos 2*fi (0*, 0*)= 0f2 (0*, 02)|_。类似代码X=inv(A)*B,可解出d0 i, d02T，总之按 上述步骤16,可以解出任意时刻的角度。在某些情况下，雅可比矩阵的行列式可能趋近 于零或等于零，意味着方程(4)多解或无解。多解可 能的原因是：机械系统由于设计导致在某个特定位 置本身就是多解的，如图2中的曲柄-滑块机构，OA 运动带动B运动通常不会有问题；但反过来B运 动带动OA运动(如蒸汽机)，当OAB共线时，OA 下一时刻可以向上也可以向下(考虑动力学时OA杆由于惯性继续运动，解是确定的，但仅考虑运动学 时OA杆位置是多解的）

11、。无解的可能是：机械系统 由于尺寸关系导致运动发生冲突，例如图2中lr 时OA杆就不能连续转动，本问题中若四杆长度满 足a + a (6) (6)动画中每一帧图的特性，如颜色、线型、宽度等；set 函数是读取数据；drawnow函数是画出当前帧的图 案；pause函数表示暂停，让每一帧有一定的视觉残 留以便更好地形成动画（类似放电影，1秒钟放24 格图像）同。2函数求导获得角速度、角加速度求出角度后，对约束方程（1）求导出现角速度, 有a2 0 1 sin 91 a3 02 sin 92 a1 sin 0 = 0a2 0 1 cos 01 + a3 02 cos 02 + a1 3 cos 0

12、 = 0由于0,01,02已在前面求出，因此得到关于0 1,02 的一组线性方程组。类似X=inv（A）*B可解出角速 度，从而可以获得角速度随时间或随0变化的关系 （图 5）。对式(5)进一步求导得到角加速度的方程a2 01 sin 01 a3 02 sin 02 a2 02 cos 01 a3 02 cos 02 a1 32 s cos 0 = 0 a2 01 cos 01 + a3 02 cos 02 a2 02 sin 01 一a3 02 sin 02 a1 32 sin 0 = 0此时0,01,02,01,02都是已知值，因此得到关于 01,02的一组线性方程组。类似X=inv（A）

13、*B可解 出角加速度，从而可以获得角加速度随时间或随0 变化的关系（图6）。这只是一类典型机构的例子，采用本文介绍的 方法，任意机构的运动轨迹、刚体的角速度和角加速 度、某点的速度和加速度都可以类似求出。而这些结 果（图3、图5、图6）都是传统方法难以获得的。-210101234567O/rad传统处理方法：取AB杆为动系，O为动点 （图8）。通过速度和加速度分析（具体过程略），可 以得到AB = 7 (顺时针-210101234567O/rad传统处理方法：取AB杆为动系，O为动点 （图8）。通过速度和加速度分析（具体过程略），可 以得到AB = 7 (顺时针)，AB = V3VO/4r2

14、(逆时针)(7)图5角速度与e的关系(TS.-peJ)/败炭Mtlab处理本问题时，只需要把一般位置角度的 表达式写出来，然后利用符号求导的方法，就可以得 到角速度和角加速度。角度以逆时针方向为正，设初 始时刻AB杆垂直，根据图9有x = r + vo ttan=r/x(8)如果运动学问题中没有具体数值，全是参数，用 Matlab也可以处理。案例2：半径为尸的圆轮在水平桌面上作直线 纯滚动，轮心速度vo的大小为常数。摇杆AB与 桌面铰接，并靠在圆轮上（图7）。当摇杆与桌面夹角 9 = 60。时，试求摇杆的角速度和角加速度。根据式（8）有9 = 2 arctanr/(r + vt)(9)对式（9

15、）求一阶导数获得角速度，求二阶导数获 得角加速度。程序的源代码见图10。(9)clc;%清除屏幕syms r phi vO x t % 定义符号 %给出运动学关系x=r+vO*t;%质心的运动学关系phi=2*atan(r/x) ;% 角度表达式%求导dphi_dt=diff(phi,t);% 对巾求一阶导数ddphi_dt2=diff (dphi_dt,t); % 对巾求二阶导数 % 显示结果 clc;%清除屏幕syms r phi vO x t % 定义符号 %给出运动学关系x=r+vO*t;%质心的运动学关系phi=2*atan(r/x) ;% 角度表达式%求导dphi_dt=diff(

16、phi,t);% 对巾求一阶导数ddphi_dt2=diff (dphi_dt,t); % 对巾求二阶导数 % 显示结果 disp(, d(b/dt=, char (dphi_dt)%显示角速度disp(，d2(!/dt2=, char (ddphi_dt2)%显示角加速度图10函数求导数的源代码因此可以画出系统运动全程的角度、角速度、角 加速度变化曲线，取归一化参数r = 1，vo = 1，得 到的曲线关系见图12和图13。源代码中具体涉及几个函数：(1)符号定义是 “syms”。(2)函数求导的格式是 diff (function, variable)， 表示函数对某个变量求导。运行后得到

17、的结 果为2.01.51.00.50角度角速度刃 角加速度&00.51.01.52.0t/s-0.52rvoAB2r2 + Vq t2 + 2rv01_4rv2 (r + Vq t)AB(2r2 + Vq t2 + 2rv01)2(10)(11)图12系统运动参量与时间的关系如果要求特定位置的值，只需要把具体值代入, 格式为subs (function, old, new)，表示用新的变量替 换老的变量，这部分程序的源代码见图11。%特定角度时结果phi60=60*pi/180;%特定位置distance=r/tan (phi60/2) -r;% 特定角度时圆心位移ts=distance/vO

18、;%时间=距离/速度newl=subs (dphi_dt, t, ts) ;%把特定时间代入角速度new2=subs (ddphi_dt2, t, ts);%把特定时间代入角加速度% 显不结果disp(f when 6=60 d d /dt=, char (newl) ) % 显示角速度 dispCC when 6=60 d2 /dt2=, char (new2) )% 显示角加速度图11求特定位置的运动量运行后屏幕显示结果为when 中=60。d/dt = v0/(2 * r)when 中=60。d2/dt2 = (3A (1/2) *v0A2)/(4 * rA2)注意负号表示顺时针方向，因

19、此计算机推导出 的结果与传统方法相同，见式(7)。利用参数替换函数，还可以把式(1)和式(11) 中的时间用角度替换(具体略)。/rad图13/rad图13系统运动参量与角度的关系c(ts幅j)/3借助Matlab还可以对题目进行深入研究。例如, 用传统方法分析时，如果动点选为圆盘上的接触点, 速度还可能会正确，但是加速度没有办法分析了，这 是为什么呢？可以看看接触点相对AB杆的轨迹，这 借助计算机很容易实现。以9 = 60。为系统初始位置，设圆盘上P点与 AB杆接触，且AP = pr当圆心运动= vot 时，圆盘因为纯滚动的转角为3 = Ax/r (图14)，且 P点在动系中有xP = r (1 cos 3)yP = y/3r + Ax + r sin 3 = (12)y/3r + r(3 + sin 3)在动系中画出的轨迹如图15所示，很像是旋轮