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 "LinearViscoelasticityBase.h"
      11             : #include "Conversion.h"
      12             : 
      13             : InputParameters
      14          96 : LinearViscoelasticityBase::validParams()
      15             : {
      16         192 :   MooseEnum integration("backward-euler mid-point newmark zienkiewicz", "backward-euler");
      17             : 
      18          96 :   InputParameters params = ComputeElasticityTensorBase::validParams();
      19         192 :   params.addParam<MooseEnum>("integration_rule",
      20             :                              integration,
      21             :                              "describes how the viscoelastic behavior is integrated through time");
      22         192 :   params.addRangeCheckedParam<Real>("theta",
      23             :                                     1,
      24             :                                     "theta > 0 & theta <= 1",
      25             :                                     "coefficient for Newmark integration rule (between 0 and 1)");
      26         192 :   params.addParam<std::string>("driving_eigenstrain",
      27             :                                "name of the eigenstrain that increases the creep strains");
      28         192 :   params.addParam<std::string>(
      29             :       "elastic_strain_name", "elastic_strain", "name of the true elastic strain of the material");
      30         192 :   params.addParam<std::string>("creep_strain_name",
      31             :                                "creep_strain",
      32             :                                "name of the true creep strain of the material"
      33             :                                "(computed by LinearViscoelasticStressUpdate or"
      34             :                                "ComputeLinearViscoelasticStress)");
      35         192 :   params.addParam<bool>("force_recompute_properties",
      36         192 :                         false,
      37             :                         "forces the computation of the viscoelastic properties at each step of"
      38             :                         "the solver (default: false)");
      39         192 :   params.addParam<bool>(
      40             :       "need_viscoelastic_properties_inverse",
      41         192 :       false,
      42             :       "checks whether the model requires the computation of the inverse viscoelastic"
      43             :       "properties (default: false)");
      44          96 :   params.suppressParameter<FunctionName>("elasticity_tensor_prefactor");
      45          96 :   return params;
      46          96 : }
      47             : 
      48          72 : LinearViscoelasticityBase::LinearViscoelasticityBase(const InputParameters & parameters)
      49             :   : ComputeElasticityTensorBase(parameters),
      50          72 :     _integration_rule(getParam<MooseEnum>("integration_rule").getEnum<IntegrationRule>()),
      51         144 :     _theta(getParam<Real>("theta")),
      52          72 :     _apparent_elasticity_tensor(
      53          72 :         declareProperty<RankFourTensor>(_base_name + "apparent_elasticity_tensor")),
      54          72 :     _apparent_elasticity_tensor_inv(
      55          72 :         declareProperty<RankFourTensor>(_base_name + "apparent_elasticity_tensor_inv")),
      56          72 :     _elasticity_tensor_inv(declareProperty<RankFourTensor>(_elasticity_tensor_name + "_inv")),
      57         144 :     _need_viscoelastic_properties_inverse(getParam<bool>("need_viscoelastic_properties_inverse")),
      58          72 :     _has_longterm_dashpot(false),
      59          72 :     _components(0),
      60          72 :     _first_elasticity_tensor(
      61          72 :         declareProperty<RankFourTensor>(_base_name + "spring_elasticity_tensor_0")),
      62          72 :     _first_elasticity_tensor_inv(
      63          72 :         _need_viscoelastic_properties_inverse
      64         126 :             ? &declareProperty<RankFourTensor>(_base_name + "spring_elasticity_tensor_0_inv")
      65             :             : nullptr),
      66          72 :     _apparent_creep_strain(declareProperty<RankTwoTensor>(_base_name + "apparent_creep_strain")),
      67          72 :     _apparent_creep_strain_old(
      68          72 :         getMaterialPropertyOld<RankTwoTensor>(_base_name + "apparent_creep_strain")),
      69          72 :     _elastic_strain_old(
      70          72 :         getMaterialPropertyOld<RankTwoTensor>(getParam<std::string>("elastic_strain_name"))),
      71          72 :     _creep_strain_old(
      72          72 :         getMaterialPropertyOld<RankTwoTensor>(getParam<std::string>("creep_strain_name"))),
      73         144 :     _has_driving_eigenstrain(isParamValid("driving_eigenstrain")),
      74          90 :     _driving_eigenstrain_name(
      75          72 :         _has_driving_eigenstrain ? getParam<std::string>("driving_eigenstrain") : ""),
      76         144 :     _driving_eigenstrain(_has_driving_eigenstrain
      77          72 :                              ? &getMaterialPropertyByName<RankTwoTensor>(_driving_eigenstrain_name)
      78             :                              : nullptr),
      79         144 :     _driving_eigenstrain_old(_has_driving_eigenstrain
      80          72 :                                  ? &getMaterialPropertyOld<RankTwoTensor>(_driving_eigenstrain_name)
      81             :                                  : nullptr),
      82         144 :     _force_recompute_properties(getParam<bool>("force_recompute_properties")),
      83         216 :     _step_zero(declareRestartableData<bool>("step_zero", true))
      84             : {
      85          72 :   if (_theta < 0.5)
      86           0 :     mooseWarning("theta parameter for LinearViscoelasticityBase is below 0.5; time integration may "
      87             :                  "not converge!");
      88             : 
      89             :   // force material properties to be considered stateful
      90          72 :   getMaterialPropertyOld<RankFourTensor>(_base_name + "apparent_elasticity_tensor");
      91         144 :   getMaterialPropertyOld<RankFourTensor>(_base_name + "apparent_elasticity_tensor_inv");
      92             :   getMaterialPropertyOld<RankFourTensor>(_elasticity_tensor_name);
      93          72 :   getMaterialPropertyOld<RankFourTensor>(_elasticity_tensor_name + "_inv");
      94          72 :   getMaterialPropertyOld<RankFourTensor>(_base_name + "spring_elasticity_tensor_0");
      95          72 :   if (_need_viscoelastic_properties_inverse)
      96         108 :     getMaterialPropertyOld<RankFourTensor>(_base_name + "spring_elasticity_tensor_0_inv");
      97          72 : }
      98             : 
      99             : void
     100          72 : LinearViscoelasticityBase::declareViscoelasticProperties()
     101             : {
     102         198 :   for (unsigned int i = 0; i < _components; ++i)
     103             :   {
     104         126 :     std::string ith = Moose::stringify(i + 1);
     105             : 
     106         126 :     if (!_has_longterm_dashpot || (_components > 0 && i < _components - 1))
     107             :     {
     108         108 :       _springs_elasticity_tensors.push_back(
     109         108 :           &declareProperty<RankFourTensor>(_base_name + "spring_elasticity_tensor_" + ith));
     110         216 :       getMaterialPropertyOld<RankFourTensor>(_base_name + "spring_elasticity_tensor_" + ith);
     111             :     }
     112             : 
     113         126 :     _dashpot_viscosities.push_back(&declareProperty<Real>(_base_name + "dashpot_viscosity_" + ith));
     114         126 :     _dashpot_viscosities_old.push_back(
     115         252 :         &getMaterialPropertyOld<Real>(_base_name + "dashpot_viscosity_" + ith));
     116             : 
     117         126 :     _viscous_strains.push_back(
     118         126 :         &declareProperty<RankTwoTensor>(_base_name + "viscous_strain_" + ith));
     119         126 :     _viscous_strains_old.push_back(
     120         252 :         &getMaterialPropertyOld<RankTwoTensor>(_base_name + "viscous_strain_" + ith));
     121             : 
     122         126 :     if (_need_viscoelastic_properties_inverse)
     123             :     {
     124          99 :       _springs_elasticity_tensors_inv.push_back(&declareProperty<RankFourTensor>(
     125          99 :           _base_name + "spring_elasticity_tensor_" + ith + "_inv"));
     126          99 :       _springs_elasticity_tensors_inv_old.push_back(&getMaterialPropertyOld<RankFourTensor>(
     127         198 :           _base_name + "spring_elasticity_tensor_" + ith + "_inv"));
     128             :     }
     129             :   }
     130          72 : }
     131             : 
     132             : void
     133         320 : LinearViscoelasticityBase::initQpStatefulProperties()
     134             : {
     135         320 :   if (_components != _viscous_strains.size())
     136           0 :     mooseError(
     137             :         "inconsistent numbers of dashpots and viscous strains in LinearViscoelasticityBase;"
     138             :         " Make sure declareViscoelasticProperties has been called in the viscoelastic model");
     139             : 
     140         320 :   _apparent_creep_strain[_qp].zero();
     141         320 :   _apparent_elasticity_tensor[_qp].zero();
     142         320 :   _apparent_elasticity_tensor_inv[_qp].zero();
     143         320 :   _elasticity_tensor_inv[_qp].zero();
     144         320 :   _first_elasticity_tensor[_qp].zero();
     145         320 :   if (_need_viscoelastic_properties_inverse)
     146         256 :     (*_first_elasticity_tensor_inv)[_qp].zero();
     147             : 
     148         896 :   for (unsigned int i = 0; i < _components; ++i)
     149             :   {
     150         576 :     if (!_has_longterm_dashpot || (_components > 0 && i < _components - 1))
     151             :     {
     152         496 :       (*_springs_elasticity_tensors[i])[_qp].zero();
     153         496 :       if (_need_viscoelastic_properties_inverse)
     154         400 :         (*_springs_elasticity_tensors_inv[i])[_qp].zero();
     155             :     }
     156             : 
     157         576 :     (*_dashpot_viscosities[i])[_qp] = 0.0;
     158         576 :     (*_viscous_strains[i])[_qp].zero();
     159             :   }
     160         320 : }
     161             : 
     162             : void
     163        3648 : LinearViscoelasticityBase::recomputeQpApparentProperties(unsigned int qp)
     164             : {
     165        3648 :   unsigned int qp_prev = _qp;
     166        3648 :   _qp = qp;
     167             : 
     168        3648 :   if (_t_step >= 1)
     169        3648 :     _step_zero = false;
     170             : 
     171             :   // 1. we get the viscoelastic properties and their inverse if needed
     172        3648 :   computeQpViscoelasticProperties();
     173        3648 :   if (_need_viscoelastic_properties_inverse)
     174        2656 :     computeQpViscoelasticPropertiesInv();
     175             : 
     176             :   // 2. we update the internal viscous strains from the previous time step
     177        3648 :   updateQpViscousStrains();
     178             : 
     179             :   // 3. we compute the apparent elasticity tensor
     180        3648 :   computeQpApparentElasticityTensors();
     181             : 
     182             :   // 4. we transform the internal viscous strains in an apparent creep strain
     183        3648 :   if (!_step_zero)
     184        3648 :     computeQpApparentCreepStrain();
     185             : 
     186        3648 :   _qp = qp_prev;
     187        3648 : }
     188             : 
     189             : void
     190       60888 : LinearViscoelasticityBase::computeQpElasticityTensor()
     191             : {
     192       60888 :   if (_force_recompute_properties)
     193           0 :     recomputeQpApparentProperties(_qp);
     194       60888 : }
     195             : 
     196             : void
     197           0 : LinearViscoelasticityBase::computeQpViscoelasticPropertiesInv()
     198             : {
     199           0 :   if (MooseUtils::absoluteFuzzyEqual(_first_elasticity_tensor[_qp].L2norm(), 0.0))
     200           0 :     (*_first_elasticity_tensor_inv)[_qp].zero();
     201             :   else
     202           0 :     (*_first_elasticity_tensor_inv)[_qp] = _first_elasticity_tensor[_qp].invSymm();
     203             : 
     204           0 :   for (unsigned int i = 0; i < _springs_elasticity_tensors.size(); ++i)
     205             :   {
     206           0 :     if (MooseUtils::absoluteFuzzyEqual((*_springs_elasticity_tensors[i])[_qp].L2norm(), 0.0))
     207           0 :       (*_springs_elasticity_tensors_inv[i])[_qp].zero();
     208             :     else
     209           0 :       (*_springs_elasticity_tensors_inv[i])[_qp] = (*_springs_elasticity_tensors[i])[_qp].invSymm();
     210             :   }
     211           0 : }
     212             : 
     213             : Real
     214       17232 : LinearViscoelasticityBase::computeTheta(Real dt, Real viscosity) const
     215             : {
     216       17232 :   if (MooseUtils::absoluteFuzzyEqual(dt, 0.0))
     217           0 :     mooseError("linear viscoelasticity cannot be integrated over a dt of ", dt);
     218             : 
     219       17232 :   switch (_integration_rule)
     220             :   {
     221             :     case IntegrationRule::BackwardEuler:
     222             :       return 1.;
     223           0 :     case IntegrationRule::MidPoint:
     224           0 :       return 0.5;
     225           0 :     case IntegrationRule::Newmark:
     226           0 :       return _theta;
     227           0 :     case IntegrationRule::Zienkiewicz:
     228           0 :       return 1. / (1. - std::exp(-dt / viscosity)) - viscosity / dt;
     229             :     default:
     230             :       return 1.;
     231             :   }
     232             :   return 1.;
     233             : }