Line data Source code
1 : /****************************************************************/
2 : /* DO NOT MODIFY THIS HEADER */
3 : /* BlackBear */
4 : /* */
5 : /* (c) 2017 Battelle Energy Alliance, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /* */
8 : /* Prepared by Battelle Energy Alliance, LLC */
9 : /* Under Contract No. DE-AC07-05ID14517 */
10 : /* With the U. S. Department of Energy */
11 : /* */
12 : /* See COPYRIGHT for full restrictions */
13 : /****************************************************************/
14 :
15 : #include "ConcreteExpansionEigenstrainBase.h"
16 : #include "RankTwoTensor.h"
17 :
18 : InputParameters
19 1996 : ConcreteExpansionEigenstrainBase::validParams()
20 : {
21 1996 : InputParameters params = ComputeEigenstrainBase::validParams();
22 3992 : MooseEnum expansion_type("Isotropic Anisotropic");
23 3992 : params.addParam<MooseEnum>(
24 : "expansion_type", expansion_type, "Type of expansion resulting from volumetric strain");
25 3992 : params.addRangeCheckedParam<Real>(
26 : "compressive_strength", "compressive_strength > 0", "Compressive strength of concrete");
27 3992 : params.addRangeCheckedParam<Real>(
28 : "expansion_stress_limit",
29 : "expansion_stress_limit > 0",
30 : "Upper bound compressive stress beyond which no expansion occurs");
31 3992 : params.addRangeCheckedParam<Real>(
32 : "tensile_strength", "tensile_strength > 0", "Tensile strength of concrete");
33 1996 : return params;
34 1996 : }
35 :
36 1527 : ConcreteExpansionEigenstrainBase::ConcreteExpansionEigenstrainBase(
37 1527 : const InputParameters & parameters, const std::string volumetric_expansion_name)
38 : : ComputeEigenstrainBase(parameters),
39 1527 : _expansion_type(getParam<MooseEnum>("expansion_type").getEnum<ExpansionType>()),
40 5238 : _f_compress(isParamValid("compressive_strength") ? getParam<Real>("compressive_strength")
41 : : 0.0),
42 4986 : _sigma_u(isParamValid("expansion_stress_limit") ? getParam<Real>("expansion_stress_limit")
43 : : 0.0),
44 5238 : _f_tensile(isParamValid("tensile_strength") ? getParam<Real>("tensile_strength") : 0.0),
45 1527 : _eigenstrain_old(getMaterialPropertyOld<RankTwoTensor>(_eigenstrain_name)),
46 1527 : _volumetric_strain(declareProperty<Real>(volumetric_expansion_name + "_volumetric_strain")),
47 1527 : _volumetric_strain_old(
48 1527 : getMaterialPropertyOld<Real>(volumetric_expansion_name + "_volumetric_strain")),
49 6108 : _stress(getMaterialPropertyOld<RankTwoTensor>("stress"))
50 : {
51 1527 : if (_expansion_type == ExpansionType::Anisotropic)
52 : {
53 1854 : if (!parameters.isParamSetByUser("compressive_strength"))
54 3 : paramError("compressive_strength",
55 : "compressive_strength is required for expansion_type = Anisotropic");
56 1848 : if (!parameters.isParamSetByUser("expansion_stress_limit"))
57 3 : paramError("expansion_stress_limit",
58 : "expansion_stress_limit is required for expansion_type = Anisotropic");
59 1842 : if (!parameters.isParamSetByUser("tensile_strength"))
60 3 : paramError("tensile_strength",
61 : "tensile_strength is required for expansion_type = Anisotropic");
62 : }
63 :
64 : // Initialize triaxial weight table
65 : const Real third = 1.0 / 3.0;
66 1518 : _triaxial_weights = {{{{third, third, 0.0, 0.0}}, // 0 -> 13
67 : {{third, third, 0.0, 0.0}}, // 1 -> 14
68 : {{0.5, 0.5, 0.0, 0.0}}, // 2 -> 15
69 : {{0.5, 0.5, 0.0, 0.0}}, // 3 -> 16
70 : {{third, third, 0.0, 0.0}}, // 4 -> 12
71 : {{third, third, 0.0, 0.0}}, // 5 -> 1
72 : {{0.5, 0.5, 0.0, 0.0}}, // 6 -> 2
73 : {{0.5, 0.5, 0.0, 0.0}}, // 7 -> 5
74 : {{0.5, 0.5, 0.0, 0.0}}, // 8 -> 11
75 : {{0.5, 0.5, 0.0, 0.0}}, // 9 -> 4
76 : {{1.0, 1.0, third, 0.0}}, // 10 -> 3
77 : {{1.0, 1.0, 0.5, 0.0}}, // 11 -> 6
78 : {{0.5, 0.5, 0.0, 0.0}}, // 12 -> 10
79 : {{0.5, 0.5, 0.0, 0.0}}, // 13 -> 9
80 : {{1.0, 1.0, 0.5, 0.0}}, // 14 -> 8
81 : {{1.0, 1.0, 1.0, third}}}}; // 15 -> 7
82 1518 : }
83 :
84 : void
85 13136 : ConcreteExpansionEigenstrainBase::initQpStatefulProperties()
86 : {
87 13136 : _volumetric_strain[_qp] = 0.0;
88 13136 : }
89 :
90 : void
91 967034 : ConcreteExpansionEigenstrainBase::computeQpEigenstrain()
92 : {
93 967034 : if (_expansion_type == ExpansionType::Anisotropic || needStressEigenvalues())
94 816954 : _stress[_qp].symmetricEigenvaluesEigenvectors(_stress_eigenvalues, _stress_eigenvectors);
95 :
96 967034 : _volumetric_strain[_qp] = computeQpVolumetricStrain();
97 :
98 967032 : if (_expansion_type == ExpansionType::Isotropic)
99 : {
100 514502 : const Real eigenstrain_comp = computeVolumetricStrainComponent(_volumetric_strain[_qp]);
101 514502 : _eigenstrain[_qp].zero();
102 514502 : _eigenstrain[_qp].addIa(eigenstrain_comp);
103 : }
104 452530 : else if (_expansion_type == ExpansionType::Anisotropic)
105 : {
106 452530 : const Real inc_volumetric_strain = _volumetric_strain[_qp] - _volumetric_strain_old[_qp];
107 :
108 452530 : RankTwoTensor inc_weighted_strain;
109 :
110 452530 : Real sig_l = _stress_eigenvalues[0];
111 452530 : Real sig_m = _stress_eigenvalues[1];
112 452530 : Real sig_k = _stress_eigenvalues[2];
113 :
114 : // Weights
115 : // W_k
116 452530 : Real W_3 = weight(sig_l, sig_m, sig_k);
117 : // W_m
118 452530 : Real W_2 = weight(sig_k, sig_l, sig_m);
119 : // W_l
120 452530 : Real W_1 = weight(sig_m, sig_k, sig_l);
121 :
122 452530 : inc_weighted_strain(0, 0) = inc_volumetric_strain * W_1;
123 452530 : inc_weighted_strain(1, 1) = inc_volumetric_strain * W_2;
124 452530 : inc_weighted_strain(2, 2) = inc_volumetric_strain * W_3;
125 :
126 : // Rotate from principal coordinate system to Cartesian coordinate system
127 452530 : inc_weighted_strain.rotate(_stress_eigenvectors);
128 :
129 452530 : _eigenstrain[_qp] = _eigenstrain_old[_qp] + inc_weighted_strain;
130 : }
131 : else
132 0 : mooseError("Invalid expansion type");
133 967032 : }
134 :
135 : Real
136 1357590 : ConcreteExpansionEigenstrainBase::weight(Real sig_l, Real sig_m, Real sig_k)
137 : {
138 : // Inputs
139 1357590 : const Real a1 = _f_tensile, a2 = -_sigma_u, a3 = -_f_compress + _sigma_u;
140 : const Real b1 = _f_tensile, b2 = -_sigma_u, b3 = -_f_compress + _sigma_u;
141 :
142 : // Lower and upper bound for l
143 1357590 : const unsigned int pbound_l = findNeighborIndex(sig_l);
144 :
145 : // Lower and upper bound for m
146 1357590 : const unsigned int pbound_m = findNeighborIndex(sig_m);
147 :
148 : // Neighboring nodes of k
149 1357590 : const unsigned int pbound_k = findNeighborIndex(sig_k);
150 :
151 : // Neighboring nodes of l, m, k
152 1357590 : const unsigned int N1 = (pbound_m)*4 + pbound_l;
153 1357590 : const unsigned int N2 = (pbound_m)*4 + (pbound_l + 1);
154 1357590 : const unsigned int N3 = (pbound_m + 1) * 4 + pbound_l;
155 1357590 : const unsigned int N4 = (pbound_m + 1) * 4 + (pbound_l + 1);
156 :
157 : const unsigned int N5 = pbound_k;
158 1357590 : const unsigned int N6 = pbound_k + 1;
159 :
160 : // Calculate a, b, sig_l, sig_m
161 1357590 : const Real a = computeAB(a1, a2, a3, pbound_l);
162 1357590 : const Real b = computeAB(b1, b2, b3, pbound_m);
163 1357590 : sig_l = computeSigma(sig_l, pbound_l);
164 1357590 : sig_m = computeSigma(sig_m, pbound_m);
165 :
166 : // Calculate weights
167 1357590 : return computeW(N1, N2, N3, N4, N5, N6, a, b, sig_l, sig_m, sig_k);
168 : }
169 :
170 : unsigned int
171 4072770 : ConcreteExpansionEigenstrainBase::findNeighborIndex(Real sig)
172 : {
173 4072770 : if (sig <= -_sigma_u)
174 : return 2;
175 3345162 : else if (sig > -_sigma_u && sig <= 0)
176 : return 1;
177 : else
178 666057 : return 0;
179 : }
180 :
181 : Real
182 2715180 : ConcreteExpansionEigenstrainBase::computeAB(const Real ab1,
183 : const Real ab2,
184 : const Real ab3,
185 : const unsigned int pbound)
186 : {
187 : mooseAssert("pbound <= 2", "pbound outside allowed range");
188 2715180 : if (pbound == 0)
189 : return ab1;
190 2271142 : else if (pbound == 1)
191 : return ab2;
192 : else
193 485072 : return ab3;
194 : }
195 :
196 : Real
197 2715180 : ConcreteExpansionEigenstrainBase::computeSigma(const Real sig, const unsigned int pbound)
198 : {
199 2715180 : if (pbound == 2)
200 : {
201 485072 : if (sig < -_f_compress) // stress conditions beyond _f_compress
202 30240 : return -_f_compress + _sigma_u;
203 : else
204 454832 : return sig + _sigma_u;
205 : }
206 : else
207 : {
208 2230108 : if (sig > _f_tensile)
209 : return _f_tensile; // stress conditions beyond _f_tensile
210 : else
211 2081916 : return sig;
212 : }
213 : }
214 :
215 : Real
216 1357590 : ConcreteExpansionEigenstrainBase::computeW(const unsigned int N1,
217 : const unsigned int N2,
218 : const unsigned int N3,
219 : const unsigned int N4,
220 : const unsigned int N5,
221 : const unsigned int N6,
222 : const Real a,
223 : const Real b,
224 : const Real sig_l,
225 : const Real sig_m,
226 : Real sig_k)
227 : {
228 1357590 : if (sig_k < -_f_compress)
229 : sig_k = -_f_compress;
230 1342470 : else if (sig_k > _f_tensile)
231 : sig_k = _f_tensile;
232 :
233 : // Calculate weights
234 : // N1
235 1357590 : const Real W1 = computeWi(N1, N5, N6, sig_k);
236 : // N2
237 1357590 : const Real W2 = computeWi(N2, N5, N6, sig_k);
238 : // N3
239 1357590 : const Real W3 = computeWi(N3, N5, N6, sig_k);
240 : // N4
241 1357590 : const Real W4 = computeWi(N4, N5, N6, sig_k);
242 :
243 : // local node numebers
244 : // 3-----4
245 : // | |
246 : // | |
247 : // 1-----2
248 1357590 : const Real WN1 = (1.0 / (a * b)) * (a - sig_l) * (b - sig_m);
249 1357590 : const Real WN2 = (1.0 / (a * b)) * sig_l * (b - sig_m);
250 1357590 : const Real WN3 = (1.0 / (a * b)) * (a - sig_l) * sig_m;
251 1357590 : const Real WN4 = (1.0 / (a * b)) * sig_l * sig_m;
252 :
253 1357590 : return W1 * WN1 + W2 * WN2 + W3 * WN3 + W4 * WN4;
254 : }
255 :
256 : Real
257 5430360 : ConcreteExpansionEigenstrainBase::computeWi(const unsigned int N,
258 : const unsigned int N5,
259 : const unsigned int N6,
260 : const Real sig_k)
261 : {
262 5430360 : const Real farr[4] = {_f_tensile, 0, -_sigma_u, -_f_compress};
263 5430360 : const Real W_1 = _triaxial_weights[N][N5];
264 5430360 : const Real W_2 = _triaxial_weights[N][N6];
265 5430360 : return W_1 + (W_2 - W_1) * (sig_k - farr[N5]) / (farr[N6] - farr[N5]);
266 : }
|
