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