FLAC3D Theory and Background • Constitutive Models

Ubiquitous-Anisotropic Model

This model combines the logic of the anisotropic elastic model and Mohr-Coulomb model with that of the ubiquitous-joint model (note that “weak plane” and “joint” are used interchangeably in the text). The new model has a single orientation of weakness, which matches the orientation of the plane of isotropy. The criterion for failure on the plane, whose orientation is given, consists of a local Mohr-Coulomb yield criterion with tension cutoff. The model can be used to simulate the behavior of layered material, such as clay shales, and account for slip conditions in the direction of layering (weak plane direction). The criterion for failure on the matrix is the conventional Mohr-Coulomb model.

In the FLAC3D/3DEC implementation, new stresses are evaluated for each step using the elastic-anisotropic incremental laws. The new stresses are updated for yielding on the matrix accordingly using Mohr-Coulomb yield and flow rules and on the weak plane and updated accordingly using the joint flow rule. The local shear flow rule is non-associated, and the local tension flow rule is associated.

Consider a local system of reference axes, with \(x'\) and \(y'\) in the “weak” plane, and the \(z'-\) axis normal to it. This axis is a principal direction of elasticity. Also, any two perpendicular directions \(x'\) and \(y'\), which are principal directions of elasticity, can be selected in the isotropic plane. The incremental elastic strain-stress relations used in the local axes are as follows:

(1)\[\begin{split}\begin{matrix} \Delta \sigma_{1'1'} = a_{11} \Delta \varepsilon^e_{1'1'} + a_{12} \Delta \varepsilon^e_{2'2'} + a_{13} \Delta \varepsilon^e_{3'3'} \\ \\ \Delta \sigma_{2'2'} = a_{12} \Delta \varepsilon^e_{1'1'} + a_{11} \Delta \varepsilon^e_{2'2'} + a_{13} \Delta \varepsilon^e_{3'3'} \\ \\ \Delta \sigma_{3'3'} = a_{13} \Delta \varepsilon^e_{1'3'} + a_{12} \Delta \varepsilon^e_{2'2'} + a_{33} \Delta \varepsilon^e_{3'3'} \\ \\ \Delta \sigma_{1'2'} = 2G \Delta \varepsilon^e_{1'2'} \\ \\ \Delta \sigma_{1'3'} = 2G' \Delta \varepsilon^e_{1'3'} \\ \\ \Delta \sigma_{2'3'} = 2G' \Delta \varepsilon^e_{2'3'} \end{matrix}\end{split}\]

where

(2)\[\begin{split}\begin{matrix} a_{11} = {{E'-{\nu'}^2 E} \over {(1+\nu)[(1-\nu)E' - 2{\nu'}^2 E]}} E \\ \\ a_{12} = {{\nu E'+{\nu'}^2 E} \over {(1+\nu)[(1-\nu)E' - 2{\nu'}^2 E]}} E \\ \\ a_{13} = {{\nu E'} \over {(1-\nu)E' - 2{\nu'}^2 E}} E \\ \\ a_{33} = {{(1-\nu) E'} \over {(1-\nu)E' - 2{\nu'}^2 E}} E \end{matrix}\end{split}\]

and

\(E = E_{1'}= E_{2'}\)	Young’s moduli in the plane of isotropy
\(E'= E_{3'}\)	Young’s moduli in the direction normal to the plane of isotropy
\(\nu = \nu_{1'2'}\)	Poisson’s ratio characterizing lateral contraction in the plane of isotropy when tension is applied in this plane
\(\nu' = \nu_{1'3'} = \nu_{2'3'}\)	Poisson’s ratio characterizing lateral contraction in the plane of isotropy when tension is applied in the direction normal to it
\(G = G_{1'2'}\)	shear modulus for the plane of isotropy, and \(G_{1'2'} = E_{1'}/[2(1+\nu_{1'2'})]\).
\(G' = G_{1'3'} = G_{2'3'}\)	shear modulus for any plane normal to the plane of isotropy.

The conversion of stress and strain tensors from local to global axes and vice-versa is obtained by the application of the usual matrix rotation operations.

Let \(\tau\) be defined as the magnitude of the tangential traction component on the weak plane:

(3)\[\tau = \sqrt{\sigma^2_{1'3'} + \sigma^2_{2'3'}}\]

Also, \(\gamma\) is the strain variable associated with \(\tau\):

(4)\[\gamma = \sqrt{\varepsilon^2_{1'3'} + \varepsilon^2_{2'3'}}\]

The local incremental stress-strain equations (Equation (1)) expressed using the above generalized stress and strain components are:

(5)\[\begin{split}\begin{matrix} \Delta \sigma_{1'1'} = a_{11} \Delta \varepsilon^e_{1'1'} + a_{12} \Delta \varepsilon^e_{2'2'} + a_{13} \Delta \varepsilon^e_{3'3'} \\ \\ \Delta \sigma_{2'2'} = a_{12} \Delta \varepsilon^e_{1'1'} + a_{11} \Delta \varepsilon^e_{2'2'} + a_{13} \Delta \varepsilon^e_{3'3'} \\ \\ \Delta \sigma_{3'3'} = a_{13} \Delta \varepsilon^e_{1'3'} + a_{12} \Delta \varepsilon^e_{2'2'} + a_{33} \Delta \varepsilon^e_{3'3'} \\ \\ \Delta \tau = 2G \Delta \gamma \end{matrix}\end{split}\]

Yield and Potential Functions, Plastic Corrections, and Large-Strain Update to Orientation

The yield and potential functions, plastic corrections, and large-strain update to orientation are similar to those joint-related sections described in the ubiquitous-joint model.

Implementation Procedure

The coefficients of the global elasticity matrix are computed and stored in the initialization phase. New stresses are estimated using the elasticity matrix and the total strain increments for the step. The global stress tensor then is resolved in the local axes of the weak plane, and the local yield conditions are tested. If yielding is detected, a relevant local stress correction (derived using the flow rule) is calculated as described above. The stress correction then is resolved into global axes and added to the elastic stress estimate. Finally, in large strain, adjustment of the weak plane orientation is performed to account for rigid body rotations.

Examples

Borehole Behavior in Elastic Anisotropic Materials

ubiquitous-anisotropic Model Properties

Use the following keywords with the zone property (FLAC3D) or block zone property (3DEC) command to set these properties of the anisotropic-elasticity ubiquitous-joint model.

ubiquitous-anisotropic

cohesion f

cohesion of the matrix, \(c\)

dilation f

dilation angle of the matrix, \(\psi\). The default is 0.0.

dip f (3D ONLY)

dip angle [degrees] of weakness plane

dip-direction f (3D ONLY)

dip direction [degrees] of weakness plane

angle f (2D ONLY)

angle [degrees] of weakness plane, taken counterclockwise from the x-axis

friction f

internal angle of friction of the matrix, \(\phi\)

joint-cohesion f

joint cohesion, \(c_j\)

joint-dilation f

joint dilation angle, \(\psi_j\). The default is 0.0.

joint-friction f

joint friction angle, \(\phi_j\)

joint-tension f

joint tension limit, \(\sigma^t_j\). The default is 0.0.

normal v

normal direction of the weakness plane, (2D or 3D vector)

normal-x f

\(x\)-component of the normal direction to the weakness plane, \(n_x\)

normal-y f

\(y\)-component of the normal direction to the weakness plane, \(n_y\)

normal-z f

\(z\)-component of the normal direction to the weakness plane, \(n_z\) (3D ONLY)

poisson-normal f

Poisson’s ratio characterizing lateral contraction in the plane of isotropy when tension is applied normal to the plane, \({\nu}'\) = \(\nu_{13}\) = \(\nu_{23}\)

poisson-plane f

Poisson’s ratio characterizing lateral contraction in the plane of isotropy when tension is applied in the plane, \(\nu\) = \(\nu_{12}\)

shear-normal f

shear modulus for any plane normal to the plane of isotropy, \(G'\) = \(G_{13}\) = \(G_{23}\)

tension f

tension limit of the matrix, \({\sigma}^t\). The default is 0.0.

young-plane f

Young’s modulus in the plane of isotropy, \(E\) = \(E_1\) = \(E_2\)

young-normal f

Young’s modulus normal to the plane of isotropy, \(E'\) = \(E_3\)

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 \(\nu\). When choosing the latter, Young’s modulus \(E\) must be assigned in advance of Poisson’s ratio \(\nu\).

• Only one of the three options is required to define the orientation of the weakness plane: dip and dip-direction; a norm vector (\(n_x, n_y, n_z\)); or three norm components: \(n_x\), \(n_y\), and \(n_z\).

• The tension cut-off is \({\sigma}^t = min({\sigma}^t, c/\tan \phi)\).

• The joint tension limit used in the model is the minimum of the input \(\sigma^t\) and \({c_j}/{\tan \phi_j}\).