FLAC3D Theory and Background • Constitutive Models

Double-Yield Model

Permanent volume changes caused by the application of isotropic pressure are taken into account in this model by including, in addition to the shear and tensile failure envelopes in the strain-softening/hardening model, a volumetric yield surface (or “cap”). For simplicity, the cap surface, defined by the “cap pressure” pc> 0, is independent of shear stress; it consists of a vertical line on a plot of shear stress versus mean stress. The hardening behavior of the cap pressure is activated by volumetric plastic strain, and follows a piecewise-linear law prescribed in a user-supplied table. The tangential bulk and shear moduli evolve as plastic volumetric strain takes place according to a special law defined in terms of a factor, R, assumed to be constant, and defined as the ratio of elastic bulk modulus to plastic bulk modulus.

Formulations

Only two additional material parameters and a table are required, in addition to those associated with the strain-softening model:

  1. the initial value of pc, which corresponds to the maximum mean pressure that the material has experienced in the past;
  2. the value of R, greater than unity, which controls the slope of the stress-strain curve on volumetric unloading (the “swelling” line, in soil mechanics terms); and
  3. the table representation of the “hardening curve,” which relates cap pressure, pc, to plastic volume strain, epv.

Hence, any laboratory-determined hardening behavior may be modeled within the constraints imposed by a two-parameter model.

Incremental Elastic Law

In the FLAC3D/3DEC implementation of this model, principal stresses σ1, σ2, σ3 are used. The principal stresses and principal directions are evaluated from the stress tensor components, and ordered so that (recall that compressive stresses are negative)

(1)σ1σ2σ3

The corresponding principal strain increments Δe1,Δe2,Δe3 are decomposed:

(2)Δei=Δeei+Δepii=1,3

where the superscripts e and p refer to elastic and plastic parts, respectively, and the plastic components are nonzero only during plastic flow. (Note that extensional strains are positive.) It is assumed that the plastic contributions of shear, tensile, and volumetric yielding are additive, so we may write

(3)Δepi=Δepsi+Δepti+Δepvi

where the superscripts ps, pt, and pv stand for plastic shear, plastic tensile, and plastic volumetric strain. By convention, in this section, the symbol Δe is used to refer to the minus volumetric strain increment (Δe1+Δe2+Δe3) with plastic part Δep and elastic part Δee. The symbol Δepv refers to minus the value of the plastic volumetric strain (Δepv1+Δepv2+Δepv3).

The incremental expression of Hooke’s law in terms of principal stress and strain has the form

(4)Δσ1=α1Δee1+α2(Δee2+Δee3)Δσ2=α1Δee2+α2(Δee1+Δee3)Δσ3=α1Δee3+α2(Δee1+Δee2)

where α1=Kc+4Gc/3, α2=Kc2Gc/3, and Kc and Gc are the current tangential bulk and shear moduli defined according to the following considerations.

Consider an isotropic compression test with increasing pressure, pc. As the material becomes more compact, its plastic stiffness (dpc/depv) usually increases. It seems reasonable that the elastic stiffness will also increase, since the grains are being forced closer together. A simple rule is adopted in this model whereby, under general loading conditions, the incremental elastic stiffness, Kc, is a constant factor, R, multiplied by the current incremental plastic stiffness. The values of bulk and shear modulus, K and G, supplied by the user are taken as upper limits to Kc and Gc, and it is assumed that the ratio Kc/Gc remains constant and equal to K/G. Using incremental notation, this law is defined by the relations

(5)Kc=RΔpcΔepvKc:=min(Kc,K)
(6)Gc=GKcK

where the factor R is given, and Δpc/Δepv is the current slope of the table of pc values.

The type of behavior exhibited by the double-yield model is illustrated in Figure 1, which shows a nonlinear volumetric loading curve with several unloading excursions; these excursions are elastic, with slope related by R to the plastic stiffness at the point of unloading.


../../../../_images/modeldoubley-1.png

Figure 1: Elastic volumetric loading/unloading paths.


Yield and Potential Functions

The shear and tensile yield functions, referred to as fs and ft, have the form

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

where

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

and ϕ is the friction angle, c is the cohesion, and σt is the tensile strength.

The volumetric yield function fv is defined as

(10)fv=13(σ1+σ2+σ3)+pc

where pc is the cap pressure.

The shear potential function, gs, corresponds to a nonassociated flow rule; the tensile and volumetric potential functions, gt and gv, correspond to associated laws. They have the form

(11)gs=σ1σ3Nψ
(12)gt=σ3
(13)gv=13(σ1+σ2+σ3)

where

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

and ψ is the dilation angle.

Hardening/Softening Parameters

The shear and volume yield surfaces can harden (or soften), and the tensile yield surface can soften, according to hardening rules that are specified by user-defined tables. Entry to the tables is by hardening parameters that record some measure of accumulated plastic strain. In shear and tension, the hardening parameter incremental forms are

(15)Δeps={ 12(Δeps1Δepsm)2+12(Δepsm)2+12(Δeps3Δepsm)2 }12
(16)Δept=ept3

where

(17)Δepsm=13(Δeps1+Δeps3)

Δepsj,j = 1,3, and Δept3 are plastic shear and tensile strain increments in the principal directions.

In the volumetric direction, the hardening parameter increment is

(18)Δepv=|Δepv1+Δepv2+Δepv3|

where Δepvj,j = 1,3 are plastic volumetric strain increments in the principal directions.

These hardening parameters are used in the tables to determine new values of friction, cohesion, dilation, tensile strength, and cap pressure. The current bulk and shear moduli are also calculated from the table values as per Equations (5) and (6).

Plastic Corrections

Let the superscript I be used to represent the elastic guess obtained by adding to the old stresses σOij, elastic increments computed using the total strain increments. In principal axes, we then have

(19)σI1=σO1+α1Δe1+α2(Δe2+Δe3)σI2=σO2+α1Δe2+α2(Δe1+Δe3)σI3=σO3+α1Δe3+α2(Δe1+Δe2)

In the FLAC3D/3DEC implementation, shear yield is detected if fs(σI1,σI3)< 0, volumetric yield if fv(σI1,σI2,σI3)< 0, and tensile yield if ft(σI3)< 0. Corresponding plastic corrections are evaluated using the following techniques.

We first consider the case where tensile failure is not detected for the step, but both shear and volumetric yield conditions are exceeded. Using Equations (2) and (3), the principal strain increments may be expressed as

(20)Δei=Δeei+Δepsi+Δepvii=1,3

The flow rules for shear and volumetric yielding are

(21)Δepsi=λsgsσi
(22)Δepvi=λvgvσi

where i = 1,3. Using Equations (11) to (13), these expressions become, after differentiation:

(23)Δeps1=λsΔeps2=0Δeps3=λsNψ

and

(24)Δepv1=λv/3Δepv2=λv/3Δepv3=λv/3

Substituting in Equation (15), we obtain

(25)Δee1=Δe1λsλv/3Δee2=Δe2λv/3Δee3=Δe3+λsNψλv/3

With these expressions for the elastic strain increments, Hooke’s incremental equations yield (see Equation (4))

(26)σN1=σI1λs(α1α2Nψ)λvKσN2=σI2α2λs(1Nψ)λvKσN3=σI3λs(α2α1Nψ)λvK

where σIi, i = 1,3, are the initial trial stresses in Equation (19), and σNi=σOi+Δσi, i = 1,3 are the new principal stresses for the step.

To determine the multipliers λs and λv, we require that if shear and volumetric yielding occur, the new stresses lie on both yield surfaces and we must have fs(σN1,σN3) = 0 and fσ(σN1,σN2,σN3) =0. Substituting Equation (25) for σi,i = 1,3 in Equations (7) and (10), and solving for λs, we obtain

(27)λs=fsIfvI(1Nϕ)α1α2Nψα2Nϕ+α1NϕNψK(1Nϕ)(1Nψ)

Hence,

(28)λv=fvIKλs(1Nψ)

In these equations, the notation fI stands for the function f evaluated for the initial trial stresses.

Equations (27) and (28) can now be used to evaluate the new stresses from Equation (26); these stresses simultaneously satisfy both yield conditions and both flow rules.

If the element is only yielding in shear, then

(29)λv=0
(30)λs=fsIα1α2Nψα2Nϕ+α1NϕNψ

If the element is only yielding in volume, then

(31)λs=0
(32)λv=fvIK

Equation (29) to (32) may be used in Equation (26), as appropriate, to compute new stresses.

We now consider the case where tensile failure is detected by the condition ft(σI3)< 0. If volumetric failure is not detected, we use the same technique and stress corrections as described in the Mohr-Coulomb model. If volumetric failure is detected in addition to tensile failure, then either fs(σI1,σI3) 0 or fs(σI1,σI3)> 0. We begin by assuming that all three yield conditions are exceeded. We assume that the plastic contributions of shear, volumetric, and tensile yielding are additive:

(33)Δei=Δeei+Δepsi+Δepvi+Δeptii=1,3

The flow rule for tensile yielding has the form

(34)Δepti=λtgtσii=1,3

Using Equation (8), these expressions become, after partial differentiation:

(35)Δept1=0Δept2=0Δept3=λt

Using the same reasoning as above, and Equations (23) and (24) for the shear and volumetric flow rule, we obtain

(36)σN1=σI1λs(α1α2Nψ)λvK+λtα2σN2=σI2α2λs(1Nψ)λvK+λtα2σN3=σI3λs(α2α1Nψ)λvK+λtα1

The multipliers λs, λv, and λt are determined by solving the system of three equations, fs(σN1,σN3) = 0, fv(σN1,σN2,σN3) = 0, and ft(σN3) = 0. This gives

(37)λs=ftI(1+2Nϕ)+3fvI2fsIα2α1λv=fvIK+3(1+Nϕ)ftI6fvI+3fsIα2α1λt=ftI[Nψ(1+2Nϕ)+2+Nϕ]3(1+Nψ)fvI+(1+2Nψ)fsIα2α1

Substitution of those expressions in Equation (36) yields

(38)σN1=σtNϕ2cNϕσN2=3pcσt(1+Nϕ)+2cNϕσN3=σt

If only tensile and volumetric yield are detected, then λs = 0 in Equation (37). The constants λv and λt are determined by requiring that conditions fv(σN1,σN2,σN3) = 0 and ft(σN3) = 0 both be fulfilled. After some manipulation, we obtain

(39)λv=α1fvI+KftIK(α1K)λt=fvI+ftIα1K

Substitution of those expressions in Equation (36) gives

(40)σN1=σI13fvI+ftI2σN2=σI23fvI+ftI2σN3=σt

Implementation Procedure

Hardening and softening behaviors for the cohesion, friction, and dilation in terms of the shear parameter eps (see Equation (15)) are provided by the user in the form of tables. Softening of the tensile strength is described in a similar manner using the parameter ept (see Equation (16)). In turn, the variation of cap pressure is specified in a table in terms of the parameter epv (see Equation (18)). Each table contains pairs of values: one for the parameter, and one for the corresponding property value. It is assumed that the property varies linearly between two consecutive parameter entries in the table.

In the implementation of the double-yield model in FLAC3D/3DEC, new stresses for the step are computed using the current values of the model properties. In this process, an elastic guess, σIij, is first computed by adding to the old stress components, increments calculated by application of Hooke’s law to the total strain increment for the step. Principal stresses σI1,σI2,σI3, and corresponding principal directions, are calculated and ordered. If these stresses violate the composite yield criterion, corrections are applied to the elastic guess to give the new stress state. The stress tensor components in the system of reference axes are then calculated from the principal values by assuming that the principal directions have not been affected by the occurrence of a plastic correction.

Plastic strain increments are evaluated from Equations (23), (24), and (35), using relevant expressions of λs, λt and λv for the mode of failure taking place. Zone hardening increments are then calculated as the surface average of values obtained from Equations (15), (18), and (16) for all triangles involved in the zone. The hardening parameters are updated, and new zone properties for cohesion, friction, dilation, tensile strength, and cap pressure are evaluated by linear interpolation in the tables. New elastic constants are derived from the cap pressure table using Equations (36). All these properties are stored for use in the next step. The hardening or softening lags one timestep behind the corresponding plastic deformation. In an explicit code, this error is small because the steps are small.

For a material with friction, the maximum value of the tensile strength is evaluated from c/tanϕ, using the new cohesion and friction angle. This value is retained by the code if it is smaller than the tensile strength updated from the table.

Choice of Volumetric Properties

The “hardening curve” and ratio, R, of elastic bulk modulus to plastic bulk modulus are volumetric properties that may be derived from the results of a triaxial test in which axial stress and confining pressure, p, are kept equal. This test (for which dep=depv) is recommended because it is best to determine the parameters related to a particular mode of failure from a test that only involves that failure mode.


../../../../_images/modeldoubley-2.png

Figure 2: Isotropic consolidation test.


Consider the experimental graph of minus mean stress (pressure) versus minus volumetric strain for an increasing stress level with a small unloading excursion, obtained from such a test and presented in Figure 2. The volumetric strain increment de at a point of the main loading path (assuming that we are above any initial pre-consolidation stress level) is composed of an elastic part, dee, and a plastic part, dep. (Recall that, in this section, de, dee, and dep refer to minus the value of the volumetric strain.) The observed tangent modulus may be expressed as

(41)dpde=hKch+Kc

where h is the plastic modulus, Kc is the elastic modulus, and, by definition:

(42)h=dpcdep
(43)Kc=dpdee

With the preceding notation convention, the volumetric property, R, may be defined as

(44)R=Kc/h

and Equation (41) becomes

(45)dpde=Kc1+R

Expressing R from this relation, we obtain

(46)R=Kcdp/de1

Values for dp/de and Kc can be estimated from main loading and unloading increments on the graph. Hence, R can be calculated from Equation (44). Note that in the context of this model, the ratio R is assumed to be constant. Using Kc=Rh and h=dpc/dep in Equation (41), we can write, after some manipulation,

(47)dpcdep=h=(1+RR)dpde

From this, it follows that values of pc for a particular ep can be obtained to the first approximation by multiplying the value p on the graph corresponding to e=ep by the ratio (1+R)/R. For example, if R = 5, then the graph curve must be scaled by a factor of 1.2 to convert it to table values, assuming no over-consolidation.

The double-yield model was initially developed to represent the behavior of mine backfill material, for which pre-consolidation pressures are low. The modified Cam-clay model is more applicable to soils such as soft clays for which pre-consolidation pressures can have a significant effect on material behavior.

Comparison of the Double-Yield and Modified Cam-Clay Model

In the double-yield model:

  • Elastic moduli remain constant during elastic loading and unloading.
  • Shear and tensile failure are not coupled to plastic volumetric change, due to volumetric yielding. The shear yield function corresponds to the Mohr-Coulomb criterion, and the tensile yield is evaluated based on a tensile strength.
  • Material hardening or softening, upon shear or tensile failure, is defined by tables relating friction angle and cohesion to plastic shear strain, and tensile strength to plastic tensile strain.
  • Volumetric yielding occurs; the isotropic stress at which yield occurs (the cap pressure) is not influenced by the amount of shear or tensile plastic deformation, but it increases (or hardens) with volumetric plastic strain.
  • A tensile strength limit and tensile softening can be defined.

In the modified Cam-clay model:

  • The elastic deformation is nonlinear, with the elastic moduli depending on mean stress.
  • hear failure is affected by the occurrence of plastic volumetric deformation: the material can harden or soften depending on the degree of pre-consolidation.
  • As shear loading increases, the material evolves toward a critical state at which unlimited shear strain occurs with no accompanying change in specific volume or stress.
  • There is no resistance to tensile mean stress.


double-yield model properties

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

double-yield
bulk-maximum f

maximum elastic bulk modulus, Kmax

cohesion f

cohesion, c

dilation f

dilation angle, ψ. The default is 0.0.

friction f

angle of internal friction, ϕ

multiplier f

multiplier on current plastic cap modulus to give elastic bulk and shear moduli, R. The default is 5.0.

pressure-cap f

current intersection of volumetric yield surface with pressure axis, pc

shear-maximum f

maximum elastic shear modulus, Gmax

tension f

tension limit, σt. The default is 0.0.

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 cohesion to plastic shear strain. The default is a null string.

table-dilation s (a)

name of the table relating dilation angle to plastic shear strain. The default is a null string.

table-friction s (a)

name of the table relating friction angle to plastic shear strain. The default is a null string.

table-pressure-cap s (a)

name of the table relating cap pressure to plastic volume strain. The default is a null string.

table-tension s (a)

name of the table relating tensile limit to plastic tensile strain. The default is a null string.

bulk f (r)

current elastic bulk modulus, K

poisson f (r)

current Poisson’s ratio, ν

shear f (r)

current elastic shear modulus, G

strain-shear-plastic f (r)

accumulated plastic shear strain

strain-tensile-plastic f (r)

accumulated plastic tensile strain

strain-volumetric-plastic f (r)

accumulated plastic volumetric strain

young f (r)

current elastic Young’s modulus, E

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

  • If table-pressure-cap is not assigned, the model moduli will be constant, so that K = Kmax and G = Gmax. In any case, Kmax and Gmax are upper limits of moduli.
  • The tension cut-off is σt=min(σt,c/tanϕ).
  • The tension table and flag-brittle should not be active at the same time.