小球运动及碰撞3D仿真模型

物体运动及碰撞是一个非常基础但是很普遍物理过程，在许多场景中都会发生，小到微观世界中研究气体反应的分子碰撞理论，大到宇宙中星体运动及星系演化。不管是游戏领域，机器人领域，还是计算化学领域，都会涉及到物体运动和碰撞过程。

关于小球运动及碰撞的仿真及代码，网站资源非常多，但绝大多数以2维平面上运动及碰撞为主，很少有研究三维的，所以这次分享内容是小球在方盒中运动及碰撞3D仿真。

废话少讲，直入主题吧。

一、前期准备工作（基础物理知识+Matlab知识点）

1 运动小球描述

描述运动的小球，我们需要的参数包括：小球的半径，球心的位置，球心速度，球心加速度。当然如果要研究小球的自旋，这里还需要引入旋转轴和角速度这两个参数，在本实例中忽略球体自旋，实际在真实碰撞过程中，如果不是小球正对相撞，会导致小球发生自旋，类似乒乓球上旋和下旋，考虑到问题复杂程度，这里忽略小球自旋过程。

因此在三维下，描述小球参数有：[x, y, z, r, vx, vy, vz, ax, ay, az];忽略小球运动系统受外力影响，则描述小球参数可以简化为：[x, y, z, r, vx, vy, vz]。

2 小球运动和碰撞

忽略小球运动系统受外力影响，则小球的运动过程通过微元法可以直接写出运动方程：

位置方程：

 速度方程：

 小球碰撞过程包含两部分：a.球与边界的碰撞，b.球与球之间碰撞

a.小球与边界的碰撞

下面以z方向上边界和下边界碰撞为例，分析碰撞过程判定条件及碰撞后小球参数演变。

 如图所示，如果在某一时刻，小球越界判定条件为：

下面分析碰撞后小球位置，速度等信息演变，如果在某一时刻，小球越过z平面上界，即

说明在时刻之前，小球刚好接触到z平面，碰撞后，小球x，y轴方向速度不变，z轴方向速度反向，因此碰撞后，

b.球与球之间碰撞

小球碰撞的判定条件是球心之间距离小于小球半径之和。

 由于碰撞发生在两小球球心连线的方向上，因此，我们需要确定某一时刻下，在球心相连的方向上两个小球的速度分量，来确定两小球刚刚接触的时刻，并基于此计算小球碰撞时刻的位置和速度信息及碰撞后位置和速度信息。具体思路与边界碰撞相似。

小球之间的碰撞过程由于无外力作用，一定满足动量守恒定律。碰撞有三种不同类型：完全弹性碰撞，完全非弹性碰撞，非完全弹性碰撞。根据小球碰撞的类型不同，考虑能量守恒定律。通过碰撞原理可以求解膨胀后小球速度演变。在这里不赘述了。

    1.3 Matlab知识点（rand函数介绍）

小球初始参数信息主要通过随机数生成器来生成，因此需要使用到rand函数。rand函数用来产生随机数，内部实现是用线性同余法实现的，是伪随机数，由于周期较长，因此在一定范围内可以看成是随机的。

rand生成均匀分布的伪随机数，分布在（0~1）之间；randn生成标准正态分布的伪随机数（均值为0，方差为1）；randi 生成均匀分布的伪随机整数。

这里，我们用rand函数生成N个小球初始位置和速度，用randn生成小球的半径。

ok，有了这些基本知识，我们可以开始码代码了。

二，Matlab代码

第一步：方盒尺寸及边界设定，小球数量，位置，半径等初始参数产生;这里通过编写particlegenerate函数，基本思路是随机生成一个颗粒，与设定边界以及其他小球进行比较，确定是否超出边界或者与其他小球重叠，如果存在，则去掉，并进行下一次循环，直到获得所需颗粒数量或者超过设定循环次数。

function [pos,rad,vel] = particlegenerate(box,N,radius,sigma,type,T) %RANDPARTICLE 此处显示有关此函数的摘要 % 此处显示详细说明 if nargin < 6T=273.15; %%开尔文温度，反映粒子运动if nargin <5type='normal';%% type 可选三种颗粒分布：正态分布normal，均值分布mean，定值constantif nargin < 4sigma=0;if nargin < 3radius=1;if nargin < 2N=10;if nargin < 1box=[10 10 10];endendendendend enda0=box(1);b0=box(2);c0=box(3); pos=zeros(N,3); rad=zeros(N,1);switch typecase 'constant'rad=ones(N,1).*radius;case 'normal'rad=radius.* ones(N,1) + randn(N,1).*sigma;case 'mean'rad=radius.* ones(N,1) + rand(N,1).*sigma; end r=rad; pos(1,:)=[rand*(a0-2*r(1))+r(1) rand*(b0-2*r(1))+r(1) rand*(c0-2*r(1))+r(1)]; i=2;k=1; %%颗粒随机填充 while i<=N&&k<=100000temp=[rand*(a0-2*r(i))+r(i) rand*(b0-2*r(i))+r(i) rand*(c0-2*r(i))+r(i)];if all(vecnorm(temp-pos(1:i-1,1:3),2,2)-r(i)-r(1:i-1)>=0)pos(i,:)=temp;i=i+1;k=1;elsek=k+1;end end %%剩余空间补填颗粒 %%网格划分结合精度和算力，设定网格尺寸为netsize=[0.01 0.01 0.01]； netsize=[0.1 0.1 0.1]; %%那么整个长方体对应网格节点大小为ceiling([a0 b0 c0]./netsize)+[1 1 1]; abc=ceil([a0 b0 c0]./netsize)+[1 1 1]; ii=i; while ii<=Nfor i=1:abc(1)for j=1:abc(2)for k=1:abc(3)temp=[i j k].*[a0 b0 c0];if all((temp-[r(ii) r(ii) r(ii)])>=0)&&all((temp-[a0-r(ii) b0-r(ii) c0-r(ii)])<=0)&&all(vecnorm(temp-pos(1:i-1,1:3),2,2)-r(i)-r(1:i-1)>=0)pos(ii,:)=temp;ii=ii+1;break;elser(ii)=0;endendendendii=ii+1; endvel=rand(N,3).*radius.*T./T; vel(pos(:,1)==0,:)=[]; pos(pos(:,1)==0,:)=[]; r(r(:,1)==0,:)=[];rad=r; end

第二步，编写小球碰撞函数collision.m和长方体绘制函数plotcuboid.m；

function [v1,v2] = collision(mas1,vel1,mas2,vel2,f) %MASSENERGYEQUATION 此处显示有关此函数的摘要 % 此处显示详细说明 % 此函数暂时只支持f=1或者0情形 if f==1%% 完全弹性碰撞v1 = vel1+2*(vel2-vel1)/(1+mas2/mas1);v2 = vel2+2*(vel1-vel2)/(1+mas1/mas2); elseif f==0%% 完全非弹性碰撞v1 = (mas1*vel1+mas2*vel2)/(mas1+mas2);v2 = v1; else%% 非完全弹性碰撞 待定v1=0;v2=0; endend function PlotCuboid(startPoint,Size) %% 绘制长方体Opoint=[0 0 0;0 0 1;0 1 0;0 1 1;1 0 0;1 0 1;1 1 0;1 1 1]; cornerpoint=startPoint+Opoint.*Size;%% 定义6个平面分别对应的顶点 face=[1 2 4 3;1 2 6 5;1 3 7 5;2 4 8 6;3 4 8 7;5 6 8 7];%% 定义顶点的颜色 color=[1;2;3;4;5;6;7;8];%% 绘制图像patch('Vertices',cornerpoint,'Faces',face,'FaceVertexCData',color,'FaceColor','interp','FaceAlpha',0.5); view(3); xlabel('X'); ylabel('Y'); zlabel('Z'); title('Cubic'); fig=gcf; fig.Color=[1 1 1]; fig.Name='cuboid'; fig.NumberTitle='off';

第三步，确定仿真时间间隔以及总时长，构建元胞数组，存储所有信息，包括位置坐标，小球半径，速度，任意两个小球之间距离等信息；

clc;clear; %%最小运动间隔0.1；总仿真时间100 Dt=0.01;t=100;S=floor(t/Dt); %%初始颗粒生成，包含位置，半径，速度参数； %%限定颗粒运动区域 box=[10 10 10]; %%颗粒生成参数 N=20;r=0.3；sigma=0.05； box=[10 10 10]; N=50;r=0.3;sigma=0.05;type='normal'; %%定义碰撞因子 f=0.5；重力加速度 g=0.98； f=1;g=0; Bound=[0 box(1) 0 box(2) 0 box(3)]; [pos,rad,vel]=particlegenerate(box,N,r,sigma,type); allinfo=cell(S,3); allinfo{1,1}=pos; allinfo{1,2}=rad; allinfo{1,3}=vel; N=length(rad); allinfo{1,4}=zeros(N,N);%%Edge distance matrix

第四步，开始仿真，首先基于运动方程，确定下一时刻每个小球位置，速度，然后开始判定是否存在越界小球，如果存在，则按照上文讨论思路，给出新的位置，速度等信息；在判断是否存在小球碰撞，如果存在则按照上文讨论思路，给出碰撞后小球新的位置，速度等信息。

for k=2:Sallinfo{k,1}=allinfo{k-1,1}+allinfo{k-1,3}.*Dt;allinfo{k,2}=allinfo{k-1,2};allinfo{k,3}=allinfo{k-1,3}-ones(N,1)*[0 0 g].*Dt;allinfo{k,4}=zeros(N,N);for j=1:Nfor i=1:3d=allinfo{k,1}(j,i)+allinfo{k,2}(j)-Bound(2*i);if d>=0dt=d/allinfo{k,3}(j,i);allinfo{k,1}(j,i)=allinfo{k,1}(j,i)-allinfo{k,3}(j,i)*dt;allinfo{k,3}(j,i)=-f*allinfo{k,3}(j,i);endendendfor j=1:Nfor i=1:3d=allinfo{k,1}(j,i)-allinfo{k,2}(j)-Bound(2*i-1);if d<=0dt=d/allinfo{k,3}(j,i);allinfo{k,1}(j,i)=allinfo{k,1}(j,i)-allinfo{k,3}(j,i)*dt;allinfo{k,3}(j,i)=-f*allinfo{k,3}(j,i);endendendfor i=1:N-1for j=i+1:Nallinfo{k,4}(i,j)=norm(allinfo{k,1}(i,:)-allinfo{k,1}(j,:))-allinfo{k,2}(i)-allinfo{k,2}(j);endendallinfo{k,4}=allinfo{k,4}+tril(abs(-1+eye(N,N)))+eye(N,N);if find(allinfo{k,4}<0)[indI,indJ]=find(allinfo{k,4}<0);for i=1:length(indI)vecIJ=normr(allinfo{k,1}(indI(i),:)-allinfo{k,1}(indJ(i),:));vecVI=allinfo{k,3}(indI(i),:)*vecIJ';vecVJ=allinfo{k,3}(indJ(i),:)*vecIJ';dt=abs(allinfo{k,4}(indI(i),indJ(i)))/(abs(vecVI+vecVJ));allinfo{k,1}([indI(i),indJ(i)],:)=allinfo{k,1}([indI(i),indJ(i)],:)-allinfo{k,3}([indI(i),indJ(i)],:)*dt;endfor i=1:length(indI)vecIJ=normr(allinfo{k,1}(indI(i),:)-allinfo{k,1}(indJ(i),:));vecVI=allinfo{k,3}(indI(i),:)*vecIJ';vecVJ=allinfo{k,3}(indJ(i),:)*vecIJ';restVI=allinfo{k,3}(indI(i),:)-vecVI.*vecIJ;restVJ=allinfo{k,3}(indJ(i),:)-vecVJ.*vecIJ;massI=4*pi/3*allinfo{k,2}(indI(i)).^3;massJ=4*pi/3*allinfo{k,2}(indJ(i)).^3;[VI,VJ]=collision(massI,vecVI,massJ,vecVJ,f);allinfo{k,3}(indI(i),:)=restVI+VI.*vecIJ;allinfo{k,3}(indJ(i),:)=restVJ+VJ.*vecIJ;allinfo{k,1}([indI(i),indJ(i)],:)=allinfo{k,1}([indI(i),indJ(i)],:)+allinfo{k,3}([indI(i),indJ(i)],:)*dt;endend end

第五步，绘制仿真动态视图，根据元胞数组存储信息绘制每个时间步下小球，通过getframe函数录制每一帧下图像，采取电影动画方式绘制动态过程。

%%绘制动画 close all; figure(1) plotcuboid([0,0,0],box); light; set(gca,'xtick',[],'xticklabel',[]); set(gca,'ytick',[],'yticklabel',[]); set(gca,'ztick',[],'zticklabel',[]); axis equal; axis(Bound); axis off hold on for i=1:10:Sfor j=1:Nh(j)=plotsphere(allinfo{i,1}(j,:),allinfo{i,2}(j),'r');hold onenddrawnow;CF=getframe(gcf);imind=frame2im(CF);[imind,cm] = rgb2ind(imind,256);if i==1imwrite(imind,cm,'test.gif','gif', 'Loopcount',inf,'DelayTime',1e-6);elseimwrite(imind,cm,'test.gif','gif','WriteMode','append','DelayTime',1e-6);endfor j=1:Ndelete(h(j));end end

三，代码测试及演示

最后展示一下结果。

总结

以上是生活随笔为你收集整理的小球运动及碰撞3D仿真模型的全部内容，希望文章能够帮你解决所遇到的问题。

如果觉得生活随笔网站内容还不错，欢迎将生活随笔推荐给好友。