Line data Source code
1 : /*************************************************/
2 : /* DO NOT MODIFY THIS HEADER */
3 : /* */
4 : /* MASTODON */
5 : /* */
6 : /* (c) 2015 Battelle Energy Alliance, LLC */
7 : /* ALL RIGHTS RESERVED */
8 : /* */
9 : /* Prepared by Battelle Energy Alliance, LLC */
10 : /* With the U. S. Department of Energy */
11 : /* */
12 : /* See COPYRIGHT for full restrictions */
13 : /*************************************************/
14 :
15 : // This code was implemented in colloboration with Ozgun Numanoglu
16 : // ([email protected]) and Omar Baltaji ([email protected]) from UIUC.
17 :
18 : #include "ADComputeISoilStress.h"
19 :
20 : #include "MooseMesh.h"
21 : #include "Function.h"
22 : #include "MastodonUtils.h"
23 : #include "ISoilUtils.h"
24 : #include "Conversion.h"
25 : #include "FEProblem.h"
26 : #include "MooseUtils.h"
27 :
28 : #include "metaphysicl/numberarray.h"
29 : #include "metaphysicl/dualnumber.h"
30 :
31 : registerMooseObject("MastodonApp", ADComputeISoilStress);
32 :
33 : InputParameters
34 168 : ADComputeISoilStress::validParams()
35 : {
36 168 : InputParameters params = ADComputeFiniteStrainElasticStress::validParams();
37 168 : params.addClassDescription(
38 : "Compute total stress for the nonlinear material "
39 : "model I-Soil using a backbone curve and automatically computes the Jacobian.");
40 336 : params.addRequiredCoupledVar("layer_variable",
41 : "The auxvariable providing the soil layer identification.");
42 336 : params.addRequiredParam<std::vector<unsigned int>>(
43 : "layer_ids",
44 : "Vector of layer ids that map one-to-one to the rest of the "
45 : "soil layer parameters provided as input.");
46 336 : params.addRequiredParam<std::vector<Real>>("poissons_ratio",
47 : "Poissons's ratio for the soil layers. The "
48 : "size of the vector should be same as the size of "
49 : "layer_ids.");
50 336 : params.addParam<Real>("b_exp",
51 336 : 0.0,
52 : "The exponential factors for pressure "
53 : "dependent stiffness for all the soil "
54 : "layers.");
55 336 : params.addParam<std::vector<Real>>("p_ref",
56 : {},
57 : "The reference pressure at which "
58 : "the parameters are defined for "
59 : "each soil layer. If 'soil_type = "
60 : "darendeli', then the reference "
61 : "pressure must be input in kilopascals.");
62 336 : params.addParam<Real>("a0",
63 336 : 1.0,
64 : "The first coefficient for pressure dependent yield strength "
65 : "calculation for all the soil layers. If a0 = 1, a1 = 0 and "
66 : "a2=0 for one soil layer, then the yield strength of that "
67 : "layer is independent of pressure.");
68 336 : params.addParam<Real>("a1",
69 336 : 0.0,
70 : "The second coefficient for pressure dependent yield "
71 : "strength calculation for all the soil layers. If a0 = "
72 : "1, a1 = 0, a2 = 0 for one soil layer, then the yield "
73 : "strength of that layer is independent of pressure.");
74 336 : params.addParam<Real>("a2",
75 336 : 0.0,
76 : "The third coefficient for pressure dependent yield "
77 : "strength calculation for all the soil layers. If a0 = "
78 : "1, a1=0 and a2=0 for one soil layer, then the yield "
79 : "strength of that layer is independent of pressure.");
80 336 : params.addParam<Real>("tension_pressure_cut_off",
81 336 : -1.0,
82 : "The tension cut-off for all the soil layers. If the "
83 : "pressure becomes lower than this value, then the "
84 : "stiffness of the soil reduces to zero. A negative "
85 : "pressure indicates tension. The default "
86 : "value is -1.0 for all the soil layers.");
87 336 : params.addParam<bool>("pressure_dependency",
88 336 : false,
89 : "Set to true to turn on pressure dependent stiffness "
90 : "and yield strength calculation.");
91 336 : params.addParam<bool>(
92 336 : "wave_speed_calculation", true, "Set to false to turn off P and S wave speed calculation.");
93 336 : params.addParam<std::vector<FunctionName>>(
94 : "initial_soil_stress",
95 : {},
96 : "A list of functions describing the initial stress. There "
97 : "must be 9 functions, corresponding to the xx, yx, zx, xy, yy, zy, xz, yz, "
98 : "zz components respectively. If not provided, all components of the "
99 : "initial stress will be zero.");
100 : // params for specific backbone types
101 336 : MooseEnum soil_type("user_defined darendeli gqh thin_layer");
102 336 : params.addRequiredParam<MooseEnum>(
103 : "soil_type",
104 : soil_type,
105 : "This parameter determines the type of backbone curve used. Use 'user_defined' "
106 : "for a user defined backbone curve provided in a data file, 'darendeli' "
107 : "if the backbone curve is to be determined using Darandeli equations, 'gqh' "
108 : "if the backbone curve is determined using the GQ/H approach and 'thin_layer' if the soil "
109 : "is "
110 : "being used to simulate a thin-layer friction interface.");
111 : // params required for user_defined backbone curve: soil_type = 'user_defined'
112 336 : params.addParam<std::vector<FileName>>(
113 : "backbone_curve_files",
114 : "The vector of file names of the files containing "
115 : "stress-strain backbone curves for the different soil layers. The "
116 : "size of the vector should be same as the size of layer_ids. All files "
117 : "should contain the same number of stress-strain points. Headers are not "
118 : "expected and it is assumed that the first column corresponds to strain values "
119 : "and the second column corresponds to the stress values. Additionally, two "
120 : "segments of a backbone curve cannot have the same slope.");
121 : // params required for soil_type = 'darendeli', 'GQ/H' and 'thin_layer'
122 336 : params.addParam<std::vector<Real>>(
123 : "initial_shear_modulus",
124 : "The initial shear modulus of the soil layers. "
125 : "This is required if Darandeli or GQ/H type backbone curves are used.");
126 : // params required for soil_type = 'darendeli' and 'GQ/H'
127 336 : params.addParam<unsigned int>("number_of_points",
128 : "The total number of data points in which the "
129 : "backbone curve needs to be split for all soil "
130 : "layers (required for Darandeli or GQ/H type backbone curves).");
131 : // params required for Darandeli backbone curve: soil_type = 'darendeli'
132 336 : params.addParam<std::vector<Real>>("over_consolidation_ratio",
133 : "The over consolidation ratio of the soil "
134 : "layers. Required for Darandeli backbone curve.");
135 336 : params.addParam<std::vector<Real>>(
136 : "plasticity_index",
137 : "The plasticity index of the soil layers. Required for Darandeli backbone curve.");
138 : // params required for GQ/H backbone curve: soil_type = "gqh"
139 336 : params.addParam<std::vector<Real>>("theta_1",
140 : "The curve fit coefficient for "
141 : "GQ/H model"
142 : "for each soil layer.");
143 336 : params.addParam<std::vector<Real>>("theta_2",
144 : "The curve fit coefficient for "
145 : "GQ/H model"
146 : "for each soil layer.");
147 336 : params.addParam<std::vector<Real>>("theta_3",
148 : "The curve fit coefficient for "
149 : "GQ/H model"
150 : "for each soil layer.");
151 336 : params.addParam<std::vector<Real>>("theta_4",
152 : "The curve fit coefficient for "
153 : "GQ/H model"
154 : "for each soil layer.");
155 336 : params.addParam<std::vector<Real>>("theta_5",
156 : "The curve fit coefficient for "
157 : "GQ/H model"
158 : "for each soil layer.");
159 336 : params.addParam<std::vector<Real>>("taumax",
160 : "The ultimate shear strength of "
161 : "the soil layers. Required for "
162 : "GQ/H model");
163 : // params required for thin_layer contact model soil_type = "thin_layer"
164 336 : params.addParam<std::vector<Real>>("friction_coefficient",
165 : "Friction coefficients of the thin layers.");
166 336 : params.addParam<std::vector<Real>>("hardening_ratio",
167 : "Post-yield hardening ratios of the layers.");
168 168 : return params;
169 168 : }
170 :
171 126 : ADComputeISoilStress::ADComputeISoilStress(const InputParameters & parameters)
172 : : ADComputeFiniteStrainElasticStress(parameters),
173 : _base_models(),
174 : _stress_model(),
175 : _stress_model_old(),
176 : _yield_stress(), // *** YIELD STRAIN NOT YIELD STRESS ***
177 : _youngs(),
178 126 : _soil_layer_variable(coupledValue("layer_variable")),
179 252 : _layer_ids(getParam<std::vector<unsigned int>>("layer_ids")),
180 252 : _wave_speed_calculation(getParam<bool>("wave_speed_calculation")),
181 252 : _poissons_ratio(getParam<std::vector<Real>>("poissons_ratio")),
182 252 : _density(_wave_speed_calculation ? &getADMaterialProperty<Real>("density") : nullptr),
183 252 : _b_exp(getParam<Real>("b_exp")),
184 252 : _p_ref(getParam<std::vector<Real>>("p_ref")),
185 252 : _a0(getParam<Real>("a0")),
186 252 : _a1(getParam<Real>("a1")),
187 252 : _a2(getParam<Real>("a2")),
188 252 : _p0(getParam<Real>("tension_pressure_cut_off")),
189 252 : _pressure_dependency(getParam<bool>("pressure_dependency")),
190 126 : _strength_pressure_correction(1.0),
191 126 : _stiffness_pressure_correction(1.0),
192 126 : _shear_wave_speed(_wave_speed_calculation ? &declareADProperty<Real>("shear_wave_speed")
193 : : nullptr),
194 126 : _P_wave_speed(_wave_speed_calculation ? &declareADProperty<Real>("P_wave_speed") : nullptr),
195 126 : _pos(0),
196 126 : _initial_soil_stress_provided(
197 252 : getParam<std::vector<FunctionName>>("initial_soil_stress").size() ==
198 252 : LIBMESH_DIM * LIBMESH_DIM)
199 : {
200 : // checking that density, and Poisson's ratio are the same size as layer_ids
201 126 : if (_poissons_ratio.size() != _layer_ids.size())
202 0 : mooseError("Error in ADComputeISoilStress. Poisson's ratio should be of the same "
203 : "size as layer_ids.");
204 : // checks for pressure dependency
205 126 : if (_pressure_dependency && _b_exp == 0.0)
206 0 : mooseWarning("Warning in ADComputeISoilStress. Pressure dependency is set to true "
207 : "but b_exp is set to 0.0. Stiffness "
208 : "pressure dependency is NOT "
209 : "turned on.");
210 126 : if (_pressure_dependency && (_a0 == 1.0 && _a1 == 0.0 && _a2 == 0.0))
211 0 : mooseWarning(
212 : "Warning in ADComputeISoilStress. Pressure dependency is set to true but a0, a1 and a2 are "
213 : "set to 1.0, 0.0 and 0.0, respectively. Strength "
214 : "pressure dependency is NOT turned on.");
215 126 : if (_pressure_dependency && _p_ref.size() != _layer_ids.size())
216 0 : mooseError("Error in ADComputeISoilStress. When pressure dependency is turned on, "
217 : "a positive reference pressure "
218 : "(compressive) has to be defined for all "
219 : "the soil layers and the same number of reference "
220 : "pressures as soil layers should be provided.");
221 126 : if (_pressure_dependency && MastodonUtils::isNegativeOrZero(_p_ref))
222 0 : mooseError("Error in ADComputeISoilStress. Please provide positive (compressive) values for "
223 : "reference pressure.");
224 :
225 : // Initializing backbone curve
226 252 : const MooseEnum & soil_type = getParam<MooseEnum>("soil_type");
227 126 : std::vector<std::vector<Real>> backbone_stress(_layer_ids.size());
228 126 : std::vector<std::vector<Real>> backbone_strain(_layer_ids.size());
229 : // Calculating backbone curve for soil_type = user_defined
230 126 : if (soil_type == "user_defined")
231 : {
232 : std::vector<FileName> backbone_curve_files =
233 189 : getParam<std::vector<FileName>>("backbone_curve_files");
234 63 : if (backbone_curve_files.size() != _layer_ids.size())
235 0 : mooseError("Error in ADComputeISoilStress. A vector of file names needs to "
236 : "be provided for `backbone_curve_files` and the size of this vector "
237 : "should be same as that of `layers_ids`.");
238 63 : ISoilUtils::computeUserDefinedBackbone(
239 : backbone_stress, backbone_strain, _layer_ids, backbone_curve_files, "ADComputeISoilStress");
240 63 : }
241 : // Calculating backbone curve for soil_type = darendeli
242 63 : else if (soil_type == "darendeli")
243 : {
244 : const std::vector<Real> initial_shear_modulus =
245 36 : getParam<std::vector<Real>>("initial_shear_modulus");
246 : const std::vector<Real> over_consolidation_ratio =
247 36 : getParam<std::vector<Real>>("over_consolidation_ratio");
248 36 : const std::vector<Real> plasticity_index = getParam<std::vector<Real>>("plasticity_index");
249 36 : unsigned int number_of_points = getParam<unsigned int>("number_of_points");
250 18 : if (initial_shear_modulus.size() != _layer_ids.size() ||
251 36 : over_consolidation_ratio.size() != _layer_ids.size() ||
252 : plasticity_index.size() != _layer_ids.size())
253 0 : mooseError("Error in ADComputeISoilStress. initial_shear_modulus, "
254 : "over_consolidation_ratio and plasticity_index must be "
255 : "of the same size as layer_ids.");
256 36 : if (MastodonUtils::isNegativeOrZero(over_consolidation_ratio) || number_of_points <= 0 ||
257 18 : MastodonUtils::isNegativeOrZero(initial_shear_modulus))
258 0 : mooseError("Error in ADComputeISoilStress. Positive values have to be provided "
259 : "for over_consolidation_ratio, "
260 : "number_of_points, "
261 : "and initial_shear_modulus.");
262 18 : ISoilUtils::computeDarendeliBackbone(backbone_stress,
263 : backbone_strain,
264 : _layer_ids,
265 : initial_shear_modulus,
266 : over_consolidation_ratio,
267 : plasticity_index,
268 : _p_ref,
269 : number_of_points,
270 : "ADComputeISoilStress");
271 18 : }
272 : // Calculating backbone curve for soil_type = gqh
273 45 : else if (soil_type == "gqh")
274 : {
275 : const std::vector<Real> initial_shear_modulus =
276 54 : getParam<std::vector<Real>>("initial_shear_modulus");
277 54 : unsigned int number_of_points = getParam<unsigned int>("number_of_points");
278 54 : const std::vector<Real> theta_1 = getParam<std::vector<Real>>("theta_1");
279 54 : const std::vector<Real> theta_2 = getParam<std::vector<Real>>("theta_2");
280 54 : const std::vector<Real> theta_3 = getParam<std::vector<Real>>("theta_3");
281 54 : const std::vector<Real> theta_4 = getParam<std::vector<Real>>("theta_4");
282 54 : const std::vector<Real> theta_5 = getParam<std::vector<Real>>("theta_5");
283 81 : const std::vector<Real> taumax = getParam<std::vector<Real>>("taumax");
284 27 : const std::vector<std::vector<Real>> modulus(_layer_ids.size());
285 27 : if (initial_shear_modulus.size() != _layer_ids.size() || theta_1.size() != _layer_ids.size() ||
286 27 : theta_2.size() != _layer_ids.size() || theta_3.size() != _layer_ids.size() ||
287 54 : theta_4.size() != _layer_ids.size() || theta_5.size() != _layer_ids.size() ||
288 : taumax.size() != _layer_ids.size())
289 0 : mooseError("Error in ADComputeISoilStress. initial_shear_modulus, theta_1, "
290 : "theta_2, theta_3, theta_4, theta_5 and taumax should be of the "
291 : "same size as layer_ids.");
292 54 : if (number_of_points <= 0 || MastodonUtils::isNegativeOrZero(taumax) ||
293 27 : MastodonUtils::isNegativeOrZero(initial_shear_modulus))
294 0 : mooseError(
295 : "Error in ADComputeISoilStress. Please provide positive values for number of points, "
296 : "taumax and initial_shear_modulus.");
297 27 : ISoilUtils::computeGQHBackbone(backbone_stress,
298 : backbone_strain,
299 : _layer_ids,
300 : initial_shear_modulus,
301 : number_of_points,
302 : theta_1,
303 : theta_2,
304 : theta_3,
305 : theta_4,
306 : theta_5,
307 : taumax);
308 27 : }
309 : // Calculating backbone curve for soil_type = thin_layer
310 18 : else if (soil_type == "thin_layer")
311 : {
312 18 : _pressure_dependency = true;
313 18 : _a0 = 0.0;
314 18 : _a1 = 0.0;
315 18 : _a2 = 1.0;
316 36 : std::vector<Real> initial_shear_modulus = getParam<std::vector<Real>>("initial_shear_modulus");
317 36 : std::vector<Real> friction_coefficient = getParam<std::vector<Real>>("friction_coefficient");
318 54 : std::vector<Real> hardening_ratio = getParam<std::vector<Real>>("hardening_ratio");
319 18 : if (initial_shear_modulus.size() != _layer_ids.size() ||
320 36 : friction_coefficient.size() != _layer_ids.size() ||
321 : hardening_ratio.size() != _layer_ids.size())
322 0 : mooseError("Error in ADComputeISoilStress. initial_shear_modulus, friction_coefficient, "
323 : "hardening_ratio and p_ref must be "
324 : "of the same size as layer_ids.");
325 36 : if (MastodonUtils::isNegativeOrZero(initial_shear_modulus) ||
326 36 : MastodonUtils::isNegativeOrZero(friction_coefficient) ||
327 18 : MastodonUtils::isNegativeOrZero(hardening_ratio))
328 0 : mooseError("Error in ADComputeISoilStress. Positive values have to be provided "
329 : "for `initial_shear_modulus`, "
330 : "`friction_coefficient` and `hardening_ratio`");
331 18 : ISoilUtils::computeCoulombBackbone(backbone_stress,
332 : backbone_strain,
333 : _layer_ids,
334 : initial_shear_modulus,
335 : friction_coefficient,
336 : hardening_ratio,
337 : _p_ref,
338 : "ADComputeISoilStress");
339 18 : }
340 : else
341 0 : mooseError("Error in ADComputeISoilStress. The parameter soil_type is invalid.");
342 : // Deconstructing the backbone curves for all the soil layers into
343 : // elastic-perfectly-plastic components. Each backbone curve is split up into
344 : // a set of youngs modulus and yield stress pairs.
345 126 : _youngs.resize(_layer_ids.size());
346 126 : _yield_stress.resize(_layer_ids.size());
347 126 : computeSoilLayerProperties(_youngs,
348 : _yield_stress, // *** CALCULATES YIELD STRAIN NOT YIELD STRESS ***
349 : backbone_stress,
350 : backbone_strain,
351 : _layer_ids,
352 : _poissons_ratio,
353 : "ADComputeISoilStress");
354 126 : _stress_model.resize(_youngs[0].size());
355 126 : _stress_model_old.resize(_youngs[0].size());
356 126 : _base_models.resize(_youngs[0].size());
357 1836 : for (std::size_t i = 0; i < _youngs[0].size(); i++)
358 : {
359 1710 : _base_models[i] = Moose::stringify(i);
360 1710 : _stress_model[i] = &declareADProperty<RankTwoTensor>(_base_models[i] + "_stress_model");
361 1710 : _stress_model_old[i] =
362 3420 : &getMaterialPropertyOldByName<RankTwoTensor>(_base_models[i] + "_stress_model");
363 : }
364 :
365 : const std::vector<FunctionName> & fcn_names(
366 252 : getParam<std::vector<FunctionName>>("initial_soil_stress"));
367 126 : const unsigned num = fcn_names.size();
368 :
369 126 : if (!(num == 0 || num == LIBMESH_DIM * LIBMESH_DIM))
370 0 : mooseError("Either zero or ",
371 0 : LIBMESH_DIM * LIBMESH_DIM,
372 : " initial soil stress functions must be provided. You supplied ",
373 : num,
374 : "\n");
375 :
376 126 : _initial_soil_stress.resize(num);
377 1017 : for (unsigned i = 0; i < num; ++i)
378 891 : _initial_soil_stress[i] = &getFunctionByName(fcn_names[i]);
379 :
380 126 : _stress_new.zero();
381 126 : _individual_stress_increment.zero();
382 126 : _deviatoric_trial_stress.zero();
383 126 : }
384 :
385 : void
386 1088 : ADComputeISoilStress::initQpStatefulProperties()
387 : {
388 1088 : ADComputeStressBase::initQpStatefulProperties();
389 1088 : if (_initial_soil_stress_provided)
390 : {
391 3840 : for (unsigned i = 0; i < LIBMESH_DIM; ++i)
392 11520 : for (unsigned j = 0; j < LIBMESH_DIM; ++j)
393 8640 : _stress[_qp](i, j) = _initial_soil_stress[i * LIBMESH_DIM + j]->value(_t, _q_point[_qp]);
394 : }
395 68032 : for (std::size_t i = 0; i < _base_models.size(); i++)
396 66944 : (*_stress_model[i])[_qp].zero();
397 :
398 : // Determine the current id for the soil. The variable which is a Real must be
399 : // converted to a unsigned int for lookup, so first
400 : // it is rounded to avoid Real values that are just below the desired value.
401 1088 : _current_id = static_cast<unsigned int>(std::round(_soil_layer_variable[_qp]));
402 : // Get the position of the current id in the layer_ids array
403 1088 : _pos = find(_layer_ids.begin(), _layer_ids.end(), _current_id) - _layer_ids.begin();
404 :
405 1088 : if (_wave_speed_calculation)
406 : {
407 1088 : ADReal initial_youngs = 0.0;
408 68032 : for (std::size_t i = 0; i < _base_models.size(); i++)
409 66944 : initial_youngs += _youngs[_pos][i];
410 :
411 : // shear wave speed is sqrt(shear_modulus/density)
412 1088 : (*_shear_wave_speed)[_qp] =
413 2176 : std::sqrt(initial_youngs / (2.0 * (1.0 + _poissons_ratio[_pos])) / (*_density)[_qp]);
414 :
415 : // P wave speed is sqrt(P wave modulus/density)
416 1088 : (*_P_wave_speed)[_qp] =
417 2176 : std::sqrt(initial_youngs * (1.0 - _poissons_ratio[_pos]) / (1.0 + _poissons_ratio[_pos]) /
418 2176 : (1.0 - 2.0 * _poissons_ratio[_pos]) / (*_density)[_qp]);
419 : }
420 :
421 : // determine the lateral and vertical stresses
422 1088 : ADReal residual_vertical = _stress[_qp](2, 2);
423 1088 : ADReal residual_xx = _stress[_qp](0, 0);
424 1088 : ADReal residual_yy = _stress[_qp](1, 1);
425 1088 : ADReal mean_stress = _stress[_qp].trace() / (-3.0);
426 :
427 1088 : if (_pressure_dependency)
428 : {
429 384 : if (!MooseUtils::absoluteFuzzyEqual(mean_stress - _p0, 0.0))
430 384 : _stiffness_pressure_correction = std::pow((mean_stress - _p0) / _p_ref[_pos], _b_exp);
431 : else
432 0 : _stiffness_pressure_correction = 0.0;
433 :
434 384 : if (!MooseUtils::absoluteFuzzyEqual(
435 768 : _a0 + _a1 * (mean_stress - _p0) + _a2 * (mean_stress - _p0) * (mean_stress - _p0), 0.0))
436 : _strength_pressure_correction =
437 128 : std::sqrt(_a0 + _a1 * (mean_stress - _p0) +
438 0 : _a2 * (mean_stress - _p0) * (mean_stress - _p0)) /
439 256 : std::sqrt(_a0 + _a1 * (_p_ref[_pos]) + _a2 * (_p_ref[_pos]) * (_p_ref[_pos]));
440 : else
441 64 : _strength_pressure_correction = 0.0;
442 : }
443 :
444 : // Calculate the K0 consistent stress distribution
445 1088 : ADRankTwoTensor dev_model;
446 68032 : for (std::size_t i = 0; i < _base_models.size(); i++)
447 : {
448 66944 : ADReal _mean_pressure = 0.0;
449 66944 : if (residual_vertical != 0.0)
450 : {
451 66688 : ADReal sum_youngs = 0.0;
452 :
453 3312960 : for (std::size_t j = i; j < _base_models.size(); j++)
454 3246272 : sum_youngs += _youngs[_pos][j];
455 :
456 133376 : (*_stress_model[i])[_qp](2, 2) = residual_vertical * _youngs[_pos][i] / sum_youngs;
457 133376 : (*_stress_model[i])[_qp](0, 0) = residual_xx * _youngs[_pos][i] / sum_youngs;
458 133376 : (*_stress_model[i])[_qp](1, 1) = residual_yy * _youngs[_pos][i] / sum_youngs;
459 66688 : dev_model = ((*_stress_model[i])[_qp]).deviatoric() / _youngs[_pos][i];
460 133376 : _mean_pressure = (*_stress_model[i])[_qp].trace() / 3.0;
461 66688 : ADReal J2_model = dev_model.doubleContraction(dev_model);
462 66688 : ADReal dev_stress_model = 0.0;
463 66688 : if (!MooseUtils::absoluteFuzzyEqual(J2_model, 0.0))
464 133376 : dev_stress_model = std::sqrt(3.0 / 2.0 * J2_model);
465 133376 : if (dev_stress_model > _yield_stress[_pos][i] * _strength_pressure_correction)
466 87296 : dev_model *= (_yield_stress[_pos][i] * _strength_pressure_correction) / dev_stress_model;
467 :
468 66688 : (*_stress_model[i])[_qp] = dev_model * _youngs[_pos][i]; // stress_model contains only the
469 : // deviatoric part of the stress
470 : }
471 133888 : residual_vertical = residual_vertical - (*_stress_model[i])[_qp](2, 2) - _mean_pressure;
472 133888 : residual_xx = residual_xx - (*_stress_model[i])[_qp](0, 0) - _mean_pressure;
473 133888 : residual_yy = residual_yy - (*_stress_model[i])[_qp](1, 1) - _mean_pressure;
474 : }
475 1088 : }
476 :
477 : void
478 956432 : ADComputeISoilStress::computeQpStress()
479 : {
480 : // Nothing to update during the first time step, return immediately
481 956432 : if (_t_step == 0)
482 : return;
483 :
484 955888 : _stress_new.zero();
485 955888 : computeStress();
486 955888 : _stress[_qp] = _rotation_increment[_qp] * _stress_new * _rotation_increment[_qp].transpose();
487 : }
488 :
489 : void
490 955888 : ADComputeISoilStress::computeStress()
491 : {
492 955888 : if (_t_step == 0)
493 0 : return;
494 :
495 : // Determine the current id for the soil. The variable which is a Real must be
496 : // converted to a unsigned int for lookup, so first
497 : // it is rounded to avoid Real values that are just below the desired value.
498 955888 : _current_id = static_cast<unsigned int>(std::round(_soil_layer_variable[_qp]));
499 :
500 : // Get the position of the current id in the layer_ids array
501 955888 : _pos = find(_layer_ids.begin(), _layer_ids.end(), _current_id) - _layer_ids.begin();
502 :
503 955888 : _individual_stress_increment.zero();
504 :
505 955888 : _deviatoric_trial_stress.zero();
506 :
507 : // current pressure calculation
508 955888 : ADReal mean_stress = _stress_old[_qp].trace() / (-3.0);
509 :
510 955888 : if (mean_stress < _p0)
511 17280 : mean_stress = 0.0;
512 :
513 955888 : if (_pressure_dependency)
514 : {
515 689472 : if (!MooseUtils::absoluteFuzzyEqual(mean_stress - _p0, 0.0))
516 689472 : _stiffness_pressure_correction = std::pow((mean_stress - _p0) / _p_ref[_pos], _b_exp);
517 : else
518 0 : _stiffness_pressure_correction = 0.0;
519 :
520 689472 : if (!MooseUtils::absoluteFuzzyEqual(
521 1378944 : _a0 + _a1 * (mean_stress - _p0) + _a2 * (mean_stress - _p0) * (mean_stress - _p0), 0.0))
522 : _strength_pressure_correction =
523 315872 : std::sqrt(_a0 + _a1 * (mean_stress - _p0) +
524 315872 : _a2 * (mean_stress - _p0) * (mean_stress - _p0)) /
525 631744 : std::sqrt(_a0 + _a1 * (_p_ref[_pos]) + _a2 * (_p_ref[_pos]) * (_p_ref[_pos]));
526 : else
527 28864 : _strength_pressure_correction = 0.0;
528 : }
529 :
530 955888 : _mean_pressure = 0.0;
531 :
532 : // compute trial stress increment - note that _elasticity_tensor here
533 : // assumes youngs modulus = 1
534 :
535 955888 : ADReal lame_1 = _elasticity_tensor[_qp](0, 0, 1, 1);
536 1911776 : ADReal lame_2 = (_elasticity_tensor[_qp](0, 0, 0, 0) - lame_1) / 2;
537 1911776 : ADReal initial_youngs = (lame_2 * (3 * lame_1 + 2 * lame_2)) / (lame_1 + lame_2);
538 :
539 : _individual_stress_increment =
540 955888 : _elasticity_tensor[_qp] * (_strain_increment[_qp]) / initial_youngs;
541 :
542 1911776 : _mean_pressure_tmp = _individual_stress_increment.trace() / 3.0 * _stiffness_pressure_correction;
543 :
544 : _deviatoric_trial_stress_tmp =
545 955888 : _individual_stress_increment.deviatoric() * _stiffness_pressure_correction;
546 :
547 955888 : _dev_trial_stress_squared = 0.0;
548 :
549 955888 : _effective_trial_stress = 0.0;
550 :
551 955888 : _yield_condition = 0.0;
552 :
553 8273232 : for (std::size_t i = 0; i < _base_models.size(); i++)
554 : {
555 :
556 : // calculate pressure for each element
557 14634688 : _mean_pressure += _mean_pressure_tmp * _youngs[_pos][i];
558 :
559 : // compute the deviatoric trial stress normalized by non pressure dependent
560 : // youngs modulus.
561 : _deviatoric_trial_stress =
562 7317344 : _deviatoric_trial_stress_tmp + (*_stress_model_old[i])[_qp] / (_youngs[_pos][i]);
563 :
564 : // compute the effective trial stress
565 : _dev_trial_stress_squared =
566 7317344 : _deviatoric_trial_stress.doubleContraction(_deviatoric_trial_stress);
567 :
568 7317344 : if (!MooseUtils::absoluteFuzzyEqual(_dev_trial_stress_squared, 0.0))
569 13976896 : _effective_trial_stress = std::sqrt(3.0 / 2.0 * _dev_trial_stress_squared);
570 :
571 : // check yield condition and calculate plastic strain
572 : _yield_condition =
573 14634688 : _effective_trial_stress - _yield_stress[_pos][i] * _strength_pressure_correction;
574 :
575 7317344 : if (_yield_condition > 0.0)
576 : _deviatoric_trial_stress *=
577 6694232 : _yield_stress[_pos][i] * _strength_pressure_correction / _effective_trial_stress;
578 :
579 14634688 : (*_stress_model[i])[_qp] = _youngs[_pos][i] * (_deviatoric_trial_stress);
580 :
581 7317344 : _stress_new += (*_stress_model[i])[_qp];
582 : }
583 955888 : _stress_new(0, 0) += _mean_pressure - mean_stress;
584 955888 : _stress_new(1, 1) += _mean_pressure - mean_stress;
585 955888 : _stress_new(2, 2) += _mean_pressure - mean_stress;
586 : }
587 :
588 : void
589 126 : ADComputeISoilStress::computeSoilLayerProperties(
590 : std::vector<std::vector<ADReal>> & youngs,
591 : std::vector<std::vector<ADReal>> & yield_stress, // *** YIELD STRAIN, NOT YIELD STRESS***
592 : const std::vector<std::vector<Real>> & backbone_stress,
593 : const std::vector<std::vector<Real>> & backbone_strain,
594 : const std::vector<unsigned int> & layer_ids,
595 : const std::vector<Real> & poissons_ratio,
596 : const std::string & name)
597 : {
598 423 : for (std::size_t k = 0; k < layer_ids.size(); k++)
599 : {
600 297 : unsigned int number = backbone_strain[k].size();
601 :
602 : // Calculating the Youngs modulus for each of the n curves for soil layer k
603 297 : std::vector<std::vector<ADReal>> A(number);
604 297 : std::vector<std::vector<ADReal>> InvA(number);
605 :
606 19107 : for (std::size_t i = 0; i < number; i++)
607 : {
608 18810 : A[i].resize(number);
609 1826190 : for (std::size_t j = 0; j < number; j++)
610 : {
611 1807380 : if (i <= j)
612 913095 : A[i][j] = backbone_strain[k][i];
613 : else
614 894285 : A[i][j] = backbone_strain[k][j];
615 : }
616 : }
617 :
618 : // create upper triangular matrix and modified stress
619 297 : youngs[k].resize(number);
620 297 : std::vector<ADReal> G0_component(number);
621 297 : yield_stress[k].resize(number);
622 297 : std::vector<ADReal> modified_backbone_stress(number);
623 297 : InvA = A;
624 18810 : for (std::size_t i = 1; i < number; i++)
625 : {
626 1807083 : for (std::size_t j = 0; j < number; j++)
627 : {
628 3577140 : InvA[i][j] = A[i][j] - A[i - 1][j];
629 1788570 : modified_backbone_stress[i] = backbone_stress[k][i] - backbone_stress[k][i - 1];
630 : }
631 : }
632 297 : modified_backbone_stress[0] = backbone_stress[k][0];
633 :
634 : // backward substitution
635 594 : G0_component[number - 1] = modified_backbone_stress[number - 1] / InvA[number - 1][number - 1];
636 297 : yield_stress[k][number - 1] = backbone_strain[k][number - 1];
637 297 : int i = number - 2;
638 18810 : while (i >= 0)
639 : {
640 18513 : G0_component[i] = 0.0;
641 912798 : for (std::size_t j = i + 1; j < number; j++)
642 894285 : G0_component[i] += InvA[i][j] * G0_component[j];
643 :
644 18513 : G0_component[i] = (modified_backbone_stress[i] - G0_component[i]) / InvA[i][i];
645 18513 : yield_stress[k][i] = backbone_strain[k][i];
646 18513 : i = i - 1;
647 : }
648 :
649 : // scaling
650 19107 : for (std::size_t i = 0; i < number; i++)
651 : {
652 37620 : youngs[k][i] = G0_component[i] * 2.0 * (1.0 + poissons_ratio[k]);
653 37620 : yield_stress[k][i] = yield_stress[k][i] * std::sqrt(3.0) / (2.0 * (1.0 + poissons_ratio[k]));
654 : }
655 297 : }
656 423 : for (std::size_t k = 0; k < layer_ids.size(); k++)
657 : {
658 19107 : for (std::size_t j = 0; j < backbone_strain[k].size(); j++)
659 : {
660 18810 : if (youngs[k][j] == 0)
661 0 : mooseError("In " + name + "Found two segments in the backbone curve with the same slope. ",
662 : "Please recheck input backbone curves.");
663 : }
664 : }
665 126 : }
|
