Commit 4b0b9cd5 authored by Laurent STAINIER's avatar Laurent STAINIER

Add themo-plasticity model with mixed hardening + several fixes

parent 97b616c1
......@@ -4,7 +4,7 @@
* This file is part of ZorgLib, a computational simulation framework
* for thermomechanics of solids and structures (systems in general).
*
* Copyright (c) 2001-2019, L. Stainier.
* Copyright (c) 2001-2020, L. Stainier.
* See file LICENSE.txt for license information.
* Please report all bugs and problems to <Laurent.Stainier@ec-nantes.fr>.
*/
......@@ -26,7 +26,7 @@ USING_MATLIB_NAMESPACE
// check consistency of properties
void ElastoPlasticBar::checkProperties(MaterialProperties& material,std::ostream *os) {
if (os) (*os) << "\nElasto-plasticity model with isotropic hardening (pure 1D bar):" << std::endl;
// density
try {
z_real_t rho = material.getRealProperty("MASS_DENSITY");
......@@ -39,7 +39,7 @@ void ElastoPlasticBar::checkProperties(MaterialProperties& material,std::ostream
catch (NoSuchPropertyException) {
if (os) (*os) << "\n\tmass density is not defined" << std::endl;
}
// bar cross section
try {
z_real_t A = material.getRealProperty("SECTION");
......@@ -56,10 +56,10 @@ void ElastoPlasticBar::checkProperties(MaterialProperties& material,std::ostream
// elastic potential
this->potential->checkProperties(material,os);
// dilatancy model
if (this->dilatancy) this->dilatancy->checkProperties(material,os);
// viscoplastic model
viscoPlasticity->checkProperties(material,os);
}
......@@ -78,10 +78,10 @@ z_real_t ElastoPlasticBar::incrementalPotential(const MaterialProperties& materi
bool update,bool tangent) {
// get cross-section
z_real_t A = material.getRealProperty("SECTION");
// get strain
z_real_t eps = state.grad[0];
// get plastic strain
z_real_t epsP0 = state0.internal[1];
......@@ -113,12 +113,12 @@ z_real_t ElastoPlasticBar::plasticUpdate(const MaterialProperties& material,
// get equivalent plastic strain
z_real_t ePl0 = intV0[0];
z_real_t ePl = intV[0];
// extract internal parameters
unsigned int nIntPar = intV.size()-2;
const MatLibArray intPar0(intV0,nIntPar,2);
MatLibArray intPar(intV,nIntPar,2);
// update
viscoPlasticity->initialize = true;
viscoPlasticity->finalize = false;
......@@ -128,17 +128,17 @@ z_real_t ElastoPlasticBar::plasticUpdate(const MaterialProperties& material,
// perform update (radial return)
radialReturn(material,extPar,intPar0,intPar,
eps,epsPl0,ePl0,ePl,dTime);
// update internal variables
intV[0] = ePl;
// update plastic strain
dEPl = ePl-ePl0;
epsPl = epsPl0+std::copysign(dEPl,eps-epsPl0);
viscoPlasticity->finalize = true;
}
// elastic free energy
z_real_t epsEl = eps-epsPl;
z_real_t We = this->storedEnergy(material,extPar,epsEl,sig,M,
......@@ -150,14 +150,14 @@ z_real_t ElastoPlasticBar::plasticUpdate(const MaterialProperties& material,
z_real_t dummy,Hp;
z_real_t Wp = viscoPlasticity->irreversibleEnergy(material,extPar,intPar0,intPar,ePl0,ePl,
dummy,Hp,dTime,false,computeTangent);
// tangents
dEPl = ePl-ePl0;
if (computeTangent && dEPl > 0.0) {
// (visco)plastic correction
M -= M*M/(M+Hp);
}
return We+Wp-intV0[1];
}
......@@ -169,17 +169,17 @@ unsigned int ElastoPlasticBar::radialReturn(const MaterialProperties& material,
z_real_t ePl0,z_real_t& ePl,z_real_t dTime) {
static const unsigned int ITMAX = 30;
static const z_real_t MULT = 0.9;
static const z_real_t PREC = 1.0e-14;
static const z_real_t PREC = 1.0e-12;
static const z_real_t TOLE = 1.0e-07;
static const z_real_t THRSHLD = 0.1*std::numeric_limits<z_real_t>::max();
// get algorithmic parameter
unsigned int maxIt;
if (material.checkProperty("RR_MAX_ITER_PARAMETER"))
maxIt = material.getIntegerProperty("RR_MAX_ITER_PARAMETER");
else
maxIt = ITMAX;
// compute elastic predictor
z_real_t epsEl = eps-epsPl;
z_real_t sig,M;
......@@ -192,7 +192,7 @@ unsigned int ElastoPlasticBar::radialReturn(const MaterialProperties& material,
sigPl,Hp,dTime,true,true);
z_real_t fct0 = std::fabs(sig)-sigPl;
if (fct0 <= 0.0) return 0;
// apply plastic corrector
z_real_t dEPl = 0.0;
z_real_t ePl00 = ePl0;
......@@ -213,9 +213,9 @@ unsigned int ElastoPlasticBar::radialReturn(const MaterialProperties& material,
dEPl = mult*(ePl-ePl00);
}
if (std::fabs(dEPl) < PREC) break;
epsEl -= std::copysign(dEPl,sig);
this->storedEnergy(material,extPar,epsEl,sig,M,true,true);
epsEl -= std::copysign(1.0,sig)*dEPl;
ePl += dEPl;
this->storedEnergy(material,extPar,epsEl,sig,M,true,true);
viscoPlasticity->irreversibleEnergy(material,extPar,intPar0,intPar,ePl0,ePl,
sigPl,Hp,dTime,true,true);
fct = std::fabs(sig)-sigPl;
......@@ -229,7 +229,7 @@ unsigned int ElastoPlasticBar::radialReturn(const MaterialProperties& material,
if (iter == maxIt) {
throw UpdateFailedException("no convergence in radial return");
}
return iter;
}
......@@ -240,7 +240,7 @@ unsigned int ElastoPlasticBar::radialReturn(const MaterialProperties& material,
// check consistency of properties
void ElastoPlasticBarMixed::checkProperties(MaterialProperties& material,std::ostream *os) {
if (os) (*os) << "\nElasto-plasticity model with mixed (isotropic+linear kinematic) hardening (pure 1D bar):" << std::endl;
// density
try {
z_real_t rho = material.getRealProperty("MASS_DENSITY");
......@@ -253,7 +253,7 @@ void ElastoPlasticBarMixed::checkProperties(MaterialProperties& material,std::os
catch (NoSuchPropertyException) {
if (os) (*os) << "\n\tmass density is not defined" << std::endl;
}
// bar cross section
try {
z_real_t A = material.getRealProperty("SECTION");
......@@ -267,13 +267,13 @@ void ElastoPlasticBarMixed::checkProperties(MaterialProperties& material,std::os
if (os) (*os) << "\n\tbar cross-section set to 1" << std::endl;
material.setProperty("SECTION",1.0);
}
// elastic potential
this->potential->checkProperties(material,os);
// dilatancy model
if (this->dilatancy) this->dilatancy->checkProperties(material,os);
// kinematic hardening modulus
z_real_t Hk;
try {
......@@ -303,39 +303,39 @@ z_real_t ElastoPlasticBarMixed::plasticUpdate(const MaterialProperties& material
// get equivalent plastic strain
z_real_t ePl0 = intV0[0];
z_real_t ePl = intV[0];
// extract internal parameters
unsigned int nIntPar = intV.size()-2;
const MatLibArray intPar0(intV0,nIntPar,2);
MatLibArray intPar(intV,nIntPar,2);
// update
viscoPlasticity->initialize = true;
viscoPlasticity->finalize = false;
z_real_t dEPl=0.0;
if (update) {
// perform update (radial return)
radialReturn(material,extPar,intPar0,intPar,
eps,epsPl0,ePl0,ePl,dTime);
// update internal variables
intV[0] = ePl;
// update plastic strain
dEPl = ePl-ePl0;
epsPl = epsPl0+std::copysign(dEPl,eps-epsPl0);
viscoPlasticity->finalize = true;
}
// elastic free energy
z_real_t epsEl = eps-epsPl;
z_real_t We = this->storedEnergy(material,extPar,epsEl,sig,M,
update || computeTangent,
computeTangent);
if (update) intV[1] = We;
// plastic free energy increment + dissipated energy
z_real_t dummy,Hp;
z_real_t Wp = viscoPlasticity->irreversibleEnergy(material,extPar,intPar0,intPar,ePl0,ePl,
......@@ -352,7 +352,7 @@ z_real_t ElastoPlasticBarMixed::plasticUpdate(const MaterialProperties& material
// (visco)plastic correction
M -= M*M/(M+Hp+Hk);
}
return We-intV0[1]+WpKin-0.5*Hk*epsPl0*epsPl0+Wp;
}
......@@ -364,17 +364,17 @@ unsigned int ElastoPlasticBarMixed::radialReturn(const MaterialProperties& mater
z_real_t ePl0,z_real_t& ePl,z_real_t dTime) {
static const unsigned int ITMAX = 30;
static const z_real_t MULT = 0.9;
static const z_real_t PREC = 1.0e-14;
static const z_real_t PREC = 1.0e-12;
static const z_real_t TOLE = 1.0e-07;
static const z_real_t THRSHLD = 0.1*std::numeric_limits<z_real_t>::max();
// get algorithmic parameter
unsigned int maxIt;
if (material.checkProperty("RR_MAX_ITER_PARAMETER"))
maxIt = material.getIntegerProperty("RR_MAX_ITER_PARAMETER");
else
maxIt = ITMAX;
// compute elastic predictor
z_real_t epsEl = eps-epsPl;
z_real_t sig,M;
......@@ -390,7 +390,7 @@ unsigned int ElastoPlasticBarMixed::radialReturn(const MaterialProperties& mater
sigPl,Hp,dTime,true,true);
z_real_t fct0 = std::fabs(sig-Hk*epsPl)-sigPl;
if (fct0 <= 0.0) return 0;
// apply plastic corrector
z_real_t dEPl = 0.0;
z_real_t ePl00 = ePl0;
......@@ -411,9 +411,9 @@ unsigned int ElastoPlasticBarMixed::radialReturn(const MaterialProperties& mater
dEPl = mult*(ePl-ePl00);
}
if (std::fabs(dEPl) < PREC) break;
epsEl -= std::copysign(dEPl,sig);
epsEl -= std::copysign(1.0,sig)*dEPl;
this->storedEnergy(material,extPar,epsEl,sig,M,true,true);
epsPl += std::copysign(dEPl,sig);
epsPl += std::copysign(1.0,sig)*dEPl;
ePl += dEPl;
viscoPlasticity->irreversibleEnergy(material,extPar,intPar0,intPar,ePl0,ePl,
sigPl,Hp,dTime,true,true);
......@@ -428,7 +428,7 @@ unsigned int ElastoPlasticBarMixed::radialReturn(const MaterialProperties& mater
if (iter == maxIt) {
throw UpdateFailedException("no convergence in radial return");
}
return iter;
}
......
This diff is collapsed.
/*
* $Id$
*
* This file is part of ZorgLib, a computational simulation framework
* for thermomechanics of solids and structures (systems in general).
*
* Copyright (c) 2001-2020, L. Stainier.
* See file LICENSE.txt for license information.
* Please report all bugs and problems to <Laurent.Stainier@ec-nantes.fr>.
*/
#include "J2ThermoPlasticityMixed.h"
#ifdef MATLIB_USE_NAMESPACE
USING_MATLIB_NAMESPACE
#endif
/*
* Methods for class LinearMixedJ2ThermoPlasticityBuilder.
*/
// the instance
LinearMixedJ2ThermoPlasticityBuilder const* LinearMixedJ2ThermoPlasticityBuilder::BUILDER
= new LinearMixedJ2ThermoPlasticityBuilder();
// constructor
LinearMixedJ2ThermoPlasticityBuilder::LinearMixedJ2ThermoPlasticityBuilder() {
ModelDictionary::add("LINEAR_MIXED_J2_THERMO_PLASTICITY",*this);
}
// build model
ConstitutiveModel* LinearMixedJ2ThermoPlasticityBuilder::build(unsigned int d) const {
switch(d) {
case 3:
return new LinearMixedJ2ThermoPlasticity3D();
break;
case 2:
return new LinearMixedJ2ThermoPlasticity2D();
break;
case 1:
return new LinearMixedJ2ThermoPlasticity1D();
break;
default:
return 0;
break;
}
}
/*
* Methods for class NonLinearMixedJ2ThermoPlasticityBuilder.
*/
// the instance
NonLinearMixedJ2ThermoPlasticityBuilder const* NonLinearMixedJ2ThermoPlasticityBuilder::BUILDER
= new NonLinearMixedJ2ThermoPlasticityBuilder();
// constructor
NonLinearMixedJ2ThermoPlasticityBuilder::NonLinearMixedJ2ThermoPlasticityBuilder() {
ModelDictionary::add("NONLINEAR_MIXED_J2_THERMO_PLASTICITY",*this);
}
// build model
ConstitutiveModel* NonLinearMixedJ2ThermoPlasticityBuilder::build(unsigned int d) const {
switch(d) {
case 3:
return new NonLinearMixedJ2ThermoPlasticity3D();
break;
case 2:
return new NonLinearMixedJ2ThermoPlasticity2D();
break;
case 1:
return new NonLinearMixedJ2ThermoPlasticity1D();
break;
default:
return 0;
break;
}
}
This diff is collapsed.
......@@ -238,6 +238,7 @@ class J2ThermoPlasticitySimple : virtual public ThermoElastoPlasticity<ALG> {
bool update,bool computeTangent) {
static const z_real_t ONE_THIRD = 1./3.;
static const SYM_TENSOR I = SYM_TENSOR::identity();
z_real_t N;
SYM_TENSOR epsEl,sigDev,Mp;
......@@ -265,7 +266,6 @@ class J2ThermoPlasticitySimple : virtual public ThermoElastoPlasticity<ALG> {
M,dSig,C,true,false);
// compute stress deviator
static const SYM_TENSOR I = SYM_TENSOR::identity();
z_real_t p = trace(sig);
sigDev = sig-ONE_THIRD*p*I;
......@@ -299,8 +299,8 @@ class J2ThermoPlasticitySimple : virtual public ThermoElastoPlasticity<ALG> {
epsEl = eps-epsPl;
// elastic free energy
z_real_t W = this->storedEnergy(material,extPar,epsEl,Th1,sig,
N,M,dSig,C,
z_real_t W = this->storedEnergy(material,extPar,epsEl,Th1,
sig,N,M,dSig,C,
update || computeTangent,
computeTangent);
if (update) intV[2] = W;
......@@ -377,7 +377,8 @@ class J2ThermoPlasticitySimple : virtual public ThermoElastoPlasticity<ALG> {
z_real_t coef3 = 4*ONE_THIRD*mu*(1.5*coef2-coef1);
M -= ((coef1*mu2)*(II-KK)+coef3*outerProd(Mp,Mp));
z_real_t dEPldT = dSigPl-innerProd2(dSig,Mp);
SYM_TENSOR dSigDev = dSig-(ONE_THIRD*trace(dSig))*I;
z_real_t dEPldT = dSigPl-innerProd2(dSigDev,Mp);
dSig += coef2*dEPldT*Mp;
C -= dEPldT*dEPldT/(mu3+Hp);
......
......@@ -53,7 +53,7 @@ add_test( NAME ConsistencyTest-nonlinear_isotropic_elasticity WORKING_DIRECTORY
add_test( NAME ConsistencyTest-isotropic_hyperelasticity WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
COMMAND ${EXECUTABLE_OUTPUT_PATH}/MatlConsistencyTest test-isotropic-hyperelasticity.dat )
add_test( NAME ConsistencyTest-neohookean_hyperelasticity WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
add_test( NAME ConsistencyTest-neohookean_hyperelasticity WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
COMMAND ${EXECUTABLE_OUTPUT_PATH}/MatlConsistencyTest test-neohookean.dat )
add_test( NAME ConsistencyTest-compressible_neohookean_hyperelasticity WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
COMMAND ${EXECUTABLE_OUTPUT_PATH}/MatlConsistencyTest test-compressible-neohookean.dat )
......@@ -110,6 +110,11 @@ add_test( NAME ConsistencyTest-linear_isotropic_hardening_J2_finite_plasticity
add_test( NAME ConsistencyTest-nonlinear_isotropic_hardening_J2_finite_plasticity WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
COMMAND ${EXECUTABLE_OUTPUT_PATH}/MatlConsistencyTest test-nonlinear-isotropic-j2-finite-plasticity.dat )
# thermo-elasticity
add_test( NAME ConsistencyTest-isotropic_thermo_elasticity WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
COMMAND ${EXECUTABLE_OUTPUT_PATH}/MatlConsistencyTest test-isotropic-thermo-elasticity.dat )
# thermo-plasticity (small strains)
#add_test( NAME ConsistencyTest-linear_isotropic_hardening_J2_thermo_plasticity WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
......
3
ISOTROPIC_THERMO_ELASTICITY
../materials/mild-steel.mat
1.0 0
0.00 0.10 0.00 0.00 0.00 0.01 10.0
1.0e-7 1.0e-7
......@@ -2,5 +2,5 @@
LINEAR_ISOTROPIC_ELASTO_PLASTIC_BAR
../materials/mild-steel.mat
1.0 0
0.10
1.0e-8 1.0e-5
\ No newline at end of file
0.01
1.0e-7 1.0e-8
......@@ -2,5 +2,5 @@
LINEAR_ISOTROPIC_J2_FINITE_PLASTICITY
../materials/mild-steel.mat
1.0 0
1.00 0.10 0.00 -0.20 1.00 0.00 0.00 0.00 1.10
1.0e-5 1.0e-4
\ No newline at end of file
1.00 0.10 0.00 -0.05 1.00 0.00 0.00 0.00 1.10
1.0e-5 1.0e-4
......@@ -3,4 +3,4 @@ LINEAR_MIXED_ELASTO_PLASTIC_BAR
../materials/mild-steel.mat
1.0 0
0.10
1.0e-6 1.0e-4
\ No newline at end of file
1.0e-6 1.0e-5
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment