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