Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       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("SolidMechanicsApp", EshelbyTensor);
      15             : registerMooseObject("SolidMechanicsApp", ADEshelbyTensor);
      16             : 
      17             : template <bool is_ad>
      18             : InputParameters
      19        1432 : EshelbyTensorTempl<is_ad>::validParams()
      20             : {
      21        1432 :   InputParameters params = Material::validParams();
      22        1432 :   params.addClassDescription("Computes the Eshelby tensor as a function of "
      23             :                              "strain energy density and the first "
      24             :                              "Piola-Kirchhoff stress");
      25        2864 :   params.addRequiredCoupledVar(
      26             :       "displacements",
      27             :       "The displacements appropriate for the simulation geometry and coordinate system");
      28        2864 :   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        2864 :   params.addParam<bool>(
      33             :       "compute_dissipation",
      34        2864 :       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        2864 :   params.addCoupledVar("temperature", "Coupled temperature");
      38        1432 :   return params;
      39           0 : }
      40             : 
      41             : template <bool is_ad>
      42        1074 : EshelbyTensorTempl<is_ad>::EshelbyTensorTempl(const InputParameters & parameters)
      43             :   : DerivativeMaterialInterface<Material>(parameters),
      44        1074 :     _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
      45        2148 :     _compute_dissipation(getParam<bool>("compute_dissipation")),
      46        2148 :     _sed(getMaterialPropertyByName<Real>(_base_name + "strain_energy_density")),
      47        2148 :     _serd(_compute_dissipation
      48        1116 :               ? &getMaterialPropertyByName<Real>(_base_name + "strain_energy_rate_density")
      49             :               : nullptr),
      50        1074 :     _eshelby_tensor(declareProperty<RankTwoTensor>(_base_name + "Eshelby_tensor")),
      51        1074 :     _eshelby_tensor_dissipation(
      52        1074 :         _compute_dissipation
      53        1074 :             ? &declareProperty<RankTwoTensor>(_base_name + "Eshelby_tensor_dissipation")
      54             :             : nullptr),
      55        1074 :     _stress(getGenericMaterialProperty<RankTwoTensor, is_ad>(_base_name + "stress")),
      56        2148 :     _stress_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "stress")),
      57        1074 :     _grad_disp(3),
      58        1074 :     _grad_disp_old(3),
      59        1074 :     _J_thermal_term_vec(declareProperty<RealVectorValue>("J_thermal_term_vec")),
      60        1074 :     _grad_temp(coupledGradient("temperature")),
      61        1074 :     _has_temp(isCoupled("temperature")),
      62        3222 :     _total_deigenstrain_dT(getOptionalMaterialProperty<RankTwoTensor>("total_deigenstrain_dT"))
      63             : {
      64        1074 :   unsigned int ndisp = coupledComponents("displacements");
      65             : 
      66             :   // Checking for consistency between mesh size and length of the provided displacements vector
      67        1074 :   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        3705 :   for (unsigned int i = 0; i < ndisp; ++i)
      73        2631 :     _grad_disp[i] = &coupledGradient("displacements", i);
      74             : 
      75             :   // set unused dimensions to zero
      76        1665 :   for (unsigned i = ndisp; i < 3; ++i)
      77         591 :     _grad_disp[i] = &_grad_zero;
      78             : 
      79             :   // Need previous step's displacements to compute deformation gradient time rate
      80        1074 :   if (_compute_dissipation)
      81             :   {
      82             :     // fetch coupled gradients previous step
      83         126 :     for (unsigned int i = 0; i < ndisp; ++i)
      84          84 :       _grad_disp_old[i] = &coupledGradientOld("displacements", i);
      85             : 
      86             :     // set unused dimensions to zero
      87          84 :     for (unsigned i = ndisp; i < 3; ++i)
      88          42 :       _grad_disp_old[i] = &_grad_zero;
      89             :   }
      90        1074 : }
      91             : 
      92             : template <bool is_ad>
      93             : void
      94        1068 : EshelbyTensorTempl<is_ad>::initialSetup()
      95             : {
      96        1068 :   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        1068 : }
     101             : 
     102             : template <bool is_ad>
     103             : void
     104           0 : EshelbyTensorTempl<is_ad>::initQpStatefulProperties()
     105             : {
     106           0 : }
     107             : 
     108             : template <bool is_ad>
     109             : void
     110    21691680 : EshelbyTensorTempl<is_ad>::computeQpProperties()
     111             : {
     112             :   // Deformation gradient
     113    21691680 :   auto F = RankTwoTensor::initializeFromRows(
     114             :       (*_grad_disp[0])[_qp], (*_grad_disp[1])[_qp], (*_grad_disp[2])[_qp]);
     115             : 
     116    21691680 :   RankTwoTensor H(F);
     117    21691680 :   F.addIa(1.0);
     118    21691680 :   Real detF = F.det();
     119    21691680 :   RankTwoTensor FinvT(F.inverse().transpose());
     120             : 
     121             :   // 1st Piola-Kirchoff Stress (P):
     122    21691680 :   RankTwoTensor P = detF * MetaPhysicL::raw_value(_stress[_qp]) * FinvT;
     123             : 
     124             :   // HTP = H^T * P = H^T * detF * sigma * FinvT;
     125    21691680 :   RankTwoTensor HTP = H.transpose() * P;
     126             : 
     127    21691680 :   RankTwoTensor WI = RankTwoTensor(RankTwoTensor::initIdentity);
     128    21691680 :   WI *= (_sed[_qp] * detF);
     129             : 
     130    21691680 :   _eshelby_tensor[_qp] = WI - HTP;
     131             : 
     132             :   // Compute deformation gradient rate
     133    21691680 :   if (_compute_dissipation == true)
     134             :   {
     135      262800 :     auto H_old = RankTwoTensor::initializeFromRows(
     136             :         (*_grad_disp_old[0])[_qp], (*_grad_disp_old[1])[_qp], (*_grad_disp_old[2])[_qp]);
     137             : 
     138      262800 :     RankTwoTensor Wdot = RankTwoTensor(RankTwoTensor::initIdentity);
     139      262800 :     Wdot *= ((*_serd)[_qp] * detF);
     140             : 
     141             :     // F_dot = (F - F_old)/dt
     142      262800 :     RankTwoTensor F_dot = (H - H_old) / _dt;
     143             : 
     144             :     // FdotTP = Fdot^T * P = Fdot^T * detF * sigma * FinvT;
     145      262800 :     RankTwoTensor FdotTP = F_dot.transpose() * P;
     146             : 
     147      262800 :     (*_eshelby_tensor_dissipation)[_qp] = Wdot - FdotTP;
     148             :   }
     149             : 
     150    21691680 :   if (_has_temp)
     151             :   {
     152             :     const Real sigma_alpha =
     153      761600 :         MetaPhysicL::raw_value(_stress[_qp]).doubleContraction(_total_deigenstrain_dT[_qp]);
     154      761600 :     _J_thermal_term_vec[_qp] = sigma_alpha * _grad_temp[_qp];
     155             :   }
     156             :   else
     157    20930080 :     _J_thermal_term_vec[_qp].zero();
     158    21691680 : }