球体碰撞算法总体速度保持递增(C++，OpenGL）

我的"装满球体粒子的盒子"项目遇到了问题。我创建了一个spheres[number_of_spheres][position，radius，velocity..]数组。边界框（与正交轴平行的边）。与长方体碰撞的球体的速度在碰撞的维度上反转。这还可以。对于球体到球体的碰撞，它非常简单，我通过比较距离和半径的总和来检测碰撞。现在我试图让它们简单地交换碰撞轴上的速度。但模拟运行了一段时间，然后速度不断增加，直到失控。

bool CollisionDetect(int i, int j)
{
    if(i==j)
    {
        return false;
    }
    float xx = (spheres[i][0]-spheres[j][0])*(spheres[i][0]-spheres[j][0]);
    float yy = (spheres[i][1]-spheres[j][1])*(spheres[i][1]-spheres[j][1]);
    float zz = (spheres[i][2]-spheres[j][2])*(spheres[i][2]-spheres[j][2]);
    if( (xx + yy + zz) <= (spheres[i][3] + spheres[j][3])*(spheres[i][3] + spheres[j][3]) ) 
    {
        return true;
    }
    else
    {
        return false;
    }
}
void CollisionSolve(int i, int j)
{
    float m1,m2,m21;
    Vec3 vel1,vel2,pos1,pos2;   //spheres info
    Vec3 nv1n,nv1b,nv1t; 
    Vec3 nv2n,nv2b,nv2t;
    Vec3 N,T,B,X,Y,Z;           //collision plane and world orthonormal basis
    X.x=1;
    X.y=0;
    X.z=0;
    Y.x=0;
    Y.y=1;
    Y.z=0;
    Z.x=0;
    Z.y=0;
    Z.z=1;
    pos1.x = spheres[i][0];
    pos2.x = spheres[j][0];
    pos1.y = spheres[i][1];
    pos2.y = spheres[j][1];
    pos1.z = spheres[i][2];
    pos2.z = spheres[j][2];
    vel1.x = spheres[i][4];
    vel2.x = spheres[j][4];
    vel1.y = spheres[i][5];
    vel2.y = spheres[j][5];
    vel1.z = spheres[i][6];
    vel2.z = spheres[j][6];
    m1 =  spheres[i][8];
    m2 =  spheres[j][8];                        //mass (for later)
    N = minus(pos2,pos1);                       //get N vector (connecting centers of spheres)
    N = Normalize(N);
    T=X;
    B = crossProduct(N,T);                      //find first perpendicular axis to N 
    if (B.x==0 && B.y==0 && B.z==0)             //then vector parallel to X axis - 
    {
        T=Y;                                    //try Y axis
        B = crossProduct(N,T);  
    }
    T = crossProduct(N,B);                      //find second perpendicular axis to N   
    T = Normalize(T);
    B = Normalize(B);

    if (simplespherecollision)
    {
        nv1n = projectUonV (vel1 , N);
        nv2n = projectUonV (vel2 , N);
        vel1 = minus (vel1,nv1n);
        vel2 = minus (vel2,nv2n);
        vel1 = plus (vel1,nv2n);
        vel2 = plus (vel2,nv1n);
                    //simply switch speed (for test)
    }
    /*/---THIS IS COMMENTED OUT (FIRST METHOD USED - DIDNT WORK)--------------------------------------------
    nv1n = projectUonV (vel1 , N);
    nv2n = projectUonV (vel2 , N);
    nv1t = projectUonV (vel1 , T);
    nv2t = projectUonV (vel2 , T);
    nv1b = projectUonV (vel1 , B);
    nv2b = projectUonV (vel2 , B);              //project velocities on new orthonormal basis
    vel1 = plus(nv1t , plus(nv1b , nv2n));      //project velocities back to world basis 
    vel2 = plus(nv2t , plus( nv2b , nv1n));     //by adding the sub vectors with swiched Xn
    /*/----------------------------------------------------------------------
    spheres[i][4] = vel1.x;
    spheres[i][5] = vel1.y;
    spheres[i][6] = vel1.z;
    spheres[j][4] = vel2.x;
    spheres[j][5] = vel2.y;
    spheres[j][6] = vel2.z;                     //reasign velocities to spheres
}
}

这是Vec3结构（以防万一）

struct Vec3
{
   float x, y, z;
};
Vec3 crossProduct(const Vec3& v1, const Vec3& v2)
{
    Vec3 r;
    r.x = (v1.y*v2.z) - (v1.z*v2.y);
    r.y = (v1.z*v2.x) - (v1.x*v2.z);
    r.z = (v1.x*v2.y) - (v1.y*v2.x);
    return r;
}
Vec3 minus(const Vec3& v1, const Vec3& v2) 
{
    Vec3 r;
    r.x = v1.x - v2.x;
    r.y = v1.y - v2.y;
    r.z = v1.z - v2.z;
    return r;
}
Vec3 plus(const Vec3& v1, const Vec3& v2) 
{
    Vec3 r;
    r.x = v1.x + v2.x;
    r.y = v1.y + v2.y;
    r.z = v1.z + v2.z;
    return r;
}
double dotProduct(const Vec3& v1, const Vec3& v2) 
{
    return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
}
Vec3 scale(const Vec3& v, double a) 
{
    Vec3 r;
    r.x = v.x * a;
    r.y = v.y * a;
    r.z = v.z * a;
    return r;
}
Vec3 projectUonV(const Vec3& u, const Vec3& v) 
{
    Vec3 r;
    r = scale(v,dotProduct(u,v));
    return r;
}
int distanceSquared(const Vec3& v1, const Vec3& v2) 
{
    Vec3 delta = minus(v2, v1);
    return dotProduct(delta, delta);
}
Vec3 Normalize (const Vec3& v)
{
    Vec3 r;
    r=v;
    int mag = sqrt(dotProduct(v,v));
    r.x = v.x/mag;
    r.y = v.y/mag;
    r.z = v.z/mag;
    return r;
}

这是我的渲染函数，它绘制球体和边界框，并调用其他函数

void Render()
{    
  //CLEARS FRAME BUFFER ie COLOR BUFFER& DEPTH BUFFER (1.0)
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);  // Clean up the colour of the window
                                                       // and the depth buffer
  glMatrixMode(GL_MODELVIEW); 
  glLoadIdentity();
  glTranslatef(-MAX_X/2,-MAX_Y/2,-MAX_Z*4);
  glTranslatef (dx,0,dz);
  //------------------------------------------------------------------------------------------------
    //-----------------------------BOUNDING BOX----------------------------------------------------
    glPushMatrix();
            glTranslatef(MAX_X/2,MAX_Y/2,MAX_Z/2);
            glColor4f(0.9,0.6,0.8,0.3);
            glScalef(MAX_X,MAX_Y,MAX_Z);
            glutSolidCube(1);
            glutWireCube(1);
    glPopMatrix();
  //-----------------------------WALL COLLISION DETECTION-----ALWAYS ON--------------------------------
    for (int i = 0; i<PART_NUM; i++)
    {
        if(spheres[i][0] > (MAX_X - spheres[i][3]) && spheres[i][4] > 0)
        {
            spheres[i][4] = spheres[i][4]*-1;
            wallcollcount++;
        }
        if(spheres[i][1] > (MAX_Y - spheres[i][3]) && spheres[i][5] > 0)
        {
            spheres[i][5] = spheres[i][5]*-1;
            wallcollcount++;
        }
        if(spheres[i][2] > (MAX_Z - spheres[i][3]) && spheres[i][6] > 0)
        {
            spheres[i][6] = spheres[i][6]*-1;
            wallcollcount++;
        }
        if(spheres[i][0] < spheres[i][3] && spheres[i][4] < 0)
        {
            spheres[i][4] = spheres[i][4]*-1;
            wallcollcount++;
        }
        if(spheres[i][1] < spheres[i][3] && spheres[i][5] < 0)
        {
            spheres[i][5] = spheres[i][5]*-1;
            wallcollcount++;
        }
        if(spheres[i][2] < spheres[i][3] && spheres[i][6] < 0)
        {
            spheres[i][6] = spheres[i][6]*-1;
            wallcollcount++;
        }
//--------------------------------------Sphere ColDit -------------------------
        if (spherecollision || simplespherecollision)
        {
            for(int j=i+1; j<PART_NUM; j++)
            {
                if(CollisionDetect(i,j))
                {
                    spherecollcount++;
                    CollisionSolve(i,j);
                }
            }
        }
//-----------------------------------------------DRAW--------------------------------
        glPushMatrix();
            glTranslatef(spheres[i][0],spheres[i][1],spheres[i][2]);
            glColor3f( (float) (i+1)/(PART_NUM)   , 1-(float)(i+1)/(PART_NUM)  ,  0.5);
            glutSolidSphere(spheres[i][3], 18,18);
        glPopMatrix();
    }
 //---------------------------------------------------------------------------------------------------------------/
  glutSwapBuffers();             // All drawing commands applied to the 
                                 // hidden buffer, so now, bring forward
                                 // the hidden buffer and hide the visible one
}