本文通过Euler方法和四阶Runge-Kutta方法，研究了忽略和考虑空气阻力情况下加农炮弹的运动轨迹。分析了发射角对轨迹、水平射程的影响，探讨了数值解的稳定性和收敛性，以及空气阻力对运动轨迹的显著影响。结果表明，空气阻力可降低射程，且改变最大射程对应的发射角。
摘要
加农炮在军事上有着广泛的应用，研究加农炮弹出膛后的运动轨迹具有重要意义。炮弹运动的实质是抛体运动。抛体运动是我们熟知的经典问题，在不计及空气阻力的情况下，其运动方程十分简单，可以解析求解，但在空气阻力不可忽略的情况下，抛体运动方程将变得复杂，一般只能数值求解。本文利用Euler方法和Runge-Kutta方法对忽略空气阻力和计及空气阻力两种情况下加农炮弹的运动轨迹进行了求解，研究了炮弹在不同发射角下的炮弹运动轨迹，考察了数值计算中时间步长对计算结果的影响。我们还研究了炮弹水平射程对发射角的依赖关系，并估计了达到最大水平射程时的发射角。进一步，我们研究了在考虑阻力情况下的中靶问题，给出了对空中目标实现精确打击所需发射角的计算。
研究背景
抛体运动
加农炮弹出膛后的运动本质上是抛体运动，为便于处理，忽略炮弹在飞行过程中的转动和形变，将其视为质点。质点的抛体运动是我们所熟知的，可以利用牛顿力学对其进行描述。在不计及空气阻力的情况下，质点只受到重力作用，因而其运动方程由牛顿运动定律给出
其中
为炮弹的质量，
为重力。将上式写成分量形式
其中
为重力加速度。这是一组容易解析求解的微分方程，若给定初始条件，炮弹的运动情况可以完全地确定。但在考虑空气阻力的情况下，上述方程右侧需要增加与阻力相关的项，这将使得炮弹的实际运动情况变得复杂。可以将炮弹的运动方程一般性地写成如下形式
其中
分别为单位质量的物体在
方向和
方向受到的除重力以外的力，
和
的形式与具体情况有关。我们仅考察空气阻力，并考虑简单的阻力形式
其中
为阻力系数，是不显著依赖于物体性质和状态的常数。在许多情况下，仅考虑空气阻力对速度二次方的依赖项，从而简化模型，得到描述炮弹运动的方程
上述方程结合一定的初始条件可以完全确定炮弹的运动状态，但方程
的解析求解是困难的，因此我们将采用数值方法给出其数值解。
抛体运动方程的数值解法
- Euler方法
我们利用Euler方法来求解上面的方程，将两个二阶微分方程降阶，得到两组一阶微分方程组，按照Euler法的思想，将微分近似为差分，这里采用向前差分，得到方程
的显式Euler法计算格式
特别地，当
时，上式退化为忽略空气阻力的情形。由上述计算格式可以看出，只要给定初始条件
，即可逐步递推后续各时刻抛体的位置、速度，也就是说，只要知道前一时刻
抛体的位置、速度，即可求得下一时刻
抛体的位置、速度，如此反复迭代，最终可以得到一系列时刻抛体的运动情况。由于微分方程的数值解法是存在误差的，上述显式Euler法的整体截断误差与时间步长
相关，因此适当缩小步长可以提高计算精度，但这往往会以增加计算时间作为代价，因此需要权衡两者的利弊，选定合适的时间步长进行数值求解，后面我们将会考察步长对计算结果的影响。
- Runge-Kutta方法
作为对比，我们采用四阶Runge-Kutta方法对方程
进行求解，同样地，先将二阶微分方程降阶，得到一阶微分方程组
为便于表述，上面并未写出微分方程组的具体形式。按照Runge-Kutta法的基本思想，得到计算格式
将
的具体形式代入上式，即可得到具体的计算格式表达式，当
时，上式退化为忽略空气阻力的情形。同样地，在给定初始条件
后，可以逐步递推得到后续各时刻
抛体的位置、速度，从而确定抛体的运动情况。四阶Runge-Kutta法的整体截断误差与时间步长
相关，是比Euler法计算精度更高的一种方法，适当缩小时间步长同样可以进一步提高该方法的精度，因此选择合适的步长是很重要的，后面也将考察步长对计算结果的影响。
加农炮弹运动的研究
运动方程的解
在忽略空气阻力的情况下，炮弹的运动方程如式
所示。给定初始条件：炮弹的初始位置
，初始速度
，通过分离变量法，可以得到解析解
其中，
是初始速度
在
方向和
方向的分量。式
表明炮弹的运动轨迹是一抛物线，这是我们所熟知的。对于方程
，我们也可以采用数值方法进行求解，通过Matlab编程，得到计算结果如
在考虑空气阻力的情况下，炮弹的运动状态变得复杂，如前所述，解析地求解运动方程
变得十分困难，这时适合采用数值计算的方法给出方程的数值解。这里我们分别利用Euler法和Runge-Kutta法求出炮弹运动方程的数值解。选取发射角
，绘制炮弹的运动轨迹如
从图中可以看到，在忽略空气阻力的情况下，数值解的各离散点较好地符合于式
所给出的解析解，可以明显地看出，Runge-Kutta法的精度要优于Euler法，前者的计算结果全部落在解析解的图线上，而通过Euler法得到的解在自变量
增大时对解析解存在小的偏离。总的来说，上述结果验证了我们采用的数值方法的正确性和可行性，这为我们在考虑阻力情况下的数值计算提供了保障。
现在我们对有无空气阻力这两种情况下炮弹运动状态的差别作一简要分析。在图1中，我们选定了发射角为
，初始速度为
的情形进行考察。可以看到，在忽略阻力时，炮弹的抛体运动遵从式
所给出的运动规律，其运动轨迹是一标准的抛物线，关于最高点确定的竖直线左右对称；在考虑阻力时，炮弹的运动轨迹不再是标准的抛物线，没有对称的特性，并且炮弹的水平射程缩短，所能达到的最高点降低，体现出空气阻力对抛体运动状态的影响。关于这一问题，我们在后面还会详细介绍。
如前所述，我们初步验证了数值方法Euler法和Runge-Kutta法在求解微分方程时的可靠性，为了考察在不同发射角下炮弹的运动状态以及进一步验证上述两种数值计算方法的正确性，我们选取了四个发射角（分别为
和
）对忽略空气阻力情况下炮弹的运动轨迹进行了数值计算，并与该抛射角下的解析解进行了对比，通过Matlab编程，我们得到如
由图可知，通过Euler法和Runge-Kutta法得到的数值解都较好地符合于精确解，后者的精度更高。炮弹在忽略空气时的运动轨迹呈现出标准的抛物线，且水平射程随着抛射角的增大先逐渐增大，后逐渐减小，在抛射角
时达到最大射程，这是我们熟知的结果。
值得说明的是，这里我们选取
和
这四个角度仅作为代表，不具有任何特殊性。若希望考察其他抛射角度下炮弹的运动情况，可以运行笔者GitHub上的程序，并手动输入
到
之间的任意角度及其它参数，程序将给出该发射角下的运行结果。
数值解的稳定性与收敛性
我们在前面提到，利用Euler法或Runge-Kutta法求解微分方程得到的数值解是存在误差的。根据两种方法的计算格式，不难发现，解的误差与计算时所采用的时间步长相关，因此步长会在一定程度上影响数值解的好坏。过长的步长可能导致解偏离真实情况，而过短的步长会显著增加计算时间，因此选择一个合适的步长对数值计算是极为重要的。
为了研究数值计算的结果与时间步长的关系，我们选取一个初始步长，进行计算，之后逐步减小步长取值，考察数值解对时间步长的依赖性，我们期待当步长减小到某一较小值时，计算结果将不显著依赖于步长的选取。取步长分别为
，并选定抛射角
。利用Euler法和Runge-Kutt法求得考虑阻力情况下炮弹轨迹的数值解，如
从图中可以发现，随着时间步长的减小，炮弹运动轨迹的数值解发生形态上的变化，当步长取到
和
时，在这两个步长下计算得到的结果没有显著的差别。据此，我们做出如下论断：
- 随着时间步长的缩小，数值解的精度显著提升，没有出现振荡的情况，因此可以相信解是稳定的；
- 随着时间步长的缩小，Runge-Kutta法的计算结果较Euler法更快地收敛到高精度的数值解，体现出Runge-Kutta法的优越性；
- 当时间步长取到
时，继续减小步长已经不能显著提高解的精确性，但会增加计算耗费的时间。因此，我们断定
是一个合适的步长，在后面的计算中，我们将一律使用此步长。
空气阻力对抛体运动的影响
在许多情况下，我们在研究抛体运动时认为物体只受到重力作用，这称为理想抛体运动，而在实际情况下，物体还可能受到除重力以外的其他外力，情况要复杂得多。正如我们现在研究的加农炮弹出膛后的运动，炮弹在飞行过程中，除受到重力外，一般还受到空气阻力的作用，其运动方程为
其中
为空气阻力。在忽略阻力时，炮弹的运动规律由式
描述，这是熟知的结果，但在考虑阻力时，炮弹的运动将遵从方程
，其运动行为将由该方程的解给出。前面提到，方程
的解析求解是困难的，其解也可能是复杂的，一般地写为
考虑到上述原因，我们利用数值方法对其进行研究。
为了对比忽略空气阻力和考虑空气阻力两种情况下炮弹的运动状态，我们将数值方法得到的解绘制于同一幅图中，并选取了多个抛射角。通过Matlab编程，得到结果如图4所示。
从图中可以看到，空气阻力的存在会显著影响炮弹的运动状态。当发射角取不同值时，炮弹的轨迹虽有所差别，但都具有一致的本质特征，即相较于无阻力的情形，空气阻力的存在使得炮弹的运动轨迹不再是标准的抛物线，水平射程和最高点都有大幅度的降低，直观地看，空气阻力将使炮弹的水平射程大约减小为无阻力时的一半，因此在实际应用时，要提高炮弹的射程，一个有效的办法就是减小空气阻力，这可以通过对炮弹形状作出适当的设计而做到。此外，空气阻力还将影响炮弹达到最大射程时的发射角，后面的计算将表明，当空气阻力存在时，炮弹的最大射程并不是在发射角
时取到，这与无阻力时的情形不同。这意味着空气阻力显著影响了炮弹射程对发射角的依赖关系。
抛体运动的最大射程
我们知道，抛体在同一发射速度情况下，所能达到的最大射程与抛射角密切相关。由于我们给定了炮弹的初始速度
，若要得到最大射程，则需要研究炮弹运动轨迹与抛射角的关系。在忽略空气阻力的情况下，抛体在抛射角
时达到最大射程，这是我们熟知的结果；在考虑空气阻力的情况下，抛体受到与其速度二次方成正比的阻力作用，其最大射程将与发射角呈现非线性关系，使得求极值的过程变得非常困难，难以用显式形式给出最大射程及对应的射角。如前所述，炮弹在受到空气阻力作用时的运动方程为
其中
为阻尼系数。方程组
为强耦合方程组，一般而言不能求出轨迹方程的显式形式，需要采用数值方法求解。为了给出考虑空气阻力情况下炮弹的最大射程及对应的抛射角，我们将在得到数值解的基础上，近似求得炮弹的落地点，再通过搜索法给出达到最大水平射程时的抛射角。
由于数值计算的离散性，一般来说我们不能直接求得抛体高度刚好为零的点（即落地点），因此我们可以考虑采用近似的方法。正如前面提到的，在利用Euler法或Runge-Kutta法进行数值计算时，只要知道前一时刻
抛体的位置
抛体的位置
方向高度
式
表明落地点必定出现在
和
之间的轨迹上。我们可以利用落地点附近的数值解进行插值，求得插值多项式函数，进而通过求该函数的零点给出炮弹落地点的近似坐标。
考虑到线性插值的精度不高，我们采用Lagrange二次插值。选取落地点附近的三个数值解
，
和
，构造Lagrange基函数
进而得二次插值多项式函数
通过求该函数的零点即可得到落地点的近似坐标。我们利用上述方法从
开始逐渐改变抛射角（角增量取为
）并计算各抛射角下落地点的横坐标（即射程），顺序搜索，直至某次计算得到的射程小于上一次的射程，结束搜索，得到最大射程及对应的抛射角。根据这一思想，计算得到达到最大射程时炮弹的轨迹如
从图中可以看出，在经Euler法和Runge-Kutta法算得的数值解的基础上得到的最大射程和对应的抛射角存在一些微小的差别，在误差的允许范围内，我们姑且认为Runge-Kutta法得到的结果更接近真实情况，因此在考虑空气阻力的情况下，当炮弹的抛射角
时达到最大射程
。利用同样的方法，我们还可以算得忽略空气阻力时的最大射程为
，可见空气阻力的存在使得最大射程减小约一半左右。这一结果对实际情况具有指导性的意义，由于炮弹在飞行过程中受到阻力，而阻力会显著影响射程对抛射角的依赖关系，因此在实际中需要根据空气阻力的具体条件来确定发射角以使炮弹能飞行更远的距离。
炮弹对目标的精确打击
在实际应用中，人们利用炮弹打击的目标往往不是固定于地面的死靶，目标的竖直高度可能是任意的，因此不能只考虑打击同高度的地面目标的情况，而是希望对不同高度的目标予以精确打击。我们在这里考察一种简单的情况，即炮弹的发射速度固定为
，加农炮处于坐标系原点，假设目标位于
。下面我们将通过数值计算给出炮弹以何种发射角出射才能命中这一目标。
考虑存在空气阻力的情况，利用Euler法和Runge-Kutta法求解运动方程
。当发射角较大时，炮弹将能够打到
附近的位置，由于数值计算的离散性，一般来说我们不能直接得到
刚好为
的点（即目标点的
方向坐标），仿照前面求射程时的办法，我们采用Lagrange线性插值得到炮弹打到
时对应的竖直高度的近似值。利用上述方法从
开始逐渐改变抛射角（角增量取为
）并计算各抛射角下炮弹打到
时对应的竖直高度和目标点的
方向坐标的偏差
，如此顺序搜索，直至某次计算得到的偏差
的值大于上一次的偏差值，结束搜索，得到击中位于
的目标所需的发射角。根据上述思想，得到炮弹对目标实现精确打击时的运动轨迹如
已标注于图中。
我们在计算时采取时间步长
以提高计算精度。可以发现，两种数值方法给出基本一致的结果，如果要打击位于
的目标，应选取发射角
左右。值得说明的是，我们在这里把约化空气阻力系数取为
而不是前面一直采用的
，我们发现在空气阻力过大时，无论以何种发射角出射，炮弹都不能击中目标。这是因为过大的空气阻力过多地耗散了炮弹的能量，只有增大炮弹的发射速度才有可能击中目标。由于我们固定了发射速度为
，那么在这一恒定初始速度下，如果空气阻力增大到一定值时，无论发射角为何，炮弹都将无法命中目标。
结论
我们利用微分方程的数值解法Euler法和Runge-Kutta法研究了加农炮弹出膛后的运动轨迹，考察了忽略空气阻力和考虑空气阻力两种情况下炮弹的运动状态。我们的研究表明，空气阻力的存在会显著影响炮弹运动轨迹的形态，使其水平射程和最高点大幅度降低，空气阻力还会影响炮弹的射程对抛射角的依赖关系，使其达到最大射程时对应的抛射角不同于无阻力的情形。此外，我们还考察了数值解的收敛性问题，研究了数值计算结果与时间步长的关系。我们指出，Euler法和Runge-Kutta法在处理炮弹运动问题时得到的解是收敛的和稳定的，这保证了我们前面作出的一系列结论的可靠性和正确性。最后我们研究了炮弹对目标的打击问题，给出了一个计算例子。我们期望本文得到的研究结果能对炮弹发射相关的实际问题提供一定的帮助。
参考文献
[1] 刘金远, 段萍, 鄂鹏. 计算物理学[M]. 北京: 科学出版社, 2012.
[2] 郝亚非. 抛体运动轨迹的数值分析[J]. 物理通报. 2015(12): 10-12.
注：本文中的全部数值计算Matlab程序代码将于近期上传至GitHub上，可供学习参考。（未经作者允许，不得私自转载本文，一经发现，严厉追究）
扫一扫
4823
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
