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 [Bolander1998a] based on the work of [Kawai1978], and the subsequent work of [Asahina2015]. The formulation provides a spring model that matches a target Young’s modulus directly, without calibration. [Rasmussen2021] was the first to adapt the model to the explicit DEM scheme and named it the Hybrid Lattice / Discrete Element Method or Hybrid LDEM. In addition, [Rasmussen2021] demonstrated the effective use of this methodology and its advantages when applied to DEM simulations of rock materials, including the effective removal of particle-scale heterogeneity in elastic response, allowing the user to control heterogeneity. [Asahina2017], [Rasmussen2018] and [Rasmusssen2022] show the ability of both RBSN and Hybrid LDEM to match a target Poisson’s ratio for isotropic elastic materials, all five independent elastic parameters of transversely isotropic materials, and both tensile and compressive young moduli and Poisson’s ratios for bi-modular elastic materials by means of a fictitious stress scheme applied to the inter-particles’ bonds, which overcomes the need for trial-and-error calibrations.
In an unbonded state, the spring network model behaves similar to 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 its bonded formulation, the model behavior is similar to that of a linear parallel bond or soft-bond model, with the possibility for the bond to fail if the bond strength is exceeded either in shear or in tension. 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
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 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:
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
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
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
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 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. [Rasmussen2018] 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:
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 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 this equation of the i Contact Resolution section.
Note:
- the unbonded behavior cannot sustain tensile normal force,
- if either the
compute-stiffnessorassign-stiffnessmethods are used then \(M_l=1\), and - there 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 shear-displacement increment of this equation 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 twist-rotation increment of this equation 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 bend-rotation increment of this equation 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 ) \Delta \pmb{δ}_\mathbf{s}\]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.
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.
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 this equation 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 bend-rotation and twist-rotation increments of this equation 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 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;
- dashpot energy, \(E_{\beta }\), defined as the total energy dissipated by the dashpots; and
| 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.
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 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.
| 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
| 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
| 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.
References
| [Asahina2015] | (1, 2) Asahina, D., K. Ito, J. Houseworth, J. Birkholzer and J. Bolander. “Simulating the Poisson Effect in Lattice Models of Elastic Continua” in Computers and Geotechnics, 70, pp. 60–67, 2015. |
| [Asahina2017] | (1, 2) Asahina, D., K. Aoyagi, K. Kim, J. Birkholzer and J. Bolander. “Elastically-homogeneous lattice models of damage in geomaterials” in Computers and Geotechnics, 81, pp. 195-206, 2017. |
| [Kawai1978] | Kawai, T. “New discrete models and their application to seismic response analysis of structures” in Nuclear Engineering and Design, 48, pp. 207–229, 1978. |
| [Rasmussen2018] | (1, 2) Rasmussen, L. L. and A. P. de Assis. “Elastically-homogeneous lattice modeling of transversely isotropic rocks” in Computers and Geotechnics, 104, pp. 96-108, 2018. |
| [Rasmussen2021] | (1, 2) Rasmussen, L. L. “Hybrid lattice/discrete element method for bonded block modeling of rocks” in Computers and Geotechnics, 130, pp. 103907, 2021. |
| [Rasmusssen2022] | Rasmussen, L. L. “A breakable grain-based model for bi-modular rocks” in International Journal of Rock Mechanics and Mining Sciences, 151, pp. 105028, 2022. |
| Was this helpful? ... | 3DEC © 2019, Itasca | Updated: Feb 25, 2024 |