FLAC3D Theory and Background • Constitutive Models
PlasticHardening (PH) Model*
Note
*Not available in 3DEC.
The PlasticHardening (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 wellknown “hyperbolicsoil” model, which is formulated as a nonlinear 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 nonlinear elastic model formulation (e.g., detection of loading/unloading pattern, nonphysical bulk modulus).
The main features of the PH model are:
hyperbolic stressstrain relationship during axial drained compression;
plastic strain in mobilizing friction (shear hardening);
plastic strain in primary compression (volumetric hardening);
stressdependent elastic stiffness according to a power law;
elastic unloading/reloading compared to virgin loading;
memory of preconsolidation stress; and
MohrCoulomb failure criterion.
The model is straightforward to calibrate using either conventional lab tests or insitu tests. It is well established for soilstructure interaction problems, excavations, tunneling, and settlements analysis, among many other applications.
Formulations
In this section, principal stress, \(\sigma_i\), and strain, \(\varepsilon_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 \(\sigma_1 \le \sigma_2 \le \sigma_3\), and \(\sigma_1\) is the most compressive stress.
Incremental Elastic Law
The PH model adopts hypoelasticity for the description of elastic behavior,
where \(p\) is the mean pressure defined as \(p = − \sigma_{ii}/3\), \(\epsilon^e_v\) is the volumetric elastic strain defined as \(\epsilon^e_v = \varepsilon^e_{ii}\), \(s_{ij} = \sigma_{ij} + p \delta_{ij}\) and \(\epsilon_{ij} = \varepsilon_{ij}  (\epsilon_v/3)\delta_{ij}\) are the deviatoric stress tensor and deviatoric elastic strain tensor, respectively. \(\delta_{ij}\) is defined as \(\delta_{ij} = 1\) if \(i = j\) and \(\delta_{ij} = 0\) if \(i \ne j\). \(K\) and \(G\) are the elastic bulk and shear moduli, which can be derived from the elastic unloadingreloading Young’s modulus, \(E_{ur}\), and the unloadingreloading Poisson’s ratio, \(\nu\), are obtained using the relations:
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 stressdependent parameter:
where
and \(E^{ref}_{ur}\), \(m\), \(c\), \(p^{ref}\), and \(\phi\) are userdefined constant parameters. \(E^{ref}_{ur}\) is the reference unloadingreloading stiffness modulus at the reference pressure \(p^{ref}\). The current unloadingreloading stiffness modulus \(E_{ur}\) depends on the maximum (minimum compressive) principal stress, \(\sigma_3\), the cohesion, \(c\), and the ultimate friction angle, \(\phi\), 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 cutoff factor \(f_{cut}\) is 0.1.
The PH model also employs an additional stiffness measure, \(E_{50}\), which defines the initial slope of the hyperbolic stressstrain curve (see Equation (7) and Figure 1). Parameter \(E_{50}\) obeys the following power law:
Here \(E^{ref}_{50}\) 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
where \(\gamma^p\) is a shear hardening parameter (one of the internal variables) defined in (13), \(E_i = 2 E_{50} / (2 − R_f)\), \(q = \sigma_3 − \sigma_1\), and \(q_a\) is given as
where
The failure ratio \(R_f = q_f / q_a\) has a value smaller than 1 (typically, \(R_f\) = 0.9 is used). Note that the ultimate deviatoric stress, \(q_f\), is defined as
is consistent with the MohrCoulomb failure criterion (Figure 1).
For a standard drained triaxial test, the connection between the axial (vertical compressional) strain, \(\varepsilon_1\), and deviatoric stress, \(q\), can be described by a hyperbolic relation:
Equation (12) is graphically represented in Figure 1 with the cutoff at \(q_f\).
The shear hardening parameter \(\gamma^p\) is defined so that its incremental form is:
Due to the increase of \(\gamma^p\), the shear yield surface will expand up to the ultimate surface, which is defined by the conventional MohrCoulomb failure criterion. In case the material experiences ultimate shear failure, the conventional MohrCoulomb failure criterion applies.
The PH model uses the following flow rule between volumetric and shear plastic strains:
where \(\sin \psi_m\) is the mobilized dilation angle, which should be smaller or equal to the userdefined ultimate dilation angle \(\psi\). The mobilized dilation angle is based on Rowe’s dilation law (1962):
where the parameter \(F_c\) 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 \(\phi_{cv}\) is defined as
The mobilized friction \(\phi_m\) is defined in terms of the current stress state
A nonassociated flow rule consistent with the MohrCoulomb yield criterion is used in the model. The shear potential function is defined as
where \(m_1 = (1+\sin \psi_m)/2\) and \(m_3 = (1+\sin \psi_m)/2\).
In order to avoid overdilatancy when the soil reaches its critical void state at \(e_{max}\), the dilation angle needs a minor modification. One way proposed by Schanz et al. (1999) is to set a cutoff rule, so that
Here we introduce a smoothing technique to avoid a sudden change of dilation angle:
The dilation rules with cutoff, smoothing technique, and without dilation cutoff are schematically compared in Figure 2.
Volumetric Cap Criterion and Flow Rule
The volumetric (cap) yield function is defined as
where \(\alpha\) is a constant derived internally from other material parameters based on a virtual (numerical) oedometer test, \(\tilde{q}\) is a shear stress measure defined as \(\tilde{q} = − [\sigma_1 + (\delta  1)\sigma_2 − \delta \sigma_3]\), and \(\delta = (3 + \sin \phi)/(3 − \sin \phi)\). The initial value of hardening parameter \(p_c\), which denotes the preconsolidation pressure, can be determined using the initial stress state and overconsolidation ratio, \(OCR\), so that
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 \(\gamma^v\) (one of the internal variables) is defined incrementally as
where \(\Delta \varepsilon^p_v\) is the volumetric plastic strain increment.
Evolution of the hardening parameter, \(p_c\), is given by the relation:
where \(H_c\) is a constant parameter that can be derived internally from other material parameters. Instead of taking \(\alpha\) and \(H_c\) as input material parameters, another two parameters, \(K_{nc}\) and \(E^{ref}_{oed}\) are specified as input. Parameter \(K_{nc}\) denotes normal consolidation coefficient and \(E^{ref}_{oed}\) stands for the tangent oedometer stiffness at the reference pressure \(p^{ref}\). If \(K_{nc}\) is not provided by the user, it is taken as \(K_{nc} = 1 − \sin \phi\) (default value). The PH model reserves the flexibility for the user to input userdefined \(\alpha\) and \(H_c\) instead of internal calculation for \(\alpha\) and \(H_c\). A typical value for \(\alpha\) is between 0.6 and 2.0 for most soils. Parameter \(H_c\) is calculated as \(H_c = k K_p\), where \(k\) is a correction factor with a typical value of \(0 < k \le 1\), \(K_p = K_1 K_2 / (K_1  K_2)\), \(K_1 = E^{ref}_{ur} / [3(12\nu)]\), and \(K_2 = E^{ref}_{oed} (1 + 2 K_{nc})/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
where \(\sigma^t\) is the tension limit. By default, \(\sigma^t\) is zero and the user can provide a value with an upper limit of \(c/\tan \phi\). The model does not consider tension hardening.
Plastic Corrections
In the implementation of the PH model, elastic trial stresses, \(\sigma^I_{ij}\), are first computed by adding to the old stresses, \(\sigma^o_{ij}\), elastic increments computed using the total strain increment \(\Delta \varepsilon_{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:
where \(\alpha_1 = K +4G/3\), \(\alpha_2 = K −2G/3\), \(K\) and \(G\) are the tangent elastic bulk and shear modulus, respectively, and (\(\Delta \varepsilon_1, \Delta \varepsilon_2, \Delta \varepsilon_3\)) is the set of principal strain increments.
The trial internal variables, \(\gamma^p\) and \(\gamma^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, \(f^s > 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
The shear hardening parameter increment is obtained from Equation (13),
Thus, the corrected (new) stress components become:
and the corrected internal variable is
Using Equation (30) and definitions \(m_1\) and \(m_2\) from Equation (19), updated shear stress and parameter \(q_a\) (see Equation (9)) are calculated as
After substituting the corrected parameters \(q\), \(q_a\), \(\gamma^p\) into Equation (8), parameter \(\lambda^s\) can be obtained by solving a nonlinear equation \(f^s(\lambda^s)\) = 0 using an iteration method. Note that \(\gamma^p\) in Equation (31) uses the current value of \(\lambda^s\) updated in each iteration step, which implies that the correction algorithm is semiimplicit.
If the trial shear stress exceeds the MohrCoulomb shear failure surface, the corrections used in MohrCoulomb model implementation apply.
Now consider the case when the elastic guess (Equation (22)) exceeds the volumetric yield criterion, \(f^v > 0\), and the volumetric hardening occurs. The plastic strain increments are related to the gradient of the potential function in the stress space as
Evolution of the hardening parameter \(\gamma^v\) is given by the relation:
In the above formulation, the plastic strain and volumetric hardening parameter increment are related to the current stress measurement (\(p\) or \(\tilde{q}\)), so the correction algorithm again uses a semiimplicit approach. The corrected stress components are obtained by analogy with shear hardening,
Noting that \(p^{I} = −(\sigma_1 + \sigma_2 + \sigma_3)/3\) and \(\tilde{q}^{I} = − [\sigma_1 + (\delta  1)\sigma_2 − \delta \sigma_3]\), and using Equation (36) we get
where \(M = 4G(\delta^2 − \delta + 1) / \alpha^2\). After substituting the corrected \(p\) and \(\tilde{q}\) into the yield function \(f^v\), (Equation (22)), parameter \(\lambda^v\) can be obtained by solving the nonlinear equation \(f^v(\lambda^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:
Parameters \(\gamma^s\) can \(\gamma^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 MohrCoulomb failure surface, the MohrCoulomb failure function should be used instead of the shear hardening function.
Implementation Procedure
In the initialization step, the initial stress, evolution parameters \(\gamma^s\) can \(\gamma^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 subzone 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 tensioncorrected 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 subzones 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 zoneaveraged stresses.
Details on the implementation can be found in Cheng & Detournay (2016).
Comparison of the PH model and CYSoil model
The PH and CYSoil models share many similar features. For example, both models have stressdependent 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 MohrCoulomb 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 (\(\gamma^p\)), while the definition for \(\gamma^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 \(E_{50}\) 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 \(\alpha\) is userdefined, 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 unloadingreloading 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/cm^{2}. The initial void ratios are 0.783, 0.786, and 0.781, respectively.
break
Calibration of Friction Angle, \(\phi\), and Cohesion, \(c\):
The calibration of these two parameter has no difference from the MohrCoulomb model.
To estimate values of friction angle \(\phi\) and cohesion \(c\) from triaxial test data using the relation \(q = \sigma_3  \sigma_1\) vs. \(p = (\sigma_1 + 2\sigma_3)/3\), follow the procedure outlined below.
Plot deviatoric stress \(q\) vs. \(p\) using triaxial compression test lab data.
Use a trend line to fit the MohrCoulomb envelope.
The slope of the trend line is \(6 \sin \phi / (3 − \sin \phi)\), which follows from an alternative form of the MohrCoulomb failure criterion. This determines friction angle \(\phi\).
The intercept of the trend line is \(6c \cos \phi / (3 − \sin \phi)\), which determines cohesion \(c\), as the friction angle \(\phi\) is already known from the previous step.
Alternatively, these two parameters can be estimated from triaxial test data using the relation \(\hat{q} = (\sigma_3  \sigma_1)/2\) vs. \(\hat{p} = (\sigma_1 + \sigma_3)/2\) by following the procedure outlined below.
Plot deviatoric stress \(\hat{q}\) vs. \(\hat{p}\) using triaxial compression test lab data.
Use a trend line to fit the MohrCoulomb envelope.
The slope of the trend line is \(\sin \phi\), which follows from an alternative form of the MohrCoulomb failure criterion. This determines friction angle \(\phi\).
The intercept of the trend line is \(c \cos \phi\), which determines cohesion \(c\), as the friction angle \(\phi\) is already known from the previous step.
For the Monterey sand, the slope of the trend line of the MohrCoulomb 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.
Calibration of \(R_f\), \(m\), and \(E^{ref}_{50}\):
To estimate values of \(R_f\), \(m\), and \(E^{ref}_{50}\) from triaxial test data, follow the procedure outlined below.
Plot the curve of \(\varepsilon_1/q\) vs. \(\varepsilon_1/q_f\) using the triaxial compression test data.
Use a trend line to fit the data. The slope of the line is \(R_f\) , the intercept is \(1/E_i\), based on the relation \(\varepsilon_1/q = R_f (\varepsilon_1/q_f) + 1/E_i\) (derived from Equation (12) and \(R_f = q_f / q_a\)).
Multiple sets of triaxial compression tests with various confining pressures can be used to produce a set of pairs of (\(R_f\), \(E_{50}\)). The final \(R_f\) can be estimated as the averaged one, and the pairs of (\(E_i\), \(\sigma_3\)) determine \(m\) and \(E^{ref}_{50}\) based on Equation (7).
Figure 6 plots the curves of \(\varepsilon_1/q\) vs. \(\varepsilon_1/q_a\) using triaxial compression test data of the Monterey sand with confining pressures 1.2, 0.6, and 0.3 kgf/cm^{2}. The slopes of these lines are \(R_f\), and the intercepts are \(100/E_i\) (as strain is given in percentage).
This figure determines three pairs of (\(R_f\), \(E_{50}\)), summarized in Table 1. The average value for \(R_f\) is 0.957. Parameter \(p^{ref}\) needs no calibration and its value is assumed to be 0.1 kgf/cm^{2}. Finally, plotting parameters \(\ln (−\sigma_3/p^{ref})\) vs. \(\ln (E50)\), as shown in Figure 7, and using Equation (7) allows for determination of \(m\) and \(E^{ref}_{50}\) (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 \(E^{ref}_{50}\) = \(\exp{}\) (4.63) = 102.5 kgf/cm^{2}.
\(\sigma_3\) 
\(R_f\) 
\(1/E_i\) 
\(E_i\) 
\(E_{50}\) 
\(\ln (E_{50})\) 
\(\ln \left( {c \cot \phi  \sigma_3} \over {c \cot \phi + p^{ref}} \right)\) 
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
break
Calibration of \(E^{ref}_{ur}\):
After the calibration of parameter \(m\), it is straightforward to obtain \(E^{ref}_{ur}\) from Equation (5) using the unloadingreloading modulus \(E_{ur}\) obtained from the original \(q\) vs. \(\varepsilon_1\) data in a triaxial compression test. The final value of \(E^{ref}_{ur}\) can be averaged using data for different confining pressures. For the example of Monterey sand, \(E^{ref}_{ur}\) is determined to be 320.0 kgf/cm^{2}. If the unloadingreloading moduli are not available, a value in the range of \((3 \sim 5) E^{ref}_{50}\) can be used for most soils. The PH model uses a default value of \(4 E^{ref}_{50}\) if no input is provided for \(E^{ref}_{ur}\).
Calibration of Dilation Angle, \(\psi\):
The ultimate dilation angle \(\psi\) can be estimated as \(\sin \psi = s_m / (s_m+2)\), where \(s_m\) is the maximum slope (\(s_m\)) of the \(\varepsilon_v\) vs. \(\varepsilon_1\) curve in the triaxial compression tests. The slope \(s_m\) is found as: \(s_m = \Delta \varepsilon_v/\Delta \varepsilon_1 \approx \Delta \varepsilon^p_v/\Delta \varepsilon^p_1 = 2\sin \psi/(1  \sin \psi)\). 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 \(K_{nc}\) and \(E^{ref}_{oed}\):
Parameters \(K_{nc}\) and \(E^{ref}_{oed}\) can be calibrated from the oedometer tests. The ultimate value of \(\sigma_3 / \sigma_1\) from an oedometer test is \(K_{nc}\). When \(\sigma_1 = p^{ref}\), the tangent modulus of the \(\sigma_1\) vs. \(\varepsilon_1\) curve is \(E^{ref}_{oed}\).
Calibration of Poisson’s Ratio, \(\nu\):
Poisson’s ratio can be estimated from the initial slope \(s_0\) (usually a negative value) of the \(\varepsilon_v\) vs. \(\varepsilon_1\) curve, so that \(\nu = (1+\tan s_0)/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 \(e_{max}\):
The material parameter \(e_{max}\) is needed if the dilation smoothing technique (Figure 2) is required. The value of \(e_{max}\) 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 trialanderror 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 \(e_{max}\) = 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 \(\sigma^{ini}_1\), \(\sigma^{ini}_2\), \(\sigma^{ini}_3\), \(e_{ini}\), and \(OCR\) are known initial parameters and should be consistent with the initial conditions. Some parameters can be estimated through the trialanderror method.
A summary of all material properties determined for the Monterey sand is provided in Table 2.
Parameters 
\(\sigma_3\) 
\(\sigma_3\) 
\(\sigma_3\) 
\(E^{ref}_{50}\) (kgf/cm^{2}) 
102.5 

\(E^{ref}_{ur}\) (kgf/cm^{2}) 
320.0 

\(m\) 
0.707 

\(R_f\) 
0.957 

\(p^{ref}\) (kgf/cm^{2}) 
0.1 

\(\nu\) 
0.3 

\(\phi\) (degrees) 
34.65 

\(c\) (kgf/cm^{2}) 
0.0 

\(e_{max}\) 
0.803 

\(K_{nc}\) 
0.5 

\(\psi\) (degrees) 
6.1 
6.4 
7.0 
\(\sigma^{ini}_{1 \sim 3}\) (kgf/cm^{2}) 
1.2 
0.6 
0.3 
\(OCR\) 
1.0 
2.0 
4.0 
\(e_{ini}\) 
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
SmallStrain Stiffness Option
Using the PlasticHardening model to simulate the unloadingreloading paths, the reloading path will be the reverse of the unloading path. However, this elastic behavior assumes strains that are in the very small strain range. Experiments illustrate that soil stiffness modulus is straindependent — it decreases nonlinearly with increase in shear strain. The strain can be categorized in three levels: very small strain, where the stiffness modulus is assumed constant (\(G_0\) or sometimes called \(G_{max}\)); small strain, where the stiffness modulus reduces nonlinearly as the strain increases; and large strain, where the soil is in or close to failure and the stiffness modulus is relatively small. Figure 10 illustrates, showing the normalized stiffness modulus degradation curve with reference to geotechnical applications and measurement accuracy from laboratory or in situ tests. It is important that the soil stiffness modulus selected for simulation should not be associated with strain levels at the end of geotechnical construction, or from laboratory tests with well disturbed samples, because the soil stiffness modulus may have decreased much from its initial value. The additional option of the smallstrain stiffness modulus for the PlasticHardening model is expected to take straindependency into account.
To activate smallstrain stiffness in the PlasticHardening model, flagsmallstrain (the default is off) should be turned on. Only two extra properties require assignment: stiffness0reference (\(E_0^{ref}\)) and strain70 (\(\gamma_{70}\)).
The initial stiffness modulus at the reference pressure, \(E_0^{ref}\), will determine the shear stiffness modulus (\(G_0\)) at the very small strain level. It can be estimated by the relations of
and
where \(Z\), \(m\), and \(\nu\) are defined in the standard PlasticHardening model. The default value of \(E_0^{ref}\) is taken as \(3 E_{ur}^{ref}\) if no specific data is available. \(E_0^{ref}\) is not allowed to be equal or smaller than \(E_{ur}^{ref}\) in this model. Together with the above two equations, \(E_0^{ref}\) could be estimated from
where \(\rho\) is the soil density, and \(V_s\) is the measured in situ soil shear wave velocity. \(G_0\) can be estimated from other insitu test data, e.g., SPT, CPT & DMT, or from some accepted empirical relations.
The reference shear strain (\(\gamma_{70}\)) is selected as the value at which the secant shear stiffness modulus is about 70% (to be exact it should be 72.2%) of \(G_0\). The default value for \(\gamma_{70}\) is 2e4. Typical values of \(\gamma_{70}\) are listed in Table 3:
Sand: 

At depth 0 to 6 m 
1.1e4 
At depth 6 to 15 m 
1.6e4 
At depth 15 to 37 m 
2.3e4 
At depth 37 to 76 m 
3.3e4 
At depth 76 to 150 m 
4.5e4 
At depth 150 to 300 m 
7.0e4 
Clay: 

Plastic Index = 30 
3.5e4 
Plastic Index = 50 
7.7e4 
Plastic Index = 100 
1.8e3 
Once \(G_0\) and \(\gamma_{70}\) are known, the smallstrain secant stiffness modulus is calculated as
and its tangent stiffness modulus is
The tangent shear stiffness modulus sets a lowerbound value of \(G_{ur} = E_{ur} / [2(1+\nu)]\) in this model.
The ‘brick’ algorithm (Cudny and Truty, 2020) to consider the unloadingreloading straindependent stiffness modulus is implemented in this model. The previous shortcoming based on the algorithm by Benz (2007) of not good performance on nested hysteretic loops (Niemunis and Cudny, 2018) is removed.
References
Atkinson, J. H., & Sallfors, G. “Experimental determination of stress/straintime characteristics in laboratory and in situ tests,” in Proceedings, 10th ECSMFE, Florence, Rotterdam: Balkema, 959999 (1991).
Benz, T. Smallstrain stiffness of soils and its numerical consequences (Vol. 5). Stuttgart: Univ. Stuttgart, Inst. f. Geotechnik (2007).
Cudny, M., and Truty, A. Refinement of the Hardening Soil model within the small strain range. Acta Geotechnica, 15(8), 20312051 (2020).
Cheng, Z., and C. Detournay. Plastic hardening model I: Implementation in FLAC3D. In Proc. 4th Itasca Symposium on Applied Numerical Modeling. 2016.
Cheng, Z., and A. Lucarelli. Plastic hardening model II: Calibration and validation. In Proc. 4th Itasca Symposium on Applied Numerical Modeling. 2016.
Duncan, J. M., and C. Y. Chang. “Nonlinear Analysis of Stress and Strain in Soils,” Soil Mechanics, 96 (SM5), 16291653 (1970).
Niemunis, A., and M. Cudny. Discussion on “Dynamic soilstructure interaction: A threedimensional numerical approach and its application to the Lotung case study”. Poor performance of the HSS model. Computer and Geotechnics, 98(6), 243245 (2018).
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,” in Beyond 2000 in computational geotechnics, 281296 (1999).
plastichardening Model Properties
Use the following keywords with the zone property
command to set these properties of the PlasticHardening model.
 cohesion f
cohesion, \(c\). The default and cutoff (allowable minimum) value is \(10^{5} \cdot p^{ref}\).
 exponent f
exponent for elastic moduli, \(m\). Allowable value is \(0 \le m \le 0.999\).
 friction f
ultimate friction angle in degrees, \(\phi\). The default and cutoff (allowable minimum) value is 0.001 degrees.
 stress1effective f
initial minimum principal effective stress, \(\sigma^{ini}_1\). See the note*.
 stress2effective f
initial median principal effective stress, \(\sigma^{ini}_2\). See the note*.
 stress3effective f
initial maximum principal effective stress, \(\sigma^{ini}_3\). See the note*.
 pressurereference f
reference pressure, \(p^{ref}\) . A nonzero value should be specified based on the unit of stress/pressure adopted in the model. The most commonly used and highly recommended reference pressure is the standard atmospheric pressure, e.g., 100.0 kPa/1.0e5 if in Pa (or other value compatible to the adopted stress/pressure unit).
 stiffness50reference f
secant stiffness, \(E^{ref}_{50}\), at 50% of the ultimate deviatoric stress, \(q_f\), when \(\sigma_3 =  p^{ref}\) in a triaxial test. A nonzero value must be specified.
 flagbrittle b (a)
If true, the tension limit is set to 0 in the event of tensile failure. The default is false.
 stiffness0reference f (a)
Initial stiffness, \(E^{ref}_{0}\), at zero strain when \(\sigma_3 =  p^{ref}\). It is required that \(E^{ref}_{0} > E^{ref}_{ur}\). If not specified, the default value is \(E^{ref}_{0} = 3 \cdot E^{ref}_{ur}\). This property is needed only when the flag for smallstrain stiffness is on, otherwise its value is useless.
 stiffnessoedometerreference f (a)
tangent stiffness, \(E^{ref}_{oed}\), when \(\sigma_1 =  p^{ref}\) in an oedometer test. By default, \(E^{ref}_{oed} = E^{ref}_{50}\) if \(E^{ref}_{oed}\) is not specified while \(E^{ref}_{50}\) has been specified. If \(E^{ref}_{oed} < 0.5 E^{ref}_{50}\), the soil is highly compressible, it is recommended to use the SoftSoil (SS) model instead.
 stiffnessurreference f (a)
unloadingreloading stiffness, \(E^{ref}_{ur}\), when \(\sigma_3 =  p^{ref}\) in a triaxial test. This value must be greater than \(2 \cdot E^{ref}_{50}\). By default, \(E^{ref}_{ur} = 4 \cdot E^{ref}_{50}\).
 coefficientnormallyconsolidation f (a)
normal consolidation coefficient, \(K_{nc}\). It is not allowed to be less than \(\nu / (1\nu)\), a range 0.5 to 0.7 is common. The default is \(K_{nc} = 1  \sin \phi\).
 constantalpha f (a)
dimensionless parameter in the cap yield function, \(\alpha\). Initialized internally. See the note* for an alternative option.
 factordilationlaw f (a)
factor of the dilation law, \(F_c\), allowable range is 0.0 to 0.25. The default is 0.0.
 flaginitialization s (a)
flag for initialization. If set to on, the internal variables will be recalculated. The default is off.
 flagsmallstrain s (a)
flag for smallstrain stiffness. If set on, smallstrain stiffness is activated. The default is off.
 overconsolidationratio f (a)
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” noncap shear plastic hardening model.
 stiffnesscaphardening f (a)
hardening modulus for cap hardening, \(H_c\). Initialized internally. See the note* for an alternative option.
 strain70 f (a)
Shear strain at which \(G_s / G_0\) = 72.2%, If not specified, the default value is \(\gamma_{70} = 2 \times 10^{4}\). This property is needed only when the flag for smallstrain stiffness is on, otherwise its value is useless.
 voidinitial f (a)
initial void ratio, \(e_{ini}\). 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.
 voidmaximum f (a)
allowable maximum void ratio, \(e_{max}\) . The default is 999.0, which is not realistic but it assures that dilation smoothing effectively will not be activated.
 bulk f (r)
current loadingunloading bulk modulus, \(K\)
 shear f (r)
current unloadingreloading shear modulus, \(G\)
 plastichardeningshear f (r)
accumulated shear plastic hardening parameter, \(\gamma^p\)
 plastichardeningvolume f (r)
accumulated volumetric plastic hardening parameter, \(\varepsilon^p_v\)
 pressurecap f (r)
cap (preconsolidation) pressure, \(p_c\). Determined by the current stress, \(\alpha\), and \(OCR\).
 void f (r)
current void ratio, \(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) Readonly property.
This property cannot be set by the user. Instead, it can be listed, plotted, or accessed through FISH.
Notes
The tension cutoff is \({\sigma}^t = min({\sigma}^t, c/\tan \phi)\).
The parameters \(\sigma^{ini}_1\), \(\sigma^{ini}_2\), \(\sigma^{ini}_3\) must be specified via commands or a FISH function the first time the PH model is assigned because it is a stressdependent 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 \(\alpha\) and \(H_c\) are specified nonzero values by the user, the PH model will use these userdefined parameters \(\alpha\) and \(H_c\), and the internal calculation for \(\alpha\) and \(H_c\) will not be performed.
Was this helpful? ...  Itasca Software © 2024, Itasca  Updated: May 07, 2024 