v0.5.86
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
*/
internalVariables.resize(7,false);
}
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++) {
}
a_w = 0.5*lambda*tR*tR;
int dd = 0;
for(;dd<3;dd++) {
a_w += mu*t*t;
}
for(;dd<6;dd++) {
a_w += 0.5*mu*t*t;
}
//plastic part
//isotropic
//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$.
*/
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++) {
f += c*t*t;
}
for(int dd = 3;dd<6;dd++) {
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;
a_y = sqrt(f) - a_internalFluxes[0];
PetscFunctionReturn(0);
}
/** \brief Associated flow potential
See \ref yieldFunction.
*/
virtual PetscErrorCode flowPotential() {
PetscFunctionBegin;
a_h = sqrt(f) - a_internalFluxes[0];
PetscFunctionReturn(0);
}
};
/** \brief Small strain plasticity with paraboloidal yield criterion (Isotropic Hardening)
* \ingroup small_strain_plasticity
*/
internalVariables.resize(2,false);
}
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++) {
}
a_w = 0.5*lambda*tR*tR;
int dd = 0;
for(;dd<3;dd++) {
a_w += mu*t*t;
}
for(;dd<6;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;
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;
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;
double alpha= (1-2*nup)/(1+nup); //relation between alpha and plastic Poission's ratio
PetscFunctionReturn(0);
}
};
#endif //__SMALL_STRAIN_PLASTICITY_MATERIAL_MODELS_HPP