A rigid sphere floating in an hydrostatic tank

This example demonstrates the API of running a rigid fluid coupling problem in PySPH. To run it one may do:

$ cd ~/pysph/pysph/examples/rigid_body/
$ python sphere_in_vessel_akinci.py

There are many command line options that this example provides, check them out with:

$ python sphere_in_vessel.py -h

The example source can be seen at sphere_in_vessel.py.

This example demonstrates:

  • Setting up a simulation involving rigid bodies and fluid
  • Discuss mainly about rigid fluid coupling

It is divided in to three parts:

  • Create particles
  • Create equations
  • Run the application

Create particles

In this example, we have a tank with a resting fluid and a sphere falling into the tank. Create three particle arrays, tank, fluid and cube. tank and fluid has to obey wcsph scheme, where as cube has to obey rigid body equations.

def create_particles(self):
    # elided
    fluid = get_particle_array_wcsph(x=xf, y=yf, h=h, m=m, rho=rho,
                                     name="fluid")

    # elided
    tank = get_particle_array_wcsph(x=xt, y=yt, h=h, m=m, rho=rho,
                                    rad_s=rad_s, V=V, name="tank")
    for name in ['fx', 'fy', 'fz']:
        tank.add_property(name)

    cube = get_particle_array_rigid_body(x=xc, y=yc, h=h, m=m, rho=rho,
                                         rad_s=rad_s, V=V, cs=cs,
                                         name="cube")

    return [fluid, tank, cube]

We will discuss the reason for adding the properties \(fx\), \(fy\), \(fz\) to the tank particle array. The next step is to setup the equations.

Create equations

def create_equations(self):
  equations = [
      Group(equations=[
          BodyForce(dest='cube', sources=None, gy=-9.81),
      ], real=False),
      Group(equations=[
          SummationDensity(
              dest='fluid',
              sources=['fluid'], ),
          SummationDensityBoundary(
              dest='fluid', sources=['tank', 'cube'], fluid_rho=1000.0)
      ]),

      # Tait equation of state
      Group(equations=[
          TaitEOSHGCorrection(dest='fluid', sources=None, rho0=self.ro,
                              c0=self.co, gamma=7.0),
      ], real=False),
      Group(equations=[
          MomentumEquation(dest='fluid', sources=['fluid'],
                           alpha=self.alpha, beta=0.0, c0=self.co,
                           gy=-9.81),
          AkinciRigidFluidCoupling(dest='fluid',
                                   sources=['cube', 'tank']),
          XSPHCorrection(dest='fluid', sources=['fluid', 'tank']),
      ]),
      Group(equations=[
          RigidBodyCollision(dest='cube', sources=['tank'], kn=1e5)
      ]),
      Group(equations=[RigidBodyMoments(dest='cube', sources=None)]),
      Group(equations=[RigidBodyMotion(dest='cube', sources=None)]),
  ]
  return equations

A few points to note while dealing with Akinci formulation,

  1. As a first point, while computing the density of the fluid due to solid, make sure to use SummationDensityBoundary, because usual SummationDensity computes density by considering the mass of the particle, where as SummationDensityBoundary will compute it by considering the volume of the particle. This makes a lot of difference while dealing with heavy density variation flows.

  2. Apply TaitEOSHGCorrection so that there is no negative pressure.

  3. The force from the boundary (here it is tank) on fluid is computed using AkinciRigidFluidCoupling equation, but in a usual case we do it using the momentum equation. There are a few advantages by doing this. If we are computing the boundary force using the momentum equation, then one should compute the density of the boundary, then compute the pressure. Using such pressure we will compute the force. But using AkinciRigidFluidCoupling we don’t need to compute the pressure of the boundary because the force is dependent only on the fluid particle’s pressure.

    def loop(self, d_idx, d_m, d_rho, d_au, d_av, d_aw,  d_p,
             s_idx, s_V, s_fx, s_fy, s_fz, DWIJ, s_m, s_p, s_rho):
        # elide
        d_au[d_idx] += -psi * _t1 * DWIJ[0]
        d_av[d_idx] += -psi * _t1 * DWIJ[1]
        d_aw[d_idx] += -psi * _t1 * DWIJ[2]
    
        s_fx[s_idx] += d_m[d_idx] * psi * _t1 * DWIJ[0]
        s_fy[s_idx] += d_m[d_idx] * psi * _t1 * DWIJ[1]
        s_fz[s_idx] += d_m[d_idx] * psi * _t1 * DWIJ[2]
    

    Since in AkinciRigidFluidCoupling (more in next point) we compute both force on fluid by solid particle and force on solid by fluid particle, which makes our sources to hold the properties fx, fy and fz.

  4. Here first few equations deal with the simulation of fluid in hydrostatic tank. The equation dealing with rigid fluid coupling is AkinciRigidFluidCoupling . Coupling equation will deal with forces exerted by fluid on solid body, and forces exerted by solid on fluid. We find the force on fluid by solid and force on the solid by fluid in a singe equation.

    Usually in an SPH equation, we tend to change properties only of a destination particle array, but in this case, both destination and sources properties are manipulated.

  5. The final equations deal with the dynamics of rigid bodies, which are discussed in other example files.

Run the application

Finally run the application by

if __name__ == '__main__':
    app = RigidFluidCoupling()
    app.run()