Simulation of Elastic Collisions

[MatinSAR]

Homework Statement

Create a simulation to model elastic collisions among particles in a 2D square area. Ensure the particles reflect elastically off the boundaries.

Relevant Equations

$\vec P_i = \vec P_f$
$T_i = T_f$

Hello everyone,

Can we use the conversion of kinetic energy (T) and momentum ($\vec P$) to find the velocities after an elastic collision in a 2D problem? Can we treat objects like particles?

 

[kuruman]

Sure. Start with the equations giving the final velocities of the two particles in terms of the initial velocities in one dimension. (See http://hyperphysics.phy-astr.gsu.edu/hbase/elacol2.html). Note that this "one" dimension is along a line joining the centers of the masses and the ignored, "second", dimension is the direction perpendicular to that. In a two dimensional simulation and in the absence of friction, the component of the velocity of each mass in the perpendicular direction does not change. What changes is the parallel component of the velocity of each mass according to the equations provided in the reference.

Can you find the parallel component for an arbitrary direction of the velocity of each mass? Hint: Find a unit vector in the direction from mass 1 to mass 2.

 

[MatinSAR]

Can you find the parallel component for an arbitrary direction of the velocity of each mass? Hint: Find a unit vector in the direction from mass 1 to mass 2.

Given that we know the positions of each particle, we can use $\dfrac{\vec{r}_2 - \vec{r}_1}{|\vec{r}_2 - \vec{r}_1|}$ as the unit vector pointing from particle 1 to particle 2.

I believe I understand your concept. However, there's a problem. In my code, I can't predict which two particles will collide to use the line connecting them as our one-dimensional reference, where velocity changes only in that specific dimension.
Additionally, I don't need an analytical solution. As long as I have the correct equations that provide the magnitude and direction of each particle's speed after a collision, I can solve the system of equations numerically.

 

[MatinSAR]

This was my intended approach.

$$m_1u_1=m_1v_1\cos\theta_1 + m_2v_2\cos\theta_2$$$$0=m_1v_1\sin\theta_1 - m_2v_2\sin\theta_2$$$$m_1u^2_1=m_1v^2_1+m_2v^2_2$$
There are two issues here. First, we have four unknowns and only three equations. Second, is this really a 2D problem? I mean, in 2D, can we consider particles to only have head-on collisions even if they have orthogonal velocities?

 

[erobz]

This was my intended approach.

View attachment 355282​

$$m_1u_1=m_1v_1\cos\theta_1 + m_2v_2\cos\theta_2$$$$0=m_1v_1\sin\theta_1 - m_2v_2\sin\theta_2$$$$m_1u^2_1=m_1v^2_1+m_2v^2_2$$
There are two issues here. First, we have four unknowns and only three equations. Second, is this really a 2D problem? I mean, in 2D, can we consider particles to only have head-on collisions even if they have orthogonal velocities?

How are they flying off in 2 directions if they are particles without dimension. I could see this scenario as something that could happen with things colliding that have some physical geometry (where a glancing blow is possible), but point masses...I don't get it.

To me, there are three scenarios $m_1$ stops - $m_2$ moves directly to the right, $m_1$ moves to the left - $m_2$ to the right, or they both continue to the right.

 

[MatinSAR]

How are they flying off in 2 directions if they are particles without dimension. I could see this scenario as something that could happen with things colliding that have some physical geometry (where a glancing blow is possible), but point masses...I don't get it.

To me, there are three scenarios $m_1$ stops, $m_2$ moves directly to the right, $m_1$ moves to the left, $m_2$ to the left, and they both continue to the right.

Third problem. Perhaps I should visit my professor's office tomorrow to discuss these issues with him.

 

[Steve4Physics]

Homework Statement: Create a simulation to model elastic collisions among particles in a 2D square area. Ensure the particles reflect elastically off the boundaries.

The problem statement seems very vague Some questions that spring to mind are:
- are the particles to be modelled as points or as discs?
- will this be a real-time simulation?
- how many (roughly) particles are in the square?
- what are the rough relative sizes of the particles, box and average particle-separation?
Or are you allowed to decide these sorts of parameters for yourself?

Knowing the key parameters could help you to decide what simplifications/approximations (if any) you might want to make. This might have an impact on your approach.

The level of sophistication of the simulation will be limited by how many person-hours of effort are allocated for the work – and the level of your own programming skills! Is there a ball-park figure?

A good start point (IMO) would be to establish answers to the above (and any other relevant) questions before diving into the detail. Then formulate an overall strategy for how to tackle the job.

So maybe there are some additional things to discuss with your professor!

 

[erobz]

Third problem. Perhaps I should visit my professor's office tomorrow to discuss these issues with him.

I could be missing something; someone will straighten this out. Spend some time thinking about it today, and then go to your prof. I was just throwing in my two cents about specifically that scenario you proposed. An oblique collision could send particles (without physical geometry) off in different directions.

 

[kuruman]

How are they flying off in 2 directions if they are particles without dimension. I could see this scenario as something that could happen with things colliding that have some physical geometry (where a glancing blow is possible), but point masses...I don't get it.

To me, there are three scenarios $m_1$ stops - $m_2$ moves directly to the right, $m_1$ moves to the left - $m_2$ to the right, or they both continue to the right.

There is no mention in the original post about point masses. As you say, there can be no 2d collision between point masses. I would assume discs of radii $R_1$ and $R_2$. If $\mathbf r_1$ and $\mathbf r_2$ are their positions, the criterion for a collision to occur is $|\mathbf r_2-\mathbf r_1|\leq(R_1+R_2).$

To @MatinSAR:
I am in the process of answering your posts #3 and #4. Please be patient.

 

[MatinSAR]

Or are you allowed to decide these sorts of parameters for yourself?

Initially, all the decisions are mine. Once I complete my work based on my preferences, I will discuss it with him, and he will add more details. I recall he mentioned that he might ask me to apply it to random shapes like squares or pentagons, and I would need to consider rotation and rotational inertia as well. I really hope that won't be the case since it's complicated enough already.

 

[Steve4Physics]

... I recall he mentioned that he might ask me to apply it to random shapes like squares or pentagons, and I would need to consider rotation and rotational inertia as well.

That seems to imply that you are expected (at least initially) to model the 'particles' as discs rather than as points.

Squares or pentagons? - aagh!

really hope that won't be the case since it's complicated enough already.

Yes indeed!

 

[kuruman]

Initially, all the decisions are mine. Once I complete my work based on my preferences, I will discuss it with him, and he will add more details. I recall he mentioned that he might ask me to apply it to random shapes like squares or pentagons, and I would need to consider rotation and rotational inertia as well. I really hope that won't be the case since it's complicated enough already.

I think you should first simulate two identical colliding discs and no spin to see what the issues are. Then consider two discs of unequal radii, then consider many disks all of unequal radii. Then you will see that squares and polygons can become a nightmare.

The main point of view to adopt is that you don't need an analytical expression to describe the mass trajectories. Since you are doing a simulation, you can use the computer to calculate the next position in equal intervals $\Delta t$ using the straight line trajectory equation $$\mathbf r_{\text{new,before}} =\mathbf r_{\text{old,before}}+\mathbf v_{\text{before}}~\Delta t.\tag{1}$$ If the collision criterion is satisfied within such an interval, you assume that the position remains the same the over that interval but the velocities change to their after the collision values then keep calculating straight line motion trajectories using $$\mathbf r_{\text{new,after}} =\mathbf r_{\text{old,before}}+\mathbf v_{\text{after}}~\Delta t.\tag{2}$$

 

[MatinSAR]

Squares or pentagons? - aagh!

At first, I was about to accept it as my project, but I'm really glad that didn't happen.

Then you will see that squares and polygons can become a nightmare.

It's already clear enough at this point.

The main point of view to adopt is that you don't need an analytical expression to describe the mass trajectories. Since you are doing a simulation, you can use the computer to calculate the next position in equal intervals $\Delta t$ using the straight line trajectory equation $$\mathbf r_{\text{new,before}} =\mathbf r_{\text{old,before}}+\mathbf v_{\text{before}}~\Delta t.\tag{1}$$ If the collision criterion is satisfied within such an interval, you assume that the position remains the same the over that interval but the velocities change to their after the collision values then keep calculating straight line motion trajectories using $$\mathbf r_{\text{new,after}} =\mathbf r_{\text{old,before}}+\mathbf v_{\text{after}}~\Delta t.\tag{2}$$

That's exactly what I did for my first simulation with just two disks. It's an effective method, but it doesn't work for disks with high speed because the position intervals become too large, causing the code to miss collisions. Reducing the time step (dt) to a very small value introduces its own problems.

Thank you, @kuruman, @Steve4Physics, and @erobz, for your time. I'll discuss this with him tomorrow. I hope this issue will be resolved with my professor's assistance tomorrow. However, if any problems persist regarding the physics of the problem, I will reach out again.

 

[haruspex]

Depending on what you want to do with the result of the simulation, there may be an alternative to taking tiny steps of time. For each pair of particles and each particle-boundary pair you can calculate if and when they will next collide, then set your $\Delta t$ as the minimum. This can greatly speed up the simulation if you are not, e.g., creating a video. The result is also more accurate usually.

 

[MatinSAR]

if you are not, e.g., creating a video.

Thank you for your reply. I am creating a gif.

 

[kuruman]

You can do what you want because it's your simulation. If I were simulating this, I would use the following piece of flowchart that requires only the 1-d equations linking the final velocity components to the initial velocity components.
Start with initial positions $\mathbf r_{1i},~\mathbf r_{2i}$ and velocities $\mathbf v_{1i},~\mathbf v_{2i}.$
Calculate the tentative final positions using $\mathbf r_{1f}= \mathbf r_{1i}+\mathbf v_{1i}\Delta t$ and $\mathbf r_{2f}= \mathbf r_{1i}+\mathbf v_{2i}\Delta t$.
If the tentative final positions do not imply a collision, then they become the actual final positions in which case the become the initial positions and you go back to step 1.
If the tentative final positions imply a collision, then the initial positions are kept, and you go back to step 1 with $\mathbf v_{1f},~\mathbf v_{2f}$ becoming $\mathbf v_{1i},~\mathbf v_{2i}.$ I show below how to get these.

A unit vector parallel to the center-to-center direction from mass 1 towards mass 2 is, in terms of the Cartesian axes that you have adopted, $$ \mathbf{\hat n}_{\parallel} =\frac{\mathbf r_{2}-\mathbf r_{1}}{|\mathbf r_{2}-\mathbf r_{1}|} =\frac{ (x_2-x_1)~\mathbf {\hat x}+ (y_2-y_1)~\mathbf {\hat y} } {\left[(x_2-x_1)^2+ (y_2-y_1)^2 \right]^{1/2}}$$ Once you have the parallel unit vector, the perpendicular unit vector can be written as $$\mathbf{\hat n}_{\perp}=\frac{ (y_2-y_1)~\mathbf {\hat x}- (x_2-x_1)~\mathbf {\hat y} } {\left[(x_2-x_1)^2+ (y_2-y_1)^2 \right]^{1/2}}.$$ The velocities before the collision have components in the parallel and perpendicular directions given by, $$ v_{k,\parallel}^{\text{before}} = \mathbf{v}_k^{\text{before}}\cdot\mathbf{\hat n}_{\parallel}~; ~~
v_{k,\perp}^{\text{before}} = \mathbf{v}_k^{\text{before}}\cdot\mathbf{\hat n}_{\perp}~~~~(k=1,2).
$$After the collision, the perpendicular components are the same as before. The parallel components are related to the "before" components as given in the reference in post #2.

What we have done here is calculate the final velocities in a rotated frame where one has a 1d collision along the center-to-center distance which is the direction along which the momentum transfer takes place.

 

[MatinSAR]

You can do what you want because it's your simulation. If I were simulating this, I would use the following piece of flowchart that requires only the 1-d equations linking the final velocity components to the initial velocity components.
Start with initial positions $\mathbf r_{1i},~\mathbf r_{2i}$ and velocities $\mathbf v_{1i},~\mathbf v_{2i}.$
Calculate the tentative final positions using $\mathbf r_{1f}= \mathbf r_{1i}+\mathbf v_{1i}\Delta t$ and $\mathbf r_{2f}= \mathbf r_{1i}+\mathbf v_{2i}\Delta t$.
If the tentative final positions do not imply a collision, then they become the actual final positions in which case the become the initial positions and you go back to step 1.
If the tentative final positions imply a collision, then the initial positions are kept, and you go back to step 1 with $\mathbf v_{1f},~\mathbf v_{2f}$ becoming $\mathbf v_{1i},~\mathbf v_{2i}.$ I show below how to get these.

A unit vector parallel to the center-to-center direction from mass 1 towards mass 2 is, in terms of the Cartesian axes that you have adopted, $$ \mathbf{\hat n}_{\parallel} =\frac{\mathbf r_{2}-\mathbf r_{1}}{|\mathbf r_{2}-\mathbf r_{1}|} =\frac{ (x_2-x_1)~\mathbf {\hat x}+ (y_2-y_1)~\mathbf {\hat y} } {\left[(x_2-x_1)^2+ (y_2-y_1)^2 \right]^{1/2}}$$ Once you have the parallel unit vector, the perpendicular unit vector can be written as $$\mathbf{\hat n}_{\perp}=\frac{ (y_2-y_1)~\mathbf {\hat x}- (x_2-x_1)~\mathbf {\hat y} } {\left[(x_2-x_1)^2+ (y_2-y_1)^2 \right]^{1/2}}.$$ The velocities before the collision have components in the parallel and perpendicular directions given by, $$ v_{k,\parallel}^{\text{before}} = \mathbf{v}_k^{\text{before}}\cdot\mathbf{\hat n}_{\parallel}~; ~~
v_{k,\perp}^{\text{before}} = \mathbf{v}_k^{\text{before}}\cdot\mathbf{\hat n}_{\perp}~~~~(k=1,2).
$$After the collision, the perpendicular components are the same as before. The parallel components are related to the "before" components as given in the reference in post #2.

What we have done here is calculate the final velocities in a rotated frame where one has a 1d collision along the center-to-center distance which is the direction along which the momentum transfer takes place.

Impressive insights! Thank you, @kuruman. I will proceed with this method and see where I get with it.

 

[haruspex]

If the tentative final positions

What if the overlap of the discs is only for part of the way between the two positions? Are you assuming $\Delta t$ is so small we can ignore such glancing blows?

then they become the actual final positions

… except that there are all the other possible collisions these two particles might have before that. Having thought about it some more, I think the approach I described in post #14 avoids some pitfalls. Given that a video is to be generated, I would:

• choose an interframe time, $t_{frame}$
• for each frame {
  • $starttime=clock()$
  • initialise $t_{next}=t_{frame}$
  • while $t_{next}>0$ {
    • $t_{collide}=t_{next}$
    • go through all possible collisions to see if any happen before $t_{collide}$, reducing $t_{collide}$ accordingly
    • compute all positions and velocities after $t_{collide}$ elapses
    • subtract $t_{collide}$ from $t_{next}$
  • }
  • generate frame
  • $now=clock()$
  • if $(starttime+t_{frame}>now)$ sleep($starttime+t_{frame}-now$)
  • // watch out for timestamp wrapping
• }

 

[MatinSAR]

Thank you for your reply @haruspex ...

What if the overlap of the discs is only for part of the way between the two positions? Are you assuming Δt is so small we can ignore such glancing blows?

This problem existed in my initial code too, as mentioned in post #13. If I start with velocities smaller than 100, then vΔt for Δt=0.01 would be smaller than 1. So, I use disks with a radius greater than 1. I understand this is a limitation, but since this is only a university project, I will ask my professor if this velocity limit is a negative point or not.

except that there are all the other possible collisions these two particles might have before that.

I do not understand this part. Could you please explain it more?

Having thought about it some more, I think the approach I described in post #14 avoids some pitfalls

But isn't it more complicated to predict which two particles are going to collide?

 

[haruspex]

I do not understand this part. Could you please explain it more?

But isn't it more complicated to predict which two particles are going to collide?

You have to take the particles in pairs to check if they have collided. According to post #16, if the first pair you pick, A and B, have not collided in time $\Delta t$ then you update their positions on the basis of non collision. But then you find A would have collided with C, so the update is wrong. Likewise, if you find A and B have collided and update them you could then find one of those would have collided with C sooner.
The only way to avoid this is to calculate for every pair (including particle and wall pairs) if and when they would next collide and recompute the whole ensemble at the earliest collision time you find.
The rest of the logic in post #18 is to arrange that the frames are equally spaced in elapsed time.

 

[MatinSAR]

You have to take the particles in pairs to check if they have collided. According to post #16, if the first pair you pick, A and B, have not collided in time Δt then you update their positions on the basis of non collision. But then you find A would have collided with C, so the update is wrong.

Why check for just a pair? We can update the positions of all particles and then check them all in pairs. The issue arises when three or more particles collide. In such cases, I'll print an error message and break the code. Am I missing something?

The only way to avoid this is to calculate for every pair (including particle and wall pairs) if and when they would next collide and recompute the whole ensemble at the earliest collision time you find.

Isn't it difficult to predict when they will collide? Each particle has the chance to collide with every other particle.

 

[haruspex]

Why check for just a pair? We can update the positions of all particles and then check them all in pairs.

If I understand your algorithm, you update all positions as though no collisions occur then check for overlapping objects. Immediately there is the possibility that you overlook a collision because the overlap occurred only briefly during the interval.
Next, how do you update the positions? Do you leave them overlapping or backtrack to where they would have made contact?

The issue arises when three or more particles collide.

How do you detect that?

In such cases, I'll print an error message and break the code. Am I missing something?

So the run is aborted? You could find that happening most of the time, and the runs that succeed will not be a representative sample. Besides, that's poor design when you can so easily make it always work.

Isn't it difficult to predict when they will collide?

It is just a matter of solving a quadratic, no? If the quadratic has no solution then they do not collide.

Each particle has the chance to collide with every other particle.

Sure, but you can calculate if and when any given pair of objects (particles and walls) would collide in the absence of all other objects then determine which happens first.

 

[MatinSAR]

If I understand your algorithm, you update all positions as though no collisions occur then check for overlapping objects. Immediately there is the possibility that you overlook a collision because the overlap occurred only briefly during the interval.

Yes this problem with my algorithm was mentioned before:

This problem existed in my initial code too, as mentioned in post #13. If I start with velocities smaller than 100, then vΔt for Δt=0.01 would be smaller than 1. So, I use disks with a radius greater than 1. I understand this is a limitation, but since this is only a university project, I will ask my professor if this velocity limit is a negative point or not.

Next, how do you update the positions? Do you leave them overlapping or backtrack to where they would have made contact?

For particles that collide, I calculate their velocities after the collision using their past positions and take another step using these velocities. For particles that don’t collide, I use their current positions and velocities and take another step.

How do you detect that?

The code check if any two pairs collide, It might detect more than one, and we cannot choose one pair.

So the run is aborted? You could find that happening most of the time, and the runs that succeed will not be a representative sample. Besides, that's poor design when you can so easily make it always work.

Is that really easy?
What formulas you suggust using to calculate velocities after each collision?

It is just a matter of solving a quadratic, no? If the quadratic has no solution then they do not collide.

By difficult, I meant we have many particles... But I don't understand why quadratic? The velocities are constant.

Sure, but you can calculate if and when any given pair of objects (particles and walls) would collide in the absence of all other objects then determine which happens first.

Great!!! I did not think of it this way. Thank you! I will think on this more.

 

[MatinSAR]

@haruspex I think I understand your point now, or at least I understand a better way to do it.
As you suggusted we predict which particles (or a particle and wall) are going to collide and when. Then we use formulas mentioned by @kuruman for velocities after collisions. Right? Because we know which particles are going to collide we can use their initial positions to use the coordinate in which only in one direction velocity changes. Then after each collision using updated position and velocity of those two particles and we calcluate which particels are going to collide first again.

 

[haruspex]

@haruspex I think I understand your point now, or at least I understand a better way to do it.
As you suggusted we predict which particles (or a particle and wall) are going to collide and when. Then we use formulas mentioned by @kuruman for velocities after collisions. Right? Because we know which particles are going to collide we can use their initial positions to use the coordinate in which only in one direction velocity changes. Then after each collision using updated position and velocity of those two particles and we calcluate which particels are going to collide first again.

Yes.

By difficult, I meant we have many particles...

Whether it is three or 1000, the logic (software) is the same, just takes longer to execute. Ultimately, there would be a limit on how fast the particles could move in the video.

But I don't understand why quadratic? The velocities are constant.

You can write the x(t) and y(t) functions for each particle in the form $x_i(t)=x_{i0}+v_{xi0}t$. The equation $(x_1(t)-x_2(t))^2+(y_1(t)-y_2(t))^2=4R^2$ gives a quadratic in t.

 

[MatinSAR]

Yes.

Thank you for your time and valuable help @haruspex ...
So I'm going to try this algorithm and see where I get with it..

Whether it is three or 1000, the logic (software) is the same, just takes longer to execute. Ultimately, there would be a limit on how fast the particles could move in the video.

You can write the x(t) and y(t) functions for each particle in the form $x_i(t)=x_{i0}+v_{xi0}t$. The equation $(x_1(t)-x_2(t))^2+(y_1(t)-y_2(t))^2=4R^2$ gives a quadratic in t.

I get it.

 

[MatinSAR]

Hello everyone,

I have one more question. For the equations in post #16, do I need to consider one particle stationary and compute the other particle's velocity with respect to it, then use those equations? Because only in this frame, one particle is stationary... When I calculate velocities after the collision using the conservation law, should I then transfer the equations back to the lab frame?

Please let me know if this is correct.

 

[MatinSAR]

Another challenge I encountered is that I considered velocity changes only along the line connecting the moving and stationary particles during a collision. However, in many references, I see that velocity changes occur in two directions.

https://phys.libretexts.org/Bookshe...ory/15.06:_Two_Dimensional_Elastic_Collisions

 

[haruspex]

Another challenge I encountered is that I considered velocity changes only along the line connecting the moving and stationary particles during a collision. However, in many references, I see that velocity changes occur in two directions.
View attachment 355372
https://phys.libretexts.org/Bookshelves/Classical_Mechanics/Classical_Mechanics_(Dourmashkin)/15:_Collision_Theory/15.06:_Two_Dimensional_Elastic_Collisions

Yes, you need to treat the collisions as oblique.
In using the equations at that link, there are ambiguities you need to resolve: square root signs and arctan values. Do you understand how to do that?

 

[MatinSAR]

Yes, you need to treat the collisions as oblique.
In using the equations at that link, there are ambiguities you need to resolve: square root signs and arctan values. Do you understand how to do that?

In my code, I assumed all collisions to be head-on, meaning the only direction in which the velocity could change was along the line connecting the colliding particles. This simplification leads to equations similar to those found in this link:
15.4: One-Dimensional Collisions Between Two Objects

Is this approach to simulation incorrect?

 

[MatinSAR]

Do you understand how to do that?

I'm having trouble understanding how to use those equations to determine the four unknowns (two angles and two speeds). IMO we need to measure at least one angle in the lab to be able to calculate the magnitudes of the velocities after the collision.

For example, the equation above don't provide any information about the velocities post-collision because we don't know the angles $\theta_{1,f}$ and $\theta_{2,f}$ and we do not know $v_{1,f}$.

Edit: By reviewing the equations in post #4, I am now 90% sure that I cannot use the equations from that link-15.6: Two Dimensional Elastic Collisions- in my code...

We need to measure one of the angles in the lab frame. Don't you agree, @haruspex?

 

[kuruman]

Isn't it difficult to predict when they will collide? Each particle has the chance to collide with every other particle.

In such problems, one disc is modeled as a point disc and the other a disc of radius $2R$. The momentum and energy transfer are the same. See figure below.

You don't need to know well in advance when two particles will collide. You need to know whether they will collide in the next interval $\Delta t$. To figure that out, first you calculate where they will be if their velocities were the same between $t=t_n$ and $t=t_n+\Delta t$. If it turns out that the point disc is inside the other, $(|\mathbf r_2-\mathbf r_1|<2R)$, then you conclude that the collision (assumed instantaneous) will be completed within the interval $\Delta t$ therefore you must discontinuously change the velocities before the collision to the velocities after the collision as described in post #16. You know where each particle is because you have to tell the machine to put it there after each$\Delta r$. Before you do that, you test if the new positions imply that one disc is inside the other. If it turns out that this will the case, you keep the positions the same and change the velocities before the collision happens. The figure below shows what I mean. The "Before" figure on the left and the "Accepted" figure on the right can be considered as two consecutive frames in the simulation.

In my code, I assumed all collisions to be head-on, meaning the only direction in which the velocity could change was along the line connecting the colliding particles. This simplification leads to equations similar to those found in this link:
15.4: One-Dimensional Collisions Between Two Objects

Is this approach to simulation incorrect?

You can use the equations in that reference (they are the same as the ones in the hyperphysics link I gave you) but you have to understand that subscript $x$ in them is the same as subscript $\parallel$ in post #16. The scattering angles are not needed if you choose to adopt my method for the simulation, which you don't have to do. I am just pointing out how approached such simulations in the past.

 

[MatinSAR]

You can use the equations in that reference (they are the same as the ones in the hyperphysics link I gave you) but you have to understand that subscript $x$ in them is the same as subscript $\parallel$ in post #16. The scattering angles are not needed if you choose to adopt my method for the simulation, which you don't have to do. I am just pointing out how approached such simulations in the past.

Thanks for the clarification with the algorithm. The current problem is related to the physics of the problem. I'm not sure if considering head-on collisions is a good simplification or not, as it makes the final velocity equations after collision too simple.


The most important lesson for me here is:
The cost of choosing a project topic later than other students is that you'll end up with the most difficult one. I'm really glad I didn't accept the professor's offer to enter the rotation in simulation.

 

[MatinSAR]

I appreciate your help, but I think considering head-on collisions is enough for my project. Thank you, @haruspex and @kuruman, for your time and assistance.

 

[haruspex]

In my code, I assumed all collisions to be head-on, meaning the only direction in which the velocity could change was along the line connecting the colliding particles.

The velocities only change along the direction of the line joining the mass centres, whether head-on or oblique. @kuruman explained that in post #16. The process is:

• calculate the direction of that line
• resolve the velocities into coordinates parallel to and normal to that line
• convert velocities along the line into being relative to one of the discs
• compute the new velocities along that line
• convert back to the lab frame
• take the velocities normal to the line of centres to be unchanged
• resolve the velocities back into your standard xy coordinates

You also need to watch out for overflow when doing divisions or using the tan function. E.g. when coding the use of a division you will need to think about the physical meaning of the denominator being extremely small. Is it possible? Is there a way to cope with it?
E.g. if the library arctan function is of the form
double arctan(double ratio)
then calling it as arctan(y/x) might blow up in the division.
A solution would be to front-end it with a dyadic version, solving the arctan ambiguity at the same time. The library function would always return a value in the range $(-\pi/2, +\pi/2)$. I think this is right:

double arctan(double x, double y) {
if (abs(x) < MinInvertible) {
if (y > 0) return Pi/2;
return -Pi/2;
}
if (y > 0) return arctan(y/x);
return arctan(y/x) + Pi;
}

But maybe the library already has a dyadic version.

An interesting sanity check you can do is to keep track of the total KE. Because of rounding errors it won't be constant, but it should only drift slowly. If it drifts persistently up or down then you could artificially tweak all the velocities to keep it in bounds.
(This is trickier in simulations involving other energy modes, such as PE, or other conservations, such as momentum. What should be tweaked?)