0

我正在尝试创建一个太阳系模拟，但在试图找出我放置在模拟中的随机物体的初始速度向量时遇到了问题。

假设： - 我使用的是高斯重力常数，所以我的所有单位都是 AU/太阳质量/天 - 使用 x、y、z 作为坐标 - 一颗星，固定在 0,0,0。为它确定准随机质量 - 我将一颗行星放置在随机 x、y、z 坐标上，并确定其自身的准随机质量。

在我开始 nbody 循环（使用 RK4）之前，我希望行星的初始速度使其具有围绕恒星的圆形轨道。一旦模拟开始，其他放置的行星当然会拉上它，但我想让它有机会拥有一个稳定的轨道......

所以，最后，我需要有一个行星的初始速度矢量 (x,y,z)，这意味着它在 1 个时间步长后将围绕恒星形成一个圆形轨道。

帮助？几个星期以来我一直在反对这个问题，但我认为我还没有任何合理的解决方案......

4

3 回答 3

4

如果你假设恒星的M质量比所有行星的总质量大得多，那就很简单了sum(m[i])。这简化了问题，因为它允许您将星星固定到坐标系的中心。此外，更容易假设所有行星的运动是共面的，这进一步将问题的维数降低到 2D。

  1. r[i]首先确定给定半径矢量（轨道半径）大小的圆形轨道速度的大小。由于上述假设，它仅取决于恒星的质量：v[i] = sqrt(mu / r[i])，其中mu是恒星的标准引力参数，mu = G * M。

  2. phi[i]通过从 中均匀采样来选择一个随机轨道相位参数[0, 2*pi)。那么行星在笛卡尔坐标中的初始位置为：
    x[i] = r[i] * cos(phi[i])
    y[i] = r[i] * sin(phi[i])

  3. 对于圆形轨道，速度矢量总是垂直于径向矢量，即它的方向是phi[i] +/- pi/2（+pi/2逆时针（CCW）旋转和-pi/2顺时针旋转）。让我们以CCW旋转为例。行星速度的笛卡尔坐标为：
    vx[i] = v[i] * cos(phi[i] + pi/2) = -v[i] * sin(phi[i])
    vy[i] = v[i] * sin(phi[i] + pi/2) = v[i] * cos(phi[i])

z[i] = 0这很容易通过添加和扩展到共面 3D 运动vz[i] = 0，但这没有任何意义，因为在 Z 方向上没有力，因此z[i]并且vz[i]永远保持等于0（即您将解决完整 3D 空间的 2D 子空间问题）。

通过完整的 3D 模拟，每个行星在随机倾斜的初始轨道上移动，可以这样工作：

  1. 此步骤等于 2D 案例中的步骤 1。

  2. 您需要在单位球体表面上选取一个初始位置。有关如何以统一随机方式执行此操作的示例，请参见此处。然后将单位球体坐标按 的大小缩放r[i]。

  3. 在 3D 情况下，行星速度所在的不是两个可能的垂直向量，而是一个完整的切平面。切平面的法向量与半径向量和 共线dot(r[i], v[i]) = 0 = x[i]*vx[i] + y[i]*vy[i] + z[i]*vz[i]。r[i]例如，可以选择任何垂直于 的向量e1[i] = (-y[i], x[i], 0)。这会在两极产生一个空向量，因此可以选择e1[i] = (0, -z[i], y[i])。r[i]然后通过取和的叉积可以找到另一个垂直向量e1[i]：
    e2[i] = r[i] x e1[i] = (r[2]*e1[3]-r[3]*e1[2], r[3]*e1[1]-r[1]*e1[3], r[1]*e1[2]-r[2]*e1[1])。现在e1[i]和e2[i]可以通过将它们除以它们的规范来归一化：
    n1[i] = e1[i] / ||e1[i]||
    n2[i] = e2[i] / ||e2[i]||
    where ||a|| = sqrt(dot(a, a)) = sqrt(a.x^2 + a.y^2 + a.z^2)。omega既然你在切平面上有一个正交向量基础，你可以选择一个随机角度[0, 2*pi)并将速度向量计算为v[i] = cos(omega) * n1[i] + sin(omega) * n2[i]，或作为笛卡尔分量：
    vx[i] = cos(omega) * n1[i].x + sin(omega) * n2[i].x
    vy[i] = cos(omega) * n1[i].y + sin(omega) * n2[i].y
    vz[i] = cos(omega) * n1[i].z + sin(omega) * n2[i].z.

请注意，通过构造，步骤 3 中的基础取决于半径矢量，但这并不重要，因为omega添加了随机方向 ( )。

至于单位的选择，在模拟科学中，我们总是倾向于将事物保持在自然单位中，即所有计算量都是无量纲的并且保持在[0, 1]或至少在 1-2 个数量级之内的单位，因此有限浮动的全分辨率可以使用点表示。如果您将恒星质量以太阳质量为单位，距离以 AU 为单位，时间以年为单位，那么对于在1AU 处的类地行星围绕类太阳恒星，轨道速度的大小为2*pi(AU/yr)，半径矢量的大小为1(AU)。

于 2013-02-13T10:18:49.420 回答
1

只要让向心加速度等于重力加速度。

m 1 v 2 / r = G m 1 m 2 / r 2

v = sqrt( G m 2 / r )

当然，恒星质量 m 2必须比行星质量 m 1大得多，否则你就没有真正的一体问题。

设置物理问题时，单位是一件很痛苦的事情。我花了几天时间解决以秒为单位的错误与时间步长单位。你对 AU/Solar Masses/Day 的选择是完全疯狂的。先解决这个问题。

并且，请记住，计算机本身的精度有限。nbody 模拟会累积积分误差，因此经过一百万步或十亿步后，无论步长如何，您肯定不会有一个圆圈。我对那个数学知之甚少，但我认为稳定的 n 体系统通过吸收微小变化的共振来保持自身稳定，无论是由附近的恒星还是由 FPU 引入。因此，对于稳定的 5 体问题，该设置可能工作得很好，但对于 1 体问题仍然会失败。

于 2013-02-13T02:16:01.677 回答
0

As Ed suggested, I would use the mks units, rather than some other set of units.

For the initial velocity, I would agree with part of what Ed said, but I would use the vector form of the centripetal acceleration:

m1v2/r r(hat) = G m1 m2 / r2 r(hat)

Set z to 0, and convert from polar coordinates to cartesian coordinates (x,y). Then, you can assign either y or x an initial velocity, and compute what the other variable is to satisfy the circular orbit criteria. This should give you an initial (Vx,Vy) that you can start your nbody problem from. There should also be quite a bit of literature out there on numerical recipes for nbody central force problems.

于 2013-02-13T05:29:37.447 回答