OGS 6.2.1-372-g67eca524a.dirty.20190919230349
BHE_1U.cpp
Go to the documentation of this file.
1 
11 #include "BHE_1U.h"
12 
13 #include <boost/math/constants/constants.hpp>
15 #include "Physics.h"
17 
18 namespace ProcessLib
19 {
20 namespace HeatTransportBHE
21 {
22 namespace BHE
23 {
24 std::array<double, BHE_1U::number_of_unknowns> BHE_1U::pipeHeatCapacities()
25  const
26 {
27  double const rho_r = refrigerant.density;
29  double const rho_g = grout.rho_g;
30  double const porosity_g = grout.porosity_g;
31  double const heat_cap_g = grout.heat_cap_g;
32 
33  return {{/*i1*/ rho_r * specific_heat_capacity,
34  /*o1*/ rho_r * specific_heat_capacity,
35  /*g1*/ (1.0 - porosity_g) * rho_g * heat_cap_g,
36  /*g2*/ (1.0 - porosity_g) * rho_g * heat_cap_g}};
37 }
38 
39 std::array<double, BHE_1U::number_of_unknowns> BHE_1U::pipeHeatConductions()
40  const
41 {
42  double const lambda_r = refrigerant.thermal_conductivity;
43  double const rho_r = refrigerant.density;
44  double const Cp_r = refrigerant.specific_heat_capacity;
45  double const alpha_L = _pipes.longitudinal_dispersion_length;
46  double const porosity_g = grout.porosity_g;
47  double const lambda_g = grout.lambda_g;
48 
49  double const velocity_norm = std::abs(_flow_velocity) * std::sqrt(2);
50 
51  // Here we calculate the laplace coefficients in the governing
52  // equations of BHE. These governing equations can be found in
53  // 1) Diersch (2013) FEFLOW book on page 952, M.120-122, or
54  // 2) Diersch (2011) Comp & Geosci 37:1122-1135, Eq. 19-22.
55  return {{// pipe i1, Eq. 19
56  (lambda_r + rho_r * Cp_r * alpha_L * velocity_norm),
57  // pipe o1, Eq. 20
58  (lambda_r + rho_r * Cp_r * alpha_L * velocity_norm),
59  // pipe g1, Eq. 21
60  (1.0 - porosity_g) * lambda_g,
61  // pipe g2, Eq. 22
62  (1.0 - porosity_g) * lambda_g}};
63 }
64 
65 std::array<Eigen::Vector3d, BHE_1U::number_of_unknowns>
67 {
68  double const& rho_r = refrigerant.density;
69  double const& Cp_r = refrigerant.specific_heat_capacity;
70 
71  return {{// pipe i1, Eq. 19
72  {0, 0, -rho_r * Cp_r * _flow_velocity},
73  // pipe o1, Eq. 20
74  {0, 0, rho_r * Cp_r * _flow_velocity},
75  // grout g1, Eq. 21
76  {0, 0, 0},
77  // grout g2, Eq. 22
78  {0, 0, 0}}};
79 }
80 
81 double compute_R_gs(double const chi, double const R_g)
82 {
83  return (1 - chi) * R_g;
84 }
85 
86 double compute_R_gg(double const chi, double const R_gs, double const R_ar,
87  double const R_g)
88 {
89  double const R_gg = 2.0 * R_gs * (R_ar - 2.0 * chi * R_g) /
90  (2.0 * R_gs - R_ar + 2.0 * chi * R_g);
91  if (!std::isfinite(R_gg))
92  {
93  OGS_FATAL(
94  "Error!!! Grout Thermal Resistance is an infinite number! The "
95  "simulation will be stopped!");
96  }
97 
98  return R_gg;
99 }
100 
106 std::pair<double, double> thermalResistancesGroutSoil(double chi,
107  double const R_ar,
108  double const R_g)
109 {
110  double R_gs = compute_R_gs(chi, R_g);
111  double R_gg =
112  compute_R_gg(chi, R_gs, R_ar, R_g); // Resulting thermal resistances.
113 
114  auto constraint = [&]() {
115  return 1.0 / ((1.0 / R_gg) + (1.0 / (2.0 * R_gs)));
116  };
117 
118  std::array<double, 3> const multiplier{chi * 0.66, chi * 0.5 * 0.66, 0.0};
119  for (double m_chi : multiplier)
120  {
121  if (constraint() >= 0)
122  {
123  break;
124  }
125  DBUG(
126  "Warning! Correction procedure was applied due to negative thermal "
127  "resistance! Chi = %f.\n",
128  m_chi);
129 
130  R_gs = compute_R_gs(m_chi, R_g);
131  R_gg = compute_R_gg(m_chi, R_gs, R_ar, R_g);
132  }
133 
134  return {R_gg, R_gs};
135 }
136 
137 void BHE_1U::updateHeatTransferCoefficients(double const flow_rate)
138 
139 {
140  auto const tm_flow_properties = calculateThermoMechanicalFlowPropertiesPipe(
142 
143  _flow_velocity = tm_flow_properties.velocity;
145  calcThermalResistances(tm_flow_properties.nusselt_number);
146 }
147 
149 std::array<double, BHE_1U::number_of_unknowns> BHE_1U::calcThermalResistances(
150  double const Nu)
151 {
152  constexpr double pi = boost::math::constants::pi<double>();
153 
154  double const lambda_r = refrigerant.thermal_conductivity;
155  double const lambda_g = grout.lambda_g;
156  double const lambda_p = _pipes.inlet.wall_thermal_conductivity;
157 
158  // thermal resistances due to advective flow of refrigerant in the _pipes
159  // Eq. 36 in Diersch_2011_CG
160  double const R_adv_i1 = 1.0 / (Nu * lambda_r * pi);
161  double const R_adv_o1 = 1.0 / (Nu * lambda_r * pi);
162 
163  // thermal resistance due to thermal conductivity of the pipe wall material
164  // Eq. 49
165  double const R_con_a =
166  std::log(_pipes.outlet.diameter / _pipes.inlet.diameter) /
167  (2.0 * pi * lambda_p);
168 
169  // the average outer diameter of the _pipes
170  double const d0 = _pipes.outlet.diameter;
171  double const D = borehole_geometry.diameter;
172  // Eq. 51
173  double const chi = std::log(std::sqrt(D * D + 2 * d0 * d0) / 2 / d0) /
174  std::log(D / std::sqrt(2) / d0);
175  // Eq. 52
176  // thermal resistances of the grout
177  double const R_g =
178  std::acosh((D * D + d0 * d0 - _pipes.distance * _pipes.distance) /
179  (2 * D * d0)) /
180  (2 * pi * lambda_g) * (1.601 - 0.888 * _pipes.distance / D);
181  // thermal resistance due to the grout transition.
182  double const R_con_b = chi * R_g;
183  // Eq. 29 and 30
184  double const R_fig = R_adv_i1 + R_con_a + R_con_b;
185  double const R_fog = R_adv_o1 + R_con_a + R_con_b;
186 
187  // thermal resistance due to inter-grout exchange
188  double const R_ar =
189  std::acosh((2.0 * _pipes.distance * _pipes.distance - d0 * d0) / d0 /
190  d0) /
191  (2.0 * pi * lambda_g);
192 
193  double R_gg;
194  double R_gs;
195  std::tie(R_gg, R_gs) = thermalResistancesGroutSoil(chi, R_ar, R_g);
196 
197  return {{R_fig, R_fog, R_gg, R_gs}};
198 
199  // keep the following lines------------------------------------------------
200  // when debugging the code, printing the R and phi values are needed--------
201  // std::cout << "Rfig =" << R_fig << " Rfog =" << R_fog << " Rgg =" <<
202  // R_gg << " Rgs =" << R_gs << "\n"; double phi_fig = 1.0 / (R_fig *
203  // S_i); double phi_fog = 1.0 / (R_fog * S_o); double phi_gg = 1.0 / (R_gg
204  // * S_g1); double phi_gs = 1.0 / (R_gs * S_gs); std::cout << "phi_fig ="
205  // << phi_fig << " phi_fog =" << phi_fog << " phi_gg =" << phi_gg << "
206  // phi_gs =" << phi_gs << "\n";
207  // -------------------------------------------------------------------------
208 }
209 
210 double BHE_1U::updateFlowRateAndTemperature(double const T_out,
211  double const current_time)
212 {
213  auto values =
214  visit([&](auto const& control) { return control(T_out, current_time); },
216  updateHeatTransferCoefficients(values.flow_rate);
217  return values.temperature;
218 }
219 } // namespace BHE
220 } // namespace HeatTransportBHE
221 } // namespace ProcessLib
std::array< double, number_of_unknowns > pipeHeatConductions() const
Definition: BHE_1U.cpp:39
double const wall_thermal_conductivity
Definition: Pipe.h:29
std::array< double, number_of_unknowns > calcThermalResistances(double const Nu)
Nu is the Nusselt number.
Definition: BHE_1U.cpp:149
double updateFlowRateAndTemperature(double T_out, double current_time)
Return the inflow temperature for the boundary condition.
Definition: BHE_1U.cpp:210
ThermoMechanicalFlowProperties calculateThermoMechanicalFlowPropertiesPipe(Pipe const &pipe, double const length, RefrigerantProperties const &fluid, double const flow_rate)
BoreholeGeometry const borehole_geometry
Definition: BHECommon.h:41
double compute_R_gs(double const chi, double const R_g)
Definition: BHE_1U.cpp:81
std::pair< double, double > thermalResistancesGroutSoil(double chi, double const R_ar, double const R_g)
Definition: BHE_1U.cpp:106
void updateHeatTransferCoefficients(double const flow_rate)
Definition: BHE_1U.cpp:137
FlowAndTemperatureControl const flowAndTemperatureControl
Definition: BHECommon.h:44
std::array< double, number_of_unknowns > _thermal_resistances
Definition: BHE_1U.h:175
std::array< double, number_of_unknowns > pipeHeatCapacities() const
Definition: BHE_1U.cpp:24
std::array< Eigen::Vector3d, number_of_unknowns > pipeAdvectionVectors() const
Definition: BHE_1U.cpp:66
#define OGS_FATAL(fmt,...)
Definition: Error.h:64
double _flow_velocity
Flow velocity inside the pipes. Depends on the flow_rate.
double compute_R_gg(double const chi, double const R_gs, double const R_ar, double const R_g)
Definition: BHE_1U.cpp:86
RefrigerantProperties const refrigerant
Definition: BHECommon.h:42