Spring Network Model

The spring network model is referred to in commands and FISH by the name springnetwork.

Introduction

The spring network model can be used to simulate unbonded and bonded systems. 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.

In an unbonded state, the spring network model behaves like the contact model proposed by [Jiang2015a], providing the ability to transmit both a force and a moment at the contact point, with frictional strength parameters limiting the shear force, bending moment, and twisting moment. Dilation and evolving shear strength as a function of slip can also be specified.

In a bonded state, the spring network model behaves like a linear parallel bond or soft-bond model, with the possibility for the bond to fail, if the bond strength is exceeded either in tension or in shear. Unlike those models, contact properties can be computed based on lattice theory to provide a theoretical match to a target macroscopic Young’s modulus and Poisson’s ratio. Similar to the soft-bond model, the bond may not be removed upon failure in tension. Instead, it may enter into a user-specified softening regime until the bond stress reaches a threshold value, at which time the bond is removed and considered broken.

Behavior Summary

A spring network bond combines two concepts. 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 contact point. Relative motion at the contact causes development of a linear force and moment that act on the two contacting pieces. This behavior is the same as the linear parallel bond and the soft-bond model. But unlike those models, a method exists 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. If the maximum normal stress exceeds the bond tensile strength, the bond may enter a user-defined force-displacement regime before breaking. If the maximum shear stress exceeds the bond shear strength, then it breaks and may enter a user-defined force-displacement regime in the shear direction. If the bond is inactive but the contact is active, frictional strength parameters govern slip with potential dilation and/or healing.

Activity-Deletion Criteria

A contact with the spring network model is active if it is bonded or if the surface gap is less than or equal to zero. The force-displacement law is skipped for inactive contacts.

Contact Stiffness Computation

The spring network model relies on lattice theory to compute stiffnesses for bonded assemblies. As outlined in [Bolander1998a], one can compute stiffnesses consistent with a target Young’s modulus and zero Poisson’s ratio elastic material. Though their formulation focuses specifically on Voronoi tessellations, it is applied here to all contact types and for all types of tessellations.

Suppose that a contact exists between two pieces, delineating a well-defined interface surface that is {a line segment in 2D; a planar polygon in 3D}. [Bolander1998a] computes a single translational stiffness, applied to all components of the incremental displacement in the contact local $(nst)$ coordinate system, as

(1) $$k_n = k_s = k_{n_c} = k_{s_c} = k_{t_c} = E A / d$$

where $E$ is the target Young’s modulus, $A$ is the interface area, and $d$ is the distance from the body centroids. The subscript designates the stiffness in the local coordinate system axes. The rotational stiffnesses are given by

(2) $$\begin{split} \begin{array}{rcl} k_{r_{n_c}} = E J / d \\ k_{r_{s_c}} = E I_{s} / d \\ k_{r_{t_c}} = E I_{t} / d \\ \end{array}\end{split}$$

where $J$ is the polar moment of inertia of the interface surface and $I_{s}$, $I_{t}$ are the two principal moments of inertia of the interface surface in 3D. In 3D, discussion below defines $\mathbf{k_{rs}} = (0,k_{r_{s_c}},k_{r_{t_c}})$ as the bending stiffness vector. Note that for 3D ball-ball contacts the interface polygon is a disk. In 2D there is only one rotational stiffness. The compute-stiffness method computes the stiffnesses as described above by setting the user_area property to the corresponding contact area; the area does not change once this method is called, regardless of changes in the piece positions/orientations. The assign-stiffness method does a similar computation, except that the Poisson’s ratio cannot be specified—though one can explicitly set the ratio of the normal to shear stiffnesses. That method is primarily meant to apply to preexisting joint surfaces.

For zero-porosity packings of {polygons in 2D; polyhedra in 3D} with coincident facets, this formulation (where contacts exist between all adjacent facets) results in a continuum-like macroscopic response corresponding to an isotropic, elastic material with zero Poisson’s ratio with the specified Young’s modulus. Importantly, detailed inspection of displacements/stresses at the particle level demonstrate that, in spite of particle size variations within the packing, the elastic response is homogeneous. In other words, using this lattice-based approach to compute stiffnesses removes elastic heterogeneity at the particle scale, allowing the user to explicitly specify this heterogeneity. For ball packings, the same effect can be achieved if the contact areas are scaled consistent with the volume they enclose. One way of achieving this scaling is by using the facet areas of a Voronoi tessellation of the ball packing to set the contact areas, where each Voronoi cell corresponds to a ball (see the ball tractions command for more details).

For ball-ball contacts, take $R$ to be the minimum ball radius or the {length of the interface line segment in 2D; radius of the equivalent disk in 3D} if the user_area property has been set. In this case:

(3) $$\begin{split} \begin{array}{l} A = \left\{ \begin{array}{rl} 2 R , & \mbox{2D} \\ \pi R^2, & \mbox{3D} \end{array} \right. \\[3mm] I_s = I_t = \left\{ \begin{array}{rl} \tfrac{2}{3} R^3, & \mbox{2D} \\ \tfrac{1}{4} \pi R^4, & \mbox{3D} \end{array} \right. \\[3mm] J = \left\{ \begin{array}{rl} 0, & \mbox{2D} \\ \tfrac{1}{2} \pi R^4, & \mbox{3D} \end{array} \right. \end{array}\end{split}$$

These expressions apply to 2D rblock-rblock contacts where $R$ is half the length of the interface line segment.

For 3D rblock-rblock contacts, denote the $N$ interface surface vertices as $(s_0,t_0)$, $(s_1,t_1)$, … $(s_N,t_N)$ that are rotated onto the the $(st)$ plane with coordinate directions in the $(st)$ axes directions. The inertial values can be computed as

(4) $$\begin{split} \begin{array}{l} I_s = {{1}\over {12}}\sum\limits_{N} (s_i t_{i+1} - s_{i+1} t_i)({t_i}^2 + t_i t_{i+1} + {t_{i+1}}^2) \\[3mm] I_t = {{1}\over {12}}\sum\limits_{N} (s_i t_{i+1} - s_{i+1} t_i)({s_i}^2 + s_i s_{i+1} + {s_{i+1}}^2) \\[3mm] J = I_s + I_t \end{array}\end{split}$$

Fictitious Stress Correction

Though the above formulation can effectively match the target Young’s modulus and remove particle-scale elastic heterogeneity, it does not allow for a realistic representation of the Poisson effect (i.e., wherein an object expands or contracts perpendicular to the direction of loading). The Poisson effect may have significant impacts on stress re-distribution in rock masses around discontinuities, or when spatially heterogeneous stresses are applied. To overcome the lack of Poisson effect in RBSN models, [Asahina2015] and [Asahina2017] introduced the fictitious stress approach. Suppose that the stress of each body is computed and stored from the previous timestep. (As outlined here, the stress in a body is computed via the sum of the outer products of the contact forces and branch vectors for all contacts with the body. See the rblock contact-resolution, ball accumulate-stress, and clump accumulate-stress commands to activate stress computation/storage.) The volumetric average of the body stresses is used to define the stress associated with the contact. Volumetric averaging of the stresses is important when the sizes of neighboring particles differ substantially. From linear elasticity theory, the stress-strain response, including Poisson effect, is given by

(5) $$\begin{split}\begin{bmatrix} {\epsilon}_{11} \\ {\epsilon}_{22} \\ {\epsilon}_{33} \\ {\epsilon}_{12} \\ {\epsilon}_{13} \\ {\epsilon}_{23} \end{bmatrix} = \begin{bmatrix} {{1}\over{E}} & -{{\nu}\over{E}} & -{{\nu}\over{E}}& & & \\ -{{\nu}\over{E}} & {{1}\over{E}} & -{{\nu}\over{E}}& & & \\ -{{\nu}\over{E}}& -{{\nu}\over{E}}& {{1}\over{E}} & & & \\ & & & {{1+\nu}\over{E}} & & \\ & & & & {{1+\nu}\over{E}} & \\ & & & & & {{1+\nu}\over{E}} \end{bmatrix} \begin{bmatrix} {\sigma}_{11} \\ {\sigma}_{22} \\ {\sigma}_{33} \\ {\sigma}_{12} \\ {\sigma}_{13} \\ {\sigma}_{23} \end{bmatrix}\end{split}$$

where $\epsilon$ denotes strain, $\nu$ denotes Poisson’s ratio, and $\sigma$ denotes stress. This is the target stress-strain response. As the formulation above provides stiffnesses corresponding to a zero Poisson’s ratio material, the actual stress-strain response is given by

(6) $$\begin{split}\begin{bmatrix} {\epsilon}_{11} \\ {\epsilon}_{22} \\ {\epsilon}_{33} \\ {\epsilon}_{12} \\ {\epsilon}_{13} \\ {\epsilon}_{23} \end{bmatrix} = \begin{bmatrix} {{1}\over{E}} & & & & & \\ & {{1}\over{E}} & & & & \\ & & {{1}\over{E}} & & & \\ & & & {{1}\over{E}} & & \\ & & & & {{1}\over{E}} & \\ & & & & & {{1}\over{E}} \end{bmatrix} \begin{bmatrix} {\sigma}_{11} \\ {\sigma}_{22} \\ {\sigma}_{33} \\ {\sigma}_{12} \\ {\sigma}_{13} \\ {\sigma}_{23} \end{bmatrix}\end{split}$$

In order to reach the target response, the difference in the previous two equations produces the required strain at the contact to account for the Poisson effect. This difference produces the fictitious stresses

(7) $$\begin{split}\begin{bmatrix} {{\sigma}_{11}}^{f} \\ {{\sigma}_{22}}^{f} \\ {{\sigma}_{33}}^{f} \\ {{\sigma}_{12}}^{f} \\ {{\sigma}_{13}}^{f} \\ {{\sigma}_{23}}^{f} \end{bmatrix} = \begin{bmatrix} {-\nu \sigma_{22} -\nu \sigma_{33}} \\ {-\nu \sigma_{11} -\nu \sigma_{33}} \\ {-\nu \sigma_{11} -\nu \sigma_{22}} \\ {\nu \sigma_{12}} \\ {\nu \sigma_{13}} \\ {\nu \sigma_{23}} \end{bmatrix}\end{split}$$

The fictitious stress tensor, computed per contact, is converted to a fictitious force ($\mathbf{F^{f}}$) in the local $(nst)$ axis system by multiplication with the contact area and each axis of the local system. By default, the off-diagonal terms of the stresses that produce shearing strains (i.e., ${\nu \sigma_{12}}$, ${\nu \sigma_{13}}$, and ${\nu \sigma_{23}}$) are not used to compute the fictitious forces because the shear strain increments are found to contribute negligibly to the macro Poisson effect. In addition, for very high aspect ratio rigid blocks, the shear strain increments of the fictitious force can feedback with the angular rotation, leading to instability. This behavior can be changed with the sn_non_diag property. [Rasmussen2018b] demonstrates how to extend this formulation to transversely isotropic elastic materials, though this is not part of the current implementation.

Tensile and Cohesion Tables

In order to accommodate more complex post-peak behavior, the spring network model allows one to specify the post-peak force-displacement response in tension and in shear. The peak strengths are specified with the sn_ten and sn_coh properties. In the simplest form, one may specify critical distances over which the tensile/cohesive stresses degrade linearly from peak values to residual values using the sn_tendc and sn_cohdc properties for the critical elongation in tension and critical slip distance in shear, respectively. A residual cohesive strength may also be specified with the sn_cohres property. For a non-linear post-peak response in tension/shear, one can add entries with varying elongation/slip using the sn_tentab and sn_cohtab properties. The values entered into these tables are scaling factors of the peak tensile/shear strength at failure, and a linear relationship between scaling factor and elongation/displacement is assumed between each entry.

Force-Displacement Law

The force-displacement law for the spring network model updates the contact force and moment:

(8) $$\mathbf{F_{c}} =\mathbf{F^{l}} +\mathbf{F^{d}} +\mathbf{F^{f}} ,\quad \mathbf{M_{c}} =\mathbf{M}$$

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

The elastic force is resolved into a normal and shear force, and the linear moment is resolved into a twisting and bending moment:

(9) $$\begin{split}\begin{array}{rcl} {\mathbf{F}} & {=} & {- F_n {\kern 1pt} \hat{\mathbf{n}}_\mathbf{c} +\mathbf{F_s}} \\[3mm] {\mathbf{M} } & {=} & {M_t {\kern 1pt} \hat{\mathbf{n}}_\mathbf{c} +\mathbf{M_b} \; \; \left({\rm 2D\; model:\; } M_t \equiv 0\right).} \end{array}\end{split}$$

where $F_n > 0$ is tension. The shear force and bending moment lie on the contact plane and are expressed in the contact plane coordinate system:

(10) $$\begin{split}\begin{array}{rcl} {\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)} \\[3mm] {\mathbf{M_b} } & {=} & {M_bs {\kern 1pt} \hat{\mathbf{s}}_\mathbf{c} +M_{bt} {\kern 1pt} \hat{\mathbf{t}}_\mathbf{c} \; \; \left({\rm 2D\; model:\; } M_{bt} \equiv 0\right).} \end{array}\end{split}$$

The forces and moment are updated as described below for the unbonded or bonded spring network contact model.

Unbonded Behavior

When unbonded, the fictitious force $\mathbf{F^f}$ is zero. In the equations below in 3D, $I$ and $J$ are computed using the radius of a disk with an area equivalent to the contact area regardless of contact type. The force and moment are updated with the following steps.

1. Update $F_n$ based on the normal-force update mode $M_l$:

   If $M_l=0$ (absolute update), the normal force is computed as:

   (11) $$\begin{split}F_n = \left\{ \begin{array}{rl} k_n A g_s, & g_s < 0 \\ 0 , & \mbox{otherwise} \end{array} \right.\end{split}$$

   where $g_s$ is the surface gap (opposite of the overlap).

   If $M_l=1$ (incremental update), the normal force is computed as:

   (12) $$F_n := \max\left(F_n + k_n {\kern 1pt} A \, \Delta \delta _n, 0\right)$$

   where $\Delta \delta _n$ is the relative normal-displacement increment of Equation (12) of the i Contact Resolution section.

   Note:

   1. the unbonded behavior cannot sustain tensile normal force,

   2. if either the compute-stiffness or assign-stiffness methods are used then $M_l=1$, and

   3. there are no fictitious stresses when unbonded.

2. Update $\mathbf{F_s}$:

   The shear force $\mathbf{F_s}$ is first updated with:

   (13) $$\mathbf{F_s} := \mathbf{F_s} - k_s {\kern 1pt} A{\kern 1pt} \Delta \pmb{δ} _\mathbf{s}$$

   where $\Delta \pmb{δ}_\mathbf{s}$ is the relative shear-displacement increment of Equation (12) of the i Contact Resolution section.

3. Update $M_t$:

   (14) $$\begin{split}M_t := \left\{ \begin{array}{rl} 0, & \mbox{2D} \\ M_t - k_{r_{n_c}} \Delta \theta_t, & \mbox{3D} \end{array} \right.\end{split}$$

   where $\Delta \theta _{t}$ is the relative twist-rotation increment of Equation (12) of the i Contact Resolution section.

4. Update $\mathbf{M_b}$:

   (15) $$\mathbf{M_b} := \mathbf{M_b} - \mathbf{k_{rs}} {\kern 1pt} \Delta \pmb{θ}_\mathbf{b}$$

   where $\Delta \pmb{θ}_\mathbf{b}$ is the relative bend-rotation increment of Equation (12) of the i Contact Resolution section and $\mathbf{k_{rs}}$ is the bending stiffness vector.

5. Enforce the 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 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 slip state is updated as:

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

   If the slip state is true, then the contact is 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 contact is in a bonded state where it cannot sustain tension but can sustain shear.

   If slipping, 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 in the form

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

   where $\psi$ is the dilation angle in degrees.

   Similar to the Coulomb friction criteria, the maximum bending and twisting moment criteria are enforced as:

   (19) $$\begin{split} \begin{array}{l} \mathbf{M_b} = \left\{ \begin{array}{rl} \mathbf{M_b} , & \| \mathbf{M_b} \| \le M_b^*\\ M_b^* \bigl( \mathbf{M_b} / \| \mathbf{M_b} \| \bigr), & \mbox{otherwise.} \end{array} \right. \\ M_t = \left\{ \begin{array}{rl} M_t , & \|M_t\| \le M_t^* \\ M_t^* \bigl( M_t / \|M_t\| \bigr), & \mbox{otherwise.} \end{array} \right. \end{array}\end{split}$$

   where the critical bending and twisting moments are:

   (20) $$\begin{split}\begin{array}{l} M_b^* & = - 2.1 \lambda_b F_n R {\kern 1pt} / {\kern 1pt} 4 \\ M_t^* & = - 0.65 \lambda_t \mu F_n * R \end{array}\end{split}$$

   where $\lambda_b$ and $\lambda_t$ are, respectively, the bending and twisting friction multiplier, that default to 1.0. The critical bending and twisting moments with $\lambda_b=1.0$ and $\lambda_t=1.0$ correspond to the values proposed by [Jiang2015a].

   The bending and twisting slip states are updated as:

   (21) $$\begin{split}\begin{array}{rl} s_b = \left\{ \begin{array}{rl} \mbox{true}, & \| \mathbf{M_b} \| = M_b^* \\ \mbox{false}, & \mbox{otherwise.} \end{array} \right. \\ s_t = \left\{ \begin{array}{rl} \mbox{true}, & \| \mathbf{M_t} \| = M_t^* \\ \mbox{false}, & \mbox{otherwise.} \end{array} \right. \\ \end{array}\end{split}$$

   Whenever the slip state, bending-slip state, or twisting-slip state change, the slip_change callback event occurs. The first argument is the contact pointer, the second argument is an integer ranging from 1 to 7 denoting which slip state(s) has(have) changed, and 3 additional arguments hold the current values of the shear, bend, and twist slip states.

6. Update the dashpot forces:

   The dashpot force is updated as in the linear model.

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. Tensile softening can be accommodated either by: 1) specifying the critical elongation ${\sigma_c}^{dc}$ over which the tensile stress evolves linearly from $\sigma_c$ to 0 as a function of elongation, or 2) specifying a tensile stress table where the tensile stress is interpolated linearly from table entries as a function of elongation.

When bonded, the forces and moment are updated with the following steps.

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

   (22) $$\begin{split}\begin{array}{rl} F_n & := F_n + k_n {\kern 1pt} A \Delta \delta _n + {F_n}^f \\ \mathbf{F_s} & := \mathbf{F_s} - k_s {\kern 1pt} A{\kern 1pt} \Delta \pmb{δ} _\mathbf{s} + \mathbf{{F_s}^f} \end{array}\end{split}$$

   where $\Delta \delta _n$ and $\Delta \pmb{δ}_\mathbf{s}$ are, respectively, the relative normal-displacement and shear-displacement increments of Equation (12) of the i Contact Resolution section, and ${F_n}^f$/ $\mathbf{{F_s}^f}$ are the normal and shear fictitious forces, respectively, described above.

2. Update $\mathbf{M_b}$ and $M_t$:

   (23) $$\begin{split}\begin{array}{rl} \mathbf{M_b} & := \mathbf{M_b} - k_{r_{n_c}} {\kern 1pt} \Delta \pmb{θ}_\mathbf{b} \\ M_t & := \left\{ \begin{array}{rl} 0, & \mbox{2D} \\ M_t - \mathbf{k_{rs}} \Delta \theta_t, & \mbox{3D} \end{array} \right. \end{array}\end{split}$$

   where $\Delta \pmb{θ}_\mathbf{b}$ and $\Delta \theta _{t}$ are, respectively, the relative bend-rotation and twist-rotation increments of Equation (12) of the i Contact Resolution section.

3. Update the maximum normal ($\sigma$, $\sigma > 0$ is tension) and shear ($\tau$) stresses at the bond periphery:

   (24) $$\begin{split}\begin{array}{l} \sigma = \displaystyle \frac{F_n} {A} + \beta \frac{\left\| \mathbf{M_b} \right\| R}{I} \\ \tau = \frac{\left\| \mathbf{F_s} \right\|}{A} + \left\{ \begin{array}{rl} 0, & \mbox{2D} \\ \beta \frac{\left| M_t \right| R}{J}, & \mbox{3D} \end{array} \right. \\ \mbox{with }\quad \beta \in [0, 1]. \end{array}\end{split}$$

   The terms $I$ and $J$ are computed for an equivalent disk.

   The effective normal stress at the bond periphery, $\sigma_e$, is used for subsequent failure computations, using the pore pressure $P_p$:

   (25) $$\sigma_e = \sigma - P_p$$

3. Update the bond state:

   Take the shear strength $\tau_c = c - \sigma_c tan \phi$.

   • If the bond has healed ($B=6$) 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=1$). If $\tau > \tau_c$ then the bond breaks in shear ($B=2$) and may enter a cohesive strength softening regime.

   • If the bond is intact ($B=3$) and the effective normal stress at the bond periphery exceeds the bond tensile strength ($\sigma_e > \sigma_c$), then the bond enters the softening regime if either a critical elongation ${\sigma_c}^{dc}$ or a tensile strength table is specified ($B=4$); otherwise the bond breaks in tension ($B=1$).

   • If the bond is intact ($B=3$) and the maximum shear stress at the bond periphery exceeds the shear strength ($\tau > \tau_c$), then the bond fails in shear ($B=2$) and may enter a cohesive strength softening regime.

   • If the bond is in the softening regime ($B=4$), take $\sigma^*$ as the interpolated tensile strength based on the current elongation. The effective normal stress at bond periphery is checked against the interpolated tensile strength. If $\sigma_e \ge \sigma^*$ both the normal force and bending moment are reduced to the softening envelope

     (26) $$\begin{split}\begin{array}{l} F_n & := F_n {\kern 1pt} \left(\sigma^* / \sigma_e \right) \\ \mathbf{M_b} & := \mathbf{M_b} {\kern 1pt} \left(\sigma^* / \sigma_e \right) \end{array}\end{split}$$

   The bond remains in the softening regime ($B = 4$) if this is the case. If $\sigma_e \lt \sigma^*$ then $B = 5$ and the bond enters the softening regime with compression.

   • If the bond is in the softening regime and in compression ($B=5$) and the maximum normal stress at bond periphery exceeds the bond tensile strength at which compression initiated, then the bond returns to the softening regime ($B=4$) and the tensile force and bending moment are scaled to project the tensile stress on the softening envelope as discussed above.

Energy Partitions

The spring-network model provides three energy partitions:

• 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:

   (27) $${\rm E} _{k} =\frac{1}{2} \left( \frac{F_n^{2} }{k_{n}{\kern 1pt} A} + \frac{\left\| \mathbf{F_s} \right\| ^{2} }{k_{s} {\kern 1pt} A} + \frac{\left\| \mathbf{M_b} \right\| ^{2} }{k_n {\kern 1pt} I} + \frac{M_t^{2} }{k_s {\kern 1pt} J} \right).$$

2. Update the slip energy:

   (28) $$\begin{split}\begin{array}{l} E_{\mu } := E_{\mu} + \Delta E_{\mu}^s + \Delta E_{\mu}^b + \Delta E_{\mu}^t \\ {\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}^{\mu } \\ \Delta E_{\mu}^b := - {\tfrac{1}{2}} \left(\left(\mathbf{M_{b}} \right)_{o} +\mathbf{M_{b}} \right)\cdot \Delta \pmb{\theta} _\mathbf{b}^{\mu } \\ \Delta E_{\mu}^t := - {\tfrac{1}{2}} \left(\left(M_{t} \right)_{o} + M_{t} \right)\cdot \Delta \pmb{\theta} _\mathbf{t}^{\mu } \\ \end{array} \\ {\rm with} \\ \qquad \begin{array}{l} \Delta \pmb{δ} _\mathbf{s}^{\mu } =\Delta \pmb{δ} _\mathbf{s} -\Delta \pmb{δ} _\mathbf{s}^{k} =\Delta \pmb{δ} _\mathbf{s} -\left(\frac{\mathbf{F_{s}} -\left(\mathbf{F_{s}} \right)_{o} }{k_{s} {\kern 1pt} A} \right) \\ \Delta \pmb{\theta} _\mathbf{b}^{\mu } =\Delta \pmb{\theta} _\mathbf{b} -\Delta \pmb{\theta} _\mathbf{b}^{k} =\Delta \pmb{\theta} _\mathbf{b} -\left(\frac{\mathbf{M_{b}} -\left(\mathbf{M_{b}} \right)_{o} }{k_{n} {\kern 1pt} I} \right) \\ \Delta \pmb{\theta} _\mathbf{t}^{\mu } =\Delta \pmb{\theta} _\mathbf{t} -\Delta \pmb{\theta} _\mathbf{t}^{k} =\Delta \pmb{\theta} _\mathbf{t} -\left(\frac{M_{t} -\left(M_{t} \right)_{o} }{k_{s} {\kern 1pt} J} \right) \\ \end{array} \end{array}\end{split}$$

3. Update the dashpot energy:

   (29) $$E_{\beta } :=E_{\beta } - \mathbf{F^{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
sn_bmul	$\lambda_b$	Bending Friction multiplier [-]	FLT	$[0.0,+\infty)$	1.0	YES	YES
sn_tmul	$\lambda_t$	Twisting Friction multiplier [-]	FLT	$[0.0,+\infty)$	1.0	YES	YES
sn_mode	$M_l$	Normal-force update mode [-]	INT	$\{0,1\}$	0	YES	NO
		$\;\;\;\;\;\;\begin{cases} \mbox{0: update is absolute} \\ \mbox{1: update is incremental} \end{cases}$
emod	$E^*$	Effective modulus [force/area]	FLT	$[0.0,+\infty)$	0.0	NO	N/A
kratio	$\kappa^*$	Normal-to-shear stiffness ratio [-]	FLT	$[0.0,+\infty)^*$	0.0$^*$	NO	N/A
		$\kappa^* \equiv \frac{k_n}{k_s}$
sn_rmul	$\lambda$	Radius multiplier [-]	FLT	$(0.0,+\infty)$	1.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_radius	$R$	Effective radius [length]	FLT	$(0.0,+\infty)$	N/A	NO (set via $\lambda$)	NO
sn_mcf	$\beta$	Moment-contribution factor [-]	FLT	$[0.0,1]$	1.0	YES	NO
sn_ten	$\sigma_c$	Tensile strength [stress]	FLT	$[0.0,+\infty)$	0.0	YES	NO
sn_tendc	${\sigma_c}^{dc}$	Critical elongation [length]	FLT	$[0.0,+\infty)$	0.0	YES	NO
sn_tentab		Tensile strength multiplier as a function of elongation	VEC2	$\mathbb{R}^2$		YES	NO
sn_coh	$c$	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_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_heal	$h$	Healing	INT	{0,1}	0	YES	NO
sn_state	$B$	Bond state	INT	{0,1,2,3,4,5,6}	0	NO	NO
		$\;\;\;\;\;\;\begin{cases} \mbox{0: unbonded} \\ \mbox{1: unbonded & broke in tension} \\ \mbox{2: unbonded & broke in shear} \\ \mbox{3: bonded} \\ \mbox{4: bonded & softening} \\ \mbox{5: bonded & compressing} \\ \mbox{6: healed, can sustain shear but not tension}\end{cases}$
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$	Total normal stress at bond periphery [stress]	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_ndisp	$\sum \delta _n$	Normal displacement [length]	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$	Slip state [-]	BOOL	{false,true}	false	NO	N/A
		$\;\;\;\;\;\;\begin{cases} \mbox{true: slipping} \\ \mbox{false: not slipping} \end{cases}$
sn_slipb	$sb$	Bending slip state [-]	BOOL	{false,true}	false	NO	N/A
		$\;\;\;\;\;\;\begin{cases} \mbox{true: slipping} \\ \mbox{false: not slipping} \end{cases}$
sn_slipt	$st$	Twisting slip state [-]	BOOL	{false,true}	false	NO	N/A
		$\;\;\;\;\;\;\begin{cases} \mbox{true: slipping} \\ \mbox{false: not slipping} \end{cases}$
sn_force	$\mathbf{F}$	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)$
sn_moment	$\mathbf{M}$	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)$
sn_non_diag	$nd$	Off-diagonal computation mode [-]	FLT	$[0,1]$	0	YES	NO
sn_pois	$\nu$	Poisson ratio [-]	FLT	$[0,+\infty)$	0.0	YES	NO
sn_pois_force	$\mathbf{F^{p}}$	Poisson force (contact plane coord. system)	VEC	$\mathbb{R}^3$	$\mathbf{0}$	NO	NO
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}$	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)$
$^*$ 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: spring-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				Slip state has changed
	1	C_PNT	N/A	Contact pointer
	2	INT	{0,7}	Slip change mode
				$\;\;\;\;\;\;\begin{cases} \mbox{1: shear slip state has changed} \\ \mbox{2: bend slip state has changed} \\ \mbox{3: twist slip state has changed} \\ \mbox{4: shear & bend slip states have changed} \\ \mbox{5: shear & twist slip states have changed} \\ \mbox{6: bend & twist slip states have changed} \\ \mbox{7: all slip states have changed} \end{cases}$
	3	INT	{0,1}	shear slip state
	4	INT	{0,1}	twist slip state
	5	INT	{0,1}	bend slip state
bond_break				Bond has broken
	1	C_PNT	N/A	Contact pointer
	2	INT	{1,2}	Failure mode
				$\;\;\;\;\;\;\begin{cases} \mbox{1: failed in tension} \\ \mbox{2: failed in shear} \end{cases}$
	3	FLT	$[0.0,+\infty)$	Failure strength [force] ($\overline{\sigma }_{c}$ or $\overline{\tau }_{c}$, according to the failure mode)
	4	FLT	$[0.0,+\infty)$	Bond strain energy $\overline{E}_k$ at onset of failure

Usage and Verification Examples

The verification example i Spring Network Contact Model Capabilities demonstrates the different capabilities the spring network contact model.

The example i Using the Rigid Body Spring Network Paradigm demonstrates usage of the spring network contact model for laboratory tests (Direct Tension, Unconfined Compression and Triaxial tests).

Model Summary

An alphabetical list of the spring network contact model methods is given here. An alphabetical list of the spring network contact model properties is given here.

• Spring Network Model Methods
• Spring Network Model Properties

References

[Asahina2015] (1,2)

Asahina, D., K. Ito, J. E. Houseworth, J. T. Birkholzer and J. E. Bolander. “Simulating the Poisson Effect in Lattice Models of Elastic Continua,” Computers and Geotechnics, 70, 60–67, 2015.

[Asahina2017] (1,2)

Asahina, D., K. Aoyagi, K. Kim, J. T. Birkholzer and J. E. Bolander. “Elastically-homogeneous lattice models of damage in geomaterials,” Computers and Geotechnics, 81, 195-206, 2017.

[Bolander1998]

Bolander, J. E. and S. Sato. “Fracture analyses using spring networks with random geometry,” Engineering Fracture Mechanics, 61, 569–591, 1998.

[Kawai1978]

Kawai, T. “New discrete models and their application to seismic response analysis of structures,” Nuclear Engineering and Design, 48, 207–229, 1978.

[Rasmussen2018a]

Rasmussen, L. L., M. M. de Farias and A. P. de Assis. “Extended Rigid Body Spring Network method for the simulation of brittle rocks,” Computers and Geotechnics, 99, 31-41, 2018.

[Rasmussen2018b] (1,2,3)

Rasmussen, L. L. and A. P. de Assis. “Elastically-homogeneous lattice modeling of transversely isotropic rocks,” Computers and Geotechnics, 104, 96-108, 2018.

[Rasmussen2019]

Rasmussen, L. L. and M. M. de Farias. “Lattice modelling of gravity and stress-driven failure of rock tunnels,” Computers and Geotechnics, 116, 103183, 2019.

[Rasmussen2021b]

Rasmussen, L. L. “Hybrid lattice/discrete element method for bonded block modeling of rocks,” Computers and Geotechnics, 130, 103907, 2021.

[Rasmussen2022] (1,2)

Rasmussen, L. L. “A breakable grain-based model for bi-modular rocks,” in International Journal of Rock Mechanics and Mining Sciences, 151, 105028, 2022.