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 "LinearViscoelasticityBase.h"
      11             : #include "Conversion.h"
      12             : 
      13             : InputParameters
      14         192 : LinearViscoelasticityBase::validParams()
      15             : {
      16         384 :   MooseEnum integration("backward-euler mid-point newmark zienkiewicz", "backward-euler");
      17             : 
      18         192 :   InputParameters params = ComputeElasticityTensorBase::validParams();
      19         384 :   params.addParam<MooseEnum>("integration_rule",
      20             :                              integration,
      21             :                              "describes how the viscoelastic behavior is integrated through time");
      22         384 :   params.addRangeCheckedParam<Real>("theta",
      23             :                                     1,
      24             :                                     "theta > 0 & theta <= 1",
      25             :                                     "coefficient for Newmark integration rule (between 0 and 1)");
      26         384 :   params.addParam<std::string>("driving_eigenstrain",
      27             :                                "name of the eigenstrain that increases the creep strains");
      28         384 :   params.addParam<std::string>(
      29             :       "elastic_strain_name", "elastic_strain", "name of the true elastic strain of the material");
      30         384 :   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         384 :   params.addParam<bool>("force_recompute_properties",
      36         384 :                         false,
      37             :                         "forces the computation of the viscoelastic properties at each step of"
      38             :                         "the solver (default: false)");
      39         384 :   params.addParam<bool>(
      40             :       "need_viscoelastic_properties_inverse",
      41         384 :       false,
      42             :       "checks whether the model requires the computation of the inverse viscoelastic"
      43             :       "properties (default: false)");
      44         192 :   params.suppressParameter<FunctionName>("elasticity_tensor_prefactor");
      45         192 :   return params;
      46         192 : }
      47             : 
      48         144 : LinearViscoelasticityBase::LinearViscoelasticityBase(const InputParameters & parameters)
      49             :   : ComputeElasticityTensorBase(parameters),
      50         144 :     _integration_rule(getParam<MooseEnum>("integration_rule").getEnum<IntegrationRule>()),
      51         288 :     _theta(getParam<Real>("theta")),
      52         144 :     _apparent_elasticity_tensor(
      53         144 :         declareProperty<RankFourTensor>(_base_name + "apparent_elasticity_tensor")),
      54         144 :     _apparent_elasticity_tensor_inv(
      55         144 :         declareProperty<RankFourTensor>(_base_name + "apparent_elasticity_tensor_inv")),
      56         144 :     _elasticity_tensor_inv(declareProperty<RankFourTensor>(_elasticity_tensor_name + "_inv")),
      57         288 :     _need_viscoelastic_properties_inverse(getParam<bool>("need_viscoelastic_properties_inverse")),
      58         144 :     _has_longterm_dashpot(false),
      59         144 :     _components(0),
      60         144 :     _first_elasticity_tensor(
      61         144 :         declareProperty<RankFourTensor>(_base_name + "spring_elasticity_tensor_0")),
      62         144 :     _first_elasticity_tensor_inv(
      63         144 :         _need_viscoelastic_properties_inverse
      64         252 :             ? &declareProperty<RankFourTensor>(_base_name + "spring_elasticity_tensor_0_inv")
      65             :             : nullptr),
      66         144 :     _apparent_creep_strain(declareProperty<RankTwoTensor>(_base_name + "apparent_creep_strain")),
      67         144 :     _apparent_creep_strain_old(
      68         144 :         getMaterialPropertyOld<RankTwoTensor>(_base_name + "apparent_creep_strain")),
      69         144 :     _elastic_strain_old(
      70         144 :         getMaterialPropertyOld<RankTwoTensor>(getParam<std::string>("elastic_strain_name"))),
      71         144 :     _creep_strain_old(
      72         144 :         getMaterialPropertyOld<RankTwoTensor>(getParam<std::string>("creep_strain_name"))),
      73         288 :     _has_driving_eigenstrain(isParamValid("driving_eigenstrain")),
      74         180 :     _driving_eigenstrain_name(
      75         144 :         _has_driving_eigenstrain ? getParam<std::string>("driving_eigenstrain") : ""),
      76         288 :     _driving_eigenstrain(_has_driving_eigenstrain
      77         144 :                              ? &getMaterialPropertyByName<RankTwoTensor>(_driving_eigenstrain_name)
      78             :                              : nullptr),
      79         288 :     _driving_eigenstrain_old(_has_driving_eigenstrain
      80         144 :                                  ? &getMaterialPropertyOld<RankTwoTensor>(_driving_eigenstrain_name)
      81             :                                  : nullptr),
      82         288 :     _force_recompute_properties(getParam<bool>("force_recompute_properties")),
      83         432 :     _step_zero(declareRestartableData<bool>("step_zero", true))
      84             : {
      85         144 :   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         144 :   getMaterialPropertyOld<RankFourTensor>(_base_name + "apparent_elasticity_tensor");
      91         288 :   getMaterialPropertyOld<RankFourTensor>(_base_name + "apparent_elasticity_tensor_inv");
      92             :   getMaterialPropertyOld<RankFourTensor>(_elasticity_tensor_name);
      93         144 :   getMaterialPropertyOld<RankFourTensor>(_elasticity_tensor_name + "_inv");
      94         144 :   getMaterialPropertyOld<RankFourTensor>(_base_name + "spring_elasticity_tensor_0");
      95         144 :   if (_need_viscoelastic_properties_inverse)
      96         216 :     getMaterialPropertyOld<RankFourTensor>(_base_name + "spring_elasticity_tensor_0_inv");
      97         144 : }
      98             : 
      99             : void
     100         144 : LinearViscoelasticityBase::declareViscoelasticProperties()
     101             : {
     102         396 :   for (unsigned int i = 0; i < _components; ++i)
     103             :   {
     104         252 :     std::string ith = Moose::stringify(i + 1);
     105             : 
     106         252 :     if (!_has_longterm_dashpot || (_components > 0 && i < _components - 1))
     107             :     {
     108         216 :       _springs_elasticity_tensors.push_back(
     109         216 :           &declareProperty<RankFourTensor>(_base_name + "spring_elasticity_tensor_" + ith));
     110         432 :       getMaterialPropertyOld<RankFourTensor>(_base_name + "spring_elasticity_tensor_" + ith);
     111             :     }
     112             : 
     113         252 :     _dashpot_viscosities.push_back(&declareProperty<Real>(_base_name + "dashpot_viscosity_" + ith));
     114         252 :     _dashpot_viscosities_old.push_back(
     115         504 :         &getMaterialPropertyOld<Real>(_base_name + "dashpot_viscosity_" + ith));
     116             : 
     117         252 :     _viscous_strains.push_back(
     118         252 :         &declareProperty<RankTwoTensor>(_base_name + "viscous_strain_" + ith));
     119         252 :     _viscous_strains_old.push_back(
     120         504 :         &getMaterialPropertyOld<RankTwoTensor>(_base_name + "viscous_strain_" + ith));
     121             : 
     122         252 :     if (_need_viscoelastic_properties_inverse)
     123             :     {
     124         198 :       _springs_elasticity_tensors_inv.push_back(&declareProperty<RankFourTensor>(
     125         198 :           _base_name + "spring_elasticity_tensor_" + ith + "_inv"));
     126         198 :       _springs_elasticity_tensors_inv_old.push_back(&getMaterialPropertyOld<RankFourTensor>(
     127         396 :           _base_name + "spring_elasticity_tensor_" + ith + "_inv"));
     128             :     }
     129             :   }
     130         144 : }
     131             : 
     132             : void
     133         640 : LinearViscoelasticityBase::initQpStatefulProperties()
     134             : {
     135         640 :   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         640 :   _apparent_creep_strain[_qp].zero();
     141         640 :   _apparent_elasticity_tensor[_qp].zero();
     142         640 :   _apparent_elasticity_tensor_inv[_qp].zero();
     143         640 :   _elasticity_tensor_inv[_qp].zero();
     144         640 :   _first_elasticity_tensor[_qp].zero();
     145         640 :   if (_need_viscoelastic_properties_inverse)
     146         512 :     (*_first_elasticity_tensor_inv)[_qp].zero();
     147             : 
     148        1792 :   for (unsigned int i = 0; i < _components; ++i)
     149             :   {
     150        1152 :     if (!_has_longterm_dashpot || (_components > 0 && i < _components - 1))
     151             :     {
     152         992 :       (*_springs_elasticity_tensors[i])[_qp].zero();
     153         992 :       if (_need_viscoelastic_properties_inverse)
     154         800 :         (*_springs_elasticity_tensors_inv[i])[_qp].zero();
     155             :     }
     156             : 
     157        1152 :     (*_dashpot_viscosities[i])[_qp] = 0.0;
     158        1152 :     (*_viscous_strains[i])[_qp].zero();
     159             :   }
     160         640 : }
     161             : 
     162             : void
     163        7296 : LinearViscoelasticityBase::recomputeQpApparentProperties(unsigned int qp)
     164             : {
     165        7296 :   unsigned int qp_prev = _qp;
     166        7296 :   _qp = qp;
     167             : 
     168        7296 :   if (_t_step >= 1)
     169        7296 :     _step_zero = false;
     170             : 
     171             :   // 1. we get the viscoelastic properties and their inverse if needed
     172        7296 :   computeQpViscoelasticProperties();
     173        7296 :   if (_need_viscoelastic_properties_inverse)
     174        5312 :     computeQpViscoelasticPropertiesInv();
     175             : 
     176             :   // 2. we update the internal viscous strains from the previous time step
     177        7296 :   updateQpViscousStrains();
     178             : 
     179             :   // 3. we compute the apparent elasticity tensor
     180        7296 :   computeQpApparentElasticityTensors();
     181             : 
     182             :   // 4. we transform the internal viscous strains in an apparent creep strain
     183        7296 :   if (!_step_zero)
     184        7296 :     computeQpApparentCreepStrain();
     185             : 
     186        7296 :   _qp = qp_prev;
     187        7296 : }
     188             : 
     189             : void
     190      112496 : LinearViscoelasticityBase::computeQpElasticityTensor()
     191             : {
     192      112496 :   if (_force_recompute_properties)
     193           0 :     recomputeQpApparentProperties(_qp);
     194      112496 : }
     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       34464 : LinearViscoelasticityBase::computeTheta(Real dt, Real viscosity) const
     215             : {
     216       34464 :   if (MooseUtils::absoluteFuzzyEqual(dt, 0.0))
     217           0 :     mooseError("linear viscoelasticity cannot be integrated over a dt of ", dt);
     218             : 
     219       34464 :   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             : }