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 "ConcreteASREigenstrain.h"
16 :
17 : registerMooseObject("BlackBearApp", ConcreteASREigenstrain);
18 :
19 : InputParameters
20 1038 : ConcreteASREigenstrain::validParams()
21 : {
22 1038 : InputParameters params = ConcreteExpansionEigenstrainBase::validParams();
23 1038 : params.makeParamRequired<Real>("compressive_strength");
24 1038 : params.makeParamRequired<Real>("tensile_strength");
25 :
26 2076 : params.addRequiredCoupledVar("temperature", "Coupled temperature");
27 2076 : params.setDocUnit("temperature", "unit specified in 'temperature_unit'");
28 2076 : params.addRequiredCoupledVar("relative_humidity", "Coupled relative humidity");
29 2076 : params.setDocUnit("relative_humidity", "unitless");
30 :
31 3114 : params.addRangeCheckedParam<Real>(
32 : "rh_exponent",
33 2076 : 4.0,
34 : "rh_exponent >= 0.0",
35 : "Power to which relative humidity is raised in computation of ASR volumetric strain");
36 2076 : params.setDocUnit("rh_exponent", "unitless");
37 :
38 2076 : params.addRequiredRangeCheckedParam<Real>(
39 : "max_volumetric_expansion",
40 : "max_volumetric_expansion > 0.0",
41 : "Final ansymptotic ASR volumetric expansion strain when reaction is complete");
42 2076 : params.setDocUnit("max_volumetric_expansion", "unitless");
43 2076 : params.addRequiredRangeCheckedParam<Real>(
44 : "characteristic_time",
45 : "characteristic_time > 0.0",
46 : "Chracteristic ASR time at reference temprature. (tau_C(T_0))");
47 2076 : params.setDocUnit("characteristic_time", "day");
48 2076 : params.addRequiredParam<Real>("latency_time",
49 : "Latency ASR time at reference temprature (tau_L(T_0))");
50 2076 : params.setDocUnit("latency_time", "day");
51 3114 : params.addRangeCheckedParam<Real>("characteristic_activation_energy",
52 2076 : 5400.0,
53 : "characteristic_activation_energy > 0.0",
54 : "Activation energy associated with characteristic_time (U_C)");
55 2076 : params.setDocUnit("characteristic_activation_energy", "K");
56 3114 : params.addRangeCheckedParam<Real>("latency_activation_energy",
57 2076 : 9400.0,
58 : "latency_activation_energy > 0.0",
59 : "Activation energy associated with latency_time (U_L)");
60 2076 : params.setDocUnit("latency_activation_energy", "K");
61 2076 : params.addRequiredParam<Real>("reference_temperature",
62 : "Reference temperature for ASR reaction constants.");
63 2076 : params.setDocUnit("reference_temperature", "unit specified in 'temperature_unit'");
64 :
65 : // Note that Fahrenheit is not supported because that would require different parameters for the
66 : // times and activation energies
67 2076 : MooseEnum temperature_units("Celsius Kelvin");
68 2076 : params.addRequiredParam<MooseEnum>(
69 : "temperature_unit",
70 : temperature_units,
71 : "Model-specific temperature unit, which is used by the 'temperature' variable and "
72 : "'reference_temperature' parameter.");
73 :
74 3114 : params.addRangeCheckedParam<Real>(
75 : "stress_latency_factor",
76 2076 : 4.0 / 3.0,
77 : "stress_latency_factor >= 0.0",
78 : "Constant for ASR latency time retardation under hydrostatic compression (alpha)");
79 2076 : params.setDocUnit("stress_latency_factor", "unitless");
80 :
81 3114 : params.addRangeCheckedParam<Real>(
82 : "tensile_absorption_factor",
83 2076 : 0.5,
84 : "tensile_absorption_factor >= 0.0 & tensile_absorption_factor <= 1.0",
85 : "Fraction of tensile strength beyond which ASR gel is absorbed into tensile cracks "
86 : "(gamma_t)");
87 2076 : params.setDocUnit("tensile_absorption_factor", "unitless");
88 3114 : params.addRangeCheckedParam<Real>(
89 : "tensile_retention_factor",
90 2076 : 0.5,
91 : "tensile_retention_factor >= 0.0 & tensile_retention_factor <= 1.0",
92 : "Residual ASR retention factor under tension (gamma_r)");
93 2076 : params.setDocUnit("tensile_retention_factor", "unitless");
94 :
95 2076 : params.addParam<Real>("compressive_stress_exponent",
96 2076 : 0.5,
97 : "Exponent for ASR retention factor under compressive stress state (beta)");
98 2076 : params.setDocUnit("compressive_stress_exponent", "unitless");
99 :
100 2076 : params.addParam<bool>("ASR_dependent_tensile_strength",
101 2076 : false,
102 : "Set true to turn ASR reaction dependent tensile strength");
103 2076 : params.addRequiredRangeCheckedParam<Real>(
104 : "residual_tensile_strength_fraction",
105 : "residual_tensile_strength_fraction >= 0.0 & residual_tensile_strength_fraction <= 1.0",
106 : "Residual fraction of tensile strength at full ASR reaction");
107 2076 : params.setDocUnit("residual_tensile_strength_fraction", "unitless");
108 :
109 2076 : params.addParamNamesToGroup("stress_latency_factor tensile_absorption_factor "
110 : "tensile_retention_factor compressive_stress_exponent",
111 : "Advanced");
112 :
113 2076 : params.addParam<unsigned int>(
114 2076 : "max_its", 30, "Maximum number of iterations for material solution");
115 2076 : params.addParam<bool>(
116 2076 : "output_iteration_info", false, "Set true to output material iteration information");
117 2076 : params.addParam<bool>("output_iteration_info_on_error",
118 2076 : false,
119 : "Set true to output material iteration information when a step fails");
120 2076 : params.addParam<Real>(
121 2076 : "relative_tolerance", 1e-8, "Relative convergence tolerance for material iteration");
122 2076 : params.addParam<Real>(
123 2076 : "absolute_tolerance", 1e-20, "Absolute convergence tolerance for material iteration");
124 :
125 2076 : params.addParamNamesToGroup("max_its output_iteration_info output_iteration_info_on_error "
126 : "relative_tolerance absolute_tolerance",
127 : "Advanced");
128 :
129 1038 : params.addClassDescription(
130 : "Computes the volumetric expansion eigenstrain due to alkali-silica reaction.");
131 :
132 1038 : return params;
133 1038 : }
134 :
135 798 : ConcreteASREigenstrain::ConcreteASREigenstrain(const InputParameters & parameters)
136 : : ConcreteExpansionEigenstrainBase(parameters, "ASR"),
137 798 : _temperature(coupledValue("temperature")),
138 :
139 798 : _relative_humidity(coupledValue("relative_humidity")),
140 1596 : _rh_exponent(getParam<Real>("rh_exponent")),
141 :
142 1596 : _max_vol_expansion(getParam<Real>("max_volumetric_expansion")),
143 1596 : _tau_c_T0(getParam<Real>("characteristic_time")),
144 1596 : _tau_L_T0(getParam<Real>("latency_time")),
145 1596 : _Uc(getParam<Real>("characteristic_activation_energy")),
146 1596 : _UL(getParam<Real>("latency_activation_energy")),
147 1596 : _ref_temp(getParam<Real>("reference_temperature")),
148 1596 : _alpha(getParam<Real>("stress_latency_factor")),
149 :
150 1596 : _gamma_tensile(getParam<Real>("tensile_absorption_factor")),
151 1596 : _gamma_residual(getParam<Real>("tensile_retention_factor")),
152 1596 : _beta(getParam<Real>("compressive_stress_exponent")),
153 :
154 1596 : _ASR_dependent_tensile_strength(getParam<bool>("ASR_dependent_tensile_strength")),
155 1596 : _beta_f(getParam<Real>("residual_tensile_strength_fraction")),
156 :
157 798 : _max_its(parameters.get<unsigned int>("max_its")),
158 1596 : _output_iteration_info(getParam<bool>("output_iteration_info")),
159 1596 : _output_iteration_info_on_error(getParam<bool>("output_iteration_info_on_error")),
160 798 : _relative_tolerance(parameters.get<Real>("relative_tolerance")),
161 798 : _absolute_tolerance(parameters.get<Real>("absolute_tolerance")),
162 :
163 798 : _ASR_extent(declareProperty<Real>("ASR_extent")),
164 3192 : _ASR_extent_old(getMaterialPropertyOld<Real>("ASR_extent"))
165 : {
166 1596 : const MooseEnum temperature_unit = getParam<MooseEnum>("temperature_unit");
167 :
168 798 : if (temperature_unit == "Celsius")
169 756 : _temp_offset = 273.15;
170 42 : else if (temperature_unit == "Kelvin")
171 42 : _temp_offset = 0.0;
172 : else
173 0 : mooseError("Unsupported temperature unit");
174 :
175 : // Convert input value of reference temperature to Kelvin
176 798 : _ref_temp += _temp_offset;
177 798 : }
178 :
179 : void
180 5916 : ConcreteASREigenstrain::initQpStatefulProperties()
181 : {
182 5916 : _ASR_extent[_qp] = 0.0;
183 5916 : ConcreteExpansionEigenstrainBase::initQpStatefulProperties();
184 5916 : }
185 :
186 : Real
187 564220 : ConcreteASREigenstrain::computeQpVolumetricStrain()
188 : {
189 : // Use Newton iteration to determine ASR reaction extent at new time step
190 564220 : Real scalar = _ASR_extent_old[_qp];
191 : unsigned int it = 0;
192 : Real residual = 10;
193 : Real norm_residual = 10;
194 : Real first_norm_residual = 10;
195 :
196 564220 : std::stringstream iter_output;
197 2268636 : while (it < _max_its && norm_residual > _absolute_tolerance &&
198 2266168 : (norm_residual / first_norm_residual) > _relative_tolerance)
199 : {
200 1704416 : residual = computeResidual(_qp, scalar);
201 : norm_residual = std::abs(residual);
202 1704416 : if (it == 0)
203 : {
204 : first_norm_residual = norm_residual;
205 564220 : if (first_norm_residual == 0)
206 : first_norm_residual = 1;
207 : }
208 :
209 1704416 : scalar -= residual / computeDerivative(_qp, scalar);
210 :
211 1704416 : if (_output_iteration_info == true || _output_iteration_info_on_error == true)
212 : {
213 16 : iter_output << " it=" << it << " ASR_extent=" << scalar
214 8 : << " rel_res=" << norm_residual / first_norm_residual
215 8 : << " rel_tol=" << _relative_tolerance << " abs_res=" << norm_residual
216 8 : << " abs_tol=" << _absolute_tolerance << std::endl;
217 : }
218 1704416 : ++it;
219 : }
220 :
221 564220 : if (_output_iteration_info)
222 8 : _console << iter_output.str();
223 564220 : if (it == _max_its && norm_residual > _absolute_tolerance &&
224 4 : (norm_residual / first_norm_residual) > _relative_tolerance)
225 : {
226 4 : if (_output_iteration_info_on_error)
227 4 : Moose::err << iter_output.str();
228 4 : mooseError("Max material iteration hit during nonlinear constitutive model solve!");
229 : }
230 :
231 : // new ASR reaction extent
232 564216 : _ASR_extent[_qp] = scalar;
233 564216 : Real inc_ASR_extent = _ASR_extent[_qp] - _ASR_extent_old[_qp];
234 :
235 : // stress dependent total ASR volumetric accounting for ASR gel absorption due to tensile and
236 : // compressive cracking Eigen solve - Note the eigen values are ranked from minimum to maximum
237 564216 : const Real sig_k = _stress_eigenvalues[2];
238 :
239 : Real gel_absorption_tensile = 1.0;
240 : Real gel_absorption_compress = 1.0;
241 :
242 : // Optionally calculate effect of ASR reaction on the tensile strength
243 564216 : Real f_t = _f_tensile;
244 564216 : if (_ASR_dependent_tensile_strength)
245 322184 : f_t = _f_tensile * (1.0 - (1.0 - _beta_f) * _ASR_extent[_qp]);
246 :
247 564216 : if (sig_k > _gamma_tensile * f_t)
248 51040 : gel_absorption_tensile =
249 51040 : _gamma_residual + (1.0 - _gamma_residual) * (_gamma_tensile * f_t / sig_k);
250 :
251 564216 : Real sig_effective = _stress[_qp].trace() / (3.0 * -_f_compress);
252 :
253 564216 : if (sig_effective > 0.0)
254 : {
255 290948 : gel_absorption_compress =
256 290948 : 1.0 - std::exp(_beta) * sig_effective / (1.0 + (std::exp(_beta) - 1.0) * sig_effective);
257 290948 : if (gel_absorption_compress > 1.0)
258 : gel_absorption_compress = 1.0;
259 290948 : else if (gel_absorption_compress < 0.0)
260 : gel_absorption_compress = 0.0;
261 : }
262 :
263 564216 : const Real inc_ASR_volumetric_strain = gel_absorption_tensile * gel_absorption_compress *
264 564216 : std::pow(_relative_humidity[_qp], _rh_exponent) *
265 564216 : inc_ASR_extent * _max_vol_expansion;
266 :
267 564216 : return _volumetric_strain_old[_qp] + inc_ASR_volumetric_strain;
268 564216 : }
269 :
270 : Real
271 1704416 : ConcreteASREigenstrain::computeResidual(unsigned qp, Real scalar)
272 : {
273 : Real f;
274 1704416 : if (_expansion_type == ExpansionType::Isotropic)
275 : f = 1.0;
276 715550 : else if (_expansion_type == ExpansionType::Anisotropic)
277 : {
278 715550 : const Real I_sigma = _stress[_qp].trace();
279 715550 : if (I_sigma >= 0.0) // hydrostatic tension
280 : f = 1.0;
281 : else // hydrostatic compression: retarding ASR rection
282 : {
283 370104 : f = 1.0 + _alpha * I_sigma / (3.0 * -_f_compress);
284 : mooseAssert("f >= 1.0", "Wrong retardation for ASR latency time calculation!");
285 : }
286 : }
287 : else
288 0 : mooseError("Invalid expansion type");
289 :
290 : // Convert current temperature to Kelvin
291 1704416 : const Real T = _temperature[qp] + _temp_offset;
292 :
293 : // ASR characteristic and latency times (in days)
294 1704416 : Real tau_c = _tau_c_T0 * std::exp(_Uc * (1.0 / T - 1.0 / _ref_temp));
295 1704416 : Real tau_L = f * _tau_L_T0 * std::exp(_UL * (1.0 / T - 1.0 / _ref_temp));
296 :
297 1704416 : Real resid = scalar - _ASR_extent_old[qp] -
298 1704416 : (_dt * (1.0 - scalar) * (scalar + std::exp(-tau_L / tau_c))) /
299 1704416 : (86400.0 * tau_c * (1.0 + std::exp(-tau_L / tau_c)));
300 :
301 1704416 : return resid;
302 : }
303 :
304 : Real
305 1704416 : ConcreteASREigenstrain::computeDerivative(unsigned qp, Real scalar)
306 : {
307 : Real f;
308 1704416 : if (_expansion_type == ExpansionType::Isotropic)
309 : f = 1.0;
310 715550 : else if (_expansion_type == ExpansionType::Anisotropic)
311 : {
312 715550 : const Real I_sigma = _stress[_qp].trace();
313 715550 : if (I_sigma >= 0.0) // hydrostatic tension
314 : f = 1.0;
315 : else // hydrostatic compression: retarding ASR rection
316 : {
317 370104 : f = 1.0 + _alpha * I_sigma / (3.0 * -_f_compress);
318 : mooseAssert("f >= 1.0", "Wrong retardation for ASR latency time calculation!");
319 : }
320 : }
321 : else
322 0 : mooseError("Invalid expansion type");
323 :
324 : // Convert current temperature to Kelvin
325 1704416 : const Real T = _temperature[qp] + _temp_offset;
326 :
327 : // ASR characteristic and latency times (in days)
328 1704416 : Real tau_c = _tau_c_T0 * std::exp(_Uc * (1.0 / T - 1.0 / _ref_temp));
329 1704416 : Real tau_L = f * _tau_L_T0 * std::exp(_UL * (1.0 / T - 1.0 / _ref_temp));
330 :
331 1704416 : Real jac = 1.0 - _dt / 86400.0 / tau_c * (1.0 - scalar) / (1.0 + std::exp(-tau_L / tau_c)) -
332 1704416 : _dt / 86400.0 / tau_c * (-1) / (1.0 + std::exp(-tau_L / tau_c)) *
333 1704416 : (scalar + std::exp(-tau_L / tau_c));
334 1704416 : return jac;
335 : }
|
