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 "FunctorThermalResistanceBC.h"
11 : #include "HeatConductionNames.h"
12 :
13 : registerMooseObject("HeatTransferApp", FunctorThermalResistanceBC);
14 :
15 : InputParameters
16 164 : FunctorThermalResistanceBC::validParams()
17 : {
18 164 : auto params = FVFluxBC::validParams();
19 328 : params.addParam<MooseFunctorName>("temperature", "temperature variable");
20 164 : params.addRequiredParam<Real>(HeatConduction::T_ambient, "constant ambient temperature");
21 328 : params.addRequiredParam<MooseFunctorName>("htc", "heat transfer coefficient");
22 :
23 328 : params.addRequiredRangeCheckedParam<Real>(HeatConduction::emissivity,
24 164 : HeatConduction::emissivity + " >= 0.0 & " +
25 164 : HeatConduction::emissivity + " <= 1.0",
26 : "emissivity of the surface");
27 :
28 328 : params.addRequiredParam<std::vector<Real>>(
29 : "thermal_conductivities",
30 : "vector of thermal conductivity values used for the conduction layers");
31 328 : params.addRequiredParam<std::vector<Real>>("conduction_thicknesses",
32 : "vector of conduction layer thicknesses");
33 :
34 328 : MooseEnum geometry("cartesian cylindrical", "cartesian");
35 328 : params.addParam<MooseEnum>("geometry", geometry, "type of geometry");
36 328 : params.addRangeCheckedParam<Real>("inner_radius",
37 : "inner_radius > 0.0",
38 : "coordinate corresponding to the first resistance layer");
39 :
40 492 : params.addRangeCheckedParam<Real>(
41 328 : "step_size", 0.1, "step_size > 0.0", "underrelaxation step size");
42 :
43 492 : params.addRangeCheckedParam<unsigned int>(
44 328 : "max_iterations", 100, "max_iterations >= 0", "maximum iterations");
45 :
46 492 : params.addRangeCheckedParam<Real>(
47 328 : "tolerance", 1E-3, "tolerance > 0.0", "tolerance to converge iterations");
48 164 : params.addClassDescription("Thermal resistance heat flux boundary condition for the "
49 : "fluid and solid energy equations");
50 164 : return params;
51 164 : }
52 :
53 88 : FunctorThermalResistanceBC::FunctorThermalResistanceBC(const InputParameters & parameters)
54 : : FVFluxBC(parameters),
55 88 : _geometry(getParam<MooseEnum>("geometry").getEnum<Moose::CoordinateSystemType>()),
56 88 : _inner_radius(
57 110 : _geometry == Moose::CoordinateSystemType::COORD_RZ ? getParam<Real>("inner_radius") : 1.0),
58 264 : _T(isParamValid("temperature") ? getFunctor<ADReal>("temperature")
59 176 : : getFunctor<ADReal>("variable")),
60 88 : _T_ambient(getParam<Real>(HeatConduction::T_ambient)),
61 176 : _k(getParam<std::vector<Real>>("thermal_conductivities")),
62 176 : _dx(getParam<std::vector<Real>>("conduction_thicknesses")),
63 176 : _h(getFunctor<ADReal>(getParam<MooseFunctorName>("htc"))),
64 88 : _emissivity(getParam<Real>(HeatConduction::emissivity)),
65 176 : _max_iterations(getParam<unsigned int>("max_iterations")),
66 176 : _tolerance(getParam<Real>("tolerance")),
67 176 : _alpha(getParam<Real>("step_size")),
68 88 : _T_surface(0.0),
69 88 : _outer_radius(_inner_radius),
70 88 : _conduction_resistance(0.0),
71 176 : _parallel_resistance(0.0)
72 : {
73 88 : if (_k.size() != _dx.size())
74 0 : paramError("conduction_thicknesses",
75 : "Number of specified thermal conductivities must match "
76 : "the number of conduction layers!");
77 :
78 88 : if (_geometry == Moose::CoordinateSystemType::COORD_RZ)
79 88 : for (const auto & d : _dx)
80 : _outer_radius += d;
81 :
82 : // because the thermal conductivities are constant, we only need to compute
83 : // the conduction resistance one time
84 88 : computeConductionResistance();
85 88 : }
86 :
87 : void
88 88 : FunctorThermalResistanceBC::computeConductionResistance()
89 : {
90 88 : Real r = _inner_radius;
91 :
92 352 : for (const auto i : index_range(_k))
93 : {
94 264 : switch (_geometry)
95 : {
96 198 : case Moose::CoordinateSystemType::COORD_XYZ:
97 198 : _conduction_resistance += _dx[i] / _k[i];
98 198 : break;
99 66 : case Moose::CoordinateSystemType::COORD_RZ:
100 : {
101 66 : _conduction_resistance += std::log((_dx[i] + r) / r) / _k[i];
102 66 : r += _dx[i];
103 66 : break;
104 : }
105 0 : default:
106 0 : mooseError("Unhandled 'GeometryEnum' in 'FunctorThermalResistanceBC'!");
107 : }
108 : }
109 88 : }
110 :
111 : ADReal
112 900 : FunctorThermalResistanceBC::computeQpResidual()
113 : {
114 : // Evaluate material properties on the face
115 900 : const auto face_arg = singleSidedFaceArg();
116 :
117 : // radiation resistance has to be solved iteratively, since we don't know the
118 : // surface temperature. We do know that the heat flux in the conduction layers
119 : // must match the heat flux in the parallel convection-radiation segment. For a
120 : // first guess, take the surface temperature as the average of _T and T_ambient.
121 2700 : _T_surface = 0.5 * (_T(face_arg, determineState()) + _T_ambient);
122 :
123 : // total flux perpendicular to boundary
124 : ADReal flux;
125 :
126 : // resistance component representing the sum of the convection and radiation
127 : // resistances, in parallel
128 900 : computeParallelResistance();
129 :
130 : // other iteration requirements
131 : unsigned int iteration = 0;
132 900 : ADReal norm = 2 * _tolerance;
133 : ADReal T_surface_previous;
134 :
135 : // iterate to find the approximate surface temperature needed for evaluating the
136 : // radiation resistance. We only do this iteration if we have radiation transfer.
137 900 : if (_emissivity > 1e-8)
138 54450 : while (norm > (_tolerance * _alpha))
139 : {
140 54450 : T_surface_previous = _T_surface;
141 :
142 : // compute the flux based on the conduction part of the circuit
143 108900 : flux = (_T(face_arg, determineState()) - _T_surface) / _conduction_resistance;
144 :
145 54450 : computeParallelResistance();
146 :
147 : // use the flux computed from the conduction half to update T_surface
148 108900 : _T_surface = flux * _parallel_resistance + _T_ambient;
149 163350 : _T_surface = _alpha * _T_surface + (1 - _alpha) * T_surface_previous;
150 163350 : norm = std::abs(_T_surface - T_surface_previous) / std::abs(T_surface_previous);
151 :
152 54450 : if (iteration == _max_iterations)
153 : {
154 450 : mooseWarning("Maximum number of iterations reached in 'FunctorThermalResistanceBC'!");
155 : break;
156 : }
157 : else
158 54000 : iteration += 1;
159 : }
160 :
161 : // once we have determined T_surface, we can finally evaluate the complete
162 : // resistance to find the overall heat flux. For Cartesian, dividing by the
163 : // 'inner_radius' has no effect, but it is required for correct normalization
164 : // for cylindrical geometries.
165 900 : flux = (_T(face_arg, determineState()) - _T_ambient) /
166 1800 : (_conduction_resistance + _parallel_resistance) / _inner_radius;
167 900 : return flux;
168 : }
169 :
170 : void
171 55350 : FunctorThermalResistanceBC::computeParallelResistance()
172 : {
173 55350 : const auto face_arg = singleSidedFaceArg();
174 :
175 : // compute the parallel convection and radiation resistances, assuming they
176 : // act on the same surface area size
177 55350 : ADReal hr = _emissivity * HeatConduction::Constants::sigma *
178 55350 : (_T_surface * _T_surface + _T_ambient * _T_ambient) * (_T_surface + _T_ambient);
179 :
180 : // for Cartesian, dividing by the 'outer_radius' has no effect, but it is
181 : // required for correct normalization for cylindrical geometries
182 166050 : _parallel_resistance = 1.0 / (hr + _h(face_arg, determineState())) / _outer_radius;
183 55350 : }
|
