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 : // MASTODON includes
19 : #include "ISoilUtils.h"
20 : #include "MastodonUtils.h"
21 :
22 : // MOOSE includes
23 : #include "ColumnMajorMatrix.h"
24 : #include "Conversion.h"
25 : #include "FEProblem.h"
26 : #include "Factory.h"
27 : #include "MooseEnum.h"
28 : #include "DelimitedFileReader.h"
29 :
30 : void
31 136 : ISoilUtils::computeUserDefinedBackbone(std::vector<std::vector<Real>> & backbone_stress,
32 : std::vector<std::vector<Real>> & backbone_strain,
33 : const std::vector<unsigned int> & layer_ids,
34 : const std::vector<FileName> & backbone_curve_files,
35 : const std::string & name)
36 : {
37 136 : backbone_strain.resize(layer_ids.size());
38 136 : backbone_stress.resize(layer_ids.size());
39 273 : for (std::size_t i = 0; i < layer_ids.size(); i++)
40 : {
41 137 : MooseUtils::DelimitedFileReader backbone_curve(backbone_curve_files[i]);
42 : backbone_curve.setHeaderFlag(MooseUtils::DelimitedFileReader::HeaderFlag::OFF);
43 137 : backbone_curve.read();
44 137 : backbone_strain[i] = backbone_curve.getData("column_0");
45 274 : backbone_stress[i] = backbone_curve.getData("column_1");
46 137 : if (backbone_strain[i][0] == 0 && backbone_stress[i][0] == 0)
47 0 : mooseError("In " + name +
48 0 : ". Found 0, 0 in the stress-strain curve provided through the file " +
49 0 : backbone_curve_files[i] +
50 : ". Please make sure that the backbone curve does not start at 0, 0.");
51 137 : }
52 272 : if (!MastodonUtils::checkEqualSize(backbone_strain) ||
53 136 : !MastodonUtils::checkEqualSize(backbone_stress))
54 0 : mooseError("In " + name + ". Please make sure that the backbone curves of all layers in "
55 : "backbone_curve_files have the same number of points.");
56 136 : }
57 :
58 : void
59 73 : ISoilUtils::computeDarendeliBackbone(std::vector<std::vector<Real>> & backbone_stress,
60 : std::vector<std::vector<Real>> & backbone_strain,
61 : const std::vector<unsigned int> & layer_ids,
62 : const std::vector<Real> & initial_shear_modulus,
63 : const std::vector<Real> & over_consolidation_ratio,
64 : const std::vector<Real> & plasticity_index,
65 : const std::vector<Real> & p_ref,
66 : const unsigned int number_of_points,
67 : const std::string & name)
68 : {
69 73 : backbone_strain.resize(layer_ids.size());
70 73 : backbone_stress.resize(layer_ids.size());
71 : Real phi1 = 0.0352;
72 : Real phi2 = 0.001;
73 : Real phi3 = 0.3246;
74 : Real phi4 = 0.3483;
75 : Real phi5 = 0.9190;
76 147 : for (std::size_t i = 0; i < layer_ids.size(); i++)
77 : {
78 74 : if (plasticity_index[i] < 0)
79 0 : mooseError("Error in " + name + ". plasticity_index should be >= 0.");
80 74 : backbone_strain[i].resize(number_of_points);
81 74 : backbone_stress[i].resize(number_of_points);
82 74 : Real ref_strain = (phi1 + phi2 * plasticity_index[i] * pow(over_consolidation_ratio[i], phi3)) *
83 74 : pow(p_ref[i] * 0.00987, phi4);
84 1624 : for (std::size_t j = 0; j < number_of_points; j++)
85 : {
86 1550 : backbone_strain[i][j] = pow(10.0, -6.0 + 5.0 / (number_of_points - 1) * j);
87 1550 : backbone_stress[i][j] =
88 1550 : initial_shear_modulus[i] *
89 1550 : (1.0 / (1.0 + pow(100.0 * backbone_strain[i][j] / ref_strain, phi5))) *
90 : backbone_strain[i][j];
91 : }
92 : }
93 73 : }
94 :
95 : void
96 109 : ISoilUtils::computeGQHBackbone(std::vector<std::vector<Real>> & backbone_stress,
97 : std::vector<std::vector<Real>> & backbone_strain,
98 : const std::vector<unsigned int> & layer_ids,
99 : const std::vector<Real> & initial_shear_modulus,
100 : const unsigned int & number_of_points,
101 : const std::vector<Real> & theta_1,
102 : const std::vector<Real> & theta_2,
103 : const std::vector<Real> & theta_3,
104 : const std::vector<Real> & theta_4,
105 : const std::vector<Real> & theta_5,
106 : const std::vector<Real> & taumax)
107 : {
108 109 : backbone_strain.resize(layer_ids.size());
109 109 : backbone_stress.resize(layer_ids.size());
110 109 : std::vector<std::vector<Real>> modulus(layer_ids.size());
111 750 : for (std::size_t i = 0; i < layer_ids.size(); i++)
112 : {
113 641 : backbone_strain[i].resize(number_of_points);
114 641 : backbone_stress[i].resize(number_of_points);
115 641 : modulus[i].resize(number_of_points);
116 59701 : for (std::size_t j = 0; j < number_of_points; ++j)
117 : {
118 59060 : backbone_strain[i][j] = pow(10.0, -6.0 + 5.0 / (number_of_points - 1) * j);
119 : Real theta =
120 59060 : theta_1[i] +
121 59060 : theta_2[i] * theta_4[i] *
122 59060 : pow(backbone_strain[i][j] / (taumax[i] / initial_shear_modulus[i]), theta_5[i]) /
123 59060 : (pow(theta_3[i], theta_5[i]) +
124 59060 : pow(backbone_strain[i][j] / (taumax[i] / initial_shear_modulus[i]), theta_5[i]));
125 59060 : modulus[i][j] =
126 : (theta == 0)
127 59060 : ? 1.0 / (1.0 + backbone_strain[i][j] / (taumax[i] / initial_shear_modulus[i]))
128 59060 : : (1.0 / (backbone_strain[i][j] / (taumax[i] / initial_shear_modulus[i]))) *
129 59060 : (0.5 / theta *
130 59060 : (1.0 + backbone_strain[i][j] / (taumax[i] / initial_shear_modulus[i]) -
131 59060 : sqrt(pow(1.0 + backbone_strain[i][j] / (taumax[i] / initial_shear_modulus[i]),
132 : 2.0) -
133 59060 : 4.0 * theta * backbone_strain[i][j] /
134 : (taumax[i] / initial_shear_modulus[i]))));
135 59060 : backbone_stress[i][j] = modulus[i][j] * initial_shear_modulus[i] * backbone_strain[i][j];
136 : }
137 : }
138 109 : }
139 :
140 : void
141 37 : ISoilUtils::computeCoulombBackbone(std::vector<std::vector<Real>> & backbone_stress,
142 : std::vector<std::vector<Real>> & backbone_strain,
143 : const std::vector<unsigned int> & layer_ids,
144 : const std::vector<Real> & initial_shear_modulus,
145 : const std::vector<Real> & friction_coefficient,
146 : const std::vector<Real> & hardening_ratio,
147 : const std::vector<Real> & p_ref,
148 : const std::string & name)
149 : {
150 37 : backbone_strain.resize(layer_ids.size());
151 37 : backbone_stress.resize(layer_ids.size());
152 75 : for (std::size_t i = 0; i < layer_ids.size(); i++)
153 : {
154 38 : backbone_strain[i].resize(2);
155 38 : backbone_stress[i].resize(2);
156 38 : if (hardening_ratio[i] >= 1.0)
157 0 : mooseError("Error in " + name + ". hardening_ratio cannot be greater than or equal to 1.0.");
158 38 : backbone_strain[i][0] = p_ref[i] * friction_coefficient[i] / initial_shear_modulus[i];
159 38 : backbone_stress[i][0] = p_ref[i] * friction_coefficient[i];
160 38 : backbone_strain[i][1] = 100.0;
161 38 : backbone_stress[i][1] = backbone_stress[i][0] +
162 38 : (backbone_strain[i][1] - backbone_strain[i][0]) * hardening_ratio[i] *
163 : initial_shear_modulus[i];
164 : }
165 37 : }
166 :
167 : void
168 226 : ISoilUtils::computeSoilLayerProperties(
169 : std::vector<std::vector<Real>> & youngs,
170 : std::vector<std::vector<Real>> & yield_stress, // *** YIELD STRAIN, NOT YIELD STRESS***
171 : const std::vector<std::vector<Real>> & backbone_stress,
172 : const std::vector<std::vector<Real>> & backbone_strain,
173 : const std::vector<unsigned int> & layer_ids,
174 : const std::vector<Real> & poissons_ratio,
175 : const std::string & name)
176 : {
177 813 : for (std::size_t k = 0; k < layer_ids.size(); k++)
178 : {
179 587 : Real number = backbone_strain[k].size();
180 :
181 : // Calculating the Youngs modulus for each of the n curves for soil layer k
182 587 : ColumnMajorMatrix A(number, number);
183 587 : ColumnMajorMatrix InvA(number, number);
184 :
185 43406 : for (std::size_t i = 0; i < number; i++)
186 : {
187 4196454 : for (std::size_t j = 0; j < number; j++)
188 : {
189 4153635 : if (i <= j)
190 2098227 : A(i, j) = backbone_strain[k][i];
191 : else
192 2055408 : A(i, j) = backbone_strain[k][j];
193 : }
194 : }
195 :
196 : // create upper triangular matrix and modified stress
197 587 : youngs[k].resize(number);
198 587 : std::vector<Real> G0_component(number);
199 587 : yield_stress[k].resize(number);
200 587 : std::vector<Real> modified_backbone_stress(number);
201 : InvA = A;
202 42819 : for (std::size_t i = 1; i < number; i++)
203 : {
204 4153048 : for (std::size_t j = 0; j < number; j++)
205 : {
206 4110816 : InvA(i, j) = A(i, j) - A(i - 1, j);
207 4110816 : modified_backbone_stress[i] = backbone_stress[k][i] - backbone_stress[k][i - 1];
208 : }
209 : }
210 :
211 587 : modified_backbone_stress[0] = backbone_stress[k][0];
212 :
213 : // backward substitution
214 587 : G0_component[number - 1] = modified_backbone_stress[number - 1] / InvA(number - 1, number - 1);
215 587 : yield_stress[k][number - 1] = backbone_strain[k][number - 1];
216 587 : int i = number - 2;
217 42819 : while (i >= 0)
218 : {
219 42232 : G0_component[i] = 0.0;
220 2097640 : for (std::size_t j = i + 1; j < number; j++)
221 2055408 : G0_component[i] += InvA(i, j) * G0_component[j];
222 :
223 42232 : G0_component[i] = (modified_backbone_stress[i] - G0_component[i]) / InvA(i, i);
224 42232 : yield_stress[k][i] = backbone_strain[k][i];
225 42232 : i = i - 1;
226 : }
227 :
228 : // scaling
229 43406 : for (std::size_t i = 0; i < number; i++)
230 : {
231 42819 : youngs[k][i] = G0_component[i] * 2.0 * (1.0 + poissons_ratio[k]);
232 42819 : yield_stress[k][i] = yield_stress[k][i] * std::sqrt(3.0) / (2.0 * (1.0 + poissons_ratio[k]));
233 : }
234 587 : }
235 813 : for (std::size_t k = 0; k < layer_ids.size(); k++)
236 : {
237 43406 : for (std::size_t j = 0; j < backbone_strain[k].size(); j++)
238 : {
239 42819 : if (youngs[k][j] == 0)
240 0 : mooseError("In " + name + "Found two segments in the backbone curve with the same slope. ",
241 : "Please recheck input backbone curves.");
242 : }
243 : }
244 226 : }
|
