21 空气阻力
在上一章中，我们模拟了一枚硬币在真空中下落，也就是说，没有空气阻力。但是我们使用的计算框架非常复杂。一般情况下，很容易添加额外的力，包括阻力。在本章中，我提出了一个阻力模型，并将其添加到模型中。
21.1阻力
当一个物体通过流体时，例如空气，这个物体会对空气施加力，根据牛顿第三运动定律，空气会对这个物体施加一个大小相等、方向相反的力。 (参考 http://modsimpy.com/newton)这个阻力的方向与运动的方向相反，它的大小由阻力方程给出。 (参考 http://modsimpy.com/drageq)
Fd为阻力，单位为牛顿
ρ是流体的密度，单位是kg/m3
v是速度的大小，单位是m/s
A是物体的参考面积，单位是平方米，在这种情况下，参考区域是投射的正面区域，即从物体运动线上的一点(或更远的地方)所看到的物体的可见区域。
Cd是阻力系数，这是一个无量纲的量，取决于物体的形状(包括长度但不包括前缘面积)、表面性质以及它与流体的相互作用方式。
对于在空气中以中等速度移动的物体，典型的阻力系数在0.1到1.0之间，钝体物体在范围的高端，流线型物体在低端。(参考 http://modsimpy.com/dragco)
对于简单的几何物体，我们有时可以以合理的精度猜测阻力系数;对于更复杂的物体，我们通常需要进行测量并根据数据估计Cd。然，阻力方程本身就是一个模型，基于这样的假设Cd不依赖于方程中的其他项:密度、速度和面积。对于在空中以中等速度移动的物体(低于45英里/小时或20米/秒)，这个模型可能已经足够好了，但我们应该记住重新考虑这个假设。对于下落的硬币，我们可以用测量值来估计Cd。特别是，我们可以测量终端速度，它是阻力等于重力时候的速度:
m是物体的质量g是重力加速度。求解Cd产率方程:
据《流言终结者》报道，一枚硬币的终极速度在35到65英里每小时(见http://modsimpy.com/mythbust)。如果以40英里/小时(约18米/秒)为下限，Cd的估计值为0.44，接近光滑球体的阻力系数。现在我们准备为模型添加空气阻力。
21.2阻力添加实现
随着系统参数数量的增加，以及我们需要做更多的工作来计算它们，我们将发现定义一个Params对象来包含创建系统对象所需的数量是很有用的。Params对象类似于System和State对象;事实上，这三种能力都是一样的。我给它们起了不同的名字，以记录它们扮演的不同角色。
下面是掉落硬币的Params对象:
parems = Params(height = 381 *m,
v_init = 0 * m / s,
g = 9.8 * m / s,
mass = 2.5e-3 * kg,
diameter = 19e-3 * m,
rho = 1.2 * kg/m**3,
v_term = 18 * m / s)
质量和直径来自http://modsimpy.com/penny。空气的密度取决于温度、气压(取决于海拔高度)、湿度和成分(http://modsimpy.com/density)。我选择的值可能是典型的波士顿，马萨诸塞州20◦C。
下面是make_system的一个版本，它接受一个Params对象并返回一个系统:
def make_system(params):
diameter, mass = params.diameter, params.mass
g, rho = params.g, params.rho,
v_init, v_term = params.v_init, params.v_term
height = params.height
area = np.pi * (diameter/2)**2
C_d = 2*mass*g / (rho*area*v_term**2)
init = State(y=height, v=v_init)
t_end = 30*s
dt = t_end / 100
return System(params, area=area, C_d = C_d,
init = init, t_end = t_end, dt = dt)
System的第一个参数是params，因此结果包含params、area、C_d和其他参数中的所有参数。
我们需要Params对象的原因可能并不明显，但它们很快就会被证明是有用的。我们可以这样创建一个系统
system = make_system(params)
这是一个包含阻力的斜率函数
def slop_func(state, t, system):
y, v = state
rho, C_d, g = system.rho, system.C_d, system.g
area, mass = system.area, system.mass
f_drag = rho * v**2 * C_d * area / 2
a_drag = f_drag / mass
dydt = v
dvdt = -g + a_drag
return dydt, dvdt
据阻力方程和牛顿第二定律，f_drag是由于阻力而产生的加速度。为了计算总加速度，我们添加了重力和拖动的加速度。g是负的，因为它在y减小的方向上，a_drag是正的，因为它在y增大的方向上，在下一章我们会用到矢量对象来跟踪力的方向，并以一种不太容易出错的方式将它们相加。为了在硬币撞到人行道时停止模拟，我们将使用第20.3节中的事件函数:
def event_func(state, t, system):
y, v = state
return y
现在我们可以像这样运行模拟:
results, details = run_ode_solver(system, slope_func, events = event_func)
图21.1显示了结果。硬币加速到终极速度只需要几秒钟;之后，速度是常数，所以高度作为时间的函数是一条直线。
在你继续讲下去之前，你可能想看看这一章的笔记本，第21章。ipynb，做练习。有关下载和运行代码的说明，请参见第四部分。
本书的中文翻译由南开大学医学院智能医学工程专业2018级、2019级的师生完成，方便后续学生学习《Python仿真建模》课程。翻译人员（排名不分前后）：薛淏源、金钰、张雯、张莹睿、赵子雨、李翀、慕振墺、许靖云、李文硕、尹瀛寰、沈纪辰、迪力木拉、樊旭波、商嘉文、赵旭、连煦、杨永新、樊一诺、刘志鑫、彭子豪、马碧婷、吴晓玲、常智星、陈俊帆、高胜寒、韩志恒、刘天翔、张艺潇、刘畅。
整理校订由刘畅完成，如果您发现有翻译不当或者错误，请邮件联系changliu@nankai.edu.cn。
本书中文版不用于商业用途，供大家自由使用。
未经允许，请勿转载。
