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 particlescale heterogeneity and provide elasticallyhomogeneous 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 bondedblock models that directly match the macroscale elastic properties of isotropic, transversely isotropic and bimodular materials without the need for timeconsuming 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 RBSNDEM bondedblock 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 softbond 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 softbond model, the bond may not be removed upon failure in tension. Instead, it may enter into a userspecified 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} crosssection 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 softbond 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 nonzero 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 userdefined forcedisplacement regime before breaking. If the maximum shear stress exceeds the bond shear strength, then it breaks and may enter a userdefined forcedisplacement 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.
ActivityDeletion 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 forcedisplacement 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 welldefined 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
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
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 ballball contacts the interface polygon is a disk. In 2D there is only one rotational stiffness. The computestiffness
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 assignstiffness
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 zeroporosity packings of {polygons in 2D; polyhedra in 3D} with coincident facets, this formulation (where contacts exist between all adjacent facets) results in a continuumlike 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 latticebased 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 ballball 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:
These expressions apply to 2D rblockrblock contacts where \(R\) is half the length of the interface line segment.
For 3D rblockrblock 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
Fictitious Stress Correction
Though the above formulation can effectively match the target Young’s modulus and remove particlescale 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 redistribution 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 contactresolution
, ball accumulatestress
, and clump accumulatestress
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 stressstrain response, including Poisson effect, is given by
where \(\epsilon\) denotes strain, \(\nu\) denotes Poisson’s ratio, and \(\sigma\) denotes stress. This is the target stressstrain response. As the formulation above provides stiffnesses corresponding to a zero Poisson’s ratio material, the actual stressstrain response is given by
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
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 offdiagonal 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 postpeak behavior, the spring network model allows one to specify the postpeak forcedisplacement 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 nonlinear postpeak 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.
ForceDisplacement Law
The forcedisplacement law for the spring network model updates the contact force and moment:
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:
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:
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.
Update \(F_n\) based on the normalforce 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 normaldisplacement increment of Equation (12) of the i Contact Resolution section.
Note:
the unbonded behavior cannot sustain tensile normal force,
if either the
computestiffness
orassignstiffness
methods are used then \(M_l=1\), andthere are no fictitious stresses when unbonded.
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 sheardisplacement increment of Equation (12) of the i Contact Resolution section.
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 twistrotation increment of Equation (12) of the i Contact Resolution section.
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 bendrotation increment of Equation (12) of the i Contact Resolution section and \(\mathbf{k_{rs}}\) is the bending stiffness vector.
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, bendingslip state, or twistingslip 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.
Update the dashpot forces:
The dashpot force is updated as in the linear model.
Bonded Behavior
As mentioned above, the computestiffness
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 nonzero 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.
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 normaldisplacement and sheardisplacement 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.
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 bendrotation and twistrotation increments of Equation (12) of the i Contact Resolution section.
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\]
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 springnetwork 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
Keyword 
Symbol 
Description 
Range 
Accumulated 

springnetwork: 


\(E_{k}\) 
strain energy 
\([0.0,+\infty)\) 
NO 

\(E_{\mu}\) 
total energy dissipated by slip 
\((\infty,0.0]\) 
YES 
Dashpot Group: 


\(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.
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).\]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}\]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 this equation of the i Contact Resolution section.
Properties
The properties defined by the springnetwork 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.
Keyword 
Symbol 
Description 
Type 
Range 
Default 
Modifiable 
Inheritable 



SpringNetwork Group: 


\(k_n\) 
Normal stiffness [stress/disp./area] 
FLT 
\([0.0,+\infty)\) 
0.0 
YES 
YES 

\(k_s\) 
Shear stiffness [stress/disp./area] 
FLT 
\([0.0,+\infty)\) 
0.0 
YES 
YES 

\(\mu\) 
Friction coefficient [] 
FLT 
\([0.0,+\infty)\) 
0.0 
YES 
YES 

\(g_r\) 
Reference gap [length] 
FLT 
\(\mathbb{R}\) 
0.0 
YES 
NO 

\(\lambda_b\) 
Bending Friction multiplier [] 
FLT 
\([0.0,+\infty)\) 
1.0 
YES 
YES 

\(\lambda_t\) 
Twisting Friction multiplier [] 
FLT 
\([0.0,+\infty)\) 
1.0 
YES 
YES 

\(M_l\) 
Normalforce update mode [] 
INT 
\(\{0,1\}\) 
0 
YES 
NO 
\(\;\;\;\;\;\;\begin{cases} \mbox{0: update is absolute} \\ \mbox{1: update is incremental} \end{cases}\) 


\(E^*\) 
Effective modulus [force/area] 
FLT 
\([0.0,+\infty)\) 
0.0 
NO 
N/A 

\(\kappa^*\) 
Normaltoshear stiffness ratio [] 
FLT 
\([0.0,+\infty)^*\) 
0.0\(^*\) 
NO 
N/A 
\(\kappa^* \equiv \frac{k_n}{k_s}\) 


\(\lambda\) 
Radius multiplier [] 
FLT 
\((0.0,+\infty)\) 
1.0 
YES 
NO 

\(A\) 
Bond area [length*length] 
FLT 
\((0.0,+\infty)\) 
N/A 
NO 
NO 

\(A\) 
Constant bond area [length*length] 
FLT 
\((0.0,+\infty)\) 
0.0 
YES 
NO 

\(R\) 
Effective radius [length] 
FLT 
\((0.0,+\infty)\) 
N/A 
NO (set via \(\lambda\)) 
NO 

\(\beta\) 
Momentcontribution factor [] 
FLT 
\([0.0,1]\) 
1.0 
YES 
NO 

\(\sigma_c\) 
Tensile strength [stress] 
FLT 
\([0.0,+\infty)\) 
0.0 
YES 
NO 

\({\sigma_c}^{dc}\) 
Critical elongation [length] 
FLT 
\([0.0,+\infty)\) 
0.0 
YES 
NO 

Tensile strength multiplier as a function of elongation 
VEC2 
\(\mathbb{R}^2\) 
YES 
NO 


\(c\) 
Cohesion [stress] 
FLT 
\([0.0,+\infty)\) 
0.0 
YES 
NO 

\(c^{dc}\) 
Critical slip distance [length] 
FLT 
\([0.0,+\infty)\) 
0.0 
YES 
NO 

Cohesive strength multiplier as a function of shear displacement 
VEC2 
\(\mathbb{R}^2\) 
YES 
NO 


\(\phi\) 
Friction angle [degrees] 
FLT 
\([0.0,90.0)\) 
0.0 
YES 
NO 

\(\psi\) 
Dilation angle [degrees] 
FLT 
\([0.0,90.0)\) 
0.0 
YES 
NO 

\(u_{cs}\) 
Critical dilation distance [length] 
FLT 
\([0.0,+\infty)\) 
+infty 
YES 
NO 

\(h\) 
Healing 
INT 
{0,1} 
0 
YES 
NO 

\(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}\) 


\(\tau_{c}\) 
Bond shear strength [stress] 
FLT 
\([0.0,+\infty)\) 
0.0 
NO (set via \(c\) and \(\sigma_c\)) 
N/A 

\(\sigma\) 
Total normal stress at bond periphery [stress] 
FLT 
\([0.0,+\infty)\) 
0.0 
NO 
N/A 

\(\sigma_e\) 
Effective normal stress at bond periphery [stress] 
FLT 
\([0.0,+\infty)\) 
0.0 
NO 
N/A 

\(\tau\) 
Shear stress at bond periphery [stress] 
FLT 
\((\infty,+\infty)\) 
0.0 
NO 
N/A 

\(\sum \delta _n\) 
Normal displacement [length] 
FLT 
\((\infty,+\infty)\) 
0.0 
NO 
N/A 

\(\sum \delta _s\) 
Shear displacement [length] 
VEC2 
\(\mathbb{R}^2\) 
\(\mathbf{0}\) 
NO 
N/A 

\(s\) 
Slip state [] 
BOOL 
{false,true} 
false 
NO 
N/A 
\(\;\;\;\;\;\;\begin{cases} \mbox{true: slipping} \\ \mbox{false: not slipping} \end{cases}\) 


\(sb\) 
Bending slip state [] 
BOOL 
{false,true} 
false 
NO 
N/A 
\(\;\;\;\;\;\;\begin{cases} \mbox{true: slipping} \\ \mbox{false: not slipping} \end{cases}\) 


\(st\) 
Twisting slip state [] 
BOOL 
{false,true} 
false 
NO 
N/A 
\(\;\;\;\;\;\;\begin{cases} \mbox{true: slipping} \\ \mbox{false: not slipping} \end{cases}\) 


\(\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)\) 


\(\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)\) 


\(nd\) 
Offdiagonal computation mode [] 
FLT 
\([0,1]\) 
0 
YES 
NO 

\(\nu\) 
Poisson ratio [] 
FLT 
\([0,+\infty)\) 
0.0 
YES 
NO 

\(\mathbf{F^{p}}\) 
Poisson force (contact plane coord. system) 
VEC 
\(\mathbb{R}^3\) 
\(\mathbf{0}\) 
NO 
NO 

\(P_p\) 
Pore pressure [stress] 
FLT 
\((\infty,+\infty)\) 
0.0 
YES 
NO 
Dashpot Group: 


\(\beta_n\) 
Normal critical damping ratio [] 
FLT 
\([0.0,1.0]\) 
0.0 
YES 
NO 

\(\beta_s\) 
Shear critical damping ratio [] 
FLT 
\([0.0,1.0]\) 
0.0 
YES 
NO 

\(M_d\) 
Dashpot mode [] 
INT 
{0,1,2,3} 
0 
YES 
NO 
\(\;\;\;\;\;\;\begin{cases} \mbox{0: full normal & full shear} \\ \mbox{1: notension normal & full shear} \\ \mbox{2: full normal & slipcut shear} \\ \mbox{3: notension normal & slipcut shear} \end{cases}\) 


\(\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 forcedisplacement calculation. When \(M_l = 0\), the normal component of the linear force is automatically overridden during the next forcedisplacement 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
Method 
Arguments 
Symbol 
Type 
Range 
Default 
Description 

springnetwork Group: 


Set 


Assign stiffness without Poisson correction 


\(k_n\) 
FLT 
\([0.0,+\infty)\) 
N/A 
Normal stiffness 


\(\kappa^*\) 
FLT 
\([0.0,+\infty)^*\) 
N/A 
Normaltoshear stiffness ratio 


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


\(G\) 
VEC2 
\(\mathbb{R}^2\) 
\((\infty,0]\) 
Gap interval 


Compute stiffnesses with Poisson correction 


\(E\) 
FLT 
\([0.0,+\infty)\) 
N/A 
Young’s modulus 


\(\nu\) 
FLT 
\([0.0,+\infty)^*\) 
N/A 
Poisson ratio 


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


\(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.
assignstiffness
Assign normal and shear stiffness without the Poisson’s ratio correction. This should be used for assigning notional joint stiffnesses for nonbonded, preexisting 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 bondinggap 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).
computestiffness
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
Event 
Array Slot 
Value Type 
Range 
Description 


Contact has become active 

1 
C_PNT 
N/A 
Contact pointer 

springnetwork Group: 


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 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.
References
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.
Asahina, D., K. Aoyagi, K. Kim, J. T. Birkholzer and J. E. Bolander. “Elasticallyhomogeneous lattice models of damage in geomaterials,” Computers and Geotechnics, 81, 195206, 2017.
Bolander, J. E. and S. Sato. “Fracture analyses using spring networks with random geometry,” Engineering Fracture Mechanics, 61, 569–591, 1998.
Kawai, T. “New discrete models and their application to seismic response analysis of structures,” Nuclear Engineering and Design, 48, 207–229, 1978.
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, 3141, 2018.
Rasmussen, L. L. and A. P. de Assis. “Elasticallyhomogeneous lattice modeling of transversely isotropic rocks,” Computers and Geotechnics, 104, 96108, 2018.
Rasmussen, L. L. and M. M. de Farias. “Lattice modelling of gravity and stressdriven failure of rock tunnels,” Computers and Geotechnics, 116, 103183, 2019.
Rasmussen, L. L. “Hybrid lattice/discrete element method for bonded block modeling of rocks,” Computers and Geotechnics, 130, 103907, 2021.
Was this helpful? ...  Itasca Software © 2024, Itasca  Updated: Apr 02, 2024 