Subspring Network Model

The subspring network model is referred to in commands and FISH by the name subspringnetwork.

Introduction

The subspring network model can be used to simulate unbonded and bonded systems. Like the Spring Network Model, the bonded formulation is based on the Rigid Body Spring Network (RBSN) lattice model introduced by [Bolander1998] based on the work of [Kawai1978].

RBSN is a lattice formulation that can eliminate particle-scale heterogeneity and provide elastically-homogeneous material behavior in bonded Voronoi block models. [Asahina2015], [Asahina2017], [Rasmussen2018b] and [Rasmussen2022] have expanded the capabilities of RBSN by introducing a fictitious stress scheme capable of generating bonded-block models that directly match the macroscale elastic properties of isotropic, transversely isotropic and bi-modular materials without the need for time-consuming calibrations. [Rasmussen2018a], [Rasmussen2018b] and [Rasmussen2019] applied the RBSN to model hard and anisotropic rocks as well as failure modes of deep underground excavations, including rock fall and spalling. Building on these developments, [Rasmussen2021b] introduced the Hybrid Lattice / Discrete Element Method (Hybrid LDEM), which integrates RBSN into the explicit DEM scheme. This approach is highly effective in bonded block simulations of rocks, allowing accurate control over various aspects of the nonlinear behavior of rocks under compression loading, such as crack initiation and damage stress thresholds. Subsequently, [Rasmussen2022] introduced breakable particles, further expanding the range of macroscale properties that the combined RBSN-DEM bonded-block models can reproduce while preserving their simplicity.

Unlike the Spring Network Model, {2 in 2D; 3 in 3D} subcontacts are present. The force-displacement law is run independently for each subcontact and subcontact gaps are incrementally tracked; each subcontact can fail independently of the others, similar to the Flat-Joint Model formulation. The total contact forces and moments, summed over all subcontacts, are applied to the pieces in contact.

When all subcontacts are bonded, the behavior is identical to the Spring Network Model. In other words, one can effectively remove particle-scale heterogeneity in elastic response and match bonded material elastic response without calibration. When a subcontact is not bonded, the fictitious forces are zero for that subcontact, and the unbonded force-displacement law is used provided the subcontact gap is less than or equal to zero. All subcontacts have the same stiffnesses, areas, frictional properties, strengths, etc. However, the subcontact forces/moments and gaps may differ between subcontacts.

Behavior Summary

A subspring network bond combines three concepts: incremental forces/moments due to relative motion, the RBSN approach to computing contact stiffnesses/fictitious forces, and subcontacts. The first envisions a bond as a set of elastic springs with constant normal and shear stiffnesses, uniformly distributed over a {rectangular in 2D; circular in 3D} cross-section lying on the contact plane and centered at the subcontact point. Relative motion at the subcontact causes the development of a linear force and moment. The total force and moment, summed over all subcontacts, acts on the two contacting pieces. The RBSN paradigm provides a method to compute contact stiffnesses that are consistent—when bonded—with an elastic material with a target Young’s modulus and zero Poisson’s ratio. If a non-zero Poisson’s ratio is specified then a fictitious force, derived from linear elasticity, is computed and added to the total force to produce the desired Poisson effect for each bonded subcontact. Unlike the Spring Network Model, failure can occur at each subcontact independently. Unbonded or failed subcontacts utilize frictional strength parameters to govern slip with potential dilation and/or healing if the tracked subcontact gap is less than or equal to zero.

Activity-Deletion Criteria

A contact with the subpring network model is active if any subcontact is active. A subcontact is active if it is bonded or if its gap is less than or equal to zero. The force-displacement law is skipped for inactive subcontacts.

Contact Stiffness Computation

The subspring network model uses the method outlined here to compute stiffnesses for bonded assemblies. The stiffnesses are stored per unit area, and the subcontact stiffness are uniform.

Fictitious Stress Correction

The subspring network model uses the method outlined here to compute fictitious forces consistent with non-zero Poisson ratio materials. This fictitious force is computed once per contact if any subcontact is bonded, with {1/2 in 2D, 1/3 in 3D} of the value being contributed to the total force for each bonded subcontact. If all subcontacts are unbonded, no fictitious force is computed.

One important aspect of the fictitious force correction is that it requires a stress estimate from the bodies attached to the pieces at each end of the contact. This method was initially developed using Voronoi cells in contact along shared faces. That geometrical arrangement allows for contacts at many orientations around each Voronoi cell. As the stress is computed via the outer product of the contact branch vector and the contact force, having more contacts with diverse orientations provides a more robust estimate of the stress state. Now imagine that contacts exist along the facets of a very elongate {triangle in 2D; tetrahedra in 3D}. The number of contacts will be reduced relative to the Voronoi cell case and some of the contact normals may be nearly coincident, meaning that the body stress may be poorly defined. This can potentially lead to instabilities when bodies with few contacts are bonded and the fictitious force is computed. The subspring network model can automatically detect this condition by assessing the convergence of the associate bodies. If the convergence of either body is greater than a specified value (changed with the sn_conv property) for some number of force-displacement law computations (changed with the sn_cntconv property - it is -1 by default, disabling the check), the fictitious force correction is disabled.

Subcontact Area and Geometry

The area associated with each subcontact is {1/2 in 2D; 1/3 in 3D} the total contact area. Unlike in the Flat-Joint Model, one cannot control the number of subcontacts. The subcontacts are located on the contact plane. In 2D, the two subcontacts are located on either side of the contact position, at a distance of half the contact length (i.e., assuming a unit thickness). In 3D, the three subcontacts are located at a distance of the radius of an equivalent disk from the contact position, 120 degrees apart (i.e., equally distributed around the contact position). The locations of the subcontacts in 3D depend on the local coordinate system, which is arbitrarily assigned. This means that the 3D subcontacts do not represent the actual contacting surface for touching polyhedra, but are rather a useful abstraction of this area. Numerical experiments with tetrahedral and Voronoi packings demonstrated that this abstraction is capable of producing very similar failure responses to a model where the subcontacts are defined based on the faces in contacts. One ramification of this abstraction is that the subspring network model can be used for contacts between any types of pieces.

Bond State

Each subcontact carries an associated bond state

(1)\[\begin{split}B_{sub} = \begin{cases} \mbox{0: unbonded (U)} \\ \mbox{1: broken in tension (T)} \\ \mbox{2: broken in shear (S)} \\ \mbox{3: bonded (B)} \\ \mbox{4: healed (H)} \end{cases}\end{split}\]

Individual subcontacts cannot be unbonded, meaning that \(B_{sub}=0\) and \(B_{sub}\ne0\) cannot happen for subcontacts that are part of the same contact. Above corresponds to the case when sn_index is non-zero (i.e., when one is querying the bond state of a specific subcontact). By default, sn_index is zero, meaning that \(B_{sub}\) returns the sum of all subcontact bond states. Take the letters \(U\), \(B\), \(T\) and \(S\) to stand for an unbonded subcontact, a bonded subcontact, a subcontact broken in tension and for a subcontact broken in shear, respectively, as shown in the equation above. Healing is left out of this presentation as it is less common situation. In 2D, the \(B_{0}\) (sn_index = 0) represents:

(2)\[\begin{split}B_{0} = \begin{cases} \mbox{0: (U,U)} \\ \mbox{2: (T,T)} \\ \mbox{3: (T,S)} \\ \mbox{4: (T,B) or (S,S)} \\ \mbox{5: (S,B)} \\ \mbox{6: (B,B)} \end{cases}\end{split}\]

3D differs due to the number of subcontacts:

(3)\[\begin{split}B_{0} = \begin{cases} \mbox{0: (U,U,U)} \\ \mbox{3: (T,T,T)} \\ \mbox{4: (T,T,S)} \\ \mbox{5: (T,T,B) or (T,S,S)} \\ \mbox{6: (S,S,S) or (T,S,B)} \\ \mbox{7: (T,B,B) or (S,S,B)} \\ \mbox{8: (S,B,B)} \\ \mbox{9: (B,B,B)} \end{cases}\end{split}\]

The sn_states property returns state strings as in the equations above.

Force-Displacement Law

Each subcontact has an associated gap, force, moment and failure state. All subcontacts lie on the contact plane. Suppose that the vector from the contact position to the position of subcontact \(sub\), in the local coordinate system, is given by \(\mathbf{r}_{sub}\). The relative displacement increment at the subcontact position is computed as:

(4)\[\Delta \pmb{δ}_{sub} = \Delta \pmb{δ} + \mathbf{r}_{sub} \times \Delta \pmb{θ}\]

where \(\Delta \pmb{θ}\) is the increment of relative rotation and \(\Delta \pmb{δ}\) is the increment of relative displacement of Equation (12) of the i Contact Resolution section. The subcontact relative displacement increment can be decomposed as:

(5)\[\Delta \pmb{δ}_{sub} = \Delta \delta _{n,sub} {\kern 1pt} \hat{\mathbf{n}}_\mathbf{c} +\Delta \pmb{δ} _{\mathbf{s},sub}\]

The subcontact gap is incremented as:

(6)\[g_{sub} \mathrel{+}= \Delta \delta _{n,sub}\]

The total contact force/moment is the sum of the subcontact forces/moments:

(7)\[\mathbf{F_{c}} = \sum_{sub=1}^{ns} \mathbf{F}_{sub} ,\quad \mathbf{M_{c}} =\sum_{sub=1}^{ns} \mathbf{M}_{sub}\]

where {\(ns = 2\) in 2D; \(ns = 3\) in 3D}, \(\mathbf{F}_{sub}\) is the force of subcontact \(sub\) and \(\mathbf{M}_{sub}\) is the moment of subcontact \(sub\).

The force-displacement law for the subspring network model updates the subcontact force and moment:

(8)\[\mathbf{F}_{sub} =\mathbf{F}_{sub}^{l} +\mathbf{F}_{sub}^{d} +\mathbf{F}_{sub}^{f} ,\quad \mathbf{M}_{sub} =\mathbf{M}\]

where \(\mathbf{F}_{sub}^{l}\) is a linear force, \(\mathbf{F}_{sub}^{d}\) is a dashpot force, \(\mathbf{F}_{sub}^{f}\) is the fictitious force corresponding to the Poisson effect, and \(\mathbf{M}\) is a linear moment. Take \(\mathbf{F} = \mathbf{F}_{sub}^{l} + \mathbf{F}_{sub}^{f}\) as the total elastic subcontact force, including the Poisson correction. These forces and moments are updated as described below.

The elastic subcontact force is resolved into a normal and shear force:

(9)\[{\mathbf{F}} = {- F_n {\kern 1pt} \hat{\mathbf{n}}_\mathbf{c} +\mathbf{F_s}}\]

where \(F_n > 0\) is tension. The subcontact shear force lies on the contact plane and is expressed in the contact plane coordinate system:

(10)\[ {\mathbf{F_s}} = {F_{ss} {\kern 1pt} \hat{\mathbf{s}}_\mathbf{c} +F_{st} {\kern 1pt} \hat{\mathbf{t}}_\mathbf{c} \; \; \left({\rm 2D\; model:\; }F_{ss} \equiv 0\right)}\]

The subcontact gap is zero when the contact is created. Specifying the contact force, either directly or using a stress initialization strategy (i.e., ball tractions or rblock tractions), initializes the subcontact gap:

(11)\[g_{sub} = g + F_n / \left(k_n {\kern 1pt} A_{sub} \right)\]

where \(g\) is the contact gap, \(k_n\) is the normal stiffness per unit area, \(A\) is the contact area and \(A_{sub} = A / ns\). The reference gap is likewise initialized to the same value if a stress initialization strategy is used; this ensures that the force-displacement law continues to be called even if the physical contact gap becomes positive when stresses have been installed.

Subcontact Moment Update

The subcontact moment is updated, when the subcontact is active, via:

(12)\[\mathbf{M}_{sub} \mathrel{+}= \mathbf{k_r} \Delta \pmb{θ} A_{sub}\]

where \(\mathbf{k_r}\) is the rotational stiffness per unit area. In 2D the rotational stiffness is a scalar quantity, while in 3D it is given by

(13)\[\mathbf{k_{r}} = (k_{r_{n_c}},k_{r_{s_c}},k_{r_{t_c}})\]

where these quantities are defined in Equation (2) of the Spring Network Model section. The subcontact moment update applies when either the contact is bonded or if it is unbonded with \(g_{sub} \le 0\).

Subcontact Force Update - Unbonded Behavior

When unbonded, there is no fictitious (\(\mathbf{F}_{sub}^f = 0\)). The force-displacement law is computed for all unbonded subcontacts with \(g_{sub} \le 0\). The force is updated with the following steps.

  1. The normal subcontact force is updated incrementally as:

    (14)\[F_n := \max\left(F_n + k_n {\kern 1pt} A_{sub} \Delta \delta _{n,sub} , 0\right)\]

    Note:

    1. the unbonded behavior cannot sustain tensile normal force,

    2. the normal-force update is always incremental, and

    3. there are no fictitious forces when unbonded.

  2. Update \(\mathbf{F_s}\):

    The shear subcontact force \(\mathbf{F_s}\) is first updated with:

    (15)\[\mathbf{F_s} := \mathbf{F_s} - k_s {\kern 1pt} A{\kern 1pt} \Delta \pmb{δ} _{\mathbf{s},sub}\]

  3. Enforce the subcontact slip criteria:

    If the critical slip distance, \(c^{dc}\), is not specified, nor any entries in the cohesion table, a Coulomb friction criterion is enforced on the subcontact shear force such that:

    (16)\[\begin{split}\mathbf{F_s} = \left\{ \begin{array}{rl} \mathbf{F_s} , & \| \mathbf{F_s} \| \le - \mu F_n \\ - \mu F_n \bigl( \mathbf{F_s} / \| \mathbf{F_s} \| \bigr), & \mbox{otherwise.} \end{array} \right.\end{split}\]

    The subcontact slip state is updated as:

    (17)\[\begin{split}s_{sub} = \left\{ \begin{array}{rl} \mbox{true}, & \| \mathbf{F_s} \| = - \mu F_n \\ \mbox{false}, & \mbox{otherwise.} \end{array} \right.\end{split}\]

    If the subcontact slip state is true, then the subcontact is sliding. The entire contact is sliding if all active subcontacts are sliding.

    If \(c^{dc}\) is specified then the shear displacement is tracked and the shear force is linearly interpolated from the peak value to \(- \mu F_n \bigl( \mathbf{F_s} / \| \mathbf{F_s} \| \bigr)\) when the shear displacement achieves \(c^{dc}\). If multiple entries are specified in the cohesion table, the shear force is linearly interpolated between entries based on the displacement where entries are multiplicative factors on the failure strength. If healing is specified, then once sliding ceases the shear strength instantaneously recovers to its peak value and the subcontact is in a bonded state where it cannot sustain tension but can sustain shear.

    If slipping, subcontact dilation may occur. Dilation occurs until the total shear displacement reaches \(u_{cs}\), at which time it ceases. The increment of shear displacement over a timestep adds a normal force to the subcontact in the form

    (18)\[F_n = F_n + k_n {\kern 1pt} A{\kern 1pt} \tan (\psi \pi /180 ) \lVert \Delta \pmb{δ}_{\mathbf{s},sub} \rVert\]

    where \(\psi\) is the dilation angle in degrees.

  4. Update the subcontact dashpot forces:

    The subcontact dashpot force is updated as in the linear model.

Subcontact Force Update - Bonded Behavior

As mentioned above, the compute-stiffness method computes stiffnesses consistent with a zero Poisson’s ratio material. If the target Young’s modulus is uniform throughout the specimen when bonded, elastic heterogeneity can effectively be removed. In addition, specifying a non-zero Poisson’s ratio results in a fictitious stress correction that produces a realistic Poisson effect.

When bonded, the subcontact force is updated with the following steps.

  1. Update \(F_n\) and \(\mathbf{F_s}\) using an incremental formulation:

    (19)\[\begin{split}\begin{array}{rl} F_n & := F_n + k_n {\kern 1pt} A_{sub} \Delta \delta _{n,sub} + {F_n}^f/{ns} \\ \mathbf{F_s} & := \mathbf{F_s} - k_s {\kern 1pt} A_{sub}{\kern 1pt} \Delta \pmb{δ} _{\mathbf{s},sub} + \mathbf{{F_s}^f}/{ns} \end{array}\end{split}\]

  2. Update the subcontact normal (\(\sigma\), \(\sigma > 0\) is tension) and shear (\(\tau\)) stresses:

    (20)\[\begin{split}\begin{array}{l} \sigma = \displaystyle \frac{F_n} {A_{sub}} \\ \tau = \frac{\left\| \mathbf{F_s} \right\|}{A_{sub}} \end{array}\end{split}\]

    The effective subcontact normal stress, \(\sigma_e\), is used for subsequent failure computations, using the pore pressure \(P_p\):

    (21)\[\sigma_e = \sigma - P_p\]

  3. Update the bond state:

    Take the subcontact shear strength \(\tau_c = c - \sigma_c tan \phi\).

    • If the bond has healed then it cannot sustain tensile stresses (i.e., \(\sigma_c = 0\)) but can sustain shear stresses. If \(\sigma_e > 0\) then the bond breaks in tension (\(B_{sub}=1\)). If \(\tau > \tau_c\) then the bond breaks in shear (\(B_{sub}=2\)).

    • If the bond is intact (\(B_{sub}=3\)) and the effective normal stress exceeds the bond tensile strength (\(\sigma_e > \sigma_c\)) then the bond breaks in tension (\(B_{sub}=1\)).

    • If the bond is intact (\(B_{sub}=3\)) and the maximum shear stress exceeds the shear strength (\(\tau > \tau_c\)), then the bond fails in shear (\(B_{sub}=2\)).

Energy Partitions

The subspring network model provides three energy partitions that are updated per subcontact:

  • strain energy, \(E_{k}\), stored in the linear springs plus the fictitious stress correction;

  • slip energy, \(E_{\mu }\), defined as the total energy dissipated by frictional slip; and

  • dashpot energy, \(E_{\beta }\), defined as the total energy dissipated by the dashpots

Table 1: Spring Network Model Energy Partitions

Keyword

Symbol

Description

Range

Accumulated

spring-network:

energy‑strain

\(E_{k}\)

strain energy

\([0.0,+\infty)\)

NO

energy‑slip

\(E_{\mu}\)

total energy dissipated by slip

\((-\infty,0.0]\)

YES

Dashpot Group:

energy‑dashpot

\(E_{\beta}\)

total energy dissipated by dashpots

\((-\infty,0.0]\)

YES

If energy tracking is activated (see the model energy command), the energy partitions are updated as described below.

  1. Update the strain energy:

    (22)\[{\rm E} _{k} =\frac{1}{2} \left( \frac{F_n^{2} }{k_{n}{\kern 1pt} A_{sub}} + \frac{\left\| \mathbf{F_s} \right\| ^{2} }{k_{s} {\kern 1pt} A_{sub}} + \frac{\left\| \mathbf{M_b} \right\| ^{2} }{k_n {\kern 1pt} I ns} + \frac{M_t^{2} }{k_s {\kern 1pt} J ns} \right).\]

  2. Update the slip energy:

    (23)\[\begin{split}\begin{array}{l} E_{\mu } := E_{\mu} + \Delta E_{\mu}^s \\ {\rm where} \\ \qquad \begin{array}{l} \Delta E_{\mu}^s := - {\tfrac{1}{2}} \left(\left(\mathbf{F_{s}} \right)_{o} +\mathbf{F_{s}} \right)\cdot \Delta \pmb{δ} _{\mathbf{s},sub}^{\mu } \\ \end{array} \\ {\rm with} \\ \qquad \begin{array}{l} \Delta \pmb{δ} _{\mathbf{s},sub}^{\mu } =\Delta \pmb{δ} _{\mathbf{s},sub} -\Delta \pmb{δ} _{\mathbf{s},sub}^{k} =\Delta \pmb{δ} _{\mathbf{s},sub} -\left(\frac{\mathbf{F_{s}} -\left(\mathbf{F_{s}} \right)_{o} }{k_{s} {\kern 1pt} A_{sub}} \right) \\ \end{array} \end{array}\end{split}\]

  3. Update the dashpot energy:

    (24)\[E_{\beta } :=E_{\beta } - \mathbf{F}_{sub}^{d} \cdot \left(\dot{\pmb{δ} }{\kern 1pt} {\kern 1pt} \Delta t\right)\]

    where \(\dot{\pmb{δ} }\) is the relative translational velocity of Equation (10) of the i Contact Resolution section.

Properties

The properties defined by the spring-network model are listed in the table below for a concise reference; see the i Contact Properties section for a description of the information in the table columns. The mapping from the surface inheritable properties to the contact model properties is also discussed below.

Table 2: Spring-Network Model Properties

Keyword

Symbol

Description

Type

Range

Default

Modifiable

Inheritable

springnetwork | Model name

Spring-Network Group:

kn

\(k_n\)

Normal stiffness [stress/disp./area]

FLT

\([0.0,+\infty)\)

0.0

YES

YES

ks

\(k_s\)

Shear stiffness [stress/disp./area]

FLT

\([0.0,+\infty)\)

0.0

YES

YES

fric

\(\mu\)

Friction coefficient [-]

FLT

\([0.0,+\infty)\)

0.0

YES

YES

rgap

\(g_r\)

Reference gap [length]

FLT

\(\mathbb{R}\)

0.0

YES

NO

kratio

\(\kappa^*\)

Normal-to-shear stiffness ratio [-]

FLT

\([0.0,+\infty)^*\)

0.0\(^*\)

NO

N/A

\(\kappa^* \equiv \frac{k_n}{k_s}\)

krot

\(\mathbf{k_r}\)

Rotational stiffness [stress/disp./area]

VEC

\(\mathbb{R}^3\)

\(\mathbf{0}\)

YES

NO

sn_area

\(A\)

Bond area [length*length]

FLT

\((0.0,+\infty)\)

N/A

NO

NO

user_area

\(A\)

Constant bond area [length*length]

FLT

\((0.0,+\infty)\)

0.0

YES

NO

sn_ten

\(\sigma_c\)

Tensile strength [stress]

FLT

\([0.0,+\infty)\)

0.0

YES

NO

sn_coh

\(c\)

Cohesion [stress]

FLT

\([0.0,+\infty)\)

0.0

YES

NO

sn_cohres

\(c_{resid}\)

Residual cohesion [stress]

FLT

\([0.0,+\infty)\)

0.0

YES

NO

sn_cohdc

\(c^{dc}\)

Critical slip distance [length]

FLT

\([0.0,+\infty)\)

0.0

YES

NO

sn_cohtab

Cohesive strength multiplier as a function of shear displacement

VEC2

\(\mathbb{R}^2\)

YES

NO

sn_cntconv

Number of times the convergence limit has been exceeded to de-activate the fictitious force correction [~]. A negative value disables the check.

INT

\([0,+\infty)\)

-1

YES

NO

sn_conv

Convergence limit to de-activate the fictitious force correction [~]

FLT

\([0.0,+\infty)\)

1000.0

YES

NO

sn_fa

\(\phi\)

Friction angle [degrees]

FLT

\([0.0,90.0)\)

0.0

YES

NO

sn_dil

\(\psi\)

Dilation angle [degrees]

FLT

\([0.0,90.0)\)

0.0

YES

NO

sn_dilzero

\(u_{cs}\)

Critical dilation distance [length]

FLT

\([0.0,+\infty)\)

+infty

YES

NO

sn_gap

\(g_{sub}\)

Incremental gap [disp.]

FLT

\((-\infty,+\infty)\)

0

YES

NO

If sn_index is zero, the minimum subcontact gap is returned, otherwise the subcontact gap is returned

sn_heal

\(h\)

Healing

INT

{0,1}

0

YES

NO

sn_index

Subcontact index

INT

{0,1,2,3}

0

YES

NO

\(0-2/3 \ \mbox{(2D model/3D model)}\)

sn_state

\(B_{sub}\)

Bond state

INT

{0,1,2,3,4,5,6,7,8,9}

0

NO

NO

See equations (2) and (3) of the Bond State section

sn_states

Bond state string

STR

NO

NO

See equations (2) and (3) of the Bond State section

sn_shear

\(\tau_{c}\)

Bond shear strength [stress]

FLT

\([0.0,+\infty)\)

0.0

NO (set via \(c\) and \(\sigma_c\))

N/A

sn_sigma

\(\sigma\)

( )

FLT

\([0.0,+\infty)\)

0.0

NO

N/A

sn_esigma

\(\sigma_e\)

Effective normal stress at bond periphery [stress]

FLT

\([0.0,+\infty)\)

0.0

NO

N/A

sn_tau

\(\tau\)

Shear stress at bond periphery [stress]

FLT

\((-\infty,+\infty)\)

0.0

NO

N/A

sn_sdisp

\(\sum \delta _s\)

Shear displacement [length]

VEC2

\(\mathbb{R}^2\)

\(\mathbf{0}\)

NO

N/A

sn_slip

\(s_{sub}\)

Slip state [-]

BOOL

{false,true}

false

NO

N/A

\(\;\;\;\;\;\;\begin{cases} \mbox{true: slipping} \\ \mbox{false: not slipping} \end{cases}\). If sn_index is zero, true is returned if any subcontact is slipping, otherwise the subcontact slip state is returned

sn_force

\(\mathbf{F}_{sub}\)

Force (contact plane coord. system)

VEC

\(\mathbb{R}^3\)

\(\mathbf{0}\)

YES

NO

\(\left( -F_n,F_{ss},F_{st} \right) \quad \left(\mbox{2D model: } F_{ss} \equiv 0 \right)\). If sn_index is zero, the accumulated subcontact force is returned, otherwise the subcontact force is returned

sn_moment

\(\mathbf{M}_{sub}\)

Moment (contact plane coord. system)

VEC

\(\mathbb{R}^3\)

\(\mathbf{0}\)

YES

NO

\(\left( M_t,M_{bs},M_{bt} \right) \quad \left(\mbox{2D model: } M_{t} \equiv M_{bt} \equiv 0 \right)\). If sn_index is zero, the accumulated subcontact moment is returned, otherwise the subcontact moment is returned

sn_non_diag

\(nd\)

Off-diagonal computation mode [-]

FLT

\([0,1]\)

1

YES

NO

sn_pois

\(\nu\)

Poisson ratio [-]

FLT

\([0,+\infty)\)

0.0

YES

NO

sn_pois_force

\(\mathbf{F^{p}}_{sub}\)

Poisson force (contact plane coord. system)

VEC

\(\mathbb{R}^3\)

\(\mathbf{0}\)

NO

NO

If sn_index is zero, the accumulated subcontact fictitious force is returned, otherwise the subcontact fictitious force is returned

sn_porp

\(P_p\)

Pore pressure [stress]

FLT

\((-\infty,+\infty)\)

0.0

YES

NO

Dashpot Group:

dp_nratio

\(\beta_n\)

Normal critical damping ratio [-]

FLT

\([0.0,1.0]\)

0.0

YES

NO

dp_sratio

\(\beta_s\)

Shear critical damping ratio [-]

FLT

\([0.0,1.0]\)

0.0

YES

NO

dp_mode

\(M_d\)

Dashpot mode [-]

INT

{0,1,2,3}

0

YES

NO

\(\;\;\;\;\;\;\begin{cases} \mbox{0: full normal & full shear} \\ \mbox{1: no-tension normal & full shear} \\ \mbox{2: full normal & slip-cut shear} \\ \mbox{3: no-tension normal & slip-cut shear} \end{cases}\)

dp_force

\(\mathbf{F^d}_{sub}\)

Dashpot force (contact plane coord. system)

VEC

\(\mathbb{R}^3\)

\(\mathbf{0}\)

NO

NO

\(\left( -F_n^d,F_{ss}^d,F_{st}^d \right) \quad \left(\mbox{2D model: } F_{ss}^d \equiv 0 \right)\). If sn_index is zero, the accumulated subcontact dashpot force is returned, otherwise the subcontact dashpot force is returned

\(^*\) By convention, \(\kappa^*\) equals zero if either normal or shear stiffness is zero.

Note

Modifying the contact model force will not alter forces accumulated to the bodies.

Therefore, any change to \(\mathbf{F^l}\) or \(\mathbf{M}\) may only be effective during the next force-displacement calculation. When \(M_l = 0\), the normal component of the linear force is automatically overridden during the next force-displacement calculation.

Surface Property Inheritance

The linear stiffnesses, \(k_n\) and \(k_s\), and the friction coefficient, \(\mu\), may be inherited from the contacting pieces. See this section from the linear formulation for details on property inheritance.

Methods

Table 3: Subspring-network Model Methods

Method

Arguments

Symbol

Type

Range

Default

Description

spring-network Group:

area

Set user_area to sn_area

assign-stiffness

Assign stiffness without Poisson correction

kn

\(k_n\)

FLT

\([0.0,+\infty)\)

N/A

Normal stiffness

kratio

\(\kappa^*\)

FLT

\([0.0,+\infty)^*\)

N/A

Normal-to-shear stiffness ratio

bond

Bond the contact if \(g_c \in G\)

gap

\(G\)

VEC2

\(\mathbb{R}^2\)

\((-\infty,0]\)

Gap interval

compute-stiffness

Compute stiffnesses with Poisson correction

emod

\(E\)

FLT

\([0.0,+\infty)\)

N/A

Young’s modulus

poisson

\(\nu\)

FLT

\([0.0,+\infty)^*\)

N/A

Poisson ratio

unbond

Unbond the contact if \(g_c \in G\)

gap

\(G\)

VEC2

\(\mathbb{R}^2\)

\((-\infty,0]\)

Gap interval

\(^*\) By convention, setting \(\kappa^*\) equal to zero sets the shear stiffness to zero but does not modify the normal stiffness.

area

Set the user_area property via the current sn_area. This operation means that the contact area stays constant and is fixed independent of changes to the piece sizes/geometries. In order for the stiffnesses to be recomputed accounting for this area, one should subsequently call the deformability method.

assign-stiffness

Assign normal and shear stiffness without the Poisson’s ratio correction. This should be used for assigning notional joint stiffnesses for non-bonded, pre-existing joints. This computation is outlined above, except one can specify a different shear stiffness and the rotational stiffnesses are 0.

bond

Activate the bond if the contact gap between the pieces is within the bonding-gap interval. If no gap is specified, then the bond is activated if the pieces overlap. A single value can be specified with the gap keyword corresponding to the maximum gap. One can ensure the existence of contacts between all pieces with a contact gap less than a specified bonding gap \((g_b)\) by specifying \(g_b\) with the proximity in the contact cmat default command of the Contact Model Assignment Table (CMAT). If the bond is activated, then the normal force calculation mode sn_mode is automatically set to 1 (incremental).

compute-stiffness

Compute stiffnesses to match the target Young’s modulus and Poisson’s ratio. This computation is outlined above.

unbond

Deactivate the bond if the contact gap between the pieces is within the gap interval. If no gap is specified, then the bond is deactivated if the pieces overlap. A single value can be specified with the gap keyword corresponding to the maximum gap. If the bond is deactivated, then the bond state becomes unbonded \((B=0)\). The force and moment are unaffected and will be updated during the next cycle.

Callback Events

Table 4: Spring Network Model Callback Events

Event

Array Slot

Value Type

Range

Description

contact_activated

Contact has become active

1

C_PNT

N/A

Contact pointer

spring-network Group:

slip_change

Subcontact slip state has changed

1

C_PNT

N/A

Contact pointer

2

INT

{1,2,3}

Subcontact index for use with sn_index

3

INT

{0,1}

slip state

bond_break

Subcontact bond has broken

1

C_PNT

N/A

Contact pointer

2

INT

{1,2}

Failure mode

See equations (2) and (3) of the Bond State section

3

VEC

\(\mathbb{R}^3\)

Vector, in the global coordinate system, from the contact position to the subcontact position

4

INT

{1,2,3}

Subcontact index for use with sn_index

bond_broken

All subcontacts have broken

1

C_PNT

N/A

Contact pointer

2

INT

{2,3,4,5,6}

Accumulated \(B_{sub}\)

See equations (2) and (3) of the Bond State section

Usage and Verification Examples

The example i Rigid Block Model of Tunnel Excavation demonstrates usage of the subspring network contact model for the sequential excavation of a tunnel. Note in that model, sn_cntconv is set to 10.

Model Summary