Impact simulation

Asked by Temitope on 2017-06-08

Hi there,

Please I am trying to impact a steel ball on the surface area of a cylindrical rod. I want to get the force-time history during the impact and also need an output file to view in paraview. Please I need help with the following:
1) Specification of material micro properties such as young modulus and poisson's ratio for the steel-ball and steel rod . The impactor is a steel and my rod is also a steel
2)Inserting the cylindrical (0.02m diamter and 1.5m length) rod into the domain
3)Basis of getting the normalK based on the material properties.

i was able to modify part of the tutorial code to suit my work. Here is the code and response that I am getting after attempting to run it.

#Import the appropriate Esys-Particle module
from esys.lsm import*
from esys.lsm.util import Vec3, BoundingBox
from esys.lsm.geometry import*
#instantiate a simulation object
#number of processor for simulation
sim = LsmMpi (numWorkerProcesses=1, mpiDimList=[1,1,1])
#initialise the neighbour search algorithm:
sim.initNeighbourSearch(
 particleType="NRotSphere",
 gridSpacing=0.0625, #i.e, 2.5 Xply by max. particle radius which is 0.025m
 verletDist=0.005 #i.e, 0.2 Xply by min. particle radius which is 0.025m
)
#setting number of timestep and timestep size
sim.setNumTimeSteps(1000)
sim.setTimeStepSize(0.001)

#Domain of simulation: here is the cube with 40m each size
domain = BoundingBox(Vec3(0,0,0), Vec3(0.2,0.5,0.2))
sim.setSpatialDomain (domain)

#create a steel ball with NRotSphere type, id=0 position at that vector posn, with the specified radius and mass
steel_ball=NRotSphere(id=0, posn=Vec3(0.1,0.3,0.1), radius=0.025, mass=0.51)

#initialise the velocity of the steel_ball Vx, Vy and Vz (m/s)
steel_ball.setLinearVelocity(Vec3(0,0,0))

#add the newly created steel_ball
sim.createParticle(steel_ball)

#initialise gravity in the domain:
sim.createInteractionGroup(
 GravityPrms(name="Gravitational_acceleration", acceleration=Vec3(0,-9.81,0))
)
#add a horizontal wall to act like the steel rod surface:
sim.createWall(
 name="steel_rod",
 posn=Vec3(0,0,0),
 normal=Vec3(0,1,0)
)
#specify the type of interaction between wall and steel ball
sim.createInteractionGroup(
 NRotElasticWallPrms(
  name="elastic_wall",
  wallName="steel_rod",
  normalK=10000.0
 )
)
#add local viscosity to simulate air resistance
sim.createInteractionGroup(
 LinDampingPrms(
  name="linDamping",
  viscosity=0.1,
  maxIterations=100
 )
)
#create a FieldSaver to store forces on the steelrod:
force_saver=WallVectorFieldSaverPrms(
 wallName=["steel_rod"],
 fieldName="Force",
 fileName="Steel_rod_Force.dat",
 YbeginTimeStep=0,
 endTimeStep=100,
 timeStepIncr=0.001
)
sim.createFieldSaver(force_saver)

#add check pointer to store simulation data for paraview
sim.createCheckPointer (
 CheckPointPrms (
  fileNamePrefix ="snapshot",
  beginTimeStep =0,
  endTime=1000,
  timeStepIncr=4
 )
)
#run the simulation
sim.run()

adminuser@temitope-oladele:~$ export PATH=/usr/local/bin/:$PATH
adminuser@temitope-oladele:~$ export LD_LIBRARY_PATH=/usr/local/lib/:$LD_LIBRARY_PATH
adminuser@temitope-oladele:~$ export LIBRARY_PATH=/usr/local/lib/:$LIBRARY_PATH
adminuser@temitope-oladele:~$ export PYTHONPATH=/usr/local/lib/python2.7/dist-packages/:$PYTHONPATH
adminuser@temitope-oladele:~$ export PYTHONPATH=/usr/local/lib/python2.7/dist-packages/gengeo/:$PYTHONPATH
adminuser@temitope-oladele:~$ mpirun -np 2 esysparticle calibration.py
CSubLatticeControler::initMPI()
slave started at local/global rank 0 / 1
Traceback (most recent call last):
  File "calibration.py", line 71, in <module>
    timeStepIncr=4
Boost.Python.ArgumentError: Python argument types in
    CheckPointPrms.__init__(CheckPointPrms)
did not match C++ signature:
    __init__(_object*, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fileNamePrefix, int beginTimeStep, int endTimeStep, int timeStepIncr)

Kind Regards
Tope

Question information

Language:
English Edit question
Status:
Solved
For:
ESyS-Particle Edit question
Assignee:
No assignee Edit question
Solved by:
Temitope
Solved:
2017-07-19
Last query:
2017-07-19
Last reply:
2017-07-10

This question was reopened

Vince Boros (v-boros) said : #1

You have an incorrect parameter in CheckPointPrms(). Change 'endTime' to 'endTimeStep' to avoid the error message you received.

Regards,

Vince

Vince Boros (v-boros) said : #2

I noticed another mistake in the code. In WallVectorFieldSaverPrms() there is a spurious 'Y' at the start of the parameter 'beginTimeStep'.

Vince

Temitope (temitope) said : #3

Hi Vince,

I did correct it . Here is the message that I am getting after the correction.

adminuser@temitope-oladele:~$ export PATH=/usr/local/bin/:$PATH
adminuser@temitope-oladele:~$ export LD_LIBRARY_PATH=/usr/local/lib/:$LD_LIBRARY_PATH
adminuser@temitope-oladele:~$ export LIBRARY_PATH=/usr/local/lib/:$LIBRARY_PATH
adminuser@temitope-oladele:~$ export PYTHONPATH=/usr/local/lib/python2.7/dist-packages/:$PYTHONPATH
adminuser@temitope-oladele:~$ export PYTHONPATH=/usr/local/lib/python2.7/dist-packages/gengeo/:$PYTHONPATH
adminuser@temitope-oladele:~$ mpirun -np 2 esysparticle calibration.py
CSubLatticeControler::initMPI()
slave started at local/global rank 0 / 1
Traceback (most recent call last):
  File "calibration.py", line 70, in <module>
    timeStepIncr=0.001
Boost.Python.ArgumentError: Python argument types in
    WallVectorFieldSaverPrms.__init__(WallVectorFieldSaverPrms)
did not match C++ signature:
    __init__(_object*, boost::python::list wallName, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fieldName, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fileName, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fileFormat, int beginTimeStep, int endTimeStep, int timeStepIncr)
    __init__(_object*, boost::python::tuple wallName, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fieldName, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fileName, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fileFormat, int beginTimeStep, int endTimeStep, int timeStepIncr)
    __init__(_object*, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > wallName, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fieldName, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fileName, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fileFormat, int beginTimeStep, int endTimeStep, int timeStepIncr)

Kind Regards
Tope

Dion Weatherley (d-weatherley) said : #4

Hi Tope,

Just checking that you replaced "YbeginTimeStep" with "beginTimeStep" on line number 68, as suggested by Vince.

If so, please post your new "calibration.py" script in full, so we can look for other possible errors.

Cheers,

Dion

Temitope (temitope) said : #5

Hi Dion,

Here is my script as requested:

#Import the appropriate Esys-Particle module
from esys.lsm import*
from esys.lsm.util import Vec3, BoundingBox
from esys.lsm.geometry import*
#instantiate a simulation object
#number of processor for simulation
sim = LsmMpi (numWorkerProcesses=1, mpiDimList=[1,1,1])
#initialise the neighbour search algorithm:
sim.initNeighbourSearch(
 particleType="NRotSphere",
 gridSpacing=0.0625, #i.e, 2.5 Xply by max. particle radius which is 0.025m
 verletDist=0.005 #i.e, 0.2 Xply by min. particle radius which is 0.025m
)
#setting number of timestep and timestep size
sim.setNumTimeSteps(1000)
sim.setTimeStepSize(0.001)

#Domain of simulation: here is the cube with 40m each size
domain = BoundingBox(Vec3(0,0,0), Vec3(0.2,0.5,0.2))
sim.setSpatialDomain (domain)

#create a steel ball with NRotSphere type, id=0 position at that vector posn, with the specified radius and mass
steel_ball=NRotSphere(id=0, posn=Vec3(0.1,0.3,0.1), radius=0.025, mass=0.51)

#initialise the velocity of the steel_ball Vx, Vy and Vz (m/s)
steel_ball.setLinearVelocity(Vec3(0,0,0))

#add the newly created steel_ball
sim.createParticle(steel_ball)

#initialise gravity in the domain:
sim.createInteractionGroup(
 GravityPrms(name="Gravitational_acceleration", acceleration=Vec3(0,-9.81,0))
)
#add a horizontal wall to act like the steel rod surface:
sim.createWall(
 name="steel_rod",
 posn=Vec3(0,0,0),
 normal=Vec3(0,1,0)
)
#specify the type of interaction between wall and steel ball
sim.createInteractionGroup(
 NRotElasticWallPrms(
  name="elastic_wall",
  wallName="steel_rod",
  normalK=10000.0
 )
)
#add local viscosity to simulate air resistance
sim.createInteractionGroup(
 LinDampingPrms(
  name="linDamping",
  viscosity=0.1,
  maxIterations=100
 )
)
#create a FieldSaver to store forces on the steelrod:
force_saver=WallVectorFieldSaverPrms(
 wallName=["steel_rod"],
 fieldName="Force",
 fileName="Steel_rod_Force.dat",
 beginTimeStep=0,
 endTimeStep=100,
 timeStepIncr=0.001
)
sim.createFieldSaver(force_saver)

#add check pointer to store simulation data for paraview
sim.createCheckPointer (
 CheckPointPrms (
  fileNamePrefix ="snapshot",
  beginTimeStep =0,
  endTimeStep=1000,
  timeStepIncr=4
 )
)
#run the simulation
sim.run()

Kind Regards
Tope

Dion Weatherley (d-weatherley) said : #6

Hi Tope,

OK, I missed this earlier: you need to add an extra line specifying the "fileFormat" for your force_saver. Replace your current definition of the force_saver with the following:

force_saver=WallVectorFieldSaverPrms(
 wallName=["steel_rod"],
 fieldName="Force",
 fileName="Steel_rod_Force.dat",
 fileFormat="RAW_SERIES",
 beginTimeStep=0,
 endTimeStep=100,
 timeStepIncr=0.001
)

Hopefully your script will work now. Have fun!

Cheers,

Dion

Temitope (temitope) said : #7

Hi Dion,

I have added the line of line and here is what i am getting:

adminuser@temitope-oladele:~$ export PATH=/usr/local/bin/:$PATH
adminuser@temitope-oladele:~$ export LD_LIBRARY_PATH=/usr/local/lib/:$LD_LIBRARY_PATH
adminuser@temitope-oladele:~$ export LIBRARY_PATH=/usr/local/lib/:$LIBRARY_PATH
adminuser@temitope-oladele:~$ export PYTHONPATH=/usr/local/lib/python2.7/dist-packages/:$PYTHONPATH
adminuser@temitope-oladele:~$ export PYTHONPATH=/usr/local/lib/python2.7/dist-packages/gengeo/:$PYTHONPATH
adminuser@temitope-oladele:~$ mpirun -np 2 esysparticle calibration.py
CSubLatticeControler::initMPI()
slave started at local/global rank 0 / 1
Traceback (most recent call last):
  File "calibration.py", line 71, in <module>
    timeStepIncr=0.001
Boost.Python.ArgumentError: Python argument types in
    WallVectorFieldSaverPrms.__init__(WallVectorFieldSaverPrms)
did not match C++ signature:
    __init__(_object*, boost::python::list wallName, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fieldName, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fileName, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fileFormat, int beginTimeStep, int endTimeStep, int timeStepIncr)
    __init__(_object*, boost::python::tuple wallName, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fieldName, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fileName, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fileFormat, int beginTimeStep, int endTimeStep, int timeStepIncr)
    __init__(_object*, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > wallName, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fieldName, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fileName, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > fileFormat, int beginTimeStep, int endTimeStep, int timeStepIncr)

here is my code:

#Calibration of SILC device: Steelball-steelrod impact
#Drop a 0.51kg steel ball from 0.3m height on a steel rod
#Author: Oladele T.P (2017)
#Acknowledged ESyS-particle tutorial

#Import the appropriate Esys-Particle module
from esys.lsm import*
from esys.lsm.util import Vec3, BoundingBox
from esys.lsm.geometry import*
#instantiate a simulation object
#number of processor for simulation
sim = LsmMpi (numWorkerProcesses=1, mpiDimList=[1,1,1])
#initialise the neighbour search algorithm:
sim.initNeighbourSearch(
 particleType="NRotSphere",
 gridSpacing=0.0625, #i.e, 2.5 Xply by max. particle radius which is 0.025m
 verletDist=0.005 #i.e, 0.2 Xply by min. particle radius which is 0.025m
)
#setting number of timestep and timestep size
sim.setNumTimeSteps(1000)
sim.setTimeStepSize(0.001)

#Domain of simulation: here is the cube with 40m each size
domain = BoundingBox(Vec3(0,0,0), Vec3(0.2,0.5,0.2))
sim.setSpatialDomain (domain)

#create a steel ball with NRotSphere type, id=0 position at that vector posn, with the specified radius and mass
steel_ball=NRotSphere(id=0, posn=Vec3(0.1,0.3,0.1), radius=0.025, mass=0.51)

#initialise the velocity of the steel_ball Vx, Vy and Vz (m/s)
steel_ball.setLinearVelocity(Vec3(0,0,0))

#add the newly created steel_ball
sim.createParticle(steel_ball)

#initialise gravity in the domain:
sim.createInteractionGroup(
 GravityPrms(name="Gravitational_acceleration", acceleration=Vec3(0,-9.81,0))
)
#add a horizontal wall to act like the steel rod surface:
sim.createWall(
 name="steel_rod",
 posn=Vec3(0,0,0),
 normal=Vec3(0,1,0)
)
#specify the type of interaction between wall and steel ball
sim.createInteractionGroup(
 NRotElasticWallPrms(
  name="elastic_wall",
  wallName="steel_rod",
  normalK=10000.0
 )
)
#add local viscosity to simulate air resistance
sim.createInteractionGroup(
 LinDampingPrms(
  name="linDamping",
  viscosity=0.1,
  maxIterations=100
 )
)
#create a FieldSaver to store forces on the steelrod:
force_saver=WallVectorFieldSaverPrms(
 wallName=["steel_rod"],
 fieldName="Force",
 fileName="Steel_rod_Force.dat",
 fileFormat="RAW_SERIES",
 beginTimeStep=0,
 endTimeStep=100,
 timeStepIncr=0.001
)
sim.createFieldSaver(force_saver)

#add check pointer to store simulation data for paraview
sim.createCheckPointer (
 CheckPointPrms (
  fileNamePrefix ="snapshot",
  beginTimeStep =0,
  endTimeStep=1000,
  timeStepIncr=4
 )
)
#run the simulation
sim.run()

Dion Weatherley (d-weatherley) said : #8

Hi Tope,

I finally found the problem! Try changing timeStepIncr from 0.001 to simply 1. On my PC your script executes once this is changed.

In your script, your wall field saver is given a fractional timeStepIncr. This needs to be an integer >=1.

Let me explain:
* 'beginTimeStep' : the timestep number to begin saving data
* 'endTimeStep' : the timestep number to end saving data
* 'timeStepIncr': the number of timesteps to ignore between each successive data output step

so if I set these to 0, 1000, 100 resp., I will record simulation data at timesteps 0,100,200,..,900,1000.

Note that each timestep equates to an increment of the simulation by dt = timeStepSize (in your script dt = 0.001).

Cheers,

Dion

Temitope (temitope) said : #9

Hi Dion,

it worked!

Thanks a lot

Tope

Temitope (temitope) said : #10

Thanks Dion Weatherley, that solved my question.

Temitope (temitope) said : #11

Hi Dion,

I have few questions in creating a moving object such as a wall or steel ball. Their properties (young's modulus, poisson's ratio, density etc) are different from the properties of the particles. How do I include impacting a spherical steel ball on my particles? And how do I describe their properties in ESyS particle?

Tope

Hi Tope,

This question requires quite some in-depth answers that might be best done via phone/email. I'm happy to collaborate with you to help with this issue. Can you please email me to arrange a time to discuss this?

For the benefit of others, there are currently 2 ways a spherical impactor might be simulated:
1) As a bonded assembly of spheres whose bond parameters are tuned to the mechanical properties of steel,
2) As a so-called "SphereBody" which is essentially a spherical wall (akin to a planar wall) that may be moved via a Runnable.

The first method is arguably better because deformation of the ball is simulated but it can result in simulation models with a very large number of DEM spheres. The second method has been successfully used to model the impact of a steel ball upon a rock, reproducing nicely the deformation and breakage of rocks in the so-called Short Impact Load Cell apparatus. The downside of the second method is that the steel ball is rigid so phenomena associated with deformation of the ball (e.g. elastic restitution) is not simulated. In any case, if the motion of the ball is carefully prescribed, this is a computationally efficient method to simulate a large spherical indentor's effect upon a bonded particle assembly.

Cheers,

Dion

Temitope (temitope) said : #13

Hi Dion,

I have been thinking about the second option but have not been able to include the properties. The first option seems to be a better idea since mechanical properties of the steel can be defined there. I have sent you the email as requested.

Kind Regards
Tope

Temitope (temitope) said : #14

Hi Dion and Vince,

I attempted making my impactor (steel ball) using gengeo. I tried implementing the Sphere argument in a similar way to the simple_box tutorial. I not sure of where i went wrong. Please see my script and the response i got after execution.

#particle cluster for steelball geometry

from gengeo import*

#-- parameters --

#position of the steel sphere
xdim=0
ydim=0.325
zdim=0

#Steelball radius
R= 0.05
#particle size range
minRadius = 0.0001 #1micron
maxRadius = 0.0001

centre = Vector3(xdim,ydim,zdim)
radius = R

#Steel ball volume
steel_ball = Sphere(centre, radius)

#neighbour table
mntable =MNTable3D(
 centre=centre,
 radius=radius,
 gridSize=2.5*maxRadius,
 numGroups=1
 )

#setup packer and iteration parameters
insertFails=1000
maxIter=1000
tol =1.0e-6

#packer; place particles into volume
packer=InsertGenerator3D(centre,radius,insertFails,maxIter,tol)
#generte the packerinto volume
packer.generatePacking(steel_ball, mntable)

#create bonds between between neighbouring paricles:
mntable.write("steel_ball.geo",1)
mntable.write("steel_ball.vtu",2)

adminuser@temitope-oladele:~$ python steel_ball.py
python: can't open file 'steel_ball.py': [Errno 2] No such file or directory
adminuser@temitope-oladele:~$ python steelball.py
Traceback (most recent call last):
  File "steelball.py", line 29, in <module>
    numGroups=1
Boost.Python.ArgumentError: Python argument types in
    MNTable3D.__init__(MNTable3D)
did not match C++ signature:
    __init__(_object*, Vector3 {lvalue} minPoint, Vector3 {lvalue} maxPoint, double gridSize, unsigned int numGroups=1)
    __init__(_object*, MNTable3D)
    __init__(_object*)
adminuser@temitope-oladele:~$ python steelball.py
Traceback (most recent call last):
  File "steelball.py", line 29, in <module>
    numGroups=0
Boost.Python.ArgumentError: Python argument types in
    MNTable3D.__init__(MNTable3D)
did not match C++ signature:
    __init__(_object*, Vector3 {lvalue} minPoint, Vector3 {lvalue} maxPoint, double gridSize, unsigned int numGroups=1)
    __init__(_object*, MNTable3D)
    __init__(_object*)

Kind Regards
Tope

Hi Tope,

There are 3 problems with your gengeo script:

1) the volume to be filled should be a "SphereVol" rather than a "Sphere":
steel_ball = SphereVol(centre, radius)

2) the MNTable3D is given incorrect arguments. You need to supply the minPoint and maxPoint of a rectangular prism surrounding the sphere:
mntable = MNTable3D(
        minPoint = Vector3(xdim-1.25*radius,ydim-1.25*radius,zdim-1.25*radius),
        maxPoint = Vector3(xdim+1.25*radius,ydim+1.25*radius,zdim+1.25*radius),
        gridSize=2.5*maxRadius,
        numGroups=1
        )

3) the InsertGenerator3D is also incorrectly specified. You need to supply the minRadius and maxRadius of DEM particles to pack inside the sphere:
packer=InsertGenerator3D(minRadius,maxRadius,insertFails,maxIter,tol,True)

I also recommend you increase minRadius and maxRadius to around 0.001 (10micron), at least for initial testing of your simulations. This results in a sphere comprised of approximately 85,000 particles. If you use 1 micron particles, I expect your sphere will contain upwards of 1 million particles, requiring of order 100 CPUs to conduct a simulation in a reasonable amount of time. This is far too many particles for initial design and testing of simulation scripts.

Cheers,

Dion

Temitope (temitope) said : #16

Hi Dion,

Thanks a lot it worked!

Please I also use this means to give you a kind reminder about the email I sent to you.

Kind Regards
Tope

Temitope (temitope) said : #17

Thanks Dion Weatherley, that solved my question.

Temitope (temitope) said : #18

Hi Dion,

I am happy with the gengeo output of my sphere but I tried to import this gengeo file into my impact simulation script. I'm not so sure of what i'm not doing correctly:

#Import the appropriate Esys-Particle module
from esys.lsm import*
from esys.lsm.util import Vec3, BoundingBox
from esys.lsm.geometry import*
from gengeo import *

from steel_ball import*

#instantiate a simulation object
#number of processor for simulation
sim = LsmMpi (numWorkerProcesses=1, mpiDimList=[1,1,1])
#initialise the neighbour search algorithm:
sim.initNeighbourSearch(
 particleType="NRotSphere",
 gridSpacing=0.0025, #i.e, 2.5 Xply by max. particle radius which is 0.001m
 verletDist=0.00002 #i.e, 0.2 Xply by min. particle radius which is 0.0001m
)
#setting number of timestep and timestep size
sim.setNumTimeSteps(1000)
sim.setTimeStepSize(0.001)

#Domain of simulation: here is the cube with
domain = BoundingBox(Vec3(0,0,0), Vec3(0.2,0.5,0.2))
sim.setSpatialDomain (domain)

#initialise the velocity of the steel_ball Vx, Vy and Vz (m/s)
steel_ball.setLinearVelocity(Vec3(0,0,0))

#initialise gravity in the domain:
sim.createInteractionGroup(
 GravityPrms(name="Gravitational_acceleration", acceleration=Vec3(0,-9.81,0))
)
#add a horizontal wall to act like the steel rod surface:
sim.createWall(
 name="steel_rod",
 posn=Vec3(0,0,0),
 normal=Vec3(0,1,0)
)
#specify the type of interaction between wall and steel ball
#normalK calculated from steel properties i.e function of young modulus (wikipedia)
sim.createInteractionGroup(
 NRotElasticWallPrms(
  name="elastic_wall",
  wallName="steel_rod",
  normalK= 10000
 )
)

#create rotational elastic-brittle bonds between between steel_ball particles:
pp_bonds = sim.createInteractionGroup (
 BrittleBeamPrms(
  name="Steel_ball_bonds",
  youngsModulus=210000000000.0, #210Gpa for steel
  poissonsRatio=0.29,
  cohesion=100.0, #check this again!
  tanAngle=1.0,
  tag=1
 )
)

#add local viscosity to simulate air resistance
sim.createInteractionGroup(
 LinDampingPrms(
  name="linDamping",
  viscosity=0.1,
  maxIterations=100
 )
)
#create a FieldSaver to store forces on the steelrod:
force_saver=WallVectorFieldSaverPrms(
 wallName=["steel_rod"],
 fieldName="Force",
 fileName="Steel_rod_Force.dat",
 fileFormat="RAW_SERIES",
 beginTimeStep=0,
 endTimeStep=1000,
 timeStepIncr=1
)
sim.createFieldSaver(force_saver)

#add check pointer to store simulation data for paraview
sim.createCheckPointer (
 CheckPointPrms (
  fileNamePrefix ="snapshot",
  beginTimeStep =0,
  endTimeStep=1000,
  timeStepIncr=4
 )
)
#run the simulation
sim.run()

adminuser@temitope-oladele:~$ export PATH=/usr/local/bin/:$PATH
adminuser@temitope-oladele:~$ export LD_LIBRARY_PATH=/usr/local/lib/:$LD_LIBRARY_PATH
adminuser@temitope-oladele:~$ export LIBRARY_PATH=/usr/local/lib/:$LIBRARY_PATH
adminuser@temitope-oladele:~$ export PYTHONPATH=/usr/local/lib/python2.7/dist-packages/:$PYTHONPATH
adminuser@temitope-oladele:~$ export PYTHONPATH=/usr/local/lib/python2.7/dist-packages/gengeo/:$PYTHONPATH
adminuser@temitope-oladele:~$ mpirun -np 2 esysparticle single_particle.py
CSubLatticeControler::initMPI()
nx,ny,nz: 27 , 27 , 27
InsertGenerator3D::seedParticles
bbx: -0.05 0.275 -0.05 - 0.05 0.375 0.05
total tries: 0
slave started at local/global rank 0 / 1
Traceback (most recent call last):
  File "single_particle.py", line 28, in <module>
    steel_ball.setLinearVelocity(Vec3(0,0,0))
AttributeError: 'SphereVol' object has no attribute 'setLinearVelocity'

kind regards
Tope

Hi Tope,

You need to import the geometry file ("steel_ball.geo") into your simulation using the sim.readGeometry(..) command.

Replace the line containing steel_ball.setLinearVelocity (...) with the following command:
sim.readGeometry("steel_ball.geo")

Cheers,

Dion

Temitope (temitope) said : #20

Hi Dion,

Thanks for the response.

I created a similar particle (particle_1) in gengeo which I intend to break by impact the first one and read it into the simulation script.

1) I do not know how to distinguish the properties (young's modulus, ....) of the two particles in the domain. I was able to create of the steel ball but i'm not sure if it allocated it to the steel ball.

2) I expect the steel ball to fall on particle_1 but seems my model is not running properly. Doesn't the initialisation of the velocity matter anymore after replacing it with sim.readGeometry("steel_ball.geo")?

#Import the appropriate Esys-Particle module
from esys.lsm import*
from esys.lsm.util import Vec3, BoundingBox
from esys.lsm.geometry import*
from gengeo import *
from steel_ball import*

#instantiate a simulation object
#number of processor for simulation
sim = LsmMpi (numWorkerProcesses=1, mpiDimList=[1,1,1])
#initialise the neighbour search algorithm:
sim.initNeighbourSearch(
 particleType="NRotSphere",
 gridSpacing=0.0025, #i.e, 2.5 Xply by max. particle radius which is 0.001m
 verletDist=0.00002 #i.e, 0.2 Xply by min. particle radius which is 0.0001m
)
#setting number of timestep and timestep size
sim.setNumTimeSteps(1000)
sim.setTimeStepSize(0.1)

#Domain of simulation: here is the cube with
domain = BoundingBox(Vec3(0,0,0), Vec3(0.2,0.5,0.2))
sim.setSpatialDomain (domain)

#initialise the velocity of the steel_ball Vx, Vy and Vz (m/s)
###steel_ball.setLinearVelocity(Vec3(0,0,0))

#read the steel_ball geometry
sim.readGeometry("steel_ball.geo")

#read the single particle geometry
sim.readGeometry("particle_1.geo")

#initialise gravity in the domain:
sim.createInteractionGroup(
 GravityPrms(name="Gravitational_acceleration", acceleration=Vec3(0,-9.81,0))
)
#add a horizontal wall to act like the steel rod surface:
sim.createWall(
 name="steel_rod",
 posn=Vec3(0,0,0),
 normal=Vec3(0,1,0)
)
#specify the type of interaction between wall and steel ball
#normalK calculated from steel properties i.e function of young modulus (wikipedia)
sim.createInteractionGroup(
 NRotElasticWallPrms(
  name="elastic_wall",
  wallName="steel_rod",
  normalK= 10000
 )
)

#create rotational elastic-brittle bonds between between steel_ball particles:
pp_bonds = sim.createInteractionGroup (
 BrittleBeamPrms(
  name="Steel_ball_bonds",
  youngsModulus=210000000000.0, #210Gpa for steel
  poissonsRatio=0.29,
  cohesion=100.0, #check this again!
  tanAngle=1.0,
  tag=1
 )
)

#add local viscosity to simulate air resistance
sim.createInteractionGroup(
 LinDampingPrms(
  name="linDamping",
  viscosity=0.1,
  maxIterations=100
 )
)
#create a FieldSaver to store forces on the steelrod:
force_saver=WallVectorFieldSaverPrms(
 wallName=["steel_rod"],
 fieldName="Force",
 fileName="Steel_rod_Force.dat",
 fileFormat="RAW_SERIES",
 beginTimeStep=0,
 endTimeStep=1000,
 timeStepIncr=1
)
sim.createFieldSaver(force_saver)

#add check pointer to store simulation data for paraview
sim.createCheckPointer (
 CheckPointPrms (
  fileNamePrefix ="snapshot",
  beginTimeStep =0,
  endTimeStep=1000,
  timeStepIncr=4
 )
)
#run the simulation
sim.run()

Nothing happened after waiting for about 1 hour for the simulation to run

adminuser@temitope-oladele:~$ export PATH=/usr/local/bin/:$PATH
adminuser@temitope-oladele:~$ export LD_LIBRARY_PATH=/usr/local/lib/:$LD_LIBRARY_PATH
adminuser@temitope-oladele:~$ export LIBRARY_PATH=/usr/local/lib/:$LIBRARY_PATH
adminuser@temitope-oladele:~$ export PYTHONPATH=/usr/local/lib/python2.7/dist-packages/:$PYTHONPATH
adminuser@temitope-oladele:~$ export PYTHONPATH=/usr/local/lib/python2.7/dist-packages/gengeo/:$PYTHONPATH
adminuser@temitope-oladele:~$ mpirun -np 2 esysparticle single_particle.py
CSubLatticeControler::initMPI()
nx,ny,nz: 27 , 27 , 27
InsertGenerator3D::seedParticles
bbx: -0.025 0.3 -0.025 - 0.025 0.35 0.025
total tries: 0
MNTable3D::generateBonds( 0 , 1e-05 , 0 )
slave started at local/global rank 0 / 1
Geometry info read from file is incompatible with previously set geometry (bounding box, circular boundaries) - Model may not run properly!
Geometry info read from file is incompatible with previously set geometry (bounding box, circular boundaries) - Model may not run properly!
constructing FieldMaster for field Force

Kind Regards
Tope

Temitope (temitope) said : #21

Hi Dion and Vince,

Please I am trying to construct my collection of sphere in the simulation script as shown in the code below but I seem not figure out what the challenge is when I ran the code. Please help me with this.

#Import the appropriate Esys-Particle module
from esys.lsm import*
from esys.lsm.util import Vec3, BoundingBox
from esys.lsm.geometry import*
#instantiate a simulation object
#number of processor for simulation
sim = LsmMpi (numWorkerProcesses=1, mpiDimList=[1,1,1])
#initialise the neighbour search algorithm:
sim.initNeighbourSearch(
 particleType="NRotSphere",
 gridSpacing=0.0625, #i.e, 2.5 Xply by max. particle radius
 verletDist=0.005 #i.e, 0.2 Xply by min. particle radius
)
#setting number of timestep and timestep size
sim.setNumTimeSteps(1000)
sim.setTimeStepSize(0.001)

#Domain of simulation: here is the cube with
domain = BoundingBox(Vec3(0,0,0), Vec3(0.2,0.5,0.2))
sim.setSpatialDomain (domain)

#create a steel ball with NRotSphere type, id=0 position at that vector posn, with the specified radius and mass
steel_ball=NRotSphere(id=0, posn=Vec3(0.1,0.325,0.1), radius=0.025, mass=0.51)

#initialise the velocity of the steel_ball Vx, Vy and Vz (m/s)
steel_ball.setLinearVelocity(Vec3(0,0,0))

#add the newly created steel_ball
sim.createParticle(steel_ball)

#read the single particle geometry
#sim.readGeometry("particle_1.geo")

#create a sphere of spherical particles:
geoRandomSphere = RandomSpherePacker (
 minRadius = 0.0001,
 maxRadius = 0.00025,
 cubicPackRadius = 0.0025,
 maxInsertFails = 200000,
 bSphere = BoundingSphere(
  Vec3(0.0000, 0.0000, 0.0000),
  Vec3(0.1000, 0.0025, 0.1000)
 ),
 circDimList = [False, False, False],
 tolerance = 1.0000e-06
)
geoRandomSphere.generate()
geoRandomSphere_particles = geoRandomSphere.getSimpleSphereCollection()

#add the particles to the simulation object:
sim.createParticles(geoRandomSphere_particles)

#bond particles together with bondTag = 1:
sim.createConnections(
 ConnectionFinder(
 maxDist = 0.0005,
 bondTag = 1,
 pList = geoRandomSphere_particles
 )
)

#initialise gravity in the domain:
sim.createInteractionGroup(
 GravityPrms(name="Gravitational_acceleration", acceleration=Vec3(0,-9.81,0))
)

#add a horizontal wall to act like the steel rod surface:
sim.createWall(
 name="steel_rod",
 posn=Vec3(0,0,0),
 normal=Vec3(0,1,0)
)
#specify the type of interaction between wall and steel ball
#normalK calculated from steel properties i.e function of young modulus (wikipedia)
sim.createInteractionGroup(
 NRotElasticWallPrms(
  name="elastic_wall",
  wallName="steel_rod",
  normalK= 1000000
 )
)
#add local viscosity to simulate air resistance
sim.createInteractionGroup(
 LinDampingPrms(
  name="linDamping",
  viscosity=0.1,
  maxIterations=100
 )
)
#create a FieldSaver to store forces on the steelrod:
force_saver=WallVectorFieldSaverPrms(
 wallName=["steel_rod"],
 fieldName="Force",
 fileName="Steel_rod_Force.dat",
 fileFormat="RAW_SERIES",
 beginTimeStep=0,
 endTimeStep=1000,
 timeStepIncr=1
)
sim.createFieldSaver(force_saver)

#add check pointer to store simulation data for paraview
sim.createCheckPointer (
 CheckPointPrms (
  fileNamePrefix ="snapshot",
  beginTimeStep =0,
  endTimeStep=1000,
  timeStepIncr=4
 )
)
#run the simulation
sim.run()

adminuser@temitope-oladele:~/PhD_ESyScode/Testrun$ mpirun -np 2 esysparticle single_breakage.py
CSubLatticeControler::initMPI()
slave started at local/global rank 0 / 1
Traceback (most recent call last):
  File "single_breakage.py", line 48, in <module>
    Vec3(0.1000, 0.0025, 0.1000)
Boost.Python.ArgumentError: Python argument types in
    BoundingSphere.__init__(BoundingSphere, Vec3, Vec3)
did not match C++ signature:
    __init__(_object*, esys::lsm::Vec3Py centrePt=<esys.lsm.util.FoundationPy.Vec3 object at 0x7fe85d9c2ed0>, double radius=0.0)
    __init__(_object*, boost::python::api::object seq, double radius)
    __init__(_object*, esys::lsm::BoundingSpherePy bSphere)
    __init__(_object*)

Kind Regards
Tope

Temitope (temitope) said : #22

i figured it out