How do I apply a gravitational gradient to this code? (particles in elastic collision)

1 view (last 30 days)
Hello,
I am attempting to simulate particles on an asteroid surface and whether or not they will "flyaway" from the asteroid surface given an input force.
The input force would be applied sideways horizontally.
I have managed to simulate an asteroid surface and two particles which will collide with eachother, however I'm now looking to add a gravitation field gradient to see if the initial velocities given would be overcome by the asteroid gravity (and hence the particle would return to the asteroid surface), or if the initial velocity would overcome the gravity (and hence the particle would flyaway from the asteroid surface indefinitely).
I'm currently considering doing this via multiple if statements dependant on the current position of the particles, and then subtracting from the velocity of the particle using a value that corresponds to the gravitational effect on the velocity. However I'm not sure how to quantify this effect, or where to fit this into my current code.
My current code is below. FYI, requires curve function from PCB toolbox to run.
%
clearvars
close all
%-----------------------
N=2;% number of particles
L=2;%length of the box - no
asteroid=curve(Radius=20,ReferencePoint=[0,0],ArcAngle=[70,(110)]);
dt=1;
X=zeros(2,N); %original position matrix
V=zeros(2,N); %original velocity matrix
R=0.3*ones(1,N); %radius of each particle
M=ones(1,N)+3*rand(1,N)+2*R; % mass of each particle
theta=0:2*pi/20:2*pi;%
nstep=30; % number of time steps
pit=0.1; % number of iterations between two plots
%---------- manually initilizing position and velocity
%position
X=[-3 3;23 23]; %-1.899 to match box of 2 with rad 0.3*rand(1,N)+0.1;
%velocity
V=[0.1 -0.1;-0.4 -0.4];
%Radius of each particule
R=[0.3 0.3];
%Mass of each particle
M=[2 2];
%auxiliary variables
Ax=zeros(N,N);
Ay=zeros(N,N);
Ar=zeros(N,N);
for i=1:N
for j=i:N
Ar(i,j)=R(1,i)+R(1,j);
end
end
hold on
axis equal
for i=1:N %plot particles initial
x=X(1,i)+R(1,i)*cos(theta);
y=X(2,i)+R(1,i)*sin(theta);
plot(x,y,'color','b')
end
plot(asteroid)
%plotdomain(L)
drawnow
pause(1)
for T=0:nstep
X=X+dt*V; %
%check if particules collided with eachother
for i=1:N
Ax(i,:)=X(1,i)-X(1,:);
Ay(i,:)=X(2,i)-X(2,:);
end
Ax=triu(Ax);
Ay=triu(Ay);
Nrm=(Ax.^2+Ay.^2).^(0.5)-Ar-10^-3; % calculate the distance between each two particles
Nrm=triu(Nrm(1:N-1,2:N));
[row,col]=find(Nrm<0); % find the particles that collided
l=length(row);
%if collided calculate the new velocities
if l~=0
col=col+1;
for t=1:l
i=row(t,1);
j=col(t,1);
C1=(2*M(1,j)/(M(1,i)+M(1,j)))*dot(V(:,i)-V(:,j),X(:,i)-X(:,j))/(norm(X(:,i)-X(:,j))^2);
C2=(2*M(1,i)/(M(1,i)+M(1,j)))*dot(V(:,j)-V(:,i),X(:,j)-X(:,i))/(norm(X(:,i)-X(:,j))^2);
V(:,i)=V(:,i)-C1*(X(:,i)-X(:,j));
V(:,j)=V(:,j)-C2*(X(:,j)-X(:,i));
end
end
%check if particles collide with asteroid
for i=1:N
if asteroid.ReferencePoint(1,2)+X(2,i)<= R(1,i)+asteroid.Radius
V(2,i)=-V(2,i);
end
end
%plot
if T==pit*ceil(T/pit)
clf
hold on
axis equal
for i=1:N %replot each particles beyond step 1!
x=X(1,i)+R(1,i)*cos(theta);
y=X(2,i)+R(1,i)*sin(theta);
plot(x,y,'color','b')
end
% plotdomain(L) %plots the boundaries visuals
plot(asteroid) %plots the asteroid beyond step 1!
title(T*dt)
drawnow
pause(0.01)
end
end
Thank you in advance for the help.
Kind regards
Alex

Answers (1)

David Goodmanson
David Goodmanson on 26 Dec 2024
Edited: David Goodmanson on 26 Dec 2024
Hi Alex
If you don't need the details of the actual trajectory, the method is pretty simple. Modeling the asteroid as a sphere of radius Ra and mass M in time-honored physics tradition, just calculate the total energy of the particle after the collision
E = (m/2)(v_horz^2+v_vert^2) - GmM/Ra
where the vertical velocity has to be away from the surface which you have already corrected for. The particle is unbounded and flies away, or is bound to the asteroid and will eventually return, according to E>0 or E<0.

Categories

Find more on Environment in Help Center and File Exchange

Products


Release

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!