Material class for 2D correspondence material model for finite strain: plane strain, generalized plane strain, weak plane stress. More...

#include <ComputePlaneFiniteStrainNOSPD.h>

Inheritance diagram for ComputePlaneFiniteStrainNOSPD:

[legend]

Public Member Functions
	ComputePlaneFiniteStrainNOSPD (const InputParameters &parameters)

virtual void	initQpStatefulProperties () override

Static Public Member Functions
static MooseEnum	decompositionType ()

Protected Member Functions
virtual void	computeQpFhat () override

virtual void	computeQpStrain () override
	Function to compute strain tensors. More...

void	computeQpStrainRotationIncrements (RankTwoTensor &e, RankTwoTensor &r)
	Function to compute strain and rotational increments. More...

void	subtractEigenstrainIncrementFromStrain (RankTwoTensor &strain)
	Function to compute the mechanical strain tensor by subtracting thermal strain from the total strain. More...

virtual void	computeProperties () override

virtual void	computeBondStretch () override

virtual void	computeQpDeformationGradient ()
	Function to compute bond-associated deformation gradient. More...


virtual Real	computeQpOutOfPlaneDeformationGradient ()
	Functions to compute the out-of-plane component of deformation gradient for generalized plane strain and weak plane stress. More...

virtual Real	computeQpOutOfPlaneDeformationGradientOld ()

Protected Attributes
std::vector< RankTwoTensor >	_Fhat
	'Incremental' deformation gradient More...

const bool	_plane_strain
	Plane strain problem or not, this is only used for mechanical stretch calculation. More...


MaterialProperty< RankTwoTensor > &	_strain_rate
	Material properties to store. More...

MaterialProperty< RankTwoTensor > &	_strain_increment

MaterialProperty< RankTwoTensor > &	_rotation_increment


const MaterialProperty< RankTwoTensor > &	_deformation_gradient_old
	Material properties to fetch. More...

const MaterialProperty< RankTwoTensor > &	_mechanical_strain_old

const MaterialProperty< RankTwoTensor > &	_total_strain_old

std::vector< const MaterialProperty< RankTwoTensor > * >	_eigenstrains_old


const MaterialProperty< RankFourTensor > &	_Cijkl
	Material properties to fetch. More...

std::vector< MaterialPropertyName >	_eigenstrain_names

std::vector< const MaterialProperty< RankTwoTensor > * >	_eigenstrains


MaterialProperty< RankTwoTensor > &	_shape2
	Material properties to store. More...

MaterialProperty< RankTwoTensor > &	_deformation_gradient

MaterialProperty< RankTwoTensor > &	_ddgraddu

MaterialProperty< RankTwoTensor > &	_ddgraddv

MaterialProperty< RankTwoTensor > &	_ddgraddw

MaterialProperty< RankTwoTensor > &	_total_strain

MaterialProperty< RankTwoTensor > &	_mechanical_strain

MaterialProperty< Real > &	_multi

Private Attributes

const bool	_scalar_out_of_plane_strain_coupled
	Scalar out-of-plane strain for generalized plane strain. More...

const VariableValue &	_scalar_out_of_plane_strain

const VariableValue &	_scalar_out_of_plane_strain_old


const bool	_out_of_plane_strain_coupled
	Out-of-plane strain for weak plane stress. More...

const VariableValue &	_out_of_plane_strain

const VariableValue &	_out_of_plane_strain_old

const DecompMethod	_decomposition_method

enum	DecompMethod { DecompMethod::TaylorExpansion, DecompMethod::EigenSolution }
	Method to decompose into rotation increment and strain increment. More...

Detailed Description

Material class for 2D correspondence material model for finite strain: plane strain, generalized plane strain, weak plane stress.

Definition at line 23 of file ComputePlaneFiniteStrainNOSPD.h.

Member Enumeration Documentation

◆ DecompMethod

enum ComputeFiniteStrainNOSPD::DecompMethod

strongprivateinherited

Method to decompose into rotation increment and strain increment.

Enumerator
TaylorExpansion
EigenSolution

Definition at line 63 of file ComputeFiniteStrainNOSPD.h.

   {
     TaylorExpansion,
     EigenSolution
   };

Constructor & Destructor Documentation

◆ ComputePlaneFiniteStrainNOSPD()

ComputePlaneFiniteStrainNOSPD::ComputePlaneFiniteStrainNOSPD ( const InputParameters & parameters )

Definition at line 33 of file ComputePlaneFiniteStrainNOSPD.C.

   : ComputeFiniteStrainNOSPD(parameters),
     _scalar_out_of_plane_strain_coupled(isCoupledScalar("scalar_out_of_plane_strain")),
     _scalar_out_of_plane_strain(_scalar_out_of_plane_strain_coupled
                                     ? coupledScalarValue("scalar_out_of_plane_strain")
                                     : _zero),
     _scalar_out_of_plane_strain_old(_scalar_out_of_plane_strain_coupled
                                         ? coupledScalarValueOld("scalar_out_of_plane_strain")
                                         : _zero),
     _out_of_plane_strain_coupled(isCoupled("out_of_plane_strain")),
     _out_of_plane_strain(_out_of_plane_strain_coupled ? coupledValue("out_of_plane_strain")
                                                       : _zero),
     _out_of_plane_strain_old(_out_of_plane_strain_coupled ? coupledValueOld("out_of_plane_strain")
                                                           : _zero)
 {
 }

Member Function Documentation

◆ computeBondStretch()

void ComputeStrainBaseNOSPD::computeBondStretch ( )

overrideprotectedvirtualinherited

Definition at line 156 of file ComputeStrainBaseNOSPD.C.

 {
   _total_stretch[0] = _current_length / _origin_length - 1.0;
   _total_stretch[1] = _total_stretch[0];
   _mechanical_stretch[0] = _total_stretch[0];
  
   Real factor = 1.0;
   if (_plane_strain)
     factor = 1.0 + ElasticityTensorTools::getIsotropicPoissonsRatio(_Cijkl[0]);
  
   for (auto es : _eigenstrains)
     _mechanical_stretch[0] -= 0.5 * factor * ((*es)[0](2, 2) + (*es)[1](2, 2));
  
   _mechanical_stretch[1] = _mechanical_stretch[0];
 }

Referenced by ComputeStrainBaseNOSPD::computeProperties().

◆ computeProperties()

void ComputeStrainBaseNOSPD::computeProperties ( )

overrideprotectedvirtualinherited

Definition at line 145 of file ComputeStrainBaseNOSPD.C.

 {
   setupMeshRelatedData();     // function from base class
   computeBondCurrentLength(); // current length of a bond from base class
   computeBondStretch();       // stretch of a bond
  
   for (_qp = 0; _qp < _nnodes; ++_qp)
     computeQpStrain();
 }

◆ computeQpDeformationGradient()

void ComputeStrainBaseNOSPD::computeQpDeformationGradient ( )

protectedvirtualinherited

Function to compute bond-associated deformation gradient.

Reimplemented in ComputeForceStabilizedSmallStrainNOSPD.

Definition at line 71 of file ComputeStrainBaseNOSPD.C.

 {
   _shape2[_qp].zero();
   _deformation_gradient[_qp].zero();
   _ddgraddu[_qp].zero();
   _ddgraddv[_qp].zero();
   _ddgraddw[_qp].zero();
   _multi[_qp] = 0.0;
  
   // for cases when current bond was broken, assign shape tensor and deformation gradient to unity
   if (_bond_status_var->getElementalValue(_current_elem) < 0.5)
   {
     _shape2[_qp] = RankTwoTensor::initIdentity;
     _deformation_gradient[_qp] = _shape2[_qp];
   }
   else
   {
     if (_dim == 2)
       _shape2[_qp](2, 2) = _deformation_gradient[_qp](2, 2) = 1.0;
  
     const Node * cur_nd = _current_elem->node_ptr(_qp);
     const Node * end_nd = _current_elem->node_ptr(1 - _qp); // two nodes for edge2 element
     std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(cur_nd->id());
     std::vector<dof_id_type> bonds = _pdmesh.getBonds(cur_nd->id());
  
     unsigned int nb =
         std::find(neighbors.begin(), neighbors.end(), end_nd->id()) - neighbors.begin();
     std::vector<dof_id_type> dg_neighbors = _pdmesh.getDefGradNeighbors(cur_nd->id(), nb);
  
     // calculate the shape tensor and prepare the deformation gradient tensor
     Real dgnodes_vsum = 0.0;
     RealGradient ori_vec(_dim), cur_vec(_dim);
  
     for (unsigned int j = 0; j < dg_neighbors.size(); ++j)
       if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[j]])) > 0.5)
       {
         Node * node_j = _pdmesh.nodePtr(neighbors[dg_neighbors[j]]);
         Real vol_j = _pdmesh.getPDNodeVolume(neighbors[dg_neighbors[j]]);
         dgnodes_vsum += vol_j;
         ori_vec = *node_j - *_pdmesh.nodePtr(cur_nd->id());
  
         for (unsigned int k = 0; k < _dim; ++k)
           cur_vec(k) = ori_vec(k) + _disp_var[k]->getNodalValue(*node_j) -
                        _disp_var[k]->getNodalValue(*cur_nd);
  
         Real ori_len = ori_vec.norm();
         for (unsigned int k = 0; k < _dim; ++k)
         {
           for (unsigned int l = 0; l < _dim; ++l)
           {
             _shape2[_qp](k, l) += _horiz_rad[_qp] / ori_len * ori_vec(k) * ori_vec(l) * vol_j;
             _deformation_gradient[_qp](k, l) +=
                 _horiz_rad[_qp] / ori_len * cur_vec(k) * ori_vec(l) * vol_j;
           }
           // calculate derivatives of deformation_gradient w.r.t displacements of node i
           _ddgraddu[_qp](0, k) += -_horiz_rad[_qp] / ori_len * ori_vec(k) * vol_j;
           _ddgraddv[_qp](1, k) += -_horiz_rad[_qp] / ori_len * ori_vec(k) * vol_j;
           if (_dim == 3)
             _ddgraddw[_qp](2, k) += -_horiz_rad[_qp] / ori_len * ori_vec(k) * vol_j;
         }
       }
     // finalize the deformation gradient and its derivatives
     _deformation_gradient[_qp] *= _shape2[_qp].inverse();
     _ddgraddu[_qp] *= _shape2[_qp].inverse();
     _ddgraddv[_qp] *= _shape2[_qp].inverse();
     _ddgraddw[_qp] *= _shape2[_qp].inverse();
  
     // force state multiplier
     _multi[_qp] = _horiz_rad[_qp] / _origin_length * _node_vol[0] * _node_vol[1] * dgnodes_vsum /
                   _horiz_vol[_qp];
   }
 }

Referenced by ComputeSmallStrainNOSPD::computeQpStrain(), and ComputeFiniteStrainNOSPD::computeQpStrain().

◆ computeQpFhat()

void ComputePlaneFiniteStrainNOSPD::computeQpFhat ( )

overrideprotectedvirtual

Reimplemented from ComputeFiniteStrainNOSPD.

Definition at line 51 of file ComputePlaneFiniteStrainNOSPD.C.

 {
   _deformation_gradient[_qp](2, 2) = computeQpOutOfPlaneDeformationGradient();
  
   RankTwoTensor deformation_gradient_old = _deformation_gradient_old[_qp];
   deformation_gradient_old(2, 2) = computeQpOutOfPlaneDeformationGradientOld();
  
   // Incremental deformation gradient of current step w.r.t previous step:
   // _Fhat = deformation_gradient * inv(deformation_gradient_old)
   _Fhat[_qp] = _deformation_gradient[_qp] * deformation_gradient_old.inverse();
 }

◆ computeQpOutOfPlaneDeformationGradient()

Real ComputePlaneFiniteStrainNOSPD::computeQpOutOfPlaneDeformationGradient ( )

protectedvirtual

Functions to compute the out-of-plane component of deformation gradient for generalized plane strain and weak plane stress.

Definition at line 64 of file ComputePlaneFiniteStrainNOSPD.C.

 {
   // This is consistent with the approximation of stretch rate tensor
   // D = log(sqrt(Fhat^T * Fhat)) / dt
  
   if (_scalar_out_of_plane_strain_coupled)
     return std::exp(_scalar_out_of_plane_strain[0]);
   else
     return std::exp(_out_of_plane_strain[_qp]);
 }

Referenced by computeQpFhat().

◆ computeQpOutOfPlaneDeformationGradientOld()

Real ComputePlaneFiniteStrainNOSPD::computeQpOutOfPlaneDeformationGradientOld ( )

protectedvirtual

Definition at line 76 of file ComputePlaneFiniteStrainNOSPD.C.

 {
   if (_scalar_out_of_plane_strain_coupled)
     return std::exp(_scalar_out_of_plane_strain_old[0]);
   else
     return std::exp(_out_of_plane_strain_old[_qp]);
 }

Referenced by computeQpFhat().

◆ computeQpStrain()

void ComputeFiniteStrainNOSPD::computeQpStrain ( )

overrideprotectedvirtualinherited

Function to compute strain tensors.

Implements ComputeStrainBaseNOSPD.

Definition at line 56 of file ComputeFiniteStrainNOSPD.C.

 {
   computeQpDeformationGradient();
  
   computeQpFhat();
  
   RankTwoTensor total_strain_increment;
  
   // Two ways to calculate these increments: TaylorExpansion(default) or EigenSolution
   computeQpStrainRotationIncrements(total_strain_increment, _rotation_increment[_qp]);
  
   _strain_increment[_qp] = total_strain_increment;
  
   // Remove the eigenstrain increment
   subtractEigenstrainIncrementFromStrain(_strain_increment[_qp]);
  
   if (_dt > 0)
     _strain_rate[_qp] = _strain_increment[_qp] / _dt;
   else
     _strain_rate[_qp].zero();
  
   // Update strain in intermediate configuration
   _mechanical_strain[_qp] = _mechanical_strain_old[_qp] + _strain_increment[_qp];
   _total_strain[_qp] = _total_strain_old[_qp] + total_strain_increment;
  
   // Rotate strain to current configuration
   _mechanical_strain[_qp] =
       _rotation_increment[_qp] * _mechanical_strain[_qp] * _rotation_increment[_qp].transpose();
   _total_strain[_qp] =
       _rotation_increment[_qp] * _total_strain[_qp] * _rotation_increment[_qp].transpose();
 }

◆ computeQpStrainRotationIncrements()

void ComputeFiniteStrainNOSPD::computeQpStrainRotationIncrements	(	RankTwoTensor &	e,
		RankTwoTensor &	r
	)

protectedinherited

Function to compute strain and rotational increments.

Definition at line 97 of file ComputeFiniteStrainNOSPD.C.

 {
   switch (_decomposition_method)
   {
     case DecompMethod::TaylorExpansion:
     {
       // inverse of _Fhat
       RankTwoTensor invFhat(_Fhat[_qp].inverse());
  
       // A = I - _Fhat^-1
       RankTwoTensor A(RankTwoTensor::initIdentity);
       A -= invFhat;
  
       // Cinv - I = A A^T - A - A^T;
       RankTwoTensor Cinv_I = A * A.transpose() - A - A.transpose();
  
       // strain rate D from Taylor expansion, Chat = (-1/2(Chat^-1 - I) + 1/4*(Chat^-1 - I)^2 + ...
       total_strain_increment = -Cinv_I * 0.5 + Cinv_I * Cinv_I * 0.25;
  
       const Real a[3] = {invFhat(1, 2) - invFhat(2, 1),
                          invFhat(2, 0) - invFhat(0, 2),
                          invFhat(0, 1) - invFhat(1, 0)};
  
       Real q = (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]) / 4.0;
       Real trFhatinv_1 = invFhat.trace() - 1.0;
       const Real p = trFhatinv_1 * trFhatinv_1 / 4.0;
  
       // cos theta_a
       const Real C1 =
           std::sqrt(p + 3.0 * Utility::pow<2>(p) * (1.0 - (p + q)) / Utility::pow<2>(p + q) -
                     2.0 * Utility::pow<3>(p) * (1.0 - (p + q)) / Utility::pow<3>(p + q));
  
       Real C2;
       if (q > 0.01)
         // (1-cos theta_a)/4q
         C2 = (1.0 - C1) / (4.0 * q);
       else
         // alternate form for small q
         C2 = 0.125 + q * 0.03125 * (Utility::pow<2>(p) - 12.0 * (p - 1.0)) / Utility::pow<2>(p) +
              Utility::pow<2>(q) * (p - 2.0) * (Utility::pow<2>(p) - 10.0 * p + 32.0) /
                  Utility::pow<3>(p) +
              Utility::pow<3>(q) *
                  (1104.0 - 992.0 * p + 376.0 * Utility::pow<2>(p) - 72.0 * Utility::pow<3>(p) +
                   5.0 * Utility::pow<4>(p)) /
                  (512.0 * Utility::pow<4>(p));
       const Real C3 =
           0.5 * std::sqrt((p * q * (3.0 - q) + Utility::pow<3>(p) + Utility::pow<2>(q)) /
                           Utility::pow<3>(p + q)); // sin theta_a/(2 sqrt(q))
  
       // Calculate incremental rotation. Note that this value is the transpose of that from Rashid,
       // 93, so we transpose it before storing
       RankTwoTensor R_incr;
       R_incr.addIa(C1);
       for (unsigned int i = 0; i < 3; ++i)
         for (unsigned int j = 0; j < 3; ++j)
           R_incr(i, j) += C2 * a[i] * a[j];
  
       R_incr(0, 1) += C3 * a[2];
       R_incr(0, 2) -= C3 * a[1];
       R_incr(1, 0) -= C3 * a[2];
       R_incr(1, 2) += C3 * a[0];
       R_incr(2, 0) += C3 * a[1];
       R_incr(2, 1) -= C3 * a[0];
  
       rotation_increment = R_incr.transpose();
       break;
     }
     case DecompMethod::EigenSolution:
     {
       std::vector<Real> e_value(3);
       RankTwoTensor e_vector, N1, N2, N3;
  
       RankTwoTensor Chat = _Fhat[_qp].transpose() * _Fhat[_qp];
       Chat.symmetricEigenvaluesEigenvectors(e_value, e_vector);
  
       const Real lambda1 = std::sqrt(e_value[0]);
       const Real lambda2 = std::sqrt(e_value[1]);
       const Real lambda3 = std::sqrt(e_value[2]);
  
       N1.vectorOuterProduct(e_vector.column(0), e_vector.column(0));
       N2.vectorOuterProduct(e_vector.column(1), e_vector.column(1));
       N3.vectorOuterProduct(e_vector.column(2), e_vector.column(2));
  
       RankTwoTensor Uhat = N1 * lambda1 + N2 * lambda2 + N3 * lambda3;
       RankTwoTensor invUhat(Uhat.inverse());
  
       rotation_increment = _Fhat[_qp] * invUhat;
  
       total_strain_increment =
           N1 * std::log(lambda1) + N2 * std::log(lambda2) + N3 * std::log(lambda3);
       break;
     }
  
     default:
       mooseError("ComputeFiniteStrainNOSPD Error: Invalid decomposition type! Please specify : "
                  "TaylorExpansion or "
                  "EigenSolution.");
   }
 }

Referenced by ComputeFiniteStrainNOSPD::computeQpStrain().

◆ decompositionType()

MooseEnum ComputeFiniteStrainNOSPD::decompositionType ( )

staticinherited

Definition at line 17 of file ComputeFiniteStrainNOSPD.C.

 {
   return MooseEnum("TaylorExpansion EigenSolution", "TaylorExpansion");
 }

Referenced by validParams< ComputeFiniteStrainNOSPD >().

◆ initQpStatefulProperties()

void ComputeStrainBaseNOSPD::initQpStatefulProperties ( )

overridevirtualinherited

Definition at line 57 of file ComputeStrainBaseNOSPD.C.

 {
   _mechanical_strain[_qp].zero();
   _total_strain[_qp].zero();
   _deformation_gradient[_qp].zero();
   _deformation_gradient[_qp].addIa(1.0);
  
   if (_qrule->n_points() < 2) // Gauss_Lobatto: order = 2n-3 (n is number of qp)
     mooseError(
         "NOSPD models require Gauss_Lobatto rule and a minimum of 2 quadrature points, i.e., "
         "order of FIRST");
 }

◆ subtractEigenstrainIncrementFromStrain()

void ComputeFiniteStrainNOSPD::subtractEigenstrainIncrementFromStrain ( RankTwoTensor & strain )

protectedinherited

Function to compute the mechanical strain tensor by subtracting thermal strain from the total strain.

Definition at line 199 of file ComputeFiniteStrainNOSPD.C.

 {
   for (unsigned int i = 0; i < _eigenstrains.size(); ++i)
   {
     strain -= (*_eigenstrains[i])[_qp];
     strain += (*_eigenstrains_old[i])[_qp];
   }
 }