FLAC3D Theory and Background • Constitutive Models

Bilinear Strain-Softening/Hardening Ubiquitous-Joint (SUBI) Model

The bilinear strain-softening/hardening ubiquitous-joint (or, SUBI, for simplicity) model is a generalization of the ubiquitous-joint model. In the bilinear model, the failure envelopes for the matrix and joint are the composite of two Mohr-Coulomb model criteria with a tension cutoff that can harden or soften according to specified laws. A nonassociated flow rule is used for shear-plastic flow and an associated flow rule is used for tensile-plastic flow.

The softening behaviors for the matrix and the joint are specified in tables in terms of four independent hardening parameters (two for the matrix and two for the joint) that measure the amount of plastic shear and tensile strain, respectively. In this numerical model, general failure is first detected for the step, and relevant plastic corrections are applied. The new stresses are then analyzed for failure on the weak plane and updated accordingly. The hardening parameters are incremented if plastic flow has taken place, and the parameters of cohesion, friction, dilation, and tensile strength are adjusted for the matrix and the joint using the tables.

Formulations

Failure Criterion and Flow Rule for the Matrix

The criterion for failure in the matrix used in this model is sketched in the principal stress plane (σ1,σ3) in Figure 1. (Recall that compressive stresses are negative and, by convention, σ1σ2σ3.)

The failure envelope is defined by two Mohr-Coulomb failure criteria: fs2 = 0 and fs1 = 0 for segments AB and BC, and a tension failure criterion ft = 0 for segment CD.

The shear failure criterion has the general form fs = 0. The criterion is characterized by a cohesion, c2, and a friction angle, ϕ2, for segment AB, and by a cohesion, c1, and a friction angle, ϕ1, for segment BC. The tensile failure criterion is specified by means of the tensile strength, σt (positive value); thus we have

(1)fs=σ1σ3 Nϕ+2cNϕ
(2)ft=σ3σt

where

(3)Nϕ=1+sin(ϕ)1sin(ϕ)

Note that the tensile cap acts on segment BC of the shear envelope and, for a material with nonzero friction angle ϕ1, the maximum value of the tensile strength is given by

(4)σtmax=c1tanϕ1

../../../../_images/modelsubi-1.png

Figure 1: FLAC3D bilinear matrix failure criterion.


In the model formulation, elastic guesses for the stresses are first evaluated for the step using total strain increments. Plastic yielding is detected if the corresponding stress point (σI1,σI3) lies outside the failure surface representation in Figure 1. In this case, a stress correction must be applied to the elastic guess. It is determined by allowing plastic flow to occur, in order to restore the condition fs2 = 0, fs1 = 0 or ft = 0, depending on the position of the stress point above AB, BC, or CD. (Bisectors are used at B and C to delimit the domain of failure attached to a particular segment of the yield surface.)

The usual assumption is made that total strain increments can be decomposed into elastic and plastic parts. The flow rule for plastic yielding has the form

(5)Δϵpi=λgσi

where i = 1,3. The potential function, g, for shear yielding is gs. This function corresponds to the nonassociated law,

(6)gs=σ1σ3Nψ

where ψ, the dilation angle, is equal to ψ2 for failure along AB, and ψ1 along BC, and

(7)Nψ=1+sin(ψ)1sin(ψ)

The potential function for tensile yielding is gt. It corresponds to the associated flow rule,

(8)gt=σ3

The plastic strain increments for shear failure have the form

(9)Δϵps1=λsΔϵps2=0Δϵps3=λsNψ

The stress corrections for shear failure are

(10)Δσ1=λs(α1α2Nψ)Δσ2=λsα2(1Nψ)Δσ3=λs(α1Nψ+α2)

where

(11)λs=σI1σI3 Nϕ+2cNϕ(α1α2 Nψ)(α1 Nψ+α2) Nϕ

and, by definition:

(12)α1=K+43 Gα2=K23 G

In turn, the plastic strain increments for tensile failure have the form

(13)Δϵpt1=0Δϵpt2=0Δϵpt3=λt

The stress corrections for tensile failure are

(14)Δσ1=λtα2Δσ2=λtα2Δσ3=λtα1

where

(15)λt=σI3σtα1

Failure Criterion and Flow Rule for the Weak Plane

The stresses, corrected for plastic flow in the matrix, are resolved into components parallel and perpendicular to the weak plane and tested for ubiquitous-joint failure. The failure criterion is expressed in terms of the magnitude of the tangential traction component, τ=σ132+σ232, and the normal traction component, σ33, on the weak plane.

The failure criterion is represented in Figure 1 and corresponds to two Mohr-Coulomb failure criteria (fs2 = 0 for segment AB, and fs1 = 0 for segment BC) and a tension failure criterion (ft = 0 for segment CD). Each shear criterion has the general form fs = 0 and is characterized by a cohesion and a friction angle cj,ϕj, equal to cj2,ϕj2 along segment AB, and cj1,ϕj1 along BC. The tensile criterion is specified by means of the tensile strength, σtj (positive value). Thus, we have

(16)fs=τ+σ33 tanϕjcj
(17)ft=σ33σtj

Note that for a weak plane with nonzero friction angle ϕj1, the maximum value of the tensile strength is given by

(18)σtjmax=cj1tanϕj1

../../../../_images/modelsubi-2.png

Figure 2: FLAC3D bilinear joint failure criterion.


Yield is detected and stress corrections are applied using a technique similar to the one described in the matrix context.

Here, the flow rule for plastic yielding has the form

(19)Δϵps33=λgσ33Δγps=λgτ

where γ is the strain variable associated with τ, and we have

(20)Δγps=Δϵps132+Δϵps232

The potential function, g, for shear yielding on the weak plane is gs. It corresponds to the nonassociated law,

(21)gs=τ+σ33tanψj

where ψj, the dilation angle, is equal to ψj2 for failure along AB, and ψj1 along BC.

The potential function g for tensile yielding on the weak plane is gt. It corresponds to the associated flow rule,

(22)gt=σ33

The local plastic strain increments for shear failure are such that

(23)Δϵps33=λstanψjΔγps=λs

The stress corrections for shear failure are

(24)Δσ11=λsα2tanψjΔσ22=λsα2tanψjΔσ33=λsα1tanψjΔτ=λs2G

where

(25)λs=τO+σO33 tanϕjcj2G+α1tanψjtanϕj

and the superscript O indicates values obtained just before detection of failure on the weak plane.

The plastic corrections for the local shear stress components on the weak plane are derived by scaling:

(26)Δσ13=ΔτσO13τOΔσ23=ΔτσO23τO

In turn, local plastic strain increments for tensile failure have the form

(27)Δϵpt33=λtΔγpt=0

The stress corrections for tensile failure are

(28)Δσ11=λtα2Δσ22=λtα2Δσ33=λtα1

where

(29)λt=σO33σtjα1

Large-Strain Update of Orientation

In large strain, the orientation of the weak plane is adjusted per zone to account for rigid-body rotations, and rotations due to deformations. The large-strain update to orientation is identical to those described in the ubiquitous-joint model.

Hardening Parameters

In the bilinear strain-softening/hardening ubiquitous-joint model, some or all of the zone yielding parameters (i.e., cohesion, friction, dilation and tensile strength) for the matrix and joint are modified automatically, after the onset of plasticity, according to piecewise linear laws specified on input in terms of a range of values for the hardening parameters. (See the section of user-defined functions in the strain-softening/hardening model.) One table number must be specified in the zone property (FLAC3D) or block zone property (3DEC) command for each softening parameter. (If no table property number is specified, the parameter is taken as constant.) The corresponding table data contain pairs of values for the parameter and the property between which a linear variation is assumed. The last property value is used for values of the hardening parameter beyond the last one specified in the table.

Four independent hardening parameters are used in this model:

  1. κs measures the matrix plastic shear strain and is used to update the matrix cohesion, friction, and dilation;
  2. κt measures the matrix plastic volumetric tensile strain and is used to update the matrix tensile strength;
  3. κsj estimates the joint plastic shear strain and controls the joint cohesion, friction, and dilation update; and
  4. κtj evaluates the joint plastic volumetric tensile strain and controls the joint tensile strength update.

The parameters are defined as the sum of incremental measures of plastic strain for the zone. The zone-hardening increments are calculated as the volumetric average of hardening increments over all tetrahedra involved in the zone.

The shear-hardening increment for a tetrahedron is the square root of the second invariant of the incremental plastic shear-strain deviator tensor for the step. For the matrix, it is given as

(30)Δκs=12(Δϵps1Δϵpsm)2+(Δϵpsm)2+(Δϵps3Δϵpsm)2

where Δϵpsm is the mean plastic strain increment,

(31)Δϵpsm=13(Δϵps1+Δϵps3)

and the plastic strain increments are given by Equation (9), using Equation (11) for λs.

For the joint, the formula is

(32)Δκsj=13(Δϵps33)2+(Δϵps13)2+(Δϵps23)2

where the plastic strain increments are given by Equation (23) (see Equation (20)), using the form (25) for λs.

The tetrahedron tensile-hardening increment is the plastic volumetric tensile-strain increment.

For the matrix, we have

(33)Δκt=Δϵpt3

where the plastic strain increment is given by Equation (13), using Equation (15) for λt.

For the joint, the expression is

(34)Δκtj=Δϵpt33

where the plastic strain increment is given by Equation (27), using the Equation (29) for λt.

Implementation Procedure

The implementation of the bilinear model in FLAC3D proceeds as indicated above. An elastic guess, σIij, is first computed using stress increments for the step evaluated by application of Hooke’s law to the total strain increments, Δϵij. Principal stresses are calculated, ordered such that σI1σI2σI3, and tested for failure in the matrix using the yield criteria defined by Equation (1) and (2). In principle, matrix failure is declared if the representation of the stress point (σI1,σI3) falls outside the yield surface in Figure 1. In this case, stress corrections are applied to the principal values of the elastic guess, which depend on the position of the stress point above AB, BC or CD. (Bisectors are used at B and C to delimit the domain of failure attached to a particular segment of the yield surface.)

The stress corrections for shear failure in the matrix are given by Equations (10) and (11), where the parameters of cohesion, c, friction, ϕ, and dilation, ψ, have value c2,ϕ2,ψ2 for failure along AB, and c1,ϕ1,ψ1 for failure along BC.

The stress corrections for tensile failure in the matrix are given by by Equations (14) and (15).

The stress tensor components in the system of reference axes, σOij, are then calculated from the corrected principal values by assuming that the principal directions have not changed during plastic flow.

Local traction components on the weak plane are defined as σ33 and τ, with σ33 being the normal component, and τ=σ132+σ232 being the magnitude of the tangential traction component. (The local axis 1สน is in the dip direction, 2สน the strike, and 3สน the normal to the weak plane.) These stresses are resolved from σOij and examined for ubiquitous-joint failure using the yield criteria (Equations (16) and (17)). In principle, ubiquitous-joint failure is declared if the representation of the stress point (σO33,τO) falls outside the yield surface in Figure 2. In this case, local stress corrections, which depend on the position of the stress point in the vicinity of AB, BC or CD, are applied. (Bisectors are used at B and C to delimit the domain of failure attached to a particular segment of the yield surface.)

The stress corrections for shear joint failure are given by Equations (24) to (26), where the parameters of cohesion, cj, friction, ϕj, and dilation, ψj, have values cj2,ϕj2,ψj2 for failure along AB, and cj1,ϕj1,ψj1 for failure along BC.

The stress corrections for tensile joint failure are given by Equations (28) and (29).

Finally, the local stress components are resolved back into global axes.

In large-strain mode, the unit normal to the weak plane is adjusted per zone to account for body rotations.

After determination of the new stresses for the step, the hardening parameters are incremented using Equations (30), (32), (33), and (34). These parameters are then used to determine new values of cohesion, friction, dilation and tensile strength for the matrix and the joint from the available input tables.

It is assumed that the tensile strength of the material can never increase. Also, for material with friction, the value will not exceed the maximum value σtmax or σtjmax (see Equations (4) and (18)).

Note that, by default, the yield model is linear in both the matrix and the joint, in which case only sections 1 (where fs1 = 0) and 3 (where ft = 0) of the yield curve are recognized (even if properties are assigned for section 2, where fs2 = 0). To activate the bilinear laws, the property flag-bilinear and/or the property flag-bilinear-joint must be set to 1.

Also, if the friction angles for sections 1 and 2 become equal, the model will be considered as linear, and section 2 will be ignored (for the matrix and/or the joint, as appropriate). Section 2 will also be ignored if the intersection of sections 1 and 2 corresponds to a stress point that violates the tensile criterion.



softening-ubiquitous model properties

Use the following keywords with the zone property (FLAC3D) or block zone property (3DEC) command to set these properties of the bilinear strain-softening/hardening ubiquitous-joint model.

softening-ubiquitous
bulk f

elastic bulk modulus, K

cohesion f

matrix cohesion, c = c1

cohesion-2 f

matrix cohesion, c2

dilation f

matrix dilation angle, ψ = ψ1. The default is 0.0.

dilation-2 f

matrix dilation angle, ψ2. The default is 0.0.

dip f

dip angle [degrees] of weakness plane

dip-direction f

dip direction [degrees] of weakness plane

friction f

matrix friction angle, ϕ = ϕ1

friction-2 f

matrix friction angle, ϕ2

poisson f

Poisson’s ratio, ν

shear f

elastic shear modulus, G

tension f

tension limit, σt. The default is 0.0.

young f

Young’s modulus, E

joint-cohesion f

joint cohesion, cj = cj1

joint-cohesion-2 f

joint cohesion, cj2

joint-dilation f

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

joint-dilation-2 f

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

joint-friction f

joint friction angle, ϕj = ϕj1

joint-friction-2 f

joint friction angle, ϕj2

joint-tension f.

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

normal v

normal direction of the weakness plane, (nx, ny, nz)

normal-x f

x-component of the normal direction to the weakness plane, nx

normal-y f

y-component of the normal direction to the weakness plane, ny

normal-z f

z-component of the normal direction to the weakness plane, nz

flag-bilinear i (a)

= 0 (default) for matrix linear model ;

= 1 for matrix bilinear model.

flag-bilinear-joint i (a)

= 0 (default) for joint linear model ;

= 1 for joint bilinear model ;

< 0 joint effect will be skipped, so that the model is degenerated into a bilinear Mohr-Coulomb model.

flag-brittle b (a)

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

table-cohesion s (a)

name of the table relating matrix cohesion c = c1 to matrix plastic shear strain.

table-cohesion-2 s (a)

name of the table relating matrix cohesion c2 to matrix plastic shear strain.

table-dilation s (a)

name of the table relating matrix dilation angle ψ = ψ1 to matrix plastic shear strain.

table-dilation-2 s (a)

name of the table relating matrix dilation ψ2 to matrix plastic shear strain.

table-friction s (a)

name of the table relating matrix friction ψ = ψ1 angle to matrix plastic shear strain.

table-friction-2 s (a)

name of the table relating matrix friction angle ψ2 to matrix plastic shear strain.

table-joint-cohesion s (a)

name of the table relating joint cohesion cj = cj1 to joint plastic shear strain.

table-joint-cohesion-2 s (a)

name of the table relating joint cohesion cj2 to joint plastic shear strain.

table-joint-dilation s (a)

name of the table relating joint dilation ψj = ψj1 to joint plastic shear strain.

table-joint-dilation-2 s (a)

name of the table relating joint dilation ψj2 to joint plastic shear strain.

table-joint-friction s (a)

name of the table relating joint friction angle ϕj = ϕj1 to joint plastic shear strain.

table-joint-friction-2 s (a)

name of the table relating joint friction angle ϕj2 to joint plastic shear strain.

table-joint-tension s (a)

name of the table relating joint tension limit σtj to joint plastic tensile strain.

strain-shear-plastic f (r)

accumulated plastic shear strain

strain-shear-plastic-joint f (r)

accumulated joint plastic shear strain

strain-tensile-plastic f (r)

accumulated plastic tensile strain

strain-tensile-plastic-joint f (r)

accumulated joint plastic tensile strain

Key

(a) Advanced property.
This property has a default value; simpler applications of the model do not need to provide a value for it.
(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 ν. When choosing the latter, Young’s modulus E must be assigned in advance of Poisson’s ratio ν.
  • 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 joint tension limit used in the model is the minimum of the input σt and cj/tanϕj.
  • The tension table and flag-brittle should not be active at the same time.