第 43 卷 第 5 期 力 学 与 实 践 2021 年 10 月 Matlab 求解理论力学问题系列 (四) 乒乓球的滚动 高云峰1) (清华大学航天航空学院, 北京 100084) 摘要 本文讨论了乒乓球运动的理论解和数值解。在考虑摩擦、滚动摩阻、空气阻力等条件的一般情况下， 只能获得部分解析解，而采用数值计算的方法可获得系统所有的运动规律。文章介绍了如何利用 Matlab 处理 分段函数问题，包括求解微分代数方程、利用积分中断条件确定分段的时间点，以及利用部分解析解验证数值 计算的准确性等，最后通过大量计算获得乒乓球向前运动及向后运动的分界线，以及乒乓球正好回到出发位置 的条件。 关键词 微分代数方程, 解析解, 数值解, 积分中断 中图分类号: O313.3 文献标识码: A doi: 10.6052/1000-0879-20-240 THE ROLL OF A TABLE TENNIS GAO Yunfeng1) (School of Aerospace Engineering, Tsinghua University, Beijing 100084, China) Abstract Theoretical and numerical solutions of the kinematics of table tennis are discussed in this paper. Under the general conditions of friction (sliding, rolling friction and air resistance), analytical solutions can be obtained only in limited conditions, but numerical solutions are universally available. This paper handles various numerical problems in solving the kinematics with Matlab, including using piecewise function，solving differential-algebraic equation, obtaining the interruption point of the integral equation, and verifying the accuracy of the numerical results through comparing with available analytical solutions. At last, two critical conditions are investigated, one for the determination of forward or backward motion of the ball, the other for the ball jumping back to the starting position. Key words differential-algebraic equation, analytical solution, numerical solution, integral interruption 可能很多人尝试过用手指按压乒乓球，在一定 条件下可以使它向前运动，然后又会返回 (图 1)。本 文从理论分析和数值计算的角度，对这一现象进行 分析。 乒乓球滚动问题的特点是，乒乓球可能打滑也 可能不打滑，摩擦力或滚动摩阻的方向也可能发生 变化，因此需要在不同阶段列写不同的方程再进行 处理。 如果同时考虑空气阻力或滚动摩阻，不能获得 所有的解析解，这也是本系列文章的一个目的：绝大 部分动力学问题没有解析解，但是可以利用数值计 算的方法获得系统所有参数变化的规律，加深对问 题的理解。 1 乒乓球滚动的动力学方程及部分解析解 假设乒乓球质量为 m，半径为 r，水平面摩擦 系数为 µ，滚动摩阻为 δ，空气阻力系数为 n。乒乓 球的运动和受力如图 2 所示，其中 xC 为乒乓球质心 2020–06–01 收到第 1 稿，2020–07–01 收到修改稿。 1) E-mail: gaoyunfeng@tsinghua.edu.cn 引用格式: 高云峰. Matlab 求解理论力学问题系列 (四) 乒乓球的滚动. 力学与实践, 2021, 43(5): 753-759 Gao Yunfeng. The roll of a table tennis. Mechanics in Engineering, 2021, 43(5): 753-759
754 力 学 与 实 践 2021 年 第 43 卷 ᤹࣋⭘ 图 1 乒乓球游戏 xC xC vA Mf θ θ . . y O C A mg N F C A x R 图 2 乒乓球运动和受力图 的运动，θ 为乒乓球的转角，N 为压力，F 为摩擦 力，Mf = δN 为滚动摩阻，R = nx˙ C 为空气阻力。 根据实际经验，乒乓球在初始阶段会打滑，后来 会进入纯滚动状态，因此下面分两阶段来处理。第一 阶段，乒乓球打滑前进 (图 2)，可以列出动力学方程 及补充的摩擦关系，为 mx¨C = −F − nx˙ C 0 = N − mg JC ¨θ = −δN + F r F = µN (1) 其中 JC = 2mr2/3 是薄壁圆球对质心的转动惯量。 令 β = n/(2m)，方程 (1) 中的未知力消去后可以得 到系统的运动微分方程，为 x¨C + 2βx˙ + µg = 0 ¨θ − 3 2 (µr − δ)g/r2 = 0 (2) 方程 (2) 存在解析解，注意 Matlab 除了可以数 值计算，还可以进行符号推导，包括对微分方程的积 分 [1]。利用 Matlab 对方程 (2) 进行积分，源代码见 图 3。 图 3 微分方程的积分 图 3 中 syms 为定义符号；diff(xc,t,2) 表示 xc 对时间的 2 阶导数；dsolve 表示对微分方程进行符 号积分。运行后屏幕上显示方程 (2) 的解为 xc=C1+C2*exp(-2*beta*t) - (g*mu*t)/(2*beta) theta=C4+C3*t-(tˆ2*(3*delta*g-3*g*mu*r))/ (4*rˆ2) 代入初始条件后有 xC = 1 4 ( gµ/β2 − 2 ˙x0/β) (1 − e −2βt)− 1 2 gµt/β θ = 3 4 (µr − δ)gt2/r2 + ˙θ0t (3) 这一阶段结束的标志是接触点 A 速度为零。A 点速 度根据图 2 可以分解为牵连速度 ve = ˙xC 和相对速 度 vr = r ˙θ，代入式 (3) 处理后有 ( x˙ 0 − 1 2 gµ/β) e −2βt − 3 2 (µr − δ)gt/r − 1 2 gµ/β − r ˙θ0 = 0 (4) 设 t ∗ 时刻触点 A 速度为零，在退化情况下方程 (4) 可以有解析解，但是在一般情况下没有 t ∗ 的解析表 达式 (这将导致后续分析出现问题)。不过在数值积 分中 t ∗ 可以利用积分中断获得，这表明了在一般情 况下利用数值计算的潜力。 第二阶段，乒乓球作纯滚动，但乒乓球可能向前 运动也可能向后运动，滚动摩阻的方向相应也应改 变，定义方向函数 λ = +1, if vC > 0 0, if vC = 0 −1, if vC < 0 (5)
第 5 期 高云峰：Matlab 求解理论力学问题系列 (四) 乒乓球的滚动 755 注意这一阶段 A 点速度为零，此时动力学方程及补 充运动学关系为 mx¨C = −F − nx˙ C 0 = N − mg JC ¨θ = −λδN + F r x¨C = r ¨θ (6) 方程 (6) 中消去未知的力，可以简化为 x¨C + 6 5 βx˙ C + 3 5 gλδ/r = 0 (7) 类似可以利用 Matlab 的符号推导得到解析解，为 xC = C1 + C2e −6βt/5 − λgt/(2βr) (8) 问题在于，由于第一阶段结束的时间 t ∗ 以及 xC (t ∗ ) 和 x˙ C (t ∗ ) 没有解析表达式，方程 (8) 中的积分常数 无法确定，因此无法得到更多的理论解。 由于在一般情况下很难得到动力学问题的全部 解析解，下面介绍数值计算。 2 微分方程的分段准备工作 首先把第一、二阶段的动力学方程 (1) 和 (6) 转 换为 AX = B 的形式，分别为 m 0 1 0 0 JC −r δ 0 0 0 1 0 0 1 −µ x¨C ¨θ F N = −nx˙ 0 mg 0 (9) m 0 1 0 0 JC −r λδ 0 0 0 1 1 −r 0 0 x¨C ¨θ F N = −nx˙ 0 mg 0 (10) 这是混合形式的微分代数方程。在数值计算中， 当前 t 时刻的 xC , θ, x˙ C , ˙θ 是已知量，而 x¨C , ¨θ, F, N 是未知量，因此先把 x¨C 和 ¨θ 视为代数量，调用 X =inv(A)*B 求出当前时刻的 x¨ ∗ C 和 ¨θ ∗ (求解后已 经是具体数值)，再把动力学方程转换为标准的一阶 方程组形式，设 y1 = xC , y2 = ˙xC , y3 = θ, y4 = ˙θ (11) 从而得到一阶微分方程组和初始条件为 y˙1 = y2 y˙2 = ¨x ∗ C y˙3 = y4 y˙4 = ¨θ ∗ , y1(0) = 0 y2(0) = ˙xC (0) y3(0) = 0 y4(0) = ˙θ(0) (12) 第一阶段和第二阶段的分界点时刻 t ∗ 可以由 Matlab 积分时的中断条件获得。在积分之前，先把 积分中具体的要求在 options 中说明，其格式为 options = odeset( ′RelTol′ ,number1, ′AbsTol′ , number2, ′ events′ , @event name); (13) 格式中黑正体表示固定格式，白斜体表示变量。注意 ‘events’，@ 表示在积分时会调用名为 event name 这个子程序，具体的中断条件可以在该子程序中表 达，在本问题中即 A 点速度为零。在设定积分条件 后，再调用积分函数，其格式为 [t, y] = ode45(@rg kt name, [start time : step time : end time], y0, options); (14) 格式中黑正体表示固定格式，白斜体表示变量。其 中 [t, y] 是积分的时间和结果；动力学方程在子程序 rg kt name 中描述；y 是初始条件；如果没有中断条 件，积分时按等步长 step time 从 start time 开始到 end time 结束，有中断条件后，最后一步的积分步 长会自动调整，结束时间也会自动提前。 在程序中先给参数赋值，然后分段计算，这一部 分的源代码见图 4。 图 4 带中断条件的分段积分
756 力 学 与 实 践 2021 年 第 43 卷 图 4 中 rg kt pingpang1，rg kt pingpang2 分别 根据方程 (9) 和方程 (10) 解出的结果进行积分；在 event pingpang1 中满足 A 点速度为零就中断积分； 第一阶段结束后把最后时刻的参数赋给 tf 和 new y0 作为第二阶段的初值，再进行第二阶段积分。第一阶 段具体中断的条件见图 5。 3 数值计算与解析解的比较 在正式进行大量计算前，首先对程序进行一些 测试。由于第一阶段有解析解，可以用于验证数值积 分。下面参数在所有计算中保持不变 m = 3 g, r = 2 cm, g = 9.8 m/s 2 µ = 0.3, JC = 2 3 mr2 , xC (0) = 0 x˙ C (0) = 1 m/s, θ(0) = 0 而 ˙θ0, β = n/(2m), δ 根据情况改变。 取 ˙θ0 = −50 rad/s, β = 0.1, δ = r/100，第一阶 段的位移和速度曲线见图 6 和图 7。数值解 (圆点) 与解析解 (位移用圆圈，半径与角度乘积用方框) 对 比表明，两者很吻合，只有最后一个圆圈或方框与圆 点不重合。原因是数值计算中运用了中断，最后一步 不是等步长的；而解析解的最后一步的时间没有解 析表达式，按等步长来处理。 图 6 和图 7 的结果表明数值计算的精度很好 (在 前面文章中已介绍精度 ∼ 10−10)，另一方面，把解 析解与数值解结合起来，可以取长补短，更好地处理 第二阶段的情况。 图 5 中断子程序 0 0.05 0.10 0.15 0.20 0.25 0.30 ᰦ䰤/s -0.15 -0.10 -0.05 0 0.05 0.10 0.15 ս〫/m xC〟࠶٬ xC⨶䇪٬ r *θ〟࠶٬ r *θ⨶䇪٬ 图 6 位移的数值解与解析解 0 0.05 0.10 0.15 0.20 0.25 0.30 ᰦ䰤/s -1.0 -0.5 0 0.5 1.0 䙏ᓖ m.s-1 dxC〟࠶٬ dxC⨶䇪٬ r *dθ〟࠶٬ r *dθ⨶䇪٬ 图 7 速度的数值解与解析解
第 5 期 高云峰：Matlab 求解理论力学问题系列 (四) 乒乓球的滚动 757 第二阶段也有解析解，但把数值计算与解析解 进行对比，发现第一阶段数值解与理论解吻合，第二 阶段差别很大 (图 8 和图 9)。原因是方程 (8) 中的参 数与第一阶段的中止时间有关，由于没有解析解，于 是按等步长来计算时间，导致解出的常数 C1 和 C2 有偏差 (表 1 给出了第一阶段终止时刻的理论解和 数值解)。 如果用第一阶段数值计算中断时的最后结果来 确定第二阶段的解析解中的初始条件，可以发现两 阶段的数值解 (圆点) 全部落在代表解析解的圆圈或 菱形中了 (图 10 和图 11)。这一结论既证明了数值 解的精确性，也表明了解析解与数值解结合的潜力。 0 1 2 3 0 0.05 0.10 0.15 0.20 0.25 0.30 ս〫/m ᰦ䰤/s xC〟࠶٬ ⨶䇪٬ ⨶䇪٬ 图 8 位移的数值解与解析解 0 0.5 1.0 1.5 2.0 2.5 3.0 -6 -4 -2 0 2 4 䖜䀂/(°) ᰦ䰤/s θ 〟࠶٬ θ ⨶䇪٬ θ ⨶䇪٬ 图 9 转角的数值解与解析解 表 1 积分常数及其他参数的变化 xC /m θ/rad x˙ C /(m·s−1 ) ˙θ/(rad·s−1 ) t ∗/s C1/m C2/m 解析解 0.161 5 0.085 7 −5.408 2 13.945 0.300 4.959 0 4.797 5 数值解 0.158 1 0.165 0 −5.704 8 8.247 8 0.273 3 5.616 1 5.458 0 0 1 2 3 0 0.05 0.10 0.15 0.20 0.25 0.30 ս〫/m ᰦ䰤/s xC〟࠶٬ ⨶䇪٬ ⨶䇪٬ 图 10 位移的数值解与修正的解析解 0 0 0.5 1.0 1.5 2.0 2.5 3.0 -6 -4 -2 2 4 䖜䀂/(°) ᰦ䰤/s θ 〟࠶٬ θ ⨶䇪٬ θ ⨶䇪٬ 图 11 转角的数值解与修正的解析解 有了前面的验证，计算的可靠性有了保证。 乒 乓球游戏有两种结果：模式 1 初始角速度小，乒乓 球只会前进不能后退；模式 2 初始角速度大，乒乓 球先前进再后退。 根据式 (3)，得到速度表达式为 x˙ C = 1 2 (gµ/β − 2 ˙x0) e−2βt − 1 2 gµ/β ˙θ = 3 2 (µr − δ)gt/r2 + ˙θ0 (15)
758 力 学 与 实 践 2021 年 第 43 卷 在第一阶段， 质心速度在小阻尼时， 利用 e −2βt ≈ 1 − 2βt，其斜率 k1 ≈ 2βx˙ 0 − gµ；牵连速 度 r ˙θ 的斜率是 k2 = (3/2)(µ − δ/r)g，在第二阶段， 质心速度的斜率 k3 = (3/5)gλδ/r。 但是当 A 点速度为 0 时 (牵连速度与相对速度 曲线相交)，乒乓球就开始作纯滚动，模式 1 中乒乓 球继续前进，模式 2 中乒乓球会返回 (图 12)；而在 速度图中，牵连速度与相对速度会一直相等且逐渐 变为 0 (静止状态)，模式 1 速度从正趋于 0 (一直 前进)，模式 2 速度从负趋于 0 (先前进再后退)，见 图 13。 0 1 2 3 4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 ⁑ᔿ1:অ㓟ࡽ䘋 ⁑ᔿ2:ࡽݸ䘋ਾ䘰 1:অ㓟 ࡽݸ:2 䘋ࡽ ⁑ᔿ2:ࡽݸ䘋ਾ䘰 ᰦ䰤/s 䐍/m 图 12 两种模式下质心位移曲线 0 1 2 1.0 0.5 0 -0.5 -1.0 -1.5 -2.0 3 4 ᰦ䰤/s 䍘ᗳ䙏ᓖ1 牵连速度1 䍘ᗳ䙏ᓖ2 牵连速度2 䙏ᓖ/(m.s-1 ) k k k 图 13 两种模式下速度曲线 如果关心角速度对最终位移的影响，可以把初 始速度比 (定义为 |r ˙θ0/x˙ C0|， 即牵连速度与质心 速度之比) 作为横坐标，进行大量计算得出曲线如 图 14 所示，存在一条分界线，分界线左侧 (角速度 较小) 乒乓球最终停在前方 (包括先前进再后退，但 停在出发点前方)，分界线右侧 (角速度较大) 乒乓球 最终会停在出发点的后方。如果希望乒乓球向前运 动后返回，正好返回原点，初始角速度是多少 (注意 方程 (8) 没有解析解)？Matlab 可以把图形放大，只 需要把图 14 中分界线与曲线的交点区域不断放大， 就能得到分界线处的速度比值是 1.770 665 387 040 5， 然后以这样的速度比值单独再算一次，就可以让乒 乓球先前进然后正好退回到原点 (图 15)。 0 1 2 3 4 ս〫/m 2 1 0 -1 -2 -3 -4 䙏ᓖ∄ 㓯⭼࠶ ਚ㜭ࡽ䘋 Ⲵ४ฏ ਟԕਾ䘰 Ⲵ४ฏ 图 14 不同速度比对最终位移的影响 0 1 2 3 0 0.05 0.10 0.15 0.20 ↓ྭ䘄എ⛩ ᰦ䰤/s 䐍/m 图 15 正好退回原点的情况 4 小结 本文讨论了乒乓球运动的理论解和数值解。在 退化情况下，可以获得全部解析解，而在一般情况下 (如同时考虑摩擦、滚动摩阻、空气阻力)，只能得到 很少的解析解，而采用数值计算的方法可获得一般 情况下所关心的所有参数。 (1) 乒乓球的运动分不同的阶段，本文介绍了在
第 5 期 高云峰：Matlab 求解理论力学问题系列 (四) 乒乓球的滚动 759 Matlab (R2016b 版本) 中如何处理分段计算的方法， 以及微分代数方程求解的一般过程。 (2) 把解析解与数值解相结合，证明数值精度很 高。 (3) 适当运用 Matlab 中的符号推导，可以获得 部分解析解表达式 (包括部分微分方程的解析解)。 (4) 通过计算，得到了一般情况下乒乓球运动的 规律，例如在初始速度一定时，找到了初始角速度和 质心最终位移的关系，得到了乒乓球质心最终是否 前进的分界线的解，也可以让乒乓球正好返回原处。 这些结论在一般情况下都没有解析解。 参 考 文 献 1 甘勤涛, 胡仁喜, 程政田等. Matlab 数学计算与工程分析从入 门到精通. 北京: 机械工业出版社, 2019 (责任编辑: 胡 漫)