Documentation of contact model chain elements

AssembleContactMatrix

Sets actual contact matrix element of the contact state to -M where M is the net interaction tensor that is computed differently for underdamped and overdamped equations of motion.

Implicit solvers that solve some system A x = b use this element as an i-j element of A, but must also sum M to the matrix diagonal to account for the symmetric interaction (see SumToDiagonalCommand).

AssembleContactMatrix_ComposedParticle

Sets the element of the contact matrix (contained in the contact state) to the appropriate ‘sub-matrix’ for any combination of (complex) composite particles. A composite particle is a particle for which the degrees of freedom do not overlap with the contact geometry. As such, the contact matrix element for the contact primitive is only ‘virtual’ and an actual friction element for the degrees of freedom needs to be properly set. To assemble a contact matrix element M to the degrees of freedom, an additional DOF contact chain element needs to be present before this assembly. It knows how many degrees of freedom each primitive should assemble to, and it computes a vector with ‘weights’ to indicate how much of this friction should be distributed to each subsequent dof. For dof i:

M_{i12} = w_{i1}*w_{i2}\,M.

This formulation ensures that the sum of all contact matrix elements remains conserved as M for both p1 and p2. For example, for deformable-cell-like models, these weights will be identical to the force weights computed using DistributeForcesAndMomentsToTriangle_1_Nodes and DistributeForcesAndMomentsToTriangle_2_Nodes.

AssembleForces

Adds up the total contact force to the force degrees of freedom of particles p1 and p2. If the total force is F, p1 receives -F and p2 receives F. If the rejectContact flag has been set to non-zero, no forces are assembled.

Note

An atomic add is performed to enable parallel compatibility.

AssembleForces_1

Adds up the total contact force to the force degrees of freedom of particle p1. If the total force is F, p1 receives -F. If the rejectContact flag has been set to non-zero, no forces are assembled.

Note

An atomic add is performed to enable parallel compatibility.

AssembleForces_2

Adds up the total contact force to the force degrees of freedom of particle p2. If the total force is F, p2 receives F. If the rejectContact flag has been set to non-zero, no forces are assembled.

Note

An atomic add is performed to enable parallel compatibility.

AssembleForces_DeformableCapsuleWithTorque

Ensures distribution of the contact force to deformable cylinder with rotational degrees of freedom. In this case the contact force is simply added to the total force applied on the cylinder.The actual distribution (de-facto integration) of forces and torques toward the nodes velocities is achieved in the overdamped case by function ComputeVelocitiesFromSegmentForcesAndTorque.

Note

An atomic add is performed to enable parallelism.

AssembleForces_DeformableCapsuleWithTorque_1

Ensures distribution of the contact force to deformable cylinder with rotational degrees of freedom. In this case the contact force is simply added to the total force applied on the cylinder.The actual distribution (de-facto integration) of forces and torques toward the nodes velocities is achieved in the overdamped case by function ComputeVelocitiesFromSegmentForcesAndTorque.

Note

An atomic add is performed to enable parallelism.

AssembleForces_DeformableCapsuleWithTorque_2

Ensures distribution of the contact force to deformable cylinder with rotational degrees of freedom. In this case the contact force is simply added to the total force applied on the cylinder.The actual distribution (de-facto integration) of forces and torques toward the nodes velocities is achieved in the overdamped case by function ComputeVelocitiesFromSegmentForcesAndTorque.

Note

An atomic add is performed to enable parallelism.

AssembleForces_DeformableCylinder_1

Ensures the distribution of the contact force to end points’ forces. We use a lever mechanism to distribute the contact force \vec{F} and contact moment \vec{M} to the end points,

\vec{F}_1 = \left[ -(I-\hat{t}\hat{t})(\vec{x}_2 - \vec{x}_C) (\hat{t}\cdot \vec{F}) +                 \hat{t} \times \vec{M} +                 (\vec{x}_2 - \vec{x}_C) \cdot \vec{F}\right]/L,

\vec{F}_2 = \left[ (I-\hat{t}\hat{t})(\vec{x}_2 - \vec{x}_C) (\hat{t}\cdot \vec{F}) -                 \hat{t} \times \vec{M} -                 (\vec{x}_2 - \vec{x}_C) \cdot \vec{F}\right]/L,

where \vec{F}_{12} represent the forces on the end points, L is the cylinder’s length, \vec{x}_\text{C} is the contact point, and \hat{t} is the unit vector along the cylinder`s axis.This assembly ensures that:

  • The sum of forces on the capsule end point is equal to the total contact force.
  • The sum of moments generated by the end point forces w.r.t. the contact point is equal to contact moment (couple).
  • The axial forces are inversely propertional to the distance of each end point to the contat point. This assumption is consistent with the assumption of linear elasticity in the deformable segment, i.e. for equal deformations the force changes inversely with the segment length.

Note

An atomic add is performed to enable parallelism.

AssembleForces_DeformableCylinder_2

Ensures the distribution of the contact force to end points’ forces. We use a lever mechanism to distribute the contact force \vec{F} and contact moment \vec{M} to the end points,

\vec{F}_2 = \left[ (I-\hat{t}\hat{t})(\vec{x}_2 - \vec{x}_C) (\hat{t}\cdot \vec{F}) -                 \hat{t} \times \vec{M} -                 (\vec{x}_2 - \vec{x}_C) \cdot \vec{F}\right]/L,

\vec{F}_1 = \left[ -(I-\hat{t}\hat{t})(\vec{x}_2 - \vec{x}_C) (\hat{t}\cdot \vec{F}) +                 \hat{t} \times \vec{M} +                 (\vec{x}_2 - \vec{x}_C) \cdot \vec{F}\right]/L,

where \vec{F}_{12} represent the forces on the end points, L is the cylinder’s length, \vec{x}_\text{C} is the contact point, and \hat{t} is the unit vector along the cylinder`s axis.This assembly ensures that:

  • The sum of forces on the capsule end point is equal to the total contact force.
  • The sum of moments generated by the end point forces w.r.t. the contact point is equal to contact moment (couple).
  • The axial forces are inversely propertional to the distance of each end point to the contat point. This assumption is consistent with the assumption of linear elasticity in the deformable segment, i.e. for equal deformations the force changes inversely with the segment length.

Note

An atomic add is performed to enable parallelism.

AssembleForces_DeformableTriangle_1

Assembles a contact force to the triangle’s DOF nodes of p1, by using distribution vector ntnf from a previous chain element. Additionally, any pure tangential forces are assumed to contribute equally and one third of these is added to each node.

AssembleForces_DeformableTriangle_2

Assembles a contact force to the triangle’s DOF nodes of p2, by using distribution vector ntnf from a previous chain element. Additionally, any pure tangential forces are assumed to contribute equally and one third of these is added to each node.

AssembleForces_Primitives_1

Assembles the contact force to a force array Fprim1 that exists for the contact primitives (the actual particles contained in the particle container passed to the contact model). Typically these forces are not used in the equation of motion but for postprocessing reasons, e.g. to analyze local forces that act on a certain geometry.

AssembleForces_Primitives_2

Assembles the contact force to a force array Fprim2 that exists for the contact primitives (the actual particles contained in the particle container passed to the contact model). Typically these forces are not used in the equation of motion but for postprocessing reasons, e.g. to analyze local forces that act on a certain geometry.

AssembleMoments

Adds up the total contact model moment to the moment degrees of freedom of particles p1 and p2. The following moment is added for p1:

\Delta \vec{M}_1 = - ( \vec{x}_{cp} - \vec{x}_1 ) \times \vec{F}_c - \vec{M}_c

and for particle p2:

\Delta \vec{M}_2 =   ( \vec{x}_{cp} - \vec{x}_2 ) \times \vec{F}_c + \vec{M}_c

where \vec{x}_{cp} is the contact point and F_c and M_c are the contact model force and moment. If the rejectContact flag has been set to non-zero, no forces are assembled.

Note

An atomic add is performed to enable parallel compatibility.

AssembleMoments_1

Adds up the total contact model moment to the moment degrees of freedom of particle p1. The following moment is added for p1:

\Delta \vec{M}_1 = - ( \vec{x}_{cp} - \vec{x}_1 ) \times \vec{F}_c - \vec{M}_c

where \vec{x}_{cp} is the contact point and F_c and M_c are the contact model force and moment. If the rejectContact flag has been set to non-zero, no forces are assembled.

Note

An atomic add is performed to enable parallel compatibility.

AssembleMoments_2

Adds up the total contact model moment to the moment degrees of freedom of particle p2. The following moment is added for p2:

\Delta \vec{M}_2 = ( \vec{x}_{cp} - \vec{x}_2 ) \times \vec{F}_c + \vec{M}_c

where \vec{x}_{cp} is the contact point and F_c and M_c are the contact model force and moment. If the rejectContact flag has been set to non-zero, no forces are assembled.

Note

An atomic add is performed to enable parallel compatibility.

AssembleMoments_cps_1

Assemble the moment vector with respect to the contact point sphere M_cps, as computed by any contact model that includes RoundedTriangles, to the moment degrees of freedom of particle 1. If the rejectContact flag has been set to non-zero, no moments are assembled.

Note

An atomic add is performed to enable parallel compatibility.

AssembleMoments_cps_2

Assemble the moment vector with respect to the contact point sphere M_cps, as computed by any contact model that includes RoundedTriangles, to the moment degrees of freedom of particle 2. If the rejectContact flag has been set to non-zero, no moments are assembled.

Note

An atomic add is performed to enable parallel compatibility.

AssembleTorque_DeformableCapsule_1

Assemble the effective torque \vec{T} based on the contact point position \vec{x} and the contact force \vec{F}. The exact contribution is given by

\Delta \vec{T} = \left( \vec{x} - \vec{x}_\text{C} \right) \cdot \hat{t} \, \vec{F},

where \vec{x}_\text{C} is the position of the center of mass of the cylinder and \hat{t} is the unit vector parallel to the cylinder axis.

The actual distribution (de-facto integration) of forces and torques toward the nodes velocities is achieved in the overdamped case by function ComputeVelocitiesFromSegmentForcesAndTorque.

See article for details of the model.

Note

An atomic add is performed to enable parallelism.

AssembleTorque_DeformableCapsule_2

Assemble the effective torque \vec{T} based on the contact point position \vec{x} and the contact force \vec{F}. The exact contribution is given by

\Delta \vec{T} = \left( \vec{x} - \vec{x}_\text{C} \right) \cdot \hat{t} \, \vec{F},

where \vec{x}_\text{C} is the position of the center of mass of the cylinder and \hat{t} is the unit vector parallel to the cylinder axis.

The actual distribution (de-facto integration) of forces and torques toward the nodes velocities is achieved in the overdamped case by function ComputeVelocitiesFromSegmentForcesAndTorque.

See article for details of the model.

Note

An atomic add is performed to enable parallelism.

Assemble_ContactArea_1

Adds the contact area computed in this chain to the contact_area array of each respective particle.

Assemble_ContactArea_2

Adds the contact area computed in this chain to the contact_area array of each respective particle.

Assemble_SpringTension

For a spring-like connection within a triangulated mesh, determines the two-dimensional tension (in N/m) based on the interaction force based on the interaction force and the mesh dimensions. For node i, the contribution to the tension from the spring force is approximated as:

T_i = - \frac{ f_{ij} }{ N_i \tan( \pi / N_i ) d_{ij} }

where N_i is the number of triangles connected to node i.

Note

This chain element assumes that an array with Scalar values for each ‘node’ is present and re-set to zero at each timestep of a simulation.

Assemble_SpringTension1D

For a spring-like connection within a one-dimensional string, determines the tension (in N) based on the interaction force based on the interaction force. Only half of the force is added because the tensile force is added on both sides of the node:

T_i = - \frac{f_{ij}}{2}

Note

This chain element assumes that an array with Scalar values for each ‘node’ is present and re-set to zero at each timestep of a simulation.

Assemble_VirialStress

Assembles the ‘virial stress’ in a ‘stress’ array with 3x3 symmetric matrix elements that exist for particle containers pc1 and pc2. We compute the stress tensor:

\mathbf{\sigma}_1 = \frac{ - 3 }{ 4 \pi r_1^2 }\vec{F}_{12} \otimes \vec{\hat{n}}_{12}

and

\mathbf{\sigma}_1 = \frac{ - 3 }{ 4 \pi r_2^2 }\vec{F}_{12} \otimes \vec{\hat{n}}_{12}

Assemble_VirialStress_2

Assembles the ‘virial stress’ in a ‘stress’ array with 3x3 symmetric matrix elements that exist for particle container pc2 onl. We compute the stress tensor:

\mathbf{\sigma}_1 = \frac{ - 3 }{ 4 \pi r_2^2 }\vec{F}_{12} \otimes \vec{\hat{n}}_{12}

Capsule_Capsule

Contact geometry resolution for interaction between two capsules. A capsule is defined as a cylinder with spheres of the same radius connected to both end points. For capsules, the contact geometry can be easily resolved by determining the closest points between the two cylinders’ line segments. The following steps are performed:

  1. Determine the closest distance d_{12} between both line segments, and compute the parameterized coordinates of the respective closest points.
  2. Compute the overlap from this distance as \delta = r_1 + r_2 - d_{12}.
  3. Compute the normal unit vector from the vector connecting the closest points.
  4. Compute the contact point at the point on the connecting line where the penetration is maximal.
  5. We set rejectContact = (overlap < 0).

Note

The effective contact radius is determined as: \hat{R} = \frac{R1\,R_2}{R_1 + R_2}.

Capsule_Sphere

Contact geometry resolution between a capsule and a sphere. A capsule is a cylinder with spheres connected to the end points with the same radius as the cylinder. This resolver computes the parameterized coordinate:

c_{12} = (\vec{x}_s-\vec{x}_{c1})\cdot(\vec{x}_{c2}-\vec{x}_{c1}),

If this coordinate is below zero or above the cylinder’s length, we just resolve contact with one of the end spheres. Otherwise, we resolve contact with the cylinder hull.

ComputeContactMatrixElementOverdamped

Sets the contact matrix element between particle p1 and p2 for an overdamped system in which a linear system \Gamma \cdot \vec{v} = \vec{F} is solved.

The following contact matrix element is computed:

M_{12} = \nabla_v \vec{F}_{12} + \lambda \Delta t \nabla_x \vec{F}_{12}.

where \lambda (implicitness) modulates the kind of semi-implicit integration method (0 for explicit Euler, 0.5 for Cranck-Nicholson and 1 for implicit Euler).

ComputeContactMatrixElementStatic

Sets the contact matrix element between p1 and p2 for an implicitly solved system K x = F that represents a static linear spring network. This computes the contact matrix element:

M_{12} = \nabla_x \vec{F}_{12}.

.

ComputeContactMatrixElementUnderdamped

Sets the contact matrix element between particle p1 and p2 for an underdamped system in which a linear system M \cdot \vec{v} = \vec{F} is solved.

A contact matrix element is computed as:

M_{12} = \lambda \Delta t ( \nabla_v \vec{F}_{12} + \Delta t \nabla_x \vec{F}_{12}  ).

where \lambda (implicitness) modulates the kind of semi-implicit integration method (0 for explicit Euler, 0.5 for Cranck-Nicholson and 1 for implicit Euler).

DMT_Capsule

Special extension to the N_DMT contact model chain element that adds an additional contact moment for the interaction between two capsules. This moment arises for capsule-capsule contacts since the tractions in the contact ‘ellipse’ for such a case are not radially symmetric, thus resulting in a net contact moment. The exact computation of the contact ellipse between two arbitrary capsules is a complex and computationally demanding problem, so this chain element provides a simplified and heuristic approximation. It assumes that the moment is generated by the DMT contact force \vec{F}_{DMT} = 2 \pi \hat{R} w, and the ‘average lever arm’ in the DMT contact circle. The lever arm L is:

L = \frac{2}{3} \sqrt{ \hat{R} \delta }.

and the net contact moment:

\vec{M}_{DMT} = - 2 f\,F_{DMT}\,L\,\frac{ \vec{t}_1 \cdot \vec{t}_2}{ (\vec{t}_1 \cdot \vec{t}_1)\cdot(\vec{t}_2 \cdot \vec{t}_2) }\,(\vec{t}_1 \times \vec{t}_2).

where f is the weighting property torque_factor and where \vec{t}_i = \vec{x}_{i,1} - \vec{x}_{i,0} are the axis vectors of both capsules. No moment is computed when the contact point does not project within one of the two capsule segments, since end-sphere end-sphere contact is symmetric and does not generate any moment.

DeformableTriangleDof_1

Chain element that provides access to degrees of freedom (DOF) of a deformable triangle, i.e. its nodes. A weights member is kept that determines the relative contribution of each node. These weights are computed by a previous contact model chain element, e.g. using DistributeForcesAndMomentsToTriangle_1_Nodes and DistributeForcesAndMomentsToTriangle_2_Nodes.

DeformableTriangleDof_2

Chain element that provides access to degrees of freedom (DOF) of a deformable triangle, i.e. its nodes. A weights member is kept that determines the relative contribution of each node. These weights are computed by a previous contact model chain element, e.g. using DistributeForcesAndMomentsToTriangle_1_Nodes and DistributeForcesAndMomentsToTriangle_2_Nodes.

Deformable_Capsule_1_Rotational_Data

Provide binding to arrays representing rotational data of deformable capsule.

  • t - unit vector representing the axis of the cylinder,

  • l - length of the cylinder,

  • F - force on the center of the cylinder,

  • T - effective torque on the cylinder

    \vec{T} = \sum\limits_i r_i \vec{F}_i,

    where r_i represent distance from the center of the cylinder to the point where force F_i acts, see article for further details.

  • gamma - vector of friction coefficients for translational degrees of freedom, where the first component is for the movement along the axis and the other two are for the movement perpendicular to the axis.

  • gamma_rot - vector of friction coefficients for rotational degrees of freedom, where the first component is for the deformation along the axis and the other two are for rotations perpendicular to the axis, see article.

Attention

Some commands assume that last two components of gamma vectors are equal.

Deformable_Capsule_2_Rotational_Data

Provide binding to arrays representing rotational data of deformable capsule.

  • t - unit vector representing the axis of the cylinder,

  • l - length of the cylinder,

  • F - force on the center of the cylinder,

  • T - effective torque on the cylinder

    \vec{T} = \sum\limits_i r_i \vec{F}_i,

    where r_i represent distance from the center of the cylinder to the point where force F_i acts, see article for further details.

  • gamma - vector of friction coefficients for translational degrees of freedom, where the first component is for the movement along the axis and the other two are for the movement perpendicular to the axis.

  • gamma_rot - vector of friction coefficients for rotational degrees of freedom, where the first component is for the deformation along the axis and the other two are for rotations perpendicular to the axis, see article.

Attention

Some commands assume that last two components of gamma vectors are equal.

Deformable_Cylinder_1_Data

Provide binding to arrays representing deformable cylinder.

  • x - positions of endpoints,
  • r - radius of cylinder,
  • v - velocities of endpoints,
  • F - forces on endpoints,
  • m - masses of endpoints,
  • vertexIndices - indeces of verteces representing the endpoints,
  • parentIndex - index of particle which we are part of.

Deformable_Cylinder_2_Data

Provide binding to arrays representing deformable cylinder.

  • x - positions of endpoints,
  • r - radius of cylinder,
  • v - velocities of endpoints,
  • F - forces on endpoints,
  • m - masses of endpoints,
  • vertexIndices - indeces of verteces representing the endpoints,
  • parentIndex - index of particle which we are part of.

Deformable_Line_1_Data

Provide binding to arrays representing deformable line.

  • x - positions of endpoints,
  • v - velocities of endpoints,
  • F - forces on endpoints,
  • m - masses of endpoints,
  • vertexIndices - indeces of verteces representing the endpoints,
  • parentIndex - index of particle which we are part of.

Deformable_Line_2_Data

Provide binding to arrays representing deformable line.

  • x - positions of endpoints,
  • v - velocities of endpoints,
  • F - forces on endpoints,
  • m - masses of endpoints,
  • vertexIndices - indeces of verteces representing the endpoints,
  • parentIndex - index of particle which we are part of.

Disk_Sphere

Contact geometry resolution for interaction between a ‘disk’ and a sphere. The disk is defined by two ‘control points’, and a radius. Depending on a template parameter Disk_Selector, the disk center is located at one of either control points. This mechanism allows us to easily modulate between the top and bottom disk of hollow or solid cylinders. The normal of the disk points either towards the other point if the radius of the disk is negative (i.e. it delimits a hollow cylinder, in which the spheres are situated) or it points in the opposite direction if the radius of the disk is negative (i.e. it delimits a solid cylinder, where the spheres are around. This resolver does not compute any contact_fraction and solely computes a valid contact if the contact point is within the disk and the overlap is positive, or invalidates the contact otherwise.

DistributeForcesAndMomentsToTriangle_1_Nodes

Computation of a distribution vector d2 which contains the amount of normal force that should be added to each of the three triangle nodes of p2. This vector is computed based on a total normal force, an arbitrary point (e.g. the contact point) and a vector of total moments wrt this arbitrary point (e.g. computed by a pressure integrator). With this information a unique set of three forces can be computed in such a way that both the total sum of forces as well as the sum of moments stays conserved. Adding this chain element is essential for any of the contact models that compute contact forces with a pressure integrator and that have the mesh nodes as degrees of freedom (e.g. the deformable cell model).

Note

For more information about this computation, see https://doi.org/10.1371/journal.pcbi.1003267 equation 19.

DistributeForcesAndMomentsToTriangle_2_Nodes

Computation of a distribution vector d2 which contains the amount of normal force that should be added to each of the three triangle nodes of p2. This vector is computed based on a total normal force, an arbitrary point (e.g. the contact point) and a vector of total moments wrt this arbitrary point (e.g. computed by a pressure integrator). With this information a unique set of three forces can be computed in such a way that both the total sum of forces as well as the sum of moments stays conserved. Adding this chain element is essential for any of the contact models that compute contact forces with a pressure integrator and that have the mesh nodes as degrees of freedom (e.g. the deformable cell model).

Note

For more information about this computation, see https://doi.org/10.1371/journal.pcbi.1003267 equation 19.

ElasticBendingForce_Triangles

Model for bending rigidity of a closed elastic sheet described by the macroscopic model of Helfrich, for which the energy is given as:

E_c = \frac{k_c}{2} \int_A { (C_1 + C_2 - 2 C_0)^2 }dA + E_0

with local principal curvatures C_1 and C_2, spontaneous curvature C_0, and macroscopic bending rigidity k_c. The microscopic model for bending between two triangles is given by:

E_b = k_b [ 1-\cos( \theta - \theta_0 ) ].

with bending rigidity k_b (k_bend), where it can be shown that for a sphere:

k_b = \frac{2}{\sqrt{3}} k_c

For a triangle-triangle common edge contact, this force model computes 4 nodal forces: Fl1, Fl2, Fc1 and Fc2 for lever l and common c nodes of two adjacent triangles 1 and 2. The derivation of these forces is based on the change of bending energy for small displacements s, i.e. the following expression:

f(s) = -k_b \sin( \theta - \theta_0) \frac{ \partial \theta }{\partial s}

The exact expressions for the four nodal forces can be found in Multiscale Modeling of Blood Flow and Soft Matter, by D.A. Fedosov pages 221-223.

Attention

This model should be more accurate to describe ‘macroscopic’ bending of a membrane than the more simplistic but often-used linear bending model in LinearBendingForce_Triangles, where a moment is distributed to nodal forces in such a way that the ‘common’ points have identical forces. It should be noted that generally, the difference with this simplistic model is very small. Moreover, numerically, the simpler LinearBendingForce_Triangles is usually better behaved e.g. at preventing ‘wrinkling’ at large flattened regions of a membrane.

ElasticShellBendingForce_Triangles

Model that computes elastic bending forces in a triangulated shell. The shell is characterized by a thickness t, a Young’s modulus E and a poisson ratio \nu. For a pair of two adjacent triangles, a bending coefficient k_b is computed as:

k_b = \frac{ E t_{12}^3 }{ 12 (1-\nu^2) }

where t_{12} denotes the effective thickness which is provided by a previous chain element e.g. TriangleShellData. Given k_b, bending forces for four nodes of the two adjacent triangles are computed as described in ElasticBendingForce_Triangles.

Note

The current value of k_b can be consulted as a read-only property with name k_bend.

ElasticShellSpring

In-plane stretching and compression model for a linearly elastic shell with small thickness. This computes spring-like interaction between nodes in a triangulated mesh, where the stiffness of the springs is dependent on the shell material properties, and the area of the two triangles that are adjacent to the node. The instantaneous stiffness k of each node-node connection is computed as (Van Gelder model):

k_{12} = E\,\frac{ A_{t1}\,t_{t1} + A_{t2}\,t_{t2} }{L_0^2}

where E is the Young’s modulus of the thin shell material, A and t are the thickness and area of adjacent triangle t1 and t2, and L_0 is the spring resting length. The equations in this model are adapted from Lloyd and Harders and Kot, Nagahashi and Szymczak.

Warning

It has been shown that this model alone is only valid for materials with Poisson number 1/3. For other Poisson factors, this model should be supplemented with the so-called Kirkwood-model, which adds angular springs between triangle edges of the mesh. It should be noted that the error made for other values of \nu is typically reasonably small.

Feedback

Returns the overlap computed by previous chain-elements as ‘feedback’ information for the contactdetector. The usual interpretation of this value is ‘distance until actual contact’ when the feedback has a negative sign. The magnitude of this feedback is compared to the keep-distance of the contactdetector to decide whether or not the contact needs to be evaluated in the next iteration. Positive values indicate that the contact is ‘real’ and should not be removed regardless of keep-distance of the contactdetector.

FeedbackDistance

Returns minus the distance computed by previous chain-elements as ‘feedback’ information for the contactdetector. This feedback is compared to the keep-distance of the contactdetector to decide whether or not the contact needs to be evaluated in the next iteration.

FeedbackDividing

Sets the ‘feedback’ for the contact model to zero (indicating to always keep the contact regardless of the interparticle distance) if p1 and p2 are particles belonging to one dividing cell.

FrictionMatrixForce

Explictly computes a force that corresponds with given friction \nabla_v \vec{F} based on the relative contact velocity \vec{v}_{12}, and adds it to the total contact model force. The following contact force is computed by this contact model chain element:

\vec{F}_{12} = -\nabla_v \vec{F} \vec{v}_{12}

FrictionMatrixForce_1

Explictly computes the partial force contribution that corresponds to the current contact matrix element M and the velocity of the particle from pc2, and adds it to particle 1. This is not a regular conservative force, so it requires another ‘assymmetric’ contact model chain element, for example, one that adds the equivalent friction element to particle 2. The following force is assembled only to the degree of freedom *j$ of particle 1:

vec{F}_{1,j} = w_{1,j} , M ,vec{v}_{2,c}

where w_{1,j} is the weight of the j-th dof of p1 in this contact, and where \vec{v}_{2,c} is the velocity of particle 2 at the contact point c.

FrictionMatrixForce_2

Explictly computes the partial force contribution that corresponds to the current contact matrix element M and the velocity of the particle from pc1, and adds it to particle 2. This is not a regular conservative force, so it requires another ‘assymmetric’ contact model chain element, for example, one that adds the equivalent friction element to particle 1. The following force is assembled only to the degree of freedom *j$ of particle 1:

vec{F}_{2,j} = w_{2,j} , M ,vec{v}_{1,c}

where w_{2,j} is the weight of the j-th dof of p2 in this contact, and where \vec{v}_{1,c} is the velocity of particle 2 at the contact point c.

KeepTimeStep

Provides access to the simulation timestep \Delta t to all subsequent contact model chain elements.

LinearBendingForce_Triangles

Compute a linear bending moment for two triangles from a mesh that share a common edge. These contacts are typically provided by a FixedListContactDetector and keep a contact state with a reference angle, which indicates the resting position. In a linearized version around the resting angle, we compute a bending moment as:

M = k_b ( \theta - \theta_0 ).

with angle and reference angle \theta and \theta_0. Subsequently, bending forces at the four nodes possessed by the two triangles are computed as:

\vec{F}_{l1} =  \frac{ M }{h_1}\,\vec{\hat{n}}_1,\hspace{2cm} \vec{F}_{l2} =  \frac{ M }{h_2}\,\vec{\hat{n}}_2,

for the non-shared ‘lever’ nodes of triangle 1 and triangle 2, where \vec{\hat{n}} are the triangle normals and h are the (orthogonal) triangle ‘heights’ from the common edge towards the lever nodes. For the shared nodes located at the common edge we compute

\vec{F}_{c1} = \vec{F}_{c2} = -\frac{1}{2}(\vec{F}_{l1}+\vec{F}_{l2})

The later formula ensures that the sum of all forces is always zero.

Attention

From a deviation analysis based on the energetics, a slightly different formula would be computed, for example to compute the bending rigidity of a solid thin ‘sheet’. The bending moment would contain a sinus (which is linearized for small angle deviations) and the distribution to node forces would be slightly different. The macroscopic net result of these changes will be very minor in almost all cases (and vanish for perfect equilateral triangles). A very detailed derivation of the proper expressions for thin sheet bending rigidity of a meshed triangulated surface can be found in Fedosov et al., 2010.

To use this approach, select a model that has the chain element ElasticBendingForce_Triangles or a derivative of this, such as ElasticShellBendingForce_Triangles

LinearDamper

Linear dashpot with damping coefficient c. A force in the direction of the normal unit vector \vec{\hat{n}} is computed based on the normal relative contact velocity v_n as:

\vec{F}_c = - c v_n \vec{\hat{n}}.

This force is always computed for each contact offered by the contact detector, so it is ideally suited for putting dashpots in spring networks with FixedListContactDetectors.

LinearDamperOverlap

Linear dashpot with damping coefficient c. A force in the direction of the normal unit vector \vec{\hat{n}} is computed based on the normal relative contact velocity v_n as:

\vec{F}_c = - c v_n \vec{\hat{n}}.

Contrary to LinearDamper, this force is only computed if the overlap distance is larger than zero. Also, the force is weighted by the contact_fraction to ensure mesh independence of the total computed force for meshed particles.

MaxwellFluidShell

Maxwell relaxation for spring-like contact models. Stresses (contact forces) are relaxed by dynamically changing the equilibrium distance (reference length) of internal springs with a relaxation time. This contact model chain element depends on a previous chain element that describes the elastic stress as a function of a Young’s modulus E, and a previous chain element that provides a function effectiveThickness, which computes the thickness of the spring-like shell element for which the relaxation is to be computed. For each connection, a dashpot coefficient c is then computed as:

c = \frac{ \tau_R\,t_{12}\,E }{\sqrt{3}}

where \tau_R is the Maxwell relaxation time. This dashpot coefficient is passed to SerialDashpot which actually relaxes the spring element.

NGon_CapsuleSphere12

Contact geometry resolution for interaction between a capsule and an n-gon, explicitly taking into account the contact of the spheres at the capsule’s end points. As a result of this, this resolver computes and stores additional data that should be separately handled by the contact force model. In essence this contact model should handle the contact with the two end spheres as well as the submerged part of the cylinder’s hull. This resolver is extremely involved, but roughly, the required geometrical information is computed in the following steps:

  1. Project the cylinder axis on the n-gon
  2. Resolve the contact with the end spheres and the plane, including values modify_factor1 and modify_factor2, which indicate the fraction of the spheres that are actually overlapping (similar to contact_fraction for ‘regular’ resolvers).
  3. Compute a ‘reduced projected cylinder, which is a projection that discounts for the part of the projected cylinder axis that is already taking into account by the spheres.
  4. Compute the intersection of this projected cylinder axis and the n-gon.
  5. Compute the ‘min_overlap’ and ‘max_overlap’ from the overlaps at the endpoints of this intersection.
  6. Compute the ‘angle’ of the cylinder axis wrt the n-gon, and store it for later use as \alpha.
  7. Compute a contact point as the properly weighted (with overlap_sphere1, overlap_sphere2, min_overlap and max_overlap) of the intersection points.
  8. Compute the overlap from the sub-overlaps. Physical contact is only present when this net overlap is larger than zero.

Attention

This resolver is scheduled for a major overhaul.

NGon_NGon

Contact geometry resolution for interaction between two n-gons. This resolver computes all necessary information to prepare for the explicit integration of a pressure formulation over their projected intersection. The general methodology for this is described in Smeets et al. 2016. Briefly, the resolver works in the following steps:

  1. Compute a contact plane by subtracting and normalizing the respective n-gon normals.
  2. Compute the unique intersection line between the two planes to which the n-gons belong. save its direction \vec{\hat{l}} and a point \vec{s}.
  3. Determine the angle (\alpha) between these planes and save it for later use.
  4. Project both n-gons to the shared contact plane.
  5. Compute the counter-clockwise sorted intersection polygon between the two projected n-gons.
  6. Compute the subset of this intersection polygon that is on the contact side of the intersection line. This subset is saved for later use by further chain elements.
  7. Compute an overlap distance, an ‘overlap-weighted’ contact point, and a contact area based on the computed physical overlap area.

Crucially, this resolver provides a function OverlapAtPosition for further chain elements that allow them to probe the local overlap at arbitrary locations within intersection, for example for integrating a contact pressure formulation. This local overlap is computed as:

\delta(\vec{x}) = 2 \tan{\alpha}\left[ (\vec{x}-\vec{s})\cdot (\vec{\hat{n}} \times \vec{\hat{l}}) \right].

In the same vein, a function to compute the relative contact velocity at any arbitrary point in the intersection polygon is provided as:

\vec{v}_r(x) = \vec{v}_{r2}(\vec{x}) - \vec{v}_{r1}(\vec{x}),

where \vec{v}_{r1} and \vec{v}_{r2} are provided by the underlying data class.

NGon_Sphere

Contact geometry resolution for interaction between an n-gon (pc1) and a sphere (pc2). The resolver determines overlap, contactPoint, contact_fraction and normalUnitVector in the following steps:

  1. Calculation of contact point and contact circle with the infinite plane
  2. Compute the intersection between the contact circle and the n-gon
  3. If the intersection is empty, reject. Else, compute the proper value contactPoint, overlap and contact_fraction.

Property reject_large_overlap allows the rejection of contacts with overlap larger than the radius. This is essential to prevent overlap with the ‘opposite’ sides of (partially) convex geometries. Property flip_normals will flip the normals if the sphere is more than one radius submerged.

Note

The effective contact radius is simply \hat{R} = r_2.

N_Budding

Computes a cell division force that drives two particles, one the a ‘mother’ and the other a ‘bud’, apart. The separating force is obtained by computing a ‘division error’ i.e. the difference in the desired separation at time t after initiation of the budding, and the actual separation of the two sub-particles. A linear controller with stiffness divstiffrep in repulsion and divstiffatt in attraction is used to come up with a normal contact force between two particles that compose a budding cell.

Note

A division model has special status within a contact model chain, since it usually is not ‘additive’ with other chain elements, but rather seeks to replace their actions specifically for p1-p2 pairs that are each others division partner. Hence it sets the force vector contained in model data rather than adds up to it, and it sets the overlap to - max_double while simultaneously ensuring that rejectContact = 0, to avoid other contact force models to come up with non-zero forces. This contact model chain element relies on the presence of a special FeedbackDividing that ensures that division partners are always retained by the contact detector, and a specialized friction matrix assembler, that, similar to this force model, overwrites the friction matrix element gamma, rather than adding up to it.

N_ConstantForce

Normal contact force with magnitude given by property f.

Note

To preserve the ‘constant’ force total magnitude, we keep into account the contact fraction taken up by this particular contact of the complete (infinite) geometry primitive, if applicable.

Note

Depending on template argument CheckOverlap, we either always compute the constant force, or only when the overlap distance is larger than or equal to zero.

N_Contractile

Contact force between 2D disk-like cells with a repulsive force due to (implicit) cell-substrate adhesion and a attractive force due to cell-cell contractility. This is the contact model that was used in the 2016 PNAS paper where cells are described as colloidal active particles with contact inhibtion of locomotion (see Smeets et al. 2016). The parameter Wa represents substrate adhesion which effectively results in cell-cell repulsion, and the parameter Wc represents adhesive tension.

Note

If allow_dewet is set to True, we optionally allow for a very crude way of ‘dewetting’ which will cause the interaction force to drop to zero. This parameter was set to True in the original paper. If False, the adhesion energy is obviously actually higher, so for most cases this should be set to its default.

N_Contractile_dewet

Adaptation of N_Contractile which models ‘dewetting’ of 2D disk-like particles in a bit more elaborate manner. Each cell contains a ‘dewet-level’ which is the ‘layer’ index of stacking that this dewetted cell is situated in. A de-wet event will randomly increase the dewet level of one of the two interacting cells. In the absence of cells ‘below’, a higher cell will drop its dewet level as far as it can without colliding in a dewet status again.

Note

To make this work properly, the command ProcessDewetLevelsCommand should be added to the simulation, which properly administrates the dewet levels.

N_Contractile_line

Simple adaptation of N_Contractile which models contractile entities (spring-like ) in a line. This model was used for the supplementary information in Smeets et al. 2016

N_DMT

Derjaguin-Muller-Toporov (DMT) model for the adhesive contact between curved elastic asperities. The DMT model is an opposite limit to the JKR limit and holds for low values of the Tabor coefficient (i.e. small, stiff spheres). The total normal interaction force is given by:

F_{DMT} = \frac{4}{3} \hat{E} \sqrt{\hat{R}} \delta^{3/2} - 2 \pi \sigma \hat{R},

where \sigma (attrConst) specifies the strength of adhesion and where

\hat{E} = \frac{E_1 E_2}{E_2 (1-\nu_1^2)+E_1 (1-\nu_2^2)},

and \hat{R} is the effective contact radius (see the appropriate resolver chain element).

Note

For more details on adhesion models in contact mechanics, see Johnson and Greenwood.

N_DampedHertz

Hertz contact model with additional normal viscous damping (often called Kuwabara and Kono model). We compute a normal contact force:

F_H = \frac{4}{3}\hat{E}\sqrt{\hat{R}}(\delta^{3/2}+\hat{A}\dot{\delta}\sqrt{\delta}),

where \hat{A} = (A_1 + A_2) / 2 is a normal damping constant (in units of time) and \dot{\delta} is the normal relative contact velocity. As in the Hertz model, we define :

\hat{E} = \frac{E_1 E_2}{E_2 (1-\nu_1^2)+E_1 (1-\nu_2^2)},

and the effective contact radius \hat{R} (see the appropriate resolver chain element). The computed force has the direction of the normal unit vector as obtained from the resolver.

Note

A contact force is only computed if the overlap distance \delta is larger than or equal to zero.

Note

For detailed explanation, see ‘Computational Granular Dynamics’ of Pöschel, page 19.

N_DampedHertz_Int

Integrated normal damped Hertz contact force for interactions with ‘rounded’ meshed particles. This chain element is identical in parameters and physics as the regular N_DampedHertz which is used for more simple contact geometries. However, we explicitly integrate (using quadrature rules) the Hertz contact pressure over the intersection area between two contacting particles, if the contact is sufficiently large. The Hertz contact pressure including viscous damping at radius r is defined as

p_H(r) = \frac{\hat{E}}{\pi} \sqrt{a^2 - r^2} \left( \frac{3 \hat{A} \dot{\delta}}{a^2} + \frac{2}{\hat{R}} \right),

for a contact radius a = \sqrt{\delta \hat{R}} and where \hat{E}, \hat{R} and \hat{A} are as defined in N_DampedHertz. This method relies on the evaluation of r which is computed as the in-contact-plane distance from the contact point of the particles’ encompassing spheres. Hence, this contact force model can only be used for any of the ‘Roundedtriangle’ contact geometries. Moreover, it relies on a numerical integration over an intersection polygon, which should be set by an appropriate contact geometry resolver. During this numerical integration, each computed pressure is translated into a point force. This point force is summed up to a total sum of forces and a sum of moments wrt contactPointSphere. For very small overlaps (i.e. a is small), numerical integration would be inaccurate. Instead, we compute the ‘full-sphere’ Hertz contact force directly and weight it with the fraction of the total ‘full-sphere’ contact area that is taken up by intersection. Doing so neglects the spatial heterogeneity of the contact pressure, but this can be justified since the contact area is comparatively small anyway. The precise (heuristic) criterion for directly computing the force is

0.1 (N_i-2) N_Q \pi a^2 < A_i

where N_i and A_i are the number of corners and the area of intersection, and N_Q=16 is the number of quadrature points (per triangle) for numerical integration. Care should be taken that an appropriate contact force/moment assembler is chosen that properly takes into account the total sum of forces as well as the sum of moments wrt the sphere-contact-point.

Note

For numerical integration, the (convex) polygon intersection is subdivided into triangles over which the pressures are integrated separately.

Note

An elaborate description of this method can be found in Smeets et al. 2015.

N_Delile

Normal contact force model for contact between biological cells as described by Delile et al.. The full mathematical description can be found in the supplementary information. In short, we compute an effective contact area and a weighted repulsive and contractile force based on the overlap distance between two spheres. The default constants are as described by the original paper. The contact area:

A_c = a ( r_{12} - r_{12}^\mathrm{max} )^2,

where r_{12}^\mathrm{max} is defined as:

r_{12}^\mathrm{max} = c_\mathrm{max} ( R_1 + R_2 ),

The contact force reads:

\vec{F}_{12} = - w\,A_c\,( r_{12} - r_{12}^\mathrm{eq} ) \vec{n}_{12},

with resting distance:

r_{12}^\mathrm{eq} = c_\mathrm{eq} (R_1 + R_2)

.The constant $c_mathrm{eq}$ is by default the area ratio of a disk to a hexagon $pi / (2sqrt{3}) approx 0.9069$.

N_Dividing

Computes a cell division force that drives two particles, which together compose a cell undergoing cytokinesis, apart. The separating force is obtained by computing a ‘division error’ i.e. the difference in the desired separation at time t after initiation of the cytokinesis, and the actual separation of the two sub-particles. A linear controller with stiffness divstiffrep in repulsion and divstiffatt in attraction is used to come up with a normal contact force between two particles that compose a dividing cell.

Note

A division model has special status within a contact model chain, since it usually is not ‘additive’ with other chain elements, but rather seeks to replace their actions specifically for p1-p2 pairs that are each others division partner. Hence it sets the force vector contained in model data rather than adds up to it, and it sets the overlap to - max_double while simultaneously ensuring that rejectContact = 0, to avoid other contact force models to come up with non-zero forces. This contact model chain element relies on the presence of a special FeedbackDividing that ensures that division partners are always retained by the contact detector, and a specialized friction matrix assembler, that, similar to this force model, overwrites the friction matrix element gamma, rather than adding up to it.

N_DividingResistance

Contact friction between two particles that represent a cell undergoing cytokinesis. This friction is added as a diagonal 3x3 matrix with identical elements gamma_dividing.

Note

To ensure compatibility with other chain elements that could potentiall add to friction as well, this chain element does not add to but set the contact data gamma element. This also implies that N_DividingResistance should be the last friction model in the contact model chain.

N_Hertz

Hertz repulsive force for the elastic contact between two bodies. We compute a normal contact force that depends on overlap distance \delta as

F_H = \frac{4}{3} \hat{E} \sqrt{\hat{R}} \delta^{3/2},

where

\hat{E} = \frac{E_1 E_2}{E_2 (1-\nu_1^2)+E_1 (1-\nu_2^2)},

and \hat{R} is the effective contact radius (see the appropriate resolver chain element). The computed force has the direction of the normal unit vector as obtained from the resolver.

Note

A contact force is only computed if the overlap distance is larger than or equal to zero.

Tip

For a Hertzian contact force that includes normal damping, look up N_DampedHertz.

N_Hertz_Int

Integrated Hertz contact force for interactions with ‘rounded’ meshed particles. This chain element is identical in parameters and physics as the regular N_Hertz which is used for more simple contact geometries. However, we explicitly integrate (using quadrature rules) the Hertz contact pressure over the intersection area between two contacting particles, if the contact is sufficiently large. The (elastic) Hertz contact pressure at radius r is defined as

p_H(r) = \frac{2\hat{E}}{\pi \hat{R}} \sqrt{a^2 - r^2},

for a contact radius a = \sqrt{\delta \hat{R}} and where \hat{E} and \hat{R} are as defined in N_Hertz. This method relies on the evaluation of r which is computed as the in-contact-plane distance from the contact point of the particles’ encompassing spheres. Hence, this contact force model can only be used for any of the ‘RoundedTriangle’ contact geometries. Moreover, it relies on a numerical integration over an intersection polygon, which should be set by an appropriate contact geometry resolver. During this numerical integration, each computed pressure is translated into a point force. This point force is summed up to a total sum of forces and a sum of moments wrt contactPointSphere. For very small overlaps (i.e. a is small), numerical integration would be inaccurate. Instead, we compute the ‘full-sphere’ Hertz contact force directly and weight it with the fraction of the total ‘full-sphere’ contact area that is taken up by intersection. Doing so neglects the spatial heterogeneity of the contact pressure, but this can be justified since the contact area is comparatively small anyway. The precise (heuristic) criterion for directly computing the force is

0.1 (N_i-2) N_Q \pi a^2 < A_i

where N_i and A_i are the number of corners and the area of intersection, and N_Q=16 is the number of quadrature points (per triangle) for numerical integration. Care should be taken that an appropriate contact force/moment assembler is chosen that properly takes into account the total sum of forces as well as the sum of moments wrt the sphere-contact-point.

Note

For numerical integration, the (convex) polygon intersection is subdivided into triangles over which the pressures are integrated separately.

Note

An elaborate description of this method can be found in Smeets et al. 2015.

N_Hertz_TruncatedTraction_geo_Int

A new ‘integrated’ contact model that combines a Hertzian repulsive pressure with a linearly increasing ‘truncated’ traction which vanishes beyond a rupture length. The contact state keeps track of a ‘stick radius’ (defined in the contact plane) which is reduced when two bodies are pulled out of contact, and can never increase. This way, plastic ‘failure’ of bonded entities can be simulated.

N_JKR

Johnson-Kendall-Roberts model for the adhesive contact between elastic curved asperities. The JKR model is an extension of the Hertz model for contact between curved bodies, in the limit of high values of the Tabor coefficient (i.e. large, compliant spheres with short-range attractive interactions) and is often used to model cell-cell adhesion. The relationship between the (adhesion-modified) contact radius a and the interaction force F_{JKR} is:

F_{JKR} = \frac{4 \hat{E}}{3 \hat{R}} a^3 - \sqrt{ 8 \pi \sigma \hat{E} a^3},

with specific adhesion energy \sigma (attrConst) and \hat{E} and \hat{R} as defined for Hertz. The contact radius is related to the overlap \delta as

\delta = \frac{a^2}{\hat{R}} - \sqrt{\frac{2 \pi \sigma}{\hat{E}} a}.

The solution for the (square root of the) contact radius requires solving a polynomial equation of degree 4, which is done iteratively in a Newton-Raphson scheme. Furthermore, N_JKR uses a special contact state with a bool contactEstablished_, which can be found in JKRContactState.

Note

Depending on the template parameter Param_Type, either the properties are ALL given as single values, or they are all expected to be in an array that is kept per particle, (i.e. an array for E1,E2,nu1,nu2,attrConst1 and attrConst2). Please verify the documentation of the complete contact model chain to check which version is used there.

Note

For more details on adhesion models in contact mechanics, see Johnson and Greenwood.

N_LinearAdhesionForce

Simple model for adhesion with an attractive force that linearly increases with separation between the particles. Detachment between the particles happens at user-specified rupture length d_detach. At this separation, the maximal attractive force occurs which is f_detach. The adhesion energy W is then trivially:

W = \frac{F_d d_d}{2}.

Note

This contact model keeps a contact state with a bool contactEstablished, which is set to True if positive overlap is reached and remains True as long as the separation between the particles is less than d_detach.

N_LinearForce

Force that increases linearly for a positive overlap with stiffness k. An optional linear viscous damping force is added with damping coefficient c.

Note

N_LinearForce takes into account contact_fraction to scale the force for partial contributions from a larger (super)particle geometry.

N_LinearPressure

Contact force model with a normal pressure that increases linearly with overlap distance. The contact force is given as:

\vec{F}_n = A_c \, (k \delta - c\,v_r) \vec{\hat{n}},

with linear layer stiffness k and damping c, and normal unit vector \vec{\hat{n}}.

Note

We ensure that this pressure is always positive. I.e. if a ‘negative’ damping force is higher than the positive repulsive force, we set the pressure to zero.

N_LinearPressure_Int

Integrated linear pressure contact force for interactions with ‘rounded’ meshed particles. This chain element is identical in parameters and physics as the regular N_LinearPressure which is used for more simple contact geometries. However, we explicitly integrate (using quadrature rules) the linear contact pressure over the intersection area between two contacting particles, if the contact is sufficiently large. For sphere-like particles, the overlap increases linearly with the contact area. Hence, we obtain a pressure that linearly increases with overlap if we integrate a uniform pressure p that scales with overlapSphere over the surface of the intersection polygon.This contact force model can only be used for any of the ‘RoundedTriangle’ contact geometries. and a geometric contact radius a should be available to delimit the contact area. Moreover, it relies on a numerical integration over an intersection polygon, which should be set by an appropriate contact geometry resolver. During this numerical integration, each computed pressure is translated into a point force. This point force is summed up to a total sum of forces and a sum of moments wrt contactPointSphere. For very small overlaps (i.e. a is small), numerical integration would be inaccurate. Instead, we compute the ‘full-sphere’ contact force directly and weight it with the fraction of the total ‘full-sphere’ contact area that is taken up by intersection. The precise (heuristic) criterion for directly computing the force is

0.1 (N_i-2) N_Q \pi a^2 < A_i

where N_i and A_i are the number of corners and the area of intersection, and N_Q=16 is the number of quadrature points (per triangle) for numerical integration. Care should be taken that an appropriate contact force/moment assembler is chosen that properly takes into account the total sum of forces as well as the sum of moments wrt the sphere-contact-point.

Note

For numerical integration, the (convex) polygon intersection is subdivided into triangles over which the pressures are integrated separately.

N_LubricationResistance

Add to the normal friction coefficient the resistance due to lubrication forces. For simple normal lubrication forces between spherical particles, the following approximation is used for the normal resistance \gamma_n:

\gamma_n = -\frac{6\pi\mu \hat{R}^2}{ \delta_{ij} }

with effective contact radius \hat{R}, and only for (negative) values of \delta < -d_\mathrm{min}.

Note

For numerical reasons, also a maximal distance d_\mathrm{max} must be specified, beyond which the resistance should be sufficiently small as to be considered negligible. This should be done to prevent artifacts from contact detection.

N_MaugisDugdale

Maugis-Dugdale (MD) model for the adhesive contact between curved elastic asperities. The MD model describes adhesive interactions in the full range of interaction types captured by the Tabor coefficient, ranging from DMT (low values) to JKR (high values). For a curved asperity, it defines a narrow zone of adhesive traction \sigma_0, related to the specific adhesion energy (attrConst) W by an effective_range of interaction h_0:

W = \sigma_0 h_0.

An integration of a formulation for the traction distribution yields following expression for the adhesive force:

F_{MD} = - 2 \sigma_0 \left\{ c^2 \arccos\left(\frac{a}{c}\right) + a\,\sqrt{c^2 - a^2} \right\}.

In this expression a and c are the radii of resp. intimate and adhesive contact with c \geq a. These radii are not constant but vary dependent on the overlap. Their computation is based on an iterative solution of two expressions that can be found in https://doi.org/10.1371/journal.pcbi.1003267 (Eq. 11 and 12). It can be easily shown that the DMT (N_DMT) and JKR (N_JKR) formulas follow directly when applying the proper limits of the Tabor coefficient.

Note

The total contact force is scaled with contact_fraction to ensure that the total force is conserved for meshed particles.

N_MaugisDugdale_Int

Integrated Maugis-Dugdale (MD) contact force for interactions with ‘rounded’ meshed particles. This chain element is identical in parameters and physics as the regular N_MaugisDugdale which is used for more simple contact geometries. However, we explicitly integrate (using quadrature rules) the Maugis-Dugdale contact pressure (and traction) over the intersection area between two contacting particles, if the contact is sufficiently large. Maugis-Dugdale theory is based on the combination of a Hertz repulsive pressure with a formulation for adhesive traction, which considers a narrow adhesive ‘ring’ at the edge (gap) of the contact area where adhesive forces are highest. The radius of the adhesive ring is c while the radius of the contact circle is a. c is always larger than or equal to a. The (elastic) Hertz contact pressure at radius r is defined as

p_H(r) = \frac{2\hat{E}}{\pi \hat{R}} \sqrt{a^2 - r^2},

where \hat{E} and \hat{R} are as defined in N_Hertz. The (adhesive) Maugis-Dugdale traction is given by

p_{MD}(r) = -\frac{\sigma_0}{\pi} \, \arccos{ \left\{ \frac{2a^2 - c^2 - r^2}{c^2 - r^2 } \right\} }

for r\leq a and

p_{MD}(r)= -\sigma_0

for a < r \leq c. The maximal adhesive traction \sigma_0 is related to the specific adhesion energy (attrConst) W by an effective_range of interaction h_0:

W = \sigma_0 h_0.

c and a are interdependent variable and their determination is based on an iterative solution of two expressions that can be found in https://doi.org/10.1371/journal.pcbi.1003267 (Eq. 11 and 12). Furthermore, this method relies on the evaluation of r which is computed as the in-contact-plane distance from the contact point of the particles’ encompassing spheres. Hence, this contact force model can only be used for any of the ‘RoundedTriangle’ contact geometries. Moreover, it relies on a numerical integration over an intersection polygon, which should be set by an appropriate contact geometry resolver. During this numerical integration, each computed pressure is translated into a point force. This point force is summed up to a total sum of forces and a sum of moments wrt contactPointSphere. For very small overlaps (i.e. a is small), numerical integration would be inaccurate. Instead, we compute the ‘full-sphere’ Maugis-Dugdale contact force (see N_MaugisDugdale) directly and weight it with the fraction of the total ‘full-sphere’ contact area that is taken up by intersection. Doing so neglects the spatial heterogeneity of the contact pressure, but this can be justified since the contact area is comparatively small anyway. The precise (heuristic) criterion for directly computing the force is

0.1 (N_i-2) N_Q \pi a^2 < A_i

where N_i and A_i are the number of corners and the area of intersection, and N_Q=16 is the number of quadrature points (per triangle) for numerical integration. Care should be taken that an appropriate contact force/moment assembler is chosen that properly takes into account the total sum of forces as well as the sum of moments wrt the sphere-contact-point.

Note

For numerical integration, the (convex) polygon intersection is subdivided into triangles over which the pressures are integrated separately.

Warning

This contact model is approximative of the adhesive behavior of solid rounded elastic and adhesive surfaces that are in contact. It uses the intimate contact radius that is based on the implicit elastic deformation of the solid body. As such, for two triangles at a ‘gap’ there is no guarantee that the edge of contact will be situated at their intersection. When the body is not actually a solid elastic body (which are usually modeled as rigid-bodies in DEM), this assumption is problematic. In this case, it is much preferable to use a modified formulation of this contact force model which takes the ‘geometric’ contact radius as the intimate radius of contact a. For this, see N_MaugisDugdale_geo_Int. Any of the deformable-cell-like contact models should use this ‘geo’ version of the same algorithm.

Note

An elaborate description of this method can be found in Odenthal et al. 2013.

N_MaugisDugdale_geo_Int

Integrated Maugis-Dugdale (MD) contact force for interactions with ‘rounded’ meshed particles that uses the geometric contact radius instead of the actual Maugis-Dugdale contact radius to delimit the contact circle. Apart from this different contact radius this chain element is almost identical to N_MaugisDugdale_Int. Maugis-Dugdale theory is based on the combination of a Hertz repulsive pressure with a formulation for adhesive traction, which considers a narrow adhesive ‘ring’ at the edge (gap) of the contact area where adhesive forces are highest. The radius of the adhesive ring is c while the radius of the geometric contact circle is a. c is always larger than or equal to a. The (elastic) Hertz contact pressure at radius r is defined as

p_H(r) = \frac{2\hat{E}}{\pi \hat{R}} \sqrt{a^2 - r^2},

where \hat{E} and \hat{R} are as defined in N_Hertz. The (adhesive) Maugis-Dugdale traction is given by

p_{MD}(r) = -\frac{\sigma_0}{\pi} \, \arccos{ \left\{ \frac{2a^2 - c^2 - r^2}{c^2 - r^2 } \right\} }

for r\leq a and

p_{MD}(r)= -\sigma_0

for a < r \leq c. The maximal adhesive traction \sigma_0 is related to the specific adhesion energy (attrConst) W by an effective_range of interaction h_0:

W = \sigma_0 h_0.

In Maugis-Dugdale theory, c and a are interdependent variable and their determination is based on an iterative solution of two expressions that can be found in https://doi.org/10.1371/journal.pcbi.1003267 (Eq. 11 and 12). Importantly, this ‘geo’-modified algorithm fixes a at the geometric contact radius and, using MD theory, computes c from this. In general, this method relies on the evaluation of r which is computed as the in-contact-plane distance from the contact point of the particles’ encompassing spheres. Hence, this contact force model can only be used for any of the ‘RoundedTriangle’ contact geometries. Moreover, it relies on a numerical integration over an intersection polygon, which should be set by an appropriate contact geometry resolver. During this numerical integration, each computed pressure is translated into a point force. This point force is summed up to a total sum of forces and a sum of moments wrt contactPointSphere. For very small overlaps (i.e. c is small), numerical integration would be inaccurate. Instead, we compute the ‘full-sphere’ Maugis-Dugdale contact force (see N_MaugisDugdale) directly and weight it with the fraction of the total ‘full-sphere’ contact area that is taken up by intersection. Doing so neglects the spatial heterogeneity of the contact pressure, but this can be justified since the contact area is comparatively small anyway. The precise (heuristic) criterion for directly computing the force is

0.1 (N_i-2) N_Q \pi c^2 < A_i

where N_i and A_i are the number of corners and the area of intersection, and N_Q=16 is the number of quadrature points (per triangle) for numerical integration. Care should be taken that an appropriate contact force/moment assembler is chosen that properly takes into account the total sum of forces as well as the sum of moments wrt the sphere-contact-point.

Note

For numerical integration, the (convex) polygon intersection is subdivided into triangles over which the pressures are integrated separately.

Note

The precise criterion for deciding whether the pressure should be integrated is slightly different than the one used e.g. in N_Hertz_Int: it uses the adhesive radius c instead of the geometric radius a to determine if the ‘full-sphere’ contact is sufficiently large.

Note

An elaborate description of this method can be found in Odenthal et al. 2013.

N_T_AreaWeightedResistance

Friction matrix model that adds a contact matrix friction to the contact model chain. The main parameters are gamma_normal and gamma_tangential. A friction element \Gamma is computed as (T=transpose):

\Gamma = A_c \left[\gamma_n \hat{\vec{n}}\hat{\vec{n}}^\mathrm{T}+\gamma_t\left( \vec{I} - \hat{\vec{n}}\hat{\vec{n}}^\mathrm{T} \right)\right]

where A_c is the contactArea provided by a function in one of the previous chain elements. Generally, \Gamma is a symmetric 3x3 sub-matrix which is to be stored in the ‘off-diagonal’ part of the friction matrix.

Attention

Since the area is used to weight the friction element, the units of the provided gamma_normal and gamma_tangential must be (Pa*s/m).

N_T_ConstantResistance

Friction matrix model that adds a contact matrix friction to the contact model chain. The main parameters are gamma_normal and gamma_tangential. A friction element \Gamma is computed as (T=transpose):

\Gamma = \gamma_n \hat{\vec{n}}\hat{\vec{n}}^\mathrm{T}+\gamma_t\left( \vec{I} - \hat{\vec{n}}\hat{\vec{n}}^\mathrm{T} \right)

Generally, \Gamma is a symmetric 3x3 sub-matrix which is to be stored in the ‘off-diagonal’ part of the friction matrix.

N_T_ConstantStiffness

Computes a constant stiffness tensor for linearized elastic collisions. The parameter is the constant linear stiffness k that must exist in a previous chain element. Based on k, a stiffness tensor is computed as:

\nabla_x \vec{F} = k \hat{\vec{n}}\hat{\vec{n}}^\mathrm{T}+ k \left( \vec{I} - \hat{\vec{n}}\hat{\vec{n}}^\mathrm{T} \right)

N_T_MaskedConstantResistance

Friction matrix model that adds a contact matrix friction to the contact model chain. The main parameters are gamma_normal and gamma_tangential, which are to be provided as symmetric matrices that contain friction parameters gamma at i and j for the interaction for particles with mask value i and j. This allows the user to differentiate between different ‘masked’ particle types. A friction element \Gamma is computed as (T=transpose):

\Gamma = \gamma_n \hat{\vec{n}}\hat{\vec{n}}^\mathrm{T}+\gamma_t\left( \vec{I} - \hat{\vec{n}}\hat{\vec{n}}^\mathrm{T} \right)

Generally, \Gamma is a symmetric 3x3 sub-matrix which is to be stored in the ‘off-diagonal’ part of the friction matrix.

N_T_MaterialSpecificAreaWeightedResistance

Friction matrix model that adds a contact matrix friction to the contact model chain. The main parameters are gamma_normal1, gamma_normal2, gamma_tangential1, gamma_tangential2 which can refer to either single values or arrays. A friction element \Gamma is computed as (T=transpose):

\Gamma = A_c \left[\gamma_n \hat{\vec{n}}\hat{\vec{n}}^\mathrm{T}+\gamma_t\left( \vec{I} - \hat{\vec{n}}\hat{\vec{n}}^\mathrm{T} \right)\right]

where A_c is the contactArea provided by a function in one of the previous chain elements. Generally, \Gamma is a symmetric 3x3 sub-matrix which is to be stored in the ‘off-diagonal’ part of the friction matrix.

Attention

Since the area is used to weight the friction element, the units of the provided friction coefficients must be (Pa*s/m).

N_T_MaterialSpecificConstantResistance

Friction matrix model that adds a contact matrix friction to the contact model chain. The main parameters are gamma_normal1, gamma_normal2, gamma_tangential1, gamma_tangential2 which can refer to either single values or arrays. A friction element \Gamma is computed as (T=transpose):

\Gamma = \gamma_n \hat{\vec{n}}\hat{\vec{n}}^\mathrm{T}+\gamma_t\left( \vec{I} - \hat{\vec{n}}\hat{\vec{n}}^\mathrm{T} \right)

Generally, \Gamma is a symmetric 3x3 sub-matrix which is to be stored in the ‘off-diagonal’ part of the friction matrix.

N_T_MeshWeightedResistance

Friction matrix model that adds a contact matrix friction to the contact model chain. The main parameters are gamma_normal and gamma_tangential. A friction element \Gamma is computed as (T=transpose):

\Gamma = w_{12} \left[ \gamma_n \hat{\vec{n}}\hat{\vec{n}}^\mathrm{T}\gamma_t\left( \vec{I} - \hat{\vec{n}}\hat{\vec{n}}^\mathrm{T} \right)\right]

Generally, \Gamma is a symmetric 3x3 sub-matrix which is to be stored in the ‘off-diagonal’ part of the friction matrix. The weights w are computed as

w_{12} = \frac{1}{2}( w_1 + w_2 )

with individual weights for node i computed as:

w_i = \frac{ \sin(\pi / N_i) }{ \sin(2\pi / N_i)  }

where N refers to the number of triangles (and nodes) connected to node i.

Warning

Maybe a bit counter-intuitively, gamma_normal represents the IN-plane viscous damping while gamma_tangential represents the ‘out-of-plane’ or ‘normal-to-the-plane’ viscous damping. This is to keep consistency with the rest of the contact model, which is effectively between two nodes, and where the ‘normal’ direction is the direction vector connecting the two nodes.

Note

For the most typical case, a 6-fold connectivity, or roughly equilateral triangles, the computed weight can be easily shown to be \frac{1}{\sqrt{3}}.

N_T_SPH_Resistance

Friction matrix model that computes contact matrix friction for SPH-based contact models. The parameters gamma_normal and gamma_tangential cannot be provided by the user but are computed in previous contact model chain elements, where they are based on a viscosity model. A friction element \Gamma is computed as (T=transpose):

\Gamma = \gamma_n \hat{\vec{n}}\hat{\vec{n}}^\mathrm{T}+\gamma_t\left( \vec{I} - \hat{\vec{n}}\hat{\vec{n}}^\mathrm{T} \right)

Generally, \Gamma is a symmetric 3x3 sub-matrix which is to be stored in the ‘off-diagonal’ part of the friction matrix.

N_T_ShellViscosity

Friction matrix model that adds a contact matrix friction to the contact model chain, given a ‘viscosity’ of the thin shell that this contact model describes. For a triangulated shell, a friction element \Gamma is computed as (T=transpose):

\Gamma =  \frac{\mu\,t_{12}}{\sqrt{3}}\left[ \hat{\vec{n}}\hat{\vec{n}}^\mathrm{T} + \left( \vec{I} - \hat{\vec{n}}\hat{\vec{n}}^\mathrm{T} \right)\right]

where \mu is the viscosity and t is the shell thickness, which is to be provided by a preceding contact model chain element. Generally, \Gamma is a symmetric 3x3 sub-matrix which is to be stored in the ‘off-diagonal’ part of the friction matrix. The pre-factor should be 1/\sqrt{3} for a regular triangulated mesh, and the effective thickness t_{12} is computed by a previous chain element TriangleShellData.

NoFeedback

No rejection feedback is given from the contact model to the contact detector. In practice, we set feedback to zero, ensuring that all executed contacts will remain within keep_distance of the contact detector, and none will be removed based on information from the contact model.

NodeShellData

Data class for node-node contacts to for shell models. Provides data access to triangle area, normal and thickness, and vertex indices of nodes and triangles required for contact resolution.

Node_Node_Shell

Contact resolution between two nodes that are part of a mesh defining a triangulated shell. This resolver is based on the resolution of a ‘spring’-like connection, and computes a distance and a normalUnitVector between two points in space p1 and p2. For new contacts, this resolver determines the indices of the two adjacent triangles of the edge connecting p1 and p2, as well as determines the indices of the indices of the two ‘lever’ points (points in triangle 1 and 2 which are not p1 or p2). A function is provided that computes the effective ‘thickness’ of the shell element as:

t_{12} = \frac{ (A_{t1}+A_{t2}) t_{t1}\,t_{t2} }{A_{t1}\,t_{t2} + A_{t2}\,t_{t1}},

with A and t the area and thickness of triangles t_1 and t_2.

Furthermore, a function computeElementVolume is provided which computes the volume of our shell element as:

V_{12} = A_{t1}\,t_{t1} + A_{t2}\,t_{t2}.

NormalStiffnessMatrixCorrectionForce

For semi-implicit integration schemes that solve a system M a = F, an additional correction contact force must be computed for linear stiffness elements based on the relative contact velocity. This chain element adds this contact force for only the normal component of the relative contact velocity:

vec{F}_{12} = - Delta t lambda nabla_x vec{F} vec{v}^n_{12}.

where \vec{v}^n_{12} is the normal relative contact velocity and \lambda is the implicitness (0, 0.5, or 1).

Point_1_Data

Provides access to the necessary arrays for point-like particles. The following arrays are opened in pc2: x, v and F.

Point_2_Data

Provides access to the necessary arrays for point-like particles. The following arrays are opened in pc2: x, v and F.

Point_Point

Contact geometry resolution for interaction between two points. Since points are considered, overlap has not meaning and the only thing this resolver does is determining the distance d_{12} as:

d_{12} = ||\vec{x}_2 - \vec{x}_1||,

and the normal unit vector

\vec{\hat{n}} = \frac{\vec{x}_2 - \vec{x}_1}{d_{12}}.

PresetForce

Contact force model that always adds a fixed applied force F_applied as long as the contact state contactEstablished is true. contactEstablished is set to true the first time that the overlap is larger than zero and is only set to false the first time that the separation between the particles is larger than property release_distance.

RemoveFromContactMatrix

Set the contact model state contact matrix element to a ‘zero’ element.

ResetForceData

Resets all forces and moments contained in ForceData to zero vectors. effectively cancelling the effects of any previous chain elements on contact forces and moments. Useful e.g. when these quantities will be later obtained by an implicit solver.

Rigid_NGon_1_Data

Provides access to the necessary arrays for n-gon geometries and their rigid-body parent. The template argument specifies the number of corners of the n-gon that is used.

  • Arrays opened in the parent of pc2: x, v, m, F, w, M and q.
  • Arrays opened in pc2: vertexIndices, parentIndex and normal.

Rigid_NGon_2_Data also provides a function to compute the relative contact velocity of particle 2 given a contact point (\vec{x}_c):

\vec{v}_{r2} = \vec{v}_2 + \vec{w}_2 \times (\vec{x}_c - \vec{x}_2).

Rigid_NGon_2_Data

Provides access to the necessary arrays for n-gon geometries and their rigid-body parent. The template argument specifies the number of corners of the n-gon that is used.

  • Arrays opened in the parent of pc2: x, v, m, F, w, M and q.
  • Arrays opened in pc2: vertexIndices, parentIndex and normal.

Rigid_NGon_2_Data also provides a function to compute the relative contact velocity of particle 2 given a contact point (\vec{x}_c):

\vec{v}_{r2} = \vec{v}_2 + \vec{w}_2 \times (\vec{x}_c - \vec{x}_2).

RoundedTriangle_Cylinder

Contact geometry resolution for interaction between a rounded triangle and a cylinder. The triangle must have a unique encompassing sphere computed for it, defined by a radius (inverse_curvature) and a center (center_sphere). Ensure that these data, as well as the triangle’s normal is computed before this contact model chain element is executed (either before in the contact model chain or - preferably - in separate commands). This resolver computes all necessary information to prepare for the explicit integration of a pressure formulation for curved asperities (e.g. Hertz-like). The general methodology for this is described in Smeets et al. 2015. Briefly the contact geometry is resolved in the following steps:

  1. Compute the overlap distance, contact point and normal unit vector based on the interaction between the triangle’s encompassing sphere and the cylinder. All these data are saved for use by further contact model chain elements.
  2. Compute the contact plane as the contact plane between the encompassing sphere and the cylinder.
  3. Project the triangle onto the contact plane.
  4. The intersection polygon intersection is simply this projected triangle.
  5. Compute the overlap distance and contact point based on the closest point within intersection to the sphere-cylinder contact point.

Note

Valid contacts are possible as long as computeEffectiveContactRadius yields a positive result. In other words, one of the two radii can be negative (concave) as long as \hat{R}=\frac{R_1 R_2 }{R_1 + R_2} > 0. For example, the cylinder can be hollow as long as the magnitude of its radius is larger than the radius of the encompassing sphere. In the previous formula, we have used the radius of the cylinder (which is not correct, since the curvature is not isotropic) and thereby introduced a small error in the simulation. In Smeets et al. 2015 we have quantified this error to be sufficiently small for most reasonable DEM simulations.

RoundedTriangle_NGon

Contact geometry resolution for interaction between a rounded triangle and an n-gon. The triangle must have a unique encompassing sphere computed for it, defined by a radius (inverse_curvature) and a center (center_sphere). Ensure that these data, as well as the triangle’s normal is computed before this contact model chain element is executed (either before in the contact model chain or - preferably - in separate commands). This resolver computes all necessary information to prepare for the explicit integration of a pressure formulation for curved asperities (e.g. Hertz-like). The general methodology for this is described in Smeets et al. 2015. Briefly the contact geometry is resolved in the following steps:

  1. Compute the overlap distance, contact point and normal unit vector based on the interaction between the triangle’s encompassing sphere and the n-gon. All these data are saved for use by further contact model chain elements.
  2. The contact plane is always the plane of the n-gon.
  3. Project the triangle onto the contact plane.
  4. Compute the intersection polygon of the projected triangle and the n-gon, and save it as intersection.
  5. Compute the overlap distance and contact point based on the closest point within intersection to the sphere-n-gon contact point.

Note

The effective contact radius is simply the inverse curvature of the triangle.

RoundedTriangle_RoundedTriangle

Contact geometry resolution for interaction between rounded triangles. These triangles must have a unique encompassing sphere computed for them, defined by a radius (inverse_curvature) and a center (center_sphere). Ensure that these data, as well as the triangle’s normal is computed before this contact model chain element is executed (either before in the contact model chain or - preferably - in separate commands). This resolver computes all necessary information to prepare for the explicit integration of a pressure formulation for curved asperities (e.g. Hertz-like). The general methodology for this is described in Smeets et al. 2015. Briefly the contact geometry is resolved in the following steps:

  1. Compute the overlap distance, contact point and normal unit vector based on the contact geometry of both triangles’ encompassing spheres. All these data are saved for use by further contact model chain elements.
  2. Compute the contact plane as the contact plane of the encompassing spheres.
  3. Project both triangles onto the contact plane.
  4. Compute the intersection polygon of both projected triangles, and save it as intersection.
  5. Compute the overlap distance and contact point based on the closest point within intersection to the sphere-sphere contact point.

Note

The effective contact radius has the same shape as the one for two spheres, with the radii and position determined by inverse_curvature and center_sphere.

Note

Valid contacts are possible as long as computeEffectiveContactRadius yields a positive result. In other words, one of the two radii can be negative (concave) as long as \hat{R}=\frac{R_1 R_2 }{R_1 + R_2} > 0.

RoundedTriangle_Sphere

Contact geometry resolution for interaction between a rounded triangle and a sphere. The triangle must have a unique encompassing sphere computed for it, defined by a radius (inverse_curvature) and a center (center_sphere). Ensure that these data, as well as the triangle’s normal is computed before this contact model chain element is executed (either before in the contact model chain or - preferably - in separate commands). This resolver computes all necessary information to prepare for the explicit integration of a pressure formulation for curved asperities (e.g. Hertz-like). The general methodology for this is described in Smeets et al. 2015. Briefly the contact geometry is resolved in the following steps:

  1. Compute the overlap distance, contact point and normal unit vector based on the interaction between the triangle’s encompassing sphere and the sphere. All these data are saved for use by further contact model chain elements.
  2. Compute the contact plane as the contact plane of these spheres.
  3. Project the triangle onto the contact plane.
  4. The intersection polygon intersection is simply this projected triangle.
  5. Compute the overlap distance and contact point based on the closest point within intersection to the sphere-sphere contact point.

Note

The effective contact radius has the same shape as the one for two spheres, with one radius and position determined by inverse_curvature and center_sphere.

Note

Valid contacts are possible as long as computeEffectiveContactRadius yields a positive result. In other words, one of the two radii can be negative (concave) as long as \hat{R}=\frac{R_1 R_2 }{R_1 + R_2} > 0.

SerialDashpot

Serial dashpot for contact models with spring-like elements. This contact model chain element relies on a previous chain element that computes an elastic contact force based on the deviation from a reference distance. This chain element dynamically relaxes the reference distance of the spring such that the equilibrium distance L_0 changes as:

\frac{d L_0}{ d t } = - \frac{f_{12}}{ c }

if L/L_0 is larger than slip_threshold (which by default is zero, so always true).

Sphere0_1_Data

Provides access to the necessary arrays for massless sphere particles with no rotational degrees of freedom. The following arrays are opened in pc2: x, r, v and F.

Sphere0_2_Data

Provides access to the necessary arrays for massless sphere particles with no rotational degrees of freedom. The following arrays are opened in pc2: x, r, v and F.

Sphere_1_Data

Provides access to the necessary arrays for sphere particles in pc2. Next to the data provided by Sphere0_2_Data (x, r, v and F), arrays associated to mass and rotational degrees of freedom are opened: m, w, M, and q.

Sphere_2_Data

Provides access to the necessary arrays for sphere particles in pc2. Next to the data provided by Sphere0_2_Data (x, r, v and F), arrays associated to mass and rotational degrees of freedom are opened: m, w, M, and q.

Sphere_Point

Contact geometry resolution for interactions between a sphere and a point. We compute a distance (not stored) as:

d_{12} = ||\vec{x}_2 - \vec{x}_1||,

a normal unit vector:

\vec{\hat{n}} = \frac{\vec{x}_2 - \vec{x}_1}{d_{12}},

and an overlap:

\delta = r_1 - d_{12}

Sphere_Sphere

Contact geometry resolution for interaction between two spheres. Given two spheres p1 and p2 with radius r1 and r2 and position x1 and x2, this resolver computes a unit length vector that points from p1 towards p2, computes an overlap distance which is negative if the spheres are separated and positive if the spheres collide, and computes a contactPoint, which is located at the point in space where the penetration of the spheres is maximal. If the overlap distance is negative, the contactPoint is set to (0,0,0). Sphere_Sphere also sets a rejectContact flag which is zero if the contact is a valid physical contact and is one if the contact should be rejected by other elements later in the contact model chain. The following functions are provided for use by later chain elements: computeEffectiveContactRadius, computeNormalRelativeContactVelocity, computeGeometricContactRadius and computeOverlapArea.

Note

Valid sphere-sphere contacts are possible as long as computeEffectiveContactRadius yields a positive result. In other words, one of the two radii can be negative (concave) as long as \hat{R}=\frac{R_1 R_2 }{R_1 + R_2} > 0.

Spring

Linear spring with force proportional to deviation from a reference distance. The magnitude of the force is:

F = k ( d_0 - d ),

for a distance d (obtained from the resolver) and a reference distance d_0. Two separate values k_stretch and k_compress can be given to enact different stiffnesses in the compressed or the stretched states. If the equilibrium distance is not provided, it will be set to the current distance when the contact model is first executed. A rupture length d_max can be provided. The spring will break (permanently) if it elongates beyond d_max.

SpringStiffness

Computes a contact stiffness tensor for linear spring systems. The parameter is the constant linear spring stiffness k that must exist in a previous chain element. Based on k, a stiffness tensor is computed as:

\nabla_x \vec{F} = k [ \vec{n} \otimes \vec{n} + (1 - \frac{ L_0 }{ L } ) \cdot ( \mathbf{I} - \vec{n} \otimes \vec{n}  ) ].

for springs with resting length L_0 and length L, and where \vec{n} is the contact normal unit vector.

Spring_Damper

Combined spring and dasphot interaction. The magnitude of the force is:

F = k ( d_0 - d ) - c v_r.

for a distance d (obtained from the resolver), a reference distance d_0 and a dasphot coefficient c. The normal relative contact velocity is computed by a function computeNormalRelativeContactVelocity that must be provided by the contact geometry resolution chain element.

Note

If the reference distance was not user-initalized, it is set to the current distance whenever this given contact is first executed.

StiffnessMatrixCorrectionForce

For semi-implicit integration schemes that solve a system M a = F, an additional correction contact force must be computed for linear stiffness elements based on the relative contact velocity. This chain element adds the contact force

vec{F}_{12} = - Delta t lambda nabla_x vec{F} vec{v}_{12}.

where \vec{v}_{12} is the relative contact velocity and \lambda is the implicitness (0, 0.5, or 1).

StiffnessMatrixReferenceDistanceForce

When implicitly solving spring networks, the springs’ reference distance adds a partial force term to the right hand side of the system K x = F. This chain element computes this force based on a reference distance L_0 as:

\vec{F}_{12} = L_0\,\nabla_x \vec{F}\,\vec{n}_{12}

with \vec{n}_{12} the spring’s normal unit vector.

T_CoulombFriction

Haff and Werner model for tangential Coulomb friction with friction coefficient mu. The static regime is enforced with a linear dashpot with damping coefficient c_t. Its value has only numerical meaning and should be chosen as high as possible, while retaining numerical stability. The Coulomb friction force can be expressed as

\vec{F}_c = -\mathrm{min}( c_t||\vec{v}_t||, \mu||\vec{F}_n||)\frac{\vec{v}_t}{||\vec{v}_t||},

where \vec{F}_n is the total normal force, and \vec{v}_t the tangential relative contact velocity (see the appropriate resolver).

Note

The viscous dashpot can never maintain a truly static regime. Therefore, this model is not suited for simulating the stability of static piles/packed beds at long timescales. Rather it should be used for cases where connectivity changes dynamically in the timescale of interest. For a simple friction model that better enforces the static regime, look up T_CundallStrack.

Attention

Since T_CoulombFriction relies on the normal force F_n, it is important that F_n has been completely computed at this point. All contact model chain elements that add to normal forces must be positioned before T_CoulombFriction in the chain!

T_CoulombFriction_Int

Haff and Werner model for tangential Coulomb friction interactions with ‘rounded’ meshed particles. This chain element is identical in parameters and physics as the regular T_CoulombFriction which is used for more simple contact geometries. However, we properly weight the friction force given the fraction of the ‘full-sphere’ contact area that is taken up by this ‘part’ of of the composing bodies’ mesh. Very importantly, the ‘full-sphere’ total normal force is used to decide if the particle is in the stick or in the slipping regime. The criterion for slipping:

c_t ||\vec{v}_t|| > \mu F_{n,\mathrm{sphere}}

only uses ‘global’ and no ‘local’ or intersection-specific information. This contact model chain element only works if this fnSphere has been computed before.

Note

Please look up T_CoulombFriction for information about the model parameters and its limitations.

Note

Since the tangential friction force is assumed to be constant over the contact area anyway, this contact force model does not need to do any actual numerical integration. However, its suffix ‘Int’ indicates that it needs an integration-family normal contact force model in front of it.

Note

An elaborate description of this method can be found in Smeets et al. 2015.

T_CundallStrack

Cundall and Strack model for tangential Coulomb friction. The static regime is enforced with a linear dashpot with damping coefficient c_t, as well as a linear spring with stiffness k_t. The spring only exists implicitly and its elongation \vec{\zeta} (kept in a contact state) is integrated as:

\vec{\zeta}(t) = \int_{t_K}^{t} \vec{v}_t(t')dt',

for tangential relative contact velocity \vec{v}_t. Moreover we have ensured that the spring’s resting length is reinitialized at the time of first contact t_K. The tangential force is then:

\vec{F}_t = -\mathrm{min}( k_t ||\vec{\zeta}|| + c_t ||\vec{v}_t||,\mu ||\vec{F}_n||)\frac{\vec{v}_t}{||\vec{v}_t||}.

The constants c_t and k_t have only numerical meaning and should be determined as high as possible while retaining numerical stability. For more information see ‘Computational Granular Dynamics’ of Pöschel, page 24.

Note

This Chain element needs a special contact state that keeps an R3::Vector with elongation.

Note

Our implementation has a few minor modifications w.r.t. the reference given in Pöschel and the original paper. Firstly, we allow different friction coefficients for the static and dynamic regime mu_s and mu_d and secondly we scale the contact force with the contact_fraction to ensure mesh independence.

Attention

Since T_CundallStrack relies on the normal force F_n, it is important that F_n has been completely computed at this point. All contact model chain elements that add to normal forces must be positioned before T_CundallStrack in the chain!

TorsionalCoulombFriction

Computation of torsional Coulomb friction which results from the relative rotation of two interacting particles around their connecting axis (normal unit vector). This model is derived for spheres and is only an approximation for other simple shapes. For complex (triangulated) shapes where contact forces are computed in an integration scheme, this torsional force is added naturally and does not need a dedicated chain element. In the static regime, the torsional moment is:

\vec{M}_t = -\frac{c_t w_{12} A_c}{ 2 \pi } \vec{n}_{12},

with w_{12} the relative angular velocity around normal unit vector \vec{n}_{12}. In the sliding regime:

\vec{M}_t = -\frac{2}{3}\,\mu\,F_N\,\hat{R}\,\vec{n}_{12}

The slip criterion is:

\left| \frac{ c w_{12} A_c }{ 2 \pi}\right| > \left| \frac{2}{3}\,\mu\,F_N\,\hat{R}\right|.

TriangleShellData

Data class for shell data of triangle geometries. Provides access to the ‘thickness’ array which should be present as a Scalar for each triangle. This chain element also provides the function effectiveThickness, which gives the current value of the shell thickness based on the thickness and area of both triangles that share a common edge:

t_{12} = \frac{ (A_1 + A_2) t_1 t_2 } { A_1 t_2 + A_2 t_1 }

This expression is an area-weighted geometric mean that preserves the thickness to t if constant, and gives bigger ‘weight’ to the contribution of a triangle with larger area.

Triangle_Triangle_CommonEdge_Angle

Contact geometry resolution for the bending interaction between two triangles that share a common edge. Given the shared edge (two ‘common’ nodes) and two ‘lever points’, this resolver computes the angle between the triangles. If the contact has not been executed before and the reference angle was not set by the user, it will be initialized to the current angle. This resolver also provides the function getAngleNormal which computes the properly signed angle given the two normals of the triangles (which must be computed at this point).

Triangle_Triangle_CommonEdge_CosAngle

Contact geometry resolution for the bending interaction between two triangles that share a common edge. Given the shared edge (two ‘common’ nodes) and two ‘lever points’, this resolver computes the cosine and sine of the angle between the triangles. If the contact has not been executed before and the reference angle was not set by the user, it will be initialized to the current angle. This resolver is written using the same symbols and definitions as defined in Multiscale Modeling of Blood Flow and Soft Matter, by D.A. Fedosov pages 221-223., keeping the non-normalized triangle normals \xi and \zeta, and storing both sine and cosine (but not the angle itself) as sin_angle and cos_angle.


In order to be able to use this module import it like this:

import mpacts.contact.models.chain_element_documenter
#or assign it to a shorter name
import mpacts.contact.models.chain_element_documenter as cha

Contents

doc() -> str :
Returns the docstring of this module (which is the documentation of all chain elements described here