www.mooseframework.org
LinearViscoelasticityBase.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
8 #include "libmesh/quadrature.h"
9 
10 template <>
11 InputParameters
13 {
14  MooseEnum integration("backward-euler mid-point newmark zienkiewicz", "backward-euler");
15 
16  InputParameters params = validParams<ComputeElasticityTensorBase>();
17  params.addParam<MooseEnum>("integration_rule",
18  integration,
19  "describes how the viscoelastic behavior is integrated through time");
20  params.addRangeCheckedParam<Real>("theta",
21  1,
22  "theta > 0 & theta <= 1",
23  "coefficient for newmark integration rule (between 0 and 1)");
24  params.addParam<std::string>("driving_eigenstrain",
25  "name of the eigenstrain that increases the creep strains");
26  params.addParam<std::string>(
27  "elastic_strain_name", "elastic_strain", "name of the true elastic strain of the material");
28  params.addParam<std::string>("stress_name", "stress", "name of the true stress of the material");
29  params.addParam<std::string>("creep_strain_name",
30  "creep_strain",
31  "name of the true creep strain of the material"
32  "(computed by LinearViscoelasticStressUpdate or"
33  "ComputeLinearViscoelasticStress");
34  params.addParam<bool>("force_recompute_properties",
35  false,
36  "forces the computation of the viscoelastic properties at each step of"
37  "the solver (default: false)");
38  params.addParam<bool>(
39  "need_viscoelastic_properties_inverse",
40  false,
41  "checks whether the model requires the computation of the inverse viscoelastic"
42  "properties (default: false)");
43  params.suppressParameter<FunctionName>("elasticity_tensor_prefactor");
44  return params;
45 }
46 
47 LinearViscoelasticityBase::LinearViscoelasticityBase(const InputParameters & parameters)
48  : ComputeElasticityTensorBase(parameters),
49  _integration_rule(getParam<MooseEnum>("integration_rule").getEnum<IntegrationRule>()),
50  _theta(getParam<Real>("theta")),
51  _apparent_elasticity_tensor(
52  declareProperty<RankFourTensor>(_base_name + "apparent_elasticity_tensor")),
53  _apparent_elasticity_tensor_inv(
54  declareProperty<RankFourTensor>(_base_name + "apparent_elasticity_tensor_inv")),
55  _instantaneous_elasticity_tensor(
56  declareProperty<RankFourTensor>(_base_name + "instantaneous_elasticity_tensor")),
57  _instantaneous_elasticity_tensor_inv(
58  declareProperty<RankFourTensor>(_base_name + "instantaneous_elasticity_tensor_inv")),
59  _need_viscoelastic_properties_inverse(getParam<bool>("need_viscoelastic_properties_inverse")),
60  _has_longterm_dashpot(false),
61  _components(0),
62  _first_elasticity_tensor(
63  declareProperty<RankFourTensor>(_base_name + "first_elasticity_tensor")),
64  _first_elasticity_tensor_inv(
65  _need_viscoelastic_properties_inverse
66  ? &declareProperty<RankFourTensor>(_base_name + "first_elasticity_tensor_inv")
67  : nullptr),
68  _springs_elasticity_tensors(
69  declareProperty<std::vector<RankFourTensor>>(_base_name + "springs_elasticity_tensors")),
70  _springs_elasticity_tensors_inv(_need_viscoelastic_properties_inverse
71  ? &declareProperty<std::vector<RankFourTensor>>(
72  _base_name + "springs_elasticity_tensors_inv")
73  : nullptr),
74  _dashpot_viscosities(declareProperty<std::vector<Real>>(_base_name + "dashpot_viscosities")),
75  _viscous_strains(declareProperty<std::vector<RankTwoTensor>>(_base_name + "viscous_strains")),
76  _viscous_strains_old(
77  getMaterialPropertyOld<std::vector<RankTwoTensor>>(_base_name + "viscous_strains")),
78  _apparent_creep_strain(declareProperty<RankTwoTensor>(_base_name + "apparent_creep_strain")),
79  _apparent_creep_strain_old(
80  getMaterialPropertyOld<RankTwoTensor>(_base_name + "apparent_creep_strain")),
81  _elastic_strain_old(
82  getMaterialPropertyOld<RankTwoTensor>(getParam<std::string>("elastic_strain_name"))),
83  _creep_strain_old(
84  getMaterialPropertyOld<RankTwoTensor>(getParam<std::string>("creep_strain_name"))),
85  _stress_old(getMaterialPropertyOld<RankTwoTensor>(getParam<std::string>("stress_name"))),
86  _has_driving_eigenstrain(isParamValid("driving_eigenstrain")),
87  _driving_eigenstrain_name(
88  _has_driving_eigenstrain ? getParam<std::string>("driving_eigenstrain") : ""),
89  _driving_eigenstrain(_has_driving_eigenstrain
90  ? &getMaterialPropertyByName<RankTwoTensor>(_driving_eigenstrain_name)
91  : nullptr),
92  _force_recompute_properties(getParam<bool>("force_recompute_properties")),
93  _step_zero(declareRestartableData<bool>("step_zero", true))
94 {
95  if (_theta < 0.5)
96  mooseWarning("theta parameter for LinearViscoelasticityBase is below 0.5; time integration may "
97  "not converge!");
98 
99  // force material properties to be considered stateful
100  getMaterialPropertyOld<RankFourTensor>(_base_name + "apparent_elasticity_tensor");
101  getMaterialPropertyOld<RankFourTensor>(_base_name + "apparent_elasticity_tensor_inv");
102  getMaterialPropertyOld<RankFourTensor>(_base_name + "instantaneous_elasticity_tensor");
103  getMaterialPropertyOld<RankFourTensor>(_base_name + "instantaneous_elasticity_tensor_inv");
104  getMaterialPropertyOld<RankFourTensor>(_base_name + "first_elasticity_tensor");
105  getMaterialPropertyOld<std::vector<RankFourTensor>>(_base_name + "springs_elasticity_tensors");
106  getMaterialPropertyOld<std::vector<Real>>(_base_name + "dashpot_viscosities");
107 
109  {
110  getMaterialPropertyOld<RankFourTensor>(_base_name + "first_elasticity_tensor_inv");
111  getMaterialPropertyOld<std::vector<RankFourTensor>>(_base_name +
112  "springs_elasticity_tensors_inv");
113  }
114 }
115 
116 void
118 {
119  _apparent_elasticity_tensor[_qp].zero();
123  _first_elasticity_tensor[_qp].zero();
124  _apparent_creep_strain[_qp].zero();
125 
126  _springs_elasticity_tensors[_qp].resize(
127  (_has_longterm_dashpot && _components > 0) ? _components - 1 : _components, RankFourTensor());
128  _dashpot_viscosities[_qp].resize(_components, 0);
129  _viscous_strains[_qp].resize(_components, RankTwoTensor());
130 
132  {
133  (*_first_elasticity_tensor_inv)[_qp].zero();
134  (*_springs_elasticity_tensors_inv)[_qp].resize(_components, RankFourTensor());
135  }
136 }
137 
138 void
140 {
141  unsigned int qp_prev = _qp;
142  _qp = qp;
143 
144  if (_t_step >= 1)
145  _step_zero = false;
146 
147  // 1. we get the viscoelastic properties and their inverse if needed
151 
152  // 2. we update the internal viscous strains from the previous time step
154 
155  // 3. we compute the apparent elasticity tensor
157 
158  // 4. we transform the internal viscous strains in an apparent creep strain
159  if (!_step_zero)
161 
162  _qp = qp_prev;
163 }
164 
165 void
167 {
170 
172 }
173 
174 void
176 {
177  if (MooseUtils::absoluteFuzzyEqual(_first_elasticity_tensor[_qp].L2norm(), 0.0))
178  (*_first_elasticity_tensor_inv)[_qp].zero();
179  else
181 
182  for (unsigned int i = 0; i < _springs_elasticity_tensors[_qp].size(); ++i)
183  {
184  if (MooseUtils::absoluteFuzzyEqual(_springs_elasticity_tensors[_qp][i].L2norm(), 0.0))
185  (*_springs_elasticity_tensors_inv)[_qp][i].zero();
186  else
187  (*_springs_elasticity_tensors_inv)[_qp][i] = _springs_elasticity_tensors[_qp][i].invSymm();
188  }
189 }
190 
191 Real
193 {
194  if (MooseUtils::absoluteFuzzyEqual(dt, 0.0))
195  mooseError("linear viscoelasticity cannot be integrated over a dt of ", dt);
196 
197  switch (_integration_rule)
198  {
200  return 1.;
202  return 0.5;
204  return _theta;
206  return 1. / (1. - std::exp(-dt / viscosity)) - viscosity / dt;
207  default:
208  return 1.;
209  }
210  return 1.;
211 }
ComputeElasticityTensorBase the base class for computing elasticity tensors.
LinearViscoelasticityBase(const InputParameters &parameters)
void recomputeQpApparentProperties(unsigned int qp)
Compute the apparent properties at a quadrature point.
MaterialProperty< std::vector< RankTwoTensor > > & _viscous_strains
The internal strain variables required by the time-stepping procedure (generally, on a one-on-one bas...
virtual void initQpStatefulProperties() override
bool _has_longterm_dashpot
Indicates if the spring-dashpot assembly has a single dashpot not associated with a spring...
bool _force_recompute_properties
If activated, the time-stepping scheme will be re-initialized at each step of the solver...
MaterialProperty< std::vector< RankFourTensor > > & _springs_elasticity_tensors
List of elasticity tensor of each subsequent spring in the chain.
bool _need_viscoelastic_properties_inverse
If active, indicates that we need to call computeQpViscoelasticPropertiesInv()
MaterialProperty< RankFourTensor > & _instantaneous_elasticity_tensor
Instantaneous elasticity tensor. This IS the real elasticity tensor of the material.
theta automatically adjusted as a function of the time step and the viscosity
double FORTRAN_CALL() viscosity(double &rho, double &T)
MaterialProperty< std::vector< RankFourTensor > > * _springs_elasticity_tensors_inv
MaterialProperty< RankFourTensor > & _apparent_elasticity_tensor
Apparent elasticity tensor. This is NOT the elasticity tensor of the material.
IntegrationRule
Determines how theta is calculated for the time-integration system.
MaterialProperty< RankFourTensor > & _first_elasticity_tensor
Elasticity tensor of a stand-alone elastic spring in the chain.
MaterialProperty< RankFourTensor > * _first_elasticity_tensor_inv
MaterialProperty< std::vector< Real > > & _dashpot_viscosities
List of viscosities of each subsequent spring in the chain.
InputParameters validParams< LinearViscoelasticityBase >()
Real computeTheta(Real dt, Real viscosity) const
Provides theta as a function of the time step and a viscosity.
Real _theta
User-defined value for theta.
virtual void computeQpApparentElasticityTensors()=0
This method computes the apparent elasticity tensor used in the internal time-stepping scheme...
InputParameters validParams< ComputeElasticityTensorBase >()
virtual void computeQpViscoelasticProperties()=0
This method assigns the mechanical properties of each spring and dashpot in the system.
unsigned int _components
This is the number of internal variables required by the model.
virtual void computeQpElasticityTensor() final
Inherited from ComputeElasticityTensorBase.
MaterialProperty< RankTwoTensor > & _apparent_creep_strain
The apparent creep strain resulting from the internal viscous strains.
virtual void computeQpApparentCreepStrain()=0
This method computes the apparent creep strain corresponding to the current viscous_strain of each da...
MaterialProperty< RankFourTensor > & _instantaneous_elasticity_tensor_inv
Inverse of the instaneous elasticity tensor.
Real L2norm(const RankTwoTensor &r2tensor)
IntegrationRule _integration_rule
Determines how theta is computed.
virtual void computeQpViscoelasticPropertiesInv()
This method computes the inverse elasticity tensor of each spring in the system (if required)...
virtual void updateQpViscousStrains()=0
Update the internal viscous strains at a quadrature point.
MaterialProperty< RankFourTensor > & _elasticity_tensor
MaterialProperty< RankFourTensor > & _apparent_elasticity_tensor_inv
Inverse of the apparent elasticity tensor.