units system definition

Hi Diego,

Thanks for your important question as it raises the somewhat confusing issue of scaling between model units and real units in ESyS-Particle simulations. In this answer I will first give you a short answer followed by a quite long, detailed answer explaining how to scale between model units and real units in simulations.

1) Short Answer:

You are correct to surmise that the problem lies with particles having incorrect masses. Unless you explicitly set the mass or density of particles in ESyS-Particle, the software will assume particles have unit density. A particle with a radius = 1 will be assigned a mass = 4/3 pi. According to the formula (Kwall > m.g/r) you need to set the stiffness Kwall > 9810 if you choose to set gravity = 9810 and measure lengths in mm. If you increase Kwall to about 10000, your particles should stop falling through the floor. You may also need to reduce your timestep increment to less than 0.001 to ensure numerical stability of the solution.

Although this is a very simple answer to your specific question, there are a number of important issues when choosing values for model parameters. In the second answer below, I will use your example to try to explain some of the issues with scaling between model units and real units. This is an important topic that all ESyS-Particle users need to keep in mind.

2) Long Answer:

The mathematical implementation of ESyS-Particle is non-dimensional – we make no assumptions on the units of measure used in the model. Consequently the user has the freedom to use the system of units that best suits the problem they wish to solve. If you decide that a particle of unit radius in the model will represent an object in the real world that is 1mm in radius then you have chosen to measure lengths in millimeters. Another way of saying this is that the unit of length in the model is equivalent to 1mm. When you analyse results from the model (such as displacement) these will be reported in mm. If your particle moves a distance 2.35 in a simulation, this is equivalent to 2.35mm in the real world.

The next model parameter you set in your question was the acceleration due to gravity. In this case you chose to set gravity = 9810 in the model. Since we have decided to measure lengths in the model in mm, by setting gravity = 9810 you have implicitly decided that time will be measured in seconds. You are correct to suppose that the unit of time in the model is equivalent to 1 second but only because you chose that to be the case when you set the value for gravity to be 9810. You could have equally well set gravity = 9.81 but by doing so you would have changed your unit of time to be equivalent to 0.0316sec i.e. each tick of the clock in the simulation is equivalent to 0.0316sec (a bit inconvenient!)

The third model parameter you set was the wall stiffness, Kwall = 10 and it would appear from your question that you wish this to be equivalent to 10 N/mm. In other orders you have decided that the unit of force in the model will be equivalent to Newtons. This is where things can get confusing. By selecting millimeters, seconds and Newtons as your units of measure for length, time and force respectively, you have implicitly set the units of measure of all other physical quantities, including (most importantly) the unit of mass in the model. In this case, a particle with unit mass in the model is equivalent to an object weighing 1000kg! It should now be obvious why your particles just fall straight through the wall...they are far too heavy! As I said above, unless you explicitly set the mass (or density) of a particle in ESyS-Particle, the software will assume particles have unit density (so a particle of radius = 1 will have a mass = 4/3 pi).

For your simulations, it would probably be easier to choose an appropriate unit of mass and rescale the wall stiffness according to the formula Kwall > m.g/r. Suppose you decide that a particle of unit mass is equivalent to an object weighing 1kg in the real world. Then, according to the formula, you should set Kwall >= 9810. However, by choosing kilograms, millimeters and seconds as measures of mass, length and time, the simulation will now report forces in milliNewtons rather than Newtons. Whenever you take a measurement of force from your simulations, you will need to divide by 1000 to convert model forces to Newtons.

Likewise, all other physical quantities will be reported in units of kilograms, millimeters and seconds so you will need to be careful to convert these as well. For example, suppose you measure the velocity of a particle in a simulation to be 37.6. This will be equivalent to 37.6 mm/s so if you wish to plot the velocity in m/s you need to remember to divide by 1000. As a second example, suppose you measure the stress experienced by a particle and it is reported as 234.5. This will be equivalent to 234.5 mN/mm^2 = 234500 Pa or if you prefer, 234.5 kPa. In other words, the unit of stress in the model is equivalent to 1kPa.

You are free to choose any system of units you wish when setting up a simulation but you must be careful to ensure you set values for all model parameters such that they are all consistent with your chosen system of units. In all cases, you need to select three physical quantities to define your fundamental units. These three quantities do not need to be mass, length and time but the chosen three quantities must include each of mass, length and time in their units, in some combination. Hence you can choose length, force and velocity to define your fundamental units but not time, velocity and acceleration (because the unit of mass will be undefined in this later case).

I hope this answer is helpful both for yourself Diego and for other ESyS-Particle users. Have fun with your simulations!

Cheers,

Dion.

Thanks Dion.
I think I almost understood the units issue
I will use meters, seconds and kilograms, that is Length, Time and Mass based system of units. The radius will be 1e-3 meters, g= 9.81. Also the Kn is 245 N/m.
Here is the difficult point. As I understood, the mass is then 4/3*pi*(1e-3)^3 = 4.19e-09 (the density is taken as 1 kg/m^3 ? or as 1000 kg/m^3?) if this is so, then m.g/r = 4.11e-5 << 245. Even if density is 1000 kg/m^3 is 4.11e-2 << 245.
The timestep I take as dt=1e-5.

The problem I run is

********************************************************************************
********************************************************************************
********************************************************************************

# -*- coding: utf-8 -*-
from esys.lsm import *
from esys.lsm.util import Vec3, BoundingBox
from esys.lsm.geometry import HexagBlock, SimpleSphereNeighbours
import cPickle as pickle
import math

#
# Generate the rectangular block of close-packed particles
#

radio = 1e-3

blk = HexagBlock(dimCount = [1,1,1], radius = radio)
blkBBox = blk.getBBox()
# Calculate the axis aligned bounding-box for use as spatial domain

minX= blkBBox.getMinPt()[0]
maxX= blkBBox.getMaxPt()[0]
minY = blkBBox.getMinPt()[1]
maxY = blkBBox.getMaxPt()[1]
minZ = blkBBox.getMinPt()[2]
maxZ = blkBBox.getMaxPt()[2]

print minX,maxX,minY,maxY,minZ,maxZ
blah = raw_input("press enter to continue")

domain = BoundingBox(Vec3(minX,minY,minZ), Vec3(maxX,100*maxY,maxZ))

#domain = BoundingBox(Vec3(-5,-5,0),Vec3(5,20,0))

# Initialise the ESyS_Particle simulation object
lsm = LsmMpi(numWorkerProcesses=1, mpiDimList=[1,1,1])
lsm.initVerletModel(\
  particleType="NRotSphere",
  gridSpacing=2.2*radio,
  verletDist=0.2*radio)

# Initialise spatial domain of simulation.
lsm.setSpatialDomain(domain)

# Make ESyS_Particle do two-dimensional calculations.
lsm.force2dComputations(True)

# Create Non rotational spherical particles within the
# model from the collection of SimpleSphere objects.
lsm.createParticles(blk)

lsm.createInteractionGroup(
   NRotElasticPrms(
      name = "elasticSpheres",
      normalK = 245
   )
)

# Create viscosity damping to prevent build up of kinetic energy
#lsm.createInteractionGroup(
# LinDampingPrms(
# name="linDamping",
# viscosity=0.05,
# maxIterations=100
# )
#)

lsm.createWall(
    name = "bottomWall",
    posn = blkBBox.getMinPt()-Vec3(0,1e-3,0),
    normal = Vec3(0, 1, 0)
)
lsm.createWall(
    name = "leftWall",
    posn = blkBBox.getMinPt(),
    normal = Vec3(1, 0, 0)
)
lsm.createWall(
    name = "rightWall",
    posn = blkBBox.getMaxPt(),
    normal = Vec3(-1, 0, 0)
)

# Create contact interaction between bottom confining wall and particles
lsm.createInteractionGroup(
    NRotElasticWallPrms(
        name = "bottomWallGroup",
        wallName = "bottomWall",
        normalK = 245
    )
)
# Create contact interaction between left confining wall and particles
lsm.createInteractionGroup(
    NRotElasticWallPrms(
        name = "leftWallGroup",
        wallName = "leftWall",
        normalK = 245
    )
)
# Create contact interaction between right confining wall and particles
lsm.createInteractionGroup(
    NRotElasticWallPrms(
        name = "rightWallGrp",
        wallName = "rightWall",
        normalK = 245
    )
)

lsm.createInteractionGroup(
   GravityPrms(name="earth-gravity", acceleration=Vec3(0,-9.810,0))
)

class DataSaverRunnable(Runnable):
    #def __init__(self, lsm, timeStepIncr,WBR, fileNameTemplate="fracData_%05d.pkl"):
    def __init__(self, lsm, timeStepIncr, fileNameTemplate="fracData_%05d.pkl"):
        Runnable.__init__(self) # NB: must call Runnable init method in child
        self.lsm = lsm
        self.timeStepIncr = timeStepIncr
        self.fileNameTemplate = fileNameTemplate
        self.saveNum = 0
 #self.WBR=WBR
    def run(self):
        """
        Saves current state of all model particles as well as
        the non-broken bond id pairs at every Nth
        time-step (where N = self.timeStepIncr).
        """
        print "Time-step = %5d of %5d" % (self.lsm.getTimeStep(), self.lsm.getNumTimeSteps())
        if (self.lsm.getTimeStep()==1 or (self.lsm.getTimeStep() % self.timeStepIncr) == 0):
            dataList = []
            dataList.append(self.lsm.getParticleList())
     print "Y = ", dataList[0][0].getPosn()[1]
     #posW=self.WBR.getDispl()
     #dataList.append(posW)
            fileName = self.fileNameTemplate % self.saveNum
            print "Writing data to file:", fileName
            print "num particles = %5d" % len(dataList[0])
     #print "Wall position: ", posW
            f = file(fileName, "w")
            pickle.dump(obj=dataList, file=f, protocol=pickle.HIGHEST_PROTOCOL)
            f.close()
            self.saveNum += 1

# Add a data saver to the model, save data every 50 time-steps.
#dataSaverRunnable = DataSaverRunnable(lsm=lsm, timeStepIncr=1, WBR=wallVibratingRunnable)
dataSaverRunnable = DataSaverRunnable(lsm=lsm, timeStepIncr=10)
lsm.addPreTimeStepRunnable(dataSaverRunnable)

# Now run the simulation
dt = 1e-5
maxTime = dt*10000
lsm.setTimeStepSize(dt)
lsm.setNumTimeSteps(int(maxTime/dt))
lsm.run()

# Remove references to the lsm object from the runnables
# for a clean MPI-program exit.
dataSaverRunnable.lsm = None
#wallVibratingRunnable.lsm = None

********************************************************************************
********************************************************************************
********************************************************************************

The results I have are very strange. The particle remains still for sometime, and then suddenly goes up at high velocity. Probably the particle has gone down a very large gap in the timestep, so the restoring force is so big that the particle is thrown away.

The question is:

What I'm doing wrong?
I expect to have the particle resting on the floor.

Thanks in advance and best regards,

Thanks Dion Weatherley, that solved my question.