机动

下面的教程详细介绍了机动的一些基本用法。包括简单的瞬时机动和更复杂的连续推力机动。

脉冲机动

脉冲机动是对航天器速度进行离散变化的简单操作。它们主要用于两种情况:

  • 用于保持轨道稳定,由于所需机动的规模较小,因此它们仍然是现实的模型
  • 用于初始化优化或搜索算法,该算法将在内部使用更现实的模型进行大规模机动

本教程展示了如何实现一系列机动,逐渐改变轨道的倾角。

让我们设置一个初始状态:

final Frame eme2000 = FramesFactory.getEME2000();
final Orbit initialOrbit =
                new KeplerianOrbit(8000000.0, 0.01,
                                   FastMath.toRadians(50.0), // ← 这是初始倾角
                                   FastMath.toRadians(140.0),
                                   FastMath.toRadians(12.0),
                                   FastMath.toRadians(-60.0), PositionAngleType.MEAN,
                                   eme2000,
                                   new AbsoluteDate(new DateComponents(2008, 6, 23),
                                                    new TimeComponents(14, 0, 0),
                                                    TimeScalesFactory.getUTC()),
                                   Constants.EIGEN5C_EARTH_MU);

机动将在航天器坐标系中定义。我们需要确保Z轴与轨道动量对齐,因此我们选择与LVLH局部轨道坐标系对齐的姿态

final AttitudeProvider attitudeProvider = new LofOffset(eme2000, LOFType.LVLH);

我们想要执行一系列3次倾角减小的机动。由于它们只修改倾角,所以它们必须发生在升交点,但并非所有升交点都适合,我们希望选择沿着-Z方向的升交点。当发生Action.STOP事件时触发机动(并将其过滤掉)

final NodeDetector ascendingNodeStopper =
                new NodeDetector(FramesFactory.getEME2000()).
                withMaxCheck(300.0).
                withThreshold(1.0e-6).
                withHandler(new StopOnIncreasing());

我们只允许在前3个轨道上进行机动

final AbsoluteDate lastAllowedDate =
                initialOrbit.getDate().shiftedBy(3 * initialOrbit.getKeplerianPeriod());
final EnablingPredicate<EventDetector> predicate =
                (state, detector, g) -> state.getDate().isBefore(lastAllowedDate);
final EventDetector trigger =
                new EventEnablingPredicateFilter<>(ascendingNodeStopper, predicate);

使用升交点探测器作为触发器创建机动

final ImpulseManeuver<EventDetector> maneuver =
                new ImpulseManeuver<>(trigger,
                                      new Vector3D(0.0, 0.0, -122.25), // ← 沿着-Z方向的122.25 m/s
                                      350.0);

将所有内容包装在一个传播器中。请注意,ImpulseManeuver是一个事件检测器,而不是一个力模型!这使得它可以用于所有的传播器,包括在这里使用的开普勒传播器。

final KeplerianPropagator propagator = new KeplerianPropagator(initialOrbit, attitudeProvider);
propagator.addEventDetector(maneuver);

进度监控

propagator.getMultiplexer().add(900.0, (state, isLast) -> {
    final Vector3D pos = state.getPVCoordinates(eme2000).getPosition();
    System.out.format(Locale.US, "%s %s hemisphere inclination = %5.3f%n",
                      state.getDate(),
                      pos.getZ() > 0 ? "Northern" : "Southern",
                      FastMath.toDegrees(state.getOrbit().getI()));
});

运行模拟

propagator.propagate(initialOrbit.getDate().shiftedBy(5 * initialOrbit.getKeplerianPeriod()));

打印结果如下所示。我们可以看到,当我们穿越降交点(即从北半球切换到南半球)时,倾角保持不变,并且在穿越前三个升交点时发生变化。

2008-06-23T14:00:00.000Z 北半球倾角 = 50.000
2008-06-23T14:15:00.000Z 北半球倾角 = 50.000
2008-06-23T14:30:00.000Z 北半球倾角 = 50.000
2008-06-23T14:45:00.000Z 南半球倾角 = 50.000
2008-06-23T15:00:00.000Z 南半球倾角 = 50.000
2008-06-23T15:15:00.000Z 南半球倾角 = 50.000
2008-06-23T15:30:00.000Z 南半球倾角 = 50.000
2008-06-23T15:45:00.000Z 北半球倾角 = 49.000
2008-06-23T16:00:00.000Z 北半球倾角 = 49.000
2008-06-23T16:15:00.000Z 北半球倾角 = 49.000
2008-06-23T16:30:00.000Z 北半球倾角 = 49.000
2008-06-23T16:45:00.000Z 南半球倾角 = 49.000
2008-06-23T17:00:00.000Z 南半球倾角 = 49.000
2008-06-23T17:15:00.000Z 南半球倾角 = 49.000
2008-06-23T17:30:00.000Z 南半球倾角 = 49.000
2008-06-23T17:45:00.000Z 北半球倾角 = 48.001
2008-06-23T18:00:00.000Z 北半球倾角 = 48.001
2008-06-23T18:15:00.000Z 北半球倾角 = 48.001
2008-06-23T18:30:00.000Z 北半球倾角 = 48.001
2008-06-23T18:45:00.000Z 南半球倾角 = 48.001
2008-06-23T19:00:00.000Z 南半球倾角 = 48.001
2008-06-23T19:15:00.000Z 南半球倾角 = 48.001
2008-06-23T19:30:00.000Z 南半球倾角 = 48.001
2008-06-23T19:45:00.000Z 北半球倾角 = 47.001
2008-06-23T20:00:00.000Z 北半球倾角 = 47.001
2008-06-23T20:15:00.000Z 北半球倾角 = 47.001
2008-06-23T20:30:00.000Z 南半球倾角 = 47.001
2008-06-23T20:45:00.000Z 南半球倾角 = 47.001
2008-06-23T21:00:00.000Z 南半球倾角 = 47.001
2008-06-23T21:15:00.000Z 南半球倾角 = 47.001
2008-06-23T21:30:00.000Z 北半球倾角 = 47.001
2008-06-23T21:45:00.000Z 北半球倾角 = 47.001
2008-06-23T22:00:00.000Z 北半球倾角 = 47.001
2008-06-23T22:15:00.000Z 北半球倾角 = 47.001
2008-06-23T22:30:00.000Z 南半球倾角 = 47.001
2008-06-23T22:45:00.000Z 南半球倾角 = 47.001
2008-06-23T23:00:00.000Z 南半球倾角 = 47.001
2008-06-23T23:15:00.000Z 南半球倾角 = 47.001
2008-06-23T23:30:00.000Z 北半球倾角 = 47.001
2008-06-23T23:45:00.000Z 北半球倾角 = 47.001

这个示例的完整代码可以在教程的源代码树中找到,文件名为src/main/java/org/orekit/tutorials/maneuvers/ImpulseAtNode.java

连续机动

连续机动是一种真实的模型,考虑了机动持续时间、机动期间的姿态变化和质量消耗。它们只能与基于积分的传播器一起使用。

本教程展示了如何实现一个远地点机动,可以使用传播器级别设置的姿态,也可以仅覆盖机动加速度方向计算。我们只使用基于日期的触发器和恒定推力推进系统,但也可以使用其他不同的推进系统。例如,我们可以将BasicConstantThrustPropulsionModel替换为ScaledConstantThrustPropulsionModel,并考虑一些校准因子(或者如果我们将此模型用于轨道确定而不是模拟中,则将这些缩放因子配置为估计值)。我们还可以将基于日期的触发器替换为基于事件的触发器,可以模拟多次燃烧机动。

让我们设置一个具有GTO轨道和2500kg的航天器的初始状态:

final Frame eme2000 = FramesFactory.getEME2000();
final AbsoluteDate date = new AbsoluteDate(new DateComponents(2004, 01, 01),
                                           new TimeComponents(23, 30, 00.000),
                                           TimeScalesFactory.getUTC());
final Orbit orbit =
                new KeplerianOrbit(24396159, 0.72831215, FastMath.toRadians(7),
                                   FastMath.toRadians(180), FastMath.toRadians(261),
                                   FastMath.toRadians(0), PositionAngleType.TRUE,
                                   eme2000, date,
                                   Constants.EIGEN5C_EARTH_MU);
final SpacecraftState initialState = new SpacecraftState(orbit, 2500.0);

准备传播器,使用自适应步长积分器。传播器将使用与VNC对齐的姿态模式,即其X轴始终沿轨道速度方向

final OrbitType orbitType = OrbitType.EQUINOCTIAL;
final double[][] tol = NumericalPropagator.tolerances(1.0, orbit, orbitType);
final AdaptiveStepsizeIntegrator integrator =
                new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]);
integrator.setInitialStepSize(60);
final NumericalPropagator propagator = new NumericalPropagator(integrator);
propagator.setOrbitType(orbitType);
propagator.setInitialState(initialState);
propagator.setAttitudeProvider(new LofOffset(eme2000, LOFType.VNC));

首先,我们希望将机动计算为惯性机动,因此不能依赖于上面配置的姿态模式,我们需要一个覆盖姿态的定律,使X轴指向特定方向

final Vector3D direction = new Vector3D(FastMath.toRadians(-7.4978),
                                        FastMath.toRadians(351));
final AttitudeProvider attitudeOverride =
                new InertialProvider(new Rotation(direction, Vector3D.PLUS_I), eme2000);

Maneuver(机动)将在已知日期开始,并在已知持续时间后停止

final AbsoluteDate firingDate = new AbsoluteDate(new DateComponents(2004, 1, 2),
                                                 new TimeComponents(4, 15, 34.080),
                                                 TimeScalesFactory.getUTC());
final double duration = 3653.99;
final ManeuverTriggers triggers = new DateBasedManeuverTriggers(firingDate, duration);

Maneuver(机动)具有恒定的推力

final double thrust = 420;
final double isp    = 318;
final PropulsionModel propulsionModel =
               new BasicConstantThrustPropulsionModel(thrust, isp,
                                                      Vector3D.PLUS_I,
                                                      "apogee-engine");

构建机动并将其作为新的力模型添加到传播器中

propagator.addForceModel(new Maneuver(attitudeOverride, triggers, propulsionModel));

进度监控

propagator.getMultiplexer().add(120.0, (state, isLast) ->
    System.out.format(Locale.US, "%s a = %12.3f m, e = %11.9f, m = %8.3f kg%n",
                      state.getDate(), state.getA(), state.getE(), state.getMass())
);

传播轨道,包括机动

propagator.propagate(fireDate.shiftedBy(-900), fireDate.shiftedBy(duration + 900));

下面显示了打印的结果。我们可以看到,在机动之前,半长轴、离心率和倾角保持不变,在机动期间它们持续变化,并在机动之后再次保持不变

2004-01-02T04:00:34.080 a = 24396159.000 m, e = 0.728312150, m = 2500.000 kg
2004-01-02T04:02:34.080 a = 24396159.000 m, e = 0.728312150, m = 2500.000 kg
2004-01-02T04:04:34.080 a = 24396159.000 m, e = 0.728312150, m = 2500.000 kg
2004-01-02T04:06:34.080 a = 24396159.000 m, e = 0.728312150, m = 2500.000 kg
2004-01-02T04:08:34.080 a = 24396159.000 m, e = 0.728312150, m = 2500.000 kg
2004-01-02T04:10:34.080 a = 24396159.000 m, e = 0.728312150, m = 2500.000 kg
2004-01-02T04:12:34.080 a = 24396159.000 m, e = 0.728312150, m = 2500.000 kg
2004-01-02T04:14:34.080 a = 24396159.000 m, e = 0.728312150, m = 2500.000 kg
2004-01-02T04:16:34.080 a = 24442112.400 m, e = 0.725043403, m = 2491.919 kg
2004-01-02T04:18:34.080 a = 24536037.252 m, e = 0.718404119, m = 2475.758 kg
2004-01-02T04:20:34.080 a = 24632720.409 m, e = 0.711627339, m = 2459.596 kg
2004-01-02T04:22:34.080 a = 24732246.157 m, e = 0.704710905, m = 2443.435 kg
2004-01-02T04:24:34.080 a = 24834702.625 m, e = 0.697652633, m = 2427.273 kg
2004-01-02T04:26:34.080 a = 24940182.000 m, e = 0.690450311, m = 2411.112 kg
2004-01-02T04:28:34.080 a = 25048780.751 m, e = 0.683101699, m = 2394.950 kg
2004-01-02T04:30:34.080 a = 25160599.875 m, e = 0.675604529, m = 2378.788 kg
2004-01-02T04:32:34.080 a = 25275745.160 m, e = 0.667956507, m = 2362.627 kg
2004-01-02T04:34:34.080 a = 25394327.460 m, e = 0.660155308, m = 2346.465 kg
2004-01-02T04:36:34.080 a = 25516463.000 m, e = 0.652198582, m = 2330.304 kg
2004-01-02T04:38:34.080 a = 25642273.694 m, e = 0.644083951, m = 2314.142 kg
2004-01-02T04:40:34.080 a = 25771887.491 m, e = 0.635809009, m = 2297.981 kg
2004-01-02T04:42:34.080 a = 25905438.748 m, e = 0.627371325, m = 2281.819 kg
2004-01-02T04:44:34.080 a = 26043068.626 m, e = 0.618768440, m = 2265.658 kg
2004-01-02T04:46:34.080 a = 26184925.522 m, e = 0.609997872, m = 2249.496 kg
2004-01-02T04:48:34.080 a = 26331165.531 m, e = 0.601057114, m = 2233.335 kg
2004-01-02T04:50:34.080 a = 26481952.946 m, e = 0.591943638, m = 2217.173 kg
2004-01-02T04:52:34.080 a = 26637460.795 m, e = 0.582654892, m = 2201.012 kg
2004-01-02T04:54:34.080 a = 26797871.426 m, e = 0.573188310, m = 2184.850 kg
2004-01-02T04:56:34.080 a = 26963377.135 m, e = 0.563541304, m = 2168.688 kg
2004-01-02T04:58:34.080 a = 27134180.848 m, e = 0.553711279, m = 2152.527 kg
2004-01-02T05:00:34.080 a = 27310496.862 m, e = 0.543695627, m = 2136.365 kg
2004-01-02T05:02:34.080 a = 27492551.643 m, e = 0.533491736, m = 2120.204 kg
2004-01-02T05:04:34.080 a = 27680584.703 m, e = 0.523096997, m = 2104.042 kg
2004-01-02T05:06:34.080 a = 27874849.544 m, e = 0.512508806, m = 2087.881 kg
2004-01-02T05:08:34.080 a = 28075614.691 m, e = 0.501724576, m = 2071.719 kg
2004-01-02T05:10:34.080 a = 28283164.817 m, e = 0.490741747, m = 2055.558 kg
2004-01-02T05:12:34.080 a = 28497801.970 m, e = 0.479557796, m = 2039.396 kg
2004-01-02T05:14:34.080 a = 28719846.917 m, e = 0.468170253, m = 2023.235 kg
2004-01-02T05:16:34.080 a = 28937941.941 m, e = 0.457162297, m = 2007.882 kg
2004-01-02T05:18:34.080 a = 28937941.941 m, e = 0.457162297, m = 2007.882 kg
2004-01-02T05:20:34.080 a = 28937941.941 m, e = 0.457162297, m = 2007.882 kg
2004-01-02T05:22:34.080 a = 28937941.941 m, e = 0.457162297, m = 2007.882 kg
2004-01-02T05:24:34.080 a = 28937941.941 m, e = 0.457162297, m = 2007.882 kg
2004-01-02T05:26:34.080 a = 28937941.941 m, e = 0.457162297, m = 2007.882 kg
2004-01-02T05:28:34.080 a = 28937941.941 m, e = 0.457162297, m = 2007.882 kg
2004-01-02T05:30:34.080 a = 28937941.941 m, e = 0.457162297, m = 2007.882 kg

如果我们不想覆盖姿态,而是想使用推进器中配置的姿态(这里是VNC对齐的姿态),那么在构建机动时,我们只需将姿态覆盖参数设置为null:

propagator.addForceModel(new Maneuver(null, triggers, propulsionModel));

使用这种配置的结果将变为:

2004-01-02T04:00:34.080 a = 24396159.000 m, e = 0.728312150, m = 2500.000 kg
2004-01-02T04:02:34.080 a = 24396159.000 m, e = 0.728312150, m = 2500.000 kg
2004-01-02T04:04:34.080 a = 24396159.000 m, e = 0.728312150, m = 2500.000 kg
2004-01-02T04:06:34.080 a = 24396159.000 m, e = 0.728312150, m = 2500.000 kg
2004-01-02T04:08:34.080 a = 24396159.000 m, e = 0.728312150, m = 2500.000 kg
2004-01-02T04:10:34.080 a = 24396159.000 m, e = 0.728312150, m = 2500.000 kg
2004-01-02T04:12:34.080 a = 24396159.000 m, e = 0.728312150, m = 2500.000 kg
2004-01-02T04:14:34.080 a = 24396159.000 m, e = 0.728312150, m = 2500.000 kg
2004-01-02T04:16:34.080 a = 24445841.001 m, e = 0.724984968, m = 2491.919 kg
2004-01-02T04:18:34.080 a = 24547017.613 m, e = 0.718218653, m = 2475.758 kg
2004-01-02T04:20:34.080 a = 24650693.763 m, e = 0.711301532, m = 2459.596 kg
2004-01-02T04:22:34.080 a = 24756972.407 m, e = 0.704231757, m = 2443.435 kg
2004-01-02T04:24:34.080 a = 24865960.944 m, e = 0.697007471, m = 2427.273 kg
2004-01-02T04:26:34.080 a = 24977771.472 m, e = 0.689626812, m = 2411.112 kg
2004-01-02T04:28:34.080 a = 25092521.078 m, e = 0.682087911, m = 2394.950 kg
2004-01-02T04:30:34.080 a = 25210332.139 m, e = 0.674388891, m = 2378.788 kg
2004-01-02T04:32:34.080 a = 25331332.649 m, e = 0.666527867, m = 2362.627 kg
2004-01-02T04:34:34.080 a = 25455656.575 m, e = 0.658502951, m = 2346.465 kg
2004-01-02T04:36:34.080 a = 25583444.233 m, e = 0.650312244, m = 2330.304 kg
2004-01-02T04:38:34.080 a = 25714842.703 m, e = 0.641953845, m = 2314.142 kg
2004-01-02T04:40:34.080 a = 25850006.270 m, e = 0.633425850, m = 2297.981 kg
2004-01-02T04:42:34.080 a = 25989096.901 m, e = 0.624726354, m = 2281.819 kg
2004-01-02T04:44:34.080 a = 26132284.765 m, e = 0.615853455, m = 2265.658 kg
2004-01-02T04:46:34.080 a = 26279748.789 m, e = 0.606805258, m = 2249.496 kg
2004-01-02T04:48:34.080 a = 26431677.269 m, e = 0.597579879, m = 2233.335 kg
2004-01-02T04:50:34.080 a = 26588268.522 m, e = 0.588175453, m = 2217.173 kg
2004-01-02T04:52:34.080 a = 26749731.602 m, e = 0.578590139, m = 2201.012 kg
2004-01-02T04:54:34.080 a = 26916287.072 m, e = 0.568822132, m = 2184.850 kg
2004-01-02T04:56:34.080 a = 27088167.851 m, e = 0.558869672, m = 2168.688 kg
2004-01-02T04:58:34.080 a = 27265620.125 m, e = 0.548731056, m = 2152.527 kg
2004-01-02T05:00:34.080 a = 27448904.347 m, e = 0.538404659, m = 2136.365 kg
2004-01-02T05:02:34.080 a = 27638296.331 m, e = 0.527888947, m = 2120.204 kg
2004-01-02T05:04:34.080 a = 27834088.441 m, e = 0.517182503, m = 2104.042 kg
2004-01-02T05:06:34.080 a = 28036590.891 m, e = 0.506284055, m = 2087.881 kg
2004-01-02T05:08:34.080 a = 28246133.176 m, e = 0.495192508, m = 2071.719 kg
2004-01-02T05:10:34.080 a = 28463065.636 m, e = 0.483906983, m = 2055.558 kg
2004-01-02T05:12:34.080 a = 28687761.178 m, e = 0.472426869, m = 2039.396 kg
2004-01-02T05:14:34.080 a = 28920617.165 m, e = 0.460751878, m = 2023.235 kg
2004-01-02T05:16:34.080 a = 29149754.286 m, e = 0.449481220, m = 2007.882 kg
2004-01-02T05:18:34.080 a = 29149754.286 m, e = 0.449481220, m = 2007.882 kg
2004-01-02T05:20:34.080 a = 29149754.286 m, e = 0.449481220, m = 2007.882 kg
2004-01-02T05:22:34.080 a = 29149754.286 m, e = 0.449481220, m = 2007.882 kg
2004-01-02T05:24:34.080 a = 29149754.286 m, e = 0.449

正如预期的那样,我们可以看到当姿态始终与速度保持一致时,对于相同的燃料消耗,半长轴的增加更多,离心率的减少也更多,而当姿态是惯性时则不然。

此示例的完整代码可以在教程的源代码树中找到,文件名为src/main/java/org/orekit/tutorials/maneuvers/ApogeeManeuver.java