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 "ADComputeSmearedCrackingStress.h"
11 : #include "ElasticityTensorTools.h"
12 : #include "StressUpdateBase.h"
13 : #include "Conversion.h"
14 :
15 : registerMooseObject("SolidMechanicsApp", ADComputeSmearedCrackingStress);
16 :
17 : InputParameters
18 209 : ADComputeSmearedCrackingStress::validParams()
19 : {
20 209 : InputParameters params = ADComputeMultipleInelasticStress::validParams();
21 209 : params.addClassDescription(
22 : "Compute stress using a fixed smeared cracking model. Uses automatic differentiation");
23 418 : params.addRequiredParam<std::vector<MaterialName>>(
24 : "softening_models",
25 : "The material objects used to compute softening behavior for loading a crack."
26 : "Either 1 or 3 models must be specified. If a single model is specified, it is"
27 : "used for all directions. If 3 models are specified, they will be used for the"
28 : "3 crack directions in sequence");
29 418 : params.addRequiredCoupledVar(
30 : "cracking_stress",
31 : "The stress threshold beyond which cracking occurs. Negative values prevent cracking.");
32 209 : MultiMooseEnum direction("x y z");
33 418 : params.addParam<MultiMooseEnum>(
34 : "prescribed_crack_directions", direction, "Prescribed directions of first cracks");
35 418 : params.addParam<unsigned int>(
36 418 : "max_cracks", 3, "The maximum number of cracks allowed at a material point.");
37 418 : params.addRangeCheckedParam<Real>("cracking_neg_fraction",
38 : 0,
39 : "cracking_neg_fraction <= 1 & cracking_neg_fraction >= 0",
40 : "The fraction of the cracking strain at which "
41 : "a transition begins during decreasing "
42 : "strain to the original stiffness.");
43 418 : params.addParam<Real>(
44 : "max_stress_correction",
45 418 : 1.0,
46 : "Maximum permitted correction to the predicted stress as a ratio of the "
47 : "stress change to the predicted stress from the previous step's damage level. "
48 : "Values less than 1 will improve robustness, but not be as accurate.");
49 :
50 627 : params.addRangeCheckedParam<Real>(
51 : "shear_retention_factor",
52 418 : 0.0,
53 : "shear_retention_factor>=0 & shear_retention_factor<=1.0",
54 : "Fraction of original shear stiffness to be retained after cracking");
55 209 : params.set<std::vector<MaterialName>>("inelastic_models") = {};
56 :
57 418 : MooseEnum crackedElasticityType("DIAGONAL FULL", "DIAGONAL");
58 418 : crackedElasticityType.addDocumentation(
59 : "DIAGONAL", "Zero out terms coupling with directions orthogonal to a crack (legacy)");
60 418 : crackedElasticityType.addDocumentation(
61 : "FULL", "Consistently scale all entries based on damage (recommended)");
62 418 : params.addParam<MooseEnum>(
63 : "cracked_elasticity_type",
64 : crackedElasticityType,
65 : "Method to modify the local elasticity tensor to account for cracking");
66 :
67 209 : return params;
68 209 : }
69 :
70 157 : ADComputeSmearedCrackingStress::ADComputeSmearedCrackingStress(const InputParameters & parameters)
71 : : ADComputeMultipleInelasticStress(parameters),
72 157 : _cracking_stress(adCoupledValue("cracking_stress")),
73 314 : _max_cracks(getParam<unsigned int>("max_cracks")),
74 314 : _cracking_neg_fraction(getParam<Real>("cracking_neg_fraction")),
75 314 : _shear_retention_factor(getParam<Real>("shear_retention_factor")),
76 314 : _max_stress_correction(getParam<Real>("max_stress_correction")),
77 157 : _cracked_elasticity_type(
78 157 : getParam<MooseEnum>("cracked_elasticity_type").getEnum<CrackedElasticityType>()),
79 157 : _crack_damage(declareADProperty<RealVectorValue>(_base_name + "crack_damage")),
80 314 : _crack_damage_old(getMaterialPropertyOld<RealVectorValue>(_base_name + "crack_damage")),
81 157 : _crack_flags(declareADProperty<RealVectorValue>(_base_name + "crack_flags")),
82 157 : _crack_rotation(declareADProperty<RankTwoTensor>(_base_name + "crack_rotation")),
83 314 : _crack_rotation_old(getMaterialPropertyOld<RankTwoTensor>(_base_name + "crack_rotation")),
84 157 : _crack_initiation_strain(
85 157 : declareADProperty<RealVectorValue>(_base_name + "crack_initiation_strain")),
86 157 : _crack_initiation_strain_old(
87 157 : getMaterialPropertyOld<RealVectorValue>(_base_name + "crack_initiation_strain")),
88 157 : _crack_max_strain(declareADProperty<RealVectorValue>(_base_name + "crack_max_strain")),
89 471 : _crack_max_strain_old(getMaterialPropertyOld<RealVectorValue>(_base_name + "crack_max_strain"))
90 : {
91 : MultiMooseEnum prescribed_crack_directions =
92 314 : getParam<MultiMooseEnum>("prescribed_crack_directions");
93 157 : if (prescribed_crack_directions.size() > 0)
94 : {
95 39 : if (prescribed_crack_directions.size() > 3)
96 0 : mooseError("A maximum of three crack directions may be specified");
97 126 : for (unsigned int i = 0; i < prescribed_crack_directions.size(); ++i)
98 : {
99 156 : for (unsigned int j = 0; j < i; ++j)
100 69 : if (prescribed_crack_directions[i] == prescribed_crack_directions[j])
101 0 : mooseError("Entries in 'prescribed_crack_directions' cannot be repeated");
102 87 : _prescribed_crack_directions.push_back(
103 87 : static_cast<unsigned int>(prescribed_crack_directions.get(i)));
104 : }
105 :
106 : // Fill in the last remaining direction if 2 are specified
107 39 : if (_prescribed_crack_directions.size() == 2)
108 : {
109 6 : std::set<unsigned int> available_dirs = {0, 1, 2};
110 18 : for (auto dir : _prescribed_crack_directions)
111 12 : if (available_dirs.erase(dir) != 1)
112 0 : mooseError("Invalid prescribed crack direction:" + Moose::stringify(dir));
113 6 : if (available_dirs.size() != 1)
114 0 : mooseError("Error in finding remaining available crack direction");
115 6 : _prescribed_crack_directions.push_back(*available_dirs.begin());
116 : }
117 : }
118 314 : if (!isParamSetByUser("cracked_elasticity_type"))
119 31 : paramWarning(
120 : "cracked_elasticity_type",
121 : "Defaulting to the legacy option of 'DIAGONAL', but the 'FULL' option is preferred");
122 :
123 156 : _local_elastic_vector.resize(9);
124 156 : }
125 :
126 : void
127 616 : ADComputeSmearedCrackingStress::initQpStatefulProperties()
128 : {
129 616 : ADComputeMultipleInelasticStress::initQpStatefulProperties();
130 :
131 616 : _crack_damage[_qp] = 0.0;
132 :
133 616 : _crack_initiation_strain[_qp] = 0.0;
134 616 : _crack_max_strain[_qp](0) = 0.0;
135 :
136 616 : switch (_prescribed_crack_directions.size())
137 : {
138 440 : case 0:
139 : {
140 440 : _crack_rotation[_qp] = ADRankTwoTensor::Identity();
141 440 : break;
142 : }
143 64 : case 1:
144 : {
145 64 : _crack_rotation[_qp].zero();
146 64 : switch (_prescribed_crack_directions[0])
147 : {
148 32 : case 0:
149 : {
150 32 : _crack_rotation[_qp](0, 0) = 1.0;
151 32 : _crack_rotation[_qp](1, 1) = 1.0;
152 32 : _crack_rotation[_qp](2, 2) = 1.0;
153 32 : break;
154 : }
155 0 : case 1:
156 : {
157 0 : _crack_rotation[_qp](1, 0) = 1.0;
158 0 : _crack_rotation[_qp](0, 1) = 1.0;
159 0 : _crack_rotation[_qp](2, 2) = 1.0;
160 0 : break;
161 : }
162 32 : case 2:
163 : {
164 32 : _crack_rotation[_qp](2, 0) = 1.0;
165 32 : _crack_rotation[_qp](0, 1) = 1.0;
166 32 : _crack_rotation[_qp](1, 2) = 1.0;
167 32 : break;
168 : }
169 : }
170 : break;
171 : }
172 0 : case 2:
173 : {
174 0 : mooseError("Number of prescribed crack directions cannot be 2");
175 : break;
176 : }
177 : case 3:
178 : {
179 448 : for (unsigned int i = 0; i < _prescribed_crack_directions.size(); ++i)
180 : {
181 : ADRealVectorValue crack_dir_vec;
182 336 : crack_dir_vec(_prescribed_crack_directions[i]) = 1.0;
183 336 : _crack_rotation[_qp].fillColumn(i, crack_dir_vec);
184 : }
185 : }
186 : }
187 616 : }
188 :
189 : void
190 154 : ADComputeSmearedCrackingStress::initialSetup()
191 : {
192 154 : ADComputeMultipleInelasticStress::initialSetup();
193 :
194 308 : if (!hasGuaranteedMaterialProperty(_elasticity_tensor_name, Guarantee::ISOTROPIC))
195 0 : mooseError("ADComputeSmearedCrackingStress requires that the elasticity tensor be "
196 : "guaranteed isotropic");
197 :
198 462 : std::vector<MaterialName> soft_matls = getParam<std::vector<MaterialName>>("softening_models");
199 351 : for (auto soft_matl : soft_matls)
200 : {
201 : ADSmearedCrackSofteningBase * scsb =
202 197 : dynamic_cast<ADSmearedCrackSofteningBase *>(&getMaterialByName(soft_matl));
203 197 : if (scsb)
204 197 : _softening_models.push_back(scsb);
205 : else
206 0 : paramError("softening_models", "Model " + soft_matl + " is not a softening model");
207 : }
208 154 : if (_softening_models.size() == 1)
209 : {
210 : // Reuse the same model in all 3 directions
211 132 : _softening_models.push_back(_softening_models[0]);
212 132 : _softening_models.push_back(_softening_models[0]);
213 : }
214 22 : else if (_softening_models.size() != 3)
215 1 : paramError("softening_models", "Either 1 or 3 softening models must be specified");
216 153 : }
217 :
218 : void
219 463820 : ADComputeSmearedCrackingStress::computeQpStress()
220 : {
221 :
222 463820 : if (!previouslyCracked())
223 32808 : computeQpStressIntermediateConfiguration();
224 : else
225 : {
226 431012 : _elastic_strain[_qp] = _elastic_strain_old[_qp] + _strain_increment[_qp];
227 :
228 : // Propagate behavior from the (now inactive) inelastic models
229 431012 : _inelastic_strain[_qp] = _inelastic_strain_old[_qp];
230 431012 : for (auto model : _models)
231 : {
232 0 : model->setQp(_qp);
233 0 : model->propagateQpStatefulProperties();
234 : }
235 :
236 : // Since the other inelastic models are inactive, they will not limit the time step
237 431012 : _material_timestep_limit[_qp] = std::numeric_limits<Real>::max();
238 :
239 : // update _local_elasticity_tensor based on cracking state in previous time step
240 431012 : updateLocalElasticityTensor();
241 :
242 : // Calculate stress in intermediate configuration
243 431012 : _stress[_qp] = _local_elasticity_tensor * _elastic_strain[_qp];
244 : }
245 :
246 : // compute crack status and adjust stress
247 463820 : updateCrackingStateAndStress();
248 :
249 463820 : if (_perform_finite_strain_rotations)
250 : {
251 463820 : finiteStrainRotation();
252 463820 : _crack_rotation[_qp] = _rotation_increment[_qp] * _crack_rotation[_qp];
253 : }
254 463820 : }
255 :
256 : void
257 431012 : ADComputeSmearedCrackingStress::updateLocalElasticityTensor()
258 : {
259 : const ADReal youngs_modulus =
260 431012 : ElasticityTensorTools::getIsotropicYoungsModulus(_elasticity_tensor[_qp]);
261 :
262 : bool cracking_locally_active = false;
263 :
264 431012 : const ADReal cracking_stress = _cracking_stress[_qp];
265 :
266 431012 : if (cracking_stress > 0)
267 : {
268 431012 : ADRealVectorValue stiffness_ratio_local(1.0, 1.0, 1.0);
269 431012 : const ADRankTwoTensor & R = _crack_rotation_old[_qp];
270 431012 : ADRankTwoTensor ePrime(_elastic_strain_old[_qp]);
271 431012 : ePrime.rotate(R.transpose());
272 :
273 1724048 : for (unsigned int i = 0; i < 3; ++i)
274 : {
275 : // Update elasticity tensor based on crack status of the end of last time step
276 1293036 : if (_crack_damage_old[_qp](i) > 0.0)
277 : {
278 1197960 : if (_cracking_neg_fraction == 0.0 && MooseUtils::absoluteFuzzyLessThan(ePrime(i, i), 0.0))
279 57936 : stiffness_ratio_local(i) = 1.0;
280 0 : else if (_cracking_neg_fraction > 0.0 &&
281 1082088 : ePrime(i, i) < _crack_initiation_strain_old[_qp](i) * _cracking_neg_fraction &&
282 541044 : ePrime(i, i) > -_crack_initiation_strain_old[_qp](i) * _cracking_neg_fraction)
283 : {
284 0 : const ADReal etr = _cracking_neg_fraction * _crack_initiation_strain_old[_qp](i);
285 : const ADReal Eo = cracking_stress / _crack_initiation_strain_old[_qp](i);
286 0 : const ADReal Ec = Eo * (1.0 - _crack_damage_old[_qp](i));
287 0 : const ADReal a = (Ec - Eo) / (4 * etr);
288 0 : const ADReal b = (Ec + Eo) / 2;
289 : // Compute the ratio of the current transition stiffness to the original stiffness
290 0 : stiffness_ratio_local(i) = (2.0 * a * etr + b) / Eo;
291 : cracking_locally_active = true;
292 : }
293 : else
294 : {
295 541044 : stiffness_ratio_local(i) = (1.0 - _crack_damage_old[_qp](i));
296 : cracking_locally_active = true;
297 : }
298 : }
299 : }
300 :
301 431012 : if (cracking_locally_active)
302 : {
303 : // Update the elasticity tensor in the crack coordinate system
304 373508 : if (_cracked_elasticity_type == CrackedElasticityType::DIAGONAL)
305 : {
306 1104 : const bool c0_coupled = MooseUtils::absoluteFuzzyEqual(stiffness_ratio_local(0), 1.0);
307 1104 : const bool c1_coupled = MooseUtils::absoluteFuzzyEqual(stiffness_ratio_local(1), 1.0);
308 1104 : const bool c2_coupled = MooseUtils::absoluteFuzzyEqual(stiffness_ratio_local(2), 1.0);
309 :
310 2208 : const ADReal c01 = (c0_coupled && c1_coupled ? 1.0 : 0.0);
311 1776 : const ADReal c02 = (c0_coupled && c2_coupled ? 1.0 : 0.0);
312 1536 : const ADReal c12 = (c1_coupled && c2_coupled ? 1.0 : 0.0);
313 :
314 : const ADReal c01_shear_retention =
315 1104 : (c0_coupled && c1_coupled ? 1.0 : _shear_retention_factor);
316 : const ADReal c02_shear_retention =
317 1104 : (c0_coupled && c2_coupled ? 1.0 : _shear_retention_factor);
318 : const ADReal c12_shear_retention =
319 1104 : (c1_coupled && c2_coupled ? 1.0 : _shear_retention_factor);
320 :
321 1104 : _local_elastic_vector[0] = (c0_coupled ? _elasticity_tensor[_qp](0, 0, 0, 0)
322 1104 : : stiffness_ratio_local(0) * youngs_modulus);
323 2208 : _local_elastic_vector[1] = _elasticity_tensor[_qp](0, 0, 1, 1) * c01;
324 2208 : _local_elastic_vector[2] = _elasticity_tensor[_qp](0, 0, 2, 2) * c02;
325 1104 : _local_elastic_vector[3] = (c1_coupled ? _elasticity_tensor[_qp](1, 1, 1, 1)
326 1104 : : stiffness_ratio_local(1) * youngs_modulus);
327 2208 : _local_elastic_vector[4] = _elasticity_tensor[_qp](1, 1, 2, 2) * c12;
328 1104 : _local_elastic_vector[5] = (c2_coupled ? _elasticity_tensor[_qp](2, 2, 2, 2)
329 1104 : : stiffness_ratio_local(2) * youngs_modulus);
330 2208 : _local_elastic_vector[6] = _elasticity_tensor[_qp](1, 2, 1, 2) * c12_shear_retention;
331 2208 : _local_elastic_vector[7] = _elasticity_tensor[_qp](0, 2, 0, 2) * c02_shear_retention;
332 2208 : _local_elastic_vector[8] = _elasticity_tensor[_qp](0, 1, 0, 1) * c01_shear_retention;
333 : }
334 : else // _cracked_elasticity_type == CrackedElasticityType::FULL
335 : {
336 : const ADReal & c0 = stiffness_ratio_local(0);
337 : const ADReal & c1 = stiffness_ratio_local(1);
338 : const ADReal & c2 = stiffness_ratio_local(2);
339 :
340 : const ADReal c01 = c0 * c1;
341 : const ADReal c02 = c0 * c2;
342 : const ADReal c12 = c1 * c2;
343 :
344 372692 : const ADReal c01_shear_retention = std::max(c01, _shear_retention_factor);
345 372404 : const ADReal c02_shear_retention = std::max(c02, _shear_retention_factor);
346 372404 : const ADReal c12_shear_retention = std::max(c12, _shear_retention_factor);
347 :
348 744808 : _local_elastic_vector[0] = _elasticity_tensor[_qp](0, 0, 0, 0) * c0;
349 744808 : _local_elastic_vector[1] = _elasticity_tensor[_qp](0, 0, 1, 1) * c01;
350 744808 : _local_elastic_vector[2] = _elasticity_tensor[_qp](0, 0, 2, 2) * c02;
351 744808 : _local_elastic_vector[3] = _elasticity_tensor[_qp](1, 1, 1, 1) * c1;
352 744808 : _local_elastic_vector[4] = _elasticity_tensor[_qp](1, 1, 2, 2) * c12;
353 744808 : _local_elastic_vector[5] = _elasticity_tensor[_qp](2, 2, 2, 2) * c2;
354 744808 : _local_elastic_vector[6] = _elasticity_tensor[_qp](1, 2, 1, 2) * c12_shear_retention;
355 744808 : _local_elastic_vector[7] = _elasticity_tensor[_qp](0, 2, 0, 2) * c02_shear_retention;
356 744808 : _local_elastic_vector[8] = _elasticity_tensor[_qp](0, 1, 0, 1) * c01_shear_retention;
357 : }
358 :
359 : // Filling with 9 components is sufficient because these are the only nonzero entries
360 : // for isotropic or orthotropic materials.
361 373508 : _local_elasticity_tensor.fillFromInputVector(_local_elastic_vector,
362 : ADRankFourTensor::symmetric9);
363 :
364 : // Rotate the modified elasticity tensor back into global coordinates
365 373508 : _local_elasticity_tensor.rotate(R);
366 : }
367 : }
368 : if (!cracking_locally_active)
369 57504 : _local_elasticity_tensor = _elasticity_tensor[_qp];
370 431012 : }
371 :
372 : void
373 463820 : ADComputeSmearedCrackingStress::updateCrackingStateAndStress()
374 : {
375 : const ADReal youngs_modulus =
376 463820 : ElasticityTensorTools::getIsotropicYoungsModulus(_elasticity_tensor[_qp]);
377 :
378 : const ADReal poissons_ratio =
379 463820 : ElasticityTensorTools::getIsotropicPoissonsRatio(_elasticity_tensor[_qp]);
380 :
381 463820 : ADReal cracking_stress = _cracking_stress[_qp];
382 :
383 463820 : if (cracking_stress > 0)
384 : {
385 : // Initializing crack states
386 463784 : _crack_rotation[_qp] = _crack_rotation_old[_qp];
387 :
388 1855136 : for (unsigned i = 0; i < 3; ++i)
389 : {
390 1391352 : _crack_max_strain[_qp](i) = _crack_max_strain_old[_qp](i);
391 1391352 : _crack_initiation_strain[_qp](i) = _crack_initiation_strain_old[_qp](i);
392 1391352 : _crack_damage[_qp](i) = _crack_damage_old[_qp](i);
393 : }
394 :
395 : // Compute crack orientations: updated _crack_rotation[_qp] based on current strain
396 : ADRealVectorValue strain_in_crack_dir;
397 463784 : computeCrackStrainAndOrientation(strain_in_crack_dir);
398 :
399 1855136 : for (unsigned i = 0; i < 3; ++i)
400 : {
401 1391352 : if (strain_in_crack_dir(i) > _crack_max_strain[_qp](i))
402 487657 : _crack_max_strain[_qp](i) = strain_in_crack_dir(i);
403 : }
404 :
405 : // Check for new cracks.
406 : // Rotate stress to cracked orientation.
407 463784 : const ADRankTwoTensor & R = _crack_rotation[_qp];
408 463784 : ADRankTwoTensor sigmaPrime(_stress[_qp]);
409 463784 : sigmaPrime.rotate(R.transpose()); // stress in crack coordinates
410 :
411 : unsigned int num_cracks = 0;
412 1855136 : for (unsigned int i = 0; i < 3; ++i)
413 : {
414 1391352 : if (_crack_damage_old[_qp](i) > 0.0)
415 598980 : ++num_cracks;
416 : }
417 :
418 : bool cracked(false);
419 : ADRealVectorValue sigma;
420 : mooseAssert(_softening_models.size() == 3, "Must have 3 softening models");
421 1855136 : for (unsigned int i = 0; i < 3; ++i)
422 : {
423 1391352 : sigma(i) = sigmaPrime(i, i);
424 :
425 1391352 : ADReal stiffness_ratio = 1.0 - _crack_damage[_qp](i);
426 :
427 1391352 : const bool pre_existing_crack = (_crack_damage_old[_qp](i) > 0.0);
428 : const bool met_stress_criterion = (sigma(i) > cracking_stress);
429 1391352 : const bool loading_existing_crack = (strain_in_crack_dir(i) >= _crack_max_strain[_qp](i));
430 1391352 : const bool allowed_to_crack = (pre_existing_crack || num_cracks < _max_cracks);
431 : bool new_crack = false;
432 :
433 1391352 : cracked |= pre_existing_crack;
434 :
435 : // Adjustments for newly created cracks
436 1391352 : if (met_stress_criterion && !pre_existing_crack && allowed_to_crack)
437 : {
438 : new_crack = true;
439 4684 : ++num_cracks;
440 :
441 : // Assume Poisson's ratio drops to zero for this direction. Stiffness is then Young's
442 : // modulus.
443 4684 : _crack_initiation_strain[_qp](i) = cracking_stress / youngs_modulus;
444 :
445 4684 : if (_crack_max_strain[_qp](i) < _crack_initiation_strain[_qp](i))
446 408 : _crack_max_strain[_qp](i) = _crack_initiation_strain[_qp](i);
447 : }
448 :
449 : // Update stress and stiffness ratio according to specified crack release model
450 1386668 : if (new_crack || (pre_existing_crack && loading_existing_crack))
451 : {
452 : cracked = true;
453 386968 : _softening_models[i]->computeCrackingRelease(sigma(i),
454 : stiffness_ratio,
455 : strain_in_crack_dir(i),
456 386968 : _crack_initiation_strain[_qp](i),
457 386968 : _crack_max_strain[_qp](i),
458 : cracking_stress,
459 : youngs_modulus,
460 : poissons_ratio);
461 773936 : _crack_damage[_qp](i) = 1.0 - stiffness_ratio;
462 : }
463 :
464 820376 : else if (cracked && _cracking_neg_fraction > 0 &&
465 2008768 : _crack_initiation_strain[_qp](i) * _cracking_neg_fraction > strain_in_crack_dir(i) &&
466 1004384 : -_crack_initiation_strain[_qp](i) * _cracking_neg_fraction < strain_in_crack_dir(i))
467 : {
468 0 : const ADReal etr = _cracking_neg_fraction * _crack_initiation_strain[_qp](i);
469 0 : const ADReal Eo = cracking_stress / _crack_initiation_strain[_qp](i);
470 0 : const ADReal Ec = Eo * (1.0 - _crack_damage_old[_qp](i));
471 0 : const ADReal a = (Ec - Eo) / (4.0 * etr);
472 0 : const ADReal b = 0.5 * (Ec + Eo);
473 0 : const ADReal c = 0.25 * (Ec - Eo) * etr;
474 0 : sigma(i) = (a * strain_in_crack_dir(i) + b) * strain_in_crack_dir(i) + c;
475 : }
476 : }
477 :
478 463784 : if (cracked)
479 432688 : updateStressTensorForCracking(_stress[_qp], sigma);
480 : }
481 :
482 927640 : _crack_flags[_qp](0) = 1.0 - _crack_damage[_qp](2);
483 927640 : _crack_flags[_qp](1) = 1.0 - _crack_damage[_qp](1);
484 927640 : _crack_flags[_qp](2) = 1.0 - _crack_damage[_qp](0);
485 463820 : }
486 :
487 : void
488 463784 : ADComputeSmearedCrackingStress::computeCrackStrainAndOrientation(
489 : ADRealVectorValue & strain_in_crack_dir)
490 : {
491 : // The rotation tensor is ordered such that directions for pre-existing cracks appear first
492 : // in the list of columns. For example, if there is one existing crack, its direction is in the
493 : // first column in the rotation tensor.
494 463784 : const unsigned int num_known_dirs = getNumKnownCrackDirs();
495 :
496 463784 : if (num_known_dirs == 0)
497 : {
498 28500 : std::vector<ADReal> eigval(3, 0.0);
499 28500 : ADRankTwoTensor eigvec;
500 :
501 28500 : _elastic_strain[_qp].symmetricEigenvaluesEigenvectors(eigval, eigvec);
502 :
503 : // If the elastic strain is beyond the cracking strain, save the eigen vectors as
504 : // the rotation tensor. Reverse their order so that the third principal strain
505 : // (most tensile) will correspond to the first crack.
506 28500 : _crack_rotation[_qp].fillColumn(0, eigvec.column(2));
507 28500 : _crack_rotation[_qp].fillColumn(1, eigvec.column(1));
508 28500 : _crack_rotation[_qp].fillColumn(2, eigvec.column(0));
509 :
510 28500 : strain_in_crack_dir(0) = eigval[2];
511 28500 : strain_in_crack_dir(1) = eigval[1];
512 28500 : strain_in_crack_dir(2) = eigval[0];
513 28500 : }
514 435284 : else if (num_known_dirs == 1)
515 : {
516 : // This is easily the most complicated case.
517 : // 1. Rotate the elastic strain to the orientation associated with the known
518 : // crack.
519 : // 2. Extract the lower 2x2 block into a separate tensor.
520 : // 3. Run the eigen solver on the result.
521 : // 4. Update the rotation tensor to reflect the effect of the 2 eigenvectors.
522 :
523 : // 1.
524 194620 : const ADRankTwoTensor & R = _crack_rotation[_qp];
525 194620 : RankTwoTensor ePrime(raw_value(_elastic_strain[_qp]));
526 389240 : ePrime.rotate(raw_value(R.transpose())); // elastic strain in crack coordinates
527 :
528 : // 2.
529 194620 : ColumnMajorMatrix e2x2(2, 2);
530 194620 : e2x2(0, 0) = ePrime(1, 1);
531 194620 : e2x2(1, 0) = ePrime(2, 1);
532 194620 : e2x2(0, 1) = ePrime(1, 2);
533 194620 : e2x2(1, 1) = ePrime(2, 2);
534 :
535 : // 3.
536 194620 : ColumnMajorMatrix e_val2x1(2, 1);
537 194620 : ColumnMajorMatrix e_vec2x2(2, 2);
538 194620 : e2x2.eigen(e_val2x1, e_vec2x2);
539 :
540 : // 4.
541 : ADRankTwoTensor eigvec(
542 194620 : 1.0, 0.0, 0.0, 0.0, e_vec2x2(0, 1), e_vec2x2(1, 1), 0.0, e_vec2x2(0, 0), e_vec2x2(1, 0));
543 :
544 194620 : _crack_rotation[_qp] = _crack_rotation_old[_qp] * eigvec; // Roe implementation
545 :
546 194620 : strain_in_crack_dir(0) = ePrime(0, 0);
547 194620 : strain_in_crack_dir(1) = e_val2x1(1, 0);
548 194620 : strain_in_crack_dir(2) = e_val2x1(0, 0);
549 : }
550 240664 : else if (num_known_dirs == 2 || num_known_dirs == 3)
551 : {
552 : // Rotate to cracked orientation and pick off the strains in the rotated
553 : // coordinate directions.
554 240664 : const ADRankTwoTensor & R = _crack_rotation[_qp];
555 240664 : ADRankTwoTensor ePrime(_elastic_strain[_qp]);
556 240664 : ePrime.rotate(R.transpose()); // elastic strain in crack coordinates
557 :
558 240664 : strain_in_crack_dir(0) = ePrime(0, 0);
559 : strain_in_crack_dir(1) = ePrime(1, 1);
560 240664 : strain_in_crack_dir(2) = ePrime(2, 2);
561 : }
562 : else
563 0 : mooseError("Invalid number of known crack directions");
564 463784 : }
565 :
566 : unsigned int
567 463784 : ADComputeSmearedCrackingStress::getNumKnownCrackDirs() const
568 : {
569 : unsigned int num_known_dirs = 0;
570 1855136 : for (unsigned int i = 0; i < 3; ++i)
571 : {
572 1391352 : if (_crack_damage_old[_qp](i) > 0.0 || _prescribed_crack_directions.size() >= i + 1)
573 859972 : ++num_known_dirs;
574 : }
575 463784 : return num_known_dirs;
576 : }
577 :
578 : void
579 432688 : ADComputeSmearedCrackingStress::updateStressTensorForCracking(ADRankTwoTensor & tensor,
580 : const ADRealVectorValue & sigma)
581 : {
582 : // Get transformation matrix
583 432688 : const ADRankTwoTensor & R = _crack_rotation[_qp];
584 : // Rotate to crack frame
585 432688 : tensor.rotate(R.transpose());
586 :
587 : // Reset stress if cracked
588 1730752 : for (unsigned int i = 0; i < 3; ++i)
589 1298064 : if (_crack_damage[_qp](i) > 0.0)
590 : {
591 601960 : const ADReal stress_correction_ratio = (tensor(i, i) - sigma(i)) / tensor(i, i);
592 601960 : if (stress_correction_ratio > _max_stress_correction)
593 0 : tensor(i, i) *= (1.0 - _max_stress_correction);
594 601960 : else if (stress_correction_ratio < -_max_stress_correction)
595 0 : tensor(i, i) *= (1.0 + _max_stress_correction);
596 : else
597 601960 : tensor(i, i) = sigma(i);
598 : }
599 :
600 : // Rotate back to global frame
601 432688 : tensor.rotate(R);
602 432688 : }
603 :
604 : bool
605 463820 : ADComputeSmearedCrackingStress::previouslyCracked()
606 : {
607 652772 : for (unsigned int i = 0; i < 3; ++i)
608 619964 : if (_crack_damage_old[_qp](i) > 0.0)
609 : return true;
610 : return false;
611 : }
|
