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 : // Moose includes
11 : #include "CHTHandler.h"
12 : #include "LinearFVFluxKernel.h"
13 : #include "LinearFVDiffusion.h"
14 : #include "LinearFVAnisotropicDiffusion.h"
15 : #include "LinearFVBoundaryCondition.h"
16 : #include "LinearFVCHTBCInterface.h"
17 : #include "FEProblemBase.h"
18 :
19 : namespace NS
20 : {
21 : namespace FV
22 : {
23 :
24 : InputParameters
25 1784 : CHTHandler::validParams()
26 : {
27 1784 : auto params = emptyInputParameters();
28 3568 : params.addParam<std::vector<BoundaryName>>(
29 : "cht_interfaces",
30 : {},
31 : "The interfaces where we would like to add conjugate heat transfer handling.");
32 :
33 5352 : params.addRangeCheckedParam<unsigned int>(
34 : "max_cht_fpi",
35 3568 : 1,
36 : "max_cht_fpi >= 1",
37 : "Number of maximum fixed point iterations (FPI). Currently only applied to"
38 : " conjugate heat transfer simulations. The default value of 1 essentially keeps"
39 : " the FPI feature turned off. CHT iteration ends after this number of iteration even if the "
40 : "tolerance is not met.");
41 :
42 5352 : params.addRangeCheckedParam<Real>(
43 : "cht_heat_flux_tolerance",
44 3568 : 1e-5,
45 : "cht_heat_flux_tolerance > 0 & cht_heat_flux_tolerance <= 1.0",
46 : "The relative tolerance for terminating conjugate heat transfer iteration before the maximum "
47 : "number of CHT iterations. Relative tolerance is ignore if the maximum number of CHT "
48 : "iterations is reached.");
49 :
50 3568 : params.addParam<std::vector<Real>>(
51 : "cht_fluid_temperature_relaxation",
52 : {},
53 : "The relaxation factors for the boundary temperature when being updated on the fluid side.");
54 3568 : params.addParam<std::vector<Real>>(
55 : "cht_solid_temperature_relaxation",
56 : {},
57 : "The relaxation factors for the boundary temperature when being updated on the solid side.");
58 3568 : params.addParam<std::vector<Real>>(
59 : "cht_fluid_flux_relaxation",
60 : {},
61 : "The relaxation factors for the boundary flux when being updated on the fluid side.");
62 3568 : params.addParam<std::vector<Real>>(
63 : "cht_solid_flux_relaxation",
64 : {},
65 : "The relaxation factors for the boundary flux when being updated on the solid side.");
66 :
67 3568 : params.addParamNamesToGroup(
68 : "cht_interfaces max_cht_fpi cht_heat_flux_tolerance cht_fluid_temperature_relaxation "
69 : "cht_solid_temperature_relaxation cht_fluid_flux_relaxation "
70 : "cht_solid_flux_relaxation",
71 : "Conjugate Heat Transfer");
72 :
73 1784 : return params;
74 0 : }
75 :
76 892 : CHTHandler::CHTHandler(const InputParameters & params)
77 : : MooseObject(params),
78 1784 : _problem(*getCheckedPointerParam<FEProblemBase *>(
79 : "_fe_problem_base", "This might happen if you don't have a mesh")),
80 892 : _mesh(_problem.mesh()),
81 1784 : _cht_boundary_names(getParam<std::vector<BoundaryName>>("cht_interfaces")),
82 892 : _cht_boundary_ids(_mesh.getBoundaryIDs(_cht_boundary_names)),
83 1784 : _max_cht_fpi(getParam<unsigned int>("max_cht_fpi")),
84 2676 : _cht_heat_flux_tolerance(getParam<Real>("cht_heat_flux_tolerance"))
85 : {
86 2676 : if (isParamSetByUser("cht_interfaces") && !_cht_boundary_names.size())
87 0 : paramError("cht_interfaces", "You must declare at least one interface!");
88 892 : }
89 :
90 : void
91 36 : CHTHandler::linkEnergySystems(SystemBase * solid_energy_system, SystemBase * fluid_energy_system)
92 : {
93 36 : _energy_system = fluid_energy_system;
94 36 : _solid_energy_system = solid_energy_system;
95 :
96 36 : if (!_energy_system || !_solid_energy_system)
97 0 : paramError("cht_interfaces",
98 : "You selected to do conjugate heat transfer treatment, but it needs two energy "
99 : "systems: a solid and a fluid. One of these systems is missing.");
100 36 : }
101 :
102 : void
103 36 : CHTHandler::deduceCHTBoundaryCoupling()
104 : {
105 36 : if (_solid_energy_system->nVariables() != 1)
106 0 : mooseError("We should have only one variable in the solid energy system: ",
107 0 : _energy_system->name(),
108 : "! Right now we have: ",
109 0 : Moose::stringify(_solid_energy_system->getVariableNames()));
110 36 : if (_energy_system->nVariables() != 1)
111 0 : mooseError("We should have only one variable in the fluid energy system: ",
112 0 : _energy_system->name(),
113 : "! Right now we have: ",
114 0 : Moose::stringify(_energy_system->getVariableNames()));
115 108 : const std::vector<std::string> solid_fluid({"solid", "fluid"});
116 :
117 : // We do some setup at the beginning to make sure the container sizes are good
118 : _cht_system_numbers =
119 36 : std::vector<unsigned int>({_solid_energy_system->number(), _energy_system->number()});
120 36 : _cht_conduction_kernels = std::vector<LinearFVFluxKernel *>({nullptr, nullptr});
121 36 : _cht_boundary_conditions.clear();
122 36 : _cht_boundary_conditions.resize(_cht_boundary_names.size(), {nullptr, nullptr});
123 :
124 : const auto flux_relaxation_param_names =
125 108 : std::vector<std::string>({"cht_solid_flux_relaxation", "cht_fluid_flux_relaxation"});
126 : const auto temperature_relaxation_param_names = std::vector<std::string>(
127 108 : {"cht_solid_temperature_relaxation", "cht_fluid_temperature_relaxation"});
128 36 : _cht_flux_relaxation_factor.clear();
129 36 : _cht_flux_relaxation_factor.resize(2, std::vector<Real>(_cht_boundary_names.size(), 1.0));
130 36 : _cht_temperature_relaxation_factor.clear();
131 36 : _cht_temperature_relaxation_factor.resize(2, std::vector<Real>(_cht_boundary_names.size(), 1.0));
132 :
133 108 : for (const auto region_index : index_range(solid_fluid))
134 : {
135 : // First thing, we fetch the relaxation parameter values
136 : const auto & flux_param_value =
137 : getParam<std::vector<Real>>(flux_relaxation_param_names[region_index]);
138 72 : if (flux_param_value.empty() || (flux_param_value.size() != _cht_boundary_names.size()))
139 0 : paramError(flux_relaxation_param_names[region_index],
140 : "The number of relaxation factors is not the same as the number of interfaces!");
141 :
142 72 : _cht_flux_relaxation_factor[region_index] = flux_param_value;
143 : // We have to do the range check here because the intput parameter check errors if the vector is
144 : // empty
145 144 : for (const auto param : _cht_flux_relaxation_factor[region_index])
146 72 : if (param <= 0 || param > 1.0)
147 0 : paramError(flux_relaxation_param_names[region_index],
148 : "The relaxation parameter should be between 0 and 1!");
149 :
150 : const auto & temperature_param_value =
151 : getParam<std::vector<Real>>(temperature_relaxation_param_names[region_index]);
152 72 : if (temperature_param_value.empty() ||
153 : (temperature_param_value.size() != _cht_boundary_names.size()))
154 0 : paramError(temperature_relaxation_param_names[region_index],
155 : "The number of relaxation factors is not the same as the number of interfaces!");
156 :
157 72 : _cht_temperature_relaxation_factor[region_index] = temperature_param_value;
158 : // We have to do the range check here because the intput parameter check errors if the vector is
159 : // empty
160 144 : for (const auto param : _cht_temperature_relaxation_factor[region_index])
161 72 : if (param <= 0 || param > 1.0)
162 0 : paramError(temperature_relaxation_param_names[region_index],
163 : "The relaxation parameter should be between 0 and 1!");
164 :
165 : // We then fetch the conduction kernels
166 : std::vector<LinearFVFluxKernel *> flux_kernels;
167 72 : _app.theWarehouse()
168 72 : .query()
169 72 : .template condition<AttribSystem>("LinearFVFluxKernel")
170 144 : .template condition<AttribVar>(0)
171 72 : .template condition<AttribSysNum>(_cht_system_numbers[region_index])
172 : .queryInto(flux_kernels);
173 :
174 180 : for (auto kernel : flux_kernels)
175 : {
176 108 : auto check_diff = dynamic_cast<LinearFVDiffusion *>(kernel);
177 108 : auto check_aniso_diff = dynamic_cast<LinearFVAnisotropicDiffusion *>(kernel);
178 108 : if (_cht_conduction_kernels[region_index] && (check_diff || check_aniso_diff))
179 0 : mooseError("We already have a kernel that describes the heat conduction for the ",
180 : solid_fluid[region_index],
181 : " domain: ",
182 : _cht_conduction_kernels[region_index]->name(),
183 : " We found another one with the name: ",
184 : (check_diff ? check_diff->name() : check_aniso_diff->name()),
185 : " Make sure that you have only one conduction kernel on the ",
186 : solid_fluid[region_index],
187 : " side!");
188 :
189 108 : if (check_diff || check_aniso_diff)
190 72 : _cht_conduction_kernels[region_index] = kernel;
191 : }
192 :
193 : // Then we check the boundary conditions, to make sure at least there is something defined
194 : // from both sides
195 144 : for (const auto bd_index : index_range(_cht_boundary_names))
196 : {
197 : const auto & boundary_name = _cht_boundary_names[bd_index];
198 72 : const auto boundary_id = _cht_boundary_ids[bd_index];
199 :
200 : std::vector<LinearFVBoundaryCondition *> bcs;
201 72 : _problem.getMooseApp()
202 : .theWarehouse()
203 72 : .query()
204 72 : .template condition<AttribSystem>("LinearFVBoundaryCondition")
205 144 : .template condition<AttribVar>(0)
206 72 : .template condition<AttribSysNum>(_cht_system_numbers[region_index])
207 72 : .template condition<AttribBoundaries>(boundary_id)
208 : .queryInto(bcs);
209 :
210 72 : if (bcs.size() != 1)
211 0 : mooseError("We found multiple or no boundary conditions for solid energy on boundary ",
212 : boundary_name,
213 : " (ID: ",
214 : boundary_id,
215 : "). Make sure you define exactly one for conjugate heat transfer applications!");
216 72 : _cht_boundary_conditions[bd_index][region_index] = bcs[0];
217 :
218 72 : if (!dynamic_cast<LinearFVCHTBCInterface *>(_cht_boundary_conditions[bd_index][region_index]))
219 0 : mooseError("The selected boundary condition cannot be used with CHT problems! Make sure it "
220 : "inherits from LinearFVCHTBCInterface!");
221 72 : }
222 72 : }
223 216 : }
224 :
225 : void
226 36 : CHTHandler::setupConjugateHeatTransferContainers()
227 : {
228 : // We already error in initialSetup if we have more variables
229 : const auto * fluid_variable =
230 36 : dynamic_cast<const MooseLinearVariableFVReal *>(&_energy_system->getVariable(0, 0));
231 : const auto * solid_variable =
232 36 : dynamic_cast<const MooseLinearVariableFVReal *>(&_solid_energy_system->getVariable(0, 0));
233 :
234 36 : _cht_face_info.clear();
235 36 : _cht_face_info.resize(_cht_boundary_ids.size());
236 36 : _boundary_heat_flux.clear();
237 36 : _boundary_temperature.clear();
238 36 : _integrated_boundary_heat_flux.clear();
239 :
240 72 : for (const auto bd_index : index_range(_cht_boundary_ids))
241 : {
242 36 : const auto bd_id = _cht_boundary_ids[bd_index];
243 : const auto & bd_name = _cht_boundary_names[bd_index];
244 :
245 : // We populate the face infos for every interface
246 : auto & bd_fi_container = _cht_face_info[bd_index];
247 74242 : for (auto & fi : _problem.mesh().faceInfo())
248 74206 : if (fi->boundaryIDs().count(bd_id))
249 576 : bd_fi_container.push_back(fi);
250 :
251 : // We do this because the coupling functors should be evaluated on both sides
252 : // of the interface and there are rigorous checks if the functors don't support a subdomain
253 : std::set<SubdomainID> combined_set;
254 36 : std::set_union(solid_variable->blockIDs().begin(),
255 36 : solid_variable->blockIDs().end(),
256 36 : fluid_variable->blockIDs().begin(),
257 36 : fluid_variable->blockIDs().end(),
258 : std::inserter(combined_set, combined_set.begin()));
259 :
260 : // We instantiate the coupling fuctors for heat flux and temperature
261 : FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>> solid_bd_flux(
262 36 : _problem.mesh(), combined_set, "heat_flux_to_solid_" + bd_name);
263 : FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>> fluid_bd_flux(
264 36 : _problem.mesh(), combined_set, "heat_flux_to_fluid_" + bd_name);
265 :
266 : _boundary_heat_flux.push_back(
267 144 : std::vector<FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>>>(
268 72 : {std::move(solid_bd_flux), std::move(fluid_bd_flux)}));
269 : auto & flux_container = _boundary_heat_flux.back();
270 :
271 72 : _integrated_boundary_heat_flux.push_back(std::vector<Real>({0.0, 0.0}));
272 :
273 : FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>> solid_bd_temperature(
274 36 : _problem.mesh(), combined_set, "interface_temperature_solid_" + bd_name);
275 : FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>> fluid_bd_temperature(
276 36 : _problem.mesh(), combined_set, "interface_temperature_fluid_" + bd_name);
277 :
278 : _boundary_temperature.push_back(
279 144 : std::vector<FaceCenteredMapFunctor<Real, std::unordered_map<dof_id_type, Real>>>(
280 72 : {std::move(solid_bd_temperature), std::move(fluid_bd_temperature)}));
281 : auto & temperature_container = _boundary_temperature.back();
282 :
283 : // Time to register the functors on all of the threads
284 72 : for (const auto tid : make_range(libMesh::n_threads()))
285 : {
286 36 : _problem.addFunctor("heat_flux_to_solid_" + bd_name, flux_container[NS::CHTSide::SOLID], tid);
287 36 : _problem.addFunctor("heat_flux_to_fluid_" + bd_name, flux_container[NS::CHTSide::FLUID], tid);
288 36 : _problem.addFunctor(
289 36 : "interface_temperature_solid_" + bd_name, temperature_container[NS::CHTSide::SOLID], tid);
290 36 : _problem.addFunctor(
291 72 : "interface_temperature_fluid_" + bd_name, temperature_container[NS::CHTSide::FLUID], tid);
292 : }
293 :
294 : // Initialize the containers, they will be filled with correct values soon.
295 : // Before any solve happens.
296 108 : for (const auto region_index : make_range(2))
297 1224 : for (auto & fi : bd_fi_container)
298 : {
299 1152 : flux_container[region_index][fi->id()] = 0.0;
300 1152 : temperature_container[region_index][fi->id()] = 0.0;
301 : }
302 36 : }
303 108 : }
304 :
305 : void
306 36 : CHTHandler::initializeCHTCouplingFields()
307 : {
308 72 : for (const auto bd_index : index_range(_cht_boundary_ids))
309 : {
310 : const auto & bd_fi_container = _cht_face_info[bd_index];
311 : auto & temperature_container = _boundary_temperature[bd_index];
312 :
313 108 : for (const auto region_index : make_range(2))
314 : {
315 : // Can't be const considering we will update members from here
316 72 : auto bc = _cht_boundary_conditions[bd_index][region_index];
317 1224 : for (const auto & fi : bd_fi_container)
318 : {
319 1152 : bc->setupFaceData(fi, fi->faceType(std::make_pair(0, _cht_system_numbers[region_index])));
320 1152 : temperature_container[1 - region_index][fi->id()] = bc->computeBoundaryValue();
321 : }
322 : }
323 : }
324 36 : }
325 :
326 : void
327 15012 : CHTHandler::updateCHTBoundaryCouplingFields(const NS::CHTSide side)
328 : {
329 : // Well we can just use the face that this enum casts into int very nicely
330 : // we can use it to get the index of the other side
331 15012 : const NS::CHTSide other_side = static_cast<NS::CHTSide>(1 - side);
332 :
333 30024 : for (const auto bd_index : index_range(_cht_boundary_ids))
334 : {
335 15012 : auto & other_bc = _cht_boundary_conditions[bd_index][other_side];
336 : auto & other_kernel = _cht_conduction_kernels[other_side];
337 :
338 : // We get the relaxation from the other side, so if we are fluid side we get the solid
339 : // relaxation
340 15012 : const auto temperature_relaxation = _cht_flux_relaxation_factor[other_side][bd_index];
341 15012 : const auto flux_relaxation = _cht_temperature_relaxation_factor[other_side][bd_index];
342 :
343 : // Fetching the right container here, if side is fluid we fetch "heat_flux_to_fluid"
344 15012 : auto & flux_container = _boundary_heat_flux[bd_index][side];
345 : // Fetching the other side's contaienr here, if side is fluid we fetch the solid temperature
346 : auto & temperature_container = _boundary_temperature[bd_index][other_side];
347 : // We will also update the integrated flux for output info
348 : auto & integrated_flux = _integrated_boundary_heat_flux[bd_index][side];
349 : // We are recomputing this so, time to zero this out
350 15012 : integrated_flux = 0.0;
351 :
352 : const auto & bd_fi_container = _cht_face_info[bd_index];
353 :
354 : // We enter the face loop to update the coupling fields
355 255204 : for (const auto & fi : bd_fi_container)
356 : {
357 240192 : other_kernel->setupFaceData(fi);
358 : // We will want the flux in W/m2 for the coupling so no face integral for now,
359 : // this can cause issues if we start using face area in the kernels
360 : // for more than just face integral multipliers.
361 : // Also, if we decide to not require overlapping meshes on the boundary
362 : // this will probably have to change.
363 240192 : other_kernel->setCurrentFaceArea(1.0);
364 240192 : other_bc->setupFaceData(fi, fi->faceType(std::make_pair(0, _cht_system_numbers[other_side])));
365 :
366 : // T_new = relaxation * T_boundary + (1-relaxation) * T_old
367 480384 : temperature_container[fi->id()] =
368 240192 : temperature_relaxation * other_bc->computeBoundaryValue() +
369 240192 : (1 - temperature_relaxation) * temperature_container[fi->id()];
370 :
371 : // Flux_new = relaxation * Flux_boundary + (1-relaxation) * Flux_old,
372 : // minus sign is due to the normal differences
373 240192 : const auto flux = flux_relaxation * other_kernel->computeBoundaryFlux(*other_bc) +
374 240192 : (1 - flux_relaxation) * flux_container[fi->id()];
375 :
376 240192 : flux_container[fi->id()] = flux;
377 :
378 : // We do the integral here
379 240192 : integrated_flux += flux * fi->faceArea() * fi->faceCoord();
380 : }
381 : }
382 15012 : }
383 :
384 : void
385 7506 : CHTHandler::sumIntegratedFluxes()
386 : {
387 15012 : for (const auto i : index_range(_integrated_boundary_heat_flux))
388 : {
389 : auto & integrated_fluxes = _integrated_boundary_heat_flux[i];
390 7506 : _problem.comm().sum(integrated_fluxes[NS::CHTSide::SOLID]);
391 7506 : _problem.comm().sum(integrated_fluxes[NS::CHTSide::FLUID]);
392 : }
393 7506 : }
394 :
395 : void
396 7506 : CHTHandler::printIntegratedFluxes() const
397 : {
398 15012 : for (const auto i : index_range(_integrated_boundary_heat_flux))
399 : {
400 : auto & integrated_fluxes = _integrated_boundary_heat_flux[i];
401 7506 : _console << " Iteration " << _fpi_it << " Boundary " << _cht_boundary_names[i]
402 7506 : << " flux on solid side " << integrated_fluxes[NS::CHTSide::SOLID]
403 7506 : << " flux on fluid side: " << integrated_fluxes[NS::CHTSide::FLUID] << std::endl;
404 : }
405 7506 : }
406 :
407 : void
408 4212 : CHTHandler::resetIntegratedFluxes()
409 : {
410 8424 : for (const auto i : index_range(_integrated_boundary_heat_flux))
411 8424 : _integrated_boundary_heat_flux[i] = std::vector<Real>({0.0, 0.0});
412 4212 : }
413 :
414 : bool
415 67764 : CHTHandler::converged() const
416 : {
417 67764 : if (_fpi_it >= _max_cht_fpi)
418 : return true;
419 :
420 43917 : for (const auto & boundary_flux : _integrated_boundary_heat_flux)
421 : {
422 9594 : const Real f1 = boundary_flux[0];
423 9594 : const Real f2 = boundary_flux[1];
424 :
425 : // Special case: both are zero at startup not converged yet
426 9594 : if (_fpi_it != 0 && (f1 == 0.0 && f2 == 0.0))
427 : return true;
428 :
429 : // These fluxes should be of opposite sign
430 9594 : const Real diff = std::abs(f1 + f2);
431 9594 : const Real denom = std::max({std::fabs(f1), std::fabs(f2), Real(1e-14)});
432 9594 : const Real rel_diff = diff / denom;
433 :
434 9594 : if (rel_diff >= _cht_heat_flux_tolerance)
435 : return false;
436 : }
437 :
438 34323 : return _fpi_it;
439 : }
440 :
441 : } // End FV namespace
442 : } // End Moose namespace
|
