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 non-fragile (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:

  1. \(\bullet_n\) for normal sense (translations along the local \(x\)-axis);

  2. \(\bullet_t\) for tangential sense (translations in the local \(yz\)-plane; may be 2d vectors if indicated);

  3. \(\bullet_w\) for twisting sense (rotations around the local \(x\)-axis);

  4. \(\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 geometry-dependent 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 2-vector 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\))

\begin{align*} k_n&=\frac{EA}{L}, & k_t&=\frac{GA}{L}, \\ k_w&=\frac{GJ}{L}, & k_r&=\frac{EI}{L} \end{align*}

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

\[ \begin{align}\begin{aligned}\alpha_w=\frac{1}{2\pi},\\\alpha_r=\frac{1}{4\pi}.\end{aligned}\end{align} \]

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 per-sense 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

  1. Bond breaks (i.e. the contact becomes non-cohesive) 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 and brkN bits of bonds set (rather than zero).

  2. When breakage in any sense occurs, the bond is broken in all senses at once and becomes fully unboded.

  3. 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 (size-independent stress) is computed as

(26)\[c_n=\underbrace{l\left(\frac{l_1}{E_1}+\frac{l_2}{E_2}\right)^{-1}}_{E'}\eps_{bn}\]

\(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 3-vector in IceMat.beta, in this order) as dimensionless scaling parameters:

\begin{align*} F_{n}& > F_{nb}, & F_{nb}&=c_n A, \\ |\vec{F}_{t}|& > F_{tb}, & F_{tb}&=\beta_t F_{nb}=\beta_t c_n A, \\ |T_{w}|& > T_{wb}, & T_{wb}&=\beta_w F_{tb}A^{\frac{1}{2}}=\beta_w \beta_t c_n A^{\frac{3}{2}}, \\ |\vec{T}_{b}|& > T_{rb}, & T_{rb}&=\beta_r F_{nb}A^{\frac{1}{2}}=\beta_r c_n A^{\frac{3}{2}}. \end{align*}

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 out-of-box. 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:

  1. 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.

  2. 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.

\begin{align*} F_{ty}&=-\min(0,F_n\tan\phi)\\ T_{wy}&=-\sqrt{A/\pi} \min(0,F_n\tan\phi) \\ T_{ry}&=-\sqrt{A/\pi} \min(0,F_n\mu) \end{align*}

Note that all these values are non-negative, 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:

\begin{align*} |\vec{F}_t|&>F_{ty} & \Longrightarrow && \vec{F}_t&\leftarrow \normalized{\vec{F}_t} F_{ty}, \\ |T_w|&>T_{wy} & \Longrightarrow && T_w&\leftarrow \normalized{T_w} T_{wy}, \\ |\vec{T}_r|&>T_{ry} & \Longrightarrow && \vec{T}_r&\leftarrow \normalized{\vec{T}_r} T_{ry}. \end{align*}

Nomenclature

Symbol

Unit

Variable

Algorithm

Meaning

\(E\)

Pa

ElastMat.young

geometry-weighted

Young’s modulus

\(A\)

\(\mathrm{m^2}\)

L6Geom.contA

computed

contact area

\(k_t/k_n\)

–

FrictMat.ktDivKn

averaged

factor to compute \(k_t\) from \(k_n\)

\(\tan\phi\)

–

FrictMat.tanPhi, FrictPhys.tanPhi

minimum

friction angle

\((\alpha_w,\alpha_r)\)

–

IceMat.alpha

averaged

factors for computing \(k_w\), \(k_r\) from \(k_t\), \(k_n\)

\(\eps_{bn}\)

–

IceMat.breakN

averaged

normal strain where cohesion stress is reached

\(c_n\)

Pa

temporary

computed from (26)

normal cohesion value

\((\beta_t,\beta_w,\beta_r)\)

–

IceMat.beta

averaged

factors for computing cohesions from \(c_n\)

\((F_{nb},F_{tb}), (T_{wb},T_{rb})\)

(N ,N, Nm, Nm)

IcePhys.brkNT, IcePhys.brkWR

computed

limit forces for bond breakage

\(\mu\)

–

IceMat.mu, IcePhys.mu

averaged

kinetic (rolling) friction coefficient

\(k_n, k_t, k_w, k_r\)

N/m, N/m, N, N

FrictPhys.kn, FrictPhys.kt, IcePhys.kWR

computed

stifffnesses in all 4 senses

Tip

Report issues or inclarities to github.