模拟盒子中的点质量（3D 空间）

Question

我想在一个封闭的盒子里模拟一个点质量。没有摩擦，点质量服从冲击定律。所以只有与盒子壁的弹性碰撞。程序的输出是时间、位置 (rx,ry ,rz) 和速度 (vx,vy,vz)。我使用 GNUplot 绘制轨迹。

我现在遇到的问题是，点质量从某个地方获取能量。所以他们的跳跃每次都变得更加激烈。

有人可以检查我的代码吗？

/* Start of the code */

#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

struct pointmass
{
  double m; // mass
  double r[3]; // coordinates
  double v[3]; // velocity
};

// Grav.constant
const double G[3] = {0, -9.81, 0};

int main()
{
  int Time = 0; // Duration
  double Dt = 0; // Time steps
  pointmass p0;

  cerr << "Duration: ";
  cin >> Time;

  cerr << "Time steps: ";
  cin >> Dt;

  cerr << "Velocity of the point mass (vx,vy,vz)? ";
  cin >> p0.v[0];
  cin >> p0.v[1];
  cin >> p0.v[2];

  cerr << "Initial position of the point mass (x,y,z)? ";
  cin >> p0.r[0];
  cin >> p0.r[1];
  cin >> p0.r[2];

  for (double i = 0; i<Time; i+=Dt)
  {
      cout << i << setw(10);

      for (int j = 0; j<=2; j++)
      {

        ////////////position and velocity///////////
        p0.r[j] = p0.r[j] + p0.v[j]*i + 0.5*G[j]*i*i;
        p0.v[j] = p0.v[j] + G[j]*i;

       ///////////////////reflection/////////////////

        if(p0.r[j] >= 250)
        {
          p0.r[j] = 500 - p0.r[j];
          p0.v[j] = -p0.v[j];
        }
        else if(p0.r[j] <= 0)
        {
          p0.r[j] = -p0.r[j];
          p0.v[j] = -p0.v[j];
        }
       //////////////////////////////////////////////
    }
       /////////////////////Output//////////////////
    for(int j = 0; j<=2; j++)
    {
       cout << p0.r[j] << setw(10);
    }
    for(int j = 0; j<=2; j++)
    {
       cout << p0.v[j] << setw(10);
    }
    ///////////////////////////////////////////////
 cout << endl;
 }

}

Answer 1

F = 马

a = F / m

a dt = F / m dt

a dt 是固定时间内的加速度 - 该帧的速度变化。 您将其设置为 F / m i 正如 cmets 所建议的那样，i 是错误的。它需要是一帧的持续时间，而不是到目前为止整个模拟的持续时间。

Answer 2

我和其他评论者有点担心时间循环 - 确保它代表时间的增量，而不是持续时间的增长。

不过，我认为主要问题是您正在更改速度的所有三个分量的符号 反思。

这与边界处的物理定律（线性动量和能量守恒）不一致。

要了解这一点，请考虑以下情况：您的粒子仅在 x-y 平面中移动（z 方向的速度为零）并即将在 x= L 处撞击墙壁。

碰撞看起来像这样：

墙施加在点质量上的力垂直于墙。所以平行于壁面的粒子的动量分量没有变化。

应用线性动量和动能守恒，并假设完美的弹性碰撞，你会发现

• 垂直于壁面的速度分量确实会改变符号
• 平行于墙壁的速度分量不改变符号

在三个维度上，要进行准确的模拟，您必须计算出碰撞时平行和垂直于墙壁的动量分量，并对由此产生的速度变化进行编码。

换句话说，这段代码：

///////////////////reflection/////////////////

        if(p0.r[j] >= 250)
        {
          p0.r[j] = 500 - p0.r[j];
          p0.v[j] = -p0.v[j];
        }
        else if(p0.r[j] <= 0)
        {
          p0.r[j] = -p0.r[j];
          p0.v[j] = -p0.v[j];
        }
       //////////////////////////////////////////////

没有正确模拟反射的物理特性。在这里修复它是要做什么的大纲：

• 将反射检查从 x、y、z 坐标上的循环中取出（但仍在时间循环内）
• 需要检查所有六面墙的碰撞情况， 根据法向量到墙壁的方向。

例如对于由 X=250, 0

因此，在从这两个墙壁反射时，速度的 X 分量会改变符号，因为它垂直（垂直）于墙壁，但 Y 和 Z 分量不会改变符号，因为它们平行于墙壁。

在立方体左侧的顶壁和底壁（常数 Y）以及前后壁（常数 Z）应用类似的考虑因素，作为计算这些表面法线的练习。

最后，您不应该在反射时更改位置矢量分量的符号，而应该只更改速度矢量。而是在给定新速度的情况下重新计算位置向量的下一个值。

Answer 3

好的，所以有一些问题。其他人指出需要使用Dt 而不是i 进行集成步骤。

但是，您说反射和能量守恒存在问题是正确的。我在下面添加了一个明确的轨道。

请注意，除了能量问题之外，反射的分量计算实际上很好。

问题在于，在反射过程中，重力引起的加速度会发生变化。在粒子撞击地板的情况下，它获得的动能等于它一直下落时的动能，但新位置的势能更高。因此，能量将增加恰好是地板和新位置之间势能差的两倍。从屋顶反弹会产生相反的效果。

如下所述，曾经的策略是计算反射的实际时间。然而，实际上直接使用能量工作更简单，也更强大。但是请注意，虽然下面的简单能量版本保证了速度和位置是一致的，但它实际上并没有正确的位置。对于大多数可能并不重要的目的。如果你真的需要正确的位置，我认为我们需要解决反弹时间。

/* Start of the code */

#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

struct pointmass
{
    double m; // mass
    double r[3]; // coordinates
    double v[3]; // velocity
};

// Grav.constant
const double G[3] = { 0, -9.81, 0 };


int main()
{
    // I've just changed the initial values to speed up unit testing; your code worked fine here.
    int Time = 50; // Duration
    double Dt = 1; // Time steps
    pointmass p0;
    p0.v[0] = 23;
    p0.v[1] = 40;
    p0.v[2] = 15;

    p0.r[0] = 100;
    p0.r[1] = 200;
    p0.r[2] = 67;


    for (double i = 0; i<Time; i += Dt)
    {
        cout << setw(10) << i << setw(10);

        double energy = 0;
        for (int j = 0; j <= 2; j++)
        {
            double oldR = p0.r[j];
            double oldV = p0.v[j];
            ////////////position and velocity///////////

            p0.r[j] = p0.r[j] + p0.v[j] * Dt + 0.5*G[j] * Dt*Dt;
            p0.v[j] = p0.v[j] + G[j] * Dt;

            ///////////////////reflection/////////////////
            if (G[j] == 0)
            {
                if (p0.r[j] >= 250)
                {
                    p0.r[j] = 500 - p0.r[j];
                    p0.v[j] = -p0.v[j];
                }
                else if (p0.r[j] <= 0)
                {
                    p0.r[j] = -p0.r[j];
                    p0.v[j] = -p0.v[j];
                }
            }
            else
            {
                // Need to capture the fact that the acceleration switches direction relative to velocity half way through the timestep.
                // Two approaches, either 
                //    Try to compute the time of the bounce and work out the detail.
                // OR
                //    Use conservation of energy to get the right speed - much easier!
                if (p0.r[j] >= 250)
                {
                    double energy = 0.5*p0.v[j] * p0.v[j] - G[j] * p0.r[j];
                    p0.r[j] = 500 - p0.r[j];
                    p0.v[j] = -sqrt(2 * (energy + G[j] * p0.r[j])); 
                }
                else if (p0.r[j] <= 0)
                {
                    double energy = 0.5*p0.v[j] * p0.v[j] - G[j] * p0.r[j];
                    p0.r[j] = -p0.r[j];
                    p0.v[j] = sqrt(2*(energy + G[j] * p0.r[j]));
                }
            }
            energy += 0.5*p0.v[j] * p0.v[j] - G[j] * p0.r[j];
        }
        /////////////////////Output//////////////////
        cout << energy << setw(10);
        for (int j = 0; j <= 2; j++)
        {
            cout << p0.r[j] << setw(10);
        }
        for (int j = 0; j <= 2; j++)
        {
            cout << p0.v[j] << setw(10);
        }
        ///////////////////////////////////////////////
        cout << endl;
    }
}