Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "GrayLambertSurfaceRadiationBase.h"
11 : #include "MathUtils.h"
12 : #include "Function.h"
13 : #include "libmesh/quadrature.h"
14 :
15 : #include <cmath>
16 :
17 : InputParameters
18 264 : GrayLambertSurfaceRadiationBase::validParams()
19 : {
20 264 : InputParameters params = SideUserObject::validParams();
21 528 : params.addParam<Real>(
22 : "stefan_boltzmann_constant",
23 528 : 5.670367e-8,
24 : "The Stefan-Boltzmann constant. Default value is in units of [W / m^2 K^4].");
25 528 : params.addRequiredCoupledVar("temperature", "The coupled temperature variable.");
26 528 : params.addRequiredParam<std::vector<FunctionName>>("emissivity",
27 : "Emissivities for each boundary.");
28 528 : params.addParam<std::vector<BoundaryName>>(
29 : "fixed_temperature_boundary",
30 : {},
31 : "The list of boundary IDs from the mesh with fixed temperatures.");
32 528 : params.addParam<std::vector<FunctionName>>(
33 : "fixed_boundary_temperatures", {}, "The temperatures of the fixed boundary.");
34 528 : params.addParam<std::vector<BoundaryName>>(
35 : "adiabatic_boundary", {}, "The list of boundary IDs from the mesh that are adiabatic.");
36 :
37 264 : params.addClassDescription(
38 : "This object implements the exchange of heat by radiation between sidesets.");
39 264 : return params;
40 0 : }
41 :
42 142 : GrayLambertSurfaceRadiationBase::GrayLambertSurfaceRadiationBase(const InputParameters & parameters)
43 : : SideUserObject(parameters),
44 142 : _sigma_stefan_boltzmann(getParam<Real>("stefan_boltzmann_constant")),
45 142 : _n_sides(boundaryIDs().size()),
46 142 : _temperature(coupledValue("temperature")),
47 142 : _radiosity(_n_sides),
48 142 : _heat_flux_density(_n_sides),
49 142 : _side_temperature(_n_sides),
50 142 : _side_type(_n_sides),
51 142 : _areas(_n_sides),
52 142 : _beta(_n_sides),
53 284 : _surface_irradiation(_n_sides)
54 : {
55 : // get emissivity functions
56 142 : auto & eps_names = getParam<std::vector<FunctionName>>("emissivity");
57 142 : _emissivity.resize(eps_names.size());
58 1092 : for (unsigned int j = 0; j < _emissivity.size(); ++j)
59 950 : _emissivity[j] = &getFunctionByName(eps_names[j]);
60 :
61 : // set up the map from the side id to the local index & check
62 : // note that boundaryIDs is not in the right order anymore!
63 : {
64 284 : std::vector<BoundaryName> boundary_names = getParam<std::vector<BoundaryName>>("boundary");
65 :
66 1098 : for (unsigned int j = 0; j < boundary_names.size(); ++j)
67 : {
68 956 : if (boundary_names[j] == "ANY_BOUNDARY_ID")
69 0 : paramError("boundary", "boundary must be explicitly provided.");
70 :
71 956 : _side_id_index[_mesh.getBoundaryID(boundary_names[j])] = j;
72 956 : _side_type[j] = VARIABLE_TEMPERATURE;
73 : }
74 :
75 : // consistency check on emissivity, must be as many entries as boundary
76 142 : if (boundary_names.size() != _emissivity.size())
77 2 : paramError("emissivity", "The number of entries must match the number of boundary entries.");
78 140 : }
79 :
80 : // get the fixed boundaries of the system if any are provided
81 280 : if (isParamValid("fixed_temperature_boundary"))
82 : {
83 : // if fixed_temperature_boundary is valid, then fixed_side_temperatures must be
84 : // valid, too
85 280 : if (!isParamValid("fixed_boundary_temperatures"))
86 0 : paramError("fixed_boundary_temperatures",
87 : "fixed_temperature_boundary is provided, but fixed_boundary_temperatures is not.");
88 :
89 420 : auto fst_fn = getParam<std::vector<FunctionName>>("fixed_boundary_temperatures");
90 296 : for (auto & fn : fst_fn)
91 156 : _fixed_side_temperature.push_back(&getFunctionByName(fn));
92 :
93 : // get the fixed boundaries and temperatures
94 : std::vector<BoundaryName> boundary_names =
95 420 : getParam<std::vector<BoundaryName>>("fixed_temperature_boundary");
96 :
97 140 : if (boundary_names.size() != _fixed_side_temperature.size())
98 2 : paramError(
99 : "fixed_boundary_temperatures",
100 : "fixed_boundary_temperatures and fixed_temperature_boundary must have the same length.");
101 :
102 : unsigned int index = 0;
103 290 : for (auto & name : boundary_names)
104 : {
105 152 : _fixed_side_id_index[_mesh.getBoundaryID(name)] = index;
106 152 : index++;
107 : }
108 :
109 : // check that fixed side ids is a subset of boundary ids
110 : // and update _side_type info
111 288 : for (auto & p : _fixed_side_id_index)
112 : {
113 152 : if (_side_id_index.find(p.first) == _side_id_index.end())
114 2 : paramError("fixed_temperature_boundary",
115 : "fixed_temperature_boundary must be a subset of boundary.");
116 150 : _side_type[_side_id_index.find(p.first)->second] = FIXED_TEMPERATURE;
117 : }
118 136 : }
119 :
120 : // get the fixed boundaries of the system if any are provided
121 272 : if (isParamValid("adiabatic_boundary"))
122 : {
123 : // get the adiabatic boundaries and temperatures
124 : std::vector<BoundaryName> boundary_names =
125 408 : getParam<std::vector<BoundaryName>>("adiabatic_boundary");
126 :
127 480 : for (auto & name : boundary_names)
128 344 : _adiabatic_side_ids.insert(_mesh.getBoundaryID(name));
129 :
130 : // check that adiabatic side ids is a subset of boundary ids
131 : // and update _side_type info
132 476 : for (auto & id : _adiabatic_side_ids)
133 : {
134 342 : if (_side_id_index.find(id) == _side_id_index.end())
135 2 : paramError("adiabatic_boundary", "adiabatic_boundary must be a subset of boundary.");
136 340 : _side_type[_side_id_index.find(id)->second] = ADIABATIC;
137 : }
138 :
139 : // make sure that adiabatic boundaries are not already fixed boundaries
140 474 : for (auto & id : _adiabatic_side_ids)
141 340 : if (_fixed_side_id_index.find(id) != _fixed_side_id_index.end())
142 0 : paramError("adiabatic_boundary", "Isothermal boundary cannot also be adiabatic boundary.");
143 134 : }
144 134 : }
145 :
146 : void
147 133068 : GrayLambertSurfaceRadiationBase::execute()
148 : {
149 : mooseAssert(_side_id_index.find(_current_boundary_id) != _side_id_index.end(),
150 : "Current boundary id not in _side_id_index.");
151 133068 : unsigned int index = _side_id_index.find(_current_boundary_id)->second;
152 :
153 465636 : for (unsigned int qp = 0; qp < _qrule->n_points(); qp++)
154 : {
155 332568 : _areas[index] += _JxW[qp] * _coord[qp];
156 :
157 : Real temp = 0;
158 332568 : if (_side_type[index] == ADIABATIC)
159 129632 : continue;
160 202936 : else if (_side_type[index] == VARIABLE_TEMPERATURE)
161 186324 : temp = _temperature[qp];
162 16612 : else if (_side_type[index] == FIXED_TEMPERATURE)
163 : {
164 16612 : unsigned int iso_index = _fixed_side_id_index.find(_current_boundary_id)->second;
165 16612 : temp = _fixed_side_temperature[iso_index]->value(_t, _q_point[qp]);
166 : }
167 :
168 202936 : _beta[index] += _JxW[qp] * _coord[qp] * _sigma_stefan_boltzmann *
169 202936 : _emissivity[index]->value(temp, Point()) * MathUtils::pow(temp, 4);
170 202936 : _side_temperature[index] += _JxW[qp] * _coord[qp] * temp;
171 : }
172 133068 : }
173 :
174 : void
175 9070 : GrayLambertSurfaceRadiationBase::initialize()
176 : {
177 : // view factors are obtained here to make sure that another object had
178 : // time to compute them on exec initial
179 9070 : _view_factors = setViewFactors();
180 :
181 : // initialize areas, beta, side temps
182 71914 : for (unsigned int j = 0; j < _n_sides; ++j)
183 : {
184 62844 : _areas[j] = 0;
185 62844 : _beta[j] = 0;
186 62844 : _side_temperature[j] = 0;
187 : }
188 9070 : }
189 :
190 : void
191 7923 : GrayLambertSurfaceRadiationBase::finalize()
192 : {
193 : // need to do some parallel communiction here
194 7923 : gatherSum(_areas);
195 7923 : gatherSum(_beta);
196 7923 : gatherSum(_side_temperature);
197 :
198 : // first compute averages from the totals
199 62556 : for (unsigned int j = 0; j < _n_sides; ++j)
200 : {
201 54633 : _beta[j] /= _areas[j];
202 54633 : _side_temperature[j] /= _areas[j];
203 : }
204 :
205 : // matrix and rhs vector for the view factor calculation
206 7923 : DenseMatrix<Real> matrix(_n_sides, _n_sides);
207 7923 : DenseVector<Real> rhs(_n_sides);
208 7923 : DenseVector<Real> radiosity(_n_sides);
209 62556 : for (unsigned int i = 0; i < _n_sides; ++i)
210 : {
211 54633 : rhs(i) = _beta[i];
212 54633 : matrix(i, i) = 1;
213 471684 : for (unsigned int j = 0; j < _n_sides; ++j)
214 : {
215 417051 : if (_side_type[i] == ADIABATIC)
216 132724 : matrix(i, j) -= _view_factors[i][j];
217 : else
218 284327 : matrix(i, j) -=
219 284327 : (1 - _emissivity[i]->value(_side_temperature[i], Point())) * _view_factors[i][j];
220 : }
221 : }
222 :
223 : // compute the radiosityes
224 7923 : matrix.lu_solve(rhs, radiosity);
225 :
226 : // store the radiosity, temperatures and heat flux density for each surface
227 62556 : for (unsigned int i = 0; i < _n_sides; ++i)
228 : {
229 54633 : _radiosity[i] = radiosity(i);
230 :
231 : // _heat_flux_density is obtained from a somewhat cumbersome relation
232 : // but it has the advantage that we do not divide by 1 - emissivity
233 : // which blows up for black bodies
234 54633 : _heat_flux_density[i] = radiosity(i);
235 471684 : for (unsigned int j = 0; j < _n_sides; ++j)
236 417051 : _heat_flux_density[i] -= _view_factors[i][j] * radiosity(j);
237 :
238 54633 : if (_side_type[i] == ADIABATIC)
239 : {
240 18592 : Real e = _emissivity[i]->value(_side_temperature[i], Point());
241 18592 : _side_temperature[i] = std::pow(
242 18592 : (radiosity(i) + (1 - e) / e * _heat_flux_density[i]) / _sigma_stefan_boltzmann, 0.25);
243 : }
244 :
245 : // compute the surface irradiation into i from the radiosities
246 54633 : _surface_irradiation[i] = 0;
247 471684 : for (unsigned int j = 0; j < _n_sides; ++j)
248 417051 : _surface_irradiation[i] += _view_factors[i][j] * radiosity(j);
249 : }
250 7923 : }
251 :
252 : void
253 1147 : GrayLambertSurfaceRadiationBase::threadJoin(const UserObject & y)
254 : {
255 : const auto & pps = static_cast<const GrayLambertSurfaceRadiationBase &>(y);
256 :
257 9358 : for (unsigned int j = 0; j < _n_sides; ++j)
258 : {
259 8211 : _areas[j] += pps._areas[j];
260 8211 : _side_temperature[j] += pps._side_temperature[j];
261 8211 : _beta[j] += pps._beta[j];
262 : }
263 1147 : }
264 :
265 : std::set<BoundaryID>
266 9 : GrayLambertSurfaceRadiationBase::getSurfaceIDs() const
267 : {
268 : std::set<BoundaryID> surface_ids;
269 45 : for (auto & p : _side_id_index)
270 36 : surface_ids.insert(p.first);
271 9 : return surface_ids;
272 : }
273 :
274 : unsigned int
275 1877658 : GrayLambertSurfaceRadiationBase::getSideIDIndex(BoundaryID id) const
276 : {
277 : const auto it = _side_id_index.find(id);
278 1877658 : if (it == _side_id_index.end())
279 0 : mooseError("Sideset ID ", id, " (", _mesh.getBoundaryName(id), ") not found in map.");
280 1877658 : return it->second;
281 : }
282 :
283 : Real
284 633352 : GrayLambertSurfaceRadiationBase::getSurfaceIrradiation(BoundaryID id) const
285 : {
286 633352 : return _surface_irradiation[getSideIDIndex(id)];
287 : }
288 :
289 : Real
290 243403 : GrayLambertSurfaceRadiationBase::getSurfaceHeatFluxDensity(BoundaryID id) const
291 : {
292 243403 : return _heat_flux_density[getSideIDIndex(id)];
293 : }
294 :
295 : Real
296 82 : GrayLambertSurfaceRadiationBase::getSurfaceTemperature(BoundaryID id) const
297 : {
298 82 : return _side_temperature[getSideIDIndex(id)];
299 : }
300 :
301 : Real
302 73 : GrayLambertSurfaceRadiationBase::getSurfaceRadiosity(BoundaryID id) const
303 : {
304 73 : return _radiosity[getSideIDIndex(id)];
305 : }
306 :
307 : Real
308 1000748 : GrayLambertSurfaceRadiationBase::getSurfaceEmissivity(BoundaryID id) const
309 : {
310 1000748 : const auto ind = getSideIDIndex(id);
311 1000748 : return _emissivity[ind]->value(_side_temperature[ind], Point());
312 : }
313 :
314 : Real
315 0 : GrayLambertSurfaceRadiationBase::getViewFactor(BoundaryID from_id, BoundaryID to_id) const
316 : {
317 0 : return _view_factors[getSideIDIndex(from_id)][getSideIDIndex(to_id)];
318 : }
|
