Ice contact model¶
Note
Developement of this model is commercially sponsored (identity undisclosed).
The Ice
model was conceived for ice simulations originally, but may be suitable for a wide range of materials. It considers elasticity, cohesion and friction in 4 senses (normal, shear, twisting, rolling); each of the senses can be independently fragile (“unbonding” upon reaching critical force/torque value – see below) or nonfragile (plastic slip taking place when critical force/torque is reached).
The model consists of a set of classes for material parameters (IceMat
, IceMatState
, IcePhys
, Cp2_IceMat_IcePhys
, Law2_L6Geom_IcePhys
) but consumes L6Geom
, generic contact geometry object. Therefore, it works for any combination of shapes
handled by L6Geom
.
Values related to respective degrees of freedom are denoted with the following indices:
\(\bullet_n\) for normal sense (translations along the local \(x\)axis);
\(\bullet_t\) for tangential sense (translations in the local \(yz\)plane; may be 2d vectors if indicated);
\(\bullet_w\) for twisting sense (rotations around the local \(x\)axis);
\(\bullet_r\) for rolling sense (rotations around local \(y\) and \(z\) axes; may be 2d vectors if indicated).
Stiffnesses¶
There are 4 stiffness values:
normal and tangential stiffnesses \(k_n\) (
kn
) and \(k_t\) (kt
) are computed the same as in Linear (Cundall) contact model.twisting stiffness \(k_w=\alpha_w k_t A\), where \(\alpha_w\) is material parameter and \(A\) is geometrydependent
contact area
.rolling stiffness \(k_r=\alpha_r k_n A\), where \(\alpha_r\) is material parameter.
\(\alpha_w\) and \(\alpha_r\) are stored together (in this order) as 2vector in IceMat.alpha
.
Note
Twisting stiffness is proportional to tangential stiffness and rolling stiffness is proportional to normal stiffness. (This is unlike for limit breakage forces, where normal–twist and tangential–rolling are proportional.) The rationale is that analogy with isotropic circular beam, where the rescpective stiffnesses are related to the same moduli (Young’s modulus \(E\) and shear modulus \(G\))
with \(A=\pi R^2\), \(J=\frac{\pi R^4}{2}\), \(I=\frac{\pi R^4}{4}\), \(G=\frac{E}{2(1+\nu)}\). Thus, to emulate isotropic beam on the contact, the coefficients would be
Note
Stiffness values are the same for bonded and unbonded contacts. Those concepts are orthogonal: stiffness is \(\frac{\partial F}{\partial u}\) (tangent stiffness, in this model identical to secant stiffness) while bonds determine the \(u\) range over which the \(F(u)\) function is linear (or elastic, physically speaking).
Bonds¶
The contact state can be bonded (cohesive) or unbonded in any of the 4 senses, independently. The initial bonding state (stored persense as bitmask in bonds
) is determined by the Cp2_IceMat_IcePhys
functor, and the possibility for breakage in 4 different senses is also initially set there.
Breakage¶
Bond breaks (i.e. the contact becomes noncohesive) when two conditions are met simultaneously:
force/torque in any sense exceeds limit value (see below);
the bond is breakable (fragile) in that sense, as indicated by flag for given sense (stored again in
bonds
).
E.g. for the normal sense, those conditions would be \(F_{n}>F_{nb}\) (see below) for the first point, and having both
bondN
andbrkN
bits ofbonds
set (rather than zero).When breakage in any sense occurs, the bond is broken in all senses at once and becomes fully unboded.
The transition from unbonded to bonded state never occurs naturally (though it can be forced by hand).
Limit force values depend on cohesion parameters; the normal cohesion (sizeindependent stress) is computed as
\(E'\) being equivalent Young’s modulus and \(\eps_{bn}\) being “strain” (as in \(\eps_{bn}=\Delta l/l\)) at breakage in the normal sense (material parameter IceMat.breakN
), \(l=l_1+l_2\) is total center distance and radii, as used in Linear (Cundall) contact model (or equivalent measures used to distribute stiffness, when contacting particles are not spheres).
Cohesion values are only useful for senses which are both bonded and breakable, and the breakage condition is slightly different for different senses. The values are all computed from \(c_n\) using \(A\) for correct dimension and \(\beta_t\), \(\beta_w\), \(\beta_r\) (stored as 3vector in IceMat.beta
, in this order) as dimensionless scaling parameters:
Note that there is no absolute value in the first equation, as there is no breakage in compression (\(F_n<0\)).
Bond setup¶
Simulations often need to create different bonding for contacts existing in the starti configuration (at the very beginning of the simulation) and different bonding for particles which only meet later, as an effect of their movement during the simulation.
This scenario is supported outofbox. Starting configuration’s IcePhys.bonds
are specified by Cp2_IceMat_IcePhys.bonds0
, and this value is assigned only during the first few steps of the simulation – precisely until S.step
reaches step01
.
After that, all new contacts are assigned Cp2_IceMat_IcePhys.bonds1
instead; contacts created before retain their old (bonds0
) value, it is only updated when the contact is encountered for the very first time.
Bond syntax¶
bonds
are represented by bits of an integeral number. Each bit represents whether the contact is bonded or breakable in one sense. The order is as follows:
First 4 bits (counting from the least significant ones, i.e. from the right) represent bonds in the 1. normal (
bondN
), 2. tangential (bondT
), 3. twisting (bondW
), 4. rolling (bondR
) senses.Second 4 bits represent breakability in those senses again, i.e.
brkN
,brkT
,brkW
,brkR
. It is understood that breakability bits are only useful in the sense where there is a bond.
Numbers can be written as binary literal in Python; thus, for instance, if the starting configuration were to bond in all senses, and break in the twisting sense only, one would set:
# the bits are (read from the left here):
# brkR=0, brkW=0, brkT=1, brkN=0
# bondR=1, bondW=1, bondT=1, bondN=1
bonds0=0b00101111
Plasticity¶
Plastic force limiters (yield values, noted \(\bullet_y\)) apply only for senses which are currently not bonded (be they broken, or simply never bonded at all). If force/torque exceeds respective yield force/torque, it is limited to that yield value (retaining its direction).
There are two plastic parameters, friction angle \(\phi\) (tanPhi
) and kinetic friction \(\mu\) (mu
), used to compute yield values. Note that the use of \(\min(0,F_n \dots)\) implies that the yield values are always zero in tension, therefore the behavior is ideally plastic in that case.
Note that all these values are nonnegative, since \(\min(0,F_n)\leq0\) and \(\mu\geq0\), \(\tan\phi\geq0\).
For each unbonded sense, when the yield condition is satisfied, the corresponding force is reduced to return to the yield value:
Nomenclature¶
Symbol 
Unit 
Variable 
Algorithm 
Meaning 

\(E\) 
Pa 
geometryweighted 
Young’s modulus 

\(A\) 
\(\mathrm{m^2}\) 
computed 
contact area 

\(k_t/k_n\) 
– 
averaged 
factor to compute \(k_t\) from \(k_n\) 

\(\tan\phi\) 
– 
minimum 
friction angle 

\((\alpha_w,\alpha_r)\) 
– 
averaged 
factors for computing \(k_w\), \(k_r\) from \(k_t\), \(k_n\) 

\(\eps_{bn}\) 
– 
averaged 
normal strain where cohesion stress is reached 

\(c_n\) 
Pa 
temporary 
computed from (26) 
normal cohesion value 
\((\beta_t,\beta_w,\beta_r)\) 
– 
averaged 
factors for computing cohesions from \(c_n\) 

\((F_{nb},F_{tb}), (T_{wb},T_{rb})\) 
(N ,N, Nm, Nm) 
computed 
limit forces for bond breakage 

\(\mu\) 
– 
averaged 
kinetic (rolling) friction coefficient 

\(k_n, k_t, k_w, k_r\) 
N/m, N/m, N, N 
computed 
stifffnesses in all 4 senses 
Tip
Report issues or inclarities to github.