Mohr Model

The mohr model with inactive dashpots provides Mohr-Coulomb style slip with optional slip weakening, including dilation and pore pressure effects with healing. It is referred to in commands and FISH by the name mohr.

Introduction

The mohr model is similar to the Mohr-Coulomb and Softening-Healing models in 3DEC. The mohr model provides a linear stress-strain response to the point of either tensile or shear failure. The dashpot component provides viscous behavior. No contact moments are generated. The model breaks in tension once the tensile strength is exceeded. The shear strength is computed based on the cohesion, normal stress, friction angle and pore pressure. Once broken in shear, a residual shear strength is computed based on the residual cohesion, normal stress, residual friction and pore pressure. The reduction in shear strength may happen instantaneously or over a specified slip distance. The model also supports healing should conditions for slip cease.

Behavior Summary

The mohr model provides the behavior of an infinitesimal interface that does not resist relative rotation, meaning that the contact moment equals zero \((\mathbf{M_{c}} \equiv 0)\). The contact force is resolved into mohr and dashpot forces \((\mathbf{F_{c}} = \mathbf{F} + \mathbf{F^{d}})\). Prior to failure, the mohr force is produced by linear springs with normal and shear stiffnesses, \(k_n\) and \(k_s\), respectively. The surface gap, \(g_s\), is defined as the difference between the contact gap, \(g_c\), and the reference gap, \(g_r\). If the reference gap is zero, the notional surfaces coincide with the piece surfaces. Unbonded or broken contacts are active if and only if the surface gap is less than or equal to zero; the force-displacement law is skipped for inactive contacts.

The dashpot force is produced by dashpots with viscosity given in terms of the normal and shear critical-damping ratios, \(\beta_n\) and \(\beta_s\). The dashpot force is affected by the dashpot mode, \(M_{d}\), which provides four combinations of normal and shear behavior. Denote the normal and shear components of the dashpot force by \(F_{n}^{d}\) and \(F_{s}^{d}\), respectively. The normal-behavior mode can be either full or no-tension, where full means that the entire dashpot load is assigned and no-tension means that \(F_{n}^{d}\) is capped to insure that the total normal force (\(F_{n} + F_{n}^{d}\)) does not become tensile. The shear-behavior mode can be either full or slip-cut, where slip-cut means that \(F_{s}^{d}\) is set to zero if the contact is sliding.

The mohr model supports both tensile and shear strengths, with residual tensile and cohesive strengths.

Activity-Deletion Criteria

Unlike many of the other contact models, there is no bond method; the notion of bonding does not apply to this model. Thus a contact the mohr model is active if the surface gap is less than or equal to zero or if it has not failed in shear or tension. If failure has occurred and the surface gap is greater than zero, the contact will become inactive. The force-displacement law is skipped for inactive contacts.

Force-Displacement Law

The discussion below uses notions of normal/shear displacement increments and directions as outlined in the “Contact Resolution” section.

No contact moments are produced in the mohr model, meaning that the contact moment is zero

(1)\[\mathbf{M_{c}} \equiv 0\]

The force-displacement law for the mohr model updates the contact force, \(\mathbf{F_{c}}\), depending on whether or not failure has previously occurred and whether it occurs during this timestep.

The total contact force is the combination of the mohr and dashpot forces:

(2)\[\mathbf{F_{c}} =\mathbf{F} + \mathbf{F^{d}}\]

where \(\mathbf{F}\) is the mohr force and \(\mathbf{F^{d}}\) is the dashpot force. Both forces are resolved into normal and shear forces:

(3)\[\mathbf{F} =-F_{n} {\kern 1pt} \hat{\mathbf{n}}_\mathbf{c} +\mathbf{F_{s}} ,\quad \mathbf{F^{d}} =-F_{n}^{d} {\kern 1pt} \hat{\mathbf{n}}_\mathbf{c} +\mathbf{F_{s}^{d}}\]

where \(F_{n}>0\), \(F_{n}^{d}>0\) is tension and \(\hat{\mathbf{n}}_\mathbf{c}\) is the contact normal direction. Thus the shear force vector lies in the contact plane, and the shear force is expressed in the contact plane coordinate system:

(4)\[\mathbf{F_{s}} = F_{ss} {\kern 1pt} \hat{\mathbf{s}}_\mathbf{c} + F_{st} {\kern 1pt} \hat{\mathbf{t}}_\mathbf{c} \; \; \left( \mbox{2D model: } F_{ss} \equiv 0 \right).\]

The force-displacement law for the mohr model consists of the following steps:

  1. Update the normal force and tensile failure state.

    A trial normal force is computed via an incremental update

    (5)\[F_{n}^{*} = \left( F_{n} \right)_{o} + A \left( k_n \Delta\delta_n + P_p \right)\]

    where \(\left( F_{n} \right)_{o}\) is the normal force at the beginning of the timestep, \(A\) is the contact area, \(k_n\) is the normal stiffness per unit area, \(\Delta \delta _{n}\) is the relative normal-displacement increment and \(P_p\) is the pore pressure.

    The limiting tensile force is defined as

    (6)\[\begin{split}F_{n}^{l} = \left\{ \begin{array}{rl} \sigma_c A, & B = 0 \\ \sigma_c^{*} A, & \mbox{otherwise} \end{array} \right.\end{split}\]

    where \(\sigma_c\) is the tensile strength per unit area, \(\sigma_c^{*}\) is the residual tensile strength per unit area, \(A\) is the area and \(B\) is the integer failure state with \(B=0\) signifying no failure.

    The tensile force is updated as

    (7)\[\begin{split}F_{n} = \left\{ \begin{array}{rl} F_{n}^{*}, & F_{n}^{*} \le F_{n}^{l} \\ F_{n}^{l}, & \mbox{otherwise} \end{array} \right.\end{split}\]

    The second condition corresponds to tensile failure and results in the integer failure state being non-zero ([1]). The contact is in a state of total tensile failure if the residual tensile strength is zero (i.e., \(\sigma_c^{*} = 0\)). In this case, the shear force is null (\(\mathbf{F_{s}} = \mathbf{0}\)) and steps 2 and 3 below are skipped.

  1. Update the slip magnitude, shear force and shear failure state if total tensile failure has not occurred.

    If sliding is not indicated from the previous timestep and healing is active, the slip magnitude is reset (i.e., \(s = 0\)). If these conditions are not met, the slip magnitude is updated as

    (8)\[s = s_{o} + \lVert \Delta \pmb{δ} _\mathbf{s} \rVert\]

    where \(s_{o}\) is the slip at the beginning of the timestep.

    A trial shear force is computed as

    (9)\[\mathbf{F_{s}^{*}} = \left( \mathbf{F_{s}} \right)_{o} - k_{s} A \Delta \pmb{δ} _\mathbf{s}\]

    where \(\left( \mathbf{F_{s}} \right)_{o}\) is the shear force at the beginning of the timestep, \(k_{s}\) is the shear stiffness per unit area, \(A\) is the area and \(\Delta \pmb{δ} _\mathbf{s}\) is the relative shear-displacement increment.

    The trial shear force magnitude, \(F_{s}^{*} = \lVert \mathbf{F_{s}^{*}} \rVert\), will be used to test for shear failure. Define the peak and residual shear force magnitudes as

    (10)\[F_{s}^{P} = c A + F_{n} \tan \left( \phi \frac{\pi}{180} \right)\]

    and

    (11)\[F_{s}^{R} = c^{*} A + F_{n} \tan \left( \phi^{*} \frac{\pi}{180} \right)\]

    respectively, where \(c\) is the cohesion, \(c^*\) is the residual cohesion, \(\phi\) is the friction angle in degrees, and \(\phi^*\) is the residual friction in degrees.

    The mohr model supports three methods of computing the limiting shear force magnitude, \(F_{s}^{l}\):

    1. If the slip weakening distance \(d_c = 0\) and no multiplicative table entries are supplied then the limiting shear force magnitude is given by

      (12)\[\begin{split}F_{s}^{l} = \left\{ \begin{array}{rl} F_{s}^{P}, & B = 0 \\ F_{s}^{R}, & \mbox{otherwise} \end{array} \right.\end{split}\]

    2. If \(d_c > 0\) and no multiplicative table entries are supplied then

      (13)\[\begin{split}F_{s}^{l} = \left\{ \begin{array}{rl} F_{s}^{P}, & s = 0 \\ F_{s}^{P} - \left( F_{s}^P - F_{s}^{R} \right) \frac{s}{d_c}, & 0 < s \le d_c \\ F_{s}^{R}, & s > d_c \end{array} \right.\end{split}\]

      where \(s\) is the accumulated shear slip magnitude.

    3. Suppose that multiplicative table entries \(m_1, m_2, ... m_n\) have been supplied as a function of slip magnitudes \(s_1, s_2, ... s_n\) in order of increasing slip magnitude. The limiting shear force magnitude is linearly interpolated from the table entries. Suppose that \(s_i < s \le s_{i+1}\). The multiplicative factor is given by

      (14)\[m = 1-\left(\frac{m_{i+1} - m_i}{s_{i+1} - s_i} \left(s - s_i \right) + m_i \right)\]

      where \(m\) goes from 0 to 1 as \(s\) goes from 0 to \(s_n\), leading to

      (15)\[F_{s}^{l} = F_{s}^{P} - \left( F_{s}^P - F_{s}^{R} \right) m\]

      This formulation provides \(F_s^l = F_{s}^{R}\) for \(s \ge s_n\).

    With the limiting shear force magnitude computed the shear force is updated as

    (16)\[\begin{split}\mathbf{F_{s}} = \left\{ \begin{array}{rl} \mathbf{F_{s}^{*}}\frac{F_s^l}{F_s^*}, &F _s^* \ge F_s^l > 0 \\ \mathbf{F_{s}^{*}} , & \mbox{otherwise} \end{array} \right.\end{split}\]

    The first condition above corresponds to slip in which case \(B\) is updated to reflect the sliding state ([1]).

  2. Update the normal force if slip occurs (as indicated in step 2) and the dilation angle, \(\psi\), is nonzero. Dilation occurs until the total shear displacement magnitude reaches \(u_{cs}\), after which time it ceases. One can also define a residual dilation angle, \(\psi^*\). Suppose that the total shear displacement magnitude is less than \(u_{cs}\) and slip has been indicated for this timestep. Define the peak and residual dilation forces as

    (17)\[F_d^P = k_n {\kern 1pt} A{\kern 1pt} \tan (\psi \frac{\pi}{180} ) \lVert \Delta \pmb{δ}_\mathbf{s} \rVert\]

    and

    (18)\[F_d^R = k_n {\kern 1pt} A{\kern 1pt} \tan (\psi^* \frac{\pi}{180} ) \lVert \Delta \pmb{δ}_\mathbf{s} \rVert\]

    respectively.

    If no failure has occurred previous to this timestep then the normal force is updated as

    (19)\[F_{n} = F_{n} + F_d^P\]

    Suppose that failure has occurred previously. If \(d_c = 0\) and no multiplicative table entries are supplied then

    (20)\[F_{n} = F_{n} + F_d^R\]

    If \(d_c \gt 0\) and no multiplicative table entries are supplied then

    (21)\[\begin{split}F_{n} = F_n + \left\{ \begin{array}{rl} F_d^P, & s = 0 \\ F_d^P - \left(F_d^P- F_d^R \right) \frac{s}{d_c}, & 0 < s \le d_c \\ F_d^R, & s > d_c \end{array} \right.\end{split}\]

    If multiplicative table entries are supplied then

    (22)\[F_{n} = F_n + F_d^P - \left(F_d^P- F_d^R \right) m\]

    where \(m\) is the interpolated multiplicative factor defined above.

  3. Update the normal force to remove the pore pressure term.

    The pore pressure is removed from the normal force via

    (23)\[F_{n} = F_{n} - A P_p\]

  4. Update the dashpot normal force based on the dashpot behavior mode:

    (24)\[\begin{split}\begin{array}{l} F_n^d = \left\{ \begin{array}{rl} F^*, & \mbox{$M_d = \{ 0, 2 \}$ (full normal)} \\ \min\left( F^*, -F_n \right), & \mbox{$M_d = \{ 1, 3 \}$ (no-tension normal)} \end{array} \right. \\[3mm] \mbox{with }\quad F^* = \left( 2 \beta_n \sqrt{m_c k_n} \right) \dot{\delta}_n \\ \hskip{15mm} m_c = \left\{ \begin{array}{rl} \displaystyle\frac{m^{(1)} m^{(2)}} {m^{(1)} + m^{(2)}}, & \mbox{ball-ball} \\ m^{(1)}, & \mbox{ball-facet} \end{array} \right. \end{array}\end{split}\]

    where \(m^{(b)}\) is the mass of body \((b)\) and \(\dot{\delta}_n\) is the relative normal translational velocity of Equation (10) of the “Contact Resolution” section. In this expression, \(F^*\) is the entire dashpot load. For the dashpot normal-behavior mode of no-tension, \(F_n^d\) is capped to ensure that the total normal force \((F_n + F_n^d)\) does not become tensile.

  5. Update the dashpot shear force based on the dashpot behavior mode:

    (25)\[\begin{split}\begin{array}{l} \mathbf{F_s^d} = \left\{ \begin{array}{rl} \left( 2 \beta_s \sqrt{m_c k_s} \right) \dot{\pmb{δ}}_\mathbf{s}, & \mbox{$s =$ false or $M_d = \{ 0, 1 \}$ (full shear)} \\ \mathbf{0}, & \mbox{$s =$ true and $M_d = \{ 2, 3 \}$ (slip-cut)} \end{array} \right. \\[3mm] \mbox{with }\quad m_c = \left\{ \begin{array}{rl} \displaystyle\frac{m^{(1)} m^{(2)}} {m^{(1)} + m^{(2)}}, & \mbox{ball-ball} \\ m^{(1)}, & \mbox{ball-facet} \end{array} \right. \end{array}\end{split}\]

    where \(m^{(b)}\) is the mass of body \((b)\) and \(\dot{\pmb{δ}}_\mathbf{s}\) is the relative shear translational velocity of Equation (10) of the “Contact Resolution” section. For the dashpot shear-behavior mode of slip-cut, \(\mathbf{F_s^d}\) is set to zero if the linear spring is sliding \(\mbox{($s =$ true)}\).

Energy Partitions

The mohr model provides two energy partitions:

  • strain energy, \(E_{k}\) and

  • 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.

Table 1: Mohr Model Energy Partitions

Keyword

Symbol

Description

Range

Accumulated

energy‑strain

\(E_{k}\)

strain energy

\([0.0,+\infty)\)

NO

energy‑slip

\(E_{\mu}\)

total energy dissipated by slip

\((-\infty,0.0]\)

YES

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.

  1. Update the strain energy:

    (26)\[{\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} \right).\]

  2. Update the slip energy:

    (27)\[E_{\mu } := - {\tfrac{1}{2}} \left(\left(\mathbf{F_{s}} \right)_{o} +\mathbf{F_{s}} \right)\cdot \Delta \pmb{δ} _\mathbf{s}^{\mu }\]

    with

    (28)\[\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)\]

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

Properties

The properties defined by the mohr 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: Mohr Model Properties

Keyword

Symbol

Description

Type

Range

Default

Modifiable

Inheritable

mohr | Model name

area

\(A\)

Contact area [length*length]

FLT

\([0.0,+\infty)\)

0.0

YES

NO

cohesion

\(c\)

Cohesion [stress]

FLT

\([0.0,+\infty)\)

0.0

YES

NO

cohesion-residual

\(c^*\)

Residual cohesion [stress]

FLT

\([0.0,+\infty)\)

0.0

YES

NO

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

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

dashpot-ratio-normal

\(\beta_n\)

Normal critical damping ratio [-]

FLT

\([0.0,1.0]\)

0.0

YES

NO

dashpot-ratio-shear

\(\beta_s\)

Shear critical damping ratio [-]

FLT

\([0.0,1.0]\)

0.0

YES

NO

dilation

\(\psi\)

Dilation [degrees]

FLT

\([0.0,+\infty)\)

0.0

YES

NO

dilation-residual

\(\psi^*\)

Residual dilation [degrees]

FLT

\([0.0,+\infty)\)

0.0

YES

NO

dilation-zero

\(u_{cs}\)

Dilation cap [length]

FLT

\([0.0,+\infty)\)

\(+\infty\)

YES

NO

disp-normal

Normal displacement [length]

FLT

\((-\infty,+\infty)\)

0

YES

NO

disp-shear

Shear displacement (contact plane coord. system)

VEC

\(\mathbb{R}^3\)

\(\mathbf{0}\)

YES

NO

force-normal

\(F_n\)

Normal force [force]

FLT

\((-\infty,+\infty)\)

0.0

YES

NO

force-shear

\(\mathbf{F_{s}}\)

Shear force (contact plane coord. system)

VEC

\(\mathbb{R}^3\)

\(\mathbf{0}\)

YES

NO

friction

\(\phi\)

Friction angle [degrees]

FLT

\([0.0,90.0)\)

0.0

YES

NO

friction-residual

\(\phi^*\)

Residual friction [degrees]

FLT

\([0.0,90.0)\)

0.0

YES

NO

healing

Healing [-]

INT

{0,1}

0

YES

NO

pore-pressure

\(P_p\)

Pore pressure [stress]

FLT

\([0.0,+\infty)\)

0.0

YES

NO

rgap

\(g_r\)

Reference gap [length]

FLT

\([0.0,+\infty)\)

0.0

YES

NO

slip

\(s\)

Slip [length]

FLT

\([0.0,+\infty)\)

0.0

NO

NO

slip-dc

\(d_c\)

Critical slip weakening distance [length]

FLT

\([0.0,+\infty)\)

0.0

YES

NO

slip-table

Slip weakening multiplier [-,length]

VEC2

\(\mathbb{R}^2\)

YES

NO

state

\(B\)

Failure state bit mask [-]

INT

[1]

0

YES

NO

state-string

Failure state string [-]

STR

[1]

NO

NO

stiffness-normal

\(k_n\)

Normal stiffness [stress/length]

FLT

\([0.0,+\infty)\)

0.0

YES

NO

stiffness-shear

\(k_s\)

Shear stiffness [stress/length]

FLT

\([0.0,+\infty)\)

0.0

YES

NO

stress-effective

\(\sigma_e\)

Effective normal stress [stress]

FLT

\((-\infty,+\infty)\)

0.0

NO

NO

table-pos

Multiplier table position [-]

INT

1

YES

NO

tension

\(\sigma_c\)

Tensile strength [stress]

FLT

\([0.0,+\infty)\)

0.0

YES

NO

tension-residual

\(\sigma_c^*\)

Residual tensile strength [stress]

FLT

\([0.0,+\infty)\)

0.0

YES

NO

user_area

\(A\)

Constant bond area [length*length]

FLT

\([0.0,+\infty)\)

0.0

YES

NO

[1] The failure state can be plotted and queried with FISH. When querying with FISH, the state is represented by a variable with 64 bits. Each bit refers to a potential state. The potential states, and corresponding state strings, are

State (\(B\))

Value

State String

No failure

0

Failure in shear now

1

slip-n

Failure in tension now

2

tension-n

Failure in shear in the past

4

slip-p

Failure in tension in the past

8

tension-p

Combinations of states, and state strings, are possible by turning on different bits. For example, if tensile failure has occurred in the past and shearing is happening now \(B = 9 = 8 + 1\) and the state string is “slip-n tension-p”.

Methods

Table 4: Mohr Model Methods

Method

Arguments

Symbol

Type

Range

Default

Description

area

Set user_area to area

Area

Set the user_area property via the current area. This operation means that the contact area stays constant and is fixed independent of changes to the piece sizes/geometries. This method is necessary when using rigid blocks, but not when using zone joints.

Callback Events

Table 5: Mohr Model Callback Events

Event

Array Slot

Value Type

Range

Description

slip_change

Slip state has changed

1

C_PNT

N/A

Contact pointer

bond_break

Bond broken when strength (\(c\) and/or \(\sigma_c\)) is non-zero

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

Model Summary

An alphabetical list of the mohr model methods is given here. An alphabetical list of the mohr model properties is given here.