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 "ComputeLagrangianStrain.h"
      11             : 
      12             : template <class G>
      13             : InputParameters
      14        5234 : ComputeLagrangianStrainBase<G>::baseParams()
      15             : {
      16        5234 :   InputParameters params = Material::validParams();
      17             : 
      18       10468 :   params.addRequiredCoupledVar("displacements", "Displacement variables");
      19       10468 :   params.addParam<bool>(
      20       10468 :       "large_kinematics", false, "Use large displacement kinematics in the kernel.");
      21       10468 :   params.addParam<bool>("stabilize_strain", false, "Average the volumetric strains");
      22       10468 :   params.addParam<std::vector<MaterialPropertyName>>(
      23             :       "eigenstrain_names", {}, "List of eigenstrains to account for");
      24       10468 :   params.addParam<std::vector<MaterialPropertyName>>(
      25             :       "homogenization_gradient_names",
      26             :       {},
      27             :       "List of homogenization gradients to add to the displacement gradient");
      28             : 
      29       10468 :   params.addParam<std::string>("base_name", "Material property base name");
      30             : 
      31             :   // We rely on this *not* having use_displaced mesh on
      32        5234 :   params.suppressParameter<bool>("use_displaced_mesh");
      33             : 
      34        5234 :   return params;
      35           0 : }
      36             : 
      37             : template <class G>
      38        3924 : ComputeLagrangianStrainBase<G>::ComputeLagrangianStrainBase(const InputParameters & parameters)
      39             :   : Material(parameters),
      40        3924 :     _ndisp(coupledComponents("displacements")),
      41        3924 :     _disp(coupledValues("displacements")),
      42        3924 :     _grad_disp(coupledGradients("displacements")),
      43        7920 :     _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
      44        7848 :     _large_kinematics(getParam<bool>("large_kinematics")),
      45        7848 :     _stabilize_strain(getParam<bool>("stabilize_strain")),
      46        7848 :     _eigenstrain_names(getParam<std::vector<MaterialPropertyName>>("eigenstrain_names")),
      47        3924 :     _eigenstrains(_eigenstrain_names.size()),
      48        3924 :     _eigenstrains_old(_eigenstrain_names.size()),
      49        3924 :     _total_strain(declareProperty<RankTwoTensor>(_base_name + "total_strain")),
      50        7848 :     _total_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "total_strain")),
      51        3924 :     _mechanical_strain(declareProperty<RankTwoTensor>(_base_name + "mechanical_strain")),
      52        7848 :     _mechanical_strain_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "mechanical_strain")),
      53        3924 :     _strain_increment(declareProperty<RankTwoTensor>(_base_name + "strain_increment")),
      54        3924 :     _spatial_velocity_increment(
      55        3924 :         declareProperty<RankTwoTensor>(_base_name + "spatial_velocity_increment")),
      56        3924 :     _vorticity_increment(declareProperty<RankTwoTensor>(_base_name + "vorticity_increment")),
      57        3924 :     _F_ust(declareProperty<RankTwoTensor>(_base_name + "unstabilized_deformation_gradient")),
      58        3924 :     _F_avg(declareProperty<RankTwoTensor>(_base_name + "average_deformation_gradient")),
      59        3924 :     _F(declareProperty<RankTwoTensor>(_base_name + "deformation_gradient")),
      60        7848 :     _F_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "deformation_gradient")),
      61        3924 :     _F_inv(declareProperty<RankTwoTensor>(_base_name + "inverse_deformation_gradient")),
      62        3924 :     _f_inv(declareProperty<RankTwoTensor>(_base_name + "inverse_incremental_deformation_gradient")),
      63        7848 :     _homogenization_gradient_names(
      64             :         getParam<std::vector<MaterialPropertyName>>("homogenization_gradient_names")),
      65        3924 :     _homogenization_contributions(_homogenization_gradient_names.size()),
      66        7848 :     _rotation_increment(declareProperty<RankTwoTensor>(_base_name + "rotation_increment"))
      67             : {
      68             :   // Setup eigenstrains
      69        4272 :   for (auto i : make_range(_eigenstrain_names.size()))
      70             :   {
      71         348 :     _eigenstrains[i] = &getMaterialProperty<RankTwoTensor>(_eigenstrain_names[i]);
      72         348 :     _eigenstrains_old[i] = &getMaterialPropertyOld<RankTwoTensor>(_eigenstrain_names[i]);
      73             :   }
      74             : 
      75             :   // In the future maybe there is a reason to have more than one, but for now
      76        3924 :   if (_homogenization_gradient_names.size() > 1)
      77           0 :     mooseError("ComputeLagrangianStrainBase cannot accommodate more than one "
      78             :                "homogenization gradient");
      79             : 
      80             :   // Setup homogenization contributions
      81        4458 :   for (unsigned int i = 0; i < _homogenization_gradient_names.size(); i++)
      82         534 :     _homogenization_contributions[i] =
      83         534 :         &getMaterialProperty<RankTwoTensor>(_homogenization_gradient_names[i]);
      84        3924 : }
      85             : 
      86             : template <class G>
      87             : void
      88      930764 : ComputeLagrangianStrainBase<G>::initQpStatefulProperties()
      89             : {
      90      930764 :   _total_strain[_qp].zero();
      91      930764 :   _mechanical_strain[_qp].zero();
      92      930764 :   _F[_qp].setToIdentity();
      93      930764 : }
      94             : 
      95             : template <class G>
      96             : void
      97    11804890 : ComputeLagrangianStrainBase<G>::computeProperties()
      98             : {
      99             :   // Average the volumetric terms, if required
     100    11804890 :   computeDeformationGradient();
     101             : 
     102    96517238 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
     103    84712348 :     computeQpProperties();
     104    11804890 : }
     105             : 
     106             : template <class G>
     107             : void
     108    84712348 : ComputeLagrangianStrainBase<G>::computeQpProperties()
     109             : {
     110             :   // Add in the macroscale gradient contribution
     111   121748492 :   for (auto contribution : _homogenization_contributions)
     112    37036144 :     _F[_qp] += (*contribution)[_qp];
     113             : 
     114             :   // If the kernel is large deformation then we need the "actual"
     115             :   // kinematic quantities
     116    84712348 :   RankTwoTensor dL;
     117    84712348 :   if (_large_kinematics)
     118             :   {
     119    39768136 :     _F_inv[_qp] = _F[_qp].inverse();
     120    39768136 :     _f_inv[_qp] = _F_old[_qp] * _F_inv[_qp];
     121    39768136 :     dL = RankTwoTensor::Identity() - _f_inv[_qp];
     122             :   }
     123             :   // For small deformations we just provide the identity
     124             :   else
     125             :   {
     126    44944212 :     _F_inv[_qp] = RankTwoTensor::Identity();
     127    44944212 :     _f_inv[_qp] = RankTwoTensor::Identity();
     128    44944212 :     dL = _F[_qp] - _F_old[_qp];
     129             :   }
     130             : 
     131    84712348 :   computeQpIncrementalStrains(dL);
     132    84712348 : }
     133             : 
     134             : template <class G>
     135             : void
     136    84712348 : ComputeLagrangianStrainBase<G>::computeQpIncrementalStrains(const RankTwoTensor & dL)
     137             : {
     138             :   // Get the deformation increments
     139    84712348 :   _strain_increment[_qp] = (dL + dL.transpose()) / 2.0;
     140    84712348 :   _vorticity_increment[_qp] = (dL - dL.transpose()) / 2.0;
     141             : 
     142             :   // Increment the total strain
     143    84712348 :   _total_strain[_qp] = _total_strain_old[_qp] + _strain_increment[_qp];
     144             : 
     145             :   // Get rid of the eigenstrains
     146             :   // Note we currently do not alter the deformation gradient -- this will be
     147             :   // needed in the future for a "complete" system
     148    84712348 :   subtractQpEigenstrainIncrement(_strain_increment[_qp]);
     149             : 
     150             :   // Increment the mechanical strain
     151    84712348 :   _mechanical_strain[_qp] = _mechanical_strain_old[_qp] + _strain_increment[_qp];
     152             : 
     153             :   // Yes, this does make sense to do it here
     154    84712348 :   _spatial_velocity_increment[_qp] = _vorticity_increment[_qp] + _strain_increment[_qp];
     155             : 
     156             :   // Faked rotation increment for ComputeStressBase materials
     157    84712348 :   _rotation_increment[_qp] = RankTwoTensor::Identity();
     158    84712348 : }
     159             : 
     160             : template <class G>
     161             : void
     162    84712348 : ComputeLagrangianStrainBase<G>::subtractQpEigenstrainIncrement(RankTwoTensor & strain)
     163             : {
     164    85529756 :   for (auto i : make_range(_eigenstrain_names.size()))
     165      817408 :     strain -= (*_eigenstrains[i])[_qp] - (*_eigenstrains_old[i])[_qp];
     166    84712348 : }
     167             : 
     168             : template <class G>
     169             : void
     170    84712348 : ComputeLagrangianStrainBase<G>::computeQpUnstabilizedDeformationGradient()
     171             : {
     172    84712348 :   _F_ust[_qp].setToIdentity();
     173   336902992 :   for (auto component : make_range(_ndisp))
     174   252190644 :     G::addGradOp(_F_ust[_qp],
     175             :                  component,
     176             :                  (*_grad_disp[component])[_qp],
     177   252190644 :                  (*_disp[component])[_qp],
     178      599004 :                  _q_point[_qp]);
     179    84712348 : }
     180             : 
     181             : template <class G>
     182             : void
     183    11804890 : ComputeLagrangianStrainBase<G>::computeDeformationGradient()
     184             : {
     185             :   // First calculate the unstabilized deformation gradient at each qp
     186    96517238 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
     187             :   {
     188    84712348 :     computeQpUnstabilizedDeformationGradient();
     189    84712348 :     _F[_qp] = _F_ust[_qp];
     190             :   }
     191             : 
     192             :   // If stabilization is on do the volumetric correction
     193    11804890 :   if (_stabilize_strain)
     194             :   {
     195     3676644 :     const auto F_avg = StabilizationUtils::elementAverage(
     196    30941588 :         [this](unsigned int qp) { return _F_ust[qp]; }, _JxW, _coord);
     197             :     // All quadrature points have the same F_avg
     198     3676644 :     _F_avg.set().setAllValues(F_avg);
     199             :     // Make the appropriate modification, depending on small or large
     200             :     // deformations
     201    30941588 :     for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
     202             :     {
     203    27264944 :       if (_large_kinematics)
     204     9490912 :         _F[_qp] *= std::pow(F_avg.det() / _F[_qp].det(), 1.0 / 3.0);
     205             :       else
     206    17774032 :         _F[_qp] += (F_avg.trace() - _F[_qp].trace()) * RankTwoTensor::Identity() / 3.0;
     207             :     }
     208             :   }
     209    11804890 : }
     210             : 
     211             : template class ComputeLagrangianStrainBase<GradientOperatorCartesian>;
     212             : template class ComputeLagrangianStrainBase<GradientOperatorAxisymmetricCylindrical>;
     213             : template class ComputeLagrangianStrainBase<GradientOperatorCentrosymmetricSpherical>;