I have some experience with C++ programming, but I am quite new to Mathematica and I think the latter requires a different mindset which I am still not used to.

I would like to simulate the effect caused by a sphere in free-fall when it hits a liquid, which I consider to be incompressible and inviscid. I will also neglect air resistance. The problem can be reduced to a 2-dimensional problem in cylindrical coordinates (r, z), since the system is symmetric under an angle translation. The sphere is then just a circle. The final code should show the time evolution of the ripples produced by the spere in the liquid.
From what I read in the many available examples, one way to define the region over which I am going to solve the differential equations is the following (let me say that the liquid is initially at rest at z=4):

ballRadiusSquared = 0.5;
zmax = 7.;  
zmin = 0.;  
rmax = 10.;
g = 9.8;  
zBallCentre[t_] := -0.5*g*t^2 + zmax;  
region[t_] = 
    r^2 + (z - zBallCentre[t])^2 >= ballRadiusSquared, 
    {{r, 0, rmax}, {z, zmin, zmax}}];

Region of interest for t=1.

Note that I did not include the upward forced caused by the liquid pressure, that would slow the circle down as soon as it reaches the z=4 height.

I now have to solve the Laplace equation in cylindrical coordinates over this region:

D[r*D[ϕ[t, r, z], r], r]/r + D[ϕ[t, r, z], {z, 2}] == 0

From the moment the potential (ϕ[t,r,z]) is known I can immediately calculate its gradient to obtain the velocity at all points inside the region of interest. Finally, to obtain the pressure I can use the Navier-Stokes equation. The important aspect of this procedure is that I will have to calculate the derivative of the potential with respect to time.

It is clear that I have to give the right boundary conditions to NDSolve.

The Neumann boundary conditions for r=1, r=rmax and z=zmin and the Dirichlet condition for r=rmax are easy to implement:

DirCond[t_] = DirichletCondition[ϕ[t, r, z] == 0., r == rmax];  
NeumCond[t_] = 
  NeumannValue[0., r == 0] + NeumannValue[0., r == rmax] + 
    NeumannValue[0., z == zmin];

To solve the Laplace equation I imagine I would have to do something like:

DiffEq = NDSolveValue[{D[r*D[ϕ[t, r, z], r], r]/r + 
         D[ϕ[t, r, z], {z, 2}] == NeumCond[t], 
         DirCond[t]}, ϕ[t,r, z], {r, z} ∈ region[t]];

The hard part is related with the two remaining boundary conditions. Indeed, besides the boundary conditions imposed by the circle, I will also have to impose a constant pressure boundary condition in the air-water interface:

  1. The boundary condition imposed by the sphere.

    I have to impose the component of the velocity which is perpendicular to the surface of the circle. I should probably use NeumannValue again, but I could not find a good example where the condition is time-dependent.

  2. The boundary condition on the air-liquid interface.

    I will probably have to split the interface into a group of evenly separated points, which will define the interface. Then, creating some sort of an iteration method, I would check the pressure at those points and correct their vertical positions such that the pressure at that point is the same as atmospheric pressure. I would also use the pressure to calculate the upward force applied in the circle, through an integration over its boundary, that would allow me to obtain the total force applied in the circle, and so predict its future position.
    If well done, the simulation would show the well-known ripples produced in this kind of situations. Maybe one solution would be to convert the upper boundary (z=4 at t=0) into a mesh, but I could not understand how this is done in Mathematica, and I further do not know if it is possible to control each of the points individually (maybe there is an easier, higher-level way?).


The boundary condition due to the circle is always changing, and this causes the boundary condition of the air-liquid interface to change as well. I have to solve the Laplace equation inside the specified time-dependent region and then I have to use the resulting potential to calculate its gradient and its time derivative, so that it is possible to get the velocity and pressure at all points. Once the pressure is known the new region can be defined, and the whole process should repeat itself.
The goal of the code is to predict the behaviour of the air-liquid interface.

I hope I made everything clear. I need help especially in the air-liquid interface. Any suggestions and/or tips are very welcome.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.