理论力学教具DIY系列(四)一个教具的受力分析与讨论

1、理论力学教具DIY系列（四）个教具的受力分析与讨论受力分析是理论力学教学中的重要内容，而多 点摩擦问题是其中的难点。在目前的教学安排中，静 力学有少量求解多点摩擦的习题，学生可以画图、列 出方程并求解；但是在动力学中，几乎没有涉及多点 摩擦的习题，而且也不要求解方程。一个实际问题如何简化为力学模型、如何分析 列出方程、如何求解、如何验证，是工科学生必需 掌握的基本技能。当然一门课程不可能教会学生处 理实际问题的所有环节或流程，但在适当的情况下 让学生了解这一过程，让部分感兴趣的学生深入探 究，还是有可能的。因此我在教学中有意识地设计 了一些教具，让学生有机会分析处理一个“完整问 题”问题的简化

2、、分析、计算、制作模型验证。1实际问题：一个小巧的力学演示装置这是一个专门设计的装置（图1），利用激光切割 后，加上减速电动机、电池和开关，就可以拼装出一 个爬行装置，且装置尾巴长度可调（图2）。图1爬行装置的设计图图2实际的爬行装置为了突出“完整问题”的分析过程而不陷入复 杂的计算，可以通过调整重心位置让装置在平面内 运动。观察发现，装置的运动状况依赖于它的尾长: 短尾时轮子一部分时间打滑，一部分时间不打滑，装 置前进得比较慢；长尾时轮子不打滑，前进得稍快 些。定性解释很容易：尾长影响重心相对位置，导致 尾部和轮子处的压力不同。尾短时重心靠近尾部，轮 子处压力小会打滑；尾长时重心靠近轮子，轮

3、子处压 力大不会打滑。但是具体分析计算就复杂得多，轮子旋转时，轮 子和尾部受到的压力、摩擦力与装置的加速度和角 加速度有关，不容易直观判断哪一点会打滑。为了进 行具体的分析计算，首先就需要建立一个力学模型。2力学模型：如何合理简化？设OXY是地面参考坐标系，Cxy是随体坐标 系，C是装置主体（不包轮子）的质心，任意时刻，尾 巴A和轮子边缘B接触地面。装置的几何尺寸参数为：尾巴长度Q1 ；质心距离尾部如，距离底部hi, 轴心O1距离尾部+。3，距离底部h ；轮长为l(3)xoi = xC + a3 cos a (h2 h1) sin ayoi = yc + a sin a + (h2 h1) c

4、os a(3)(图 3) oxb = xoi + l sin Q yb = yoi l cos 8(4)图3装置的力学模型图3装置的力学模型(2)装置分为两部分，主体和轮子。设主体质量为 mi，轮子质量为m2；主体质心C的坐标为xc，yc； 车身俯仰角度为a (绕z轴)；观察发现车轮以匀角 速度曲顺时针转动，设t时刻8 = n/3 - !ot，因此 8 2 (-冗/3,n/3)为某轮接触地面时的角度范围。由 于A, B始终接触地面，可以得到a与轮子转角8 的关系，根据轴心Oi距离水平面的高度，有(a1 + a2 + a3) sin a + h2 cos a = l cos 8(1)从式解出a，

5、求导后a、a也都已知。类似可以 得到yc, yoi与a或8的关系。为了保证装置能正常运动，轮子长度需要满足 2个条件：l cos(n/3) h2，保证轮子总能接触地 面；(2)轮子不能太长，保证重心总落在A、B之间, 否则装置会翻倒。3数学模型装置在平面内运动时，涉及的运动参数有 xc，yc，a，力参数有Fa，NA，Fb, Nb。注意到 装置中A、B两点接触地面，哪点会打滑是多点摩 擦问题中的难点。通常假设某点不打滑，然后得到计 算结果，再判断结果是否合理。在本问题中，存在三 种模式。3.1模式一(假设A不动，B打滑)设t时刻A点坐标xa(t) = xA，yA = 0已知, 根据图3,其他各点

6、的坐标为xc = xA + (ai + a2) cos a hi sin ayc = (ai + a2) sin a + hi cos a相应各点速度、加速度都可以求出，从而各点的位 置、速度、加速度都是已知函数。根据受力图(图4)，利用动静法，可以列出X 和Y方向力的平衡方程、B点打滑的摩擦力补充方 程、对A点力矩的平衡方程，为mi xc m2 xoi Fa + Fb = 0 mi yc m2 yoi + Na + Nb 一mi g m2 g = 0Fb = Nb Jia + nb xab mi(g + yc )xac 一m2 (g + yoi )xaoi + mi xc yc+ m2 xo

7、i yoi = 0图4装置的受力分析这时方程(5)中所有的运动学参数都是已知量, 因此是代数方程，容易解出。解出后要验证是否满足(1)所有压力、摩擦力 均大于0； (2) Fa 6 N 是否成立。如果条件都满 足，说明模式一假设成立(A不动，B打滑)，然后对 t + At时刻进行计算；否则说明模式一不成立，转入 模式二。3.2模式二(假设B不动，A打滑)设t时刻B点坐标xb(t) = xB，yB = 0已知, 根据图3,其他各点的坐标为xoi = xB l sin 8； yoi = l cos 8(6)xc = xoi a3 cos a + (h2 hi) sin a (7) yc = yoi

8、 a3 sin a (h2 hi) cos a IxA = xc (ai + a2) cos a + hi sin al(8) yA = yc (al + a2) sin a hi cos a相应各点的速度、加速度都可以求出，从而各点的位 置、速度、加速度都是已知函数。所有压力、摩擦力均大于0；Fa 6 N，Fb 6g 是否成立。如果条件都成立说明假设模式三成根据受力图（图4），利用动静法，系统方程为 立（A和B都打滑），然后对t +追时刻进行计算;mi xc m2 xoi Fa + Fb = 0 m yc m2 yoi + Na + Nb 一m1 g m2 g = 0Fa =旧Na& NA

9、xAB + m1 (g + yC )xCB + m2 (g + yoi )xoiB + mi xc yc+ m2 xoi yoi = 0方程（9）也是代数方程，解出后要验证是否满足 所有压力、摩擦力均大于0；Fb 6 Nb是否 成立。如果条件都满足，说明假设模式二成立（B不 动，A打滑），然后对t + At时刻进行计算；否则说 明模式二不成立，转入模式三。3.3模式三（假设A、B均打滑）在这种模式下xc，xoi成为未知量，利用运动 学关系把xoi用xc表示，因此得到一组微分代数 方程（mi + m2）xc + Fa Fb = fi （a；如 a）Na + Nb = f2 （a； a； a）（m

10、iyc + m2yoi）xc Naxab = fs（a； a； a）（1。）Fa INa = 0Fb iNb = 0方程（10）中fi, f2, f3是a, a, a的函数，5个方程5个未知数，可以求解。求解结束后要验证是否满足 否则说明公式或程序中可能存在错误，要返回检查。 但是如何求解微分代数方程，需要一些技巧，在 下面数值计算中会介绍。4数值计算有了方程，可以算出装置运动和受力随时间变 化规律，然后把数据进行可视化处理。数值求解前通 常先要画出计算框图，考虑各种可能情况（图5）。下 面采用Matlab程序进行计算，可以在计算结束后直 观从屏幕上查看虚拟装置如何运动，方便与实际装 置进行比

11、较。Matlab程序采用矩阵和列阵计算很方便，可以 把前面的方程统一改写为AX = B的形式，例如方图5求解方程的计算框图程（10）图5求解方程的计算框图10-10(mi + m2)010100 xAB00(mi yc + m2 yoi)1mu000001mu0(Fa(fi (a, a, a)、Naf2 (a, aFb=f3(a, a, a)Nb0I xc/0(11)在模式一、模式二中，调用X = inv(A)*B的 命令格式就完成了求解。但模式三的微分代数方程 处理起来要复杂一些：先把Xc当作代数量，用 X = inv(A)*B求解出来，设解出Xc = XC (这时 XC是具体数值)，再把X

12、c = XC当作微分方程，化 为标准的一阶微分方程组后，直接调用荣格库塔法 (Runge-Kutta methods)求解常微分方程。具体过程是：设t时刻系统各个参数(位置、速 度)都已经算出，现在要通过微分方程Xc = XC计 算出t + At时刻的参数。设yi = xc，y?=加，这样 就把二阶微分方程化为标准的一阶微分方程组以及 相应的初始条件，为yi=*(12a)y? = xcyi(0)= Xc (t)(12b)y?(0) = xc (t) J积分时间段为t,t + At,求解常微分方程调用的格式为t,iy = ode45rg_kt：time, time + hh, y0, optio

13、ns); (13) 程序式(13)中参数按顺序含义为：“t”是存放积分 时间；“iy”是存放积分结果；“ode45”是求解常微分 方程的标准函数；“rg_kt”是微分方程表达式的子函 数；“time, time+hh”是积分时间段；“y”是初始 值；“options”是积分的误差选项。值得说明的是，考虑到数值计算存在误差，因0。0A A A A/产-A打滑BB打滑0。0A A A A/产-A打滑BB打滑AA和B者；不打滑L不打滑 狩丁滑X-O-%,C 八 c10080604020102030405060尾长/mm图7显示了 ai = 5 mm时主体质心加速度的变 化，一方面看出模式转换时加速度

14、会有跳变，另一方 面可以看出质心加速度比重力加速度小3个数量级, 因此可以认为装置运动是“准静态”过程，从而可以 从静力学角度近似估计出aj：如果轮子在x方向距 离质心最远时有Nb Na，则龙点会一直要打滑。 分别对图4中的A和B取矩，有Naxabmig l sin：冗 + ag cos a。一(h hi) sin qq + mglNaxab(15)Nbxab = mig (a? +ai) cos qq hi sin qq +(15)mg (ai + a? + ag) cos qq h? sin qq考虑到实际上q是小量，从式(15)中近似得到a l sin 3n - a? + (mi - m

15、2)a3/(mi + m?) (16)此在判断摩擦关系是否违约时要稍微放松一点要求, 假设积分允许误差为= 1.0 x 10-6,则把摩擦关系F 6 uN放松为代入装置的数据，得到aj x 39 mm，接近图6中的40 mm。F 6 (1 + )uN(14)数值计算得到的结果符合定性分析的结论，但 是有更多细节。下面是装置在一个周期内(0 2 n/3, n/3)各参数变化的情况。图6表明了尾长与各点打滑时长的关系： (1)尾巴越长，A点打滑时间越多；存在一个临界 尾长a*，ai ai时A点始终打滑。(2)大部分情况 下装置不会出现A和B都打滑的情况，只在临界尾 长附近才会出现，且存在的时间都很

16、短，在实际中不 容易观察到。-0.041in-质心妇-*-质心如字向加速度 亍向加速度0.020-0.020.51.51.0t/s0图7尾长与质心加速度的关系2.0为什么在临界尾长垢附近才会出现A和B均 打滑的情况？从图8和图9中可以看出，当A和B 压力曲线有交点时，会出现模式转换。当压力曲线交 点区域较窄时（斜率大），只在模式一和模式二间进 行转换，不会出现模式三；当压力曲线交点区域较宽 时（斜率小），三种模式都会出现。这可以理解为：临 界尾长由的装置在t ! 0时，A和B压力就很接 近，此时质心速度很小，数值计算时不能一步跨出这 个区域，就会出现A和B都打滑这种模式。在非临 界尾长情况下，

17、若A和B压力接近，此时质心速度 较大，数值计算时可以一步跨出这个区域，不会出现 A和B都打滑的状况。表1中的数据支持了 A和B都打滑的必要条 件：刚开始的时刻，且两者的压力、摩擦力都很接 近，以及在这种状况下如果按模式一或模式二计算, 都会出现摩擦关系违约的情况= 36 mm，t = 0.0814 s，3 = 53.9394。）。图10表明一个周期内同一尾长A点的压力随 时间变小，因为时间增加时B点更靠近重心；而同 一时刻尾巴越长A点压力越小。图11表明通过判 断摩擦是否违约后再选择适当的模式，不同尾长时 A点摩擦关系都没有出现违约情况，B点也类似。图8两种模式的情况1.21.00.80.60

18、.40.2000.51.01.5t/s图9三种模式的情况FaNaFbNbFa /NaFb /Nb模式一 B打滑0.152 20.493 80.152 30.507 80.308 2 x0.299 9 P模式二A打滑0.14780.492 60.151 80.504 80.300 0 P0.300 7 x模式三A3都打滑0.14780.492 60.151 50.504 90.300 0 P0.300 0 P表1不同模式计算出的摩擦力与压力0.35NRMISRBM图10不同尾长时Fa和Na0.300.25 0.20田 0.150.100.050/7一一/一一一无尾：0粗隹:U mm中尾:25 mm长尾：50JmmL00.51.0t/s1.5图11不同尾长时Fa