SmallStrainPlasticityMaterialModels.hpp

/** \file SmallStrainPlasticityMaterialModels.hpp
 * \ingroup small_strain_plasticity
 * \brief Small Strain plasticity material models
 * \example SmallStrainPlasticityMaterialModels.hpp
 */
 
/*
 * This file is part of MoFEM.
 * MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM 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 Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
 
#ifndef __SMALL_STRAIN_PLASTICITY_MATERIAL_MODELS_HPP
#define __SMALL_STRAIN_PLASTICITY_MATERIAL_MODELS_HPP
 
#ifndef WITH_ADOL_C
  #error "MoFEM need to be compiled with ADOL-C"
#endif
 
/** \brief J2 plasticity (Kinematic Isotropic (Linear) Hardening)
* \ingroup small_strain_plasticity
*/
struct SmallStrainJ2Plasticity: public SmallStrainPlasticity::ClosestPointProjection {
 
  SmallStrainJ2Plasticity():
  SmallStrainPlasticity::ClosestPointProjection() {
    internalVariables.resize(7,false);
    internalVariables.clear();
  }
 
  double mu;      //< Lamé parameter
  double lambda;  //< Lamé parameters
  double H;       //< Isotropic hardening
  double K;       //< Kinematic hardening
  double phi;     //< Combined isotropic/kinematic hardening
  double sIgma_y; //< Yield stress
 
  /** \brief Available energy
 
  \f[
  \psi =
  \frac{1}{2} \lambda \textrm{tr}[\varepsilon]^2 + \mu \varepsilon : \varepsilon
  +
  \sigma_y\alpha
  +
  \frac{1}{2} \phi H \alpha^2
  +
  \frac{1}{2} (1-\phi)K \beta^2
  \f]
 
  Isotropic hardening variable \f$\alpha\f$ is in first index of
  internal variables vector. Kinematic hardening variable is in
  the remaining indices of internal variables vector.
 
  */
  virtual PetscErrorCode freeHemholtzEnergy() {
    PetscFunctionBegin;
    //elastic part
    adouble tR = 0;
    for(int dd = 0;dd<3;dd++) {
      tR += a_sTrain[dd]-a_plasticStrain[dd];
    }
    a_w = 0.5*lambda*tR*tR;
    adouble t;
    int dd = 0;
    for(;dd<3;dd++) {
      t = a_sTrain[dd]-a_plasticStrain[dd];
      a_w += mu*t*t;
    }
    for(;dd<6;dd++) {
      t = a_sTrain[dd]-a_plasticStrain[dd];
      a_w += 0.5*mu*t*t;
    }
    //plastic part
    //isotropic
    a_w += a_internalVariables[0]*sIgma_y+0.5*phi*H*a_internalVariables[0]*a_internalVariables[0];
    //kinematic
    for(unsigned int dd = 1;dd<a_internalVariables.size();dd++) {
     a_w += 0.5*(1-phi)*K*a_internalVariables[dd]*a_internalVariables[dd];
    }
    PetscFunctionReturn(0);
  }
 
  /** \brief Auxiliary function
 
  \f[
  \eta=\textrm{dev}[\sigma]-\overline{\beta}
  \f]
 
  \f[
  f = \frac{3}{2} \eta:\eta
  \f]
 
  This is \f$3J_2\f$.
 
  */
  adouble t,f;
  PetscErrorCode evalF() {
    PetscFunctionBegin;
    adouble tR = 0;
    for(int dd = 0;dd<3;dd++) {
      tR += a_sTress[dd];
    }
    // Fix the constant to mach uniaxial test
    const double c = 3./2.;
    tR /= 3.;
    f = 0;
    for(int dd = 0;dd<3;dd++) {
      t = (a_sTress[dd]-tR)-a_internalFluxes[dd+1];
      f += c*t*t;
    }
    for(int dd = 3;dd<6;dd++) {
      t = a_sTress[dd]-a_internalFluxes[dd+1];
      f += c*2.*t*t; // s:s off diagonal terms need to go twice
    }
    PetscFunctionReturn(0);
  }
 
  /** \brief Evaluate yield function
 
  \f[
  y = \sqrt{f} - \overline{\alpha}
  \f]
  where \f$f\f$ is defined in \ref evalF.
 
  */
  virtual PetscErrorCode yieldFunction() {
    PetscFunctionBegin;
    ierr = evalF(); CHKERRQ(ierr);
    a_y = sqrt(f) - a_internalFluxes[0];
    PetscFunctionReturn(0);
  }
 
  /** \brief Associated flow potential
 
  See \ref yieldFunction.
 
  */
  virtual PetscErrorCode flowPotential() {
    PetscFunctionBegin;
    ierr = evalF(); CHKERRQ(ierr);
    a_h = sqrt(f) - a_internalFluxes[0];
    PetscFunctionReturn(0);
  }
 
};
 
 
 
/** \brief Small strain plasticity with paraboloidal yield criterion (Isotropic Hardening)
* \ingroup small_strain_plasticity
*/
struct SmallStrainParaboloidalPlasticity: public SmallStrainPlasticity::ClosestPointProjection {
 
  SmallStrainParaboloidalPlasticity():
  SmallStrainPlasticity::ClosestPointProjection() {
    internalVariables.resize(2,false);
    internalVariables.clear();
  }
 
  double mu;      //< Lamé parameter
  double lambda;  //< Lamé parameters
  double nup;     //< Plastic Poisson’s ratio
  double Ht, Hc;  //< Isotropic hardening for tension and compression
  double sIgma_yt, sIgma_yc; //< Yield stress in tension and compression
  double nt, nc; //< Hardening parameters for tension and compression
 
 
 
  /** \brief Available energy
 
  \f[
  \psi =
  \frac{1}{2} \lambda \textrm{tr}[\varepsilon]^2 + \mu \varepsilon : \varepsilon
  +
  \sigma_{yt}\alpha_0
  +
  \frac{1}{2} H_t \alpha_0^2
  +
 \sigma_{yc}\alpha_1
  +
  \frac{1}{2} H_c \alpha_1^2
  \f]
 
  //Energy with linear hardening
//  */
//  virtual PetscErrorCode freeHemholtzEnergy() {
//    PetscFunctionBegin;
//    //elastic part
//    adouble tR = 0;
//    for(int dd = 0;dd<3;dd++) {
//      tR += a_sTrain[dd]-a_plasticStrain[dd];
//    }
//    a_w = 0.5*lambda*tR*tR;
//    adouble t;
//    int dd = 0;
//    for(;dd<3;dd++) {
//      t = a_sTrain[dd]-a_plasticStrain[dd];
//      a_w += mu*t*t;
//    }
//    for(;dd<6;dd++) {
//      t = a_sTrain[dd]-a_plasticStrain[dd];
//      a_w += 0.5*mu*t*t;
//    }
//    //plastic part
//    //isotropic
//    a_w += a_internalVariables[0]*sIgma_yt + 0.5*Ht*a_internalVariables[0]*a_internalVariables[0];
//    a_w += a_internalVariables[1]*sIgma_yc + 0.5*Hc*a_internalVariables[1]*a_internalVariables[1];
//    PetscFunctionReturn(0);
//  }
 
 
  //Energy with exponential hardening
  //  */
  virtual PetscErrorCode freeHemholtzEnergy() {
    PetscFunctionBegin;
    //elastic part
    adouble tR = 0;
    for(int dd = 0;dd<3;dd++) {
      tR += a_sTrain[dd]-a_plasticStrain[dd];
    }
    a_w = 0.5*lambda*tR*tR;
    adouble t;
    int dd = 0;
    for(;dd<3;dd++) {
      t = a_sTrain[dd]-a_plasticStrain[dd];
      a_w += mu*t*t;
    }
    for(;dd<6;dd++) {
      t = a_sTrain[dd]-a_plasticStrain[dd];
      a_w += 0.5*mu*t*t;
    }
    //plastic part
    //isotropic
//    a_w += a_internalVariables[0]*sIgma_yt + 0.5*Ht*a_internalVariables[0]*a_internalVariables[0];
//    a_w += a_internalVariables[1]*sIgma_yc + 0.5*Hc*a_internalVariables[1]*a_internalVariables[1];
    a_w +=(sIgma_yt+Ht)*a_internalVariables[0] + (Ht/nt)*exp(-nt*a_internalVariables[0]);
    a_w +=(sIgma_yc+Hc)*a_internalVariables[1] + (Hc/nc)*exp(-nc*a_internalVariables[1]);
 
    PetscFunctionReturn(0);
  }
 
 
 
  /** \brief Auxiliary function
  \f[
   I_1 = \textrm{tr} (\boldsymbol{\sigma})
  \f]
 
  \f[
  \eta=\textrm{dev}[\boldsymbol{\sigma}]
  \f]
 
  \f[
   J_2 = \frac{1}{2} \eta:\eta
  \f]
 
  */
  adouble I1,J2;
  PetscErrorCode evalF() {
    PetscFunctionBegin;
    adouble tR = 0;
    for(int dd = 0;dd<3;dd++) {
      tR += a_sTress[dd];
    }
    I1=tR;
 
    tR /= 3.;
    J2 = 0;
    adouble t;
    for(int dd = 0;dd<3;dd++) {
      t = (a_sTress[dd]-tR);
      J2 += t*t;
    }
   for(int dd = 3;dd<6;dd++) {
      J2 += 2.*a_sTress[dd]*a_sTress[dd]; // s:s off diagonal terms need to go twice
    }
    J2=0.5*J2;
 
    PetscFunctionReturn(0);
  }
 
  /** \brief Evaluate yield function
 
  \f[
  y = 6J_2 + 2I_1\left(\overline{\alpha_1}-\overline{\alpha_0}\right) - 2\overline{\alpha_0} \,\overline{\alpha_1}
  \f]
  where
 
  \f[
  \overline{\alpha_0}=\frac{\partial \psi}{\partial \alpha_0}=\sigma_{yt} + H_t \alpha_0
  \f]
 
  \f[
  \overline{\alpha_1}=\frac{\partial \psi}{\partial \alpha_1}=\sigma_{yc} + H_c \alpha_1
  \f]
 
 
  */
  virtual PetscErrorCode yieldFunction() {
    PetscFunctionBegin;
    ierr = evalF(); CHKERRQ(ierr);
    a_y = 6*J2 + 2*I1*(a_internalFluxes[1]-a_internalFluxes[0]) - 2*a_internalFluxes[0]*a_internalFluxes[1];
    PetscFunctionReturn(0);
  }
 
  /** \brief Flow potential
 
  \f[
  \Psi = 6J_2 + 2\alpha I_1 \left(\overline{\alpha_1}-\overline{\alpha_0}\right) - 2\overline{\alpha_0} \,\overline{\alpha_1}
  \f]
  \f[
  \alpha= \frac{1-2\nu_p}{1+\nu_p}
  \f]
 
 
  */
  virtual PetscErrorCode flowPotential() {
    PetscFunctionBegin;
    ierr = evalF(); CHKERRQ(ierr);
    double alpha= (1-2*nup)/(1+nup); //relation between alpha and plastic Poission's ratio
    a_h =6*J2 + 2*alpha*I1*(a_internalFluxes[1]-a_internalFluxes[0]) - 2*a_internalFluxes[0]*a_internalFluxes[1];
    PetscFunctionReturn(0);
  }
};
 
 
#endif //__SMALL_STRAIN_PLASTICITY_MATERIAL_MODELS_HPP