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