25 转矩（力矩）
在上一章中，我们模拟了一个角速度恒定的场景。在本章中，我们将使其变得更加复杂；我们将模拟一个茶壶，在转台上以恒定的角加速度和减速度旋转。
25.1 角加速度
就像线加速度是速度的导数一样，角加速度也是角速度的导数。另外，类比于线加速度是由力引起的，角加速度也是由力、力矩的转动引起的。如果您不熟悉力矩，可以阅读http://modsimpy.com/torque进行了解。
一般来说，力矩是一个向量，定义为和
的叉积，其中
是杠杆臂，表示从旋转点到施力点的向量，
是表示力的大小和方向的向量。
然而，对于本章中的问题，我们只关注力矩的大小，而不关心方向。在这种情况下，我们可以运用公式：
其中，是力矩，
是杠杆臂的长度，
是力的大小，
是
和
之间的夹角。
由于力矩是长度和力的乘积，所以用牛米（N m）表示。
25.2
与线性加速度通过牛顿第二运动定律与力相关的方式相同，，角加速度通过另一种形式的牛顿定律与力矩相关：
其中是角加速度，
是转动惯量。正如质量使物体难以加速一样[1]，转动惯量也使得物体难以旋转。
在最一般的情况下，三维物体绕任意轴旋转，转动惯量是张量，是一个以向量为参数并返回一个向量的函数。
幸运的是，在一个系统中，所有的旋转和力矩都发生在一个单轴上，我们不必处理最一般的情况，而可以把转动惯量看作标量。
对于质量为的小物体，在距离r处绕一点旋转，转动惯量为
。对于更复杂的物体，我们可以把物体分成多个小质量，计算每个质量的转动惯量，然后把它们相加。
然而，对于大多数简单的形状，人们已经做了计算，你可以直接查找结果。例如，可以参见http://modsimpy.com/moment。
25.3 茶壶和转盘
中餐馆的桌子通常有一个旋转的托盘或转盘，以方便顾客分享菜肴。这些转盘由低摩擦轴承支撑，使它们能够轻松转动和滑动。然而，它们可能很重，特别是当它们装满食物时，因此它们有很大的转动惯量。
假设我坐在一张桌子前，转盘上放着一壶茶，坐在我正对面的人让我递茶。我用1 N的力推转台的边缘，直到它转动了0.5 rad，然后放开，转盘滑动到距离起始位置的1.5 rad处停止。那么我应该施加多大的力才能使茶壶滑到我正对面的停止位置？
我们将按照以下步骤回答这个问题：
我将使用第一次推动的结果来估算转盘的摩擦系数。
作为练习，您将使用该摩擦系数来估计旋转转盘通过剩余角度所需的力。
我们的模拟将使用以下参数：
转台半径为0.5 m，重量为7 kg。
茶壶重0.3 kg，离转盘中心0.4 m。
图25.1显示了这种情况，其中是施加在转台周边的力
，垂直于力臂
，
是产生的力矩。靠近底部的蓝色圆圈是茶壶。
图25.1：带茶壶的转盘示意图
下面是一个带有场景参数的Params对象：
params = Params(radius_disk=0.5*m,
mass_disk=7*kg,
radius_pot=0.4*m,
mass_pot=0.3*kg,
force=1*N,
torque_friction=0.2*N*m,
theta_end=0.5*radian,
t_end=20*s)
make_system创建初始状态init，并计算转盘和茶壶的总转动惯量。
def make_system(params):
mass_disk, mass_pot = params.mass_disk, params.mass_pot
radius_disk, radius_pot = params.radius_disk, params.radius_pot
init = State(theta=0*radian, omega=0*radian/s)
I_disk = mass_disk * radius_disk**2 / 2
I_pot = mass_pot * radius_pot**2
return System(params, init=init, I=I_disk+I_pot)
在初始状态下，theta代表桌子的角度，单位为rad；omega表示角速度，单位为rad/s。
I_disk是转盘的转动惯量，它是基于水平圆盘绕一个通过其中心的垂直轴旋转的转动惯量：
I_pot是将其视为一个点质量的茶壶的转动惯量：
在国际单位制中，转动惯量用kg m2表示。
现在我们可以创建一个 System对象了：
system1 = make_system(params)
这是一个采用当前状态的斜率，包含角度和角速度，并返回导数，角速度和角加速度：
def slope_func(state, t, system):
theta, omega = state
radius_disk, force = system.radius_disk, system.force
torque_friction, I = system.torque_friction, system.I
torque = radius_disk * force - torque_friction
alpha = torque / I
return omega, alpha
在这种情况下，施加在转台上的力始终垂直于杠杆臂，因此，力产生的力矩为
。
torque_friction表示摩擦产生的力矩。因为转台是在theta的正方向上旋转，所以摩擦力在theta的负方向上作用。
现在我们已经准备好了进行模拟，但我们必须先解决一个问题。
当我停止推动转盘时，角加速度突然改变。我们可以用if语句来实现斜率函数，该语句检查theta的值并相应地设置力。对于这样一个粗糙的模型，这可能是好的。但是，如果我们分两个阶段进行系统模拟，我们将得到更准确的结果：
在第一阶段，力是恒定的，我们运行到
theta为0.5 rad。在第二阶段，力为0，我们一直运行到
omega为0。
然后我们可以将两个阶段的结果合并到一个单一的时间框架中。
这是我将在第一阶段中使用的事件函数，它在theta到达 theta_end端时停止模拟，也就是我停止推动时：
def event_func1(state, t, system):
theta, omega = state
return theta - system.theta_end
现在我们可以开始运行第一阶段了。
results1, details1 = run_ode_solver(system1, slope_func,
events=event_func1)
在运行第二阶段之前，我们必须提取第一阶段的最终时间和状态。
t_0 = get_last_label(results1) * s
init2 = results1.last_row()
现在我们可以为第二阶段创建一个系统对象，初始状态来自第一阶段，力为0。
system2 = System(system1, t_0=t_0, init=init2, force=0*N)
对于第二阶段，我们需要一个事件函数，它在转盘停止时停止；也就是说，当角速度为0时。
def event_func2(state, t, system):
theta, omega = state
return omega
现在我们可以运行第二阶段了。
results2, details2 = run_ode_solver(system2, slope_func,
events=event_func2)
Pandas提供可以结合results1和results2的combine_first。
results = results1.combine_first(results2)
图25.2：施加力和摩擦力时转台的角度和角速度
图25.2显示了结果。角速度omega在我推的时候线性增大，在我放开后线性减小。角theta是角速度的积分，所以它在每个相位形成一条抛物线。
在下一节中，我们将使用这个模型来估计由于摩擦而产生的力矩。
25.4 估算摩擦
让我们从上一节中获取代码并将其包装在一个函数中。
def run_two_phases(force, torque_friction, params):
# put the parameters into the Params object
params = Params(params, force=force,
torque_friction=torque_friction)
# run phase 1
system1 = make_system(params)
results1, details1 = run_ode_solver(system1, slope_func,
events=event_func1)
# get the final state from phase 1
t_0 = results1.last_label() * s
init2 = results1.last_row()
# run phase 2
system2 = System(system1, t_0=t_0, init=init2, force=0*N)
results2, details2 = run_ode_solver(system2, slope_func,
events=event_func2)
# combine and return the results
results = results1.combine_first(results2)
return TimeFrame(results)
我们可以用run_two_phases来写一个误差函数，用root_bisect来求出由于摩擦产生的力矩，这个力矩产生了第一次推动的观察结果，总共旋转了1.5 rad。
def error_func1(torque_friction, params):
force = 1 * N
results = run_two_phases(force, torque_friction, params)
theta_final = results.last_row().theta
print(torque_friction, theta_final)
return theta_final - 1.5 * radian
现在我们可以使用root_bisect来估计由于摩擦产生的力矩了。
res = root_bisect(error_func1, [0.5, 2], params)
force = res.root
结果是0.166 N m，比最初的猜测要小一些。
现在我们知道了由于摩擦产生的力矩，我们可以计算出旋转转盘通过剩余角度所需的力，即从1.5 rad到3.14 rad。
在本章的笔记本“chap25.ipynb”中，您将有机会完成练习。有关下载和运行代码的说明，请参阅0.4节。
本书的中文翻译由南开大学医学院智能医学工程专业2018级、2019级的师生完成，方便后续学生学习《Python仿真建模》课程。翻译人员（排名不分前后）：薛淏源、金钰、张雯、张莹睿、赵子雨、李翀、慕振墺、许靖云、李文硕、尹瀛寰、沈纪辰、迪力木拉、樊旭波、商嘉文、赵旭、连煦、杨永新、樊一诺、刘志鑫、彭子豪、马碧婷、吴晓玲、常智星、陈俊帆、高胜寒、韩志恒、刘天翔、张艺潇、刘畅。
整理校订由刘畅完成，如果您发现有翻译不当或者错误，请邮件联系changliu@nankai.edu.cn。
本书中文版不用于商业用途，供大家自由使用。
未经允许，请勿转载。
