模拟分子运动.doc

弹性碰撞的模拟实验Author:hank考虑以下系统，在由多边形围成的平面区域内，有小球作无摩擦水平运动。在此过程中，小球将和多边形边界发生碰撞或者与其他球碰撞，而使自身的运动速度产生变化。本实验假定，所有的碰撞均是弹性的，并且具有相同的恢复系数e.本此试验以OpenGL为平台经行编程，由于时间仓促，水平较低，尚有许多不足，敬请见谅理论背景：对心弹性碰撞：假设现在有两个物体发生了平面内的对心碰撞，质量分别为m1,m2；碰撞前速度以v1,v2表示,碰后变为u1,u2.以经过A,B质心的直线作为x轴，垂直纸面向外为z轴正方向，建立右手直角坐标系。把v1,v2,u1,u2分解到坐标轴上，根据碰撞前后的动量守恒，可得方程：m1*v1x+m2*v2x=m1*u1x+m2*u2x1式m1*v1y+m2*v2y=m1*u1y+m2*u2y2式又在y轴上，物体未受到冲击力的作用，而且不计摩擦，以此速度在y轴上的分量，碰撞前后未发生变化。即：u1y=v1y3式u2y=y2y4式在x轴方向上，可根据实验补充方程，e=(u2x-u1x)/(v1x-v2x)5式与1式联立，可得到：u1x=(m1-e*m2)/(m1+m2)*v1x+(e+1)*m2/(m1+m2)*v2x6式u2x=(e+1)*m1/(m1+m2)*v1x+(m2-e*m1)/(m1+m2)*v2x7式程序中的碰撞分一下三种：1.对于小球与一条边（可视为一堵墙壁）碰撞：点p1和点p2分别是这条边的起点与终点，设向量a表示p1p2方向的单位向量，向量b在xy平面内与之正交，满足：b=ka（k 是 z 轴的单位向量）于是：u=(v*a)*-e*(v*b)*b碰撞条件的判断：a.小球中心到直线p1p2的距离d=r;其中，D = (op1&p1p2) / |p1p2|b.球和墙壁有相对挤压的趋势，即：（v*b）*b指向斜边的一侧可知 （v*b）*b*b0;进而 v*b0;2.小球与两条边同时碰撞这里只考虑两边正交的情况，如图所示：这种情况下，小球与两墙壁的作用可以认为互不影响，把入射速度分解到与两边正交的方向上，反射后速度的收缩系数（即与原来速度的比率）相同，因此小球将反方向返回，速度变为原来的e倍，u = e*v;碰撞条件：同时满足与两块墙壁的碰撞条件，即可。3.球与球的碰撞此时点p1和点p2是两球的质心，同样以p1p2为轴建立坐标系，向量c是其方向上的单位向量（即c=(p1-p2)/|p1-p2|）,单位向量d在图示平面内与之正交。 由对心碰撞的理论可知：碰后两球速度：u1=(v1*d)*d+(m1-e*m2)/(m1+m2)*(v1*c)+(e+1)*m2/(m1+m2)*(v2*c)*cu2=(v2*d)*d+(e+1)*m1/(m1+m2)*(v1*c) +(m2-e*m1)/(m1+m2)*(v2*c)*c碰撞条件：a. 两球中心间距p1p2应小于它们的半径之和|p1p2| = r1+r2;b. 球p2对球p1的相对速度vr要指向p1侧，以使它们有相对挤压的趋势：Vr*p1p20;其中 vr = v2 v1;效果图：使用说明：要使用的源文件：1. hant_main.cpp2. ball.h3. vector3f.h4. wall.h文件的处理：把以下文件移到相应位置：a. glut.h文件 移动到 include/glb. glut32.lib文件 移动到 libc. glut32.dll文件 移动到 windows/system32/运行说明：由于是新手上路，还没能把一些交互手段用上，这是最大遗憾。你可以通过主文件设置有关参数进行实验。通过设置头文件vector3f.h中的e值，即弹性碰撞的回复系数，可以观察不同情况下的碰撞行为。如果把e值调至大于1.0，回发现小球的速度会在很短时间内增殖很大，随之小球将脱离墙壁。这一部分因为程序本身设计不够完善，但是在一定程度上依赖于计算机的响应时间。可想如果一个物体速度极快，以至时间本身对于它也不能很好描述，那么这个物体就会脱离我们现在的世界，也许将进行时空穿梭。谁知道？附录文件：hant_main.cpp#includevector3f.h/定义矢量运算#includewall.h/关于墙壁的类#includeball.h/可通过在这个文件中设置e值，改变恢复系数，/如果设置e值大于1.0，会出现什么情况？/#pragma comment(linker, /subsystem:windows /entry:mainCRTStartup)/运行这个程序请将压缩文件包的下列文件复制到如下路径：/ glut.h 到 include/gl/glut32.lib 到 lib/glut32.dll 到 windows/system32/菜单Reset出现问题/全局变量vector3f loc,vel;/第一个小球的初位置和初速度，其他均基于它给出GLfloat speed = 50;/vector3f x(1.0,0.0,0.0),y(0.0,1.0,0.0);/x轴，y轴上的单位向量GLfloat rds = 5.,mass = 0.6;/前两球的半径和质量GLint i,j,k;const int m = 16;const int n = m*m;/指定小球个数ball bn;/小球们vector3f a4;/定义墙壁所需的四个点wall w4;GLfloat timern; /每个小球所处的时间间隔bool IsWrite = 0; /判断是否写入数据bool IsW2 = 0;/判断多少毫秒写入数据int xx; /每多少帧读入一次数据float yy = 3000;bool Reset = 0;int SpeedState = 0;int w2f,menu;/菜单float rndnp()if(rand()%10=4)return 1.*rand()/RAND_MAX;elsereturn -1.*rand()/RAND_MAX;float rndpi()float x = 1.*rand()/RAND_MAX;return 2*pi*x;void init(void);void Display()if(Reset)/重置switch(SpeedState)case 0:break;case 1:speed += 10;break;case 2:speed -= 10;break;default:break;init();Reset = 0;SpeedState = 0;glClear(GL_COLOR_BUFFER_BIT); for(i = 0;i4;i+)wi.Draw(10,10);/画出墙壁/ 两球的碰撞for(k = 0;kn;k+)if(k!=n-1)for(j = k+1;jn;j+)bk.Cld(bj);elsebk.Cld(b0);if(IsWrite)if(IsW2)if(k=0)fouttime:xx+endl;foutspeed radius number = speedendsrdsendsnendl;fout(bk.GetVel().GetVal()=0.1)bk.Draw();for(k=0;kn;k+)/与两块墙壁碰撞for(i=0;i4;i+)if(i!=3)bk.Cld(wi,wi+1);if(i!=0)bk.Cld(wi,wi-1);elsebk.Cld(w3,w0);for(k=0;kn;k+)/与一块墙壁碰撞for(i=0;i4;i+)bk.Cld(wi);if(IsWrite)if(IsW2)foutendlendlendl;IsW2 = false;glutSwapBuffers();glutPostRedisplay(); void reshape (int w, int h)/glViewport(0, 0, (GLsizei) w, (GLsizei) h);glViewport(0, 0, SCREENMAXX, SCREENMAXX);glMatrixMode(GL_PROJECTION);glLoadIdentity ();/gluOrtho2D(0.0, SCREENMAXX, 0.0, SCREENMAXY);if(wh)glOrtho(-1*SCREENMAXX,SCREENMAXX,-1*(GLsizei)h/(GLsizei)w*SCREENMAXX,(GLsizei)h/(GLsizei)w*SCREENMAXX,-1*SCREENMAXX,SCREENMAXX);elseglOrtho(-1*(GLsizei)w/(GLsizei)h*SCREENMAXX,(GLsizei)w/(GLsizei)h*SCREENMAXX,-1*SCREENMAXX,SCREENMAXX,-1*SCREENMAXX,SCREENMAXX);glMatrixMode(GL_MODELVIEW);glLoadIdentity();void init(void)glClearColor(1.0,1.0,1.0,0.0);glShadeModel(GL_FLAT);/设置小球的参数float alpha;float k1=-600,k2=600; /用以确定墙壁坐标float k=k2-2*rds;float dx=2*k/(m-1);srand(unsigned)time(NULL);for(i=0;im;i+)for(j=0;jm;j+)alpha = rndpi();vel.Set(speed*cos(alpha),speed*sin(alpha),0);/使小球均匀分布在方形框内loc.Set(-1.0*k+j*dx,-1.0*k+i*dx,0);bi*m+j.Set(loc,vel,rds,mass,i); /设置墙壁参数a0.Set(k1,k2,0);a1.Set(k1,k1,0);a2.Set(k2,k1,0);a3.Set(k2,k2,0);for(i = 0;i4;i+)if(i!=3)/设置墙壁wi.Set(ai,ai+1);elsewi.Set(a3,a0);void mytimer(int value)IsW2 = true;glutTimerFunc(yy,mytimer,1);void mymenu(int value)switch(value)case 0:coutResetendl;Reset = 1;/重置break;case 1:fout.close();exit(0);break;case 2:coutFasterendl;Reset = 1;SpeedState = 1;break;case 3:Reset = 1;SpeedState = 2;break;default:break;void mysubmenu(int value)switch(value)case 0:coutStartendl;IsWrite = 1;break;case 1:coutStopendl;IsWrite = 0;break;default:break;int main(int argc, char* argv)glutInit(&argc,argv);glutInitDisplayMode(GLUT_DOUBLE|GLUT_RGB);glutInitWindowSize(800,800);/1200,800glutInitWindowPosition(100,40);glutCreateWindow(对心碰撞!);init();/副菜单w2f = glutCreateMenu(mysubmenu);glutAddMenuEntry(Begin,0);glutAddMenuEntry(Stop,1);glutAttachMenu(GLUT_RIGHT_BUTTON);/主菜单menu = glutCreateMenu(mymenu);glutAddMenuEntry(Reset,0);glutAddMenuEntry(Faster,2);glutAddMenuEntry(Slower,3);glutAddSubMenu(Write2File,w2f);glutAddMenuEntry(Exit,1);glutSetMenu(menu);glutAttachMenu(GLUT_RIGHT_BUTTON);glutReshapeFunc(reshape);glutDisplayFunc(Display);glutIdleFunc(Display);/glutKeyboardFunc(mykey);glutTimerFunc(yy,mytimer,1);glutMainLoop();fout.close();return 0;/*extern void APIENTRY glutDisplayFunc(void (*func)(void);extern void APIENTRY glutReshapeFunc(void (*func)(int width, int height);extern void APIENTRY glutKeyboardFunc(void (*func)(unsigned char key, int x, int y);extern void APIENTRY glutMouseFunc(void (*func)(int button, int state, int x, int y);extern void APIENTRY glutMotionFunc(void (*func)(int x, int y);extern void APIENTRY glutPassiveMotionFunc(void (*func)(int x, int y);extern void APIENTRY glutEntryFunc(void (*func)(int state);extern void APIENTRY glutVisibilityFunc(void (*func)(int state);extern void APIENTRY glutIdleFunc(void (*func)(void);extern void APIENTRY glutTimerFunc(unsigned int millis, void (*func)(int value), int value);*/*hant_main.obj : error LNK2001: unresolved external symbol _glutInitWithExit12hant_main.obj : error LNK2001: unresolved external symbol _glutCreateWindowWithExit8Debug/hant_main.exe : fatal error LNK1120: 2 unresolved externals 出现上述情况，请更换glut32.dll的版本 */ball.h/ball.h/vector3f E(0,0,1);GLfloat const e = 1;/玻璃对玻璃的恢复系数GLfloat const pi = 3.1415926;class ballprivate:vector3f Lct;/球的位置GLfloat Rds;/半径GLfloat Mass;/质量GLint uID;/球的编号vector3f Vel;/速度public:ball(vector3f loc,vector3f vel,GLfloat ra,GLfloat m,GLint id);ball();ball();vector3f GetLct();/返回位置vector3f GetVel();/返回速度GLfloat GetRds();/返回半径GLfloat GetMass();/返回质量GLint GetId();/返回编号/设置球的属性void Set(vector3f loc,vector3f vel,GLfloat ra,GLfloat m,GLint id);void SetLct(vector3f loc);/设置坐标void SetVel(vector3f v);/设置速度void SetId(GLint id);/临时改变编号bool IsCld(ball &other);/判断是否与球other碰撞bool IsCld(wall &Wall);/判断是否与墙壁Wall碰撞void Draw();/画出小球void Update(GLfloat,GLfloat*);/更新小球位置void Cld(ball &);/与另一只球相撞的行为void Cld(wall &);void Cld(wall &,wall &);/与墙壁相撞的行为;/函数定义/ball:ball():Lct(0.0,0.0,0.0),Vel(0.0,0.0,0.0)Rds = 0;Mass = 0;uID = 0;ball:ball(vector3f loc,vector3f vel,GLfloat ra,GLfloat m,GLint id)Lct = loc;Vel = vel;Rds = ra;Mass = m;uID = id;ball:ball()vector3f ball:GetLct()return Lct;vector3f ball:GetVel()return Vel;/返回速度GLfloat ball:GetRds()return Rds;/返回半径GLfloat ball:GetMass()return Mass;/返回质量GLint ball:GetId()return uID;/返回编号void ball:Set(vector3f loc,vector3f vel,GLfloat ra,GLfloat m,GLint id)Lct = loc;Vel = vel;Rds = ra;Mass = m;uID = id;void ball:SetLct(vector3f loc)Lct = loc;void ball:SetVel(vector3f v)Vel = v;void ball:SetId(GLint id)uID = id;bool ball:IsCld(ball &other)/具体依据请参看压缩包中的doc文档vector3f v21 = other.Vel - Vel;vector3f p1p2 = Lct - other.Lct;float dis = p1p2.value - (Rds+other.Rds);if(dis0.00001)return true;return false;bool ball:IsCld(wall &w)vector3f a = Lct - w.p1;vector3f b = w.p2 - w.p1;GLfloat d = (a&b).value/b.value;if(dRds&Vel*w.b0)return true;else return false;void ball:Draw()/画出小球glColor3f(sin(uID*30),sin(uID*60),sin(uID*45);GLint circle_points = 12; /GLfloat centerx = Lct.vec0;GLfloat centery = Lct.vec1;GLfloat radius = Rds;GLint i;GLdouble angle;glBegin(GL_POLYGON); float tmp = 2*pi/circle_points;for (i = 0; i circle_points; i+) angle = tmp*i; glVertex2f(centerx+radius*cos(angle), centery+radius*sin(angle); glEnd();glColor3f(0.6,0.3,0.);void ball:Update(GLfloat dt,GLfloat *timer)/coutin Update()endl;/GLfloat dt = 0.005;/这里index相当单位时间间隔，即dt/GLfloat p3 ; /Lct.GetVec(p);/GLfloat s3 ; /Vel.GetVec(s);/coutendl;/coutFor balluIDendl;/coutVelendl;/Vel.Show();/coutspeed=Vel.GetVal()endl;/for(GLint i=0;i3;i+)/更新小球位置Lct.veci += Vel.veci*dt;/*timer += dt;/coutLct:GetVel();vector3f v2=other.GetVel();vector3f a = (other.Lct-Lct).to1();vector3f b = E&a;GLfloat m1 = GetMass();GLfloat m2 = other.GetMass();GLfloat k0,k1,k2,k3,k4;k0 = m1 + m2;k1 = (m1-e*m2)/k0;k2 = (e+1)*m2/k0;k3 = (e+1)*m1/k0;k4 = (m2-e*m1)/k0;vector3f t1 = (v1*b)*b + (k1*(v1*a)+k2*(v2*a)*a;vector3f t2 = (v2*b)*b + (k3*(v1*a)+k4*(v2*a)*a;/现在两球的速度都会改变吗？v1 = t1;v2 = t2;SetVel(v1);other.SetVel(v2);/*GLfloat k0,k1,k2,k3,k4;k0 = Mass + other.Mass;k1 = (Mass-e*other.Mass)/k0;k2 = (e+1)*other.Mass/k0;k3 = (e+1)*other.Mass/k0;k4 = (Mass-e*other.Mass)/k0;vector3f t1 = (Vel*b)*b + (k1*(Vel*a)+k2*(other.Vel*a)*a;other.Vel = (other.Vel*b)*b + (k3*(Vel*a)+k4*(other.Vel*a)*a;/现在两球的速度都会改变吗？/假设所有球质量相同vector3f tmp = (1-e)*(Vel*a)+(1+e)*(other.Vel*a)*a*0.5;other.Vel = (other.Vel*b)*b + tmp;/现在两球的速度都会改变吗？Vel = (Vel*b)*b + tmp;*/与另一只球相撞的行为void ball:Cld(wall &w)if(!IsCld(w)return ;Vel = (Vel*w.a)*w.a - e*(Vel*w.b)*w.b;/与墙壁相撞的行为*/void ball:Cld(wall &w1,wall &w2)if(!(IsCld(w1)&IsCld(w2)return;/vector3f b1 = w1.GetB(), b2 = w2.GetB();Vel = -1*e*(Vel*w1.b)*w1.b - e*(Vel*w2.b)*w2.b;vector3f.h#include /the windows include file, required by all windows applications#include / included for the cos and sin functions#include /the glut file for windows operations#include / it also includes gl.h and glu.h for the openGL library calls#include #include #include #include using namespace std;#define SCREENMAXX 900#define SCREENMAXY 800ofstream fout(data.txt);/ vector3f.hclass ball;class vector3ffriend class ball;private:float vec3;float value;public:vector3f();vector3f(float,float,float);vector3f(float *p);vector3f(vector3f &p);vector3f();float * GetVec();void GetVec(float p);float GetVal();void Set(float,float,float);void Set(float *p);/定义向量运算/重载运算符void operator =(vector3f &);vector3f operator +(vector3f &);vector3f operator -(vector3f &);float operator *(vector3f &);vector3f operator *(float);vector3f operator /(float);friend vector3f operator *(float ,vector3f &);friend vector3f operator &(vector3f &,vector3f &);/单位化函数vector3f to1()