J2 (Von Misses) plasticity. More...

#include <users_modules/adolc-plasticity/src/ADOLCPlasticityMaterialModels.hpp>

Inheritance diagram for ADOLCPlasticity::J2Plasticity< 3 >:

Collaboration diagram for ADOLCPlasticity::J2Plasticity< 3 >:

Public Types
enum	ModelParams { LAMBDA, MU, ISOH, KINK, PHI, SIGMA_Y, LAST_PARAM }

using	ParamsArray = std::array< double, LAST_PARAM >

Public Types inherited from ADOLCPlasticity::ClosestPointProjection
enum	TypesTags { W = 0, Y, H, LAST_TAPE }

Public Member Functions
	J2Plasticity ()

MoFEMErrorCode	getDefaultMaterialParameters ()

MoFEMErrorCode	addMatBlockOps (MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string block_name, Sev sev)
	Get material parameters from block. More...

MoFEMErrorCode	setParams (short tag, bool &recalculate_elastic_tangent)
	Set material parameters at integration point. More...

MoFEMErrorCode	freeHemholtzEnergy ()
	[Free energy J2] More...

MoFEMErrorCode	evalF ()
	[Yield function J2] More...

MoFEMErrorCode	yieldFunction ()
	Set yield function. More...

MoFEMErrorCode	flowPotential ()
	[Yield function J2] More...

Public Member Functions inherited from ADOLCPlasticity::ClosestPointProjection
VectorAdaptor	getPlasticStrain ()

VectorAdaptor	getTotalStrain ()

VectorAdaptor	getInternalVariables ()

VectorAdaptor	getActiveVariablesYH ()

VectorAdaptor	getStress ()

VectorAdaptor	getInternalFluxes ()

	ClosestPointProjection ()

MoFEMErrorCode	recordW ()
	Record strain energy. More...

MoFEMErrorCode	recordY ()
	Record yield function. More...

MoFEMErrorCode	recordH ()
	Record flow potential. More...

MoFEMErrorCode	playW ()

MoFEMErrorCode	playW_NoHessian ()

MoFEMErrorCode	playY ()

MoFEMErrorCode	playY_NoGradient ()

MoFEMErrorCode	playH ()

MoFEMErrorCode	playH_NoHessian ()

MoFEMErrorCode	createMatAVecR ()

MoFEMErrorCode	evaluatePotentials ()

MoFEMErrorCode	playPotentials ()

MoFEMErrorCode	playPotentials_NoHessian ()

MoFEMErrorCode	calculateR (Vec R)

MoFEMErrorCode	calculateA ()

MoFEMErrorCode	snesCreate ()
	Create nested snes. More...

MoFEMErrorCode	solveClosetProjection ()
	Solve nonlinear system of equations to find stress, internal fluxes, and Lagrange plastic multiplier. More...

MoFEMErrorCode	consistentTangent ()
	Calculate consistent tangent matrix. More...

MoFEMErrorCode	recordTapes ()
	Record tapes. More...

Public Attributes
FTensor::Index< 'i', 3 >	i

FTensor::Index< 'j', 3 >	j

FTensor::Index< 'Z', 6 >	Z

FTensor::Tensor2_symmetric< adouble, 3 >	tTotalStrain

FTensor::Tensor2_symmetric< adouble, 3 >	tPlasticStrain

FTensor::Tensor2_symmetric< adouble, 3 >	tElasticStrain

FTensor::Tensor2_symmetric< adouble, 3 >	tInternalTensor

FTensor::Tensor2_symmetric< adouble, 3 >	tStress

FTensor::Tensor2_symmetric< adouble, 3 >	tBackStress

FTensor::Tensor2_symmetric< adouble, 3 >	tDeviator

adouble	tR

double	mu

double	lambda

double	H

double	K

double	phi

double	sigmaY

boost::shared_ptr< ParamsArray >	defaultParamsArrayPtr = nullptr

boost::shared_ptr< ParamsArray >	paramsArrayPtr = nullptr

boost::shared_ptr< ParamsArray >	oldParamsArrayPtr = nullptr

adouble	f
	[Free energy J2] More...

Public Attributes inherited from ADOLCPlasticity::ClosestPointProjection
int	nbInternalVariables

boost::function< int(int, int, int)>	integrationRule

VectorDouble	internalVariables0

VectorDouble	plasticStrain0

std::array< int, LAST_TAPE >	tapesTags

VectorAdaptor	activeVariablesW

VectorAdaptor	gradientW

double	w

double	y

double	h

double	deltaGamma

MatrixDouble	Ep

MatrixDouble	Cp

MatrixDouble	Cep

ublas::symmetric_matrix< double, ublas::lower >	C

ublas::symmetric_matrix< double, ublas::lower >	D

MatrixDouble	hessianW

VectorDouble	gradientY

VectorDouble	gradientH

MatrixDouble	hessianH

MatrixDouble	partialWStrainPlasticStrain

VectorAdaptor	partialYSigma

VectorAdaptor	partialYFlux

VectorAdaptor	partialHSigma

VectorAdaptor	partialHFlux

ublas::symmetric_matrix< double, ublas::lower >	partial2HSigma

ublas::symmetric_matrix< double, ublas::lower >	partial2HFlux

MatrixDouble	partial2HSigmaFlux

SmartPetscObj< Mat >	A

SmartPetscObj< Vec >	R

SmartPetscObj< Vec >	Chi

SmartPetscObj< Vec >	dChi

ublas::matrix< double, ublas::column_major >	dataA

SmartPetscObj< SNES >	sNes

ublas::vector< adouble >	a_sTrain

ublas::vector< adouble >	a_plasticStrain

ublas::vector< adouble >	a_internalVariables

ublas::vector< adouble >	a_sTress

ublas::vector< adouble >	a_internalFluxes

adouble	a_w

adouble	a_y

adouble	a_h

Detailed Description

J2 (Von Misses) plasticity.

Free energy:

\[ \psi(\varepsilon^e, \alpha, \beta) = \frac{1}{2} \lambda \textrm{tr}[\varepsilon^e]^2 + \mu \varepsilon^e : \varepsilon^e + \sigma_y\alpha + \frac{1}{2} \phi H \alpha^2 + \frac{1}{2} (1-\phi)K \beta^2 \]

Isotropic hardening variable \(\alpha\) is in first index of internal variables vector. Kinematic hardening variable is in the remaining indices of internal variables vector.

Yield function:

Deviator of actual stress:

\[ \eta=\textrm{dev}[\sigma]-\overline{\beta} \]

where \(\overline{\beta}\) is back stress.

Square of the deviator norm

\[ f = \frac{3}{2} \eta:\eta \]

Yield function:

\[ y = \sqrt{f} - \overline{\alpha} \]

where \(f\) is defined in eva

This is associated potential model, so flow potential is equal to yield

Definition at line 73 of file ADOLCPlasticityMaterialModels.hpp.

Member Typedef Documentation

◆ ParamsArray

using ADOLCPlasticity::J2Plasticity< 3 >::ParamsArray = std::array<double, LAST_PARAM>

Definition at line 103 of file ADOLCPlasticityMaterialModels.hpp.

Member Enumeration Documentation

◆ ModelParams

enum ADOLCPlasticity::J2Plasticity< 3 >::ModelParams

Enumerator
LAMBDA
MU
ISOH
KINK
PHI
SIGMA_Y
LAST_PARAM

Definition at line 101 of file ADOLCPlasticityMaterialModels.hpp.

101 { LAMBDA, MU, ISOH, KINK, PHI, SIGMA_Y, LAST_PARAM };

Constructor & Destructor Documentation

◆ J2Plasticity()

ADOLCPlasticity::J2Plasticity< 3 >::J2Plasticity ( )

inline

Examples: ADOLCPlasticityMaterialModels.hpp.

Definition at line 75 of file ADOLCPlasticityMaterialModels.hpp.

                  : ClosestPointProjection() {
     nbInternalVariables = 7;
     activeVariablesW.resize(12 + nbInternalVariables, false);
   }

Member Function Documentation

◆ addMatBlockOps()

MoFEMErrorCode ADOLCPlasticity::J2Plasticity< 3 >::addMatBlockOps	(	MoFEM::Interface &	m_field,
		boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &	pip,
		std::string	block_name,
		Sev	sev
	)

inlinevirtual

Get material parameters from block.

Parameters

m_field
pip
block_name
sev

Returns: MoFEMErrorCode

Reimplemented from ADOLCPlasticity::ClosestPointProjection.

Definition at line 175 of file ADOLCPlasticityMaterialModels.hpp.

                                                 {
     MoFEMFunctionBegin;
  
     struct OpMatBlocks : public ForcesAndSourcesCore::UserDataOperator {
  
       using OP = ForcesAndSourcesCore::UserDataOperator;
  
       OpMatBlocks(boost::shared_ptr<ParamsArray> params_array_ptr,
                   boost::shared_ptr<ParamsArray> def_params_array_otr,
                   MoFEM::Interface &m_field, Sev sev,
                   std::vector<const CubitMeshSets *> meshset_vec_ptr)
           : OP(NOSPACE, OP::OPSPACE), paramsArrayPtr(params_array_ptr),
             defParamsArray(def_params_array_otr) {
         CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
                           "Can not get data from block");
       }
  
       MoFEMErrorCode doWork(int side, EntityType type,
                             EntitiesFieldData::EntData &data) {
         MoFEMFunctionBegin;
  
         for (auto &b : blockData) {
  
           if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
             std::copy(b.paramsArray.begin(), b.paramsArray.end(),
                       paramsArrayPtr->begin());
             MoFEMFunctionReturnHot(0);
           }
         }
  
         std::copy(defParamsArray->begin(), defParamsArray->end(),
                   paramsArrayPtr->begin());
  
         MoFEMFunctionReturn(0);
       }
  
     private:
       boost::shared_ptr<ParamsArray> paramsArrayPtr;
       boost::shared_ptr<ParamsArray> defParamsArray;
  
       struct BlockData {
         ParamsArray paramsArray;
         Range blockEnts;
       };
  
       std::vector<BlockData> blockData;
  
       MoFEMErrorCode
       extractBlockData(MoFEM::Interface &m_field,
                        std::vector<const CubitMeshSets *> meshset_vec_ptr,
                        Sev sev) {
         MoFEMFunctionBegin;
  
         for (auto m : meshset_vec_ptr) {
           MOFEM_TAG_AND_LOG("ADOLCPlasticityWold", sev, "JP2 parameters") << *m;
           std::vector<double> block_data;
           CHKERR m->getAttributes(block_data);
           if (block_data.size() != 6) {
             SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                     "Expected that block has two attribute");
           }
           auto get_block_ents = [&]() {
             Range ents;
             CHKERR
             m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
             return ents;
           };
  
           blockData.push_back({*defParamsArray, get_block_ents()});
           std::copy(block_data.begin(), block_data.end(),
                     blockData.back().paramsArray.begin());
           const auto young_modulus = blockData.back().paramsArray[LAMBDA];
           const auto poisson_ratio = blockData.back().paramsArray[MU];
  
           // Is assumed that uset provide young_modulus and poisson_ratio in
           // fist two argiumens of the block
  
           MOFEM_LOG("ADOLCPlasticityWold", sev)
               << "Young modulus: " << (blockData.back().paramsArray)[LAMBDA];
           MOFEM_LOG("ADOLCPlasticityWold", sev)
               << "Poisson ratio: " << (blockData.back().paramsArray)[MU];
  
           blockData.back().paramsArray[LAMBDA] =
               LAMBDA(young_modulus, poisson_ratio);
           blockData.back().paramsArray[MU] = MU(young_modulus, poisson_ratio);
  
           MOFEM_LOG("ADOLCPlasticityWold", sev)
               << "LAMBDA: " << (blockData.back().paramsArray)[LAMBDA];
           MOFEM_LOG("ADOLCPlasticityWold", sev)
               << "MU: " << (blockData.back().paramsArray)[MU];
           MOFEM_LOG("ADOLCPlasticityWold", sev)
               << "H: " << (blockData.back().paramsArray)[ISOH];
           MOFEM_LOG("ADOLCPlasticityWold", sev)
               << "K: " << (blockData.back().paramsArray)[KINK];
           MOFEM_LOG("ADOLCPlasticityWold", sev)
               << "PHI: " << (blockData.back().paramsArray)[PHI];
           MOFEM_LOG("ADOLCPlasticityWold", sev)
               << "SIGMA_Y: " << (blockData.back().paramsArray)[SIGMA_Y];
         }
         MoFEMFunctionReturn(0);
       }
     };
  
     auto cubit_meshsets_vec_ptr =
         m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
             std::regex((boost::format("%s(.*)") % block_name).str()));
  
     pip.push_back(new OpMatBlocks(paramsArrayPtr, defaultParamsArrayPtr,
                                   m_field, sev, cubit_meshsets_vec_ptr));
  
     MoFEMFunctionReturn(0);
   }

◆ evalF()

MoFEMErrorCode ADOLCPlasticity::J2Plasticity< 3 >::evalF ( )

inline

[Yield function J2]

Definition at line 373 of file ADOLCPlasticityMaterialModels.hpp.

                          {
     MoFEMFunctionBegin;
     auto t_voight_stress_op = voight_to_stress_op();
     auto t_vioght_stress = getFTensor1FromPtr<6>(&a_sTress[0]);
     auto t_back_stress = getFTensor1FromPtr<6>(&a_internalFluxes[1]);
     tStress(i, j) = t_voight_stress_op(i, j, Z) * t_vioght_stress(Z);
     tBackStress(i, j) = t_voight_stress_op(i, j, Z) * t_back_stress(Z);
  
     constexpr double third = boost::math::constants::third<double>();
     constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<double>();
     tR = tStress(i, i) / 3.;
     tDeviator(i, j) = tStress(i, j) - tR * t_kd(i, j) - tBackStress(i, j);
  
     // Fix the constant to mach uniaxial test
     constexpr double c = 3. / 2.;
     f = c * tDeviator(i, j) * tDeviator(i, j);
  
     MoFEMFunctionReturn(0);
   }

◆ flowPotential()

MoFEMErrorCode ADOLCPlasticity::J2Plasticity< 3 >::flowPotential ( )

inlinevirtual

[Yield function J2]

Implements ADOLCPlasticity::ClosestPointProjection.

Definition at line 401 of file ADOLCPlasticityMaterialModels.hpp.

                                  {
     MoFEMFunctionBegin;
     CHKERR evalF();
     a_h = sqrt(f) - a_internalFluxes[0];
     MoFEMFunctionReturn(0);
   }

◆ freeHemholtzEnergy()

MoFEMErrorCode ADOLCPlasticity::J2Plasticity< 3 >::freeHemholtzEnergy ( )

inlinevirtual

[Free energy J2]

Implements ADOLCPlasticity::ClosestPointProjection.

Reimplemented in ADOLCPlasticity::J2Plasticity< 2 >.

Definition at line 324 of file ADOLCPlasticityMaterialModels.hpp.

                                       {
     MoFEMFunctionBegin;
  
     // ADOL-C register parameters (those can varied for integration points)
     auto p_lambda = mkparam(lambda);  // 0
     auto p_mu = mkparam(mu);          // 1
     auto p_H = mkparam(H);            // 2
     auto p_K = mkparam(K);            // 3
     auto p_phi = mkparam(phi);        // 4
     auto p_sigma_y = mkparam(sigmaY); // 5
  
     // Operator converting Voight and tensor notation
     auto t_voight_strain_op = voight_to_strain_op();
  
     // Variable in voight notation
     auto t_voight_total_strain = getFTensor1FromPtr<6>(&a_sTrain[0]);
     auto t_voight_plastic_strain = getFTensor1FromPtr<6>(&a_plasticStrain[0]);
     auto t_voight_internal_tensor =
         getFTensor1FromPtr<6>(&a_internalVariables[1]);
  
     // Convert variables to tensor notation
     tTotalStrain(i, j) = t_voight_strain_op(i, j, Z) * t_voight_total_strain(Z);
     tPlasticStrain(i, j) =
         t_voight_strain_op(i, j, Z) * t_voight_plastic_strain(Z);
     tInternalTensor(i, j) =
         t_voight_strain_op(i, j, Z) * t_voight_internal_tensor(Z);
  
     // Evaluate elastic strain
     tElasticStrain(i, j) = tTotalStrain(i, j) - tPlasticStrain(i, j);
     tR = tElasticStrain(i, i);
  
     // Strain energy
     a_w = 0.5 * p_lambda * tR * tR +
           p_mu * (tElasticStrain(i, j) * tElasticStrain(i, j));
     // Isotropic strain energy
     a_w +=
         a_internalVariables[0] * p_sigma_y +
         0.5 * p_phi * p_H * (a_internalVariables[0] * a_internalVariables[0]);
     // Kinematic strain energy
     a_w += (0.5 * (1 - p_phi) * p_K) *
            (tInternalTensor(i, j) * tInternalTensor(i, j));
  
     MoFEMFunctionReturn(0);
   }

◆ getDefaultMaterialParameters()

MoFEMErrorCode ADOLCPlasticity::J2Plasticity< 3 >::getDefaultMaterialParameters ( )

inline

Definition at line 108 of file ADOLCPlasticityMaterialModels.hpp.

                                                 {
     MoFEMFunctionBegin;
  
     double young_modulus = 1;
     double poisson_ratio = 0.25;
     sigmaY = 1;
     H = 0.1;
     K = 0;
     phi = 1;
  
     CHKERR PetscOptionsGetScalar(PETSC_NULL, "-young_modulus", &young_modulus,
                                  PETSC_NULL);
     CHKERR PetscOptionsGetScalar(PETSC_NULL, "-poisson_ratio", &poisson_ratio,
                                  0);
     CHKERR PetscOptionsGetScalar(PETSC_NULL, "-sigma_y", &sigmaY, PETSC_NULL);
     CHKERR PetscOptionsGetScalar(PETSC_NULL, "-H", &H, PETSC_NULL);
     CHKERR PetscOptionsGetScalar(PETSC_NULL, "-K", &K, PETSC_NULL);
     CHKERR PetscOptionsGetScalar(PETSC_NULL, "-phi", &phi, PETSC_NULL);
  
     MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
         << "Young modulus: " << young_modulus;
     MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
         << "Poisson ratio: " << poisson_ratio;
  
     mu = MU(young_modulus, poisson_ratio);
     lambda = LAMBDA(young_modulus, poisson_ratio);
  
     defaultParamsArrayPtr = boost::make_shared<ParamsArray>();
     (*defaultParamsArrayPtr)[LAMBDA] = lambda;
     (*defaultParamsArrayPtr)[MU] = mu;
     (*defaultParamsArrayPtr)[H] = H;
     (*defaultParamsArrayPtr)[K] = K;
     (*defaultParamsArrayPtr)[PHI] = phi;
     (*defaultParamsArrayPtr)[SIGMA_Y] = sigmaY;
  
     MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
         << "LAMBDA: " << (*defaultParamsArrayPtr)[LAMBDA];
     MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
         << "MU: " << (*defaultParamsArrayPtr)[MU];
     MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
         << "H: " << (*defaultParamsArrayPtr)[ISOH];
     MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
         << "K: " << (*defaultParamsArrayPtr)[KINK];
     MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
         << "PHI: " << (*defaultParamsArrayPtr)[PHI];
     MOFEM_LOG("ADOLCPlasticityWold", Sev::inform)
         << "SIGMA_Y: " << (*defaultParamsArrayPtr)[SIGMA_Y];
  
     paramsArrayPtr = boost::make_shared<ParamsArray>();
     std::copy(defaultParamsArrayPtr->begin(), defaultParamsArrayPtr->end(),
               paramsArrayPtr->begin());
     oldParamsArrayPtr = boost::make_shared<ParamsArray>();
     std::fill(oldParamsArrayPtr->begin(), oldParamsArrayPtr->end(), 0);
  
     MoFEMFunctionReturn(0);
   };

◆ setParams()

MoFEMErrorCode ADOLCPlasticity::J2Plasticity< 3 >::setParams	(	short	tag,
		bool &	recalculate_elastic_tangent
	)

inlinevirtual

Set material parameters at integration point.

Parameters

tag
recalculate_elastic_tangent

Returns: MoFEMErrorCode

Reimplemented from ADOLCPlasticity::ClosestPointProjection.

Definition at line 297 of file ADOLCPlasticityMaterialModels.hpp.

                                                                          {
     MoFEMFunctionBegin;
     switch (tag) {
     case TypesTags::W: {
       set_param_vec(tapesTags[tag], paramsArrayPtr->size(),
                     paramsArrayPtr->data());
       for (auto p = 0; p != LAST_PARAM; p++) {
         if ((*paramsArrayPtr)[p] != (*oldParamsArrayPtr)[p]) {
           recalculate_elastic_tangent = true;
           std::copy(paramsArrayPtr->begin(), paramsArrayPtr->end(),
                     oldParamsArrayPtr->begin());
           break;
         }
       }
  
     } break;
     case TypesTags::Y:
     case TypesTags::H:
       break;
     default:
       SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
               "Unknown tag for setParams");
     }
     MoFEMFunctionReturn(0);
   }

◆ yieldFunction()

MoFEMErrorCode ADOLCPlasticity::J2Plasticity< 3 >::yieldFunction ( )

inlinevirtual

Set yield function.

Implements ADOLCPlasticity::ClosestPointProjection.

Definition at line 393 of file ADOLCPlasticityMaterialModels.hpp.

                                  {
     MoFEMFunctionBegin;
     CHKERR evalF();
     a_y = sqrt(f) - a_internalFluxes[0];
     MoFEMFunctionReturn(0);
   }