2.2 运动方程

本节将描述用于确定再入飞行器轨迹的运动方程。首先，在计算中使用的参考系及各参考系之间的变换矩阵将在2.2.1节给出。其次，轨迹确定中需要的受力将会在2.2.2节给出。最后，一个关于再入问题的在球坐标系中的运动合成方程将在2.2.3节给出。

2.2.1 坐标系

为了描述运动方程，有必要定义多个参考系，因为外力是在不同参考系下描述的。除了参考系的定义之外，需要知道这些参考系之间的转换，以便在同一个参考系下表述所有相关分量。在之后的叙述中，所有参考系都遵循右手系原则。

①惯性平面中心坐标系I (图2.2)。

本书使用的惯性参考系原点为地球质心，z轴指向旋转方向 (北方为正), x轴指向参考方向，y轴与xz轴构成右手系。对于地球而言，x轴的参考方向通常是春分点的方向 (白羊座的第一点，Υ)。通常以J2000参考系为标准参考系。

②旋转平面中心坐标系R (图2.2)。

旋转系的原点与惯性系相同，但固定在地球上。z轴是旋转轴，x轴通常与格林尼治子午线相交 ( τ = 0°, τ是经度)。I系和R系的x轴和y轴之间的夹角称为格林威治平恒星时 (GMST)，并表示为θGMST。根据它在t=0处的值和角速度ωcb可以确定惯性和行星固定系之间的转换。然而，由于这个角度对于所有的模拟来说都是恒定的，所以它也可以被置零，并且对结果不会有影响。因此，为了简单起见，我们假设I和R系之间在t=0处共线。此外，地球的角速度及旋转轴的防线刚被认为是恒定的，因此任何章动或进动在此被忽略。

③垂直坐标系V (图2.2和图2.3)。

垂直参考系以机身为中心，原点位于机身质心。各轴的方向基于相对中心天体的位置确定。z轴指向中心天体的质心，xy平面是局部水平的平面 (对于球形中心天体而言)，其中x轴指向北，y轴与x轴和z轴构成右手系。

④轨迹参考系T (图2.3和图2.4)。

轨迹参考系以机身为中心，因此原点也在飞行器质心。轨迹参考系的x轴指向飞行方向，z轴位于垂直平面内，指向下方，y轴与x轴和z轴构成右手系。应当注意的是，由于z轴必须垂直于x轴并且位于垂直平面中，所以其方向已被定义，除了符号 (定义为正向下) 外。

⑤气动坐标系A (图2.4)。

气动参考系的原点在飞行器质心。x轴是速度矢量的方向，z轴与升力共线(图2.5)，但方向相反，y轴与x轴和z轴构成右手系。可以定义两种气动参考系，一个参考空速，一个参考地面速度。然而，由于没有考虑风速，这两种参考系重合并且两者之间没有区别。值得注意的是，当倾侧角σ为零时，A系和T系重合。

⑥机体坐标系B (见图2.4)。

机体坐标系的原点位于机体质心，但是不同于气动参考系，该原点固定在机体质心上。x轴通常指向机体前方，z轴指向下，y轴与x轴和z轴构成右手系。

使用欧拉角来处理各坐标系之间的转换。利用欧拉角，坐标系之间的旋转可以通过以此围绕对应轴线的旋转来表示。应该强调的是，旋转的顺序并不是可以忽略的，因为改变它们会导致不同的坐标转换。关于第i个轴的正定义角度θi的旋转方向余弦矩阵为

这些变换被称为单位轴变换，因为每个变换都会给出x轴、y轴或z轴的旋转。

从系统A到系统C的总变换矩阵是由这些单独矩阵构成的，如关于旋转顺序j, k, l (首先围绕轴线j旋转，然后围绕轴线k，最后围绕轴线l) 的矩阵为

这些矩阵的顺序来自一个事实，即应该连续地将矩阵与其相关的矩阵左乘。类似地，当从系统A变化到系统C时，可以使用中间系统B (如CB, A和CC, B可以直接使用)，即

尽管存在很多使系统A变换到系统C的欧拉角集合，但是这些角可能不存在很清晰的物理解释。这使得选择中间系统B成为首选方法。应该理解的是，在等式(2.16) 中的一系列三次旋转中，围绕同一轴线的非连续两次旋转可能不能简单地作为绕该轴线旋转两次角度的总和。这是由于在中间旋转过程中轴线的方向发生了改变。

Mooij (1998) 提出了一系列相关变换的轮换顺序：

式中，γ和χ为轨迹角和航向角；α和β为迎角和侧滑角。每个角度的定义已经在本节给出 (图2.2~图2.4)。应该注意的是，逆变换矩阵是很容易得到的，因为变换矩阵都是正交的，所以这个矩阵的倒数就很容易求得。

2.2.2 作用力

本节将描述飞行器平移运动的力。我们将陈述用来描述飞行器平移运动的假设，讨论所得到的方程的物理意义。

这里假定空气动力和引力是作用于机体的唯一外力，所以有

式中，下标I表示在惯性坐标系下描述这些外力。由于只考虑无动力再入，所以在外力中不包括推力。

在前文中已经定义了引力 (见2.1.2节)，即

式中gV= (gn0 gd)T，这些引力分量在V系下给出，所以为了获得I系下的各分量我们需要使用静态坐标变换：

式中，CI, V=CI, RCR, V；式 (2.18) ~式 (2.22) 中给出了该式矩阵的转置矩阵CR, I和CV, R。

空气动力参数定义在空气动力坐标系A中，并由以下公式给出：

式中，CD、CS和CL分别代表阻力、侧向力和升力系数；V代表飞行器相对于大气的速度。对于某些空气动力数据库，侧向力在正yA方向上定义为正值。然而，在这里可以观察到，如式 (2.26) 所示，以右手方式定义的空气动力坐标系中的三个力系数的公式，图2.5为示意图。或者，空气动力也可以通过以下方式在机体坐标系B中表示：

其中两者的关系来自2.2.1节中描述的转换。空气动力可以通过以下方式在坐标系I中表示：

例如，C RCR, VCV, TCT, A为一个前文定义的 (逆) 矩阵组合。

2.2.3 再入方程

平移运动方程直接遵循牛顿运动定律，这些物理定律仅仅在惯性坐标系下成立。然而，通过引入在目标坐标系和惯性坐标系之间的相对旋转并且包含由此产生的表观力 (如科里奥利，向心)，这些运动方程也可以在旋转系下表示。方程的结果公式由状态变量的选择决定。关于通过三维大气层的平移运动，需要6个状态变量来完全确定，特别是3个位置坐标和3个速度分量。在我们的数值积分中，状态由笛卡儿位置和I系下的速度给出。

在最简单的形式中，描述惯性运动的状态变量是位置参数rI=(xIyIzI)和速度参数VI=(·x I· yI· zI)。将空气动力和引力转换到I坐标系下，运动的一般方程为

为了从这组6个一阶、非线性、耦合的常微分方程组中确定飞行器的运动轨迹，由于不存在一般解析解，所以必须采用数值积分方法。数值积分方法有多种，都具有不同的精度、计算复杂性和数值稳定性。由于本书研究的概念性质，再入轨迹所需的准确性有限。这是由建模误差而产生的，如空气动力系数和飞行器质量也会导致轨迹上的错误。因为必须传播大量的轨迹，所以数值计算的时间应该保持最小。可以使用变步长积分器，在结果变化较小的区域中使用更大的积分步长。然而，使用具有一定更新频率的制导系统会使这些较大的时间步长不太适用于飞行器行为的比较分析。

在此，我们使用一个比较常见的数值积分方法，即四阶龙格-库塔法，也可以使用更加精确的变步长方法，比较典型的例子是由Dormand-Prince 4 (5) 给出的一种结合五阶方法的变步长四阶方法，在给定预定相对或/和绝对容限的情况下确定步长，或者不同精度组合的龙格-库塔-费尔贝格法。有关这些方法和其他方法的详细信息，包括方法稳定性、一致性和收敛性等有关问题，可以在各种文献中找到，如Shampine (1994)、Lambert (1991) 和Press等 (2007)。虽然式 (2.29) 可能被用来有效地模拟飞行器轨迹，但是不同的表征方式可能会更加有用，无论是从物理观点还是计算观点出发。一个很好的例子是球面位置和速度的使用，这产生了一个非常直观的公式和运动方程的表达。这些方法通常用于制导系统的开发，即我们在本节中所采用的近似值 (见2.3.2节)。

球坐标位置和速度分量如图2.6所示，这6个状态分量定义如下：

①径向位置r，地心到飞行器质心的标量距离；

②经度τ，从格林尼治线开始测量，向东为正 ( -180°≤τ≤180°)；

③纬度δ，从赤道开始沿着当地子午线测量，向北为正 ( -90°≤δ≤90°)；

④ (相对) 速度V，标量速度被测量为地面速度，因为无风在本书中等同于空气速度。注意到由于地球旋转，该速度与惯性坐标的速度不同，所以对于加速级的再入来说，V<VI。

⑤飞行路径角度γ，即地面速度矢量与局部水平平面之间的夹角 ( -90°≤γ≤90°)；

⑥航向角χ，也就是在平面上的北方方向和投影到这个平面上的地面速度矢量 ( -180°≤χ≤180°)。按照顺时针方向测量，0°航向角向北，90°向东。

这里所说的速度矢量是相对于旋转的中心体，所以方程 (2.29) 需要被调整以适合中心体的旋转。在这个过程中，速度被考虑进去。飞行器的位置矢量与地球局部曲率的切线的关系定义了垂直平面，垂直于这个平面的就叫局部水平平面，包含了飞行器的重心和水平速度分量。同样，这个平面是由速度坐标系V的x轴和y轴定义的。当飞行器移动时，该局部水平面也跟随飞行器移动，并且这个平面的旋转必须应用到方程 (2.29) 的二次修正。

当使用r, τ, δ, V, γ和χ对运动方程求偏导时，推导过程冗长而单调，我们只要陈述最终的结果就足够了。注意，我们假设S=0，分析得出β =0。欲了解更多的信息，可参考Vinh等 (1980)、Regan和Anandakrishnan (1993) 和Mooj (1998) 等的相关书籍。

应该注意的是，上面的这些方程有两个奇异点：在方程 (2.31) 中这个奇异点是在γ =± 90°时遇到，这意味着一个纯粹的垂直飞行 (向上或向下)；第二个奇异点是在δ =± 90°时遇到，它在北极或者南极。在本书的研究过程中，这两种情况都不会发生，所以这些方程将作为理论指导的基础。

通常，这些方程被改写为L/D =(CL/CD)，以及一个轨迹参数或系数，定义为

其中，质量有时会被重量 (海平面) 代替，轨迹参数表示引力和空气阻力的相对影响。它可以用来指示出现最大加热、减速等的高度。但是，一般来说，由于飞行角度、速度和高度的变化，阻力系数会发生变化，因此一般情况下轨迹参数不会是常数。