Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://www.mooseframework.org
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "EshelbyTensor.h"
      11             : #include "RankTwoTensor.h"
      12             : #include "MooseMesh.h"
      13             : 
      14             : registerMooseObject("TensorMechanicsApp", EshelbyTensor);
      15             : registerMooseObject("TensorMechanicsApp", ADEshelbyTensor);
      16             : 
      17             : template <bool is_ad>
      18             : InputParameters
      19         628 : EshelbyTensorTempl<is_ad>::validParams()
      20             : {
      21         628 :   InputParameters params = Material::validParams();
      22         628 :   params.addClassDescription("Computes the Eshelby tensor as a function of "
      23             :                              "strain energy density and the first "
      24             :                              "Piola-Kirchhoff stress");
      25        1256 :   params.addRequiredCoupledVar(
      26             :       "displacements",
      27             :       "The displacements appropriate for the simulation geometry and coordinate system");
      28        1256 :   params.addParam<std::string>("base_name",
      29             :                                "Optional parameter that allows the user to define "
      30             :                                "multiple mechanics material systems on the same "
      31             :                                "block, i.e. for multiple phases");
      32        1256 :   params.addParam<bool>(
      33             :       "compute_dissipation",
      34        1256 :       false,
      35             :       "Whether to compute Eshelby tensor's dissipation (or rate of change). This tensor"
      36             :       "yields the increase in dissipation per unit crack advanced");
      37        1256 :   params.addCoupledVar("temperature", "Coupled temperature");
      38         628 :   return params;
      39           0 : }
      40             : 
      41             : template <bool is_ad>
      42         471 : EshelbyTensorTempl<is_ad>::EshelbyTensorTempl(const InputParameters & parameters)
      43             :   : DerivativeMaterialInterface<Material>(parameters),
      44         471 :     _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
      45         942 :     _compute_dissipation(getParam<bool>("compute_dissipation")),
      46         942 :     _sed(getMaterialPropertyByName<Real>(_base_name + "strain_energy_density")),
      47         942 :     _serd(_compute_dissipation
      48         489 :               ? &getMaterialPropertyByName<Real>(_base_name + "strain_energy_rate_density")
      49             :               : nullptr),
      50         471 :     _eshelby_tensor(declareProperty<RankTwoTensor>(_base_name + "Eshelby_tensor")),
      51         471 :     _eshelby_tensor_dissipation(
      52         471 :         _compute_dissipation
      53         471 :             ? &declareProperty<RankTwoTensor>(_base_name + "Eshelby_tensor_dissipation")
      54             :             : nullptr),
      55         471 :     _stress(getGenericMaterialProperty<RankTwoTensor, is_ad>(_base_name + "stress")),
      56         942 :     _stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "stress")),
      57         471 :     _grad_disp(3),
      58         471 :     _grad_disp_old(3),
      59         471 :     _J_thermal_term_vec(declareProperty<RealVectorValue>("J_thermal_term_vec")),
      60         471 :     _grad_temp(coupledGradient("temperature")),
      61         471 :     _has_temp(isCoupled("temperature")),
      62        1413 :     _total_deigenstrain_dT(getOptionalMaterialProperty<RankTwoTensor>("total_deigenstrain_dT"))
      63             : {
      64         471 :   unsigned int ndisp = coupledComponents("displacements");
      65             : 
      66             :   // Checking for consistency between mesh size and length of the provided displacements vector
      67         471 :   if (ndisp != _mesh.dimension())
      68           0 :     mooseError(
      69             :         "The number of variables supplied in 'displacements' must match the mesh dimension.");
      70             : 
      71             :   // fetch coupled gradients
      72        1632 :   for (unsigned int i = 0; i < ndisp; ++i)
      73        1161 :     _grad_disp[i] = &coupledGradient("displacements", i);
      74             : 
      75             :   // set unused dimensions to zero
      76         723 :   for (unsigned i = ndisp; i < 3; ++i)
      77         252 :     _grad_disp[i] = &_grad_zero;
      78             : 
      79             :   // Need previous step's displacements to compute deformation gradient time rate
      80         471 :   if (_compute_dissipation)
      81             :   {
      82             :     // fetch coupled gradients previous step
      83          54 :     for (unsigned int i = 0; i < ndisp; ++i)
      84          36 :       _grad_disp_old[i] = &coupledGradientOld("displacements", i);
      85             : 
      86             :     // set unused dimensions to zero
      87          36 :     for (unsigned i = ndisp; i < 3; ++i)
      88          18 :       _grad_disp_old[i] = &_grad_zero;
      89             :   }
      90         471 : }
      91             : 
      92             : template <bool is_ad>
      93             : void
      94         468 : EshelbyTensorTempl<is_ad>::initialSetup()
      95             : {
      96         468 :   if (_has_temp && !_total_deigenstrain_dT)
      97           0 :     mooseError("EshelbyTensor Error: To include thermal strain term in Fracture integral "
      98             :                "calculation, must both couple temperature in DomainIntegral block and compute "
      99             :                "total_deigenstrain_dT using ThermalFractureIntegral material model.");
     100         468 : }
     101             : 
     102             : template <bool is_ad>
     103             : void
     104      262776 : EshelbyTensorTempl<is_ad>::initQpStatefulProperties()
     105             : {
     106      262776 : }
     107             : 
     108             : template <bool is_ad>
     109             : void
     110     7698784 : EshelbyTensorTempl<is_ad>::computeQpProperties()
     111             : {
     112             :   // Deformation gradient
     113     7698784 :   auto F = RankTwoTensor::initializeFromRows(
     114             :       (*_grad_disp[0])[_qp], (*_grad_disp[1])[_qp], (*_grad_disp[2])[_qp]);
     115             : 
     116     7698784 :   RankTwoTensor H(F);
     117     7698784 :   F.addIa(1.0);
     118     7698784 :   Real detF = F.det();
     119     7698784 :   RankTwoTensor FinvT(F.inverse().transpose());
     120             : 
     121             :   // 1st Piola-Kirchoff Stress (P):
     122     7698784 :   RankTwoTensor P = detF * MetaPhysicL::raw_value(_stress[_qp]) * FinvT;
     123             : 
     124             :   // HTP = H^T * P = H^T * detF * sigma * FinvT;
     125     7698784 :   RankTwoTensor HTP = H.transpose() * P;
     126             : 
     127     7698784 :   RankTwoTensor WI = RankTwoTensor(RankTwoTensor::initIdentity);
     128     7698784 :   WI *= (_sed[_qp] * detF);
     129             : 
     130     7698784 :   _eshelby_tensor[_qp] = WI - HTP;
     131             : 
     132             :   // Compute deformation gradient rate
     133     7698784 :   if (_compute_dissipation == true)
     134             :   {
     135      108880 :     auto H_old = RankTwoTensor::initializeFromRows(
     136             :         (*_grad_disp_old[0])[_qp], (*_grad_disp_old[1])[_qp], (*_grad_disp_old[2])[_qp]);
     137             : 
     138      108880 :     RankTwoTensor Wdot = RankTwoTensor(RankTwoTensor::initIdentity);
     139      108880 :     Wdot *= ((*_serd)[_qp] * detF);
     140             : 
     141             :     // F_dot = (F - F_old)/dt
     142      108880 :     RankTwoTensor F_dot = (H - H_old) / _dt;
     143             : 
     144             :     // FdotTP = Fdot^T * P = Fdot^T * detF * sigma * FinvT;
     145      108880 :     RankTwoTensor FdotTP = F_dot.transpose() * P;
     146             : 
     147      108880 :     (*_eshelby_tensor_dissipation)[_qp] = Wdot - FdotTP;
     148             :   }
     149             : 
     150     7698784 :   if (_has_temp)
     151             :   {
     152             :     const Real sigma_alpha =
     153      318400 :         MetaPhysicL::raw_value(_stress[_qp]).doubleContraction(_total_deigenstrain_dT[_qp]);
     154      318400 :     _J_thermal_term_vec[_qp] = sigma_alpha * _grad_temp[_qp];
     155             :   }
     156             :   else
     157     7380384 :     _J_thermal_term_vec[_qp].zero();
     158     7698784 : }