电影《隐藏人物》中卫星轨道和再进入计算的现代方法

2016年上映的电影《隐藏人物》获得了良好的评价。它涉及到美国历史上的一个重要时期，触及了许多话题，包括民权和太空竞赛。电影详细介绍了Katherine Johnson 和她的同事（Dorothy Vaughan 和Mary Jackson）在美国宇航局执行水星任务和美国早期探索载人太空飞行期间的不为人知的故事。这部电影重点关注了当时非裔美国妇女在美国宇航局的戏剧性民权斗争，而这些斗争与Johnson和她的同事们的数字计算能力形成了对比。此时，计算机正处于早期阶段，因此，Johnson和她的团队在不使用计算机的情况下完成复杂的导航轨道力学问题的能力，对早期的计算机结果提供了重要的检查。

我将谈谈影片中提到的她的科学工作的两个方面：轨道计算和再进入计算。对于轨道计算，我将首先完全遵循Johnson的做法，然后与利用Wolfram语言提供的一系列工具的更现代、直接的方法进行比较。在电影中提到的使用欧拉方法解决微分方程的地方，我将在火箭的一个重要问题上将这种方法与更现代的方法进行比较：从火箭方程和阻力项（使用直接从Wolfram语言中获得的大气模型数据得出）计算再入轨道。

电影没有过多地关注Johnson和她的团队所处理的问题类型的数学细节，但为了这个博客的目的，我希望至少提供一个在Johnson的时代与现在相比可能使用的方法。

将卫星放置在选定的位置上

Johnson与人合著的最早的论文之一，"为将卫星置于选定的地球位置而确定燃尽时的方位角"，涉及的问题是，在给定的起始位置（如佛罗里达州卡纳维拉尔角）和轨道轨迹的情况下，确保卫星能在指定的轨道数量后被置于特定的地球位置上。约翰逊团队使用的方法是根据其他轨道参数，确定火箭发射的方位角（由航天器在发动机关闭时的速度矢量与一个固定的参考方向，例如北方形成的角度）。这是确保宇航员在返回地球时处于正确位置的一个重要步骤。

常量和初始处理

在论文中，Johnson定义了一些解决手头问题所需的常数和输入参数。需要解释的一个细节是 "燃尽(burnout) "一词，它是指火箭发动机的关闭。燃尽后，轨道参数基本上被 "冻结(frozen)"，航天器只在地球引力下移动（当然是通过牛顿定律确定）。在本节中，我尽可能地遵循论文的单位惯例。

为了方便起见，一些函数被定义为以度数而不是弧度来处理角度。这使得在角度计算中可以顺利地处理时间。

Johnson接着描述了其他几个派生参数，但有趣的是，她有时会采用这些参数的值，而不是使用她的公式所返回的值。她采用的数值往往接近于公式中的数值。为了简单起见，这里使用公式中的数值。

轨道椭圆的半通径：

近地点和燃烧点之间在轨道平面上的角度：

轨道偏心率：

轨道周期：

要描述下一个参数，最简单的是引用原始论文："要求燃尽位置φ1，λ1的卫星在完成n个轨道后经过选定的位置φ2，λ2，相当于要求在第一个轨道期间，卫星经过一个等效位置，其纬度φ2与选定位置的纬度相同，但其经度λ2e从λ2向东移动，移动量足以补偿n个完整轨道期间地球的旋转，即极地时角n ωE T，因此这个等效位置的经度由关系给出"：

角度θ 的近地点的时间：

计算

最终解决方案的一部分是确定中间参数δλ1-2e和θ2e的值。这可以通过几种方式完成。首先，我可以使用ContourPlot通过论文中的方程19和20获得一个图形解：

FindRoot可以用来以数字方式找到解决方案：

当然，Johnson没有机会使用ContourPlot或FindRoot，所以她的论文描述了一种迭代技术。我把论文中描述的技术翻译成了Wolfram语言，并通过她的迭代方法解决了其他一些参数。因为基础计算是针对球形地球的，所以她的方法中包括了对钝化的修正：

对各种迭代的θ2e值作图显示了快速收敛的情况：

我可以在FindRoot命令中转换方法如下（这以完全自洽的方式考虑到了钝化效应，并计算出方程中涉及的所有九个变量的值）：

有趣的是，即使这个更复杂的系统的迭代寻根步骤也能很快收敛：

制图

在确定了轨道参数后，最好能将解决方案可视化。首先，需要从以前的解决方案中提取一些关键参数：

接下来，需要得出卫星的纬度和经度作为方位角的函数：

φs和λs是纬度和经度与θs的关系：

卫星地面轨道可以通过创建一个点的表格来构建：

Johnson的论文提出了一个轨道解决方案的草图，包括显示燃尽、选定和等效位置的标记。在这里很容易重现一个类似的素描图：

作为比较，这是她的原图：

可以用GeoGraphics构建一个更有视觉效果的版本，注意将地心坐标转换为大地坐标：

今天如何计算轨道

今天，几乎我们每个人都可以立即获得比1960年代整个美国国家航空航天局所拥有的计算资源强大得多的计算资源。现在，只需使用一台台式电脑和Wolfram语言，你就可以很容易地找到轨道力学问题的直接数字解决方案，如那些向Katherine Johnson和她的团队提出的问题。虽然与老方法相比，我们的智慧可能不那么费劲，但人们从这些探索中得到的结果却不乏趣味和实用性。

为了用更现代的方法求解方位角ψ，让我们为一个简单的圆形轨道设定参数，在佛罗里达州上空燃烧后开始，假设地球是球形对称的（我不会费力地试图精确匹配Johnson论文中的轨道，我将使用现代SI单位系统重新定义上面的某些数量）。从Johnson使用的相同的低地球轨道高度开始，并使用一点球面三角学，可以直接推导出我们轨道的初始条件：

相关的物理参数可以直接从Wolfram语言中获得：

接下来，我得到一个微分方程，在地球的引力场下，我们的航天器的运动。有几种方法可以模拟地球附近的引力势。假设是一个球形对称的星球，并且自始至终利用直角坐标系，引力势仅限：

或者，你可以使用一个更现实的地球引力模型，其中行星的形状被认为是一个旋转的扁球形椭圆体。来自这样一个椭圆体的势的确切形式（假设椭圆体外壳的质量密度不变），尽管很复杂（包含多个椭圆积分），但可以通过EntityValue获得：

对于一般的同质三轴椭圆体，引力势包含分段函数：

这里，κ是x2/（a2+κ）+y2/（b2+κ）+z2/（c2+κ）=1的最大根。在扁形椭圆体的情况下，前面的公式可以简化为只包含基本函数...

...其中κ=((2 z2 (a2-c2+x2+y2)+(-a2+c2+x2+y2)2+z4)1/2-a2-c2+x2+y2+z2) /2。

在地理和空间科学界广泛使用的一个更简单的形式，我将在这里使用，由所谓的国际重力公式（IGF）给出。IGF考虑到与球面对称势的差异，直到球面谐波的二阶，并给出与前面提到的精确势在数值上没有区别的结果。就四个测量的大地测量参数而言，IGF势可以定义如下：

我可以很容易地通过GeogravityModelData使用甚至更好的引力值。对于起始位置，IGF势力与高阶近似值的偏差只有0.06%：

有了这些引力势的函数形式，寻找轨道路径就相当于取势的梯度来得到引力场矢量，然后应用牛顿第三定律。这样做，我得到了两个重力模型的轨道运动方程：

我现在准备利用NDSolve的力量来计算轨道轨迹。不过，在这之前，将轨道路径显示为三维空间中的曲线会很好。为了给这些曲线提供背景，我将把它们绘制在地球表面的纹理图上，投影到一个球体上。在这里，我构建了所需的图形对象：

虽然在惯性框架中计算的轨道路径形成了一条周期性的封闭曲线，但当你考虑到地球的旋转时，它将导致航天器在随后的每次旋转中经过地球表面的不同点。我可以通过在我从NDSolve获得的解决方案中增加一个额外的旋转项来直观地看到这一效果。出于可视化的目的，我把轨道周期数定为三个（与John Glenn的飞行相似），我构建了以下操纵，以观察轨道路径如何受到方位角发射ψ的影响，与Johnson论文中的研究相似。我将同时绘制一个假设球形地球的路径（白色）和另一个使用IGF的路径（绿色），以了解钝化效应的大小（注意，两个路径的分歧随着每个轨道的增加而增加）：

在本博客所附的笔记本中，你可以看到这个Manipulate的动作，并注意到获得每个新解决方案的速度。你会希望Katherine Johnson和她在美国国家航空航天局的同事们会留下深刻的印象!

现在，改变燃烧时间的角度ψ，就可以直接计算出航天器的位置，比如说，三圈之后：

卫星再进入的建模

电影中还提到了与再进入阶段有关的欧拉方法。正如前几节所做的那样，在找到方位角的初始问题被解决后，是时候回到地球了。火箭被发射，使轨道上的物体减速，当飞船从空间真空过渡到大气环境时，一系列复杂的事件发生了。变化的大气密度、快速减速和摩擦加热都成为必须考虑的重要因素，以使宇航员安全返回地球。高度、速度和作为时间函数的加速度都是需要解决的问题。这组问题可以用欧拉方法来解决，就像Katherine Johnson所做的那样，或者使用Wolfram语言中的微分方程求解功能。

对于简单的微分方程，人们可以通过指定的正交方法得到详细的分步解。对于一个随时间变化的质量m(t)，牛顿著名的F=m a的等价物是所谓的理想火箭方程（在一维）...

...其中m(t)是火箭质量，ve是发动机排气速度，m'p(t)是推进剂质量的时间导数。假设m'p(t)为常数，该方程的结构相对简单，很容易以封闭形式求解：

有了质量的初始和最终条件，我得到了著名的火箭方程（Tsiolkovsky 1903）：

我可以从Wolfram|Alpha获得用具体参数值和例如经典欧拉方法解这个方程的细节。以下是这些细节以及与精确解的详细比较，还有与其他数值积分方法的比较：

按照电影情节，我现在将实现一个再入大气层过程的最小化ODE模型。我首先定义了模仿格伦飞行的参数：

我假设制动过程使用第一阶段发动机1%的推力，并且运行，比如说，持续60秒。运动方程是：

这里，Fgrav是引力，Fexhaust(t)是明确随时间变化的发动机力，Ffriction(x(t),v(t))是摩擦力。后者通过空气密度明确地取决于位置x(t)，通过摩擦规律取决于v(t)。

对于依赖高度的空气密度，我可以方便地使用StandardAtmosphereData函数。由于降落伞在离地面约8.5公里处打开，我还考虑了与高度有关的面积：

这就得到了以下一组需要解决的耦合非线性微分方程。最后的WhenEvent[...]指定在太空舱到达地球表面时结束积分。我使用矢量值的位置和速度变量X和V：

有了这些关于重量、排气和空气摩擦力条款的定义...

......总的力量可以通过以下方式找到：

在这个简单的模型中，我忽略了地球的旋转、太空舱的内在旋转、主动飞行角度的变化、超音速对摩擦力的影响等等。坐标分量中微分方程的明确形式如下。Katherine Johnson解决的方程会和这些很相似：

辅以初始位置和速度，对这个方程组进行数值求解是很直接的。今天，这只是对NDSolve的一个简单调用。我不必担心要使用的方法、步长控制、误差控制等问题，因为Wolfram语言会自动选择保证有意义的结果的值：

下面是高度、速度和加速度与时间的关系图：

绘制成高度而不是时间的函数显示，空气密度的指数增长是造成高减速的原因。这不是由于降落伞的原因，因为降落伞发生在相对较低的高度。减速的峰值发生在一个非常高的高度，因为太空舱从真空到大气环境的速度非常快：

而这里是太空舱在再入过程中的垂直和切向速度图：

现在我用固定步数的欧拉方法重复数值解：

从质量上看，该解决方案看起来与之前的解决方案相同：

对于所使用的时间积分的步长，累积误差在百分之几的数量级上。更小的步长会减少误差（见之前的Wolfram|Alpha输出）：

请注意，由欧拉方法预测的着陆时间与之前的时间仅有0.11%的偏差。(作为比较，如果我用两种现代方法来解这个方程，比如 "BDF "与 "Adams"，误差会小几个数量级)。

现在，重返大气层的过程产生了大量的热量。这就是需要热防护罩的地方。在哪个高度产生的单位面积q的热量最大？在没有详细推导的情况下，我可以从纯粹的尺寸理由来猜测：

还有很多有趣的东西可以计算出来（Hicks 2009），但就像电影必须把所有东西都塞进两小时七分钟一样，现在我的博客也将因为时间的关系而结束。Wolfram 语言可以让您畅游无阻！

更多实验、建模

数学应用案例，软件技巧

考研、竞赛真题、资讯

考研竞赛数学(ID: xwmath)