Soft-Bond Model

The soft-bond model is referred to in commands and FISH by the name softbond.

Introduction

The soft-bond model can be used to simulate both unbonded and bonded systems.

In an unbonded state, it behaves essentially 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.

In its bonded formulation, the behavior is similar to that of a linear parallel bond model, with the possibility for the bond to fail if the bond strength is exceeded either in shear or in tension. However, contrary to the linear parallel bond model, the bond is not removed upon failure. Instead, it may enter into a softening regime until the bond stress reaches a threshold value at which the bond is removed and considered broken. The slope and tensile breakage strength during softening can be specified by the user (via the softening factor and the softening tensile strength factor, respectively). Another difference with the linear parallel bond model is that only one set of stiffnesses is used for both the unbonded and bonded formulations. This behavior is essentially similar to that proposed by [Ma2018], with the difference that the bond elongation used to update the normal stress in the softening regime accounts for both the normal displacement and bending increments.

Behavior Summary

A soft-bond contact can be envisioned 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 a force and moment that act on the two contacting pieces to develop. The maximum normal and shear stresses at the bond periphery are tracked. The bond can be inactive or active. If the bond is inactive, then frictional-strength parameters limit the shear force and moment (but not the normal force). If the bond is active, then it can break in shear or soften before breaking in tension: (a) if the maximum shear stress exceeds the shear strength, then the bond breaks, (b) if the maximum normal stress exceeds the tensile strength, then the bond may soften until the maximum normal stress reaches some fraction (that may be zero) of the tensile strength, at which point, the bond breaks. The bond becomes inactive after it breaks.

Activity-Deletion Criteria

A contact with the soft bond 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.

Force-Displacement Law

The force-displacement law for the soft-bond model updates the contact force and moment:

(1)Fc=F+Fd,Mc=M

where F is the linear force, Fd is the dashpot force, and M is the moment. These forces and moments are updated as described below.

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

(2)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:

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

The cross-sectional properties of the soft-bond contact are updated as:

(4)R=λ{min(R(1),R(2)),ball-ballR(1),ball-facetA={2Rt,2D (t=1)πR2,3DI={23tR3,2D (t=1)14πR4,3DJ={0,2D12πR4,3D

where A is the cross-sectional area, I is the moment of inertia of the cross section (about the line passing through xc and in the direction of Mb), and J is the polar moment of inertia of the cross section (about the line passing through xc and in the direction of ˆnc). The cross section is {rectangular in 2D; circular in 3D}. The property user_area can be used to specify a different cross-sectional area that remains constant, independent of the piece geometries.

The forces and moment are updated as described below for the unbonded or bonded soft-bond contact model.

Unbonded Behavior

When unbonded, the forces 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:

    (5)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:

    (6)Fn:=min(Fn+knAΔδn,0)

    where Δδn is the relative normal-displacement increment of Equation (12) of the “Contact Resolution” section.

    Note that the unbonded behavior cannot sustain tensile normal force.

  2. Update Fs:

    The shear force Fs is first updated with:

    (7)Fs:=FsksAΔδδs

    where Δδδs is the relative shear-displacement increment of :Equation (12) of the “Contact Resolution” section.

  3. Update Mt:

    (8)Mt:={0,2DMtksJΔθt,3D

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

  4. Update Mb:

    (9)Mb:=MbknIΔθθb

    where Δθθb is the relative bend-rotation increment of Equation (12) of the “Contact Resolution” section.

  5. Enforce the slip criteria.

    A Coulomb friction criterion is enforced on the shear force such that

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

    The shear slip state is updated as

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

    Maximum bending and twisting moment criteria are enforced as

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

    where the critical bending and twisting moments are

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

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

    The bend and twist slip states are updated as

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

    Whenever the shear, bend, or twist slip states change, the slip_change callback event occurs where 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 holding 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

The bonded behavior of the soft-bond model is like that of the linear parallel bond model when the tensile-softening factor ζ=0 or the tensile strength reduction factor γ=1. Softening behavior in tension can occur when ζ0 and γ1.0. Softening applies only to the extension response; the shear response does not soften. Once the maximum normal stress at the bond periphery is exceeded during extension, the normal stress is not automatically set to 0 as in the linear parallel bond model. Instead, the normal stress is reduced with continued extension past the tensile strength. Softening continues until the normal stress during extension is less than the product of γ and the tensile strength. Thus, when γ=0, tensile softening occurs until the normal stress reaches 0 during extension. The rate at which softening occurs is governed by the tensile-softening factor ζ. This behavior is inspired by the model of [Ma2018], although the soft-bond behavior is different in a number of ways.

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

  1. Update Fn and Fs using an incremental formulation:

    (15)Fn:=Fn+knAΔδnFs:=FsksAΔδδs

    where Δδn and Δδδs are respectively the relative normal-displacement and shear-displacement increments of Equation (12) of the “Contact Resolution” section.

  2. Update Mb and Mt:

    (16)Mb:=MbknIΔθθbMt:={0,2DMtksJΔθt,3D

    where Δθθb and Δθt are respectively the relative bend-rotation and twist-rotation increments of Equation (12) of the “Contact Resolution” section.

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

    (17)σ=FnA+βMbRIτ=FsA+{0,2Dβ|Mt|RJ,3Dwith β[0,1].
  1. Update the bond state:

    • If the bond is intact (B=3) and the maximum normal stress at bond periphery exceeds the bond tensile strength (σ>σc), then the bond enters in softening regime (B=4), and the maximal bond elongation is set to:

      (18)l=lc(1.0+ζ)with lc=FnknA+βMbRknI

    where lc is the critical bond elongation (at peak strength) and ζ is the bond softening factor.

    • If the bond is in softening regime (B=4), then the normal stress at bond periphery is checked against the softening envelope as follows.

    The maximum stress is given by:

    (19)σ=σc(ll)ζlc

    where the current bond elongation l is given by:

    (20)l=lc+δl+R|δθδθb|

    where δl is a measure of the bond elongation since softening started:

    (21)δl:=δl+Δδn

    and |δθδθb| is a measure of the accumulated bending since softening started.

    If the current tensile stress is greater than σ, then it is projected back on the softening envelope, otherwise the bond enters in compression:

    (22)if σσ then{Fn:=Fn(σ/σ)Mb:=Mb(σ/σ)else {B=5 (compression starts)σm=σlm:=lc+δl+R|δθδθb|

    where σm and lm are respectively the tensile stress and elongation at which the bond enters in compression.

    • If the bond is in 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 state is set to be back in 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.

  1. Check for bond failure:

    • Tensile failure is first checked if the bond is in softening regime (B=4):

      (23)if B=4 and σσcγ then{B=1 (tensile failure)F=0M=0
    • If the bond has not failed in tension, then shear failure is assessed. The shear strength τc=cσtanϕ, where σ=Fn/A is the average normal stress acting on the parallel bond cross section. If the shear-strength limit is exceeded (τ>τc), then break the bond in shear:

      (24)if τ>τc thenB=2 (shear failure)

    Note that if the bond fails in shear, then F and M are not reset to zero. Instead, the frictional strength properties are enforced using equations (10) and (12) and the slip states updated according to (11) and (14).

    If the bond has broken, then the bond_break callback event is triggered, and the behavior becomes that of an unbonded soft-bond model.

Energy Partitions

The soft-bond model provides three energy partitions:

  • strain energy, Ek, stored in the linear springs;

  • 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: Soft-Bond Model Energy Partitions

Keyword

Symbol

Description

Range

Accumulated

Soft-Bond:

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:

    (25)Ek=12(F2nknA+Fs2ksA+Mb2knI+M2tksJ).
  2. Update the slip energy:

    (26)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:

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

    where ˙δδ is the relative translational velocity of Equation (10) of the “Contact Resolution” section.

Properties

The properties defined by the soft-bond model are listed in the table below for a concise reference; see the “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: Soft-Bond Model Properties

Keyword

Symbol

Description

Type

Range

Default

Modifiable

Inheritable

softbond

Model name

Soft-Bond Group:

kn

kn

Normal stiffness [stress/length]

FLT

[0.0,+)

0.0

YES

YES

ks

ks

Shear stiffness [stress/length]

FLT

[0.0,+)

0.0

YES

YES

fric

μ

Friction coefficient [-]

FLT

[0.0,+)

0.0

YES

YES

sb_bmul

λb

Bending-friction multiplier [-]

FLT

[0.0,+)

1.0

YES

YES

sb_tmul

λt

Twisting-friction multiplier [-]

FLT

[0.0,+)

1.0

YES

YES

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

rgap

gr

Reference gap [length]

FLT

R

0.0

YES

NO

sb_mode

Ml

Normal-force update mode [-]

INT

{0,1}

0

YES

NO

{0: update is absolute1: update is incremental

sb_rmul

λ

Radius multiplier [-]

FLT

(0.0,+)

1.0

YES

NO

sb_radius

R

Radius [length]

FLT

(0.0,+)

N/A

NO (set via λ)

NO

sb_area

A

Cross-sectional area

FLT

(0.0,+)

N/A

NO

NO

user_area

A

Cross-sectional area (held constant)

FLT

(0.0,+)

0.0

YES

NO

sb_ten

σc

Tensile strength [stress]

FLT

[0.0,+)

0.0

YES

NO

sb_soft

ζ

Tensile-softening factor [-]

FLT

[0.0,+)

0.0

YES

NO

sb_cut

γ

Tensile strength reduction factor [-]

FLT

[0.0,1]

1.0

YES

NO

sb_coh

c

Cohesion [stress]

FLT

[0.0,+)

0.0

YES

NO

sb_fa

ϕ

Friction angle [degrees]

FLT

[0.0,90.0)

0.0

YES

NO

sb_state

B

Bond state [-]

INT

{0,1,2,3,4,5}

0

NO

NO

{0: unbonded1: unbonded & broke in tension2: unbonded & broke in shear3: bonded4: bonded & softening5: bonded & compressing

sb_shear

τc

Shear strength [stress]

FLT

[0.0,+)

0.0

NO (set via c and ϕ)

N/A

sb_mcf

β

Moment-contribution factor [-]

FLT

[0.0,1]

1.0

YES

NO

sb_sigma

σ

Maximum normal stress at bond periphery

FLT

[0.0,+)

0.0

NO

N/A

sb_tau

τ

Maximum shear stress at bond periphery

FLT

[0.0,+)

0.0

NO

N/A

sb_slip

s

Shear slip state [-]

BOOL

{false,true}

false

NO

N/A

{true: slippingfalse: not slipping

sb_slipb

sb

Bend slip state [-]

BOOL

{false,true}

false

NO

N/A

{true: slippingfalse: not slipping

sb_slipt

st

Twist slip state [-]

BOOL

{false,true}

false

NO

N/A

{true: slippingfalse: not slipping

sb_force

F

Force (contact plane coord. system)

VEC

R3

0

YES

NO

(Fn,Fss,Fst)(2D model: Fss0)

sb_moment

M

Moment (contact plane coord. system)

VEC

R3

0

YES

NO

(Mt,Mbs,Mbt)(2D model: MtMbt0)

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: Soft-Bond Model Methods

Method

Arguments

Symbol

Type

Range

Default

Description

Soft-Bond Group:

area

Set user_area to sb_area

deformability

Set deformability

emod

E

FLT

[0.0,+)

N/A

Effective modulus

kratio

κ

FLT

[0.0,+)

N/A

Normal-to-shear stiffness ratio

bond

Bond the contact if gcG

gap

G

VEC2

R2

(,0]

Gap interval

soft

ζ

FLT

[0.0,+)

0.0

Softening parameter

cut

γ

FLT

[0.0,1.0]

1.0

Softening threshold parameter

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 sb_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 deformabilty method.

Deformability

The deformability can be specified with the deformability method, which sets

(28)kn:=EL,ks:=knκwith L={R(1)+R(2),ball-ballR(1),ball-facet

The first term in this expression is obtained by equating the normal stiffness to the axial stiffness of the volume of material shown in this figure of the linear model formulation.

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 sb_mode is automatically set to 1 (incremental).

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: Soft-Bond Model Callback Events

Event

Array Slot

Value Type

Range

Description

contact_activated

Contact has become active

1

C_PNT

N/A

Contact pointer

Soft-Bond Group:

slip_change

Slip state has changed

1

C_PNT

N/A

Contact pointer

2

INT

{1,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 [stress] (σ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 example “Genesis and Testing of a Soft-Bonded Material” demonstrates usage of the soft-bond contact model.

Model Summary

An alphabetical list of the soft-bond contact model methods is given here. An alphabetical list of the soft-bond contact model properties is given here.

References

[Jiang2015a] (1,2)

Jiang, M., Z. Shen and J. Wang. “A Novel Three-Dimensional Contact Model for Granulates Incorporating Rolling and Twisting Resistances,” Computer and Geotechnics, 65, 147-163, 2015.

[Ma2018] (1,2)

Ma, Y., and H. Huang. “DEM Analysis of Failure Mechanisms in the Intact Brazilian Test,” International Journal of Rock Mechanics and Mining Sciences, 102, 109-119, 2018.