OGS
BHE_2U.cpp
Go to the documentation of this file.
1
11#include "BHE_2U.h"
12
13#include <boost/math/constants/constants.hpp>
14
16#include "Physics.h"
18
19namespace ProcessLib
20{
21namespace HeatTransportBHE
22{
23namespace BHE
24{
26 RefrigerantProperties const& refrigerant,
27 GroutParameters const& grout,
28 FlowAndTemperatureControl const& flowAndTemperatureControl,
29 PipeConfigurationUType const& pipes,
30 bool const use_python_bcs)
31 : BHECommonUType{borehole, refrigerant, grout, flowAndTemperatureControl,
32 pipes, use_python_bcs}
33{
34 _thermal_resistances.fill(std::numeric_limits<double>::quiet_NaN());
35
36 // Initialize thermal resistances.
37 auto values = visit(
38 [&](auto const& control) {
40 0. /* initial time */);
41 },
43 updateHeatTransferCoefficients(values.flow_rate);
44}
45
46std::array<double, BHE_2U::number_of_unknowns> BHE_2U::pipeHeatCapacities()
47 const
48{
49 double const rho_r = refrigerant.density;
50 double const specific_heat_capacity = refrigerant.specific_heat_capacity;
51 double const rho_g = grout.rho_g;
52 double const porosity_g = grout.porosity_g;
53 double const heat_cap_g = grout.heat_cap_g;
54
55 return {{/*i1*/ rho_r * specific_heat_capacity,
56 /*i2*/ rho_r * specific_heat_capacity,
57 /*o1*/ rho_r * specific_heat_capacity,
58 /*o2*/ rho_r * specific_heat_capacity,
59 /*g1*/ (1.0 - porosity_g) * rho_g * heat_cap_g,
60 /*g2*/ (1.0 - porosity_g) * rho_g * heat_cap_g,
61 /*g3*/ (1.0 - porosity_g) * rho_g * heat_cap_g,
62 /*g4*/ (1.0 - porosity_g) * rho_g * heat_cap_g}};
63}
64
65std::array<double, BHE_2U::number_of_unknowns> BHE_2U::pipeHeatConductions()
66 const
67{
68 double const lambda_r = refrigerant.thermal_conductivity;
69 double const rho_r = refrigerant.density;
70 double const Cp_r = refrigerant.specific_heat_capacity;
71 double const alpha_L = _pipes.longitudinal_dispersion_length;
72 double const porosity_g = grout.porosity_g;
73 double const lambda_g = grout.lambda_g;
74
75 // Here we calculate the laplace coefficients in the governing
76 // equations of BHE. These governing equations can be found in
77 // 1) Diersch (2013) FEFLOW book on page 952, M.120-122, or
78 // 2) Diersch (2011) Comp & Geosci 37:1122-1135, Eq. 19-22.
79 return {{// pipe i1
80 (lambda_r + rho_r * Cp_r * alpha_L * _flow_velocity),
81 // pipe i2
82 (lambda_r + rho_r * Cp_r * alpha_L * _flow_velocity),
83 // pipe o1
84 (lambda_r + rho_r * Cp_r * alpha_L * _flow_velocity),
85 // pipe o2
86 (lambda_r + rho_r * Cp_r * alpha_L * _flow_velocity),
87 // pipe g1
88 (1.0 - porosity_g) * lambda_g,
89 // pipe g2
90 (1.0 - porosity_g) * lambda_g,
91 // pipe g3
92 (1.0 - porosity_g) * lambda_g,
93 // pipe g4
94 (1.0 - porosity_g) * lambda_g}};
95}
96
97std::array<Eigen::Vector3d, BHE_2U::number_of_unknowns>
98BHE_2U::pipeAdvectionVectors(Eigen::Vector3d const& /*elem_direction*/) const
99{
100 double const rho_r = refrigerant.density;
101 double const Cp_r = refrigerant.specific_heat_capacity;
102
103 return {{// pipe i1
104 {0, 0, -rho_r * Cp_r * _flow_velocity},
105 // pipe i2
106 {0, 0, -rho_r * Cp_r * _flow_velocity},
107 // pipe o1
108 {0, 0, rho_r * Cp_r * _flow_velocity},
109 // pipe o2
110 {0, 0, rho_r * Cp_r * _flow_velocity},
111 // grout g1
112 {0, 0, 0},
113 // grout g2
114 {0, 0, 0},
115 // grout g3
116 {0, 0, 0},
117 // grout g4
118 {0, 0, 0}}};
119}
120
121double compute_R_gs_2U(double const chi, double const R_g)
122{
123 return (1 - chi) * R_g;
124}
125
126double compute_R_gg_2U(double const chi, double const R_gs, double const R_ar,
127 double const R_g)
128{
129 double const R_gg = 2.0 * R_gs * (R_ar - 2.0 * chi * R_g) /
130 (2.0 * R_gs - R_ar + 2.0 * chi * R_g);
131 if (!std::isfinite(R_gg))
132 {
133 OGS_FATAL(
134 "Error!!! Grout Thermal Resistance is an infinite number! The "
135 "simulation will be stopped!");
136 }
137
138 return R_gg;
139}
140
146std::array<double, 4> thermalResistancesGroutSoil2U(double const chi,
147 double const R_ar_1,
148 double const R_ar_2,
149 double const R_g)
150{
151 double R_gs = compute_R_gs_2U(chi, R_g);
152 double R_gg_1 = compute_R_gg_2U(chi, R_gs, R_ar_1, R_g);
153 double R_gg_2 = compute_R_gg_2U(chi, R_gs, R_ar_2,
154 R_g); // Resulting thermal resistances.
155 double chi_new = chi;
156
157 auto constraint = [&]()
158 { return 1.0 / ((1.0 / R_gg_1) + (1.0 / (2.0 * R_gs))); };
159
160 std::array<double, 3> const multiplier{chi * 2.0 / 3.0, chi * 1.0 / 3.0,
161 0.0};
162 for (double m_chi : multiplier)
163 {
164 if (constraint() >= 0)
165 {
166 break;
167 }
168 DBUG(
169 "Warning! Correction procedure was applied due to negative thermal "
170 "resistance! Chi = {:f}.\n",
171 m_chi);
172 R_gs = compute_R_gs_2U(m_chi, R_g);
173 R_gg_1 = compute_R_gg_2U(m_chi, R_gs, R_ar_1, R_g);
174 R_gg_2 = compute_R_gg_2U(m_chi, R_gs, R_ar_2, R_g);
175 chi_new = m_chi;
176 }
177
178 return {chi_new, R_gg_1, R_gg_2, R_gs};
179}
180
181void BHE_2U::updateHeatTransferCoefficients(double const flow_rate)
182
183{
184 auto const tm_flow_properties = calculateThermoMechanicalFlowPropertiesPipe(
186
187 _flow_velocity = tm_flow_properties.velocity;
189 calcThermalResistances(tm_flow_properties.nusselt_number);
190}
191
193std::array<double, BHE_2U::number_of_unknowns> BHE_2U::calcThermalResistances(
194 double const Nu)
195{
196 constexpr double pi = boost::math::constants::pi<double>();
197
198 double const lambda_r = refrigerant.thermal_conductivity;
199 double const lambda_g = grout.lambda_g;
200 double const lambda_p = _pipes.inlet.wall_thermal_conductivity;
201
202 // thermal resistances due to advective flow of refrigerant in the _pipes
203 // Eq. 36 in Diersch_2011_CG
204 double const R_adv_i = 1.0 / (Nu * lambda_r * pi);
205 double const R_adv_o = 1.0 / (Nu * lambda_r * pi);
206
207 // thermal resistance due to thermal conductivity of the pipe wall material
208 // Eq. 49
209 double const R_con_a =
211 (2.0 * pi * lambda_p);
212
213 // the average outer diameter of the _pipes
214 double const d0 = _pipes.outlet.outsideDiameter();
215 double const D = borehole_geometry.diameter;
216 // Eq. 38
217 double const chi =
218 std::log(std::sqrt(D * D + 4 * d0 * d0) / 2 / std::sqrt(2) / d0) /
219 std::log(D / 2 / d0);
220 // Eq. 39
221 // thermal resistances of the grout
222 double const R_g =
223 std::acosh((D * D + d0 * d0 - 2 * _pipes.distance * _pipes.distance) /
224 (2 * D * d0)) /
225 (2 * pi * lambda_g) *
226 (3.098 - 4.432 * std::sqrt(2) * _pipes.distance / D +
227 2.364 * 2 * _pipes.distance * _pipes.distance / D / D);
228
229 // thermal resistance due to inter-grout exchange
230 double const R_ar_1 =
231 std::acosh((2.0 * _pipes.distance * _pipes.distance - d0 * d0) / d0 /
232 d0) /
233 (2.0 * pi * lambda_g);
234
235 double const R_ar_2 =
236 std::acosh((2.0 * 2.0 * _pipes.distance * _pipes.distance - d0 * d0) /
237 d0 / d0) /
238 (2.0 * pi * lambda_g);
239
240 auto const [chi_new, R_gg_1, R_gg_2, R_gs] =
241 thermalResistancesGroutSoil2U(chi, R_ar_1, R_ar_2, R_g);
242
243 // thermal resistance due to the grout transition.
244 double const R_con_b = chi_new * R_g;
245
246 // Eq. 29 and 30
247 double const R_fig = R_adv_i + R_con_a + R_con_b;
248 double const R_fog = R_adv_o + R_con_a + R_con_b;
249
250 return {{R_fig, R_fog, R_gg_1, R_gg_2, R_gs}};
251
252 // keep the following lines------------------------------------------------
253 // when debugging the code, printing the R and phi values are needed--------
254 // std::cout << "Rfig =" << R_fig << " Rfog =" << R_fog << " Rgg =" <<
255 // R_gg << " Rgs =" << R_gs << "\n"; double phi_fig = 1.0 / (R_fig *
256 // S_i); double phi_fog = 1.0 / (R_fog * S_o); double phi_gg = 1.0 / (R_gg
257 // * S_g1); double phi_gs = 1.0 / (R_gs * S_gs); std::cout << "phi_fig ="
258 // << phi_fig << " phi_fog =" << phi_fog << " phi_gg =" << phi_gg << "
259 // phi_gs =" << phi_gs << "\n";
260 // -------------------------------------------------------------------------
261}
262
263std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>
265 std::size_t const top_node_id,
266 std::size_t const /*bottom_node_id*/,
267 int const in_component_id)
268{
269 return {std::make_pair(top_node_id, in_component_id),
270 std::make_pair(top_node_id, in_component_id + 2)};
271}
272
273std::optional<
274 std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>>
276 std::size_t const bottom_node_id,
277 int const in_component_id,
278 int const out_component_id)
279{
280 return {{std::make_pair(bottom_node_id, in_component_id),
281 std::make_pair(bottom_node_id, out_component_id)}};
282}
283
284std::array<double, BHE_2U::number_of_unknowns> BHE_2U::crossSectionAreas() const
285{
286 return {{
287 _pipes.inlet.area(),
288 _pipes.inlet.area(),
295 }};
296}
297
298double BHE_2U::updateFlowRateAndTemperature(double const T_out,
299 double const current_time)
300{
301 auto values =
302 visit([&](auto const& control) { return control(T_out, current_time); },
304 updateHeatTransferCoefficients(values.flow_rate);
305 return values.temperature;
306}
307} // namespace BHE
308} // namespace HeatTransportBHE
309} // namespace ProcessLib
#define OGS_FATAL(...)
Definition Error.h:26
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
double _flow_velocity
Flow velocity inside the pipes. Depends on the flow_rate.
static std::array< std::pair< std::size_t, int >, 2 > getBHEInflowDirichletBCNodesAndComponents(std::size_t const top_node_id, std::size_t const, int const in_component_id)
Definition BHE_2U.cpp:264
std::array< double, number_of_unknowns > calcThermalResistances(double const Nu)
Nu is the Nusselt number.
Definition BHE_2U.cpp:193
void updateHeatTransferCoefficients(double const flow_rate)
Definition BHE_2U.cpp:181
std::array< double, number_of_unknowns > _thermal_resistances
Definition BHE_2U.h:231
std::array< double, number_of_unknowns > pipeHeatConductions() const
Definition BHE_2U.cpp:65
std::array< Eigen::Vector3d, number_of_unknowns > pipeAdvectionVectors(Eigen::Vector3d const &) const
Definition BHE_2U.cpp:98
double updateFlowRateAndTemperature(double T_out, double current_time)
Return the inflow temperature for the boundary condition.
Definition BHE_2U.cpp:298
std::array< double, number_of_unknowns > pipeHeatCapacities() const
Definition BHE_2U.cpp:46
BHE_2U(BoreholeGeometry const &borehole, RefrigerantProperties const &refrigerant, GroutParameters const &grout, FlowAndTemperatureControl const &flowAndTemperatureControl, PipeConfigurationUType const &pipes, bool const use_python_bcs)
Definition BHE_2U.cpp:25
std::array< double, number_of_unknowns > crossSectionAreas() const
Definition BHE_2U.cpp:284
static std::optional< std::array< std::pair< std::size_t, int >, 2 > > getBHEBottomDirichletBCNodesAndComponents(std::size_t const bottom_node_id, int const in_component_id, int const out_component_id)
Definition BHE_2U.cpp:275
std::variant< TemperatureCurveConstantFlow, TemperatureCurveFlowCurve, FixedPowerConstantFlow, FixedPowerFlowCurve, PowerCurveConstantFlow, PowerCurveFlowCurve, BuildingPowerCurveConstantFlow > FlowAndTemperatureControl
double compute_R_gg_2U(double const chi, double const R_gs, double const R_ar, double const R_g)
Definition BHE_2U.cpp:126
double compute_R_gs_2U(double const chi, double const R_g)
Definition BHE_2U.cpp:121
ThermoMechanicalFlowProperties calculateThermoMechanicalFlowPropertiesPipe(Pipe const &pipe, double const length, RefrigerantProperties const &fluid, double const flow_rate)
std::array< double, 4 > thermalResistancesGroutSoil2U(double const chi, double const R_ar_1, double const R_ar_2, double const R_g)
Definition BHE_2U.cpp:146
RefrigerantProperties const refrigerant
Definition BHECommon.h:42
FlowAndTemperatureControl const flowAndTemperatureControl
Definition BHECommon.h:44
double outsideArea() const
Area of the pipe's outside including the wall thickness.
Definition Pipe.h:35
double area() const
Area of the pipe's inside without the wall.
Definition Pipe.h:32