24 旋转
在本章中，我们对涉及旋转物体的系统进行建模。一般来说，旋转是复杂的：在三维空间中，物体可以绕三个轴旋转；物体通常比其他物体更容易绕某些轴旋转；当它们绕着某些轴旋转时，它们可能是稳定的，而不是其他轴。
如果一个物体的结构随着时间的推移而改变，它可能会变得更容易或更难旋转，这就解释了体操运动员、潜水员、溜冰运动员等令人惊讶的动态。
当你对旋转物体施加扭转力时，效果往往与直觉相反。例如，请看这段关于陀螺进动的视频http://modsimpy.com/precess。
在这一章中，我们将不讨论旋转的物理学。相反，我们将集中在简单的场景中，所有的旋转和扭转力都围绕着一个轴。在这种情况下，我们可以把一些向量量当作标量来处理（就像我们有时把速度当作一个隐含方向的标量一样）。
这种方法使模拟和分析许多有趣的系统成为可能，但是您也会遇到使用更通用的工具箱更好地处理的系统。
本章和下一章的基本概念是角速度、角加速度、力矩和转动惯量。如果你还没有
图24.1：一卷卫生纸的示意图，显示了由于小旋转而导致的纸张长度变化，dθ。
熟悉这些概念，接下来我将对它们进行定义，并指出其他阅读资料。
在下一章的最后，您将使用这些工具来模拟溜溜球的行为（参见http://modsimpy.com/yoyo). 不过，我们会慢慢来的，从厕纸开始。
24.1卫生纸的物理特性
作为一个带旋转的系统的简单示例，我们将模拟一卷卫生纸的制造。从中心的纸板管开始，我们将卷起47米长的纸，这是美国一卷厕纸的典型长度（见http://modsimpy.com/paper).
图24.1显示了一个系统图：r代表一个时间点的滚动半径。最初，r是纸板芯的半径，Rmin。当滚动完成时，r是Rmax。
我将使用θ来表示横摇的总旋转弧度。在图中，dθ表示θ的微小增加，它对应于沿r dθ的滚动圆周的距离。
最后，我将使用y来表示卷纸的总长度。最初，θ=0且y=0。θ每增加一小部分，y就会相应增加：
如果我们把两边除以时间的一个小增量，dt，我们就得到了y随时间变化的微分方程。
当我们卷起纸张时，r也增加了。假设r每转增加一个固定的量，我们可以写
其中k是一个未知的常数，我们要弄清楚。同样，我们可以将两边除以dt得到一个时间微分方程：
最后，假设θ以10rad/s的恒定速率增加（大约每分钟95转）：
这种变化率称为角速度。现在我们有一个由三个微分方程组成的系统，我们可以用来模拟这个系统。
24.2实施
在这一点上，我们有一个相当标准的过程来编写这样的模拟。首先，我们会从品脱中得到我们需要的单位：
radian = UNITS.radian
m = UNITS.meter
s = UNITS.second
并使用系统参数创建一个Params对象：
params = Params(Rmin = 0.02 * m,
Rmax = 0.055 * m,
L = 47 * m,
omega = 10 * radian / s,
t_end = 130 * s,
dt = 1*s)
Rmin和Rmax是半径的初始值和最终值，r.L是纸张的总长度。t_end是模拟的时间长度，dt是ODE解算器的时间步长。
我们使用Params对象生成一个系统对象：
def make_system(params):
init = State(theta = 0 * radian,
y = 0 * m,
r = params.Rmin)
k = estimate_k(params)
return System(params, init=init, k=k)
初始状态包含三个变量，θ、y和r。
estimate_k计算与θ和r相关的参数k。其工作原理如下：
def estimate_k(params):
Rmin, Rmax, L = params.Rmin, params.Rmax, params.L
Ravg = (Rmax + Rmin) / 2
Cavg = 2 * pi * Ravg
revs = L / Cavg
rads = 2 * pi * revs
k = (Rmax - Rmin) / rads
return k
Ravg是平均半径，介于Rmin和Rmax之间，因此Cavg是r为Ravg时的滚动周长。
revs是指如果r在Ravg下为常数，卷起长度L所需的总转数。rads只是转数转换成弧度。
最后，k是r每转一弧度的变化。对于这些参数，k约为2.8e-5 m/rad。
现在我们可以使用第24.1节中的微分方程来编写斜率函数：
def slope_func(state, t, system):
theta, y, r = state
k, omega = system.k, system.omega
dydt = r * omega
drdt = k * omega
return omega, dydt, drdt
与往常一样，slope函数接受一个状态对象、一个时间和一个系统对象。State对象包含时间t的theta、y和r的假设值。
斜率函数的作用是计算这些值的时间导数。θ通常表示为角速度的导数。
当纸卷上的纸张长度为L时，我们希望停止模拟。我们可以使用一个事件函数，当y等于L时通过0：
def event_func(state, t, system):
theta, y, r = state
return y - system.L
现在我们可以这样运行模拟：
results, details = run_ode_solver(system, slope_func,
events=event_func)
图24.2显示了结果。θ随时间线性增长，正如我们所预期的。因此，r也呈线性增长。但是由于y的导数依赖于r，并且r在增加，所以y随着斜率的增大而增大。
图24.2：纸卷模拟结果，显示随时间变化的旋转、长度和半径。
因为这个系统非常简单，所以模拟它几乎是愚蠢的。我们将在下一节中看到，用解析方法求解微分方程是很容易的。但从一个简单的模拟开始，作为探索和检验假设的一种方法，通常是有用的。
为了使模拟工作，我们必须正确地选择单元，这有助于及早发现概念错误。通过插入真实的参数，我们可以检测出导致不切实际结果的错误。例如，在该系统中，我们可以检查：
模拟的总时间约为2分钟，这似乎是合理的时间，它将需要滚动47米的纸。
θ的最终值约为1250rad，相当于大约200转，这似乎也是合理的。
r的初始值和最终值与Rmin和Rmax一致，正如我们在选择k时所预期的那样。
但是现在我们有了一个工作模拟，做一些分析也是有用的。
24.3分析
第24.1节中的微分方程非常简单，我们可以直接求解。由于角速度恒定：
我们可以通过对两边的积分来找到θ作为时间的函数：
当初始条件θ（0）=0时，我们得到C1=0。同样，
所以
当初始条件r（0）=Rmin时，我们发现C2=Rmin。然后我们可以把r的解代入y的方程中：
整合双方的利益：
$$
y(t) =
kωt2
/2 + Rmint
ω + C3
$$
所以y是一条抛物线，你可能已经猜到了。当初始条件y（0）=0时，我们得到C3=0。
我们也可以用这些方程来找出y和r之间的关系，与时间无关，我们可以用它来计算k。利用14.3节中的移动，我将方程24.1和24.2分开，得到
分离变量产生
整合双方的收益
当y=0时，r=Rmin，所以
解y，我们有
当y=L，r=Rmax时，代入这些值得到
求k产量
插入参数值得到2.8e-5 m/rad，与我们在第24.2节中计算的“估计值”相同。在这种情况下，估计结果是准确的。
在你继续之前，你可能想读一下本章的笔记本，chap24.ipynb，并练习一下。有关下载和运行代码的说明，请参阅0.4节。
本书的中文翻译由南开大学医学院智能医学工程专业2018级、2019级的师生完成，方便后续学生学习《Python仿真建模》课程。翻译人员（排名不分前后）：薛淏源、金钰、张雯、张莹睿、赵子雨、李翀、慕振墺、许靖云、李文硕、尹瀛寰、沈纪辰、迪力木拉、樊旭波、商嘉文、赵旭、连煦、杨永新、樊一诺、刘志鑫、彭子豪、马碧婷、吴晓玲、常智星、陈俊帆、高胜寒、韩志恒、刘天翔、张艺潇、刘畅。
整理校订由刘畅完成，如果您发现有翻译不当或者错误，请邮件联系changliu@nankai.edu.cn。
本书中文版不用于商业用途，供大家自由使用。
未经允许，请勿转载。
