====== EP-MOMAS ======

===== Description =====

Elasto-plastic constitutive law for solid elements at constant temperature (non-associated) with linear elasticity. Isotropic softening of cohesion is possible.

==== The model ====

This law is used for mechanical analysis of elasto-plastic isotropic porous media undergoing large strains. 

==== Files ====

Prepro: LMOMA.F \\
Lagamine: MOMA2EA.F
===== Availability =====

|Plane stress state| NO |
|Plane strain state| YES |
|Axisymmetric state| YES |
|3D state| NO |
|Generalized plane state| NO |

===== Input file =====

==== Parameters defining the type of constitutive law ====

^ Line 1 (2I5, 60A1)^^
|IL|Law number|
|ITYPE| 586|
|COMMNT| Any comment (up to 60 characters) that will be reproduced on the output listing|

==== Integer parameters ====

^ Line 1 (7I5) ^^
|NINTV| Number of sub-steps used to integrate numerically the constitutive equation in a time step |
|:::| = 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 |
|:::| ≠ 0 : Use of effective stresses in the constitutive law (See annex 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 |
|KMETH| = 2 : Actualised VGRAD integration |
|:::| = 3 : Mean VGRAD integration (Default value) |

==== Real parameters ====

^ Line 1 (7G10.0) ^^
|E| Young's elastic modulus |
|ANU| Poisson's ratio |
|PSIC| Dilatancy angle (in degrees) for compressive paths |
|PSIE| Dilatancy angle (in degrees) for extensive paths (iff 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 |
^ Line 2 (3G10.0) ^^
|PHIC| Coulomb's angle (in degrees) for compressive paths |
|PHIE| Coulomb's angle (in degrees) for extensive paths |
|AN| Van Eekelen exponent (default value=-0.229) |
^ Line 3 (4G10.0) ^^
|COH0| Initial value of cohesion |
|ALPHA| Coefficient in the softening relationship |
|GAMMARP| |
|BIOPT||
===== Stresses =====

==== Number of stresses ====

6 for 3D state \\ 4 for the other cases

==== Meaning ====

The stresses are the components of CAUCHY stress tensor in global (X,Y,Z) coordinates. \\

For the 3-D state:
|SIG(1)|$\sigma_{xx}$|
|SIG(2)|$\sigma_{yy}$|
|SIG(3)|$\sigma_{zz}$|
|SIG(4)|$\sigma_{xy}$|
|SIG(5)|$\sigma_{xz}$|
|SIG(6)|$\sigma_{yz}$|

For the other cases:
|SIG(1)|$\sigma_{xx}$|
|SIG(2)|$\sigma_{yy}$|
|SIG(3)|$\sigma_{xy}$|
|SIG(4)|$\sigma_{zz}$|

===== State variables =====

==== Number of state variables ====

= 36 for 2D plane strain analysis with bifurcation criterion (ICBIF=1) \\ = 24 in all the other cases

==== List of state variables ====

|Q(1)| = 1 : Plane strain state |
|:::| Circumferential strain rate ($\dot{\varepsilon}_{\theta}$) in axisymmetrical state |
|Q(2)| Actualised specific mass |
|Q(3)| = 0 : Current state is elastic |
|:::| = 1 : Current state is elasto-plastic |
|Q(4)| Plastic work per unit volume ($W^p$) |
|Q(5)| Volumic variation |
|Q(6)| Equivalent strain n°1 : $\varepsilon_{eq1}=\int\Delta\dot{\varepsilon}_{eq}\;\Delta t$ |
|Q(7)| Equivalent strain indicator n°1 (Villote n°1) : $\alpha_1=\frac{\Delta\dot{\varepsilon}_{eq}\;\Delta t}{\varepsilon_{eq1}}$ |
|Q(8)| X deformation |
|Q(9)| Y deformation |
|Q(10)| Z deformation |
|Q(11)| XY deformation |
|Q(12)| Equivalent strain n°2 : $\varepsilon_{eq2}=\int\Delta\varepsilon_{eq}$ |
|Q(13)| Equivalent strain indicator n°2 (Villote n°2) $\alpha_2=\frac{\Delta\varepsilon_{eq}}{\varepsilon_{eq2}}$ |
|Q(14)| Actualised value of equivalent plastic strain $\varepsilon_{eq}^p$ |
|Q(15)| Actualised value of cohesion $c$ |
|Q(16)| Actualised value of Coulomb's frictional angle in compressive path ($\phi_C$) |
|Q(17)| Actualised value of Coulomb's frictional angle in extensive path ($\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)| |
|Q(21)| |
|Q(22)| Actualised value of volumetric plastic deformations |
|Q(23)| Second deviatoric strain increment invariant |
|Q(24)| Plastic loading index |
|Q(25)| Memory of localisation calculated during the re-meshing |
|Q(26)$\rightarrow$Q(36)| Reserved for bifurcation |

===== Formulation =====

==== Yield and flow surfaces ====

The stresses and stress invariants are : 
\[I_{\sigma} = \sigma_{ij}\quad ; \quad \hat{\sigma}_{ij}=\sigma_{ij}-\frac{I_{\sigma}}{3}\delta_{ij} \]
\[II_{\hat{\sigma}}=\sqrt{\frac{1}{2}\hat{\sigma}_{ij}\hat{\sigma}_{ij}}\quad ;\quad III_{\hat{\sigma}} = \frac{1}{3}\hat{\sigma}_{ij}\hat{\sigma}_{jk}\hat{\sigma}_{ki}\]
\[\beta =-\frac{1}{3}\sin^{-1}\left(\frac{3\sqrt{3}}{2}\frac{III_{\hat{\sigma}}}{II^3_{\hat{\sigma}}}\right)\]

=== Criterion with friction angle different from 0 (Drücker Prager or Van Eekelen) ===

The regular criterion is used if $I_{\sigma}-m'II_{\hat{\sigma}}<\frac{3c}{\tan\phi_c}$ : \[f=II_{\hat{\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 functions of $\phi_C$, $\phi_E$ and $n$.\\

The apex criterion is used if $I_{\sigma}-m'II_{\hat{\sigma}}\geq\frac{3c}{\tan\phi_C}$ : \[f=I_{\sigma}-\frac{3c}{\tan\phi_c}=0\] where $m'$ is the equivalent of $m$ but for the flow surface (i.e. $\phi$ is replaced by $\psi$ )

==== Softening ====

Softening is assumed to be represented by the evolution of cohesion as a function of the Von Mises equivalent plastic strain : \[\gamma^p=\sqrt{\hat{\varepsilon}_{ij}^p\hat{\varepsilon}_{ij}^p}\] where $\hat{\varepsilon}_{ij} = \varepsilon_{ij}-\frac{I_{\varepsilon}}{3}\delta_{ij}$ is the deviatoric strain tensor.\\

The following function is used : \begin{align*}f(\gamma^p)&=\left(1-(1-\alpha)\frac{\gamma^p}{\gamma_R^p}\right)^2\;quad\text{if}\quad 0<\gamma^p<\gamma^p_R \\ &= \alpha^2 \quad\text{if}\quad \gamma^p\geq \gamma_R^p \end{align*} where $\alpha$ and $\gamma_R^p$ are two models parameters.