Simulation of two spheres colliding through classical mechanics
<p>I am trying to simulate the collision of two spheres on Fortran 90 by using the equations of classical mechanics. </p>

<p>Firstly, I'm doing this on one dimension: I am fixing a sphere on the ground and putting the other one at a certain height (both are at rest at t=0), and both of them have the same radius and mass. </p>

<p>Being so, the forces acting at the moving particle is the gravitational force and the elastic force (when the spheres start to collide). So, at first, the moving sphere is only subjected to gravitational force, that is, Fg = -g.m. I'm considering only one dimension, so since the particle falling will give me lower values than the height it was in the beginning, it is <strong>negative</strong>. Now, when the moving particle hits the fixed one, the latter will cause a force of Fe = -k.dx.n, where k is the elastic constant, dx is the overlap and n is the direction of that force (its norm is 1). The elastic force is <strong>restitutive</strong>, so it should be <strong>negative</strong> also. To be clear, dx = r1 + r2 - |x1 - x2| and n = (x1-x2)/|x1 -x2|, where r1, r2 are the radii and x1, x2 are the positions of the particles at given time.</p>

<p><strong><em>The code compiles ok</em></strong>, the only problem is when I go check the positions for the moving particle - sometimes, <strong>it keeps getting lower</strong> (as it should at first) and then it <strong><em>continues</em></strong> to decrease to negative values - something that should not happen. Or, when the moving sphere <strong><em>does</em></strong> bounce back where it should, it goes back to a height <strong>HIGHER</strong> than the initial condition (<strong>it bounces higher, violating conservation of energy</strong>).</p>

<p>The thing is that, the <strong>best situation</strong> that it seems to be working (that last one, in which it bounces higher), the formula for the elastic force is actually not right (I put dx = r1 + r2 - |x1 <strong>+</strong> x2|, instead of a minus sign). I've tried correcting the higher bounce by changing some parameters (lowering and increasing elastic constant and lowering the time step), but with no success. Also, the higher bouncing I mentioned is about four times higher, so it is not some problem of error of calculations. </p>

<p>Here are the important parts of the code: where the forces are calculated, and the time integration:</p>

<pre><code>coef1 = abs(1.0d0*x1 - 1.0d0*x2)
coef2 = (1.0d0*x1 - 1.0d0*x2)/coef1
Fx_elasticaparcial(m) = -1.0d0*k*((1.0d0*a1+1.0d0*a2) + 1.0d0*coef1)*coef2
Fx_elasticatotal(l) = 1.0d0*Fx_elasticatotal(l) +1.0d0*Fx_elasticaparcial(m)

forcax(l) = 1.0d0*flagx(l)*(1.0d0*gravidade(massadummy)+1.0d0*Fx_elasticatotal(l))
</code></pre>

<p><strong>m</strong> is the partial elastic force and <strong>l</strong> is the total elastic force exerted ON a particle (I've done it thinking of adding more particles further down the road)</p>

<pre><code>function gravidade(massa)
use parametros, only: g
implicit none
real*8 :: gravidade, massa
gravidade = -1.0d0*g*massa
end function gravidade
</code></pre>

<p><strong><em>This one is a function for gravitational force</em></strong></p>

<pre><code>subroutine integracaoEuler (xo, xn, vxo, vxn, fx, mass)
use parametros, only: N, dt
implicit none
real*8 :: xo, xn, vxo, vxn, fx, mass

xn = 1.0d0*xo + 1.0d0*dt*vxo
xo = 1.0d0*xn
vxn = 1.0d0*vxo + 1.0d0*dt*((1.0d0*fx)/(1.0d0*mass))
vxo = 1.0d0*vxn

return

end subroutine integracaoEuler
</code></pre>

<p><strong><em>Euler integration subroutine</em></strong></p>

<pre><code>do while ((t + dt) &lt; (tmax)) !time loop
    t = t+dt
    do i = 1, N !particle loop
        do j = 1, N !neighbours of i particle loop
            if (i .NE. j) then
                if (abs(SQRT(1.0d0*xold(i) - 1.0d0*xold(j))**2) .LE. 1.0d0*a(i) + 1.0d0*a(j)) then
                    call forcastotais (i, j, m(i), a(i), a(j), xold(i), xold(j))
                else
                    forcax(i) = 1.0d0*flagx(i)*gravidade(m(i))
                end if
            end if
        end do
        call integracaoEuler (xold(i), xnew(i), vxold(i), vxnew(i), forcax(i), m(i))
        write(i, *) t, xold(i)
    end do
end do
</code></pre>

<p><strong><em>Main loop of the program</em></strong></p>

<p>Sorry for making this really big, to fix it probably must be something easy. Again, I need to really use the equations of <strong>classical mechanics</strong> for this, not by putting the velocity to change direction when it hits the other sphere or that sort of stuff. I am frustrated with this, and don't know how to fix it (I have been stuck on this for 2 weeks). I've done everything else on the model I'm making (viscosity, cell-link) with a lot of particles, but this problem keeps reappearing and messing up my graphs. Any help is welcome; please make it explicit, since I am new to programming. Thank you in advance.</p>