小球运动及碰撞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平面上界，即

说明在 $\Delta t=d/v_{z}$ 时刻之前，小球刚好接触到z平面，碰撞后，小球x，y轴方向速度不变，z轴方向速度反向，因此碰撞后，

$z'=z_{ub}-r-v_{z}*\Delta t$

$v_{z}=-v_{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仿真模型相关推荐

js实现简单的小球与边框碰撞反弹改变运动方向及颜色，并且继续运动的特效
js实现简单的小球与边框碰撞反弹改变运动方向及颜色,并且继续运动的特效 (代码可以直接复制使用,只需要把body中的div的id换成对应的就行,css中可以设置小球的大小和初始位置,修改小球大小之后需 ...
QT5-实现小球运动碰撞
QT实现小球运动碰撞反弹思路结构设置运动速度画出小球 BoundingRect和Shape函数运动和碰撞反弹编程工具为vs2010+qt 5.1.2 思路由于第一次刚使用qt,所以就写写 ...
canvas小球连线碰撞效果 html+css+js
先看效果(完整代码在底部): canvas小球连线碰撞完整源代码 html+css+js 这是我的B站地址~热烈欢迎实现(看注释,可一步一步跟着实现): 0:定义标签: <h1>北极光之 ...
Windows Store App JavaScript 开发：小球运动示例
通过前面内容的学习,相信读者已经对开发基于JavaScript的Windows应用商店应用有了一定的了解,本小节通过一个小球运动的示例来介绍如何新建一个JavaScript的Windows应用商店项目 ...
Unity空间与运动（中山大学3D游戏作业3）
Unity空间与运动(中山大学3D游戏作业3) 目录 Unity空间与运动(中山大学3D游戏作业3) 一.程序验证物体运动的本质三种方法实现抛物线运动实现太阳系二.牧师与恶魔游戏代码仓库:h ...
简单html js 特效,Js实现简单的小球运动特效
废话不多说了,直接给大家贴js代码了,具体代码如下所示: //定义局部变量 var directX=;//定义x轴方向 var directY=;//定义y轴方向 var ballX=;//定义x轴坐 ...
css动画，小球运动
提示:文章写完后,目录可以自动生成,如何生成可参考右边的帮助文档前言随着用户的需求,以及为提升用户体验今天我们讲利用css制作动画提示:以下是本篇文章正文内容,下面案例可供参考一.定义关键帧 ...
canvs之多个小球运动
这个案例的原理及分析和上一个:canvas之小球动画是一样的:只不过在这里小球的运动速度.小球的半径以及小球的颜色我们都设置为了随机而已 canvas {border: 1px solid} < ...
[Java]桌球小游戏（小球任意角度碰撞）
import javax.swing.*; import java.awt.*; public class BallGame extends JFrame {/*继承swing里面的窗口类*///加载 ...
【Unity】控制小球运动
跟着B站教程,做了个简单的控制小球运动的场景,记录一下: 文章目录搭建场景小球运动脚本相机跟随小球运动脚本效果展示搭建场景建立地面Plane.小球Player和四面墙Wall. 小球运动脚 ...

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

小球运动及碰撞3D仿真模型相关推荐

最新文章

热门文章