====== ORTHOPLA ======
===== Description =====
Elasto‑plastic constitutive law for solid elements at constant temperature (non-associated) with linear anisotropic elasticity. Isotropic hardening/softening of friction angle and cohesion is possible.\\

This law is a combination of [[laws:ortho3d|ORTHO3D]] (for the elastic part of the law) and PLA3D (for the plastic behaviour). 

This law is only  used for mechanical analysis of elasto-plastic anisotropic porous media undergoing large strains.
==== Files ====
Prepro: LORTHOPLA.F \\
===== Formulation =====
==== Yield and flow surfaces ====
See [[laws:epplasol|PLASOL]] law
==== Hardening/softening ====
See [[laws:epplasol|PLASOL]] law (above)
==== Cohesion anisotropy with major principal stress orientation relative to bedding (IANISO = 0)====
The material cohesion depends on the angle  $\alpha_{\sigma_1}$ between the major compressive principal stress $\vec{\sigma_1}$ and the normal to the bedding plane  $\vec{e_3}$ :
\[
\alpha_{\sigma_1} = \arccos\left( \frac{\vec{e_3}\vec{\sigma_{1}}'}{||\vec{e_3}||\ ||\vec{\sigma_1}'||} \right)
\]

<imgcaption image1|Schematic view of the angle between the normal to bedding plane and the
direction of major principal stress>{{ 
 :laws:schematic_view_of_the_angle_between_the_normal_to_bedding_plane_and_the_direction_of_major_principal_stress.png?150  |}}</imgcaption>

Three cohesion values are defined ($c_{0^{\circ}}, c_{min}, c_{90^{\circ}}$), for major principal stress parallel $\alpha_{\sigma_1} = 0^{\circ}$ (perpendicular), perpendicular $\alpha_{\sigma_1} = 90^{\circ}$ (parallel) and with an angle of $\alpha_{\sigma_1, min}$ with respect to the normal to bedding plane (with respect to the bedding plane)  (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.)). Between those values, cohesion varies linearly with $\alpha_{\sigma_1}$. The mathematical expression of the cohesion is as follows:
\[
c = \max \left[\left( \frac{c_{min} - c_{0^{\circ}}}{\alpha_{\sigma_1, min}} \right)\alpha_{\sigma_1} + c_{0^{\circ}} ;  \left( \frac{c_{90^{\circ}} - c_{min}}{90^{\circ} - \alpha_{\sigma_1, min}} \right)\left( \alpha_{\sigma_1} - \alpha_{\sigma_1, min} \right)+ c_{0^{\circ}}  \right]
\]

<imgcaption image2|Schematic view of the cohesion evolution as a function of the angle between the
normal to bedding plane and the direction of major principal stress>{{  :laws:schematic_view_of_the_cohesion_evolution_as_a_function_of_the_angle_between_the_normal_to_bedding_plane_and_the_direction_of_major_principal_stress.png?400  |}}</imgcaption>

==== Cohesion anisotropy by microstructure fabric tensor (IANISO = 1)====
The material cohesion anisotropy is defined with a microstructure tensor $a_{ij}$, which is a measure of the material fabric. The eigenvectors of this tensor correspond to the preferred or principal material axes (orthotropic axes) $e_1, e_2, e_3$. The cohesion corresponds to the projection of this tensor on a generalized loading vector l ,  therefore the cohesion specifies the effect of load orientation relative to the material axes.
\[c = a_{ij}l_il_j\]
\[l_i = \sqrt{\frac{\sigma_{l1}^2 + \sigma_{l2}^2 +\sigma_{l3}^2}{\sigma_{ij}\sigma_{ij}}}\]

Where $\sigma_{ij}$ expressed in reference to the material axes. The cohesion can be expressed as : 
\[c = c_0\left( 1 + A_{ij}l_il_j\right)\]
\[A_{ij} = \frac{dev(a_{ij})}{c_0} =\frac{a_{ij}}{c_0} - \delta_{ij}\]
\[c_0 = \frac{a_{kk}}{3}\]

where $A_{ij}$ is a symmetric traceless operator, $A_{kk}=0$. The above expression can be generalized by considering higher order tensors :
\[c= c_0 \left( 1+A_{ij}l_il_j + b_1(A_{ij}l_il_j)^2 + b_2(A_{ij}l_il_j)^3 + … \right)\]
where $B_1,b_2,…$ are constants.\\

Considering cross-anisotropy, i.e. transverse isotropy, and refering the problem to the principal material axes implies $A_{ij} = 0$ for $i \neq j$, $A_{ii} = A_{11}+A_{22}+A_{33} = 0$, $A_{11} = A_{33}$ if the bedding plane is in ($e_1, e_3$) anisotropic plane, $A_{22} = -2A_{11}$, implying :
\[A_{ij}l_il_j = A_{l1}(1-3l_2^2)\]
where $A_{11}$ is the component of the microstructure operator $A_{ij}$ in the isotropic (bedding) plane. The late expression for cohesion becomes (Pardoen, 2015)((Pardoen, B. (2015)  Hydro-mechanical analysis of the fracturing induced by the excavation of nuclear waste repository galleries using shear banding. Thesis, Liège University.)):
\[c= c_0 \left( 1+A_{l1}(1-3l_2^2) + b_1A_{l1}^2(1-3l_2^2)^2 + b_2A_{l1}^3(1-3l_2^2)^3 + … \right)\]

The constants $c_0, A_{11}, b_1, b_2, …$ can be obtained from experimental data and laboratory tests. For uniaxial compression : $l_2 = \cos(\alpha)$ , $\alpha$ being the angle between the compression direction and the normal to the bedding plane (<imgref image1>), $\alpha=0^{\circ}$ if the loading is perpendicular to the bedding plane and $\alpha = 90^{\circ}$ if parallel. The cohesion and the uniaxial compressive straight $f_c$ are linked by $f_c = 2c\cos\phi /(1-\sin\phi)$. The constants can be obtained from results coming from tests performed for differents orientation $\alpha$ (<imgref image3>). For axial compression $\sigma_1$ with confinement $p_0$ (triaxial) : $l_2^2 = \frac{p_0^2\sin^2(\alpha)+\sigma_1^2\cos^2(\alpha)}{2p_0^2+\sigma_1^2}$.

<imgcaption image3|Schematic view of the cohesion and uniaxial compressive strength evolution as a
function of the angle between the normal to bedding plane and the direction of axial loading>{{ 
 :laws:schematic_view_of_the_cohesion_and_uniaxial_compressive_strength_evolution_as_a_function_of_the_angle_between_the_normal_to_bedding_plane_and_the_direction_of_axial_loading.png?600  |}}</imgcaption>
==== Vicoplasticity ====
See [[laws:epplasol|PLASOL]] \\
Remark : For anisotropic Biot’s coeffcient, the deviatoric stress  is calculated from the effective stresses (more details about this anisotropy are available in the definition of element CSOL2 and ISOL=9 in [[appendices:a7|Appendix 7]]).

===== 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| 608 |
|COMMENT| Any comment (up to 60 characters) that will be reproduced on the output listing|
==== Integer parameters ====
^ Line 1 (14I5) ^^
|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:a7|Appendix 7]] |
|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 deviator 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) |
|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 = c_0 + AK1.log (s + 0.0001) + AK2$) |
|:::| 3 : Young's modulus is a function of suction |
|IBEDDING| Bedding orientation (normal to the bedding plane) |
|:::| 1 : bedding plane in $e_1e_2$ anisotropic plane (normal $e_3$) |
|:::| 2 : bedding plane in $e_1e_3$ anisotropic plane (normal $e_2$) |
|:::| 3 : bedding plane in $e_2e_3$ anisotropic plane (normal $e_1$) |
|IANISO| = 0 : anisotropy of cohesion with major principal stress orientation relative to bedding |
|:::| = 1 : anisotropy of cohesion by microstructure fabric tensor |
|IVISCO| = 0 : nothing |
|:::| 1 to 3 : viscoplastic model |
|IDAM| = 0 : nothing |
|:::| 1 : Damage on elastic properties |
|IHSS| = 0 : nothing |
|:::| 1 : Coupling stiffness-deformation |

==== Real parameters ====
{{  :laws:ortho3d.png?300  |}}
^ Line 1 (3G10.0) ^^
|ALPHA|Angle of rotation of the anisotropic axis around Z axis (see figure above)|
|THETA|Angle of rotation of the anisotropic axis around e1 axis (see figure above)|
|PHI|Angle of rotation of the anisotropic axis around e2 axis (see figure above)|
^ Line 2 (9G10.0) ^^
|E1|Elastic Young modulus E($e_1$)|
|E2|Elastic Young modulus E($e_2$)|
|E3|Elastic Young modulus E($e_3$)|
|G12|Elastic shear modulus G($e_1e_2$)|
|G13|Elastic shear modulus G($e_1e_3$)|
|G23|Elastic shear modulus G($e_1e_3$)|
|AE1|E1 = E1 + AE1 * (-SIG(M)), SIG(M) is confinement pressure |
|AE2|E2 = E2 + AE2 * (-SIG(M)) |
|AE3|E3 = E3 + AE3 * (-SIG(M)) |
^ Line 3 (5G10.0) ^^
|ANU12|Poisson ratio NU($e_1e_2$)|
|ANU13|Poisson ratio NU($e_1e_3$)|
|ANU23|Poisson ratio NU($e_2e_3$)|
|RHO|Specific mass|
|DIV|Size of sub-steps for computation of NINTV (only if NINTV=0; Default value=5.D-3)|
^ Line 4 (7G10.0) ^^
|PSIC|Dilatancy angle (in degrees) for compressive paths|
|PSIE|Dilatancy angle (in degrees) for extensive paths (ssi ILODEG=2)|
|PHMPS|Constant value for definition of |
|BIOPT|Bifurcation computation parameter|
|AK1|Capillary cohesion first parameter|
|AK2|Capillary cohesion second parameter|
|DECCOH|Cohesion hardening shifting|
^ Line 5 (7G10.0) ^^
|PHICF|Final Coulomb's angle (in degrees) for compressive paths|
|PHIEF|Final Coulomb’s angle (in degrees) for extensive paths (ssi ILODEF = 2)|
|RAYPHIC|Ratio between initial and residual friction angle for compressive paths|
|BPHI|Only if there is hardening/softening|
|AN|Van Eekelen exponent (default value=-0.229)|
|DECPHI|Coulomb’s angle hardening shifting|
|RAYPHIE|Ratio between initial and residual friction angle for extensive paths (ssi ILODEF = 2)|
^ Line 6 (6G10.0) ^^
|COHF0|Residual value of cohesion for major principal stress 1 perpendicular to the bedding plane (parallel to the normal to the bedding plane) (if IANISO = 0)|
|:::|COHC0 = c0 (if IANISO = 1), see [[laws:orthopla#Cohesion anisotropy by microstructure fabric tensor|Cohesion anisotropy by microstructure fabric tensor]]|
|COHFMIN|Minimal residual value of cohesion (if IANISO = 0)|
|:::|COHAISO = A11 (if IANISO = 1), see [[laws:orthopla#Cohesion anisotropy by microstructure fabric tensor|Cohesion anisotropy by microstructure fabric tensor]]|
|COHF90|Residual value of cohesion for major principal stress 1 parallel to the bedding plane (perpendicular to the normal to the bedding plane) (if IANISO = 0)|
|:::|COHB1 = b1 (if IANISO = 1), see [[laws:orthopla#Cohesion anisotropy by microstructure fabric tensor|Cohesion anisotropy by microstructure fabric tensor]]|
|ANGLEMIN|Angle between the normal to the bedding plane and the major principal stress 1 for which the cohesion is minimum (if IANISO = 0)|
|:::|COHB2 = b2 (if IANISO = 1), see [[laws:orthopla#Cohesion anisotropy by microstructure fabric tensor|Cohesion anisotropy by microstructure fabric tensor]]|
|RAYCOH|Ratio between initial and residual cohesion|
|BCOH|Only if there is hardening/softening|
^ Line 7 (4G10.0) (Only if IHSS = 1) ^^ 
|E1F|Final elastic Young modulus E($e_{1f}$)|
|E2F|Final elastic Young modulus E($e_{2f}$)|
|E3F|Final elastic Young modulus E($e_{3f}$)|
|Gamma7|equivalent strain at which the Young's modulus has reduced to 0.7 times |
|Aa|Fitting parameter |
^ Line 8 (7G10.0) (Only if IECPS = 2 or 3) ^^ 
|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 9 (2G10.0) (Only if IDAM = 1) ^^
|P|Parameter controlling the damage evolution rate|
|YD0|Initial 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 ====
= 48 : for 2D plane strain analysis with bifurcation criterion (ICBIF=1)\\
= 36 : 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 [[laws:epplasol|PLASOL]]) |
|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 for bedding perpendicular to 1st principal stress (COHF0, if IANISO = 0) or COHCO ( if IANISO=1) |
|Q(16)| = actualized value of cohesion $c$|
|Q(17)| = actualised value of Coulomb’s friction angle for compr. paths $\phi_C$|
|Q(18)| = actualised value of Coulomb’s friction angle for ext. paths $\phi_E$ |
|Q(19)| = 0 : if the stress state is not at the criterion apex |
|:::| = 1 : if the stress state is at the criterion apex | 
|Q(20)| = number of sub-intervals used for the integration |
|Q(21)| = memory of localisation calculated during the re-meshing |
|Q(22)| = ? |
|Q(23)| = ? |
|Q(24)| = ORIENTBED |
|Q(25)| = dilatancy angle in compression |
|Q(26)| = dilatancy angle in extension |
|Q(27)| = damage variable |
|Q(28)| = x plastic deformation |
|Q(29)| = y plastic deformation |
|Q(30)| = z plastic deformation |
|Q(31)| = xy plastic deformation |
|Q(32)| = ? |
|Q(33)| = ? |
|Q(34)$\rightarrow$ Q(36)| = reserved for small strain stiffness (E1, E2, E3) |
|Q(37)$\rightarrow$ Q(48)| = reserved for bifurcation |