《理论力学》课程教学资源（文献资料）Matlab求解理论力学问...列（二）典型机构的运动分析_高云峰

第43卷第3期力学与实践2021年6月育研家Matlab 求解理论力学问题系列(二)典型机构的运动分析高云峰1)（清华大学航天航空学院，北京100084)摘要运动学计算机辅助分析可以帮助人们方便了解机械系统的整个运动过程，获得系统中任意点的运动轨迹、速度和加速度的变化规律，直观明了。Matlab软件是实现计算机辅助分析的工具之一。本文着重介绍了2个典型机构的运动分析，一个案例侧重从数值计算的角度介绍Matlab如何求解机构运动涉及的非线性方程组和如何进行动画演示：另一个案例着重介绍如何利用Matlab中的符号推导功能，求解全参数化的运动学问题。关键词非线性代数方程组，牛顿一拉普森方法，运动轨迹，动画显示，符号推导KINEMATICSANALYSISOFTYPICALMECHANISMSGAO Yunfengl)(School of Aerospace Engineering, Tsinghua University, Beijing 100084, China)AbstractComputer-aided analysis of kinematics can help understand the whole movement process of amechanical system, along with the motion trajectory, the velocity and the acceleration of the system at anytime or position.Matlab can be used to realize the computer-aided analysis.This paper focuses on thekinematics analysis of two typical mechanisms. In one case, the focus is on how to solve the nonlinear algebraicequations by Matlab and how to perform the animation demonstration.In another case, the focus is on how tosolve the fully parametric kinematics problem by using the symbolic derivation function in Matlab.Key wordsnonlinear algebraic equations, Newton-Raphson method, motion trajectory,animation demon-stration, symbolic derivation在传统的运动学教学中，为了避免复杂的求导运动量。Matlab是实现计算机辅助分析的工具之一，而陷入纯数学演算，通常强调基点法和点的复合运可以帮助学生加深理解机械运动的特点，为将来学动等具有物理含义的方法，求解系统在特定时刻或生自己设计装置提供了方便。位置的运动信息]。这样处理的优点是强调物理概本文着重介绍2个典型机构的运动轨迹、速度念，缺点是对系统的整体运动不容易了解，某些点的和加速度分析，以及利用Matlab中的符号推导功运动轨迹有时难以想象。这样训练出来的学生，对于能，从位移函数求导得到速度和加速度的表达式，并机械运动缺乏了解，也难以胜任未来的设计工作。进行画图和动画显示。利用运动学计算机辅助分析方法和相关软件，可在传统点的复合运动中，要特别注意动点、动系以很容易演示系统的整个运动过程、求出任意点或刚的选择，如果选择不当，就分析不下去，因此对很多体在任意时刻的速度、加速度、角速度、角加速度等同学而言有一定的难度。而采用Matlab，可以采用2020-06-01收到第1稿，2020-07-23收到修改稿。1) E-mail: gaoyunfeng@tsinghua.edu.cn引用格式：高云峰。Matlab求解理论力学问题系列（二）典型机构的运动分析.力学与实践，2021，43(3)：414-419Gao Yunfeng.Kinematics analysis of typical mechanisms. Mechanics in Engineering,2021, 43(3):=414-419(C)1994-2022ChinaAcademic JournalElectronicPublishingHouse.Allrights reserved.http://www.cnki.net

第 43 卷 第 3 期 力 学 与 实 践 2021 年 6 月 Matlab 求解理论力学问题系列 (二) 典型机构的运动分析 高云峰1) (清华大学航天航空学院, 北京 100084) 摘要 运动学计算机辅助分析可以帮助人们方便了解机械系统的整个运动过程，获得系统中任意点的运 动轨迹、速度和加速度的变化规律，直观明了。Matlab 软件是实现计算机辅助分析的工具之一。本文着重介绍 了 2 个典型机构的运动分析，一个案例侧重从数值计算的角度介绍 Matlab 如何求解机构运动涉及的非线性方 程组和如何进行动画演示；另一个案例着重介绍如何利用 Matlab 中的符号推导功能，求解全参数化的运动学 问题。 关键词 非线性代数方程组，牛顿 − 拉普森方法，运动轨迹，动画显示，符号推导 KINEMATICS ANALYSIS OF TYPICAL MECHANISMS GAO Yunfeng1) (School of Aerospace Engineering, Tsinghua University, Beijing 100084, China) Abstract 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 paper 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 symbolic derivation function in Matlab. Key words nonlinear algebraic equations, Newton–Raphson method, motion trajectory, animation demon￾stration, symbolic derivation 在传统的运动学教学中，为了避免复杂的求导 而陷入纯数学演算，通常强调基点法和点的复合运 动等具有物理含义的方法，求解系统在特定时刻或 位置的运动信息[1]。这样处理的优点是强调物理概 念，缺点是对系统的整体运动不容易了解，某些点的 运动轨迹有时难以想象。这样训练出来的学生，对于 机械运动缺乏了解，也难以胜任未来的设计工作。 利用运动学计算机辅助分析方法和相关软件，可 以很容易演示系统的整个运动过程、求出任意点或刚 体在任意时刻的速度、加速度、角速度、角加速度等 运动量。Matlab 是实现计算机辅助分析的工具之一， 可以帮助学生加深理解机械运动的特点，为将来学 生自己设计装置提供了方便。 本文着重介绍 2 个典型机构的运动轨迹、速度 和加速度分析，以及利用 Matlab 中的符号推导功 能，从位移函数求导得到速度和加速度的表达式，并 进行画图和动画显示。 在传统点的复合运动中，要特别注意动点、动系 的选择，如果选择不当，就分析不下去，因此对很多 同学而言有一定的难度。而采用 Matlab，可以采用 2020–06–01 收到第 1 稿，2020–07–23 收到修改稿。 1) E-mail: gaoyunfeng@tsinghua.edu.cn 引用格式: 高云峰. Matlab 求解理论力学问题系列 (二) 典型机构的运动分析. 力学与实践, 2021, 43(3): 414-419 Gao Yunfeng. Kinematics analysis of typical mechanisms. Mechanics in Engineering, 2021, 43(3): 414-419

第3期415高云峰：Matlab求解理论力学问题系列（二）典型机构的运动分析统一的方法，只要列出一般位置的表达式就可以利步骤（4)：类似一元函数的泰勒展开式，f()=用符号推导功能，得到速度、加速度等运动量的表f(o)+f(o)(-)+o(-o)，多元函数为达式。fi(a) = fi(*) + J(a*)da +o(da)1Matlab中非线性方程的求解及动画演示略去高阶小量，从而得到关于修正值da=[dr1，dr2.，drnjT的一组线性方程组案例1：如图1，已知四连杆机构ABCD，AB杆长为a1,BC杆长为a2，CD杆长为a3，AD距离J(r*)dr + fi(a*)= 0为a4。若AB杆以匀角速度wo转动，初始%=0步骤(5)：解线性方程组，得到修正值da。求BC和CD杆的角度、角速度变化规律。步骤（6)：得到新的估计值a*=a*+da，返回到步骤(2)。Oa具体回到案例1，把方程(1)改写为1fi(1, p2)=a2cos P1+ a3cos(p2+ a1cos -a42f2(1, p2) = a2 sin 1 + a3 sin (p2+a1 sin 已知1,a2,3a4和，估计值是（i,)。计算雅可比矩阵，并代入当前估计值o-a2 sinpi -ag sinp2(3)J(p1,p2) =图1四连杆机构a2 cos (p1a3 cos 4P2得到一组关于修正值[dp1,dp2]T的线性方程组这个问题中尺寸不是特别给定的值，一般学生分析起来会有困难。而用Matlab可以进行规范化处-a2 sin gt -a3 sing2dpi理：建立坐标Ary，=wot已知，广义坐标为1和a2cospia3cosp2do(2。系统的约束方程为[2][ fi(i, ) (4)[2(e1,) = 10acos+a2cos1+a3cos2-a4=0(1)a1sin+ a2sin 1+a3sin p2=0类似代码x=inv(A)*B，可解出[dp1,dp2]T，总之按上述步骤1～6，可以解出任意时刻的角度。方程（1）是关于转角1和22的非线性方程组，通常没有解析解，下面给出一般的处理方法。在某些情况下，雅可比矩阵的行列式可能趋近于零或等于零，意味着方程（4）多解或无解。多解可设有非线性方程组f(1,2，.，n)，i=1,2,：几，没有解析解，可用多种数值计算的方法求解，能的原因是：机械系统由于设计导致在某个特定位置本身就是多解的，如图2中的曲柄－滑块机构，OA这里以牛顿－拉普森法为例，其主要步骤如下：运动带动B运动通常不会有问题：但反过来B运步骤（1)：给定计算中的允许误差ε，给出一组动带动OA运动（如蒸汽机），当OAB共线时，OA粗略的估计值*=[,2…，。下一时刻可以向上也可以向下（考虑动力学时OA杆步骤(2)：计算f:(a*)，判断If;(a*)I<是否成立，若成立，当前估计值即为一组数值解，求解结BU束。若不成立，到步骤(3)。0ee0步骤(3)：计算非线性方程组的雅可比矩阵J(a)，代入当前估计值得到J(e)。多元函数f(r1,.，n）的雅可比矩阵定义为E2.UBBofi/or1afi/or2...afi/orn00f2/0r1f2/0r2...0f2/0rnJ() =:3:Ofn/oriafn/or2...Ofn/orn图2无解或多解的情况(C)1994-2022ChinaAcademicJournalElectronicPublishingHouse.Allrights reserved.http://www.cnki.net

第 3 期 高云峰：Matlab 求解理论力学问题系列 (二) 典型机构的运动分析 415 统一的方法，只要列出一般位置的表达式就可以利 用符号推导功能，得到速度、加速度等运动量的表 达式。 1 Matlab 中非线性方程的求解及动画演示 案例 1：如图 1，已知四连杆机构 ABCD，AB 杆长为 a1，BC 杆长为 a2，CD 杆长为 a3，AD 距离 为 a4。若 AB 杆以匀角速度 ω0 转动，初始 θ0 = 0。 求 BC 和 CD 杆的角度、角速度变化规律。 a a a a ϕ ϕ y A B D C x θ 图 1 四连杆机构 这个问题中尺寸不是特别给定的值，一般学生 分析起来会有困难。而用 Matlab 可以进行规范化处 理：建立坐标 Axy，θ = ω0t 已知，广义坐标为 φ1 和 φ2。系统的约束方程为[2] a1 cos θ + a2 cos φ1 + a3 cos φ2 − a4 = 0 a1 sin θ + a2 sin φ1 + a3 sin φ2 = 0 } (1) 方程 (1) 是关于转角 φ1 和 φ2 的非线性方程组，通 常没有解析解，下面给出一般的处理方法。 设有非线性方程组 fi(x1, x2, · · · , xn), i = 1, 2, · · · , n，没有解析解，可用多种数值计算的方法求解， 这里以牛顿−拉普森法为例，其主要步骤如下： 步骤 (1)：给定计算中的允许误差 ε，给出一组 粗略的估计值 x ∗ = [x ∗ 1 , x∗ 2 , · · · , x∗ n] T。 步骤 (2)：计算 fi(x ∗ )，判断 |fi(x ∗ )| < ε 是否 成立，若成立，当前估计值即为一组数值解，求解结 束。若不成立，到步骤 (3)。 步骤 (3)：计算非线性方程组的雅可比矩阵 J(x)，代入当前估计值得到 J(x ∗ )。多元函数 fi(x1, x2, · · · , xn) 的雅可比矩阵定义为 J(x) =        ∂f1/∂x1 ∂f1/∂x2 · · · ∂f1/∂xn ∂f2/∂x1 ∂f2/∂x2 · · · ∂f2/∂xn . . . . . . . . . . . . ∂fn/∂x1 ∂fn/∂x2 · · · ∂fn/∂xn        步骤 (4)：类似一元函数的泰勒展开式，f(x) = f(x0) + f ′ (x0)(x − x0) + o(x − x0)，多元函数为 fi(x) = fi(x ∗ ) + J(x ∗ )dx + o(dx) 略去高阶小量，从而得到关于修正值 dx = [dx1, dx2, · · · , dxn] T 的一组线性方程组 J(x ∗ )dx + fi(x ∗ ) = 0 步骤 (5)：解线性方程组，得到修正值 dx。 步骤 (6)：得到新的估计值 x ∗ = x ∗ + dx，返回 到步骤 (2)。 具体回到案例 1，把方程 (1) 改写为 f1(φ1, φ2) = a2 cos φ1 + a3 cos φ2 + a1 cos θ − a4 f2(φ1, φ2) = a2 sin φ1 + a3 sin φ2 + a1 sin θ } (2) 已知 a1, a2, a3, a4 和 θ，估计值是 (φ ∗ 1 , φ∗ 2 )。计算雅可 比矩阵，并代入当前估计值 J(φ ∗ 1 , φ∗ 2 ) = [ −a2 sin φ ∗ 1 −a3 sin φ ∗ 2 a2 cos φ ∗ 1 a3 cos φ ∗ 2 ] (3) 得到一组关于修正值 [dφ1, dφ2] T 的线性方程组 [ −a2 sin φ ∗ 1 −a3 sin φ ∗ 2 a2 cos φ ∗ 1 a3 cos φ ∗ 2 ] [ dφ1 dφ2 ] + [ f1(φ ∗ 1 , φ∗ 2 ) f2(φ ∗ 1 , φ∗ 2 ) ] = [ 0 0 ] (4) 类似代码 X=inv(A)*B，可解出 [dφ1, dφ2] T，总之按 上述步骤 1 ∼ 6，可以解出任意时刻的角度。 在某些情况下，雅可比矩阵的行列式可能趋近 于零或等于零，意味着方程 (4) 多解或无解。多解可 能的原因是：机械系统由于设计导致在某个特定位 置本身就是多解的，如图 2 中的曲柄−滑块机构，OA 运动带动 B 运动通常不会有问题；但反过来 B 运 动带动 OA 运动 (如蒸汽机)，当 OAB 共线时，OA 下一时刻可以向上也可以向下 (考虑动力学时 OA 杆 vB B B l l A A r r O O θ ϕ ω 图 2 无解或多解的情况

416力学与实践2021年第43卷由于惯性继续运动，解是确定的，但仅考虑运动学各铰点的位置并连接起来，就得到了四连杆机构在某一时刻的图象，延迟一定的时间后再画出下一时时OA杆位置是多解的)。无解的可能是：机械系统刻的图象，就形成了动画。本问题中动画的源代码由于尺寸关系导致运动发生冲突，例如图2中1<r见图4，其中plot函数表示画线段；hl是句柄，定义时OA杆就不能连续转动，本问题中若四杆长度满足a2+a3<a1+a4，当=元时方程（1)的确无解-（物理解释是AB杆运动范围有限制，不能达到这一位置）。利用方程无解这一特点，未来自已在设计装置时，可以发现可能的运动冲突。编程计算得到角度的变化关系后，可以算出任意时刻各较的位置，以及BC杆上不同点的运动轨D迹（图3)：很明显B点轨迹是圆，C点轨迹是圆的部分（AB杆大范围运动时，CD杆只在小范围运-动，而在BC杆上不同的点轨迹就很复杂了。知道任意时刻系统的各铰位置，在屏幕上画出图3BC杆上不同点的运动轨迹plot（gx,gy，k，Linewidth，1)：holdon；%画地面hl=line（Color'，"k，LineStyle'，，Linewidth,2,Marker'，"o"，Markersize,3)；%画图格式一动画部分for ii-l:NUMset(hl,XData'，x(ii,:)，YData，y(ii,：))：%读取数据%显示范围axis([-30,50,-30,50]);axis equal%画出当前顿drawnow;%暂停pause(0.05);end图4动画部分的源代码截图动画中每一帧图的特性，如颜色、线型、宽度等：set对式（5）进一步求导得到角加速度的方程函数是读取数据：drawnow函数是画出当前顿的图案；pause函数表示暂停，让每一顿有一定的视觉残-a2p1 sin P1 - a3p2 sin P2 - a2pr cos 1-留以便更好地形成动画（类似放电影，1秒钟放24a3p2cosp2-aiwgscos0=0(6)格图像)[3]。a21cos41+a3p2cos p2-a20i sin1-a3p sin (p2-a1w sing= 02函数求导获得角速度、角加速度求出角度后，对约束方程（1）求导出现角速度，此时9,1,21,P2都是已知值，因此得到关于有1,2的一组线性方程组。类似X=inv(A)*B可解-a2p1 sin P1 - a3p2 sin p2 -aiwo sin = 0出角加速度，从而可以获得角加速度随时间或随(5a291 cos 1 + a3p2 cos (P2 + a1wo cos 0 = 0变化的关系(图6)。这只是一类典型机构的例子，采用本文介绍的由于9,1,2已在前面求出，因此得到关于1,92方法，任意机构的运动轨迹、刚体的角速度和角加速的一组线性方程组。类似X=inv(A)*B可解出角速度、某点的速度和加速度都可以类似求出。而这些结度，从而可以获得角速度随时间或随变化的关系(图 5)。果（图3、图5、图6）都是传统方法难以获得的。(C)1994-2022ChinaAcademicJournalElectronicPublishingHouse.Allrights reserved.http://www.cnki.net

416 力 学 与 实 践 2021 年 第 43 卷 由于惯性继续运动，解是确定的，但仅考虑运动学 时 OA 杆位置是多解的)。无解的可能是：机械系统 由于尺寸关系导致运动发生冲突，例如图 2 中 l < r 时 OA 杆就不能连续转动，本问题中若四杆长度满 足 a2 + a3 < a1 + a4，当 θ = π 时方程 (1) 的确无解 (物理解释是 AB 杆运动范围有限制，不能达到这一 位置)。利用方程无解这一特点，未来自己在设计装 置时，可以发现可能的运动冲突。 编程计算得到角度的变化关系后，可以算出任 意时刻各铰的位置，以及 BC 杆上不同点的运动轨 迹 (图 3)：很明显 B 点轨迹是圆，C 点轨迹是圆的 一部分 (AB 杆大范围运动时，CD 杆只在小范围运 动)，而在 BC 杆上不同的点轨迹就很复杂了。 知道任意时刻系统的各铰位置，在屏幕上画出 各铰点的位置并连接起来，就得到了四连杆机构在 某一时刻的图象，延迟一定的时间后再画出下一时 刻的图象，就形成了动画。本问题中动画的源代码 见图 4，其中 plot 函数表示画线段；h1 是句柄，定义 A B D C 图 3 BC 杆上不同点的运动轨迹 图 4 动画部分的源代码截图 动画中每一帧图的特性，如颜色、线型、宽度等；set 函数是读取数据；drawnow 函数是画出当前帧的图 案；pause 函数表示暂停，让每一帧有一定的视觉残 留以便更好地形成动画 (类似放电影，1 秒钟放 24 格图像) [3]。 2 函数求导获得角速度、角加速度 求出角度后，对约束方程 (1) 求导出现角速度， 有 −a2φ˙ 1 sin φ1 − a3φ˙ 2 sin φ2 − a1ω0 sin θ = 0 a2φ˙ 1 cos φ1 + a3φ˙ 2 cos φ2 + a1ω0 cos θ = 0 } (5) 由于 θ, φ1, φ2 已在前面求出，因此得到关于 φ˙ 1, φ˙ 2 的一组线性方程组。类似 X=inv(A)*B 可解出角速 度，从而可以获得角速度随时间或随 θ 变化的关系 (图 5)。 对式 (5) 进一步求导得到角加速度的方程 −a2φ¨1 sin φ1 − a3φ¨2 sin φ2 − a2φ˙ 2 1 cos φ1− a3φ˙ 2 2 cos φ2 − a1ω 2 0 s cos θ = 0 a2φ¨1 cos φ1 + a3φ¨2 cos φ2 − a2φ˙ 2 1 sin φ1− a3φ˙ 2 2 sin φ2 − a1ω 2 0 sin θ = 0    (6) 此时 θ, φ1, φ2, φ˙ 1, φ˙ 2 都是已知值， 因此得到关于 φ¨1, φ¨2 的一组线性方程组。类似 X=inv(A)*B 可解 出角加速度，从而可以获得角加速度随时间或随 θ 变化的关系 (图 6)。 这只是一类典型机构的例子，采用本文介绍的 方法，任意机构的运动轨迹、刚体的角速度和角加速 度、某点的速度和加速度都可以类似求出。而这些结 果 (图 3、图 5、图 6) 都是传统方法难以获得的

第3期高云峰：Matlab求解理论力学问题系列（二）典型机构的运动分析4172传统处理方法：取AB杆为动系，O为动点（图8）。通过速度和加速度分析（具体过程略），可以得到(se)/vo0（顺时针），=AB=V3u/4r2（逆时针）（7)WAB=2rB...dedtdp/dtU.WAEdpa/dt2345o160/rad图5角速度与的关系2.. e/dt?d/dt?图8速度分析da/dt?(z8peMtlab处理本问题时，只需要把一般位置角度的表达式写出来，然后利用符号求导的方法，就可以得-2到角速度和角加速度。角度以逆时针方向为正，设初始时刻AB杆垂直，根据图9有012346+wot三0/rad(8)图6角加速度与的关系3Matlab中的符号求导如果运动学问题中没有具体数值，全是参数，用Matlab也可以处理。案例2：半径为工的圆轮在水平桌面上作直线纯滚动，轮心速度o的大小为常数。摇杆AB与桌面铰接，并靠在圆轮上（图7)。当摇杆与桌面夹角$=60°时，试求摇杆的角速度和角加速度。+图9运动学关系根据式（8)有(9)(p=2arctan[r/(r+vot)对式（9)求一阶导数获得角速度，求二阶导数获图7一类接触问题得角加速度。程序的源代码见图10。(C)1994-2022ChinaAcademic JournalElectronicPublishingHouse.Allrights reserved.http://www.cnki.net

第 3 期 高云峰：Matlab 求解理论力学问题系列 (二) 典型机构的运动分析 417 0 1 2 3 4 5 6 7 θ/rad -2 1 0 1 2 dθ/dt dϕ/dt dϕ/dt 䀂䙏ᓖ/(rad.s-1 ) 图 5 角速度与 θ 的关系 0 1 2 3 4 5 6 7 θ/rad d  θ/dt  d ϕ/dt  d ϕ/dt  䀂࣐䙏ᓖ/(rad.s-2 ) -6 -4 -2 0 2 4 6 图 6 角加速度与 θ 的关系 3 Matlab 中的符号求导 如果运动学问题中没有具体数值，全是参数，用 Matlab 也可以处理。 案例 2：半径为 r 的圆轮在水平桌面上作直线 纯滚动，轮心速度 vO 的大小为常数。摇杆 AB 与 桌面铰接，并靠在圆轮上 (图 7)。当摇杆与桌面夹角 φ = 60◦ 时，试求摇杆的角速度和角加速度。 O vO A B ϕ 图 7 一类接触问题 传统处理方法：取 AB 杆为动系，O 为动点 (图 8)。通过速度和加速度分析 (具体过程略)，可 以得到 ωAB = vO 2r (顺时针), εAB = √ 3v 2 O/4r 2 (逆时针) (7) ve vO ωAB vr O A B 图 8 速度分析 Mtlab 处理本问题时，只需要把一般位置角度的 表达式写出来，然后利用符号求导的方法，就可以得 到角速度和角加速度。角度以逆时针方向为正，设初 始时刻 AB 杆垂直，根据图 9 有 x = r + vOt tan ( 1 2 φ ) = r/x    (8) ωAB O r x A B ϕ 图 9 运动学关系 根据式 (8) 有 φ = 2 arctan[r/(r + vOt)] (9) 对式 (9) 求一阶导数获得角速度，求二阶导数获 得角加速度。程序的源代码见图 10

418力学与实践2021年第43卷clc;%清除屏幕因此可以画出系统运动全程的角度、角速度、角%定义符号symsrphivo x t加速度变化曲线，取归一化参数r=1，vo=1，得一给出运动学关系到的曲线关系见图12和图13。x=r+vo*t:%质心的运动学关系phi=2*atan(r/x);%角度表达式一求导2.0%对中求一阶导数dphi dt=diff(phi,t'):“角度(-s-per)/3ddphi_dt2-diff(dphi_dt，t）；%对中求二阶导数角速度w1.5一显示结果角加速度disp(['d/dt=",char(dphi_dt))%显示角速度disp（d2小/dt2=，char(ddphi_dt2)】）%显示角加速度1.0*(_s.pe.) /m图10函数求导数的源代码0.5源代码中具体涉及几个函数：（1）符号定义是pea0“"syms"。（2）函数求导的格式是diff（function，variable')，表示函数对某个变量求导。运行后得到的结0.50.51.01.52.0果为0t/s2rvo(10)WAB2r2+ugt2+2rot图12系统运动参量与时间的关系4ru(r + vot)(11)EAB-(2r2 +t2 +2rvot)20.4可以看到，传统方法的分析要很大的篇幅，而Matlab只需要几行就得到了结果（顺便说一下，前0.2(e-8-pe.)/a *(f-s.面式(5)和式(6)也可以利用符号推导从式(1)得0到)。如果要求特定位置的值，只需要把具体值代入，0.2spe)格式为subs(function，old,new)，表示用新的变量替换老的变量，这部分程序的源代码见图11。角速度w-0.4角加速度一特定角度时结果-0.6%特定位置0.81.01.21.4phi60=60*pi/1801.6/raddistance=r/tan(phi60/2)-r;%特定角度时圆心位移%时间=距离/速度ts=distance/vo:图13系统运动参量与角度的关系%把特定时间代入角速度newl=subs(dphi_dt,t,ts);new2=subs(ddphi_dt2,t,ts):%把特定时间代入角加速度显示结果借助Matlab还可以对题且进行深入研究。例如，disp（[when中=60°d/dt=，char(new1)J）%显示角速度用传统方法分析时，如果动点选为圆盘上的接触点，disp（[when巾=60°d2/dt2=，char（new2)）%显示角加速度速度还可能会正确，但是加速度没有办法分析了，这图11求特定位置的运动量是为什么呢？可以看看接触点相对AB杆的轨迹，这借助计算机很容易实现。运行后屏幕显示结果为以=60°为系统初始位置，设圆盘上P点与when 重=60odΦ/dt=-v0/(2*r)AB杆接触，且AP=V3r。当圆心运动△r=otwhenΦ=60°d24/dt2=(3^(1/2)*时，圆盘因为纯滚动的转角为6=△/r（图14)，且v02)/(4*^2)P点在动系中有注意负号表示顺时针方向，因此计算机推导出的结果与传统方法相同，见式(7)。p = r(1 - cos0)利用参数替换函数，还可以把式(1)和式(11)(12)yp=V3r+A+rsingV3r + r( + sin 0)中的时间用角度替换(具体略)。(C)1994-2022 China Academic Journal Electronic Publishing House.All rights reserved.http://www.cnki.net

418 力 学 与 实 践 2021 年 第 43 卷 图 10 函数求导数的源代码 源代码中具体涉及几个函数：(1) 符号定义是 “syms”。(2) 函数求导的格式是 diff (function, ’vari￾able’)，表示函数对某个变量求导。运行后得到的结 果为 ωAB = − 2rv0 2r 2 + v 2 0 t 2 + 2rv0t (10) εAB = 4rv2 0 (r + v0t) (2r 2 + v 2 0 t 2 + 2rv0t) 2 (11) 可以看到，传统方法的分析要很大的篇幅，而 Matlab 只需要几行就得到了结果 (顺便说一下，前 面式 (5) 和式 (6) 也可以利用符号推导从式 (1) 得 到)。 如果要求特定位置的值，只需要把具体值代入， 格式为 subs (function, old, new)，表示用新的变量替 换老的变量，这部分程序的源代码见图 11。 图 11 求特定位置的运动量 运行后屏幕显示结果为 when Φ = 60◦ dΦ/dt = −v0/(2 ∗ r) when Φ = 60◦ d2Φ/dt2 = (3∧(1/2) ∗ v0∧2)/(4 ∗ r ∧2) 注意负号表示顺时针方向，因此计算机推导出 的结果与传统方法相同，见式 (7)。 利用参数替换函数，还可以把式 (1) 和式 (11) 中的时间用角度替换 (具体略)。 因此可以画出系统运动全程的角度、角速度、角 加速度变化曲线，取归一化参数 r = 1，vO = 1，得 到的曲线关系见图 12 和图 13。 0 0.5 1.0 1.5 2.0 t/s φ/rad, ω/(rad.s-1), ε/(rad.s-2 ) -0.5 0 0.5 1.5 1.0 2.0 䀂ᓖφ 䀂䙏ᓖω 䀂࣐䙏ᓖε 图 12 系统运动参量与时间的关系 0.8 1.0 1.2 1.4 1.6 -0.6 -0.4 -0.2 0 0.2 0.4 䀂䙏ᓖω 䀂࣐䙏ᓖε φ/rad ω/(rad.s-1), ε/(rad.s-2 ) 图 13 系统运动参量与角度的关系 借助 Matlab 还可以对题目进行深入研究。例如， 用传统方法分析时，如果动点选为圆盘上的接触点， 速度还可能会正确，但是加速度没有办法分析了，这 是为什么呢？可以看看接触点相对 AB 杆的轨迹，这 借助计算机很容易实现。 以 φ = 60◦ 为系统初始位置，设圆盘上 P 点与 AB 杆接触，且 AP = √ 3r。当圆心运动 ∆x = vOt 时，圆盘因为纯滚动的转角为 θ = ∆x/r (图 14)，且 P 点在动系中有 xP = r(1 − cos θ) yP = √ 3r + ∆x + r sin θ = √ 3r + r(θ + sin θ)    (12)

第3期419高云峰：Matlab求解理论力学问题系列（二）典型机构的运动分析在动系中画出的轨迹如图15所示，很像是旋轮的学生自己也不能说得很清楚，只是凭感觉），碰巧线。如果P点的轨迹是旋轮线，则曲线的尖点应该是正确的；但是在加速度分析时，由于相对运动的轨接触AB杆：而实际上P点接触杆时，相对速度不迹不是太清楚，相对切向加速度未知，相对法向加速为零，但轨迹的切线方向与杆平行，不过该曲线的曲度也未知，再加上AB杆的角加速度未知，未知数率在详细分析前还不清楚。太多求解不了。注意在地面上看P点的轨迹是旋轮线，通过反y转和平移，旋轮线与图15中的轨迹完全相同，这个BP?结论是否出乎意料？4小结在运动学中引入Matlab，可以方便了解系统整00O个运动过程，获得系统中任意点的运动轨迹，以及速度和加速度随时间或某角度的变化规律，直观明了让学生更容易理解系统的运动。42Lr本篇着重介绍了非线性方程组的求解方法，以及画图和动画的格式；另外介绍了参数方程的求导图14接触点的轨迹分析和替换。画图涉及到的Matlab函数是plot（画直+ 9./m线），而句柄函数、set函数和drawon函数的组合就可以实现动画；另一个重要函数是diff（函数求导）接触点轨迹AB杆和subs（变量替换）。掌握这儿个命令后，就可以对一般的运动学问题进行全程的计算和动画演示，获得丰富的运动学信息。更进一步，借助Matlab可以更加深入研究运动学的问题，解决传统分析方法无法处理的很多问题。参考文献Aa./m1李俊峰等。理论力学，北京：清华大学出版社，2003图15接触点的轨迹2高云峰.理论力学辅导与习题集，北京：清华大学出版社，20033甘勤涛，胡仁喜，程政田等，Matlab数学计算与工程分析这就回答了前面的问题：如果选P点为动点，从入门到精通.北京：机械工业出版社，2019在分析速度时认为相对速度沿杆方向（相信这样做(贵任编辑：胡漫)(C)1994-2022ChinaAcademicJournalElectronicPublishingHouse.Allrights reserved.http://www.cnki.net

第 3 期 高云峰：Matlab 求解理论力学问题系列 (二) 典型机构的运动分析 419 在动系中画出的轨迹如图 15 所示，很像是旋轮 线。如果 P 点的轨迹是旋轮线，则曲线的尖点应该 接触 AB 杆；而实际上 P 点接触杆时，相对速度不 为零，但轨迹的切线方向与杆平行，不过该曲线的曲 率在详细分析前还不清楚。 yr Dx vO xr B A O P θ θ ϕ 图 14 接触点的轨迹分析 xr/m yr/m A ᧕䀖⛩䖘䘩 ABᵶ 图 15 接触点的轨迹 这就回答了前面的问题：如果选 P 点为动点， 在分析速度时认为相对速度沿杆方向 (相信这样做 的学生自己也不能说得很清楚，只是凭感觉)，碰巧 是正确的；但是在加速度分析时，由于相对运动的轨 迹不是太清楚，相对切向加速度未知，相对法向加速 度也未知，再加上 AB 杆的角加速度未知，未知数 太多求解不了。 注意在地面上看 P 点的轨迹是旋轮线，通过反 转和平移，旋轮线与图 15 中的轨迹完全相同，这个 结论是否出乎意料？ 4 小结 在运动学中引入 Matlab，可以方便了解系统整 个运动过程，获得系统中任意点的运动轨迹，以及速 度和加速度随时间或某角度的变化规律，直观明了， 让学生更容易理解系统的运动。 本篇着重介绍了非线性方程组的求解方法，以 及画图和动画的格式；另外介绍了参数方程的求导 和替换。画图涉及到的 Matlab 函数是 plot (画直 线)，而句柄函数、set 函数和 drawon 函数的组合就 可以实现动画；另一个重要函数是 diff (函数求导) 和 subs (变量替换)。掌握这几个命令后，就可以对 一般的运动学问题进行全程的计算和动画演示，获 得丰富的运动学信息。更进一步，借助 Matlab 可以 更加深入研究运动学的问题，解决传统分析方法无 法处理的很多问题。 参 考 文 献 1 李俊峰等. 理论力学. 北京: 清华大学出版社, 2003 2 高云峰. 理论力学辅导与习题集. 北京: 清华大学出版社, 2003 3 甘勤涛, 胡仁喜, 程政田等. Matlab 数学计算与工程分析 —— 从入门到精通. 北京: 机械工业出版社, 2019 (责任编辑: 胡 漫)