FLAC3D Theory and Background • Constitutive Models

Burgers-Mohr Model

A viscoplastic model in FLAC3D or 3DEC is characterized by a visco-elasto-plastic deviatoric behavior and an elasto-plastic volumetric behavior. The viscoelastic and viscoplastic strain-rate components are assumed to act in series. The viscoelastic constitutive law corresponds to a Burgers model (Kelvin cell in series with a Maxwell component), and the plastic constitutive law corresponds to a Mohr-Coulomb model.

As a notation convention in this section, we use the symbols \(S_{ij}\) and \(e_{ij}\) to denote deviatoric stress and strain components:

(1)\[S_{ij}=\sigma_{ij}-{\sigma}_{0}\delta_{ij}\]

(2)\[e_{ij}=\epsilon_{ij}-{e_{vol} \over 3} \delta_{ij}\]

where

(3)\[{\sigma}_{0}={\sigma_{kk} \over 3}\]

and

(4)\[e_{vol}=\epsilon_{kk}\]

Also, Kelvin, Maxwell, and plastic contributions to stresses and strains are labeled using the superscripts .^K, .^M and .^p, respectively. With those conventions, the model deviatoric behavior may be described by the following relations.

Strain rate partitioning:

(5)\[\dot{e}_{ij}=\dot{e}_{ij}^K+\dot{e}_{ij}^M+\dot{e}_{ij}^p\]

Kelvin:

(6)\[S_{ij}=2\eta ^K\dot{e}_{ij}^K+2G^Ke_{ij}^K\]

Maxwell:

(7)\[\dot{e}_{ij}^M={{\dot{S}_{ij}} \over {2G^M}}+{{S_{ij}} \over {2\eta ^M}}\]

Mohr-Coulomb:

(8)\[\dot{e}_{ij}^p = \lambda ^{*}{{\partial g} \over {\partial \sigma_{ij}}}- {1 \over 3} \dot{e}_{vol}^p\delta_{ij}\]

(9)\[\dot{e}_{vol}^p = \lambda ^{*}\left[ {{\partial g} \over {\partial \sigma_{11}} }+{{\partial g} \over {\partial \sigma_{22}}}+{{\partial g} \over {\partial \sigma_{33}}}\right]\]

In turn, the volumetric behavior is given by

(10)\[\dot{\sigma}_0=K(\dot{e}_{vol}-\dot{e}_{vol}^p)\]

In those formulas, the properties \(K\) and \(G\) are the bulk and shear moduli, and \(\eta\) is the viscosity. The Mohr-Coulomb yield envelope is a composite of shear and tensile criteria. The yield criterion is \(f\) = 0, and in the principal axes formulation we have:

Shear yielding:

(11)\[f=\sigma_1-\sigma_3N_{\phi} +2C\sqrt{N_{\phi} }\]

Tension yielding:

(12)\[f=\sigma ^t-\sigma_3\]

where \(C\) is the material cohesion, \(\phi\) the friction, \(N_{\phi}=(1+\sin \phi )/(1-\sin \phi )\), \(\sigma ^t\) is the tensile strength, and \(\sigma_1\) and \(\sigma_3\) are the minimum and maximum principal stresses (compression negative). The potential function \(g\) has the following form:

Shear failure:

(13)\[g=\sigma_1-\sigma_3N_{\psi}\]

Tension failure:

(14)\[g=-\sigma_3\]

where \(\psi\) is the material dilation, and \(N_{\psi} =(1+\sin \psi )/(1-\sin \psi )\). Finally, \(\lambda ^{*}\) is a parameter that is nonzero during plastic flow only, which is determined by application of the plastic yield condition \(f\) = 0.

The model implementation closely follows the procedures described in the manual for the Burgers-creep and Mohr-Coulomb models. The principle is to write Equation (5) to (10) in the form of finite increments:

(15)\[\Delta e_{ij}=\Delta e_{ij}^K+\Delta e_{ij}^M+\Delta e_{ij}^p\]

(16)\[\overline{S_{ij}}\Delta t=2\eta ^K\Delta e_{ij}^K+2G^K\overline{e_{ij}} ^K\Delta t\]

(17)\[\Delta e_{ij}^M={{\Delta S_{ij}} \over {2G^M}}+{{\overline{S_{ij}}} \over {2\eta ^M}}\Delta t\]

(18)\[\Delta \sigma_0=K(\Delta e_{vol}-\Delta e_{vol}^p)\]

where the overbar indicates mean value over the timestep \(\Delta t\):

(19)\[\overline{S_{ij}}={{S_{ij}^N+S_{ij}^O} \over 2}\]

(20)\[\overline{e_{ij}}={{e_{ij}^N+e_{ij}^O} \over 2}\]

and the superscripts .^N and .^O denote new and old values.

After substitution of Equations (19) and (20) in Equation (16), and solving for \(e_{ij}^{K,N}\), the Kelvin strain contribution may be expressed in the form

(21)\[e_{ij}^{K,N}= {1 \over A} \left[ Be_{ij}^{K,O}+{{\Delta t} \over {4\eta ^K}}\left( S_{ij}^N+S_{ij}^O\right) \right]\]

where

(22)\[A = 1+{{G^K\Delta t} \over {2\eta ^K}}\]

(23)\[B = 1-{{G^K\Delta t} \over {2\eta ^K}}\]

After substitution of Equations (17) and (21) in Equation (15) and solving for the new deviatoric stress component, we find (using the mean value definitions in Equations (19) and (20))

(24)\[S_{ij}^N= {1 \over a} \left[ \Delta e_{ij}-\Delta e_{ij}^p+bS_{ij}^O-({B \over A} -1)e_{ij}^{K,O}\right]\]

where

(25)\[a = {1 \over {2G^M}}+{{\Delta t} \over 4} \left( {1 \over {\eta ^M}}+ {1 \over {A\eta ^K}}\right)\]

(26)\[b = { 1 \over {2G^M}}-{{\Delta t} \over 4} \left( {1 \over {\eta ^M}}+ { 1 \over {A\eta ^K}}\right)\]

and Equation (21) is used as an evolution law to evaluate \(e_{ij}^{K,O}\) in Equations (24). For completeness, Equation (18) is written in the form

(27)\[\sigma_0^N=\sigma_0^O+K(\Delta e_{vol}-\Delta e_{vol}^p)\]

In the model implementation in FLAC3D and 3DEC, new trial stress components \(\widehat{S_{ij}^N}\) and \(\widehat{\sigma_0^N}\) are computed from Equations (24) and (27), assuming viscoelastic increments. Trial principal stress components are calculated and sorted, and the yield function is computed. As long as \(f \geq\) 0, the trial stresses are taken for new stresses. If \(f <\) 0, plastic flow is taking place and the trial stresses must be corrected by a component due to incremental plastic strain before their value is assigned to the new stresses and the evolution law is updated. Expressing as Equations (24) and (27) in principal axes, we may then write, by definition of trial stresses:

(28)\[S_i^N = \widehat{S_i^N}- { 1 \over a} \Delta e_i^p\]

(29)\[\sigma_0^N = \widehat{\sigma_0^N}-K\Delta e_{vol}^p\]

or, using the definition for deviatoric components:

(30)\[\sigma_1^N = \widehat{\sigma_1^N}-\left[ \alpha_1\Delta \epsilon _1^p+\alpha_2(\Delta \epsilon_2^p+\Delta \epsilon_3^p)\right]\]

(31)\[\sigma_2^N = \widehat{\sigma_2^N}-\left[ \alpha_1\Delta \epsilon_2^p+\alpha_2(\Delta \epsilon_1^p+\Delta \epsilon_3^p)\right]\]

(32)\[\sigma_3^N = \widehat{\sigma_3^N}-\left[ \alpha_1\Delta \epsilon_3^p+\alpha_2(\Delta \epsilon_1^p+\Delta \epsilon_2^p)\right]\]

where

(33)\[\alpha_1 = K+ {2 \over {3a}}\]

(34)\[\alpha_2 = K- { 1 \over {3a}}\]

Except for the definitions of \(\alpha_1\) and \(\alpha_2\), these formulas are similar to those obtained in the Mohr-Coulomb model derivation. The plasticity formulation may proceed along similar lines. In doing so, for shear yielding we obtain:

(35)\[\sigma_1^N = \widehat{\sigma_1^N}-\lambda (\alpha_1-\alpha_2N_{\psi} )\]

(36)\[\sigma_2^N = \widehat{\sigma_2^N}-\lambda \alpha_2(1-N_{\psi} )\]

(37)\[\sigma_3^N = \widehat{\sigma_3^N}-\lambda (\alpha_2-\alpha_1N_{\psi} )\]

with

(38)\[\lambda ={{\widehat{\sigma_1^N}-\widehat{\sigma_3^N}N_{\phi} +2C\sqrt{ N_{\phi} }} \over {\left( \alpha_1-\alpha_2N_{\psi} \right) -\left( \alpha_2-\alpha _1N_{\psi} \right) N_{\phi} }}\]

and for tensile yielding:

(39)\[\sigma_1^N = \widehat{\sigma_1^N}+\lambda \alpha_2\]

(40)\[\sigma_2^N = \widehat{\sigma_2^N}+\lambda \alpha_2\]

(41)\[\sigma_3^N = \widehat{\sigma_3^N}+\lambda \alpha_1\]

with

(42)\[\lambda ={{\sigma ^t-\widehat{\sigma_3^N}} \over {\alpha_1}}\]

Finally, new global stress components are calculated, assuming that the principal directions have not been affected by the occurrence of plastic flow.

The square root of the second invariant and the modulus of the first invariant of incremental plastic-strain tensor are used as incremental contributions to measure the amount of plastic strain associated with shear and tensile failure, respectively (see the corresponding expression in strain-softening model).

By default, both Maxwell and Kelvin viscosity properties, \(\eta ^M\) and \(\eta^K\), are infinite (although stored as zero in FLAC3D’s property arrays). Note that if the default value for \(\eta ^K\) is adopted, then the model assumes that \(G^K\) = 0, even if a different value has been assigned to that property. The default value for \(G^K\) is zero and the default value for \(G^M\) is 10^-20, irrespective of the system of units adopted.

The default value for the timestep is zero, in which case the program treats the material as elasto-plastic, with only the elastic part of the Maxwell cell active.

If the stresses are changed in a FLAC3D model with the zone initialize command, the internal Kelvin strains, \(e^K_{ij}\), will not be compatible with them and movement will occur until the strains adjust. To avoid this incompatibility, the internal strains may be set to reflect the current values of stresses. The internal Kelvin strains, \(e^K_{ij}\), are available for user inspection and modification as zone property variables: strain-kelvin-xx, strain-kelvin-yy, etc. An example FISH function to perform this step is given here. This function should be invoked immediately following initialization of stresses.

burgers-mohr Model Properties

Use the following keywords with the zone property (FLAC3D) or zone property (3DEC) command to set these properties of the Burgers-Mohr model.

burgers-mohr

bulk f

bulk modulus, \(K\)

cohesion f

cohesion, \(c\)

dilation f

dilation angle in degrees, \(ψ\)

friction f

internal friction angle in degrees, \(ϕ\)

shear-kelvin f

Kelvin shear modulus, \(G^K\)

shear-maxwell f

shear modulus, \(G^M\)

strain-kelvin-xx f

Kelvin strain, \(e^K_{xx}\)

strain-kelvin-xy f

Kelvin strain, \(e^K_{xy}\)

strain-kelvin-xz f

Kelvin strain, \(e^K_{xz}\)

strain-kelvin-yy f

Kelvin strain, \(e^K_{yy}\)

strain-kelvin-yz f

Kelvin strain, \(e^K_{yz}\)

strain-kelvin-zz f

Kelvin strain, \(e^K_{zz}\)

tension f

tension limit \(σ^t\). The default is zero.

viscosity-kelvin f

Kelvin viscosity, \(η^K\)

viscosity-maxwell f

Maxwell dynamic viscosity, \(η^M\)

strain-shear-plastic f (r)

accumulated plastic shear strain

strain-tensile-plastic f (r)

accumulated plastic tensile strain

Key

(r) Read-only property.

This property cannot be set by the user. Instead, it can be listed, plotted, or accessed through FISH.

Notes

• Only one of the two options is required to define the elasticity: bulk modulus \(K\) and shear modulus \(G\), or Young’s modulus \(E\) and Poisson’s ratio \(v\).
• The tension cut-off is \(σ^t\) = min (\(σ^t\), \(c\)/tan\(ϕ\)).
• The creep behavior is triggered by deviatoric stress, while the volumetric behavior does not consider creep.
• If the stresses are changed in a FLAC3D model with the zone initialize command, the internal Kelvin strains, \(e^K_{ij}\), will not be compatible with them, and movement will occur until the strains adjust. To avoid this incompatibility, the internal strains may be set to reflect the current values of stresses. The internal Kelvin strains, \(e^K_{ij}\), are available for user inspection and modification. An example FISH function is SetKStrain.f3fis. This function should be invoked immediately following initialization of stresses.