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