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 "ComputeISoilStress.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 :
27 : registerMooseObject("MastodonApp", ComputeISoilStress);
28 :
29 : InputParameters
30 302 : ComputeISoilStress::validParams()
31 : {
32 302 : InputParameters params = ComputeFiniteStrainElasticStress::validParams();
33 302 : params.addClassDescription("Compute total stress for the nonlinear material "
34 : "model I-Soil using a backbone curve.");
35 604 : params.addRequiredCoupledVar("layer_variable",
36 : "The auxvariable providing the soil layer identification.");
37 604 : params.addRequiredParam<std::vector<unsigned int>>(
38 : "layer_ids",
39 : "Vector of layer ids that map one-to-one to the rest of the "
40 : "soil layer parameters provided as input.");
41 604 : params.addRequiredParam<std::vector<Real>>("poissons_ratio",
42 : "Poissons's ratio for the soil layers. The "
43 : "size of the vector should be same as the size of "
44 : "layer_ids.");
45 604 : params.addParam<Real>("b_exp",
46 604 : 0.0,
47 : "The exponential factors for pressure "
48 : "dependent stiffness for all the soil "
49 : "layers.");
50 604 : params.addParam<std::vector<Real>>("p_ref",
51 : {},
52 : "The reference pressure at which "
53 : "the parameters are defined for "
54 : "each soil layer. If 'soil_type = "
55 : "darendeli', then the reference "
56 : "pressure must be input in kilopascals.");
57 604 : params.addParam<Real>("a0",
58 604 : 1.0,
59 : "The first coefficient for pressure dependent yield strength "
60 : "calculation for all the soil layers. If a0 = 1, a1 = 0 and "
61 : "a2=0 for one soil layer, then the yield strength of that "
62 : "layer is independent of pressure.");
63 604 : params.addParam<Real>("a1",
64 604 : 0.0,
65 : "The second coefficient for pressure dependent yield "
66 : "strength calculation for all the soil layers. If a0 = "
67 : "1, a1 = 0, a2 = 0 for one soil layer, then the yield "
68 : "strength of that layer is independent of pressure.");
69 604 : params.addParam<Real>("a2",
70 604 : 0.0,
71 : "The third coefficient for pressure dependent yield "
72 : "strength calculation for all the soil layers. If a0 = "
73 : "1, a1=0 and a2=0 for one soil layer, then the yield "
74 : "strength of that layer is independent of pressure.");
75 604 : params.addParam<Real>("tension_pressure_cut_off",
76 604 : -1.0,
77 : "The tension cut-off for all the soil layers. If the "
78 : "pressure becomes lower than this value, then the "
79 : "stiffness of the soil reduces to zero. A negative "
80 : "pressure indicates tension. The default "
81 : "value is -1.0 for all the soil layers.");
82 604 : params.addParam<bool>("pressure_dependency",
83 604 : false,
84 : "Set to true to turn on pressure dependent stiffness "
85 : "and yield strength calculation.");
86 604 : params.addParam<bool>(
87 604 : "wave_speed_calculation", true, "Set to false to turn off P and S wave speed calculation.");
88 604 : params.addParam<std::vector<FunctionName>>(
89 : "initial_soil_stress",
90 : {},
91 : "A list of functions describing the initial stress. There "
92 : "must be 9 functions, corresponding to the xx, yx, zx, xy, yy, zy, xz, yz, "
93 : "zz components respectively. If not provided, all components of the "
94 : "initial stress will be zero.");
95 : // params for specific backbone types
96 604 : MooseEnum soil_type("user_defined darendeli gqh thin_layer");
97 604 : params.addRequiredParam<MooseEnum>(
98 : "soil_type",
99 : soil_type,
100 : "This parameter determines the type of backbone curve used. Use 'user_defined' "
101 : "for a user defined backbone curve provided in a data file, 'darendeli' "
102 : "if the backbone curve is to be determined using Darandeli equations, 'gqh' "
103 : "if the backbone curve is determined using the GQ/H approach and 'thin_layer' if the soil is "
104 : "being used to simulate a thin-layer friction interface.");
105 : // params required for user_defined backbone curve: soil_type = 'user_defined'
106 604 : params.addParam<std::vector<FileName>>(
107 : "backbone_curve_files",
108 : "The vector of file names of the files containing "
109 : "stress-strain backbone curves for the different soil layers. The "
110 : "size of the vector should be same as the size of layer_ids. All files "
111 : "should contain the same number of stress-strain points. Headers are not "
112 : "expected and it is assumed that the first column corresponds to strain values "
113 : "and the second column corresponds to the stress values. Additionally, two "
114 : "segments of a backbone curve cannot have the same slope.");
115 : // params required for soil_type = 'darendeli', 'GQ/H' and 'thin_layer'
116 604 : params.addRequiredParam<std::vector<Real>>("initial_shear_modulus",
117 : "The initial shear modulus of the soil layers. ");
118 604 : params.addParam<MaterialPropertyName>("shear_modulus",
119 : "Name of Material Property or a constant real number "
120 : "defining the shear modulus of the materials.");
121 : // params required for soil_type = 'darendeli' and 'GQ/H'
122 604 : params.addParam<unsigned int>("number_of_points",
123 : "The total number of data points in which the "
124 : "backbone curve needs to be split for all soil "
125 : "layers (required for Darandeli or GQ/H type backbone curves).");
126 : // params required for Darandeli backbone curve: soil_type = 'darendeli'
127 604 : params.addParam<std::vector<Real>>("over_consolidation_ratio",
128 : "The over consolidation ratio of the soil "
129 : "layers. Required for Darandeli backbone curve.");
130 604 : params.addParam<std::vector<Real>>(
131 : "plasticity_index",
132 : "The plasticity index of the soil layers. Required for Darandeli backbone curve.");
133 : // params required for GQ/H backbone curve: soil_type = "gqh"
134 604 : params.addParam<std::vector<Real>>("theta_1",
135 : "The curve fit coefficient for "
136 : "GQ/H model"
137 : "for each soil layer.");
138 604 : params.addParam<std::vector<Real>>("theta_2",
139 : "The curve fit coefficient for "
140 : "GQ/H model"
141 : "for each soil layer.");
142 604 : params.addParam<std::vector<Real>>("theta_3",
143 : "The curve fit coefficient for "
144 : "GQ/H model"
145 : "for each soil layer.");
146 604 : params.addParam<std::vector<Real>>("theta_4",
147 : "The curve fit coefficient for "
148 : "GQ/H model"
149 : "for each soil layer.");
150 604 : params.addParam<std::vector<Real>>("theta_5",
151 : "The curve fit coefficient for "
152 : "GQ/H model"
153 : "for each soil layer.");
154 604 : params.addParam<std::vector<Real>>("taumax",
155 : "The ultimate shear strength of "
156 : "the soil layers. Required for "
157 : "GQ/H model");
158 : // params required for thin_layer contact model soil_type = "thin_layer"
159 604 : params.addParam<std::vector<Real>>("friction_coefficient",
160 : "Friction coefficients of the thin layers.");
161 604 : params.addParam<std::vector<Real>>("hardening_ratio",
162 : "Post-yield hardening ratios of the layers.");
163 302 : return params;
164 302 : }
165 :
166 227 : ComputeISoilStress::ComputeISoilStress(const InputParameters & parameters)
167 : : ComputeFiniteStrainElasticStress(parameters),
168 227 : _strain_increment(getMaterialProperty<RankTwoTensor>(_base_name + "strain_increment")),
169 : _base_models(),
170 : _stress_model(),
171 : _stress_model_old(),
172 : _yield_stress(), // *** YIELD STRAIN NOT YIELD STRESS ***
173 : _youngs(),
174 227 : _soil_layer_variable(coupledValue("layer_variable")),
175 454 : _layer_ids(getParam<std::vector<unsigned int>>("layer_ids")),
176 454 : _wave_speed_calculation(getParam<bool>("wave_speed_calculation")),
177 454 : _poissons_ratio(getParam<std::vector<Real>>("poissons_ratio")),
178 454 : _density(_wave_speed_calculation ? &getMaterialProperty<Real>("density") : nullptr),
179 454 : _b_exp(getParam<Real>("b_exp")),
180 454 : _p_ref(getParam<std::vector<Real>>("p_ref")),
181 454 : _a0(getParam<Real>("a0")),
182 454 : _a1(getParam<Real>("a1")),
183 454 : _a2(getParam<Real>("a2")),
184 454 : _p0(getParam<Real>("tension_pressure_cut_off")),
185 454 : _pressure_dependency(getParam<bool>("pressure_dependency")),
186 227 : _strength_pressure_correction(1.0),
187 227 : _stiffness_pressure_correction(1.0),
188 227 : _shear_wave_speed(_wave_speed_calculation ? &declareProperty<Real>("shear_wave_speed")
189 : : nullptr),
190 227 : _P_wave_speed(_wave_speed_calculation ? &declareProperty<Real>("P_wave_speed") : nullptr),
191 227 : _tangent_modulus(0.0),
192 227 : _pos(0),
193 227 : _initial_soil_stress_provided(
194 454 : getParam<std::vector<FunctionName>>("initial_soil_stress").size() ==
195 454 : LIBMESH_DIM * LIBMESH_DIM)
196 : {
197 :
198 : // checking that density, and Poisson's ratio are the same size as layer_ids
199 227 : if (_poissons_ratio.size() != _layer_ids.size())
200 0 : mooseError("Error in " + name() + ". Poisson's ratio should be of the same "
201 : "size as layer_ids.");
202 : // checks for pressure dependency
203 227 : if (_pressure_dependency && _b_exp == 0.0)
204 29 : mooseWarning("Warning in " + name() + ". Pressure dependency is set to true "
205 : "but b_exp is set to 0.0. Stiffness "
206 : "pressure dependency is NOT "
207 : "turned on.");
208 226 : if (_pressure_dependency && (_a0 == 1.0 && _a1 == 0.0 && _a2 == 0.0))
209 2 : mooseWarning("Warning in " + name() +
210 : ". Pressure dependency is set to true but a0, a1 and a2 are "
211 : "set to 1.0, 0.0 and 0.0, respectively. Strength "
212 : "pressure dependency is NOT turned on.");
213 225 : if (_pressure_dependency && (_a0 == 0.0 && _a1 == 0.0 && _a2 == 0.0))
214 0 : mooseError("Error in " + name() +
215 : ". When pressure dependency is turned on, "
216 : "all three strength coefficients, a0, a1, and a2, "
217 : "cannot simultaneously be set to 0.0. This "
218 : "combination results in division by 0.");
219 225 : if (_pressure_dependency && _p_ref.size() != _layer_ids.size())
220 0 : mooseError("Error in " + name() + ". When pressure dependency is turned on, "
221 : "a positive reference pressure "
222 : "(compressive) has to be defined for all "
223 : "the soil layers and the same number of reference "
224 : "pressures as soil layers should be provided.");
225 225 : if (_pressure_dependency && MastodonUtils::isNegativeOrZero(_p_ref))
226 0 : mooseError("Error in " + name() +
227 : ". Please provide positive (compressive) values for reference pressure.");
228 :
229 : // Initializing backbone curve
230 450 : const MooseEnum & soil_type = getParam<MooseEnum>("soil_type");
231 225 : std::vector<std::vector<Real>> backbone_stress(_layer_ids.size());
232 225 : std::vector<std::vector<Real>> backbone_strain(_layer_ids.size());
233 450 : std::vector<Real> initial_shear_modulus = getParam<std::vector<Real>>("initial_shear_modulus");
234 :
235 : // Calculating backbone curve for soil_type = user_defined
236 225 : if (soil_type == "user_defined")
237 : {
238 : std::vector<FileName> backbone_curve_files =
239 216 : getParam<std::vector<FileName>>("backbone_curve_files");
240 72 : if (backbone_curve_files.size() != _layer_ids.size())
241 0 : mooseError("Error in " + name() +
242 : ". A vector of file names needs to "
243 : "be provided for `backbone_curve_files` and the size of this vector "
244 : "should be same as that of `layers_ids`.");
245 72 : ISoilUtils::computeUserDefinedBackbone(
246 : backbone_stress, backbone_strain, _layer_ids, backbone_curve_files, name());
247 72 : }
248 : // Calculating backbone curve for soil_type = darendeli
249 153 : else if (soil_type == "darendeli")
250 : {
251 : std::vector<Real> over_consolidation_ratio =
252 108 : getParam<std::vector<Real>>("over_consolidation_ratio");
253 108 : std::vector<Real> plasticity_index = getParam<std::vector<Real>>("plasticity_index");
254 108 : unsigned int number_of_points = getParam<unsigned int>("number_of_points");
255 54 : if (initial_shear_modulus.size() != _layer_ids.size() ||
256 108 : over_consolidation_ratio.size() != _layer_ids.size() ||
257 : plasticity_index.size() != _layer_ids.size())
258 0 : mooseError("Error in " + name() + ". initial_shear_modulus, "
259 : "over_consolidation_ratio and plasticity_index must be "
260 : "of the same size as layer_ids.");
261 108 : if (MastodonUtils::isNegativeOrZero(over_consolidation_ratio) || number_of_points <= 0 ||
262 54 : MastodonUtils::isNegativeOrZero(initial_shear_modulus))
263 0 : mooseError("Error in " + name() + ". Positive values have to be provided "
264 : "for over_consolidation_ratio, "
265 : "number_of_points, "
266 : "and initial_shear_modulus.");
267 54 : ISoilUtils::computeDarendeliBackbone(backbone_stress,
268 : backbone_strain,
269 : _layer_ids,
270 : initial_shear_modulus,
271 : over_consolidation_ratio,
272 : plasticity_index,
273 : _p_ref,
274 : number_of_points,
275 : name());
276 54 : }
277 : // Calculating backbone curve for soil_type = gqh
278 99 : else if (soil_type == "gqh")
279 : {
280 162 : unsigned int number_of_points = getParam<unsigned int>("number_of_points");
281 162 : std::vector<Real> theta_1 = getParam<std::vector<Real>>("theta_1");
282 162 : std::vector<Real> theta_2 = getParam<std::vector<Real>>("theta_2");
283 162 : std::vector<Real> theta_3 = getParam<std::vector<Real>>("theta_3");
284 162 : std::vector<Real> theta_4 = getParam<std::vector<Real>>("theta_4");
285 162 : std::vector<Real> theta_5 = getParam<std::vector<Real>>("theta_5");
286 243 : std::vector<Real> taumax = getParam<std::vector<Real>>("taumax");
287 81 : std::vector<std::vector<Real>> modulus(_layer_ids.size());
288 81 : if (initial_shear_modulus.size() != _layer_ids.size() || theta_1.size() != _layer_ids.size() ||
289 81 : theta_2.size() != _layer_ids.size() || theta_3.size() != _layer_ids.size() ||
290 162 : theta_4.size() != _layer_ids.size() || theta_5.size() != _layer_ids.size() ||
291 : taumax.size() != _layer_ids.size())
292 0 : mooseError("Error in " + name() +
293 : ". initial_shear_modulus, theta_1, "
294 : "theta_2, theta_3, theta_4, theta_5 and taumax should be of the "
295 : "same size as layer_ids.");
296 162 : if (number_of_points <= 0 || MastodonUtils::isNegativeOrZero(taumax) ||
297 81 : MastodonUtils::isNegativeOrZero(initial_shear_modulus))
298 0 : mooseError("Error in " + name() + ". Please provide positive values for number of points, "
299 : "taumax and initial_shear_modulus.");
300 81 : ISoilUtils::computeGQHBackbone(backbone_stress,
301 : backbone_strain,
302 : _layer_ids,
303 : initial_shear_modulus,
304 : number_of_points,
305 : theta_1,
306 : theta_2,
307 : theta_3,
308 : theta_4,
309 : theta_5,
310 : taumax);
311 81 : }
312 : // Calculating backbone curve for soil_type = thin_layer
313 18 : else if (soil_type == "thin_layer")
314 : {
315 18 : _pressure_dependency = true;
316 18 : _a0 = 0.0;
317 18 : _a1 = 0.0;
318 18 : _a2 = 1.0;
319 36 : std::vector<Real> friction_coefficient = getParam<std::vector<Real>>("friction_coefficient");
320 54 : std::vector<Real> hardening_ratio = getParam<std::vector<Real>>("hardening_ratio");
321 18 : if (initial_shear_modulus.size() != _layer_ids.size() ||
322 36 : friction_coefficient.size() != _layer_ids.size() ||
323 : hardening_ratio.size() != _layer_ids.size())
324 0 : mooseError("Error in " + name() + ". initial_shear_modulus, friction_coefficient, "
325 : "hardening_ratio and p_ref must be "
326 : "of the same size as layer_ids.");
327 36 : if (MastodonUtils::isNegativeOrZero(initial_shear_modulus) ||
328 36 : MastodonUtils::isNegativeOrZero(friction_coefficient) ||
329 18 : MastodonUtils::isNegativeOrZero(hardening_ratio))
330 0 : mooseError("Error in " + name() + ". Positive values have to be provided "
331 : "for `initial_shear_modulus`, "
332 : "`friction_coefficient` and `hardening_ratio`");
333 18 : ISoilUtils::computeCoulombBackbone(backbone_stress,
334 : backbone_strain,
335 : _layer_ids,
336 : initial_shear_modulus,
337 : friction_coefficient,
338 : hardening_ratio,
339 : _p_ref,
340 : name());
341 18 : }
342 : else
343 0 : mooseError("Error in " + name() + ". The parameter soil_type is invalid.");
344 : // Deconstructing the backbone curves for all the soil layers into
345 : // elastic-perfectly-plastic components. Each backbone curve is split up into
346 : // a set of youngs modulus and yield stress pairs.
347 225 : _youngs.resize(_layer_ids.size());
348 225 : _yield_stress.resize(_layer_ids.size());
349 225 : ISoilUtils::computeSoilLayerProperties(
350 : _youngs,
351 : _yield_stress, // *** CALCULATES YIELD STRAIN NOT YIELD STRESS ***
352 : backbone_stress,
353 : backbone_strain,
354 : _layer_ids,
355 : _poissons_ratio,
356 : name());
357 225 : _stress_model.resize(_youngs[0].size());
358 225 : _stress_model_old.resize(_youngs[0].size());
359 225 : _base_models.resize(_youngs[0].size());
360 7038 : for (std::size_t i = 0; i < _youngs[0].size(); i++)
361 : {
362 6813 : _base_models[i] = Moose::stringify(i);
363 6813 : _stress_model[i] = &declareProperty<RankTwoTensor>(_base_models[i] + "_stress_model");
364 6813 : _stress_model_old[i] =
365 13626 : &getMaterialPropertyOld<RankTwoTensor>(_base_models[i] + "_stress_model");
366 : }
367 :
368 : const std::vector<FunctionName> & fcn_names(
369 450 : getParam<std::vector<FunctionName>>("initial_soil_stress"));
370 225 : const unsigned num = fcn_names.size();
371 :
372 225 : if (!(num == 0 || num == LIBMESH_DIM * LIBMESH_DIM))
373 0 : mooseError("Either zero or ",
374 0 : LIBMESH_DIM * LIBMESH_DIM,
375 : " initial soil stress functions must be provided. You supplied ",
376 : num,
377 : "\n");
378 :
379 225 : _initial_soil_stress.resize(num);
380 1764 : for (unsigned i = 0; i < num; ++i)
381 1539 : _initial_soil_stress[i] = &getFunctionByName(fcn_names[i]);
382 :
383 : _stress_new.zero();
384 : _individual_stress_increment.zero();
385 : _deviatoric_trial_stress.zero();
386 :
387 : // // checking that the input and the backbone shear modulus values are consistent.
388 : std::vector<Real> initial_shear;
389 225 : initial_shear.resize(_poissons_ratio.size());
390 : Real tmp = 0.0;
391 810 : for (std::size_t j = 0; j < _poissons_ratio.size(); j++)
392 : {
393 : tmp = 0.0;
394 43398 : for (std::size_t i = 0; i < _base_models.size(); i++)
395 : {
396 42813 : tmp += _youngs[j][i] / (2 * (1 + _poissons_ratio[j]));
397 : }
398 585 : initial_shear[j] = tmp;
399 : }
400 :
401 225 : bool value_bool = MastodonUtils::checkEqual(initial_shear, initial_shear_modulus, 5.0);
402 225 : if (value_bool == false)
403 9 : mooseWarning(
404 : "Shear moduli inferred from the backbone curve are different from the input values."
405 : " Using the backbone curve inferred value for further computations.");
406 225 : }
407 :
408 : void
409 1888 : ComputeISoilStress::initQpStatefulProperties()
410 : {
411 1888 : ComputeStressBase::initQpStatefulProperties();
412 1888 : if (_initial_soil_stress_provided)
413 : {
414 7040 : for (unsigned i = 0; i < LIBMESH_DIM; ++i)
415 21120 : for (unsigned j = 0; j < LIBMESH_DIM; ++j)
416 15840 : _stress[_qp](i, j) = _initial_soil_stress[i * LIBMESH_DIM + j]->value(_t, _q_point[_qp]);
417 : }
418 152192 : for (std::size_t i = 0; i < _base_models.size(); i++)
419 150304 : (*_stress_model[i])[_qp].zero();
420 :
421 : // Determine the current id for the soil. The variable which is a Real must be
422 : // converted to a unsigned int for lookup, so first
423 : // it is rounded to avoid Real values that are just below the desired value.
424 1888 : _current_id = static_cast<unsigned int>(std::round(_soil_layer_variable[_qp]));
425 : // Get the position of the current id in the layer_ids array
426 1888 : _pos = find(_layer_ids.begin(), _layer_ids.end(), _current_id) - _layer_ids.begin();
427 :
428 1888 : if (_wave_speed_calculation)
429 : {
430 : Real initial_youngs = 0.0;
431 152192 : for (std::size_t i = 0; i < _base_models.size(); i++)
432 150304 : initial_youngs += _youngs[_pos][i];
433 :
434 : // shear wave speed is sqrt(shear_modulus/density)
435 1888 : (*_shear_wave_speed)[_qp] =
436 1888 : std::sqrt(initial_youngs / (2.0 * (1.0 + _poissons_ratio[_pos])) / (*_density)[_qp]);
437 :
438 : // P wave speed is sqrt(P wave modulus/density)
439 1888 : (*_P_wave_speed)[_qp] =
440 1888 : std::sqrt(initial_youngs * (1.0 - _poissons_ratio[_pos]) / (1.0 + _poissons_ratio[_pos]) /
441 1888 : (1.0 - 2.0 * _poissons_ratio[_pos]) / (*_density)[_qp]);
442 : }
443 :
444 : // determine the lateral and vertical stresses
445 1888 : Real residual_vertical = _stress[_qp](2, 2);
446 1888 : Real residual_xx = _stress[_qp](0, 0);
447 1888 : Real residual_yy = _stress[_qp](1, 1);
448 1888 : Real mean_stress = _stress[_qp].trace() / (-3.0);
449 :
450 1888 : if (_pressure_dependency)
451 : {
452 192 : _stiffness_pressure_correction = pow((mean_stress - _p0) / _p_ref[_pos], _b_exp);
453 192 : _strength_pressure_correction =
454 192 : std::sqrt(_a0 + _a1 * (mean_stress - _p0) +
455 192 : _a2 * (mean_stress - _p0) * (mean_stress - _p0)) /
456 192 : std::sqrt(_a0 + _a1 * (_p_ref[_pos]) + _a2 * (_p_ref[_pos]) * (_p_ref[_pos]));
457 : }
458 :
459 : // Calculate the K0 consistent stress distribution
460 1888 : RankTwoTensor dev_model;
461 152192 : for (std::size_t i = 0; i < _base_models.size(); i++)
462 : {
463 : Real _mean_pressure = 0.0;
464 150304 : if (residual_vertical != 0.0)
465 : {
466 : Real sum_youngs = 0.0;
467 :
468 7268224 : for (std::size_t j = i; j < _base_models.size(); j++)
469 7124704 : sum_youngs += _youngs[_pos][j];
470 :
471 143520 : (*_stress_model[i])[_qp](2, 2) = residual_vertical * _youngs[_pos][i] / sum_youngs;
472 143520 : (*_stress_model[i])[_qp](0, 0) = residual_xx * _youngs[_pos][i] / sum_youngs;
473 143520 : (*_stress_model[i])[_qp](1, 1) = residual_yy * _youngs[_pos][i] / sum_youngs;
474 143520 : dev_model = ((*_stress_model[i])[_qp]).deviatoric() / _youngs[_pos][i];
475 143520 : _mean_pressure = (*_stress_model[i])[_qp].trace() / 3.0;
476 143520 : Real J2_model = dev_model.doubleContraction(dev_model);
477 143520 : Real dev_stress_model = std::sqrt(3.0 / 2.0 * J2_model);
478 143520 : if (dev_stress_model > _yield_stress[_pos][i] * _strength_pressure_correction)
479 92336 : dev_model *= (_yield_stress[_pos][i] * _strength_pressure_correction) / dev_stress_model;
480 :
481 143520 : (*_stress_model[i])[_qp] = dev_model * _youngs[_pos][i]; // stress_model contains only the
482 : // deviatoric part of the stress
483 : }
484 150304 : residual_vertical = residual_vertical - (*_stress_model[i])[_qp](2, 2) - _mean_pressure;
485 150304 : residual_xx = residual_xx - (*_stress_model[i])[_qp](0, 0) - _mean_pressure;
486 150304 : residual_yy = residual_yy - (*_stress_model[i])[_qp](1, 1) - _mean_pressure;
487 : }
488 1888 : }
489 :
490 : void
491 1566064 : ComputeISoilStress::computeQpStress()
492 : {
493 : // Nothing to update during the first time step, return immediately
494 1566064 : if (_t_step == 0)
495 : return;
496 :
497 : _stress_new.zero();
498 1565120 : computeStress();
499 1565120 : _stress[_qp] = _rotation_increment[_qp] * _stress_new * _rotation_increment[_qp].transpose();
500 :
501 : // Compute dstress_dstrain
502 1565120 : if (_tangent_modulus == 0.0)
503 120848 : _tangent_modulus = _youngs[_pos][_youngs.size() - 1];
504 :
505 1565120 : Real lame_1 = _elasticity_tensor[_qp](0, 0, 1, 1);
506 1565120 : Real lame_2 = (_elasticity_tensor[_qp](0, 0, 0, 0) - lame_1) / 2;
507 1565120 : Real initial_youngs = (lame_2 * (3 * lame_1 + 2 * lame_2)) / (lame_1 + lame_2);
508 :
509 1565120 : _Jacobian_mult[_qp] =
510 1565120 : _elasticity_tensor[_qp] * _tangent_modulus / initial_youngs; // This is NOT the exact jacobian
511 : }
512 :
513 : void
514 1565120 : ComputeISoilStress::computeStress()
515 : {
516 1565120 : if (_t_step == 0)
517 0 : return;
518 :
519 : // Determine the current id for the soil. The variable which is a Real must be
520 : // converted to a unsigned int for lookup, so first
521 : // it is rounded to avoid Real values that are just below the desired value.
522 1565120 : _current_id = static_cast<unsigned int>(std::round(_soil_layer_variable[_qp]));
523 :
524 : // Get the position of the current id in the layer_ids array
525 1565120 : _pos = find(_layer_ids.begin(), _layer_ids.end(), _current_id) - _layer_ids.begin();
526 :
527 : _individual_stress_increment.zero();
528 :
529 : _deviatoric_trial_stress.zero();
530 :
531 1565120 : _tangent_modulus = 0.0;
532 :
533 : // current pressure calculation
534 1565120 : Real mean_stress = _stress_old[_qp].trace() / (-3.0);
535 1565120 : if (mean_stress < _p0)
536 : mean_stress = 0.0;
537 :
538 1565120 : if (_pressure_dependency)
539 : {
540 361152 : _stiffness_pressure_correction = pow((mean_stress - _p0) / _p_ref[_pos], _b_exp);
541 361152 : _strength_pressure_correction =
542 361152 : std::sqrt(_a0 + _a1 * (mean_stress - _p0) +
543 361152 : _a2 * (mean_stress - _p0) * (mean_stress - _p0)) /
544 361152 : std::sqrt(_a0 + _a1 * (_p_ref[_pos]) + _a2 * (_p_ref[_pos]) * (_p_ref[_pos]));
545 : }
546 :
547 1565120 : _mean_pressure = 0.0;
548 :
549 1565120 : Real lame_1 = _elasticity_tensor[_qp](0, 0, 1, 1);
550 1565120 : Real lame_2 = (_elasticity_tensor[_qp](0, 0, 0, 0) - lame_1) / 2;
551 1565120 : Real initial_youngs = (lame_2 * (3 * lame_1 + 2 * lame_2)) / (lame_1 + lame_2);
552 :
553 : // compute trial stress increment - note that _elasticity_tensor here
554 1565120 : _individual_stress_increment =
555 1565120 : _elasticity_tensor[_qp] * (_strain_increment[_qp]) / initial_youngs;
556 :
557 1565120 : _mean_pressure_tmp = _individual_stress_increment.trace() / 3.0 * _stiffness_pressure_correction;
558 :
559 1565120 : _deviatoric_trial_stress_tmp =
560 1565120 : _individual_stress_increment.deviatoric() * _stiffness_pressure_correction;
561 :
562 1565120 : _dev_trial_stress_squared = 0.0;
563 :
564 1565120 : _effective_trial_stress = 0.0;
565 :
566 1565120 : _yield_condition = 0.0;
567 :
568 48474752 : for (std::size_t i = 0; i < _base_models.size(); i++)
569 : {
570 :
571 : // calculate pressure for each element
572 46909632 : _mean_pressure += _mean_pressure_tmp * _youngs[_pos][i];
573 :
574 : // compute the deviatoric trial stress normalized by non pressure dependent
575 : // youngs modulus.
576 46909632 : _deviatoric_trial_stress =
577 46909632 : _deviatoric_trial_stress_tmp + (*_stress_model_old[i])[_qp] / (_youngs[_pos][i]);
578 :
579 : // compute the effective trial stress
580 46909632 : _dev_trial_stress_squared =
581 46909632 : _deviatoric_trial_stress.doubleContraction(_deviatoric_trial_stress);
582 46909632 : _effective_trial_stress = std::sqrt(3.0 / 2.0 * _dev_trial_stress_squared);
583 :
584 : // check yield condition and calculate plastic strain
585 46909632 : _yield_condition =
586 46909632 : _effective_trial_stress - _yield_stress[_pos][i] * _strength_pressure_correction;
587 :
588 46909632 : if (_yield_condition > 0.0)
589 : _deviatoric_trial_stress *=
590 18559303 : _yield_stress[_pos][i] * _strength_pressure_correction / _effective_trial_stress;
591 : else
592 28350329 : _tangent_modulus += _youngs[_pos][i];
593 :
594 46909632 : (*_stress_model[i])[_qp] = _youngs[_pos][i] * (_deviatoric_trial_stress);
595 :
596 46909632 : _stress_new += (*_stress_model[i])[_qp];
597 : }
598 1565120 : _stress_new(0, 0) += _mean_pressure - mean_stress;
599 1565120 : _stress_new(1, 1) += _mean_pressure - mean_stress;
600 1565120 : _stress_new(2, 2) += _mean_pressure - mean_stress;
601 :
602 1565120 : _tangent_modulus *= _stiffness_pressure_correction;
603 : }
|
