圆形轨道中重力的矢量方向

Question

我目前正在使用 C# 进行一个项目，我在其中玩弄行星引力，我知道这是一个要充分利用它的核心主题，但我喜欢挑战。我一直在阅读牛顿定律和开普勒定律，但我无法弄清楚的一件事是如何获得正确的引力方向。

在我的示例中，我只有 2 个身体。一颗卫星和一颗行星。这是为了简化它，所以我可以掌握它 - 但我的计划是让多个对象动态地相互影响，并希望最终得到一个有点现实的多体系统。

当你有一个轨道时，卫星就有一个引力，这当然是在行星的方向上，但那个方向不是一个常数。为了更好地解释我的问题，我将尝试使用一个示例：

假设我们有一颗卫星以 50 m/s 的速度移动，并以 10 m/s/s 的速度向地球加速，半径为 100 m。 （所有理论数字）如果我们说帧速率为 1，那么一秒钟后对象将向前 50 个单位，向下 10 个单位。

当卫星在一帧内移动多个单元并移动大约 50% 的半径时，在这一帧期间，重力方向发生了很大变化，但施加的力只是“向下”。这会产生很大的误差，尤其是当对象移动半径的很大百分比时。

在我们的示例中，我们可能需要我们的重力方向基于我们当前位置和这一帧结束时的位置之间的平均值。

如何计算这个？

我对三角学有基本的了解，但主要关注三角形。假设我很愚蠢，因为与你们中的任何人相比，我可能是。

（我提出了一个先前的问题，但最终删除了它，因为它产生了一些敌意，而且措辞基本上没有那么好，而且都是笼统的 - 这不是一个真正的具体问题。我希望这会更好。如果不是，那么请通知我，我是来学习的：））

仅供参考，这是我现在的运动功能：

foreach (ExtTerBody OtherObject in UniverseController.CurrentUniverse.ExterTerBodies.Where(x => x != this))
{
    double massOther = OtherObject.Mass;

    double R = Vector2Math.Distance(Position, OtherObject.Position);

    double V = (massOther) / Math.Pow(R,2) * UniverseController.DeltaTime;

    Vector2 NonNormTwo = (OtherObject.Position - Position).Normalized() * V;

    Vector2 NonNormDir = Velocity + NonNormTwo;
    Velocity = NonNormDir;

    Position += Velocity * Time.DeltaTime;
}

如果我的措辞很糟糕，请让我重新措辞 - 英语不是我的母语，如果您不知道正确的技术术语，特定主题可能很难措辞。 :)

我有一种预感，这在开普勒第二定律中有所涵盖，但如果是这样，那么我不确定如何使用它，因为我没有完全理解他的定律。

感谢您的宝贵时间 - 这意味着很多！

（如果有人在我的函数中发现多个错误，请指出！）

Answer 1

我目前正在使用 C# 进行一个项目，我在其中玩弄行星引力

这是一种同时学习模拟技术、编程和物理的有趣方式。

我无法弄清楚的一件事是如何获得正确的重力方向。

我假设您不是在尝试模拟相对论引力。地球不在围绕太阳的轨道上，地球在围绕八分钟前太阳所在的轨道。纠正引力不是瞬时的这一事实可能很困难。 （更新：根据评论，这是不正确的。我知道什么；第二年牛顿动力学后我停止学习物理，对张量微积分只有最模糊的了解。）

在这个早期阶段，您最好假设引力是瞬时的，并且行星是所有质量都在中心的点。引力矢量是从一点到另一点的直线。

假设我们有一颗以 50 m/s 的速度移动的卫星......如果我们说帧速率是每秒一帧，那么一秒后对象将向右移动 50 个单位，向下移动 10 个单位。

让我们更清楚地说明这一点。力等于质量乘以加速度。你计算出身体之间的力。你知道它们的质量，所以你现在知道每个物体的加速度。每个物体都有一个位置和一个速度。加速度改变速度。速度改变位置。因此，如果粒子开始时向左的速度为 50 m/s，向下的速度为 0 m/s，然后施加一个力使其向下加速 10 m/s/s，那么我们可以计算出变化为速度，然后改变位置。正如您所注意到的，在那一秒结束时，位置和速度与它们现有的大小相比都会发生巨大的变化。

当卫星在一帧内移动多个单元并移动大约 50% 的半径时，在这一帧中，重力方向发生了很大变化，但施加的力只是“向下”。这会产生很大的误差，尤其是当对象移动半径的很大百分比时。

正确。问题在于帧速率非常太低，无法正确模拟您所描述的交互。您需要运行模拟，以便在对象快速改变方向时查看十分之一、百分之一或千分之几秒。时间步长的大小通常称为模拟的“delta t”，而你的太大了。

对于行星体，您现在所做的就像试图通过每隔几个月模拟地球的位置来模拟地球，同时假设它沿直线运动。您需要每隔几分钟而不是每隔几个月实际模拟一次它的位置。

在我们的示例中，我们可能需要我们的重力方向基于我们当前位置和这一帧结束时的位置之间的平均值。

您可以这样做，但简单地减少计算的“delta t”会更容易。那么帧开始和结束的方向之间的差异要小得多。

一旦你解决了这个问题，你就可以使用更多的技术。例如，您可以检测帧之间的位置变化太大，然后返回并以更小的时间步重做计算。如果位置几乎没有变化，则增加时间步长。

一旦您对进行了排序，您可以在物理模拟中使用许多更高级的技术，但我首先要让基本的时间步进真正扎实。更高级的技术本质上是对您“对时间步长的变化进行更智能的插值”的想法的变体——您在这里走在正确的轨道上，但您应该在跑步之前先走。

Answer 2

我将从一个几乎与您一直在使用的 Euler-Cromer 积​​分一样简单但明显更准确的技术开始。这就是越级技术。这个想法很简单：位置和速度彼此保持半个时间步长。

初始状态在时间 t₀ 具有位置和速度。要获得半步偏移，第一步需要一个特殊情况，其中速度使用间隔开始时的加速度提前半个时间步，然后位置提前一整步。在第一次特殊情况之后，代码就像您的 Euler-Cromer 积​​分器一样工作。

在伪代码中，算法看起来像

void calculate_accel (orbiting_body_collection, central_body) {
    foreach (orbiting_body : orbiting_body_collection) {
        delta_pos = central_body.pos - orbiting_body.pos;
        orbiting_body.acc =
            (central_body.mu / pow(delta_pos.magnitude(),3)) * delta_pos;
    }
}

void leapfrog_step (orbiting_body_collection, central_body, delta_t) {
    static bool initialized = false;
    calculate_accel (orbiting_body_collection, central_body);
    if (! initialized) {
        initialized = true;
        foreach orbiting_body {
            orbiting_body.vel += orbiting_body.acc*delta_t/2.0;
            orbiting_body.pos += orbiting_body.vel*delta_t;
        }
    }
    else {
        foreach orbiting_body {
            orbiting_body.vel += orbiting_body.acc*delta_t;
            orbiting_body.pos += orbiting_body.vel*delta_t;
        }
    }
}

请注意，我已将加速度添加为每个轨道物体的场。这是保持算法与您的相似的临时步骤。另请注意，我将加速度的计算移到了它自己的单独函数中。这不是一个临时步骤。这是迈向更高级集成技术的第一步。

下一个重要步骤是撤消临时添加的加速度。加速度正确地属于积分器，而不是身体。另一方面，加速度的计算属于问题空间，而不是积分器。您可能想要添加相对论修正、太阳辐射压力或行星与行星的引力相互作用。积分器应该不知道计算这些加速度的内容。函数calculate_accels是一个被积分器调用的黑盒子。

不同的积分器对于何时需要计算加速度有非常不同的概念。有些存储最近加速度的历史，有些需要额外的工作空间来计算某种平均加速度。有些人对速度做同样的事情（保留历史记录，有一些速度工作空间）。一些更高级的集成技术在内部使用多种技术，从一种技术切换到另一种技术，以在准确性和 CPU 使用率之间提供最佳平衡。如果你想模拟太阳系，你需要一个极其精确的积分器。 （而且你需要远离浮标。即使双打也不足以实现高精度的太阳系集成。使用浮标，超过 RK4 没有多大意义，甚至可能没有跨越。）

正确区分属于谁的内容（集成者与问题空间）可以优化问题域（添加相关性等），并可以轻松切换集成技术，以便您可以评估一种技术与另一种技术。

Answer 3

所以我找到了一个解决方案，它可能不是最聪明的，但它有效，并且在阅读了 Eric 的答案和阅读了 marcus 的评论之后，我想到了它，你可以说它是两者的结合：

这是新代码：

foreach (ExtTerBody OtherObject in UniverseController.CurrentUniverse.ExterTerBodies.Where(x =>  x != this))
{
double massOther = OtherObject.Mass;

double R = Vector2Math.Distance(Position, OtherObject.Position);

double V = (massOther) / Math.Pow(R,2) * Time.DeltaTime;

float VRmod = (float)Math.Round(V/(R*0.001), 0, MidpointRounding.AwayFromZero);
if(V > R*0.01f)
{
    for (int x = 0; x < VRmod; x++)
    {
        EulerMovement(OtherObject, Time.DeltaTime / VRmod);
    }
}
else
    EulerMovement(OtherObject, Time.DeltaTime);

}

public void EulerMovement(ExtTerBody OtherObject, float deltaTime)
    {

            double massOther = OtherObject.Mass;

            double R = Vector2Math.Distance(Position, OtherObject.Position);

            double V = (massOther) / Math.Pow(R, 2) * deltaTime;

            Vector2 NonNormTwo = (OtherObject.Position - Position).Normalized() * V;

            Vector2 NonNormDir = Velocity + NonNormTwo;
            Velocity = NonNormDir;



            //Debug.WriteLine("Velocity=" + Velocity);
            Position += Velocity * deltaTime;
    }

解释一下：

我得出的结论是，如果问题是卫星在一帧中的速度太大，那么为什么不将它分成多个帧呢？所以这就是“它”现在所做的。

当卫星的速度超过当前半径的 1% 时，它会将计算分成多个部分，使其更加精确。这当然会在高速工作时降低帧率，但对于像这样的项目。

仍然非常欢迎不同的解决方案。我可能会调整触发量，但最重要的是它可以工作，然后我可以担心让它更流畅！

感谢所有看过的人，以及所有帮助自己找到结论的人！ :) 人们能提供这样的帮助真是太棒了！