Merge pull request #20 from 3MAH/add_plastic_law

Add plastic law for full quadratic anisotropic and DFA equivalent stress

include/simcoon/Continuum_mechanics/Functions/criteria.hpp

@@ -134,7 +134,7 @@ arma::vec dTresca_stress(const arma::vec &v);
     P_{ani} = \left( \begin{array}{cccccc}
         1 & -1/2 & -1/2 & 0 & 0 & 0 \\
         -1/2 & 1 & -1/2 & 0 & 0 & 0 \\
-        -1/2 & -1/2 & -1/2 & 0 & 0 & 0 \\
+        -1/2 & -1/2 & 1 & 0 & 0 & 0 \\
         0 & 0 & 0 & 3 & 0 & 0 \\
         0 & 0 & 0 & 0 & 3 & 0 \\
         0 & 0 & 0 & 0 & 0 & 3 \end{array} \right)
@@ -150,24 +150,24 @@ arma::mat P_Ani(const arma::vec &P_params);
  * @brief Provides an anisotropic configurational tensor considering the quadratic Hill yield criterion \cite Hill.48 in the Voigt format (6x6 matrix), given its vector representation
  * @param P_params
  * @return The anisotropic Hill48 configurational tensor (arma::mat)
- * @details The vector of parameters must be constituted of 5 values, respectively:
-    \f$ F^*,G^*,H^*,L,M,N \f$, such that
+ * @details The vector of parameters must be constituted of 6 values, respectively:
+    \f$ F,G,H,L,M,N \f$, such that
 \f[
     P_{Hill48} = \left( \begin{array}{cccccc}
         G + H & -H & -G & 0 & 0 & 0 \\
         -H & F + H & -F & 0 & 0 & 0 \\
         -G & -F & F + G & 0 & 0 & 0 \\
-        0 & 0 & 0 & 2 \, L & 0 & 0 \\
+        0 & 0 & 0 & 2 \, N & 0 & 0 \\
         0 & 0 & 0 & 0 & 2 \, M & 0 \\
-        0 & 0 & 0 & 0 & 0 & 2 \, N \end{array} \right)
+        0 & 0 & 0 & 0 & 0 & 2 \, L \end{array} \right)
 \f]
     considering the input stress \f$ \mathbf{\sigma} \f$.
     Note that the equivalent anisotropic Hill 1948) tensor is written as : \f$ \sigma^{eq}_{ani} = \sqrt{ \mathbf{\sigma} : P : \mathbf{\sigma} } \f$
     which reduces to : 
 \f[ 
     \begin{align}    
-    \sigma^{H48} & = \left( F\, \left( \sigma_{11} - \sigma_{22} \right)^2 + G\, \left( \sigma_{11} - \sigma_{33} \right)^2 + H\, \left( \sigma_{22} - \sigma_{33} \right)^2 \right. \\
-     & + \left. 2\,L\,\sigma_{12}^2 + 2\,M\,\sigma_{13}^2 + 2\,N\,\sigma_{23}^2 \right)^{1/2}
+    \sigma^{H48} & = \left( H\, \left( \sigma_{11} - \sigma_{22} \right)^2 + G\, \left( \sigma_{11} - \sigma_{33} \right)^2 + F\, \left( \sigma_{22} - \sigma_{33} \right)^2 \right. \\
+     & + \left. 2\,N\,\sigma_{12}^2 + 2\,M\,\sigma_{13}^2 + 2\,L\,\sigma_{23}^2 \right)^{1/2}
     \end{align}
 \f]
     Considering the Mises equivalent strain 
@@ -195,13 +195,53 @@ So that \f$ F = H = G = 1/2 \f$, \f$ = L = M = N = 3/2 \f$
 */
 arma::mat P_Hill(const arma::vec &P_params);
 
+/**
+ * @brief Provides an anisotropic configurational tensor considering the quadratic Deshpande–Fleck–Ashby yield criterion \cite Deshpande, V. S., Fleck, N. A. and Ashby, M. F. (2001) in the Voigt format (6x6 matrix), given its vector representation
+ * @param P_params
+ * @return The anisotropic DFA configurational tensor (arma::mat)
+ * @details The vector of parameters must be constituted of 7 values, respectively:
+    \f$ F,G,H,L,M,N,K \f$, such that
+\f[
+    P_{DFA} = \left( \begin{array}{cccccc}
+        G + H + K & -H + K & -G + K & 0 & 0 & 0 \\
+        -H + K & F + H + K & -F + K & 0 & 0 & 0 \\
+        -G + K & -F + K & F + G + K & 0 & 0 & 0 \\
+        0 & 0 & 0 & 2 \, N & 0 & 0 \\
+        0 & 0 & 0 & 0 & 2 \, M & 0 \\
+        0 & 0 & 0 & 0 & 0 & 2 \, L \end{array} \right)
+\f]
+    considering the input stress \f$ \mathbf{\sigma} \f$.
+    Note that the equivalent anisotropic tensor is written as : \f$ \sigma^{eq}_{ani} = \sqrt{ \mathbf{\sigma} : P : \mathbf{\sigma} } \f$
+    which reduces to : 
+\f[ 
+    \begin{align}    
+    \sigma^{DFA} & = \left( H\, \left( \sigma_{11} - \sigma_{22} \right)^2 + G\, \left( \sigma_{11} - \sigma_{33} \right)^2 + F\, \left( \sigma_{22} - \sigma_{33} \right)^2 \right. \\
+     & + \left. 2\,L\,\sigma_{12}^2 + 2\,M\,\sigma_{13}^2 + 2\,N\,\sigma_{23}^2 \right)^{1/2} + K \left( \sigma_{11} + \sigma_{22} + \sigma_{33} \right)^2
+    \end{align}
+\f]
+    Considering the full anisotric formulation:
+\f[ 
+    \begin{align}
+    \sigma^{ani} & = \left( P_{11}\,\sigma_{11}^2 + P_{22}\,\sigma_{22}^2 + P_{33}\,\sigma_{33}^2 \right. \\
+     & + 2\,P_{12}\,\sigma_{11}\,\sigma_{22} + 2\,P_{13}\,\sigma_{11}\,\sigma_{33} + 2\,P_{23}\,\sigma_{22} \sigma_{33} \\
+     & + \left. 2\,P_{44}\,\sigma_{12}^2 + 2\,P_{55}\,\sigma_{13}^2 + 2\,P_{66}\,\sigma_{23}^2 \right)^{1/2}
+    \end{align}
+\f]
+    the above matrix is identified.
+ * @code 
+        vec P_params = {0.5,0.6,0.7,1.5,1.5,1.6,1.};
+        mat P = P_Hill(P_params);
+ * @endcode
+*/
+arma::mat P_DFA(const arma::vec &P_params);
+
 /**
  * @brief Provides the anisotropic equivalent stress, given the stress in a vector format and given a configurational tensor P
  * @param sigma, H
  * @return The anisotropic equivalent stress (double)
  * @details Returns anisotropic equivalent stress \f$ \sigma^{eq}_{ani} = \sqrt{ \mathbf{\sigma} : P : \mathbf{\sigma} } \f$.
  * @code 
-        vec P_params = {0.5,0.6,0.7,3.,3.,3.2};
+        vec P_params = {0.5,0.6,0.7,1.5,1.5,1.6};
         mat P = P_Hill(P_params);
         vec sigma = randu(6);
         double sigma_ani = Eq_stress_P(sigma,P_Hill);
@@ -218,7 +258,7 @@ double Eq_stress_P(const arma::vec &sigma, const arma::mat &H);
     \sigma^{eq}_{ani} = \frac{\mathbf{H} : \mathbf{\sigma}}{ \sqrt{ \mathbf{\sigma} : P : \mathbf{\sigma} } }
 \f]
  * @code 
-        vec P_params = {0.5,0.6,0.7,3.,3.,3.2};
+        vec P_params = {0.5,0.6,0.7,1.5,1.5,1.6};
         mat P = P_Hill(P_params);
         vec sigma = randu(6);
         vec dsigma_ani = dEq_stress_P(sigma,P_Hill);
@@ -231,9 +271,9 @@ arma::vec dEq_stress_P(const arma::vec &sigma, const arma::mat &H);
  * @param sigma, P_params
  * @return The Hill48 anisotropic equivalent stress (double)
  * @details Returns the Hill48 anisotropic equivalent stress \f$ \sigma^{eq}_{ani} = \sqrt{ \mathbf{\sigma} : P : \mathbf{\sigma} } \f$.
- * see the function P_hill() for more details to obtain the tensor H from the set of parameters (F,G,H,L,M,N).
+ * see the function P_Hill() for more details to obtain the tensor H from the set of parameters (F,G,H,L,M,N).
  * @code 
-        vec P_params = {0.5,0.6,0.7,3.,3.,3.2};
+        vec P_params = {0.5,0.6,0.7,1.5,1.5,1.6};
         vec sigma = randu(6);
         double sigma_Hill = Hill_stress(sigma, P_params);
  * @endcode
@@ -248,15 +288,46 @@ double Hill_stress(const arma::vec &sigma, const arma::vec &P_params);
 \f[
     \sigma^{eq}_{ani} = \frac{\mathbf{H} : \mathbf{\sigma}}{ \sqrt{ \mathbf{\sigma} : P : \mathbf{\sigma} } }
 \f]
-see the function P_hill() for more details to obtain the tensor H from the set of parameters (F,G,H,L,M,N).
+see the function P_Hill() for more details to obtain the tensor H from the set of parameters (F,G,H,L,M,N).
  * @code 
-        vec P_params = {0.5,0.6,0.7,3.,3.,3.2};
+        vec P_params = {0.5,0.6,0.7,1.5,1.5,1.6};
         vec sigma = randu(6);
-        double sigma_Hill = Hill_stress(sigma, P_params);
+        vec dsigma_Hill = dHill_stress(sigma, P_params);
  * @endcode
 */
 arma::vec dHill_stress(const arma::vec &, const arma::vec &);
 
+/**
+ * @brief Provides the \cite Deshpande, V. S., Fleck, N. A. and Ashby, M. F. (2001) (DFA) anisotropic equivalent stress, given the stress in a vector format and a vector of parameters (F,G,H,L,M,N,K)
+ * @param sigma, P_params
+ * @return The DFA anisotropic equivalent stress (double)
+ * @details Returns the DFA anisotropic equivalent stress \f$ \sigma^{eq}_{ani} = \sqrt{ \mathbf{\sigma} : P : \mathbf{\sigma} } \f$.
+ * see the function P_DFA() for more details to obtain the tensor H from the set of parameters (F,G,H,L,M,N,K).
+ * @code 
+        vec P_params = {0.5,0.6,0.7,1.5,1.5,1.6,1.2};
+        vec sigma = randu(6);
+        double sigma_DFA = DFA_stress(sigma, P_params);
+ * @endcode
+*/
+double DFA_stress(const arma::vec &sigma, const arma::vec &P_params);
+
+/**
+ * @brief Provides the derivative of the \cite Deshpande, V. S., Fleck, N. A. and Ashby, M. F. (2001) (DFA) anisotropic equivalent stress, given the stress in a vector format and a vector of parameters (F,G,H,L,M,N,K)
+ * @param sigma, P_params
+ * @return The derivative of the DFA anisotropic equivalent stress (arma::vec)
+ * @details Returns the derivative of the DFA anisotropic equivalent stress \f$ \frac{\partial \mathbf{\sigma}^{ani}}{\partial \mathbf{\sigma}} \f$, considering
+\f[
+    \sigma^{eq}_{ani} = \frac{\mathbf{H} : \mathbf{\sigma}}{ \sqrt{ \mathbf{\sigma} : P : \mathbf{\sigma} } }
+\f]
+see the function P_DFA() for more details to obtain the tensor H from the set of parameters (F,G,H,L,M,N,K).
+ * @code 
+        vec P_params = {0.5,0.6,0.7,1.5,1.5,1.6,1.2};
+        vec sigma = randu(6);
+        vec dsigma_DFA = dDFA_stress(sigma, P_params);
+ * @endcode
+*/
+arma::vec dDFA_stress(const arma::vec &, const arma::vec &);
+
 /**
  * @brief Provides the anisotropic equivalent stress, given the stress in a vector format and a vector of parameters \f$ ( P_{11},P_{22},P_{33},P_{12},P_{13},P_{23},P_{44},P_{55},P_{66} ) \f$
  * @param sigma, P_params

include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/Ani_chaboche_ccp.hpp

@@ -0,0 +1,50 @@
+/* This file is part of simcoon.
+ 
+ simcoon is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+ 
+ simcoon is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+ 
+ You should have received a copy of the GNU General Public License
+ along with simcoon.  If not, see <http://www.gnu.org/licenses/>.
+ 
+ */
+
+///@file plastic_isotropic_ccp.hpp
+///@brief User subroutine for elastic-plastic materials in 1D-2D-3D case
+///@brief This subroutines uses a convex cutting plane algorithm
+///@brief Isotropic hardening with a power-law hardenig is considered
+///@version 1.0
+
+#pragma once
+#include <armadillo>
+
+namespace simcoon{
+
+///@brief The elastic-plastic UMAT with isotropic hardening requires 6 constants:
+
+///@brief props[1] : Young modulus
+///@brief props[2] : Poisson ratio
+///@brief props[3] : CTE
+///@brief props[4] : J2 equivalent yield stress limit : sigmaY
+///@brief props[5] : hardening parameter k
+///@brief props[6] : exponent m
+
+///@brief The elastic-plastic UMAT with isotropic hardening requires 8 statev:
+///@brief statev[0] : T_init : Initial temperature
+///@brief statev[1] : Accumulative plastic parameter: p
+///@brief statev[2] : Plastic strain 11: EP(0,0)
+///@brief statev[3] : Plastic strain 22: EP(1,1)
+///@brief statev[4] : Plastic strain 33: EP(2,2)
+///@brief statev[5] : Plastic strain 12: EP(0,1)
+///@brief statev[6] : Plastic strain 13: EP(0,2)
+///@brief statev[7] : Plastic strain 23: EP(1,2)
+
+void umat_ani_chaboche_CCP(const arma::vec &, const arma::vec &, arma::vec &, arma::mat &, arma::mat &, arma::vec &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &,const double &,const double &, double &, double &, double &, double &, const int &, const int &, const bool &, const int &, double &);
+
+} //namespace simcoon

include/simcoon/Continuum_mechanics/Umat/Mechanical/Plasticity/DFA_chaboche_ccp.hpp

@@ -0,0 +1,50 @@
+/* This file is part of simcoon.
+ 
+ simcoon is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+ 
+ simcoon is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+ 
+ You should have received a copy of the GNU General Public License
+ along with simcoon.  If not, see <http://www.gnu.org/licenses/>.
+ 
+ */
+
+///@file plastic_isotropic_ccp.hpp
+///@brief User subroutine for elastic-plastic materials in 1D-2D-3D case
+///@brief This subroutines uses a convex cutting plane algorithm
+///@brief Isotropic hardening with a power-law hardenig is considered
+///@version 1.0
+
+#pragma once
+#include <armadillo>
+
+namespace simcoon{
+
+///@brief The elastic-plastic UMAT with isotropic hardening requires 6 constants:
+
+///@brief props[1] : Young modulus
+///@brief props[2] : Poisson ratio
+///@brief props[3] : CTE
+///@brief props[4] : J2 equivalent yield stress limit : sigmaY
+///@brief props[5] : hardening parameter k
+///@brief props[6] : exponent m
+
+///@brief The elastic-plastic UMAT with isotropic hardening requires 8 statev:
+///@brief statev[0] : T_init : Initial temperature
+///@brief statev[1] : Accumulative plastic parameter: p
+///@brief statev[2] : Plastic strain 11: EP(0,0)
+///@brief statev[3] : Plastic strain 22: EP(1,1)
+///@brief statev[4] : Plastic strain 33: EP(2,2)
+///@brief statev[5] : Plastic strain 12: EP(0,1)
+///@brief statev[6] : Plastic strain 13: EP(0,2)
+///@brief statev[7] : Plastic strain 23: EP(1,2)
+
+void umat_dfa_chaboche_CCP(const arma::vec &, const arma::vec &, arma::vec &, arma::mat &, arma::mat &, arma::vec &, const arma::mat &, const int &, const arma::vec &, const int &, arma::vec &, const double &, const double &,const double &,const double &, double &, double &, double &, double &, const int &, const int &, const bool &, const int &, double &);
+
+} //namespace simcoon

src/Continuum_mechanics/Functions/criteria.cpp

@@ -159,6 +159,31 @@ mat P_Hill(const vec &params) {
     return P;
 }
 
+mat P_DFA(const vec &params) { //Deshpande–Fleck–Ashby 
+    assert(params.n_elem == 7); //F,G,H,L,M,N,K
+    mat P = zeros(6,6);
+    //param(0) = F
+    //param(1) = G
+    //param(2) = H
+    //param(3) = L
+    //param(4) = M
+    //param(5) = N
+    P(0,0) = params(1) + params(2) + params(6); //P_11 = G+H+K
+    P(1,1) = params(0) + params(2) + params(6); //P_22 = F+H+K
+    P(2,2) = params(0) + params(1) + params(6); //P_33 = F+G+K
+    P(0,1) = -1.*params(2) + params(6); //P_12 = -H+K
+    P(1,0) = -1.*params(2) + params(6); //P_12 = -H+K
+    P(0,2) = -1.*params(1) + params(6); //P_13 = -G+K
+    P(2,0) = -1.*params(1) + params(6); //P_13 = -G+K
+    P(1,2) = -1.*params(0) + params(6); //P_23 = -F+K
+    P(2,1) = -1.*params(0) + params(6); //P_23 = -F+K
+    P(3,3) = 2.*params(5); //P_44 = N
+    P(4,4) = 2.*params(4); //P_55 = M
+    P(5,5) = 2.*params(3); //P_66 = L
+
+    return P;
+}
+
 double Eq_stress_P(const vec &v, const mat &H) {
 
     if (norm(v,2) > sim_iota) {
@@ -188,6 +213,16 @@ vec dHill_stress(const vec &v, const vec &params) {
    mat P = P_Hill(params);
    return dEq_stress_P(v,P);
 }
+
+double DFA_stress(const vec &v, const vec &params) {
+    mat P = P_DFA(params);
+    return Eq_stress_P(v,P);
+}
+
+vec dDFA_stress(const vec &v, const vec &params) {
+   mat P = P_DFA(params);
+   return dEq_stress_P(v,P);
+}
 
 double Ani_stress(const vec &v, const vec &params) {
     mat P = P_Ani(params);