《模拟分子运动》由会员分享，可在线阅读，更多相关《模拟分子运动（25页珍藏版）》请在金锄头文库上搜索。
1、弹性碰撞的模拟实验Author:hank考虑以下系统，在由多边形围成的平面区域内，有小球作无摩擦水平运动。在此过程中，小球将和多边形边界发生碰撞或者与其他球碰撞，而使自身的运动速度产生变化。本实验假定，所有的碰撞均是弹性的，并且具有相同的恢复系数e.本此试验以OpenGL为平台经行编程，由于时间仓促，水平较低，尚有许多不足，敬请见谅理论背景：对心弹性碰撞：假设现在有两个物体发生了平面内的对心碰撞，质量分别为m1,m2；碰撞前速度以v1,v2表示,碰后变为u1,u2.以经过A,B质心的直线作为x轴，垂直纸面向外为z轴正方向，建立右手直角坐标系。把v1,v2,u1,u2分解到坐标轴上，根据碰撞前后
2、的动量守恒，可得方程：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.对于小球
3、与一条边（可视为一堵墙壁）碰撞：点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.小球与两条边同时碰撞这里只考虑两边正交的情况，如图所示：这种情况下，小球与两墙壁的作用可以认为互不影响，把入射速度分解到与两边正交的方向上，反射后速度的收缩系数（
4、即与原来速度的比率）相同，因此小球将反方向返回，速度变为原来的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碰撞
5、条件：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/运行说明：由于是新手上路，还没能把一些交互手段用上，这是最大遗
6、憾。你可以通过主文件设置有关参数进行实验。通过设置头文件vector3f.h中的e值，即弹性碰撞的回复系数，可以观察不同情况下的碰撞行为。如果把e值调至大于1.0，回发现小球的速度会在很短时间内增殖很大，随之小球将脱离墙壁。这一部分因为程序本身设计不够完善，但是在一定程度上依赖于计算机的响应时间。可想如果一个物体速度极快，以至时间本身对于它也不能很好描述，那么这个物体就会脱离我们现在的世界，也许将进行时空穿梭。谁知道？附录文件：hant_main.cpp#includevector3f.h/定义矢量运算#includewall.h/关于墙壁的类#includeball.h/可通过在这个文件中设
7、置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,
8、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
9、 = 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;i
10、nit();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.Dr
11、aw();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,
12、(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(voi
