FLAC3D Theory and Background • Constitutive Models

Simplified Cap-Yield (CHSoil) Model

As discussed previously in the CYSoil model, when soils are subjected to deviatoric loading, they usually exhibit a decrease in stiffness accompanied by irreversible deformation. The well-known Duncan and Chang model (1970) is commonly used to simulate this hyperbolic stress-strain behavior. The Duncan and Chang model is relatively easy to use. However, as noted in in the CYSoil model, this model has several drawbacks. In addition, hyperbolic relations that rely on nonlinear elasticity are known to have significant limitations (Duncan et al. 1980), including: 1) the relations are applicable prior to failure, but may produce unrealistic behavior of soils at and after failure; 2) the hyperbolic relations do not include volume change resulting from change in shear stress (the implied dilation relies on the Poisson’s ratio), and thus may not be able to predict deformations in dilatant soils accurately, such as dense sands under low confining pressure; and 3) the relations predict an isotropic behavior in the \(\Pi\) plane, which is not always realistic (soil strengths in triaxial compression and extension usually differ in magnitude).

A simplified version of the CYSoil model, named the CHSoil model, is provided as an alternative to the Duncan and Chang model that does not have the drawbacks of this model. The CHSoil model is derived from the same strain hardening/softening logic that exists in the CYSoil model, and therefore can provide a realistic stress-strain relation at failure and post-failure. The CHSoil model provides built-in features and does not have a volumetric cap.

The CHSoil model has three specific features:

1. a built-in friction hardening law that uses hyperbolic model parameters as direct input;
2. a frictional Mohr-Coulomb shear envelope; and
3. two built-in dilation laws — one is based on Rowe’s stress dilatancy theory (Rowe 1962), and the other is a user-defined dilation hardening/softening law.

The unloading behavior of the CHSoil model is elastic. Reloading is elastic up to the outermost yield envelope reached previously. Note that in its present form, the CHSoil model is not intended to simulate cyclic loading.

Incremental Elastic Law

The elastic behavior of the CHSoil model is expressed using Hooke’s law. The incremental expression of the law in terms of principal stress and strain has the form

(1)\[\begin{split}\begin{matrix} \Delta \sigma^{\prime}_1 &= \alpha_1 \Delta e_1^e + \alpha_2 (\Delta e_2^e + \Delta e_3^e) \\ \\ \Delta \sigma^{\prime}_2 &= \alpha_1 \Delta e_2^e + \alpha_2 (\Delta e_1^e + \Delta e_3^e) \\ \\ \Delta \sigma^{\prime}_3 &= \alpha_1 \Delta e_3^e + \alpha_2 (\Delta e_1^e + \Delta e_2^e) \end{matrix}\end{split}\]

where \(\alpha_1 = K^e + 4G^e / 3\), \(\alpha_2 = K^e - 2G^e / 3\), and \(K^e\) and \(G^e\) are the tangent elastic bulk and shear modulus, respectively.

In the CHSoil model, the elastic shear modulus depends on the initial value of mean effective stress, \(p^{\prime}_m\). The mean effective stress is specified by the user, using, for example, the expression

(2)\[ p^{\prime}_m = - {{\sigma^{\prime}_1+\sigma^{\prime}_2+ \sigma^{\prime}_3}\over 3}\]

The variation of \(G^e\) with \(p^{\prime}_m\) is represented by the equation (see, for example, Byrne et al. 2003)

(3)\[ G^e = G_{ref} p_{ref} \Bigl( {p^{\prime}_m \over p_{ref}}\Bigr)^n\]

The parameter \(p_{ref}\) is the reference pressure, and \(G_{ref}\) is the shear modulus number (\(G_{ref} \times p_{ref}\) is the value of the tangent elastic shear modulus at reference pressure). \(n\) is a constant modulus exponent (\(n \le 1\)). Also, the tangent elastic bulk modulus is described by the relation

(4)\[ K^e = K_{ref} p_{ref} \Bigl( {p^{\prime}_m \over p_{ref}}\Bigr)^m\]

The parameter \(K_{ref}\) is the bulk modulus number, the product (\(K_{ref} \times p_{ref}\) is the value of the tangent elastic bulk modulus at reference pressure), and \(m\) is a constant modulus exponent (\(m \le 1\)).

Some useful relations between \(K_{ref}\), \(G_{ref}\), Young’s modulus number, \(E_{ref}\), and Poisson’s ratio at reference pressure, \(\nu_{ref}\), are listed for reference:

(5)\[ K_{ref} = {E_{ref} \over {3 (1 - 2\nu_{ref})}} \qquad G_{ref} = {E_{ref} \over {2 (1 + \nu_{ref})}}\]

(6)\[ {K_{ref} \over G_{ref}} = {{2(1 + \nu_{ref})} \over {3 (1 - 2\nu_{ref})}}\]

The user can specify either \(E_{ref}\) and \(\nu_{ref}\) or \(K_{ref}\) and \(G_{ref}\) as input properties for the model. If \(E_{ref}\) and \(\nu_{ref}\) are specified, \(K_{ref}\) and \(G_{ref}\) are calculated internally from Equation (5). and the resulting values are used with Equations (3) and (4). Values of Poisson’s ratio are restricted to positive values smaller than 0.49. Accordingly, upper and lower bounds for \(K^e\) are specified internally as

(7)\[ {2 G^e \over 3} < K^e < 49.66 G^e\]

As noted above, the initial mean effective pressure, \(p^{\prime}_m\), must be provided as an input property for the model. (This value may be evaluated using FISH and stored in the associated zone property offset.) The value of \(p^{\prime}_m\) is used to calculate the tangent elastic shear and bulk moduli, according to Equations (3) and (4). \(G^e\) and \(K^e\) stay constant in this implementation; the moduli are not updated automatically in terms of mean effective pressure.

Shear Yield Criterion and Flow Rule

Shear yielding is defined by a Mohr-Coulomb criterion. The yield envelope is expressed as

(8)\[ f^s = \sigma^{\prime}_1 - \sigma^{\prime}_3 N_{\phi_m} +2 c \sqrt{N_{\phi_m}}\]

where \(c\) is the cohesion, \(\phi_m\) is the mobilized friction angle, and, by definition,

(9)\[ N_{\phi_m}\ =\ {1 + \sin \phi_m \over {1 - \sin \phi_m}}\]

The potential function is nonassociated and has the form

(10)\[ g = \sigma^{\prime}_1 - \sigma^{\prime}_3 N_{\psi_m}\]

where

(11)\[ N_{\psi_m}\ =\ {1 + \sin \psi_m \over {1 - \sin \psi_m}}\]

and \(\psi_m\) is the mobilized dilation angle. Several different laws are available in the literature to characterize \(\psi_m\).

By default, \(\psi_m\) can be given in terms of plastic shear strain via a user-defined table of \(\psi_m\) versus \(\gamma^p\). The evolution parameter for shear yielding, \(\gamma^p\), is defined incrementally by

(12)\[ \Delta \gamma^p = \Delta \epsilon^{ps} = \sqrt{{(\Delta \varepsilon^{dp}_1)^2} + {(\Delta \varepsilon^{dp}_2)^2} + {(\Delta \varepsilon^{dp}_3)^2} } / \sqrt{2}\]

where \(\Delta \varepsilon^{dp}_i, i\) = 1,3 are the principal, deviatoric, plastic shear-strain increments.

The model parameter flag-dilation must be set to zero to activate this option. If flag-dilation = 0 and no table is provided, it is then assumed that dilation is equal to the input value for ultimate dilation property set with dilation. Note that dilation is not used if a dilation table is provided.

Two built-in dilation laws are also available. A law based on Rowe’s (1962) dilatancy theory is used if flag-dilation = 1. This is the default setting. Alternatively, a simple step function can be used in which \(\psi_m\) is zero for \(\phi_m < \phi_{cv}\), where \(\phi_{cv}\) is a constant specified by the user, and \(\psi_m\) is equal to the ultimate dilation value \(\psi_f\) (set by dilation) for values of mobilized friction larger than \(\psi_{cv}\). The model property flag-dilation must be set to 2 to activate this option. The additional constraint that \(\psi_m\) cannot exceed \(\phi_m\) (to prevent unwanted generation of energy from taking place) is enforced internally by the code.

Tensile Yield Criterion and Flow Rule

The tensile yield function is the same as that used for the CYSoil model, and is of the form

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

The tensile strength, \(\sigma^t\), is given in terms of the plastic tensile-strain measure, \(e^{pt}\), and input by means of a user-defined table. If no table is provided, it is assumed that tensile strength is constant and equal to the input value of the tensile strength property.

The evolution parameter for tensile yielding is the modulus of plastic tensile strain, \(e^{pt}\). The increment of plastic tensile strain is defined as

(14)\[ \Delta \epsilon^{pt} = \Delta \epsilon^{pt}_3\]

where \(\Delta \epsilon_3^{pt}\) is the increment of tensile plastic strain in the direction of the major principal stress (recall that tensile stresses are positive).

Friction Hardening

For most soils, the plot of deviatoric stress versus axial strain obtained in a drained triaxial test can be approximated by a hyperbola. The CHSoil model incorporates a friction strain-hardening law to capture this behavior. In this formulation, the mobilized friction angle, \(\phi_m\), is given in terms of plastic shear strain measure, \(\gamma^p\), by means of the following differential law, similar to the one implemented by Byrne et al. (2003) in their UBCSAND liquefaction model:

(15)\[ d(\sin{\phi_m}) = {G^p \over p^{\prime}_m} d(\gamma^p)\]

where the plastic shear modulus, \(G^p\), is given by

(16)\[ \phi_m \le \phi_f : G^p = G^e \Bigl(1 - {\sin{\phi_m} \over \sin{\phi_f}} R_f \Bigr)^2\]

In this formula, \(G^e\) is the elastic tangent shear modulus, \(\phi_f\) is the ultimate friction angle, and \(R_f\) (the failure ratio) is a constant smaller than 1 (0.9 in most cases) used to assign a lower bound for \(G^p\).

According to Equation (14), the elastic tangent shear modulus is a function of \(p^{\prime}_m\),

(17)\[ G^e = G_{ref} p_{ref} \biggl({p^{\prime}_m \over p_{ref}} \biggr)^n\]

After substitution of Equation (17) into Equation (16), the resulting expression in Equation (15), and rearranging terms, we obtain

(18)\[ d(\gamma^p) = {p^{\prime}_m \over G^e} {d(\sin{\phi_m}) \over \Bigl(1 - {\sin{\phi_m} \over \sin{\phi_f}} R_f \Bigr)^2}\]

Using \(\phi_m\) = 0 at \(\gamma^p\) = 0, integration of this equation gives

(19)\[ \gamma^p = {p^{\prime}_m \over G^e} {\sin{\phi_f} \over R_f} \Bigl( {1 \over {1 - {\sin{\phi_m} \over \sin{\phi_f}} R_f}} - 1 \Bigr)\]

\[\]

Solving for \(\sin{\phi_m}\), we obtain

(20)\[ \sin{\phi_m} = {\sin{\phi_f} \over R_f} \Bigl( 1 - {1 \over {1 + \gamma^p {G^e \over p^{\prime}_m} {R_f \over \sin{\phi_f}}}} \Bigr)\]

This expression is used in the CHSoil model to calibrate mobilized friction in terms of plastic shear strain. The use of this hardening law for modeling primary loading in a triaxial test produces a hyperbolic-type curve of deviatoric stress versus axial strain.

Over-Consolidation

The initial state of a normally consolidated soil or an over-consolidated soil is prescribed by specifying an initial value to the friction property that is equal to or larger than the normally consolidated value, \(\phi_{nc}\). For normal consolidation, Equations (15) and (19) give

(21)\[ \phi_{nc} = \arcsin{ {\sigma^{\prime}_1 - \sigma^{\prime}_3} \over {\sigma^{\prime}_1 + \sigma^{\prime}_3 }}\]

The material behavior is considered to be elastic for stress points below the current yield envelope.

The initial value of friction must be smaller than the ultimate value, \(\phi_f\). (An upper bound equal to \(\phi_f\) is set automatically by the model.) An initial value of the evolution parameter, \(\gamma^p\), that is consistent with the specified initial value of mobilized friction angle is also assigned automatically.

Shear-Induced Compaction and Dilation

A certain amount of irrecoverable volumetric strain, \(\epsilon^p\), is expected to take place as a result of soil shearing. Also, under small monotonic shear strains, there is a tendency for the soil skeleton to contract due to grain rearrangements. For larger shear strains, the soil skeleton may dilate if the soil is dense, as a result of grains riding over each other. A dilation strain-hardening law is incorporated in the logic to model this behavior.

The shear-hardening flow rule implemented in the CHSoil model has the form

(22)\[ \Delta \epsilon^p = \Delta \gamma^p \sin{\psi_m}\]

where \(\Delta \epsilon^p\) is the plastic volumetric strain increment, \(\Delta \gamma^p\) is the plastic shear strain increment, and \(\psi_m\) is the (mobilized) dilation angle.

One possible way to characterize \(\psi_m\) is to adopt an equation based on Rowe’s stress-dilatancy theory (Rowe 1962). According to this theory, there is a (constant-volume) friction angle, \(\phi_{cv}\), below which the material contracts (i.e., for \(\phi_m \le \phi_{cv}\)), while for higher stress ratios (i.e., for \(\phi_m > \phi_{cv}\)), the material dilates. The equation has the form

(23)\[ \sin{\psi_m} = {{\sin{\phi_m} - \sin{\phi_{cv}} \over {1 - \sin{\phi_m}\sin{\phi_{cv}}}}}\]

where

(24)\[ \sin{\phi_{cv}} = {{\sin{\phi_f} - \sin{\psi_f} \over {1 - \sin{\phi_f}\sin{\psi_f}}}}\]

and \(\phi_f\) and \(\psi_f\) are ultimate (known) values of friction and dilation, respectively.

Dilation is evaluated in terms of plastic shear strain based on the last two equations and the assumed relation between \(\phi_m\) and \(\gamma^p\) reported in Equation (19). Rowe’s dilation law is selected by setting the material property flag-dilation = 1 (default value). A simple dilation law also can be used, in which \(\psi_m\) = 0 for \(\phi_m < \phi_{cv}\), and \(\psi_m = \psi_{f}\) for \(\phi_m \ge \phi_{cv}\); the option is activated by specifying flag-dilation = 2. Finally, the user may choose to define an alternative dilation law by creating an input table of dilation values versus plastic shear strain. The model property flag-dilation must be set to zero to activate this option.

Correlation Between CHSoil Properties and Duncan and Chang Properties

A parallel can be drawn, for particular cases, between properties used in the hyperbolic stress-strain relations reported in Duncan et al. (1980) for the Duncan and Chang model and the properties of the CHSoil model. The correlations depend on property input; they are listed in the following table.

Table 1: Property Correlation for Input of \(K_{ref}\) and \(G_{ref}\)

Model	\(m\) = \(n\)	No Dilation
Duncan and Chang	\(K_{ur}\) & \(n\)	\(K_{b}\) & \(m\)
CHSoil \(K_{ref}\) & \(G_{ref}\)	\({9K_{ref}G_{ref}} \over {3K_{ref} + G_{ref}}\) & \(n\)	\(K_{ref}\) & \(m\)
CHSoil \(E_{ref}\) & \(\nu_{ref}\)	\(E_{ref}\) & \(n\)	\({E_{ref}} \over {3(1-2\nu_{ref})}\) & \(m\)

Calibration to Triaxial Test Results on Nevada Sand

Results of drained triaxial compression tests at constant mean effective stress, \(p^{\prime}\), on Nevada sand are used to provide an example of the methodology that can be used to calibrate the properties of the CHSoil model. The test results are listed in data files “40-40.txt”, “40-80.txt”, and “40-160.txt” and correspond to values of the mean stress at 40 kPa, 80 kPa, and 160 kPa, respectively. The relative density of the sand for these tests is 40%. The test results (see the drained triaxial compression tests for more details) consist of three sets of data:

1. deviatoric stress versus mean effective stress;
2. deviatoric stress versus axial strain (up to 25% strain); and
3. volumetric strain versus axial strain (up to 25% strain).

For this calibration exercise, we consider axial strains up to 5%, because strains are not expected to develop above this level for the intended (static, drained) application. The Nevada sand appears to be purely frictional in character. Also, two main features characterize the data: 1) hyperbolic behavior is observed in the plot of deviatoric stress versus axial strain; and 2) bilinear dilatant behavior is exhibited by the plot of volumetric strain versus axial strain. The CHSoil model capabilities, including friction hardening and variable dilation, are used to simulate these features.

The elastic tangent shear and bulk moduli are functions of the initial mean effective stress according to Equations (3) through (5). For friction hardening, we use the built-in law, Equation (19). The dilation law is bilinear; its value is zero below the stress ratio \(\phi=\phi_{cv}\), and a constant, \(\psi_f\), above it:

(25)\[\begin{split}\begin{matrix} \psi &= 0 \qquad \ \ \ \phi < \phi_{cv} \\ \\ \psi &= \psi_f \qquad \phi \ge \phi_{cv} \end{matrix}\end{split}\]

The failure ratio, \(R_f\), is 0.99 for this exercise. Eight properties must be defined:

\[E_{ref}, \nu, p_{ref}, \phi_f, \psi_f, \phi_{cv}, m, n\]

Calibration of the model properties is done in two steps. First, using theoretical considerations and values recorded in the literature, we derive a first estimate for the property values. Second, we improve on the estimates by modeling the triaxial experiments numerically and matching the results obtained in the laboratory.

First Estimates — The value of ultimate friction is derived from a linear fit to plots of deviatoric stress, \(q\), versus mean effective stress, \(p^{\prime}\), obtained from the laboratory tests. The linear fit to maximum \(q\) at given \(p^{\prime}\) provides a value for the maximum stress ratio, \(a\):

(26)\[ a = {q \over p^{\prime}}\]

From the definition of the purely frictional Mohr-Coulomb criterion, we obtain the relation

(27)\[ \sin{\phi_f} = {3a \over 6 + a}\]

from which the ultimate friction angle, \(\phi_f\), can be derived. The estimated value for \(D_r\) = 40% sand is \(\phi_f\) = 34°.

The value of ultimate dilation angle is derived from the slope, \(b\), of a linear fit to the laboratory curves of minus volumetric strain versus axial strain. To relate \(b\) to \(\psi_f\), we use, as a first approximation, the expression for bilinear idealization of triaxial stress results provided by Vermeer and de Borst (1984):

(28)\[ b = {2 \sin{\psi_f} \over 1 - \sin{\psi_f}}\]

The first estimate for \(\psi_f\) is 8.2° for \(D_r\) = 40% sand.

To estimate the parameter \(\phi_{cv}\) for a given \(D_r\), we first derive the maximum value of axial strain, \(\varepsilon^*_a\), at which the volumetric strain is negligible from the laboratory plots of volumetric strain versus axial strain at a given \(p^{\prime}\). From the knowledge of axial strain, \(\varepsilon^*_a\), we estimate \(q\) from the laboratory plot of deviatoric stress versus axial strain at that \(p^{\prime}\), and then the ratio \(q/p^{\prime}\). Three values of the ratio are available (one at each \(p^{\prime}\)). The mean value of \(q/p^{\prime}\) is used to calculate the corresponding friction angle using Equation (27), where \(\phi_f\) now is replaced by \(\phi_{cv}\), and \(a\) by the average \(q/p^{\prime}\) (see Equation (26)). The estimate is \(\phi_{cv}\) = 27.7° for \(D_r\) = 40%.

The first guess for \(E_{ref}\) is taken as 1200 for \(D_r\) = 40% at atmospheric pressure; this value is equivalent to that used in the triaxial numerical experiment for dense sand with the cap-yield (CYSoil) model.

The value of \(\nu\) for this exercise is arbitrarily selected as 0.35. Also, we select \(n = m = 0.5\). The elastic constant, \(E_{ref}\), is estimated by matching the initial slopes of \(q\) versus axial strain curves obtained in similar triaxial tests (under constant mean stress) performed numerically and in the laboratory. These estimates may not be very accurate; more robust estimates for the elastic constants can be obtained from laboratory results of small unloading-reloading excursions. Such results were not available for this example.

cap-yield-simplified model properties

Use the following keywords with the zone property command to set these properties of the simplified Cap-Yield (CHSoil) model.

cap-yield-simplified

bulk-reference f

dimensionless bulk modulus reference value, \(K_{ref}\)

cohesion f

ultimate cohesion, \(c\)

dilation f

ultimate dilation angle, \(\psi_f\)

friction f

ultimate friction angle, \(\phi_f\)

friction-mobilized f

mobilized friction angle, \(\phi_m\), which should be initialized no less than the value determined by this equation via a command or FISH function

poisson f

current Poisson’s ratio, \(\nu\)

pressure-initial f

initial mean effective pressure, \(p_m\), which should be initialized via a command or FISH function

pressure-reference f

reference pressure, \(p_{ref}\). A non-zero value should be specified based on the unit of stress/pressure adopted in the model.

shear-reference f

dimensionless shear modulus reference value, \(G_{ref}\)

tension f

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

young-reference f

dimensionless Young’s modulus reference value, \(E_{ref}\)

dilation-mobilized f

current mobilized dilation angle, \(\psi_m\). It is read-only, except for table-dilation ǂ 0.^(a)

flag-brittle b

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

flag-dilation i

= 0 for mobilized dilation angle, \(\psi_m\), equal to input value, dilation or a function of plastic shear strain if table is input with table-dilation

= 1 (default) for mobilized dilation angle, \(\psi_m\), characterized by Rowe’s stress-dilatancy theory

= 2 for mobilized dilation angle, \(\psi_m\) = 0 if \(\phi_m\) < \(\phi_{cv}\), and \(\psi_m\) = \(\psi_f\), if \(\phi_m \ge \phi_{cv}\)   ^(a)

exponent-bulk f

bulk modulus exponent, \(m\). The default is 0.5.^(a)

exponent-shear f

shear modulus exponent, \(n\). The default is 0.5.^(a)

friction-0 f

initial mobilized fiction angle, \(\phi_0\) associated with zero plastic shear strain. By default, \(\phi_0\) = \(\phi_m\) if strain-shear-plastic is initialized to 0; otherwise, \(\phi_0\) = 0 if not specified.^(a)

friction-critical f

constant used in dilation law (flag-dilation = 1 or 2), \(\phi_{cv}\)   ^(a)

strain-shear-plastic f

accumulated plastic shear strain, \(\epsilon^{ps}\). By default, it is initialized internally.^(a)

strain-tensile-plastic f

accumulated plastic tensile strain, \(\epsilon^{pt}\). The default initial value is 0.^(a)

failure-ratio f

failure ratio, \(R_f\). The default is 0.9.^(a)

table-cohesion s

name of the table relating cohesion, \(c\), to plastic shear strain^(a)

table-dilation s

name of the table relating mobilized dilation angle to plastic shear strain^(a)

table-tension s

name of the table relating mobilized tensile strength to plastic tensile strain^(a)

bulk f

[read only] current elastic bulk modulus, \(K^e\)

shear f

[read only] current elastic shear modulus, \(G^e\)

young f

[read only] current elastic Young’s modulus, \(E\)

Notes:

• Only one of the two options is required to define the elasticity: 1) \(K_{ref}\) and \(G_{ref}\), 2) \(E_{ref}\) and \(\nu\). When choosing the latter, \(E_{ref}\) must be assigned in advance of \(\nu\).
• The tension cut-off is \({\sigma}^t = min({\sigma}^t, c/\tan \phi)\).
• The tension table and flag-brittle should not be active at the same time.

Footnotes

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

Read only properties cannot be set by the user. However, they may be listed, plotted, or accessed through FISH.