粒子碰撞

Question

我正在尝试编写一个代码，它给出一维盒子中 N 个粒子的质量、速度、位移和半径，将以 0.03 秒的小增量返回它们的位置。 该程序符合。但是不会产生正确的值。

#include <iostream>
#include <cmath>
#include <fstream>
#include <vector>


using namespace std;



struct particle {
    double x;   //position of particle
    double im;  //inverse mass
    double v;   //velocity of particle
    double r;   //radius of particle
};

void bonkers(vector<particle> &array, int NN, double t, double ddt, int T);
double time(vector<particle> &a, int NN, int &whichn);
void collision(particle &a, particle &b);
void position(particle &a, double ddt);

void position(particle &a, double ddt)//position of particles
{
    a.x += ddt * a.v; //position = t * velocity
}

void collision(particle &a, particle &b)//velocities after collision
{
    double realativeV = a.v - b.v;
    double af = a.im / (a.im + b.im);
    double bf = b.im / (a.im + b.im);
    a.v -= 2.0 * af * realativeV;
    b.v += 2.0 * bf * realativeV;
}

double time(vector<particle> &a, int NN, int &whichn)//Finds time of first collision.
{
    double dt = 1e100; //really large number
    for (int n = 0; n < NN - 1 ; n++) 
    {
        double RelativeV = a[n].v - a[n + 1].v;
        if (RelativeV > 0)//this means dt won't be 0, no negative times
        {
            double Collisiont = ((a[n + 1].x - a[n + 1].r) - (a[n].x + a[n].r))/ RelativeV;//gives time of collision
            //time of collision worked out as time when the distance between the two balls centre is equal to their combined radi
            if (Collisiont < dt)//finds smallest possible time of collision
            {
                dt = Collisiont; //returns time of first collision
                whichn = n; //gives which particles collide. Therefore which particles velocities need to change
            }
        }
    }
    return dt;
}

void bonkers(vector<particle> &array, int NN, double t, double ddt, int T)
{

    ofstream myfile;
    myfile.open("example.txt");//sends all information to file called example.txt
    int whichn = 0;
    double k;
    double l = 0;
    for (int i = 0; i <= T; i++)//set arbitrary number of collision i.e. T
    {
        k = t;//after every collision k reset to 0 as time function works from 0
        while( k <= time(array, NN, whichn) )
        {
            myfile << l; //prints the time
            k += ddt; //incriments time
            l += ddt; //increments time, but doesn't reset time so gives total time

            for (int n = 1; n < NN -1 ; n++)
            {
                position(array[n], ddt);//calculates new position
                myfile << "\t" << array[n].x;//prints the positions of all the particles at said time
            } 
            myfile << endl;
        }
        collision(array[whichn], array[whichn + 1]);//asings new velocities to the two particles which collided
    }
    myfile.close();

}

int main()
{
    int N; //number of balls
    double m; //mass of balls
    int T = 50; //number of collisions
    double t = 0.0; //starts from 0 time
    double ddt = 0.03; //invrements by 0.03 seconds each time

    cout << "Input number of balls" << endl;
    cin >> N;

    vector<particle> a(N + 2); //extra two for the two walls, wall1 at 0 and wall2 at 20

    for (int k = 1; k <= N; k++)
    {
        cout << "Input mass of ball " << k << endl;
        cin >> m;
        a[k].im = 1.0 / m;  //finds inverse mass
        cout << "Input position of ball " << k << " between 0 and 20" << endl;
        cin >> a[k].x;
        cout << "Input speed of ball " << k << endl;
        cin >> a[k].v;
        cout << "Input radius of ball " << k << endl;
        cin >> a[k].r;
    }

    //asign wall variables
    a[0].im = 0;    //infinte mass, so walls don't move
    a[0].r = 0;
    a[0].x = 0; //wall1 at 0
    a[0].v = 0;
    a[N + 1].im = 0;    
    a[N + 1].r = 0;
    a[N + 1].x = 20;    //wall2 at 20
    a[N + 1].v = 0;

    bonkers(a, N + 2, t, ddt, T);   
    return 0;
}

当你只有一个球时有一个非常明显的问题。半径、质量和速度均为 1。 由于某种原因，球的速度在 10.03 位置发生了变化，而不是在 20 处反弹。

我相信问题出在代码的“时间”函数部分。 通过在 bonkers 函数中将“k”而不是“l”输出到 example.txt，这个问题变得更加明显。

非常感谢:)

Answer 1

问题是这一行。

while( k <= time(array, NN, whichn) )

您在 while 循环的每次迭代中重新计算碰撞时间，并增加经过的时间 k，然后比较它们。直到事件发生的时间和经过的时间没有意义进行比较。您可以计算碰撞时间并增加时间k，直到您到达那里，或者您可以重新计算碰撞时间并在每次迭代中与 0 进行比较。