SPH simulation of liquid objects colliding with a surface using PySPH
Full code available on GitHub.
Creating an SPH simulation using PySPH involves defining the particles, equations and the solver. We will simulate how a liquid object falls onto the ground.
The simulation involves two entities: the object (created from an STL file) and the boundary. We’re going to make them collide. Excited?
The STL model I chose is a wheel. At first, I didn’t know how to import an STL model into PySPH and use it in the simulation. But then I had an epiphany: I could generate a mesh for the model and place particles in it’s nodes. Genius, isn’t it?
Generating the mesh
First of all, I generated a mesh for the STL model using gmsh.
gmsh.initialize() path = os.path.dirname(os.path.abspath(__file__)) gmsh.merge(os.path.join(path, 'wheel.stl')) angle = 40 forceParametrizablePatches = True includeBoundary = True curveAngle = 180 gmsh.model.mesh.classifySurfaces(angle * math.pi / 180., includeBoundary, forceParametrizablePatches, curveAngle * math.pi / 180.) gmsh.model.mesh.createGeometry() s = gmsh.model.getEntities(2) l = gmsh.model.geo.addSurfaceLoop([s[i] for i in range(len(s))]) gmsh.model.geo.addVolume([l]) gmsh.model.geo.synchronize() f = gmsh.model.mesh.field.add("MathEval") gmsh.model.mesh.field.setString(f, "F", "4") gmsh.model.mesh.field.setAsBackgroundMesh(f) gmsh.model.mesh.generate(3) gmsh.write('wheel.msh')
Processing the node’s coordinates
The code produces a
.msh file, which contains details about the mesh. We are going to open it and extract it’s nodes.
gmsh.initialize() gmsh.open(path) nodeTags, nodesCoord, parametricCoord = gmsh.model.mesh.getNodes()
nodesCoord is what we’re interested in. Let’s scale the coordinates (I picked a scale factor so the wheel’s not too big).
scale = 600 x = nodesCoord[0::3]/scale y = nodesCoord[2::3]/scale
Next, we’ll make our wheel be above the ground by adding the lowest value along the y axis. We additionally raise it to a defined height.
y = y + abs(min(y)) + height
The mesh is dense. I mean 60k nodes dense. We’re going to pick just 6k nodes (every 10th one). Actually, you can try picking more particles. But it would take a while to compute.
x = x[0::10] y = y[0::10]
Thus, we implemented a function to extract coordinates from an msh file. Creating the particles looks like this.
liquid_x, liquid_y = self.particles_from_model('wheel.msh') liquid = gpa(name='liquid', x=liquid_x, y=liquid_y)
Let’s make the walls not too big and not too small: to surround the object, leaving enough room. I decided the room to be 500*dx (500 particles) in each direction.
liquid_x, liquid_y = self.particles_from_model('wheel.msh') min_x = min(liquid_x) - 500*dx max_x = max(liquid_x) + 500*dx min_y = 0 max_y = max(liquid_y) + 500*dx _x = np.arange(min_x, max_x, dx) _y = np.arange(min_y, max_y, dx) x, y = np.meshgrid(_x, _y) x = x.ravel() y = y.ravel()
The walls itself are extracted from our overall coordinates array.
walls_x =  walls_y =  for i in range(x.size): if (y[i] < (min_y + wall_thickness)) or (x[i] >= (max_x - wall_thickness)) or (x[i] <= (min_x + wall_thickness)): walls_x.append(x[i]) walls_y.append(y[i])
Now, we’ve got the arrays of coordinates, where wall particles should be. Let’s generate them.
walls = gpa(name='walls', x=walls_x, y=walls_y)
We can preview our setup by plotting it.
plt.scatter(liquid_x, liquid_y) plt.scatter(walls_x, walls_y) plt.gca().set_aspect('equal') plt.show()
Finally, we define properties for the particles.
# particle volume liquid.add_property('V') walls.add_property('V') # kernel sum term for boundary particles walls.add_property('wij') # advection velocities and accelerations for name in ('auhat', 'avhat', 'awhat'): liquid.add_property(name) liquid.rho[:] = rho0 walls.rho[:] = 5*rho0 liquid.rho0[:] = rho0 walls.rho0[:] = 5*rho0 # mass is set to get the reference density of rho0 volume = dx * dx # volume is set as dx^2 liquid.V[:] = 1. / volume walls.V[:] = 1. / volume liquid.m[:] = volume * rho0 walls.m[:] = volume * rho0 # smoothing lengths liquid.h[:] = hdx * dx walls.h[:] = hdx * dx
The next step is to define the equations, which describe particle’s behavior and boundary conditions. I took them straight out of the hydrostatic_tank example in the PySPH repository (check out other examples, there’s some pretty cool stuff there). I chose the REF3 formulation.
# For the multi-phase formulation, we require an estimate of the # particle volume. This can be either defined from the particle # number density or simply as the ratio of mass to density. Group(equations=[ VolumeFromMassDensity(dest='liquid', sources=None) ], ), # Equation of state is typically the Tait EOS with a suitable # exponent gamma Group(equations=[ TaitEOS( dest='liquid', sources=None, rho0=rho0, c0=c0, gamma=gamma), TaitEOS( dest='walls', sources=None, rho0=rho0, c0=c0, gamma=gamma), ], ), # Main acceleration block. The boundary conditions are imposed by # peforming the continuity equation and gradient of pressure # calculation on the solid phase, taking contributions from the # fluid phase Group(equations=[ # Continuity equation ContinuityEquation(dest='liquid', sources=['liquid', 'walls']), ContinuityEquation(dest='walls', sources=['liquid']), # Pressure gradient with acceleration damping. MomentumEquationPressureGradient( dest='liquid', sources=['liquid', 'walls'], pb=0.0, gy=gy), # artificial viscosity for stability MomentumEquationArtificialViscosity( dest='liquid', sources=['liquid', 'walls'], alpha=0.25, c0=c0), # Position step with XSPH XSPHCorrection(dest='liquid', sources=['liquid'], eps=0.5) ])
Finally, we define the solver.
kernel = CubicSpline(dim=2) integrator = EPECIntegrator(liquid=WCSPHStep(), walls=WCSPHStep()) solver = Solver(kernel=kernel, dim=2, integrator=integrator, tf=tf, dt=dt, adaptive_timestep=True, cfl=0.5, output_at_times=output_at_times)
Running the simulation took about 30 minutes on my PC.
We’ve got some pretty slick results. Check it out!