#impactLoadCell.py:
#=================
#
# An ESyS-Particle simulation script to implement breakage
#  of a cylindrical specimen of DEM rock impacted with 
#  a spherical indenter a la the laboratory Short Impact Load Cell.
#
# D.Weatherley, UQ, 2018.
#
#Import the relevant modules: 
from esys.lsm import *
from esys.lsm.util import *
from math import sqrt
from ball_mover import SphereBodyMoverRunnable

#Fundamental units of measure:
#	length = millimetres (mm)
#	time = milliseconds (ms)
#	force = Newtons (N)

#Derived units of measure:
#	velocity = metres per second (m/s)
#	stress = megaPascals (MPa)
#	mass = grams (g)
#	density = 1.0e-6 kg/m^3 (g/mm^3)
#	acceleration = 1.0e-3 m/^2 (mm/s^2)

#Model Parameters:
#================
cylRadius = 5.0 	#radius of the cylindrical specimen of rock
cylHeight = 10.0	#length of the cylindrical specimen of rock

Rmin = 0.2		#minimum DEM particle radius
Rmax = 0.5		#maximum DEM particle radius

bondModulus = 1.0e5		#elastic modulus of beam interactions
bondPoissonRatio = 0.25		#Poisson's ratio of beam interactions
bondCohesion = 100.		#cohesive strength of beam interactions
bondTanAngle = 1.0		#tan(friction angle) of beam interactions
frictionCoeff = 0.6		#friction coefficient for frictional contacts

density = 2.6e-3		#density of individual DEM particles
gravity = Vec3(0,-9.81e-3,0)	#acceleration due to gravity

dropHeight = 200.0		#initial height of the spherical indenter
ballRadius = 25.0		#radius of the spherical indenter
ballDensity = 7.5e-3		#density of the spherical indenter (i.e. steel)

impactSpeed = sqrt(2.*9.81e-3*200.0)	#calculated impact speed

wallInteractionType = "glued"	#specifies type of wall interaction
				#(valid types: "elastic" or "glued")

numT = 100000	#total number of simulation timesteps

#calculated stable timestep increment:
dt = 0.1*sqrt(4.*density*Rmin**3./3./bondModulus/Rmax)

#Construct a simulation container:
sim = LsmMpi (
   numWorkerProcesses = 1,
   mpiDimList = [1,1,1]
)

#Initialise the neighbour search algorithm (using RotSphere):
sim.initNeighbourSearch (
   particleType = "RotSphere",
   gridSpacing = 2.5*Rmax,
   verletDist = 0.2*Rmin
)

#Specify the number of simulation timesteps and timestep increment:
sim.setNumTimeSteps (numT)
sim.setTimeStepSize (dt)

#Read the initial DEM particle geometry:
sim.readGeometry("cylinder.geo")

#Set the density of the DEM particles:
sim.setParticleDensity (tag = 1, mask = -1, Density = density)
sim.setParticleDensity (tag = 2, mask = -1, Density = density)

#Add a "SphereBody" to represent the spherical indenter:
sim.createSphereBody (
   name = "steel_ball",
   posn = Vec3 (0,2.*cylRadius+ballRadius,0.5*cylHeight),
   radius = ballRadius
)

#Add a planar wall beneath the DEM rock specimen to record forces:
sim.createWall (
   name = "floor",
   posn = Vec3(0,0,0),
   normal = Vec3(0,1,0)
)

#Initialise elastic interactions between DEM particles and the indenter:
sim.createInteractionGroup (
   NRotElasticSphereBodyPrms (
      name = "pb_repel",
      sphereName = "steel_ball",
      normalK = bondModulus
   )
)

#Initialise elastic interactions between DEM particles and the wall:
if (wallInteractionType=="elastic"):
   sim.createInteractionGroup (
      NRotElasticWallPrms (
         name = "pw_floor",
         wallName = "floor",
         normalK = bondModulus
      )
   )
elif (wallInteractionType == "glued"):
   sim.createInteractionGroup (
      NRotSoftBondedWallPrms (
         name = "pw_floor",
         wallName = "floor",
         normalK = bondModulus,
         shearK = 0.0,
         particleTag = 2,
         tagMask = -1,
         scaling = True
      )
   )

#Specify brittle-elastic beam interactions between bonded particles:
sim.createInteractionGroup (
   BrittleBeamPrms (
      name = "pp_beams",
      youngsModulus = bondModulus,
      poissonsRatio = bondPoissonRatio,
      cohesion = bondCohesion,
      tanAngle = bondTanAngle,
      tag = 1
   )
)

#Specify rotational frictional interactions between unbonded particles:
sim.createInteractionGroup (
   FrictionPrms (
      name = "pp_friction",
      youngsModulus = bondModulus,
      poissonsRatio = bondPoissonRatio,
      dynamicMu = frictionCoeff,
      staticMu = frictionCoeff
   )
)

#Ensure each particle pair only undergoes beam or frictional interactions:
sim.createExclusion ("pp_beams","pp_friction")

#Specify the acceleration due to gravity:
sim.createInteractionGroup (
   GravityPrms (
      name = "gravity",
      acceleration = gravity
   )
)

#Attach a Runnable to move the indenter at constant speed:
ball_mover = SphereBodyMoverRunnable (
   sim,
   "steel_ball",
   speed = impactSpeed,
   dirn = Vec3 (0,-1,0)
)
sim.addPostTimeStepRunnable(ball_mover)

#Add a checkpointer to record simulation data for visualisation:
sim.createCheckPointer (
   CheckPointPrms (
      fileNamePrefix = "snapshot",
      beginTimeStep = 0,
      endTimeStep = numT,
      timeStepIncr = 1000
   )
)

#Add a FieldSaver to record wall forces:
sim.createFieldSaver (
   WallVectorFieldSaverPrms (
      wallName = ["floor"],
      fieldName = "Force",
      fileName = "wall_force.csv",
      fileFormat = "RAW_SERIES",
      beginTimeStep = 0,
      endTimeStep = numT,
      timeStepIncr = 1
   )
)

#Add a FieldSaver to record total particle kinetic energy each timestep:
sim.createFieldSaver (
   ParticleScalarFieldSaverPrms (
      fieldName="e_kin",
      fileName="Ekin.csv",
      fileFormat="SUM",
      beginTimeStep = 0,
      endTimeStep = numT,
      timeStepIncr = 1
   )
)

#Add a FieldSaver to record the number of bonded interactions each timestep:
sim.createFieldSaver (
   InteractionScalarFieldSaverPrms (
      interactionName="pp_beams",
      fieldName="count",
      fileName="Nbonds.csv",
      fileFormat="SUM",
      beginTimeStep = 0,
      endTimeStep = numT,
      timeStepIncr = 1
   )
)

#Add a FieldSaver to record the total stored elastic strain energy:
sim.createFieldSaver (
   InteractionScalarFieldSaverPrms (
      interactionName="pp_beams",
      fieldName="potential_energy",
      fileName="Epot.csv",
      fileFormat="SUM",
      beginTimeStep = 0,
      endTimeStep = numT,
      timeStepIncr = 1
   )
)

#execute the simulation:
sim.run()
GlossyBlue theme adapted by David Gilbert
Powered by PmWiki
www.000webhost.com