FLAC3D Theory and Background • Constitutive Models

Power-Ubiquitous Model

This model simulates the visco-elasto-plastic mechanical behavior of ubiquitous joint rock with creep occurring in a rock matrix similar to the viscoplastic two-component power model. This model combines two existing FLAC3D models: the viscoplastic two-component power model and the elasto-plastic ubiquitous-joint model.

If compared to the viscoplastic two-component power model, this model accounts for the presence of an orientation of weakness (weak plane), and, therefore, yield may occur in either the solid or along the weak plane, or both, which depends on the stress state, the orientation of the weak plane, and the material properties of the solid and weak planes.

General Failure Correction

In the model formulation, the deviatoric stress rate is visco-elasto-plastic, and is expressed as:

(1)\[ \dot{S}_{ij} = 2G (\dot{e}_{ij} - \dot{e}^c_{ij} - \dot{e}^p_{ij})\]

where \(S_{ij}\) and \(e_{ij}\) are the deviatoric parts of the stress and strain rate tensors, \(\dot{S}_{ij}\) and \(\dot{e}_{ij}\) are creep and plastic strain components, respectively, and \(G\) is tangent shear modulus.

The creep strain rate \(\dot{e}^c_{ij}\) is derived using the von Mises stress \(q=\sqrt{3J_2}\) in accordance with the Norton Power Law (where \(J_2\) is the second invariant of deviatoric stress tensor):

(2)\[ \dot{e}^c_{ij} = \dot{e}_{cr} {\partial q \over \partial S_{ij}}\]

where \(\dot{e}_{cr}\), the creep intensity, has two components by definition (The following section is not linked see Section 1.2.2), and \(\partial q / \partial S_{ij}\) is the direction of creep flow, derived from the definition of \(q\):

(3)\[ {\partial q \over \partial S_{ij}} = {3 \over 2} {S_{ij} \over q}\]

The plastic strain rate, \(\dot{e}^p_{ij}\), is defined using the Mohr-Coulomb flow rule:

(4)\[ \dot{e}^p_{ij} = \dot{e}_p {\partial g \over \partial \sigma_{ij}} - {1 \over 3} \dot{e}^p_{vol} \delta_{ij}\]

where the plastic volumetric strain rate, \(\dot{e}^p_{vol}\), is

(5)\[ \dot{e}^p_{vol} = \left[ {\partial g \over \partial \sigma_{11}} + {\partial g \over \partial \sigma_{22}} + {\partial g \over \partial \sigma_{33}} \right]\]

Plastic flow rate intensity, \(\dot{e}_p\), is derived from the Mohr-Coulomb yield criterion, and \(\partial g / \partial \sigma_{ij}\) is expressed using the Mohr-Coulomb potential function \(g^s\) (see the Mohr-Coulomb model).

The volumetric stress rate is elasto-plastic and has the form:

(6)\[ \dot{\sigma}_o = K (\dot{e}_{vol} - \dot{e}^p_{vol})\]

where \(\dot{\sigma}_o = (\dot{\sigma}_{11} + \dot{\sigma}_{22} + \dot{\sigma}_{33})/3\), \(\dot{e}_{vol} = \dot{e}_{11} + \dot{e}_{22} + \dot{e}_{33}\), and \(K\) is tangent bulk modulus.

The final stress tensor is given by the sum of the deviatoric and isotropic parts:

(7)\[ \sigma_{ij} = S_{ij} + \sigma_o \delta_{ij}\]

In the principal axes formulation, the Mohr-Coulomb yield criterion is checked and relevant plastic are applied plastic what are applied?. Note that in the limiting case when the weak plane (ubiquitous joint) is not active, the model behaves like the viscoplastic two-component power model.

Weak Plane Corrections

Stresses calculated in a previous step are analyzed for failure on the weak plane. The global stress components, \(\sigma_{ij}\) (obtained after the application of the plastic corrections), are resolved into local components (see the section on transformation between global and local coordinate system in the ubiquitous-joint model).

The criterion for failure on the plane is a local form of the Mohr-Coulomb yield criterion with tension cutoff. Corrections on the stresses are made on the local axis and then resolved into the global axis.

Note that in the limiting case when the viscoplastic response of the model is disabled, this behaves like the ubiquitous-joint model.

The model implementation closely follows the procedures for the viscoplastic two-component power model and the elasto-plastic ubiquitous-joint model. Model creep response in rock mass is first calculated. Then the procedure follows the lines presented in the ubiquitous-joint model implementation, with the viscoelastic response replacing the “elastic guess.” Finally, the weak plane logic in the ubiquitous-joint model is employed.

In large-strain mode, the orientation of the weak plane is adjusted to account for body rotations, similar to the ubiquitous-joint model.

power-ubiquitous Model Properties

Use the following keywords with the zone property command to set these properties of the power-ubiquitous model.

power-ubiquitous

bulk f

elastic bulk modulus, K

cohesion f

cohesion, c

constant-1 f

power-law constant, A₁

constant-2 f

power-law constant, A₂

dilation f

dilation angle, ψ

exponent-1 f

power-law exponent, n₁

exponent-2 f

power-law exponent, n₂

friction f

angle of internal friction, ϕ

dip f

dip angle [degrees] of weakness plane

dip-direction f

dip direction [degrees] of weakness plane

joint-cohesion f

joint cohesion, cj

joint-dilation f

joint dilation angle, ψj. The default is 0.0.

joint-friction f

joint friction angle, ϕj

joint-tension f

joint tension limit, σtj. The default is 0.0.

normal v

normal direction of the weakness plane, (\(n_x\), \(n_y\), \(n_z\))

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

poisson f

Poisson’s ratio, v

shear f

elastic shear modulus, G

stress-reference-1 f

reference stress, \(\sigma^{ref}_1\)

stress-reference-2 f

reference stress, \(\sigma^{ref}_2\)

tension f

tension limit, σt. The default is 0.0.

young f

Young’s modulus, E

flag-brittle b

[advanced] If true, the tension limit is set to 0 in the event of tensile failure. The default is false.

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.
• Only one of the three options is required to define the orientation of the weakness plane: dip and dip-direction; a norm vector (nx, ny, nz); or three norm components: nx, ny, and nz.
• 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.

Footnotes

Advanced properties have default values and do not require specification for simpler applications of the model.