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