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 "ADViscoplasticityStressUpdate.h"
11 :
12 : #include "libmesh/utility.h"
13 :
14 : registerMooseObject("SolidMechanicsApp", ADViscoplasticityStressUpdate);
15 :
16 : InputParameters
17 1104 : ADViscoplasticityStressUpdate::validParams()
18 : {
19 1104 : InputParameters params = ADViscoplasticityStressUpdateBase::validParams();
20 1104 : params += ADSingleVariableReturnMappingSolution::validParams();
21 1104 : params.addClassDescription(
22 : "This material computes the non-linear homogenized gauge stress in order to compute the "
23 : "viscoplastic responce due to creep in porous materials. This material must be used in "
24 : "conjunction with ADComputeMultiplePorousInelasticStress");
25 2208 : MooseEnum viscoplasticity_model("LPS GTN", "LPS");
26 2208 : params.addParam<MooseEnum>(
27 : "viscoplasticity_model", viscoplasticity_model, "Which viscoplastic model to use");
28 2208 : MooseEnum pore_shape_model("spherical cylindrical", "spherical");
29 2208 : params.addParam<MooseEnum>("pore_shape_model", pore_shape_model, "Which pore shape model to use");
30 2208 : params.addRequiredParam<MaterialPropertyName>(
31 : "coefficient", "Material property name for the leading coefficient for Norton power law");
32 2208 : params.addRequiredRangeCheckedParam<Real>(
33 : "power", "power>=1.0", "Stress exponent for Norton power law");
34 2208 : params.addParam<Real>(
35 : "maximum_gauge_ratio",
36 2208 : 1.0e6,
37 : "Maximum ratio between the gauge stress and the equivalent stress. This "
38 : "should be a high number. Note that this does not set an upper bound on the value, but "
39 : "rather will help with convergence of the inner Newton loop");
40 2208 : params.addParam<Real>(
41 : "minimum_equivalent_stress",
42 2208 : 1.0e-3,
43 : "Minimum value of equivalent stress below which viscoplasticiy is not calculated.");
44 2208 : params.addParam<Real>("maximum_equivalent_stress",
45 2208 : 1.0e12,
46 : "Maximum value of equivalent stress above which an exception is thrown "
47 : "instead of calculating the properties in this material.");
48 :
49 2208 : params.addParamNamesToGroup("verbose maximum_gauge_ratio maximum_equivalent_stress", "Advanced");
50 1104 : return params;
51 1104 : }
52 :
53 828 : ADViscoplasticityStressUpdate::ADViscoplasticityStressUpdate(const InputParameters & parameters)
54 : : ADViscoplasticityStressUpdateBase(parameters),
55 : ADSingleVariableReturnMappingSolution(parameters),
56 828 : _model(parameters.get<MooseEnum>("viscoplasticity_model").getEnum<ViscoplasticityModel>()),
57 828 : _pore_shape(parameters.get<MooseEnum>("pore_shape_model").getEnum<PoreShapeModel>()),
58 828 : _pore_shape_factor(_pore_shape == PoreShapeModel::SPHERICAL ? 1.5 : std::sqrt(3.0)),
59 1656 : _power(getParam<Real>("power")),
60 828 : _power_factor(_model == ViscoplasticityModel::LPS ? (_power - 1.0) / (_power + 1.0) : 1.0),
61 1656 : _coefficient(getADMaterialProperty<Real>("coefficient")),
62 828 : _gauge_stress(declareADProperty<Real>(_base_name + "gauge_stress")),
63 1656 : _maximum_gauge_ratio(getParam<Real>("maximum_gauge_ratio")),
64 1656 : _minimum_equivalent_stress(getParam<Real>("minimum_equivalent_stress")),
65 1656 : _maximum_equivalent_stress(getParam<Real>("maximum_equivalent_stress")),
66 828 : _hydro_stress(0.0),
67 828 : _identity_two(RankTwoTensor::initIdentity),
68 828 : _dhydro_stress_dsigma(_identity_two / 3.0),
69 1656 : _derivative(0.0)
70 : {
71 828 : _check_range = true;
72 828 : }
73 :
74 : void
75 244226 : ADViscoplasticityStressUpdate::updateState(ADRankTwoTensor & elastic_strain_increment,
76 : ADRankTwoTensor & inelastic_strain_increment,
77 : const ADRankTwoTensor & /*rotation_increment*/,
78 : ADRankTwoTensor & stress,
79 : const RankTwoTensor & /*stress_old*/,
80 : const ADRankFourTensor & elasticity_tensor,
81 : const RankTwoTensor & elastic_strain_old,
82 : bool /*compute_full_tangent_operator = false*/,
83 : RankFourTensor & /*tangent_operator = _identityTensor*/)
84 : {
85 : using std::sqrt;
86 :
87 : // Compute initial hydrostatic stress and porosity
88 244226 : if (_pore_shape == PoreShapeModel::CYLINDRICAL)
89 168448 : _hydro_stress = (stress(0, 0) + stress(1, 1)) / 2.0;
90 : else
91 320004 : _hydro_stress = stress.trace() / 3.0;
92 :
93 244226 : updateIntermediatePorosity(elastic_strain_increment);
94 :
95 : // Compute intermediate equivalent stress
96 244224 : const ADRankTwoTensor dev_stress = stress.deviatoric();
97 244224 : const ADReal dev_stress_squared = dev_stress.doubleContraction(dev_stress);
98 484272 : const ADReal equiv_stress = dev_stress_squared == 0.0 ? 0.0 : sqrt(1.5 * dev_stress_squared);
99 :
100 244224 : computeStressInitialize(equiv_stress, elasticity_tensor);
101 :
102 : // Prepare values
103 244224 : _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp];
104 244224 : _inelastic_strain[_qp] = _inelastic_strain_old[_qp];
105 244224 : inelastic_strain_increment.zero();
106 :
107 : // Protect against extremely high values of stresses calculated by other viscoplastic materials
108 244224 : if (equiv_stress > _maximum_equivalent_stress)
109 0 : mooseException("In ",
110 : _name,
111 : ": equivalent stress (",
112 : equiv_stress,
113 : ") is higher than maximum_equivalent_stress (",
114 : _maximum_equivalent_stress,
115 : ").\nCutting time step.");
116 :
117 : // If equivalent stress is present, calculate creep strain increment
118 244224 : if (equiv_stress > _minimum_equivalent_stress)
119 : {
120 : // Initalize stress potential
121 240048 : ADReal dpsi_dgauge(0);
122 :
123 240048 : computeInelasticStrainIncrement(_gauge_stress[_qp],
124 : dpsi_dgauge,
125 : inelastic_strain_increment,
126 : equiv_stress,
127 : dev_stress,
128 : stress);
129 :
130 : // Update elastic strain increment due to inelastic strain calculated here
131 240048 : elastic_strain_increment -= inelastic_strain_increment;
132 : // Update stress due to new strain
133 240048 : stress = elasticity_tensor * (elastic_strain_old + elastic_strain_increment);
134 :
135 : // Compute effective strain from the stress potential. Note that this is approximate and to be
136 : // used qualitatively
137 480096 : _effective_inelastic_strain[_qp] += dpsi_dgauge * _dt;
138 : // Update creep strain due to currently computed inelastic strain
139 240048 : _inelastic_strain[_qp] += inelastic_strain_increment;
140 : }
141 :
142 244224 : const ADRankTwoTensor new_dev_stress = stress.deviatoric();
143 244224 : const ADReal new_dev_stress_squared = new_dev_stress.doubleContraction(new_dev_stress);
144 : const ADReal new_equiv_stress =
145 724320 : new_dev_stress_squared == 0.0 ? 0.0 : sqrt(1.5 * new_dev_stress_squared);
146 :
147 244224 : if (MooseUtils::relativeFuzzyGreaterThan(new_equiv_stress, equiv_stress))
148 0 : mooseException("In ",
149 : _name,
150 : ": updated equivalent stress (",
151 : MetaPhysicL::raw_value(new_equiv_stress),
152 : ") is greater than initial equivalent stress (",
153 : MetaPhysicL::raw_value(equiv_stress),
154 : "). Try decreasing max_inelastic_increment to avoid this exception.");
155 :
156 244224 : computeStressFinalize(inelastic_strain_increment);
157 244224 : }
158 :
159 : ADReal
160 238096 : ADViscoplasticityStressUpdate::initialGuess(const ADReal & effective_trial_stress)
161 : {
162 238096 : return effective_trial_stress;
163 : }
164 :
165 : ADReal
166 238096 : ADViscoplasticityStressUpdate::maximumPermissibleValue(const ADReal & effective_trial_stress) const
167 : {
168 238096 : return effective_trial_stress * _maximum_gauge_ratio;
169 : }
170 :
171 : ADReal
172 238096 : ADViscoplasticityStressUpdate::minimumPermissibleValue(const ADReal & effective_trial_stress) const
173 : {
174 238096 : return effective_trial_stress;
175 : }
176 :
177 : ADReal
178 1900886 : ADViscoplasticityStressUpdate::computeResidual(const ADReal & equiv_stress,
179 : const ADReal & trial_gauge)
180 : {
181 : using std::abs, std::cosh, std::sinh;
182 :
183 1900886 : const ADReal M = abs(_hydro_stress) / trial_gauge;
184 3801772 : const ADReal dM_dtrial_gauge = -M / trial_gauge;
185 :
186 1900886 : const ADReal residual_left = Utility::pow<2>(equiv_stress / trial_gauge);
187 1900886 : const ADReal dresidual_left_dtrial_gauge = -2.0 * residual_left / trial_gauge;
188 :
189 1900886 : ADReal residual = residual_left;
190 1900886 : _derivative = dresidual_left_dtrial_gauge;
191 :
192 1900886 : if (_pore_shape == PoreShapeModel::SPHERICAL)
193 : {
194 3203307 : residual *= 1.0 + _intermediate_porosity / 1.5;
195 3203307 : _derivative *= 1.0 + _intermediate_porosity / 1.5;
196 : }
197 :
198 1900886 : if (_model == ViscoplasticityModel::GTN)
199 : {
200 1234824 : residual += 2.0 * _intermediate_porosity * cosh(_pore_shape_factor * M) - 1.0 -
201 1234824 : Utility::pow<2>(_intermediate_porosity);
202 617412 : _derivative += 2.0 * _intermediate_porosity * sinh(_pore_shape_factor * M) *
203 617412 : _pore_shape_factor * dM_dtrial_gauge;
204 : }
205 : else
206 : {
207 1283474 : const ADReal h = computeH(_power, M);
208 1283474 : const ADReal dh_dM = computeH(_power, M, true);
209 :
210 1283474 : residual += _intermediate_porosity * (h + _power_factor / h) - 1.0 -
211 2566948 : _power_factor * Utility::pow<2>(_intermediate_porosity);
212 2566948 : const ADReal dresidual_dh = _intermediate_porosity * (1.0 - _power_factor / Utility::pow<2>(h));
213 1283474 : _derivative += dresidual_dh * dh_dM * dM_dtrial_gauge;
214 : }
215 :
216 1900886 : if (_verbose)
217 : {
218 : Moose::out << "in computeResidual:\n"
219 0 : << " position: " << _q_point[_qp] << " hydro_stress: " << _hydro_stress
220 : << " equiv_stress: " << equiv_stress << " trial_grage: " << trial_gauge
221 : << " M: " << M << std::endl;
222 : Moose::out << " residual: " << residual << " derivative: " << _derivative << std::endl;
223 : }
224 :
225 1900886 : return residual;
226 : }
227 :
228 : ADReal
229 3016212 : ADViscoplasticityStressUpdate::computeH(const Real n, const ADReal & M, const bool derivative)
230 : {
231 : using std::pow;
232 :
233 6032424 : const ADReal mod = pow(M * _pore_shape_factor, (n + 1.0) / n);
234 :
235 : // Calculate derivative with respect to M
236 3016212 : if (derivative)
237 : {
238 1492690 : const ADReal dmod_dM = (n + 1.0) / n * mod / M;
239 5970760 : return dmod_dM * pow(1.0 + mod / n, n - 1.0);
240 : }
241 :
242 3047044 : return pow(1.0 + mod / n, n);
243 : }
244 :
245 : ADRankTwoTensor
246 240048 : ADViscoplasticityStressUpdate::computeDGaugeDSigma(const ADReal & gauge_stress,
247 : const ADReal & equiv_stress,
248 : const ADRankTwoTensor & dev_stress,
249 : const ADRankTwoTensor & stress)
250 : {
251 : using std::abs, std::sinh;
252 :
253 : // Compute the derivative of the gauge stress with respect to the equivalent and hydrostatic
254 : // stress components
255 240048 : const ADReal M = abs(_hydro_stress) / gauge_stress;
256 240048 : const ADReal h = computeH(_power, M);
257 :
258 : // Compute the derviative of the residual with respect to the hydrostatic stress
259 240048 : ADReal dresidual_dhydro_stress = 0.0;
260 240048 : if (_hydro_stress != 0.0)
261 : {
262 : const ADReal dM_dhydro_stress = M / _hydro_stress;
263 240048 : if (_model == ViscoplasticityModel::GTN)
264 : {
265 30832 : dresidual_dhydro_stress = 2.0 * _intermediate_porosity * sinh(_pore_shape_factor * M) *
266 61664 : _pore_shape_factor * dM_dhydro_stress;
267 : }
268 : else
269 : {
270 : const ADReal dresidual_dh =
271 627648 : _intermediate_porosity * (1.0 - _power_factor / Utility::pow<2>(h));
272 209216 : const ADReal dh_dM = computeH(_power, M, true);
273 209216 : dresidual_dhydro_stress = dresidual_dh * dh_dM * dM_dhydro_stress;
274 : }
275 : }
276 :
277 : // Compute the derivative of the residual with respect to the equivalent stress
278 240048 : ADReal dresidual_dequiv_stress = 2.0 * equiv_stress / Utility::pow<2>(gauge_stress);
279 240048 : if (_pore_shape == PoreShapeModel::SPHERICAL)
280 473352 : dresidual_dequiv_stress *= 1.0 + _intermediate_porosity / 1.5;
281 :
282 : // Compute the derivative of the equilvalent stress to the deviatoric stress
283 480096 : const ADRankTwoTensor dequiv_stress_dsigma = 1.5 * dev_stress / equiv_stress;
284 :
285 : // Compute the derivative of the residual with the stress
286 240048 : const ADRankTwoTensor dresidual_dsigma = dresidual_dhydro_stress * _dhydro_stress_dsigma +
287 240048 : dresidual_dequiv_stress * dequiv_stress_dsigma;
288 :
289 : // Compute the deritative of the gauge stress with respect to the stress
290 : const ADRankTwoTensor dgauge_dsigma =
291 480096 : dresidual_dsigma * (gauge_stress / dresidual_dsigma.doubleContraction(stress));
292 :
293 240048 : return dgauge_dsigma;
294 : }
295 :
296 : void
297 240048 : ADViscoplasticityStressUpdate::computeInelasticStrainIncrement(
298 : ADReal & gauge_stress,
299 : ADReal & dpsi_dgauge,
300 : ADRankTwoTensor & inelastic_strain_increment,
301 : const ADReal & equiv_stress,
302 : const ADRankTwoTensor & dev_stress,
303 : const ADRankTwoTensor & stress)
304 : {
305 : using std::sqrt, std::pow;
306 :
307 : // If hydrostatic stress and porosity present, compute non-linear gauge stress
308 240048 : if (_intermediate_porosity == 0.0)
309 1952 : gauge_stress = equiv_stress;
310 238096 : else if (_hydro_stress == 0.0)
311 0 : gauge_stress = equiv_stress / sqrt(1.0 - (1.0 + _power_factor) * _intermediate_porosity +
312 0 : _power_factor * Utility::pow<2>(_intermediate_porosity));
313 : else
314 238096 : returnMappingSolve(equiv_stress, gauge_stress, _console);
315 :
316 : mooseAssert(gauge_stress >= equiv_stress,
317 : "Gauge stress calculated in inner Newton solve is less than the equivalent stress.");
318 :
319 : // Compute stress potential
320 480096 : dpsi_dgauge = _coefficient[_qp] * pow(gauge_stress, _power);
321 :
322 : // Compute strain increment from stress potential and the gauge stress derivative with respect
323 : // to the stress stress. The current form is explicit, and should eventually be changed
324 : inelastic_strain_increment =
325 240048 : _dt * dpsi_dgauge * computeDGaugeDSigma(gauge_stress, equiv_stress, dev_stress, stress);
326 240048 : }
327 :
328 : void
329 0 : ADViscoplasticityStressUpdate::outputIterationSummary(std::stringstream * iter_output,
330 : const unsigned int total_it)
331 : {
332 0 : if (iter_output)
333 : {
334 0 : *iter_output << "At element " << _current_elem->id() << " _qp=" << _qp << " Coordinates "
335 0 : << _q_point[_qp] << " block=" << _current_elem->subdomain_id() << '\n';
336 : }
337 0 : ADSingleVariableReturnMappingSolution::outputIterationSummary(iter_output, total_it);
338 0 : }
339 :
340 : Real
341 1900886 : ADViscoplasticityStressUpdate::computeReferenceResidual(const ADReal & /*effective_trial_stress*/,
342 : const ADReal & gauge_stress)
343 : {
344 : // Use gauge stress for relative tolerance criteria, defined as:
345 : // abs(residual / gauge_stress) <= _relative_tolerance
346 1900886 : return MetaPhysicL::raw_value(gauge_stress);
347 : }
|
