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 : using std::abs;
115 :
116 : // Evaluate material properties on the face
117 900 : const auto face_arg = singleSidedFaceArg();
118 :
119 : // radiation resistance has to be solved iteratively, since we don't know the
120 : // surface temperature. We do know that the heat flux in the conduction layers
121 : // must match the heat flux in the parallel convection-radiation segment. For a
122 : // first guess, take the surface temperature as the average of _T and T_ambient.
123 2700 : _T_surface = 0.5 * (_T(face_arg, determineState()) + _T_ambient);
124 :
125 : // total flux perpendicular to boundary
126 : ADReal flux;
127 :
128 : // resistance component representing the sum of the convection and radiation
129 : // resistances, in parallel
130 900 : computeParallelResistance();
131 :
132 : // other iteration requirements
133 : unsigned int iteration = 0;
134 900 : ADReal norm = 2 * _tolerance;
135 : ADReal T_surface_previous;
136 :
137 : // iterate to find the approximate surface temperature needed for evaluating the
138 : // radiation resistance. We only do this iteration if we have radiation transfer.
139 900 : if (_emissivity > 1e-8)
140 54450 : while (norm > (_tolerance * _alpha))
141 : {
142 54450 : T_surface_previous = _T_surface;
143 :
144 : // compute the flux based on the conduction part of the circuit
145 108900 : flux = (_T(face_arg, determineState()) - _T_surface) / _conduction_resistance;
146 :
147 54450 : computeParallelResistance();
148 :
149 : // use the flux computed from the conduction half to update T_surface
150 108900 : _T_surface = flux * _parallel_resistance + _T_ambient;
151 163350 : _T_surface = _alpha * _T_surface + (1 - _alpha) * T_surface_previous;
152 163350 : norm = abs(_T_surface - T_surface_previous) / abs(T_surface_previous);
153 :
154 54450 : if (iteration == _max_iterations)
155 : {
156 450 : mooseWarning("Maximum number of iterations reached in 'FunctorThermalResistanceBC'!");
157 : break;
158 : }
159 : else
160 54000 : iteration += 1;
161 : }
162 :
163 : // once we have determined T_surface, we can finally evaluate the complete
164 : // resistance to find the overall heat flux. For Cartesian, dividing by the
165 : // 'inner_radius' has no effect, but it is required for correct normalization
166 : // for cylindrical geometries.
167 900 : flux = (_T(face_arg, determineState()) - _T_ambient) /
168 1800 : (_conduction_resistance + _parallel_resistance) / _inner_radius;
169 900 : return flux;
170 : }
171 :
172 : void
173 55350 : FunctorThermalResistanceBC::computeParallelResistance()
174 : {
175 55350 : const auto face_arg = singleSidedFaceArg();
176 :
177 : // compute the parallel convection and radiation resistances, assuming they
178 : // act on the same surface area size
179 55350 : ADReal hr = _emissivity * HeatConduction::Constants::sigma *
180 55350 : (_T_surface * _T_surface + _T_ambient * _T_ambient) * (_T_surface + _T_ambient);
181 :
182 : // for Cartesian, dividing by the 'outer_radius' has no effect, but it is
183 : // required for correct normalization for cylindrical geometries
184 166050 : _parallel_resistance = 1.0 / (hr + _h(face_arg, determineState())) / _outer_radius;
185 55350 : }
|
