物體運動及碰撞是一個非常基礎但是很普遍物理過程，在許多場景中都會發生，小到微觀世界中研究氣體反應的分子碰撞理論，大到宇宙中星體運動及星系演化。不管是游戲領域，機器人領域，還是計算化學領域，都會涉及到物體運動和碰撞過程。

關于小球運動及碰撞的仿真及代碼，網站資源非常多，但絕大多數以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仿真模型的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。