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 1016 : ADViscoplasticityStressUpdate::validParams()
18 : {
19 1016 : InputParameters params = ADViscoplasticityStressUpdateBase::validParams();
20 1016 : params += ADSingleVariableReturnMappingSolution::validParams();
21 1016 : 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 2032 : MooseEnum viscoplasticity_model("LPS GTN", "LPS");
26 2032 : params.addParam<MooseEnum>(
27 : "viscoplasticity_model", viscoplasticity_model, "Which viscoplastic model to use");
28 2032 : MooseEnum pore_shape_model("spherical cylindrical", "spherical");
29 2032 : params.addParam<MooseEnum>("pore_shape_model", pore_shape_model, "Which pore shape model to use");
30 2032 : params.addRequiredParam<MaterialPropertyName>(
31 : "coefficient", "Material property name for the leading coefficient for Norton power law");
32 2032 : params.addRequiredRangeCheckedParam<Real>(
33 : "power", "power>=1.0", "Stress exponent for Norton power law");
34 2032 : params.addParam<Real>(
35 : "maximum_gauge_ratio",
36 2032 : 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 2032 : params.addParam<Real>(
41 : "minimum_equivalent_stress",
42 2032 : 1.0e-3,
43 : "Minimum value of equivalent stress below which viscoplasticiy is not calculated.");
44 2032 : params.addParam<Real>("maximum_equivalent_stress",
45 2032 : 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 2032 : params.addParamNamesToGroup("verbose maximum_gauge_ratio maximum_equivalent_stress", "Advanced");
50 1016 : return params;
51 1016 : }
52 :
53 762 : ADViscoplasticityStressUpdate::ADViscoplasticityStressUpdate(const InputParameters & parameters)
54 : : ADViscoplasticityStressUpdateBase(parameters),
55 : ADSingleVariableReturnMappingSolution(parameters),
56 762 : _model(parameters.get<MooseEnum>("viscoplasticity_model").getEnum<ViscoplasticityModel>()),
57 762 : _pore_shape(parameters.get<MooseEnum>("pore_shape_model").getEnum<PoreShapeModel>()),
58 762 : _pore_shape_factor(_pore_shape == PoreShapeModel::SPHERICAL ? 1.5 : std::sqrt(3.0)),
59 1524 : _power(getParam<Real>("power")),
60 762 : _power_factor(_model == ViscoplasticityModel::LPS ? (_power - 1.0) / (_power + 1.0) : 1.0),
61 1524 : _coefficient(getADMaterialProperty<Real>("coefficient")),
62 762 : _gauge_stress(declareADProperty<Real>(_base_name + "gauge_stress")),
63 1524 : _maximum_gauge_ratio(getParam<Real>("maximum_gauge_ratio")),
64 1524 : _minimum_equivalent_stress(getParam<Real>("minimum_equivalent_stress")),
65 1524 : _maximum_equivalent_stress(getParam<Real>("maximum_equivalent_stress")),
66 762 : _hydro_stress(0.0),
67 762 : _identity_two(RankTwoTensor::initIdentity),
68 762 : _dhydro_stress_dsigma(_identity_two / 3.0),
69 1524 : _derivative(0.0)
70 : {
71 762 : _check_range = true;
72 762 : }
73 :
74 : void
75 446434 : 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 : // Compute initial hydrostatic stress and porosity
86 446434 : if (_pore_shape == PoreShapeModel::CYLINDRICAL)
87 377216 : _hydro_stress = (stress(0, 0) + stress(1, 1)) / 2.0;
88 : else
89 515652 : _hydro_stress = stress.trace() / 3.0;
90 :
91 446434 : updateIntermediatePorosity(elastic_strain_increment);
92 :
93 : // Compute intermediate equivalent stress
94 446432 : const ADRankTwoTensor dev_stress = stress.deviatoric();
95 446432 : const ADReal dev_stress_squared = dev_stress.doubleContraction(dev_stress);
96 885024 : const ADReal equiv_stress = dev_stress_squared == 0.0 ? 0.0 : std::sqrt(1.5 * dev_stress_squared);
97 :
98 446432 : computeStressInitialize(equiv_stress, elasticity_tensor);
99 :
100 : // Prepare values
101 446432 : _effective_inelastic_strain[_qp] = _effective_inelastic_strain_old[_qp];
102 446432 : _inelastic_strain[_qp] = _inelastic_strain_old[_qp];
103 446432 : inelastic_strain_increment.zero();
104 :
105 : // Protect against extremely high values of stresses calculated by other viscoplastic materials
106 446432 : if (equiv_stress > _maximum_equivalent_stress)
107 0 : mooseException("In ",
108 : _name,
109 : ": equivalent stress (",
110 : equiv_stress,
111 : ") is higher than maximum_equivalent_stress (",
112 : _maximum_equivalent_stress,
113 : ").\nCutting time step.");
114 :
115 : // If equivalent stress is present, calculate creep strain increment
116 446432 : if (equiv_stress > _minimum_equivalent_stress)
117 : {
118 : // Initalize stress potential
119 438592 : ADReal dpsi_dgauge(0);
120 :
121 438592 : computeInelasticStrainIncrement(_gauge_stress[_qp],
122 : dpsi_dgauge,
123 : inelastic_strain_increment,
124 : equiv_stress,
125 : dev_stress,
126 : stress);
127 :
128 : // Update elastic strain increment due to inelastic strain calculated here
129 438592 : elastic_strain_increment -= inelastic_strain_increment;
130 : // Update stress due to new strain
131 438592 : stress = elasticity_tensor * (elastic_strain_old + elastic_strain_increment);
132 :
133 : // Compute effective strain from the stress potential. Note that this is approximate and to be
134 : // used qualitatively
135 877184 : _effective_inelastic_strain[_qp] += dpsi_dgauge * _dt;
136 : // Update creep strain due to currently computed inelastic strain
137 438592 : _inelastic_strain[_qp] += inelastic_strain_increment;
138 : }
139 :
140 446432 : const ADRankTwoTensor new_dev_stress = stress.deviatoric();
141 446432 : const ADReal new_dev_stress_squared = new_dev_stress.doubleContraction(new_dev_stress);
142 : const ADReal new_equiv_stress =
143 1323616 : new_dev_stress_squared == 0.0 ? 0.0 : std::sqrt(1.5 * new_dev_stress_squared);
144 :
145 446432 : if (MooseUtils::relativeFuzzyGreaterThan(new_equiv_stress, equiv_stress))
146 0 : mooseException("In ",
147 : _name,
148 : ": updated equivalent stress (",
149 : MetaPhysicL::raw_value(new_equiv_stress),
150 : ") is greater than initial equivalent stress (",
151 : MetaPhysicL::raw_value(equiv_stress),
152 : "). Try decreasing max_inelastic_increment to avoid this exception.");
153 :
154 446432 : computeStressFinalize(inelastic_strain_increment);
155 446432 : }
156 :
157 : ADReal
158 433472 : ADViscoplasticityStressUpdate::initialGuess(const ADReal & effective_trial_stress)
159 : {
160 433472 : return effective_trial_stress;
161 : }
162 :
163 : ADReal
164 433472 : ADViscoplasticityStressUpdate::maximumPermissibleValue(const ADReal & effective_trial_stress) const
165 : {
166 433472 : return effective_trial_stress * _maximum_gauge_ratio;
167 : }
168 :
169 : ADReal
170 433472 : ADViscoplasticityStressUpdate::minimumPermissibleValue(const ADReal & effective_trial_stress) const
171 : {
172 433472 : return effective_trial_stress;
173 : }
174 :
175 : ADReal
176 3419804 : ADViscoplasticityStressUpdate::computeResidual(const ADReal & equiv_stress,
177 : const ADReal & trial_gauge)
178 : {
179 3419804 : const ADReal M = std::abs(_hydro_stress) / trial_gauge;
180 6839608 : const ADReal dM_dtrial_gauge = -M / trial_gauge;
181 :
182 3419804 : const ADReal residual_left = Utility::pow<2>(equiv_stress / trial_gauge);
183 3419804 : const ADReal dresidual_left_dtrial_gauge = -2.0 * residual_left / trial_gauge;
184 :
185 3419804 : ADReal residual = residual_left;
186 3419804 : _derivative = dresidual_left_dtrial_gauge;
187 :
188 3419804 : if (_pore_shape == PoreShapeModel::SPHERICAL)
189 : {
190 5359284 : residual *= 1.0 + _intermediate_porosity / 1.5;
191 5359284 : _derivative *= 1.0 + _intermediate_porosity / 1.5;
192 : }
193 :
194 3419804 : if (_model == ViscoplasticityModel::GTN)
195 : {
196 2127432 : residual += 2.0 * _intermediate_porosity * std::cosh(_pore_shape_factor * M) - 1.0 -
197 2127432 : Utility::pow<2>(_intermediate_porosity);
198 1063716 : _derivative += 2.0 * _intermediate_porosity * std::sinh(_pore_shape_factor * M) *
199 1063716 : _pore_shape_factor * dM_dtrial_gauge;
200 : }
201 : else
202 : {
203 2356088 : const ADReal h = computeH(_power, M);
204 2356088 : const ADReal dh_dM = computeH(_power, M, true);
205 :
206 2356088 : residual += _intermediate_porosity * (h + _power_factor / h) - 1.0 -
207 4712176 : _power_factor * Utility::pow<2>(_intermediate_porosity);
208 4712176 : const ADReal dresidual_dh = _intermediate_porosity * (1.0 - _power_factor / Utility::pow<2>(h));
209 2356088 : _derivative += dresidual_dh * dh_dM * dM_dtrial_gauge;
210 : }
211 :
212 3419804 : if (_verbose)
213 : {
214 : Moose::out << "in computeResidual:\n"
215 0 : << " position: " << _q_point[_qp] << " hydro_stress: " << _hydro_stress
216 : << " equiv_stress: " << equiv_stress << " trial_grage: " << trial_gauge
217 : << " M: " << M << std::endl;
218 : Moose::out << " residual: " << residual << " derivative: " << _derivative << std::endl;
219 : }
220 :
221 3419804 : return residual;
222 : }
223 :
224 : ADReal
225 5526576 : ADViscoplasticityStressUpdate::computeH(const Real n, const ADReal & M, const bool derivative)
226 : {
227 11053152 : const ADReal mod = std::pow(M * _pore_shape_factor, (n + 1.0) / n);
228 :
229 : // Calculate derivative with respect to M
230 5526576 : if (derivative)
231 : {
232 2731896 : const ADReal dmod_dM = (n + 1.0) / n * mod / M;
233 10927584 : return dmod_dM * std::pow(1.0 + mod / n, n - 1.0);
234 : }
235 :
236 5589360 : return std::pow(1.0 + mod / n, n);
237 : }
238 :
239 : ADRankTwoTensor
240 438592 : ADViscoplasticityStressUpdate::computeDGaugeDSigma(const ADReal & gauge_stress,
241 : const ADReal & equiv_stress,
242 : const ADRankTwoTensor & dev_stress,
243 : const ADRankTwoTensor & stress)
244 : {
245 : // Compute the derivative of the gauge stress with respect to the equivalent and hydrostatic
246 : // stress components
247 438592 : const ADReal M = std::abs(_hydro_stress) / gauge_stress;
248 438592 : const ADReal h = computeH(_power, M);
249 :
250 : // Compute the derviative of the residual with respect to the hydrostatic stress
251 438592 : ADReal dresidual_dhydro_stress = 0.0;
252 438592 : if (_hydro_stress != 0.0)
253 : {
254 : const ADReal dM_dhydro_stress = M / _hydro_stress;
255 435008 : if (_model == ViscoplasticityModel::GTN)
256 : {
257 59200 : dresidual_dhydro_stress = 2.0 * _intermediate_porosity * std::sinh(_pore_shape_factor * M) *
258 118400 : _pore_shape_factor * dM_dhydro_stress;
259 : }
260 : else
261 : {
262 : const ADReal dresidual_dh =
263 1127424 : _intermediate_porosity * (1.0 - _power_factor / Utility::pow<2>(h));
264 375808 : const ADReal dh_dM = computeH(_power, M, true);
265 375808 : dresidual_dhydro_stress = dresidual_dh * dh_dM * dM_dhydro_stress;
266 : }
267 : }
268 :
269 : // Compute the derivative of the residual with respect to the equivalent stress
270 438592 : ADReal dresidual_dequiv_stress = 2.0 * equiv_stress / Utility::pow<2>(gauge_stress);
271 438592 : if (_pore_shape == PoreShapeModel::SPHERICAL)
272 761376 : dresidual_dequiv_stress *= 1.0 + _intermediate_porosity / 1.5;
273 :
274 : // Compute the derivative of the equilvalent stress to the deviatoric stress
275 877184 : const ADRankTwoTensor dequiv_stress_dsigma = 1.5 * dev_stress / equiv_stress;
276 :
277 : // Compute the derivative of the residual with the stress
278 438592 : const ADRankTwoTensor dresidual_dsigma = dresidual_dhydro_stress * _dhydro_stress_dsigma +
279 438592 : dresidual_dequiv_stress * dequiv_stress_dsigma;
280 :
281 : // Compute the deritative of the gauge stress with respect to the stress
282 : const ADRankTwoTensor dgauge_dsigma =
283 877184 : dresidual_dsigma * (gauge_stress / dresidual_dsigma.doubleContraction(stress));
284 :
285 438592 : return dgauge_dsigma;
286 : }
287 :
288 : void
289 438592 : ADViscoplasticityStressUpdate::computeInelasticStrainIncrement(
290 : ADReal & gauge_stress,
291 : ADReal & dpsi_dgauge,
292 : ADRankTwoTensor & inelastic_strain_increment,
293 : const ADReal & equiv_stress,
294 : const ADRankTwoTensor & dev_stress,
295 : const ADRankTwoTensor & stress)
296 : {
297 : // If hydrostatic stress and porosity present, compute non-linear gauge stress
298 438592 : if (_intermediate_porosity == 0.0)
299 1536 : gauge_stress = equiv_stress;
300 437056 : else if (_hydro_stress == 0.0)
301 : gauge_stress =
302 10752 : equiv_stress / std::sqrt(1.0 - (1.0 + _power_factor) * _intermediate_porosity +
303 7168 : _power_factor * Utility::pow<2>(_intermediate_porosity));
304 : else
305 433472 : returnMappingSolve(equiv_stress, gauge_stress, _console);
306 :
307 : mooseAssert(gauge_stress >= equiv_stress,
308 : "Gauge stress calculated in inner Newton solve is less than the equivalent stress.");
309 :
310 : // Compute stress potential
311 877184 : dpsi_dgauge = _coefficient[_qp] * std::pow(gauge_stress, _power);
312 :
313 : // Compute strain increment from stress potential and the gauge stress derivative with respect
314 : // to the stress stress. The current form is explicit, and should eventually be changed
315 : inelastic_strain_increment =
316 438592 : _dt * dpsi_dgauge * computeDGaugeDSigma(gauge_stress, equiv_stress, dev_stress, stress);
317 438592 : }
318 :
319 : void
320 0 : ADViscoplasticityStressUpdate::outputIterationSummary(std::stringstream * iter_output,
321 : const unsigned int total_it)
322 : {
323 0 : if (iter_output)
324 : {
325 0 : *iter_output << "At element " << _current_elem->id() << " _qp=" << _qp << " Coordinates "
326 0 : << _q_point[_qp] << " block=" << _current_elem->subdomain_id() << '\n';
327 : }
328 0 : ADSingleVariableReturnMappingSolution::outputIterationSummary(iter_output, total_it);
329 0 : }
330 :
331 : Real
332 3419804 : ADViscoplasticityStressUpdate::computeReferenceResidual(const ADReal & /*effective_trial_stress*/,
333 : const ADReal & gauge_stress)
334 : {
335 : // Use gauge stress for relative tolerance criteria, defined as:
336 : // std::abs(residual / gauge_stress) <= _relative_tolerance
337 3419804 : return MetaPhysicL::raw_value(gauge_stress);
338 : }
|
