Compresion test with contant pressure in a cylinder packing

Asked by Pablo Cid on 2018-12-04

Hi everyone!
I'm working on a project where as the title says, I need to apply a constant pressure over a cylinder packing.
In order to simulate friction with the cylinder walls I seted two particle taggs, one for the ones with fixed position (Walls, Tag=12)) with the setParticleNonDynamic(..) command and another to the particles inside(tag=1). I tried the Wall Piston Runnable
http://www.esys-particle.org/wiki/pmwiki.php/Main/WallPistonRunnable
with a third tag (tag=2) but at the step when it start working my wall and particles just "explode" as when you have a big time step size.

Is there any way to make kind of a circular piston or maybe other idea that archive what I need?
Here is my script. And I apologize by my poor English, I hope it's not that bad.

from esys.lsm import *
from esys.lsm.util import *
from esys.lsm.geometry import *
from LHDRun import LHD
from PistonRun import PistonRun

sim = LsmMpi(numWorkerProcesses=4, mpiDimList=[2, 2, 1])

Rmin=0.005
Rmax=0.005
minPoint = Vec3(-2, -2, 0)
maxPoint = Vec3(2, 2, 1.2)
Time=15
TimeBCkP=0.1
TimeStepSize=1e-3

sim.initNeighbourSearch(
    particleType="RotSphere",
    gridSpacing=2.5*Rmax,
    verletDist=0.2*Rmin
)
NumTimeSteps=int(Time/TimeStepSize)
sim.setNumTimeSteps(NumTimeSteps)
sim.setTimeStepSize(TimeStepSize)
sim.setSpatialDomain(
   BoundingBox(minPoint, maxPoint)
)
sim.createInteractionGroup(
   GravityPrms(name="earth-gravity",
               acceleration=Vec3(0, 0, -9.81)
   )
)
sim.readGeometry('particles.geo')
sim.setParticleDensity(
       tag = 1,
       mask = -1,
       Density = 2600
)
sim.setParticleDensity(
       tag = 12,
       mask = -1,
       Density = 2600
)
sim.setParticleDensity(
       tag = 2,
       mask = -1,
       Density = 2600
)
sim.setParticleNonDynamic(tag=12)
sim.setParticleNonRotational(tag=12)
sim.createInteractionGroup(
    FrictionPrms(
      name="friccion",
      youngsModulus=100000,
      poissonsRatio=0.25,
      dynamicMu=0.4,
      staticMu=0.6
   )
)
sim.createWall(
   name='disco',
   posn=Vec3(0,0,0.9),
   normal=Vec3(0,0,-1)
)
sim.createInteractionGroup(
   NRotBondedWallPrms(
      name='prensa',
      wallName='disco',
      normalK=100000,
      particleTag=2
   )
)
sim.createInteractionGroup(
   LinDampingPrms(
      name="linDamping",
      viscosity=0.1,
      maxIterations=100
   )
)
sim.createInteractionGroup(
   RotDampingPrms(
      name="Damping2",
      viscosity=0.002,
      maxIterations=50
   )
)

sim.createCheckPointer(
   CheckPointPrms(
      fileNamePrefix="prensa",
      beginTimeStep=0,
      endTimeStep=NumTimeSteps,
      timeStepIncr=int(float(TimeBCkP/TimeStepSize))
   )
)

sim.createFieldSaver(
 WallVectorFieldSaverPrms(
  wallName='disco',
  fieldName='Force',
  fileName='out_force.dat',
  fileFormat='RAW_SERIES',
  beginTimeStep=0,
  endTimeStep=sim.getNumTimeSteps(),
  timeStepIncr=int(float(TimeBCkP/TimeStepSize))
 )
)
sim.createFieldSaver(
 WallVectorFieldSaverPrms(
  wallName='disco',
  fieldName='Position',
  fileName='Posicion_disco.dat',
  fileFormat='RAW_SERIES',
  beginTimeStep=0,
  endTimeStep=sim.getNumTimeSteps(),
  timeStepIncr=int(float(TimeBCkP/TimeStepSize))
 )
)
#Runnable (this one works perfectly)
lhd=LHD(sim=sim)
sim.addPostTimeStepRunnable(lhd)

#Runnable (here I have some problems)
piston=PistonRun(
   sim=sim,
   wallName='disco',
   wallInteractionName='prensa',
   rampTime=10/sim.getTimeStepSize(),
   pressure=10000,
   startTime=1/sim.getTimeStepSize()
)
sim.addPostTimeStepRunnable(piston)

sim.run()

Question information

Language:
English Edit question
Status:
Solved
For:
ESyS-Particle Edit question
Assignee:
No assignee Edit question
Solved by:
Dion Weatherley
Solved:
2020-01-16
Last query:
2020-01-16
Last reply:
2018-12-04
Best Dion Weatherley (d-weatherley) said : #1

Hi Pablo,

Firstly, it would appear that your timestep increment is still too large. I calculate that you will need to decrease it from 1e-3 to 3e-4 for the parameters you have chosen.

The specific equation I have used is this:

dt < 0.2 * sqrt (M_min/K_max)

where:
M_min = 4./3.*pi*density*R_min**3.
K_max = Y_max*R_max

density = 2600.
R_min = 0.005
R_max = 0.005
Y_max = 100000.

This equation means that dt < 3.3e-4. I would normally try dt < 1.e-4 to be safe.

I would definitely try this first. I have not looked carefully at the rest of your script.

Another reason why simulations might "explode" is if particles are on the wrong side of walls. Check the position and normal vectors of walls to make sure particles are always on the positive normal side of the wall.

Let us know if you have further problems.

Cheers,

Dion

Pablo Cid (marmolin) said : #2

Thanks Dion Weatherley, that solved my question.