Merged in prabhuramachandran/pysph/rigid_collision (pull request #151)

Better support for rigid-rigid interactions and nicer examples

examples/SPHERIC/moving_square.py

@@ -35,7 +35,7 @@
 obstacle_height = 1.0
 
 # Reynolds number and kinematic viscosity
-Re = 100; nu = Umax * obstacle_width/Re
+Re = 150; nu = Umax * obstacle_width/Re
 
 # Numerical setup
 nx = 50; dx = 0.20* Lx/nx
@@ -62,15 +62,15 @@ def _get_interior(x, y):
         if ( (x[i] > 0.0) and (x[i] < Lx) ):
             if ( (y[i] > 0.0) and (y[i] < Ly) ):
                 indices.append(i)
-                
+
     return indices
 
 def _get_obstacle(x, y):
     indices = []
     for i in range(x.size):
         if ( (1.0 <= x[i] <= 2.0) and (2.0 <= y[i] <= 3.0) ):
             indices.append(i)
-            
+
     return indices
 
 def _setup_particle_properties(particles, volume):
@@ -130,7 +130,7 @@ def _setup_particle_properties(particles, volume):
     fluid.m[:] = volume * rho0
     solid.m[:] = volume * rho0
     obstacle.m[:] = volume * rho0
-    
+
     # volume is set as dx^2
     fluid.V[:] = 1./volume
     solid.V[:] = 1./volume
@@ -147,28 +147,28 @@ def _setup_particle_properties(particles, volume):
 
     solid.set_output_arrays( ['x', 'y', 'rho', 'p'] )
     obstacle.set_output_arrays( ['x', 'y', 'u0', 'rho', 'p'] )
-            
+
     particles = [fluid, solid, obstacle]
-    return particles    
+    return particles
 
 def create_particles(hcp=False, **kwargs):
     "Initial distribution using Hexagonal close packing of particles"
     # create all particles
     global dx
     if hcp:
         x, y, dx, dy, xmin, xmax, ymin, ymax = uniform_distribution.uniform_distribution_hcp2D(
-            dx=dx, xmin=-ghost_extent, xmax=Lx+ghost_extent, 
+            dx=dx, xmin=-ghost_extent, xmax=Lx+ghost_extent,
             ymin=-ghost_extent, ymax=Ly+ghost_extent)
     else:
         x, y, dx, dy, xmin, xmax, ymin, ymax = uniform_distribution.uniform_distribution_cubic2D(
-            dx=dx, xmin=-ghost_extent, xmax=Lx+ghost_extent, 
+            dx=dx, xmin=-ghost_extent, xmax=Lx+ghost_extent,
             ymin=-ghost_extent, ymax=Ly+ghost_extent)
 
     x = x.ravel(); y = y.ravel()
-    
+
     # create the basic particle array
     solid = get_particle_array(name='solid', x=x, y=y)
-    
+
     # now sort out the interior from all particles
     indices = _get_interior(solid.x, solid.y)
     fluid = solid.extract_particles( indices )
@@ -218,7 +218,7 @@ class SPHERICBenchmarkAcceleration(Equation):
     We use scipy.optimize to fit the Gaussian:
 
     .. math::
-    
+
         a \exp( -\frac{(t-b)^2}{2c^2} ) + d
 
     to the SPHERIC Motion.dat file. The values for the parameters are
@@ -238,7 +238,7 @@ def loop(self, d_idx, d_ax, t=0.0):
         b = 0.525652151
         c = 0.14142151
         d = -2.55580905e-08
-        
+
         # compute the acceleration and set it for the destination
         d_ax[d_idx] = a*exp( -(t-b)*(t-b)/(2.0*c*c) ) + d
 
@@ -255,7 +255,7 @@ def loop(self, d_idx, d_ax, t=0.0):
     # phase. This is done for all local and remote particles. At the
     # end of this group, the fluid phase has the correct density
     # taking into consideration the fluid and solid
-    # particles. 
+    # particles.
     Group(
         equations=[
             SummationDensity(dest='fluid', sources=['fluid','solid','obstacle']),
@@ -288,17 +288,17 @@ def loop(self, d_idx, d_ax, t=0.0):
             # Pressure gradient terms
             MomentumEquationPressureGradient(
                 dest='fluid', sources=['fluid', 'solid','obstacle'], pb=p0),
-            
+
             # fluid viscosity
             MomentumEquationViscosity(
                 dest='fluid', sources=['fluid'], nu=nu),
-            
+
             # No-slip boundary condition. This is effectively a
             # viscous interaction of the fluid with the ghost
             # particles.
             SolidWallNoSlipBC(
                 dest='fluid', sources=['solid','obstacle'], nu=nu),
-            
+
             # Artificial stress for the fluid phase
             MomentumEquationArtificialStress(dest='fluid', sources=['fluid']),
 

examples/rigid_body/README.md

@@ -0,0 +1,19 @@
+This directory contains a bunch of examples that demonstrate rigid body
+motion and interaction both with other rigid bodies and rigid-fluid
+coupling.
+
+The demos here are only proofs of concept.  They need work to make sure
+that the physics is correct, the equations correct and produce the right
+numbers. In particular,
+
+ - the rigid_block in tank does not work without the pressure rigid body
+   equation which is incorrect.
+
+ - the formulation and parameters used for the rigid body collision is not
+   tested if it conserves energy and works correctly in all cases.  The choice
+   of parameters is currently ad-hoc.
+
+ - the rigid-fluid coupling should also be looked at a bit more carefully with
+   proper comparisons to well-known results.
+
+Right now, it looks pretty and is a reasonable start.

examples/rigid_body/dam_break3D_sph.py

@@ -0,0 +1,140 @@
+""" 3D Dam Break Over a block of text spelling "SPH".
+
+ """
+import numpy as np
+from os import pardir
+import sys
+sys.path.append(pardir)
+
+from db_geometry import DamBreak3DGeometry
+
+from pysph.base.kernels import WendlandQuintic
+from pysph.base.utils import get_particle_array_wcsph, get_particle_array_rigid_body
+
+from pysph.sph.equation import Group
+from pysph.sph.basic_equations import ContinuityEquation, XSPHCorrection
+from pysph.sph.wc.basic import TaitEOS, TaitEOSHGCorrection, MomentumEquation
+
+from pysph.solver.application import Application
+from pysph.solver.solver import Solver
+from pysph.sph.integrator import EPECIntegrator
+from pysph.sph.integrator_step import WCSPHStep
+from pysph.tools.gmsh import vtk_file_to_points
+
+from pysph.sph.rigid_body import (BodyForce, NumberDensity, RigidBodyCollision,
+    RigidBodyMoments, RigidBodyMotion,
+    RK2StepRigidBody, ViscosityRigidBody, PressureRigidBody)
+
+
+dim = 3
+
+dt = 1e-5
+tf = 2.0
+
+# parameter to chane the resolution
+dx = 0.02
+nboundary_layers=3
+hdx = 1.2
+rho0 = 1000.0
+
+geom = DamBreak3DGeometry(
+    dx=dx, nboundary_layers=nboundary_layers, hdx=hdx, rho0=rho0,
+    with_obstacle=False)
+
+def create_particles():
+    fluid, boundary = geom.create_particles()
+
+    x, y, z = vtk_file_to_points('sph.vtk.gz', cell_centers=False)
+    y -= 0.15
+    z += 0.05
+    m = np.ones_like(x)*fluid.m[0]
+    h = np.ones_like(x)*fluid.h[0]
+    rho = np.ones_like(x)*fluid.rho[0]
+
+    obstacle = get_particle_array_rigid_body(name='obstacle', x=x, y=y, z=z,
+        m=m, h=h, rho=rho, rho0=rho)
+    obstacle.total_mass[0] = np.sum(m)
+    obstacle.add_property('cs')
+    obstacle.add_property('arho')
+    obstacle.set_lb_props( obstacle.properties.keys() )
+    obstacle.set_output_arrays( ['x', 'y', 'z', 'u', 'v', 'w', 'fx', 'fy', 'fz',
+                                 'rho', 'm', 'h', 'p', 'tag', 'pid', 'gid'] )
+
+    boundary.add_property('V')
+    boundary.add_property('fx')
+    boundary.add_property('fy')
+    boundary.add_property('fz')
+
+    return [fluid, boundary, obstacle]
+
+h0 = dx * hdx
+co = 10.0 * geom.get_max_speed(g=9.81)
+
+gamma = 7.0
+alpha = 0.5
+beta = 0.0
+B = co*co*rho0/gamma
+
+# Create the application.
+app = Application()
+
+# Create the kernel
+kernel = WendlandQuintic(dim=dim)
+
+# Setup the integrator.
+integrator = EPECIntegrator(fluid=WCSPHStep(),
+                            obstacle=RK2StepRigidBody(),
+                            boundary=WCSPHStep())
+
+# Create a solver.
+solver = Solver(kernel=kernel, dim=dim, integrator=integrator, tf=tf, dt=dt,
+                adaptive_timestep=True, n_damp=0)
+
+# create the equations
+equations = [
+
+    Group(equations=[
+            BodyForce(dest='obstacle', sources=None, gz=-9.81),
+            NumberDensity(dest='obstacle', sources=['obstacle']),
+            NumberDensity(dest='boundary', sources=['boundary']),
+            ], ),
+
+    # Equation of state
+    Group(equations=[
+
+            TaitEOS(dest='fluid', sources=None, rho0=rho0, c0=co, gamma=gamma),
+            TaitEOSHGCorrection(dest='boundary', sources=None, rho0=rho0, c0=co, gamma=gamma),
+            TaitEOSHGCorrection(dest='obstacle', sources=None, rho0=rho0, c0=co, gamma=gamma),
+
+            ], real=False),
+
+    # Continuity, momentum and xsph equations
+    Group(equations=[
+
+            ContinuityEquation(dest='fluid', sources=['fluid', 'boundary', 'obstacle']),
+            ContinuityEquation(dest='boundary', sources=['fluid']),
+            ContinuityEquation(dest='obstacle', sources=['fluid']),
+
+            MomentumEquation(dest='fluid', sources=['fluid', 'boundary'],
+                             alpha=alpha, beta=beta, gz=-9.81, c0=co,
+                             tensile_correction=True),
+
+            PressureRigidBody(dest='fluid', sources=['obstacle'], rho0=rho0),
+
+            XSPHCorrection(dest='fluid', sources=['fluid']),
+
+            RigidBodyCollision(
+                dest='obstacle', sources=['boundary'], k=1.0, d=2.0, eta=0.1, kt=0.1
+            ),
+
+            ]),
+    Group(equations=[RigidBodyMoments(dest='obstacle', sources=None)]),
+    Group(equations=[RigidBodyMotion(dest='obstacle', sources=None)]),
+
+    ]
+
+# Setup the application and solver.  This also generates the particles.
+app.setup(solver=solver, equations=equations,
+          particle_factory=create_particles)
+
+app.run()