====== EP-PLASOL ======
===== Description =====
Elasto-plastic constitutive law for solid elements at constant temperature (non-associated) with linear elasticity. Isotropic hardening/softening of friction angle and cohesion is possible. \\
\\
Integration is performed using an implicit backward Euler scheme with a return mapping normal to the flow surface g. This law can take into account the influence of:
  * the first stress invariant, i.e. the yield surface is either parallel to the pressure axis in the p-q plane (see <imgref image1>-A)
  * the third stress invariant, i.e. the Lode angle : the trace in the deviatoric plane is either a circle or a smoothed irregular hexagon (see <imgref image1>-B).
<imgcaption image1|A) Influence of the first stress invariant in the p-q plane; B) Influence of the third stress invariant in the deviatoric plane>{{  :laws:p_q_plane_and_deviatoric_plane.png?600  |}}</imgcaption>

The von Mises, the Drücker Prager and a smoothed Mohr Coulomb yield surfaces can be represented. 
==== The model ====
This law is only  used for mechanical analysis of elasto-plastic isotropic porous media undergoing large strains.
=== Yield and flow surfaces ===
Stresses and stress invariants \\
\[ I_{\sigma} = \sigma_{ij}\delta_{ij} = \sigma_{ii}; \widehat{\sigma}_{ij} = \sigma_{ij} - \frac{I_\sigma}{3}\delta_{ij};\] 
\[ II_{\widehat{\sigma}} = \sqrt{\frac{1}{2}\widehat{\sigma}_{ij}\widehat{\sigma}_{ij}}; III_{\widehat{\sigma}} = \frac{1}{3}\widehat{\sigma}_{ij}\widehat{\sigma}_{jk}\widehat{\sigma}_{ki} ;\] 
\[\beta = -\frac{1}{3}\sin^{-1}\left( \frac{3\sqrt{3}}{2} \frac{III_{\widehat{\sigma}}}{II_{\widehat{\sigma}}^3} \right) \]

__Criterion with friction angle different from 0 (Drücker Prager or Van Eekelen):__ \\
Regular criterion used if $I_{sigma} - m’ II_{\widehat{\sigma}} < (3c/\tan \phi_c)$
\[ f = II_{\widehat{\sigma}} + m\left( I_{\sigma} - \frac{3c}{\tan\phi_c} \right) = 0\]
with 
|Drücker Prager| $m = \frac{2 \sin \phi_c}{\sqrt{3}(3 - \sin \phi_c)}$ |
|Van Eekelen| $m = a(1+b\sin 3\beta)^n$|
where $a$ and $b$ are function of $\phi_C$, $\phi_E$ and $n$. \\
Apex criterion used if $I_{\sigma} - m’ II_{\widehat{\sigma}} \geq (3c / \tan \phi_c)$
\[ f = I_{\sigma} - \frac{3c}{\tan\phi_c} = 0 \]
$m’$ is the equivalent of m but for the flow surface (i.e. $\phi$ is replaced by $\psi$)

__Criterion with friction angle egal tp 0 (Von Mises criterion):__
\[ f = II_{\widehat{\sigma}} - \frac{2c}{\sqrt{3}} = 0 \]

=== Hardening/softening ===
Hardening/softening is assumed to be represented by the evolution of friction angles and/or cohesion as a function of the Von Mises equivalent plastic strain 
\[ \varepsilon_{ep}^p = \sqrt{\frac{2}{3}\widehat{\varepsilon}_{ij}^p \widehat{\varepsilon}_{ij}^p }\]
Hyperbolic functions are used:
  - If ILODE = 1 or 2, \\
    * if $\varepsilon_{eq}^p <$ decphi : $\phi_C = \phi_{C0}$
    * if $\varepsilon_{eq}^p >$ decphi : $\phi_C = \phi_{C0}+\frac{(\phi_{Cf} - \phi_{C0})(\varepsilon_{eq}^p - decphi)}{B_p + (\varepsilon_{eq}^p -decphi)}$
    * if $\varepsilon_{eq}^p <$ decphi : $c = c_0$
    * if $\varepsilon_{eq}^p >$ decphi : $c = c_0 + \frac{(c_f - c_0)(\varepsilon_{eq}^p - deccoh)}{B_c + (\varepsilon_{eq}^p -deccoh)} $
  - ONLY if ILODE = 2, 
    * if $\varepsilon_{eq}^p <$ decphi : $\phi_E = \phi_{E0}$
    * if $\varepsilon_{eq}^p >$ decphi : $\phi_E = \phi_{E0}+\frac{(\phi_{Ef} - \phi_{E0})(\varepsilon_{eq}^p - decphi)}{B_p + (\varepsilon_{eq}^p -decphi)}$
Where coefficients $B_p$ and $B_c$ are respectively the values of equivalent plastic strain for which half of the hardening/softening on friction angles and cohesion is achieved (see <imgref image2>). 
<imgcaption image2|Hardening softening hyperbolic relation for 2 values of coefficient Bp>
{{  :laws:hardening_softening_hyperbolic_relation.png?600  |}}</imgcaption>

=== Viscoplasticity ===
When viscosity is taken into account (IVISCO=1), one can assume that the plastic strain is composed of a time independent instantaneous strain $\varepsilon_{ij}^p$, but also of a time-dependent creep strain $\varepsilon_{ij}^{vp}$. The total strain (elastic, plastic and viscoplastic) reads :
\[\varepsilon_{ij} = \varepsilon_{ij}^e +\varepsilon_{ij}^p + \varepsilon_{ij}^{vp}\]
The material viscosity implies a time-dependent strain $\varepsilon_{ij}^{vp}$ which is a delayed plastic deformation. Under soil mechanic convention (compressive stress is taken as positive), the loading surface of the viscoplastic flow $f_{vp}$ and the viscoplastic potential $Q_{vp}$ are controlled by a delayed viscoplastic hardening function $\alpha_{vp}$ and read :
\[ f_{vp} = q - \alpha_{vp} g(\theta) \bar{R}\sqrt{A\left(C_s + \frac{p'}{\bar{R}}\right)} \geq 0 \]
\[ Q_{vp} = q - g(\theta)\left( \alpha_{vp} - \beta_p \right)\left( p' + C_s \bar{R} \right)\]
\[ \alpha_{vp} = \alpha_{vp,0} + \left( 1 - \alpha_{vp,0} \right) \frac{\gamma_{vp}}{B_{vp} + \gamma_{vp}}\]
where: 
  * $q=\sqrt{3}II_{\widehat{\sigma}}$ is the deviatoric stress 
  * $p' = \frac{\sigma_{ii}'}{3}$ is the mean effective stress 
  * $\gamma_{vp}$ is the equivalent plastic shear strain, i.e. the generalized plastic distorsion 
  * $g(\theta)$ is a function allowing to take into account the influence of the Lode angle $\theta$, i.e. the third deviatoric stress invariant, on yield surface in the deviatoric plane, $g(\theta) = 1$ for the sake of simplicity 
  * $\bar{R}$ is a normalising parameter, for convenience it is taken as equal to the uniaxial compression strength $\bar{R} = R_c$ 
  * $\alpha_{vp,0}$ is the initial threshold for the viscoplastic flow (a value of $\alpha_{vp,0} = 0$ could be taken in agreement with the instantaneous plastic mechanism) 
  * $B_{vp}$ is a parameter controlling the evolution of $\alpha_{vp}$ and therefore of $f_{vp}$ 
  * $A$ is an internal friction coefficient defining the curvature of the failure surface and $C_s$ is a cohesion coefficient, it denotes the material cohesion in saturated condition. Finally, $\beta_p$ is a parameter which defines the transition from compressibility $(\alpha_{vp} < \beta_p)$ to dilatancy $(\alpha_{vp} > \beta_p)$. $A$, $C_s$ and $\beta_p$ are parameters linked to the instantaneous plastic deformation modelling. \\

The equivalent plastic shear strain increment is given by : 
\[ \dot{\gamma_{vp}} = \sqrt{\frac{2}{3}\dot{e}_{ij}^{vp}\dot{e}_{ij}^{vp}} \]
where $\dot{e}_{ij}^{vp}$ is the viscoplastic deviatoric strain :
\[ \dot{e}_{ij}^{vp} = \dot{\varepsilon}_{ij}^{vp} - \dot{\varepsilon}_{kk}^{vp}\delta_{ij} \]
where $\delta_{ij}$ is the Kronecker symbol. The viscoplastic flow rule is determined as follows : 
\[ \dot{\varepsilon}_{ij}^{vp} = A(T) \langle \frac{f_{vp}}{\bar{R}} \rangle^n \frac{\partial Q_{vp}}{\partial \sigma_{ij}} \]
where $\langle\rangle$ is the Macauley bracket, $\langle x\rangle = x$ if $x\geq 0$ and $\langle x\rangle = 0$ if $x<0$, $A(T)$ is thefluidity coefficient generally dependent of the temperature $T$ and $n$ is a parameter which describes the shape of the creep curve. The following function is used for the fluidity : 
\[ A(T) = A_0 \exp\left( -\frac{\zeta}{RT} \right) \]
where $A_0$ is the fluidity value at a reference temperature, $R$ is the perfect gas universal constant, $T$ is the absolute temperature and $\zeta$ is a parameter controlling the influence of temperature on the material viscosity. However, the temperature is generally assumed constant. \\
\\
The viscoplastic law definition and typical values of parameters for sandstone and hard clay can be found in Zhou et al. (2008)((H. Zhou, Y. Jia, and J.F. Shao. A unified elastic-plastic and viscoplastic damage model for quasi-brittle rocks. //International Journal of Rock Mechanics and Mining Sciences//, 45:1237–1251, 2008.)), Jia et al. (2008)((Y. Jia, H.B. Bian, G. Duveau, K. Su, and J.F. Shao. Hydromechanical modelling of shaft excavation in meuse/haute-marne laboratory. //Physics and Chemistry of the Earth//, 33:422–435, 2008.)). To modify parameters go to the file “LMLVP.F” in LAGAMINE code. 

==== Files ====
Prepro: LPLA.F \\
Lagamine: PLA2EA.F, PLA3D.F
===== Availability =====
|Plane stress state| NO |
|Plane strain state| YES |
|Axisymmetric state| YES |
|3D state| YES |
|Generalized plane state| NO |
===== Input file =====
==== Parameters defining the type of constitutive law ====
^ Line 1 (2I5, 60A1)^^
|IL|Law number|
|ITYPE| 72 |
|COMMENT| Any comment (up to 60 characters) that will be reproduced on the output listing|
==== Integer parameters ====
^ Line 1 (9I5) ^^
|NINTV| = number of sub-steps used to integrate numerically the constitutive equation in a time step. |
|:::| If NINTV = 0 : number of sub-steps is based on the norm of the deformation increment and on DIV|
|ISOL| = 0 : use of total stresses in the constitutive law|
|:::| $\neq$ 0 : use of effective stresses in the constitutive law. See [[appendices:a8|Appendix 8]] |
|ICBIF| = 0 : nothing|
|:::| = 1 : Rice bifurcation criterion is computed (only for 2D plane strain analysis) |
|ILODEF| Shape of the yield surface in the deviatoric plane : |
|:::| = 1 : circle in the deviatoric plane|
|:::| = 2 : smoothed irregular hexagon in the deviatoric plane|
|ILODEG| Shape of the flow surface in the deviatoric plane : |
|:::| = 1 : circle in the deviatoric plane|
|:::| = 2 : smoothed irregular hexagon in the deviatoric plane|
|IECPS| = 0 : $\psi$ is defined with PSIC and PSIE|
|:::| = 1 : $\psi$ is defined with PHMPS|
|:::| = 2 : Variable dilatancy (El Moustapha,2014) ((El Moustapha, K. (2014) ‘Identification of an enriched constitutive law for geomaterials in the presence of a strain localisation’, Thesis, Liège University.))|
|:::| = 3 : Variable dilatancy (Salehnia, 2015)((Salehnia, F. (2015) From some obscurity to clarity in Boom clay behavior: Analysis of its coupled hydro-mechanical response in the presence of strain localization. Thesis, Liège University.))|
|KMETH| = 2 : actualised VGRAD integration|
|:::| = 3 : Mean VGRAD integration (Default value) |
|IREDUC| = 0 : nothing |
|:::| 1 : Phi-C reduction method|
|ICOCA| = 0 : nothing |
|:::| 1 : Capillary cohesion formulation ($c = c_0 + AK1.s + AK2.s^2$) |
|:::| 2 : Capillary cohesion formulation ($c = c0 + AK1.log(s) + AK2$) |
|:::| 3 : Capillary cohesion formulation ($c = c_0.(1+AK1.s)$) and capillary Young’s modulus formulation ($E = E0.(1+AK2.S)$). |
|:::| Available for cohesion hardening/softening (IECROU=2). |
|IVISCO| = 0 : nothing |
|:::| 1 to 3 : viscoplastic model |
|IDAM| = 0 : nothing |
|:::| 1 : damage of elastic properties |
|:::| 2 : concrete hydration via .hydr file |
==== Real parameters ====
^ Line 1 (8G10.0) ^^
|E| YOUNG’s elastic modulus |
|ANU| Poisson ratio |
|PSIC| Dilatancy angle (in degrees) for compressive paths |
|PSIE| Dilatancy angle (in degrees) for extensive paths (ssi ILODEG=2) |
|RHO| Specific mass |
|DIV| Size of sub-steps for computation of NINTV (only if NINTV=0; Default value=5.D-3) |
|PHMPS| Constant value for definition of |
|AE| If AE$\neq$0, Linear evolution of young modulus with confinement pressure is considered $E=E0+AE \cdot \sigma$|
^ Line 2 (7G10.0) ^^
|PHIC0| Initial Coulomb's angle (in degrees) for compressive paths  |
|PHICF| Final Coulomb's angle (in degrees) for compressive paths  |
|BPHI| Only if there is hardening/softening |
|PHIE0| Initial Coulomb’s angle (in degrees) for extensive paths  |
|PHIEF| Final Coulomb’s angle (in degrees) for extensive paths (if and only if ILODEF = 2) |
|AN| Van Eekelen exponent (default value=-0.229) |
|DECPHI| Coulomb’s angle hardening shifting |
^ Line 3 (7G10.0) ^^
|COH0| Initial value of cohesion |
|COHF| Final value of cohesion |
|BCOH| Only if there is hardening/softening |
|BIOPT| Parameter for optimising the bifurcation computation based on Linear Comparison Solid (L.C.S). If BIOPT = 0 then Upper Bound L.C.S. If BIOPT = 1 then Lower Bound of L.C.S. If 0<BIOPT<1 then Intermediate L.C.S.|
|AK1|Capillary cohesion first parameter |
|AK2|Capillary cohesion second parameter |
|DECCOH| Cohesion hardening shifting |
^ Line 4 - Only If IECPS = 2 Or IECPS = 3 (7G10.0) ^^
|PSICPEAK| Peak of dilatancy angle for compressive paths (If IECPS=2 then PSICPEAK is the initial value of dilatancy angle|
|PSICLIM| Limit value of dilatancy angle for compressive paths|
|RATPSI| Ratio between initial and peak of dilatancy angle|
|BPSI| Value of EEQU for which PSIC=0.5 (PSICPEAK - PSICLIM)  |
|PSIEPEAK| Peak of dilatancy angle for extensive paths (If IECPS=2 then PSIEPEAK is the initial value of dilatancy angle)  |
|PSIELIM| Limit value of dilatancy angle for extensive paths|
|DECPSI| Value of EEQU when the dilatancy angle has been half decreased between its initial and final values|
^ Line 4 - Only if IDAM = 1 AND (IECPS = 0 or IECPS =1) (2G10.0) \\ Or Line 5 - Only If IDAM = 1 AND (IECPS = 2 or  IECPS = 3) (2G10.0) ^^
|P| Parameter controlling the damage evolution rate  |
|YD0| Initial threshold  |
^ Line 4 - Only if IDAM = 2 AND (IECPS = 0 or IECPS =1) (5G10.0) \\ Or Line 5 - Only If IDAM = 2 AND (IECPS = 2 or IECPS = 3)(5G10.0) ^^
|__Note__: The evolution of the hydratation degree must be specified in file *.hydr such that: \\ TIMES \\ (I10) number of time steps defined below \\ (2G10.0, repeated) Time step, alpha||
|FCF| Final simple compression resistance |
|ALPHATH| Hydration threshold (properties are equal to Eth/E0 |
|Pa| Parameter for the evolution of strength |
|Pb|Parameter for the evolution of stiffness|
|Eth|Elastic stiffness below threshold |
===== Stresses =====
==== Number of stresses ====
= 4 : for 2D analysis \\ 
= 6 : for 3D analysis
==== Meaning ====
The stresses are the components of CAUCHY stress tensor in global (X,Y,Z) coordinates. \\
For 2D analysis :
|SIG(1)|$\sigma_{xx}$|
|SIG(2)|$\sigma_{yy}$|
|SIG(3)|$\sigma_{xy}$|
|SIG(4)|$\sigma_{zz}$|
For the 3-D analysis :
|SIG(1)|$\sigma_{xx}$|
|SIG(2)|$\sigma_{yy}$|
|SIG(3)|$\sigma_{zz}$|
|SIG(4)|$\sigma_{xy}$|
|SIG(5)|$\sigma_{xz}$|
|SIG(6)|$\sigma_{yz}$|

===== State variables =====
==== Number of state variables ====
= 43 : for 2D plane strain analysis with bifurcation criterion (ICBIF=1)\\
= 31 : in all the other cases
==== List of state variables ====
|Q(1)| = 1 in plane strain state |
|:::| = circumferential strain rate ($\dot{\varepsilon_{\theta}}$) in axisymmetrical state|
|Q(2)| = actualised specific mass |
|Q(3)| = Reduced deviatoric stress (varies from 0 to 1) |
|Q(4)| = 0 if the current state is elastic |
|:::| = 1 if the current state is elasto-plastic|
|Q(5)| = equivalent viscoplastic shear strain, i.e. the generalized plastic distorsion, which increment is $\dot{\gamma_{vp}} = \sqrt{\frac{2}{3}\dot{e}_{ij}^{vp}\dot{e}_{ij}^{vp}}$ (see below for the viscoplastic law definition) |
|Q(6)| = equivalent strain $n^o$1 $\varepsilon_{eq1} = \int \Delta \dot{\varepsilon}_{eq}\Delta t$ |
|Q(7)| = equivalent strain indicator $n^o 1$ (Villote $n^o 1$) $\alpha_1 = (\Delta\dot{\varepsilon}_{eq}\Delta t ) / \varepsilon_{eq1}$ |
|Q(8)| = $\varepsilon_{xx}$ |
|Q(9)| = $\varepsilon_{yy}$ |
|Q(10)| = $\varepsilon_{zz}$ |
|Q(11)| = $\gamma_{xy} = 2.\varepsilon_{xy}$ |
|Q(12)| = equivalent strain $n^o 2$ $\varepsilon_{eq2} = \int \Delta \varepsilon_{eq}$ |
|Q(13)| = equivalent strain indicator $n^o 2$ (Villote $n^o 2$) $\alpha_2 = \Delta\varepsilon_{eq} / \varepsilon_{eq2}$|
|Q(14)| = actualised value of equivalent plastic strain $\varepsilon_{ep}^{p}$ |
|Q(15)| = actualised value of cohesion $c$ |
|Q(16)| = actualised value of Coulomb’s friction angle for compr. paths $\phi_C$|
|Q(17)| =  actualised value of Coulomb’s friction angle for ext. paths $\phi_E$ |
|Q(18)| = 0 : if the stress state is not at the criterion apex |
|:::| = 1 : if the stress state is at the criterion apex | 
|Q(19)| = number of sub-intervals used for the integration |
|Q(20)| = memory of localisation calculated during the re-meshing |
|Q(22)| = ? |
|Q(23)| = ? |
|Q(24)| = damage variable |
|Q(25)| = x plastic deformation |
|Q(26)| = y plastic deformation |
|Q(27)| = z plastic deformation |
|Q(28)| = xy plastic deformation |
|Q(29)| = Young modulus |
|Q(30)| = dilatancy angle in compression |
|Q(31)| = dilatancy angle in extension |
|Q(31)$\rightarrow$ Q(43)| = reserved for bifurcation |