simulation.PyGranSim.models

A module that provides contact models for running numerical experiments or DEM simulation

Created on July 1, 2016

Author: Andrew Abi-Mansour

This is the:

██████╗ ██╗   ██╗ ██████╗ ██████╗  █████╗ ███╗   ██╗
██╔══██╗╚██╗ ██╔╝██╔════╝ ██╔══██╗██╔══██╗████╗  ██║
██████╔╝ ╚████╔╝ ██║  ███╗██████╔╝███████║██╔██╗ ██║
██╔═══╝   ╚██╔╝  ██║   ██║██╔══██╗██╔══██║██║╚██╗██║
██║        ██║   ╚██████╔╝██║  ██║██║  ██║██║ ╚████║
╚═╝        ╚═╝    ╚═════╝ ╚═╝  ╚═╝╚═╝  ╚═╝╚═╝  ╚═══╝

DEM simulation and analysis toolkit http://www.pygran.org, support@pygran.org

Core developer and main author: Andrew Abi-Mansour, andrew.abi.mansour@pygran.org

PyGran is open-source, distributed under the terms of the GNU Public License, version 2 or later. It is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have received a copy of the GNU General Public License along with PyGran. If not, see http://www.gnu.org/licenses . See also top-level README and LICENSE files.

Warning

Only particle-wall collisions are supported

Todo

Support 2-particle analysis by replacing mass, radius, etc. with reduced mass, radius, etc. i.e. 1/m_{ij} = 1/m_i + 1/m_j

Functions

template_tablet(nspheres, radius, length)

This function creates a multi-sphere tablet (cylinder) of radius “radius” and height “length” constituting “nspheres” spheres.

Classes

HertzMindlin(**params)

A class that implements the linear spring model for granular materials

Model(**params)

This class implements a contact model.

SpringDashpot(**params)

A class that implements the linear spring model for granular materials

ThorntonNing(**params)

A basic class that implements the Thornton elasto-plastic model based on [TN98].

simulation.PyGranSim.models.template_tablet(nspheres, radius, length)[source]

This function creates a multi-sphere tablet (cylinder) of radius “radius” and height “length” constituting “nspheres” spheres. This is used for multisphere simulations.

Todo

Move this function elsewhere.

Parameters
  • nspheres (int) – number of spheres each tablet consists of

  • radius (float) – particle radius

  • length (float) – length of the tablet

Returns

DEM representation of the tablet

Return type

tuple

simulation.PyGranSim.models.template_multisphere(func)[source]

This function creates a generic multi-sphere particle based on a user-supplied function.

Parameters

func (function) – returns a list of tuples: [(x1,y1,z1,radius1), …]

class simulation.PyGranSim.models.Model(**params)[source]

Bases: object

This class implements a contact model. It can be used to run a DEM simulation (e.g. with LIGGGHTS) or numerical experiments for simulating particle-wall collisions.

Parameters
  • mesh (dict) – a dictionary that defines all meshes and their properties

  • box (tuple) – simulation box size boundaries (xmin, xmax, ymin, ymax, zmin, zmax)

  • restart (tuple) – specify restart options via (freq=int, dirname=str, filename=str, restart=bool, frame=int)

  • nSim (int) – number of concurrent simulations to run (default 1)

  • units (str) – unit system (default ‘si’). See here for available options.

  • comm (MPI Intracomm) – MPI communicator (default COMM_WORLD)

  • rank (int) – processor rank

  • debug (bool) – set debug mode on (default False)

  • engine (str) – DEM engine (default ‘engine_liggghts’)

  • species (tuple) – defines the number and properties of all species

Todo

Support particle-particle collisions

setupProps()[source]

Creates class attributes for all material properties

contactTime()[source]

Computes the characteristic collision time assuming for a spring dashpot model

Returns

the contact time

Return type

float

displacement(dt=None)[source]

Generator that computes (iteratively) the overlap distance (displacement) as a function of time

Parameters

dt (float) – timestep used by the ODE integrator, by default it is 0.1% of contact duration

Returns

time, displacement, and force as numpy arrays

Return type

tuple

Todo

Enable the user control the timestep resolution (number of timesteps).

contactRadius

Returns the contact radius.

Returns

contact radius

Return type

numpy array

Note

This function is not useful since self._contRadius is not updated anywhere.

_contactRadius(delta, radius)[source]

Internal function that computes the contact radius for Hertzian and JKR models

Parameters
  • delta (float) – normal displacement

  • radius (float) – effective radius

Returns

contact radius

Return type

float

numericalForce(time, delta)[source]

Returns the force used for numerical solvers

class simulation.PyGranSim.models.SpringDashpot(**params)[source]

Bases: simulation.PyGranSim.models.Model

A class that implements the linear spring model for granular materials

Parameters
  • material (dict) – a python dictionary that specifies the material properties

  • limitForce (bool) – turns on a limitation on the force, preventing it from becoming attractive at the end of a contact

springStiff(delta=None)[source]

Computes the spring constant (k_n) for the normal force F_n = k_n \delta_n.

Parameters

delta (float) – normal displacement (\delta_n)

Returns

spring stiffness (k_n)

Return type

float

dissCoef(delta=None)[source]

Computes the normal dissipative coefficient (c_n) for the dissipative force F_d = - c_n \dot{\delta_n}

Returns

dissipative coefficient (c_n)

Return type

float

dissForce(delta, deltav)[source]

Returns the dissipative (viscous) force: F_d = - c_n \dot{\delta_n} where c_n is the dissipative coefficient.

Parameters
  • delta (float) – normal displacement

  • deltav (float) – time derivative of the normal displacement

Returns

dissipative force F_d

Return type

float

displacementAnalytical(dt=None)[source]

Computes the displacement based on the analytical solution for a spring-dashpot model.

Parameters

dt (float) – timestep. By default, the timestep is 1% of the contact time.

Returns

time, displacement, and force arrays of size N, N by 2, and N, respectively, where N is the total number of steps taken by the integrator. The displacement array stores the overlap in its 1st column and the overlap velocity (time derivative) in its 2nd column.

Return type

tuple

Todo

Take cohesion (simplified JKR) into account

elasticForce(delta)[source]

Returns the elastic force based on Hooke’s law: F_n = k_n \delta_n

Parameters

delta (float) – normal displacement

Returns

elastic contact force

Return type

float

cohesiveForce(delta)[source]

Returns the cohesive force F_c (for the SJKR model)

Parameters

delta (float) – normal displacement

Returns

cohesvie force

Return type

float

normalForce(delta, deltav)[source]

Returns the total normal force F_n + F_d + F_c

Parameters
  • delta (float) – normal displacement

  • deltav (float) – time derivative of the normal displacement

Returns

total normal force

Return type

float

tangForce(delta, deltav)[source]

Returns the total tangential force based on the Coulomb model

_contactRadius(delta, radius)

Internal function that computes the contact radius for Hertzian and JKR models

Parameters
  • delta (float) – normal displacement

  • radius (float) – effective radius

Returns

contact radius

Return type

float

contactRadius

Returns the contact radius.

Returns

contact radius

Return type

numpy array

Note

This function is not useful since self._contRadius is not updated anywhere.

contactTime()

Computes the characteristic collision time assuming for a spring dashpot model

Returns

the contact time

Return type

float

displacement(dt=None)

Generator that computes (iteratively) the overlap distance (displacement) as a function of time

Parameters

dt (float) – timestep used by the ODE integrator, by default it is 0.1% of contact duration

Returns

time, displacement, and force as numpy arrays

Return type

tuple

Todo

Enable the user control the timestep resolution (number of timesteps).

numericalForce(time, delta)

Returns the force used for numerical solvers

setupProps()

Creates class attributes for all material properties

class simulation.PyGranSim.models.HertzMindlin(**params)[source]

Bases: simulation.PyGranSim.models.Model

A class that implements the linear spring model for granular materials

springStiff(delta)[source]

Computes the spring constant k_n for F_n = k_n delta_n^{3/2}

Parameters

delta (float) – normal displacement

Returns

stiffness (k_n)

Return type

float

elasticForce(delta)[source]

Returns the Hertzian elastic force F_n = k_n \delta_n^{3/2}

Parameters

delta (float) – normal displacement

Returns

elastic normal force

Return type

float

dissCoef(delta)[source]

Returns the dissipative force coefficient

Parameters

delta (float) – normal displacement

Returns

coefficient for dissipative (viscous) normal force

Return type

float

dissForce(delta, deltav)[source]

Returns the dissipative force

Parameters
  • delta (float) – normal displacement

  • deltav (float) – time derivative of the normal displacement

Returns

dissipative force

Return type

float

normalForce(delta, deltav)[source]

Returns the total normal force F_n + F_d + F_c

Parameters
  • delta (float) – normal displacement

  • deltav (float) – time derivative of the normal displacement

Returns

total normal force

Return type

float

cohesiveForce(delta)[source]

Returns the cohesive force for the simplified JKR model

Parameters

delta (float) – normal displacement

Returns

cohesive force

Return type

float

_contactRadius(delta, radius)

Internal function that computes the contact radius for Hertzian and JKR models

Parameters
  • delta (float) – normal displacement

  • radius (float) – effective radius

Returns

contact radius

Return type

float

contactRadius

Returns the contact radius.

Returns

contact radius

Return type

numpy array

Note

This function is not useful since self._contRadius is not updated anywhere.

contactTime()

Computes the characteristic collision time assuming for a spring dashpot model

Returns

the contact time

Return type

float

displacement(dt=None)

Generator that computes (iteratively) the overlap distance (displacement) as a function of time

Parameters

dt (float) – timestep used by the ODE integrator, by default it is 0.1% of contact duration

Returns

time, displacement, and force as numpy arrays

Return type

tuple

Todo

Enable the user control the timestep resolution (number of timesteps).

numericalForce(time, delta)

Returns the force used for numerical solvers

setupProps()

Creates class attributes for all material properties

class simulation.PyGranSim.models.ThorntonNing(**params)[source]

Bases: simulation.PyGranSim.models.Model

A basic class that implements the Thornton elasto-plastic model based on [TN98].

computeYieldRadius()[source]

Computes the contact radius at the yield point based on Eq. (65) in [TN98]

return: yielding contact radius rtype: float

springStiff(delta)[source]

Computes the spring constant k_n for F_n = k_n delta_n^{3/2}

Parameters

delta (float) – normal displacement

Returns

stiffness (k_n)

Return type

float

elasticForce(delta)[source]

Returns the Hertzian-like elastic force

Parameters

delta (float) – normal displacement

Returns

elastic normal force

Return type

float

dissForce(delta, deltav=None)[source]

Computes the piece-wise defined elastic force

Parameters
  • delta (float) – normal displacement

  • deltav (float) – time derivative of the normal displacement

Returns

dissipative force

Return type

float

_contactRadius(delta, radius)

Internal function that computes the contact radius for Hertzian and JKR models

Parameters
  • delta (float) – normal displacement

  • radius (float) – effective radius

Returns

contact radius

Return type

float

contactRadius

Returns the contact radius.

Returns

contact radius

Return type

numpy array

Note

This function is not useful since self._contRadius is not updated anywhere.

contactTime()

Computes the characteristic collision time assuming for a spring dashpot model

Returns

the contact time

Return type

float

displacement(dt=None)

Generator that computes (iteratively) the overlap distance (displacement) as a function of time

Parameters

dt (float) – timestep used by the ODE integrator, by default it is 0.1% of contact duration

Returns

time, displacement, and force as numpy arrays

Return type

tuple

Todo

Enable the user control the timestep resolution (number of timesteps).

normalForce(delta, deltav)[source]

Returns the total normal force

Parameters
  • delta (float) – normal displacement

  • deltav (float) – time derivative of the normal displacement

Returns

total normal force

Return type

float

numericalForce(time, delta)

Returns the force used for numerical solvers

setupProps()

Creates class attributes for all material properties

cohesiveForce(delta=None)[source]

Returns the JKR cohesive force

Parameters

delta (None) – useless variable kept here for syntax consistency

Returns

cohesive force

Return type

float

yieldVel

Returns the minimum velocity required for a colliding particle to undergo plastic deformation

Returns

yielding velocity

Return type

float