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