FLAC3D Theory and Background • Constitutive Models

Plastic-Hardening Model

The Plastic-Hardening (PH) model is a shear and volumetric hardening constitutive model for the simulation of soil behavior. When subjected to deviatoric loading (e.g., during a conventional drained triaxial test), soils usually exhibit a decrease in stiffness, accompanied by irreversible deformation. In most cases, the plot of deviatoric stress versus axial strain obtained in a drained triaxial test may be approximated by a hyperbola. This feature was discussed by Duncan and Chang (1970) in their well-known “hyperbolic-soil” model, which is formulated as a non-linear elastic model. The PH model is formulated within the framework of hardening plasticity (Schanz et al. 1999) allowing the removal of the main drawbacks of the original non-linear elastic model formulation (e.g., detection of loading/unloading pattern, non-physical bulk modulus).

The main features of the PH model are:

  1. hyperbolic stress-strain relationship during axial drained compression;
  2. plastic strain in mobilizing friction (shear hardening);
  3. plastic strain in primary compression (volumetric hardening);
  4. stress-dependent elastic stiffness according to a power law;
  5. elastic unloading/reloading compared to virgin loading;
  6. memory of pre-consolidation stress; and
  7. Mohr-Coulomb failure criterion.

The model is straightforward to calibrate using either conventional lab tests or in-situ tests. It is well established for soil-structure interaction problems, excavations, tunneling, and settlements analysis, among many other applications.

In this section, principal stress, σi, and strain, εi, i = 1,3, components are positive in tension. Also, all stresses are effective by default if not explicitly stated. The principal effective stresses are with the order σ1σ2σ3, and σ1 is the most compressive stress.

Incremental Elastic Law

The PH model adopts hypo-elasticity for the description of elastic behavior,

(1)Δp=KΔϵev
(2)Δsij=2GΔϵij

where p is the mean pressure defined as p=σii/3, ϵev is the volumetric elastic strain defined as ϵev=εeii, sij=σij+pδij and ϵij=εij(ϵv/3)δij are the deviatoric stress tensor and deviatoric elastic strain tensor, respectively. δij is defined as δij=1 if i=j and δij=0 if ij. K and G are the elastic bulk and shear moduli, which can be derived from the elastic unloading-reloading Young’s modulus, Eur, and the unloading-reloading Poisson’s ratio, ν, are obtained using the relations:

(3)K=Eur3(12ν)
(4)G=Eur2(1+ν)

In the PH model, the Poisson’s ratio is assumed to be a constant with a typical value of 0.2 (if otherwise not provided at input).

The Young’s modulus is a stress-dependent parameter:

(5)Eur=ErefurZm

where

(6)Z=ccotϕσ3ccotϕ+preffcut

and Erefur, m, c, pref, and ϕ are user-defined constant parameters. Erefur is the reference unloading-reloading stiffness modulus at the reference pressure pref. The current unloading-reloading stiffness modulus Eur depends on the maximum (minimum compressive) principal stress, σ3, the cohesion, c, and the ultimate friction angle, ϕ, as well as the power, m. For clays, m is usually close to 1. For sands, m is usually between 0.4 and 0.9. The default cut-off factor fcut is 0.1.

The PH model also employs an additional stiffness measure, E50, which defines the initial slope of the hyperbolic stress-strain curve (see Equation (7) and Figure 1). Parameter E50 obeys the following power law:

(7)E50=Eref50Zm

Here Eref50 is a material parameter, which could be estimated from a set of triaxial compression tests with various cell stresses.

Shear Yield Criterion and Flow Rule

The shear yield function determining the onset and development of shear hardening is defined as

(8)fs=EurEiqaq(qaq)qEurγp2=0

where γp is a shear hardening parameter (one of the internal variables) defined in (13), Ei=2E50/(2Rf), q=σ3σ1, and qa is given as

(9)qa=qfRf=Kqa(ccotϕσ3)

where

(10)Kqa=1Rf2sinϕ1sinϕ

The failure ratio Rf=qf/qa has a value smaller than 1 (typically, Rf = 0.9 is used). Note that the ultimate deviatoric stress, qf, is defined as

(11)qf=2sinϕ(ccotϕσ3)1sinϕ

is consistent with the Mohr-Coulomb failure criterion (Figure 1).

For a standard drained triaxial test, the connection between the axial (vertical compressional) strain, ε1, and deviatoric stress, q, can be described by a hyperbolic relation:

(12)ε1=qaqEi(qaq)

Equation (12) is graphically represented in Figure 1 with the cut-off at qf.


../../../../_images/modelph-1.png

Figure 1: Hyperbolic stress-strain relation in primary shear loading.


The shear hardening parameter γp is defined so that its incremental form is:

(13)Δγp=(Δεp1Δεp2Δεp3)

Due to the increase of γp, the shear yield surface will expand up to the ultimate surface, which is defined by the conventional Mohr-Coulomb failure criterion. In case the material experiences ultimate shear failure, the conventional Mohr-Coulomb failure criterion applies.

The PH model uses the following flow rule between volumetric and shear plastic strains:

(14)Δεpv=sinψmΔγp

where sinψm is the mobilized dilation angle, which should be smaller or equal to the user-defined ultimate dilation angle ψ. The mobilized dilation angle is based on Rowe’s dilation law (1962):

(15)sinψm=Fcsinϕmsinϕcv1sinϕmsinϕcv\;,\;\;if\;sinϕm<sinϕcv
(16)sinψm=sinϕmsinϕcv1sinϕmsinϕcv\;,\;\;if\;sinϕmsinϕcv

where the parameter Fc is a contraction scale factor, with the allowable range of 0 to 0.25 and default value of 0.0 for most soils if the contraction is not considered. The critical state friction angle ϕcv is defined as

(17)sinψcv=sinϕsinψ1sinϕsinψ

The mobilized friction ϕm is defined in terms of the current stress state

(18)sinϕm=σ1σ3σ1+σ32ccotϕ

A non-associated flow rule consistent with the Mohr-Coulomb yield criterion is used in the model. The shear potential function is defined as

(19)gs=m1σ1+m3σ3

where m1=(1+sinψm)/2 and m3=(1+sinψm)/2.

In order to avoid over-dilatancy when the soil reaches its critical void state at emax, the dilation angle needs a minor modification. One way proposed by Schanz et al. (1999) is to set a cut-off rule, so that

(20)sinψcutm=0, if eemax

Here we introduce a smoothing technique to avoid a sudden change of dilation angle:

(21)sinψsmoothm=100(1eemax)sinψm, if e0.99emax

The dilation rules with cut-off, smoothing technique, and without dilation cut-off are schematically compared in Figure 2.


../../../../_images/modelph-2.png

Figure 2: Volumetric strain curve for a standard triaxial compression test with dilation cut-off and smoothing.


Volumetric Cap Criterion and Flow Rule

The volumetric (cap) yield function is defined as

(22)fv=gv=˜q2α2+p2p2c=0

where α is a constant derived internally from other material parameters based on a virtual (numerical) oedometer test, ˜q is a shear stress measure defined as ˜q=[σ1+(δ1)σ2δσ3], and δ=(1+sinϕ)/(1sinϕ). The initial value of hardening parameter pc, which denotes the preconsolidation pressure, can be determined using the initial stress state and over-consolidation ratio, OCR, so that

(23)pc,ini=OCR˜q2iniα2+p2ini

If OCR has a large value, the model behaves as if there is no cap present. The associated flow rule is adopted for volumetric hardening, which means that the potential volumetric function is assumed to be the same as the volumetric yield function (see Equation (24)).

The volumetric hardening parameter γv (one of the internal variables) is defined incrementally as

(24)Δγv=Δεpv=(Δεp1+Δεp2+Δεp3)

where Δεpv is the volumetric plastic strain increment.

Evolution of the hardening parameter, pc, is given by the relation:

(25)Δpc=Hc(ccotϕ+pcccotϕ+pref)mΔγv

where Hc is a constant parameter that can be derived internally from other material parameters. Instead of taking α and Hc as input material parameters, another two parameters, Knc and Erefoed are specified as input. Parameter Knc denotes normal consolidation coefficient and Erefoed stands for the tangent oedometer stiffness at the reference pressure pref. If Knc is not provided by the user, it is taken as Knc=1sinϕ (default value). The PH model reserves the flexibility for the user to input user-defined α and Hc instead of internal calculation for α and Hc. A typical value for α is between 0.6 and 2.0 for most soils. Parameter Hc is calculated as Hc=kKp, where k is a correction factor with a typical value of 0<k1, Kp=K1K2/(K1K2), K1=Erefur/[3(12ν)], and K2=Erefoed(1+2Knc)/3. Specifically, k = 1 if at isotropic compression state.

Tensile Yield Criterion and Flow Rule

The model checks for the tension failure condition. The tension failure and potential functions are

(26)ft=gt=σ3σt

where σt is the tension limit. By default, σt is zero and the user can provide a value with an upper limit of c/tanϕ. The model does not consider tension hardening.

Plastic Corrections

In the implementation of the PH model, elastic trial stresses, σIij, are first computed by adding to the old stresses, σoij, elastic increments computed using the total strain increment Δεij for the step (see Eq. (1.395)). During this step, the moduli and mobilized dilation angle are assumed to be constants (for simplicity) and calculated based on the old stress components. All stresses are assumed to be effective.

The formulation for the trial stresses (denoted by a superscript I) in principal axes are:

(27)σI1=σo1+α1Δε1+α2Δε2+α3Δε3σI2=σo2+α1Δε2+α2Δε3+α3Δε1σI3=σo3+α1Δε3+α2Δε1+α3Δε2

where α1=K+4G/3, α2=K2G/3, K and G are the tangent elastic bulk and shear modulus, respectively, and (Δε1,Δε2,Δε3) is the set of principal strain increments.

The trial internal variables, γp and γv , take their values from the previous step. If the hardening yield criteria (Equation (8) and/or (22)) are exceeded, both the stress components and internal variables will be corrected as described below.

First, consider the situation when only shear hardening occurs, fs>0. The plastic strain increment is oriented in the direction of the gradient of the potential function in the principal stress space. Using Equation (19), we get

(28)Δεps1=λsgsσ1=m1λsΔεps2=λsgsσ2=0Δεps3=λsgsσ3=m3λs

The shear hardening parameter increment is obtained from Equation (13),

(29)Δγp=λs(m1m3)=λs

Thus, the corrected (new) stress components become:

(30)σ1=σI1λs(α1m1+α2m3)σ2=σI2λs(α2m1+α2m3)σ3=σI3λs(α1m3+α2m1)

and the corrected internal variable is

(31)γp=γp,I+Δγp=γp,I+λs

Using Equation (30) and definitions m1 and m2 from Equation (19), updated shear stress and parameter qa (see Equation (9)) are calculated as

(32)q=qI2Gλs
(33)qa=qIa+Kqa(α1m3+α2m1)λs

After substituting the corrected parameters q, qa, γp into Equation (8), parameter λs can be obtained by solving a nonlinear equation fs(λs) = 0 using an iteration method. Note that γp in Equation (31) uses the current value of λs updated in each iteration step, which implies that the correction algorithm is semi-implicit.

If the trial shear stress exceeds the Mohr-Coulomb shear failure surface, the corrections used in Mohr-Coulomb model implementation apply.

Now consider the case when the elastic guess (Equation (22)) exceeds the volumetric yield criterion, fv>0, and the volumetric hardening occurs. The plastic strain increments are related to the gradient of the potential function in the stress space as

(34)Δεpv1=λvgvσ1=[2˜qα223p]λvΔεpv2=λvgvσ2=[2˜qα2(δ1)23p]λvΔεpv3=λvgvσ3=[2˜qα2δ23p]λv

Evolution of the hardening parameter γv is given by the relation:

(35)Δγv=Δεpv=2pλv

In the above formulation, the plastic strain and volumetric hardening parameter increment are related to the current stress measurement (p or ˜q), so the correction algorithm again uses a semi-implicit approach.

The corrected stress components are obtained by analogy with shear hardening,

(36)σ1=σI1λv[2pK+4˜qGα2]σ2=σI2λv[2pK+4(δ1)˜qGα2]σ3=σI3λv[2pK4δ˜qGα2]

Noting that p(I)=(σ1+σ2+σ3)/3 and ˜q(I)=[σ1+(δ1)σ2δσ3], and using Equation (36) we get

(37)p=pI1+2Kλv
(38)˜q=˜qI1+2Mλv

where M=4G(δ2δ+1)/α2. After substituting the corrected p and ˜q into the yield function fv, (Equation (22)), parameter λv can be obtained by solving the nonlinear equation fv(λv) = 0 for the smallest root using an iterative method.

Finally,we consider the situation when the elastic trial stresses exceed both the shear yield hardening surface and volumetric yield hardening surface. In this case, the plastic strain increment is related to the partial derivatives of the potential functions as follows:

(39)Δεp1=λsgsσ1+λvgvσ1=m1λs+[2˜qα223p]λvΔεp2=λsgsσ2+λvgvσ2=[2˜qα2(δ1)23p]λvΔεp3=λsgsσ3+λvgvσ3=m3λs+[2˜qα2δ23p]λv

Parameters γs can γv can be found by using an iterative method similar to that adopted in the shear or volumetric hardening correction technique. If the trial stress is out of both the volumetric hardening surface and Mohr-Coulomb failure surface, the Mohr-Coulomb failure function should be used instead of the shear hardening function.

Implementation Procedure

In the initialization step, the initial stress, evolution parameters γs can γv and strain increment are defined for the zone and are assumed to be constant during this step. For simplicity, the stiffness moduli and dilation, which are dependent on the current stress and evolution parameters, are also assumed to be constant in the zone during this step. The trial elastic stresses and the possible stress corrections are defined at the sub-zone level. The trial stress increments obey the linear elastic (Hooke’s) law. If the trial stresses violate the tension limit, the tension failure procedure is called. The next step is to check whether the tension-corrected stress is in the tension or compression side. If the corrected stress is tensile, volumetric hardening will not apply. In either side, if the stress is out of the shear failure or shear yield surface, the shear failure or hardening procedure will be called, respectively. In the compression side, if the stress is out of the volumetric yield surface, the volumetric hardening procedure will be called. In particular, if the stress is also out of the shear failure or shear yield surface, the mixed procedure with both the volumetric hardening and shear failure/hardening corrections will be called. To ensure the averaged stress in the zone level is within the tension limit, a second tension failure procedure is called if any volumetric hardening or shear hardening/failure occurs. After all sub-zones complete the stress check for tension and shear failure, shear and volumetric hardening criteria, the internal variables, stiffness moduli and dilation are updated based on the zone-averaged stresses.

Comparison of the PH model and CYSoil model

The PH and CYSoil models share many similar features. For example, both models have stress-dependent moduli, both have primary shear and volumetric hardening, both employ a hyperbolic relation for the deviatoric stress versus axial strain during triaxial compression, the volumetric yield and potential functions are the same in both models, etc. However, there are some differences between the models as outlined below.

  • For primary shear, the CYSoil model uses a Mohr-Coulomb type yield function associated with the mobilized friction angle, while the PH model uses the yield function described by Equation (19).
  • The CYmodel uses the accumulated plastic shear strain as one of the evolution parameters (γp), while the definition for γp differs from plastic shear strain in the PH model (see Equation (13)).
  • The CYSoil model moduli are dependent on the mean pressure, while the PH model moduli are dependent on the maximum principal stress (minimum compressive).
  • The PH model introduces a modulus E50 to determine the initial slope of the hyperbolic curve during the triaxial compression, while the CYSoil model does not have such modulus; this implies that these two models may have different hyperbolic curves during the drained triaxial compression, although the ultimate shear stresses are the same.
  • In the CYSoil model, the parameter α is user-defined, while in the PH model, this parameter can be internally calculated using other input parameters.
  • A constant ratio between the elastic and plastic volumetric strain is assumed during isotropic compression at any pressure in the CYSoil model, while the PH model does not adopt this assumption. As a result, the hardening law provided by Equation (25) is used in the PH model.
  • Most parameters of the PH model can be calibrated through the conventional triaxial and oedometer compression tests (see the next section on calibration), while calibration of the CYSoil model parameters is more elaborate.
  • In the CYSoil model, many parameters can be input through tables. Currently, the PH model does not have this flexibility feature.

Calibration of Material Parameters Using Laboratory Data

This section presents information on calibration techniques for the PH model material parameters using the conventional geotechnical laboratory tests. The calibration example uses original test results obtained from the triaxial compression tests on Monterey sand (Lade 1972). The triaxial test consisted of loading the specimen followed by the unloading-reloading regimes. The original test data are provided in Figure 3 and 4. The data are based on three sets of triaxial compression tests with confining pressures of 1.2, 0.6, and 0.3 kgf/cm2. The initial void ratios are 0.783, 0.786, and 0.781, respectively.


../../../../_images/modelph-calibration-1.png

Figure 3: Original q vs. ε1 curves for the confining pressures of 1.2, 0.6, and 0.3 kgf/cm2 of the triaxial compression tests of Monterey sand (Lade 1972).

../../../../_images/modelph-calibration-2.png

Figure 4: Original εv vs. ε1 curves for confining pressures of 1.2, 0.6, and 0.3 kgf/cm2 of the triaxial compression tests of Monterey sand (Lade 1972).

break

Calibration of Friction Angle, ϕ, and Cohesion, c:

The calibration of these two parameter has no difference from the Mohr-Coulomb model.

To estimate values of friction angle ϕ and cohesion c from triaxial test data using the relation q=σ3σ1 vs. p=(σ1+2σ3)/3, follow the procedure outlined below.

  1. Plot deviatoric stress q vs. p using triaxial compression test lab data.
  2. Use a trend line to fit the Mohr-Coulomb envelope.
  3. The slope of the trend line is 6sinϕ/(3sinϕ), which follows from an alternative form of the Mohr-Coulomb failure criterion. This determines friction angle ϕ.
  4. The intercept of the trend line is 6ccosϕ/(3sinϕ), which determines cohesion c, as the friction angle ϕ is already known from the previous step.

Alternatively, these two parameters can be estimated from triaxial test data using the relation ˆq=(σ3σ1)/2 vs. ˆp=(σ1+σ3)/2 by following the procedure outlined below.

  1. Plot deviatoric stress ˆq vs. ˆp using triaxial compression test lab data.
  2. Use a trend line to fit the Mohr-Coulomb envelope.
  3. The slope of the trend line is sinϕ, which follows from an alternative form of the Mohr-Coulomb failure criterion. This determines friction angle ϕ.
  4. The intercept of the trend line is ccosϕ, which determines cohesion c, as the friction angle ϕ is already known from the previous step.

For the Monterey sand, the slope of the trend line of the Mohr-Coulomb envelope using the relation q vs. p, presented in Figure 5, is 1.403. Based on this, the friction angle is calculated as 34.65 degrees. The intercept is approximately zero, which implies that the cohesion is zero.

../../../../_images/modelph-calibration-3.png

Figure 5: Determination of friction angle and cohesion from the triaxial compression test data.

Calibration of Rf, m, and Eref50:

To estimate values of Rf, m, and Eref50 from triaxial test data, follow the procedure outlined below.

  1. Plot the curve of ε1/q vs. ε1/qf using the triaxial compression test data.
  2. Use a trend line to fit the data. The slope of the line is Rf , the intercept is 1/Ei, based on the relation ε1/q=Rf(ε1/qf)+1/Ei (derived from Equation (12) and Rf=qf/qa).

Multiple sets of triaxial compression tests with various confining pressures can be used to produce a set of pairs of (Rf, E50). The final Rf can be estimated as the averaged one, and the pairs of (Ei, σ3) determine m and Eref50 based on Equation (7).

Figure 6 plots the curves of ε1/q vs. ε1/qa using triaxial compression test data of the Monterey sand with confining pressures 1.2, 0.6, and 0.3 kgf/cm2. The slopes of these lines are Rf, and the intercepts are 100/Ei (as strain is given in percentage).

This figure determines three pairs of (Rf, E50), summarized in Table 1. The average value for Rf is 0.957. Parameter pref needs no calibration and its value is assumed to be 0.1 kgf/cm2. Finally, plotting parameters ln(σ3/pref) vs. ln(E50), as shown in Figure 7, and using Equation (7) allows for determination of m and Eref50 (remember that cohesion c = 0 for this example). The slope of the trend line in Figure 7 is m = 0.707 and the intercept determines Eref50 = exp (4.63) = 102.5 kgf/cm2.

Table 1: Determination of Rf, m, and Eref50
σ3 Rf 1/Ei Ei E50 ln(E50) ln(ccotϕσ3ccotϕ+pref)
-0.3 0.9558 0.002459 406.7 212.3 5.358 1.099
-0.6 0.9678 0.001286 777.6 401.3 5.995 1.792
-1.2 0.9476 0.000930 1075.3 565.8 6.338 2.485
break
../../../../_images/modelph-calibration-4.png

Figure 6: Determination of Rf and 1/Ei from three sets of triaxial compression tests with three different confining stresses.

../../../../_images/modelph-calibration-5.png

Figure 7: Determination of m and Eref50 from three sets of triaxial compression tests with three different confining stresses.

break

Calibration of Erefur:

After the calibration of parameter m, it is straightforward to obtain Erefur from Equation (5) using the unloading-reloading modulus Eur obtained from the original q vs. ε1 data in a triaxial compression test. The final value of Erefur can be averaged using data for different confining pressures. For the example of Monterey sand, Erefur is determined to be 320.0 kgf/cm2. If the unloading-reloading moduli are not available, a value in the range of (35)Eref50 can be used for most soils. The PH model uses a default value of 4Eref50 if no input is provided for Erefur.

Calibration of Dilation Angle, ψ:

The ultimate dilation angle ψ can be estimated as sinψ=sm/(sm+2), where sm is the maximum slope (sm) of the εv vs. ε1 curve in the triaxial compression tests. The slope sm is found as: sm=Δεv/Δε1Δεpv/Δεp1=2sinψ/(1sinψ). For this example of the Monterey sand, the dilation angles are 6 ~ 7° (depending on confining pressure). The values are summarized in Table 2.

Calibration of Knc and Erefoed:

Parameters Knc and Erefoed can be calibrated from the oedometer tests. The ultimate value of σ3/σ1 from an oedometer test is Knc. When σ1=pref, the tangent modulus of the σ1 vs. ε1 curve is Erefoed.

Calibration of Poisson’s Ratio, ν:

Poisson’s ratio can be estimated from the initial slope s0 (usually a negative value) of the εv vs. ε1 curve, so that ν=(1+tans0)/2. In the PH model, the default value is 0.2. In this example for Monetary Sand, the Poisson’s ratio is estimated to be 0.3.

Calibration of emax:

The material parameter emax is needed if the dilation smoothing technique (Figure 2) is required. The value of emax can be determined by conducting standard laboratory tests (ASTM D4254). If the standard ASTM D4254 laboratory test data are not available, its value can be estimated through the trial-and-error method based on the curves of the volumetric strain vs. axial strain of the triaxial compression tests. In this example, its value is estimated as 0.803. The default value of emax = 999.0 is used when it is not provided as input, which implies that the dilation smoothing technique will not be activated.

Calibration of Other Parameters:

Such material parameters as σini1, σini2, σini3, eini, and OCR are known initial parameters and should be consistent with the initial conditions. Some parameters can be estimated through the trial-and-error method.

A summary of all material properties determined for the Monterey sand is provided in Table 2.


Table 2: Calibrated Material Parameters for Loose Monterey Sand
Parameters σ3 = -1.2 (kgf/cm2) σ3 = -0.6 (kgf/cm2) σ3 = -0.3 (kgf/cm2)
Eref50 (kgf/cm2)   102.5  
Erefur (kgf/cm2)   320.0  
m   0.707  
Rf   0.957  
pref (kgf/cm2)   0.1  
ν   0.3  
ϕ (degrees)   34.65  
c (kgf/cm2)   0.0  
emax   0.803  
Knc   0.5  
ψ (degrees) 6.1 6.4 7.0
σini13 (kgf/cm2) -1.2 -0.6 -0.3
OCR 1.0 2.0 4.0
eini 0.783 0.786 0.781

Using these parameters, the triaxial compression tests can be simulated by the PH model. The plots presented in Figures 8 and 9 reveal a close match between the simulated results and laboratory test data.

break
../../../../_images/modelph-calibration-6a.png

Figure 8: Deviatoric stress vs. axial strain for consolidated drained triaxial compression tests on fine Monterey sand.

../../../../_images/modelph-calibration-6b.png

Figure 9: Volumetric strain vs. axial strain for consolidated drained triaxial compression tests on fine Monterey sand.

Small-Strain Stiffness Option

Using the Plastic-Hardening model to simulate the unloading-reloading paths, the reloading path will be the reverse of the unloading path. However, this elastic behavior assumes of very small strain. Experiments illustrate that soil stiffness modulus is strain-dependent – it decreases nonlinearly with the increase of shear strain. The strain can be categorized into three levels: the very small strain, where the stiffness modulus is assumed constant (G0 or sometimes called Gmax); the small strain, where the stiffness modulus reduces non-linearly as the strain increases; and the large strain, where the soil is in or close to failure and the stiffness modulus is relatively small. Figures 10 illustrates this explanation using the normalized stiffness modulus degradation curve by comparing with the geotechnical applications and the measurement accuracy from the laboratory or in-situ tests. It is important that the soil stiffness modulus selected for simulation should not be associated to the strain levels at the end of the geotechnical constructions or from the laboratory tests with well disturbed samples, because the soil stiffness modulus could have decreased much from its initial value. The additional option of the small-strain stiffness modulus for the Plastic-Hardening model is objected to take the strain-dependency into account.

../../../../_images/modelph-strainrange.png

Figure 10: Stiffness modulus degradation curve and typical strain ranges (Modified from Atkinson and Sallfors 1991 and Ishihara 1996).

To activate the small-strain stiffness in Plastic-Hardening model, flag-smallstrain (the default is off) should be turned into on. Only two extra properties are required to be assigned: stiffness-0-reference (Eref0) and strain-70 (γ70).

The initial stiffness modulus at the reference pressure, Eref0, will determine the shear stiffness modulus (G0) at the very small strain level. It can be estimated by the relations of

(40)E0=Eref0Zm

and

(41)G0=E02(1+ν)

where Z, m and ν are defined in the standard Plastic-Hardening model. The default value of Eref0 is taken as 3Erefur if no specific data is available. Eref0 is not allowed to be equal or smaller than Erefur in this model. Together with the above two equations, Eref0 could be estimated from

(42)G0=ρV2s

where ρ is the soil density and Vs is the measured in-situ soil shear wave velocity. G0 can be estimated from other in-situ tests data, e.g., SPT, CPT & DMT, or from some accepted empirical relations.

The reference shear strain (γ70) is selected the value at which the secant shear stiffness modulus is about 70% (accurately should be 72.2%) of G0. The default value for γ70 is 2e-4. Typical vales of γ70 is listed in Table 3:

Table 3: Typical vales of γ70
Sand:  
At depth 0 to 6 m 1.1e-4
At depth 6 to 15 m 1.6e-4
At depth 15 to 37 m 2.3e-4
At depth 37 to 76 m 3.3e-4
At depth 76 to 150 m 4.5e-4
At depth 150 to 300 m 7.0e-4
Clay:  
Plastic Index = 30 3.5e-4
Plastic Index = 50 7.7e-4
Plastic Index = 100 1.8e-3

Once G0 and γ70 are known, the small-strain secant stiffness modulus is calculated as

(43)Gs=G01+0.385γγ70

and its tangent stiffness modulus is

(44)G=G0(1+0.385γγ70)2

The tangent shear stiffness modulus is setting a lower-bound value of Gur=Eur/[3(12ν)] in this model.

The detailed algorithm to consider the unloading-reloading strain-dependent stiffness modulus is presented in Benz (2007). It should be noted that this algorithm may not preform well on nested hysteretic loops and is subjected to the so-called over-shooting problem. The plastic-Hardening model with small-strain stiffness is not recommended for fully dynamic problems with random unloading-reloading paths which leads to nested hysteretic loops.

References

Atkinson, J. H., & Sallfors, G. Experimental determination of stressstrain-time charactristics in laboratory and in situ tests. Proceedings, 10th ECSMFE, Florence, Rotterdam: Balkema, 959-999 (1991).

Benz, T. Small-strain stiffness of soils and its numerical consequences (Vol. 5). Stuttgart: Univ. Stuttgart, Inst. f. Geotechnik (2007).

Duncan, J. M., and C. Y. Chang. “Nonlinear Analysis of Stress and Strain in Soils,” Soil Mechanics, 96 (SM5), 1629-1653 (1970).

Ishihara, K. Soil behaviour in earthquake geotechnics. Clarendon Press (1996).

Schanz, T., Vermeer, P. A., & Bonnier, P. G. The hardening soil model: formulation and verification. Beyond 2000 in computational geotechnics, 281-296 (1999).

plastic-hardening Model Properties

Use the following keywords with the zone property command to set these properties of the plastic hardening model.

plastic-hardening
cohesion f

cohesion, c. The default and cut-off (allowable minimum) value is 105pref.

dilation f

ultimate dilation angle in degrees, ψ. The default is 0.0.

exponent f

exponent for elastic moduli, m. Allowable value is 0m0.999.

friction f

ultimate friction angle in degrees, ϕ. The default and cut-off (allowable minimum) value is 0.001 degrees.

stress-1-effective f

initial minimum principal effective stress, σini1. See the note*.

stress-2-effective f

initial median principal effective stress, σini2. See the note*.

stress-3-effective f

initial maximum principal effective stress, σini3. See the note*.

pressure-reference f

reference pressure, pref . A non-zero value should be specified based on the unit of stress/pressure adopted in the model. Most commonly used and highly recommended reference pressure is the standard atmospheric pressure, e.g., 100.0 if in kPa, 1.0e5 if in Pa. or other value compatible to the adopted stress/pressure unit.

stiffness-50-reference f

secant stiffness, Eref50, at 50% of the ultimate deviatoric stress, qf, when σ3=pref in a triaxial test. A non-zero value must be specified.

stiffness-0-reference f

[advanced] Initial stiffness, Eref0, at zero strain when σ3=pref. It is required that Eref0>Erefur. If not input or by default, Eref0=3Erefur. Need to be assigned only when the flag for small-strain stiffness is on, otherwise its value is useless.

stiffness-oedometer-reference f

[advanced] tangent stiffness, Erefoed, when σ1=pref in an oedometer test. By default, Erefoed=Eref50 if Erefoed is not specified while Eref50 has been specified.

stiffness-ur-reference f

[advanced] unloading-reloading stiffness, Erefur, when σ3=pref in a triaxial test. This value must be greater than 2Eref50. By default, Erefur=4Eref50.

coefficient-normally-consolidation f

[advanced] normal consolidation coefficient, Knc. It is not allowed to be less than ν/(1ν), a range 0.5 to 0.7 is common. The default is Knc=1sinϕ.

constant-alpha f

[advanced] dimensionless parameter in the cap yield function, α. Initialized internally. See the note* for an alternative option.

factor-contraction f

[advanced] contraction factor, Fc, allowable range is 0.0 to 0.25. The default is 0.0.

factor-cut f

[advanced] cut-off factor, fcut. The default is 0.1.

flag-initialization f

[advanced] flag for initialization. If set to 1, the internal variables will be recalculated. The default is 0.

flag-smallstrain f

[advanced] flag for small-strain stiffness . If set to on, the option of small-strain stiffness will be activated. The default is off.

over-consolidation-ratio f

[advanced] over consolidation ratio, OCR. For normally consolidated soil, the value should be input as 1.0. The default is 100.0, which denotes an “apparent” non-cap shear plastic hardening model.

poisson f

[advanced] Poisson’s ratio, ν. The default is 0.2.

stiffness-cap-hardening f

[advanced] hardening modulus for cap hardening, Hc. Initialized internally. See the note* for an alternative option.

strain-70 f

[advanced] Shear strain at which Gs/G0 = 72.2%, If not input or by default, γ70=2×104. Need to be assigned only when the flag for small-strain stiffness is on, otherwise its value is useless.

failure-ratio f

[advanced] failure ratio, Rf. The default is 0.9.

tension f

[advanced] tension limit, σt. The default is 0.0.

void-initial f

[advanced] initial void ratio, eini. If the current void ratio is not important, or you do not wish to activate the dilation smoothing technique, just use the default value. The default is 1.0.

void-maximum f

[advanced] allowable maximum void ratio, emax . The default is 999.0, which is not physical but it denotes that the dilation smoothing technical will nominally not activated.

bulk f

[read only] current loading-unloading bulk modulus, K

shear f

[read only] current unloading-reloading shear modulus, G

plastic-hardening-shear f

[read only] accumulated shear plastic hardening parameter, γp

plastic-hardening-volume f

[read only] accumulated volumetric plastic hardening parameter, εpv

pressure-cap f

[read only] cap (preconsolidation) pressure, pc. Determined by the current stress, α, and OCR.

void f

[read only] current void ratio, e

Notes:
  • The tension cut-off is σt=min(σt,c/tanϕ).
  • The parameters σini1, σini2, σini3 must be specified via commands or a FISH function the first time the PH model is assigned because it is a stress-dependent model. Negative values denote compressive stresses.
  • Some parameter combinations may be rejected due to: 1) one or more parameters are out of predetermined limits; 2) internal parameters cannot be correctly determined; or 3) numerical instability.
  • If both α and Hc are specified non-zero values by the user, the PH model will use these user-defined parameters α and Hc, and the internal calculation for α and Hc will not be performed.

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.