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

(1)kn=ks=knc=ksc=ktc=EA/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)krnc=EJ/dkrsc=EIs/dkrtc=EIt/d

where J is the polar moment of inertia of the interface surface and Is, It are the two principal moments of inertia of the interface surface in 3D. In 3D, discussion below defines krs=(0,krsc,krtc) 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)A={2R,2DπR2,3DIs=It={23R3,2D14πR4,3DJ={0,2D12πR4,3D

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 (s0,t0), (s1,t1), … (sN,tN) that are rotated onto the the (st) plane with coordinate directions in the (st) axes directions. The inertial values can be computed as

(4)Is=112N(siti+1si+1ti)(ti2+titi+1+ti+12)It=112N(siti+1si+1ti)(si2+sisi+1+si+12)J=Is+It

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)[ϵ11ϵ22ϵ33ϵ12ϵ13ϵ23]=[1EνEνEνE1EνEνEνE1E1+νE1+νE1+νE][σ11σ22σ33σ12σ13σ23]

where ϵ denotes strain, ν denotes Poisson’s ratio, and σ 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)[ϵ11ϵ22ϵ33ϵ12ϵ13ϵ23]=[1E1E1E1E1E1E][σ11σ22σ33σ12σ13σ23]

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)[σ11fσ22fσ33fσ12fσ13fσ23f]=[νσ22νσ33νσ11νσ33νσ11νσ22νσ12νσ13νσ23]

The fictitious stress tensor, computed per contact, is converted to a fictitious force (Ff) 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., νσ12, νσ13, and νσ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:

(8)Fc=Fl+Fd+Ff,Mc=M

where Fl is a linear force, Fd is a dashpot force, Ff is the fictitious force corresponding to the Poisson effect, and M is a linear moment. Take F=Fl+Ff 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)F=Fnˆnc+FsM=Mtˆnc+Mb(2Dmodel:Mt0).

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

(10)Fs=Fssˆsc+Fstˆtc(2Dmodel:Fss0)Mb=Mbsˆsc+Mbtˆtc(2Dmodel:Mbt0).

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 Ff 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 Fn based on the normal-force update mode Ml:

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

    (11)Fn={knAgs,gs<00,otherwise

    where gs is the surface gap (opposite of the overlap).

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

    (12)Fn:=max(Fn+knAΔδn,0)

    where Δδn is the relative normal-displacement increment of this equation 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 Ml=1, and
    3. there are no fictitious stresses when unbonded.
  2. Update Fs:

    The shear force Fs is first updated with:

    (13)Fs:=FsksAΔδδs

    where Δδδs is the relative shear-displacement increment of this equation of the i Contact Resolution section.

  3. Update Mt:

    (14)Mt:={0,2DMtkrncΔθt,3D

    where Δθt is the relative twist-rotation increment of this equation of the i Contact Resolution section.

  4. Update Mb:

    (15)Mb:=MbkrsΔθθb

    where Δθθb is the relative bend-rotation increment of this equation of the i Contact Resolution section and krs is the bending stiffness vector.

  5. Enforce the slip criteria:

    If the critical slip distance, cdc is not specified, nor any entries in the cohesion table, a Coulomb friction criterion is enforced on the shear force such that:

    (16)Fs={Fs,FsμFnμFn(Fs/Fs),otherwise.

    The slip state is updated as:

    (17)s={true,Fs=μFnfalse,otherwise.

    If the slip state is true, then the contact is sliding.

    If cdc is specified then the shear displacement is tracked and the shear force is linearly interpolated from the peak value to μFn(Fs/Fs) when the shear displacement achieves cdc. 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 ucs, at which time it ceases. The increment of shear displacement over a timestep adds a normal force in the form

    (18)Fn=Fn+knAtan(ψπ/180)Δδδs

    where ψ is the dilation angle in degrees.

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

    (19)Mb={Mb,MbMbMb(Mb/Mb),otherwise.Mt={Mt,MtMtMt(Mt/Mt),otherwise.

    where the critical bending and twisting moments are:

    (20)Mb=2.1λbFnR/4Mt=0.65λtμFnR

    where λb and λt are, respectively, the bending and twisting friction multiplier, that default to 1.0. The critical bending and twisting moments with λb=1.0 and λt=1.0 correspond to the values proposed by [Jiang2015a].

    The bending and twisting slip states are updated as:

    (21)sb={true,Mb=Mbfalse,otherwise.st={true,Mt=Mtfalse,otherwise.

    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 σcdc over which the tensile stress evolves linearly from σ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 Fn and Fs using an incremental formulation:

    (22)Fn:=Fn+knAΔδn+FnfFs:=FsksAΔδδs+Fsf

    where Δδn and Δδδs are, respectively, the relative normal-displacement and shear-displacement increments of this equation of the i Contact Resolution section, and Fnf/ Fsf are the normal and shear fictitious forces, respectively, described above.

  2. Update Mb and Mt:

    (23)Mb:=MbkrncΔθθbMt:={0,2DMtkrsΔθt,3D

    where Δθθb and Δθt are, respectively, the relative bend-rotation and twist-rotation increments of this equation of the i Contact Resolution section.

  3. Update the maximum normal (σ, σ>0 is tension) and shear (τ) stresses at the bond periphery:

    (24)σ=FnA+βMbRIτ=FsA+{0,2Dβ|Mt|RJ,3Dwith β[0,1].

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

    The effective normal stress at the bond periphery, σe, is used for subsequent failure computations, using the pore pressure Pp:

    (25)σe=σPp
  1. Update the bond state:

    Take the shear strength τc=cσctanϕ.

    • If the bond has healed (B=6) then it cannot sustain tensile stresses (i.e., σc=0) but can sustain shear stresses. If σe>0 then the bond breaks in tension (B=1). If τ>τ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 (σe>σc), then the bond enters the softening regime if either a critical elongation σcdc 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 (τ>τ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 σ 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 σeσ both the normal force and bending moment are reduced to the softening envelope

      (26)Fn:=Fn(σ/σe)Mb:=Mb(σ/σe)

    The bond remains in the softening regime (B=4) if this is the case. If σe<σ 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, Ek, stored in the linear springs plus the fictitious stress correction;
  • slip energy, Eμ, defined as the total energy dissipated by frictional slip;
  • dashpot energy, Eβ, defined as the total energy dissipated by the dashpots; and
Table 1: Spring Network Model Energy Partitions
Keyword Symbol Description Range Accumulated
spring-network:
energy‑strain Ek strain energy [0.0,+) NO
energy‑slip Eμ total energy dissipated by slip (,0.0] YES
Dashpot Group:
energy‑dashpot Eβ total energy dissipated by dashpots (,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)Ek=12(F2nknA+Fs2ksA+Mb2knI+M2tksJ).
  2. Update the slip energy:

    (28)Eμ:=Eμ+ΔEsμ+ΔEbμ+ΔEtμwhereΔEsμ:=12((Fs)o+Fs)ΔδδμsΔEbμ:=12((Mb)o+Mb)ΔθθμbΔEtμ:=12((Mt)o+Mt)ΔθθμtwithΔδδμs=ΔδδsΔδδks=Δδδs(Fs(Fs)oksA)Δθθμb=ΔθθbΔθθkb=Δθθb(Mb(Mb)oknI)Δθθμt=ΔθθtΔθθkt=Δθθt(Mt(Mt)oksJ)
  3. Update the dashpot energy:

    (29)Eβ:=EβFd(˙δδΔt)

    where ˙δδ 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.

Table 2: Spring-Network Model Properties
Keyword Symbol Description Type Range Default Modifiable Inheritable
springnetwork | Model name
Spring-Network Group:
kn kn Normal stiffness [stress/disp./area] FLT [0.0,+) 0.0 YES YES
ks ks Shear stiffness [stress/disp./area] FLT [0.0,+) 0.0 YES YES
fric μ Friction coefficient [-] FLT [0.0,+) 0.0 YES YES
rgap gr Reference gap [length] FLT R 0.0 YES NO
sn_bmul λb Bending Friction multiplier [-] FLT [0.0,+) 1.0 YES YES
sn_tmul λt Twisting Friction multiplier [-] FLT [0.0,+) 1.0 YES YES
sn_mode Ml Normal-force update mode [-] INT {0,1} 0 YES NO
    {0: update is absolute1: update is incremental          
emod E Effective modulus [force/area] FLT [0.0,+) 0.0 NO N/A
kratio κ Normal-to-shear stiffness ratio [-] FLT [0.0,+) 0.0 NO N/A
    κknks          
sn_rmul λ Radius multiplier [-] FLT (0.0,+) 1.0 YES NO
sn_area A Bond area [length*length] FLT (0.0,+) N/A NO NO
user_area A Constant bond area [length*length] FLT (0.0,+) 0.0 YES NO
sn_radius R Effective radius [length] FLT (0.0,+) N/A NO (set via λ) NO
sn_mcf β Moment-contribution factor [-] FLT [0.0,1] 1.0 YES NO
sn_ten σc Tensile strength [stress] FLT [0.0,+) 0.0 YES NO
sn_tendc σcdc Critical elongation [length] FLT [0.0,+) 0.0 YES NO
sn_tentab   Tensile strength multiplier as a function of elongation VEC2 R2   YES NO
sn_coh c Cohesion [stress] FLT [0.0,+) 0.0 YES NO
sn_cohdc cdc Critical slip distance [length] FLT [0.0,+) 0.0 YES NO
sn_cohtab   Cohesive strength multiplier as a function of shear displacement VEC2 R2   YES NO
sn_fa ϕ Friction angle [degrees] FLT [0.0,90.0) 0.0 YES NO
sn_dil ψ Dilation angle [degrees] FLT [0.0,90.0) 0.0 YES NO
sn_dilzero ucs Critical dilation distance [length] FLT [0.0,+) +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
    {0: unbonded1: unbonded & broke in tension2: unbonded & broke in shear3: bonded4: bonded & softening5: bonded & compressing6: healed, can sustain shear but not tension          
sn_shear τc Bond shear strength [stress] FLT [0.0,+) 0.0 NO (set via c and σc) N/A
sn_sigma σ Total normal stress at bond periphery [stress] FLT [0.0,+) 0.0 NO N/A
sn_esigma σe Effective normal stress at bond periphery [stress] FLT [0.0,+) 0.0 NO N/A
sn_tau τ Shear stress at bond periphery [stress] FLT (,+) 0.0 NO N/A
sn_ndisp δn Normal displacement [length] FLT (,+) 0.0 NO N/A
sn_sdisp δs Shear displacement [length] VEC2 R2 0 NO N/A
sn_slip s Slip state [-] BOOL {false,true} false NO N/A
    {true: slippingfalse: not slipping          
sn_slipb sb Bending slip state [-] BOOL {false,true} false NO N/A
    {true: slippingfalse: not slipping          
sn_slipt st Twisting slip state [-] BOOL {false,true} false NO N/A
    {true: slippingfalse: not slipping          
sn_force F Force (contact plane coord. system) VEC R3 0 YES NO
    (Fn,Fss,Fst)(2D model: Fss0)          
sn_moment M Moment (contact plane coord. system) VEC R3 0 YES NO
    (Mt,Mbs,Mbt)(2D model: MtMbt0)          
sn_non_diag nd Off-diagonal computation mode [-] FLT [0,1] 0 YES NO
sn_pois ν Poisson ratio [-] FLT [0,+) 0.0 YES NO
sn_pois_force Fp Poisson force (contact plane coord. system) VEC R3 0 NO NO
sn_porp Pp Pore pressure [stress] FLT (,+) 0.0 YES NO
Dashpot Group:
dp_nratio βn Normal critical damping ratio [-] FLT [0.0,1.0] 0.0 YES NO
dp_sratio βs Shear critical damping ratio [-] FLT [0.0,1.0] 0.0 YES NO
dp_mode Md Dashpot mode [-] INT {0,1,2,3} 0 YES NO
    {0: full normal & full shear1: no-tension normal & full shear2: full normal & slip-cut shear3: no-tension normal & slip-cut shear          
dp_force Fd Dashpot force (contact plane coord. system) VEC R3 0 NO NO
    (Fdn,Fdss,Fdst)(2D model: Fdss0)          
By convention, κ 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 Fl or M may only be effective during the next force-displacement calculation. When Ml=0, the normal component of the linear force is automatically overridden during the next force-displacement calculation.

Surface Property Inheritance

The linear stiffnesses, kn and ks, and the friction coefficient, μ, 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 kn FLT [0.0,+) N/A Normal stiffness
  kratio κ FLT [0.0,+) N/A Normal-to-shear stiffness ratio
bond Bond the contact if gcG
  gap G VEC2 R2 (,0] Gap interval
compute-stiffness Compute stiffnesses with Poisson correction
  emod E FLT [0.0,+) N/A Young’s modulus
  poisson ν FLT [0.0,+) N/A Poisson ratio
unbond Unbond the contact if gcG
  gap G VEC2 R2 (,0] Gap interval
By convention, setting κ 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 (gb) by specifying gb 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
        {1: shear slip state has changed2: bend slip state has changed3: twist slip state has changed4: shear & bend slip states have changed5: shear & twist slip states have changed6: bend & twist slip states have changed7: all slip states have changed
  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
        {1: failed in tension2: failed in shear
  3 FLT [0.0,+) Failure strength [force] (¯σc or ¯τc, according to the failure mode)
  4 FLT [0.0,+) Bond strain energy ¯Ek 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.