Spring Network Model

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

Introduction

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

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

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

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

Behavior Summary

A spring network bond combines two concepts. The first envisions a bond as a set of elastic springs with constant normal and shear stiffnesses, uniformly distributed over a {rectangular in 2D; circular in 3D} cross-section lying on the contact plane and centered at the contact point. Relative motion at the contact causes development of a linear force and moment that act on the two contacting pieces. This behavior is the same as the linear parallel bond and the soft-bond model. But unlike those models, a method exists to compute contact stiffnesses that are consistent—when bonded—with an elastic material with a target Young’s modulus and zero Poisson’s ratio. If a non-zero Poisson’s ratio is specified then a fictitious force, derived from linear elasticity, is computed and added to the total force to produce the desired Poisson effect. If the maximum normal stress exceeds the bond tensile strength, the bond may enter a user-defined force-displacement regime before breaking. If the maximum shear stress exceeds the bond shear strength, then it breaks and may enter a user-defined force-displacement regime in the shear direction. If the bond is inactive but the contact is active, frictional strength parameters govern slip with potential dilation and/or healing.

Activity-Deletion Criteria

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

Contact Stiffness Computation

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

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

(1)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. [Rasmussen2018b] demonstrates how to extend this formulation to transversely isotropic elastic materials, though this is not part of the current implementation.

Tensile and Cohesion Tables

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

Force-Displacement Law

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

(8)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 Equation (12) of the i Contact Resolution section.

    Note:

    1. the unbonded behavior cannot sustain tensile normal force,

    2. if either the compute-stiffness or assign-stiffness methods are used then 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 Equation (12) of the i Contact Resolution section.

  3. Update Mt:

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

    where Δθt is the relative twist-rotation increment of Equation (12) of the i Contact Resolution section.

  4. Update Mb:

    (15)Mb:=MbkrsΔθθb

    where Δθθb is the relative bend-rotation increment of Equation (12) 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 Equation (12) 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 Equation (12) 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; and

  • dashpot energy, Eβ, defined as the total energy dissipated by the dashpots

Table 1: Spring Network Model Energy Partitions

Keyword

Symbol

Description

Range

Accumulated

spring-network:

energy‑strain

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 Equation (10) of the i Contact Resolution section.

Properties

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

Table 2: Spring-Network Model Properties

Keyword

Symbol

Description

Type

Range

Default

Modifiable

Inheritable

springnetwork | Model name

Spring-Network Group:

kn

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. E. Houseworth, J. T. Birkholzer and J. E. Bolander. “Simulating the Poisson Effect in Lattice Models of Elastic Continua,” Computers and Geotechnics, 70, 60–67, 2015.

[Asahina2017] (1,2)

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

[Bolander1998]

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

[Kawai1978]

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

[Rasmussen2018a]

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

[Rasmussen2018b] (1,2,3)

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

[Rasmussen2019]

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

[Rasmussen2021b]

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

[Rasmussen2022] (1,2)

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