diff --git a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
index f21fa05b1cde55a59349641c267e10a3e0705d4a..8e4215cecc334825633064306fa63b61e0e827c5 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
@@ -23,11 +23,11 @@ namespace BHE
std::array<double, BHECommonCoaxial::number_of_unknowns>
BHECommonCoaxial::pipeHeatCapacities() const
{
- double const& rho_r = refrigerant.density;
- double const& specific_heat_capacity = refrigerant.specific_heat_capacity;
- double const& rho_g = grout.rho_g;
- double const& porosity_g = grout.porosity_g;
- double const& heat_cap_g = grout.heat_cap_g;
+ double const rho_r = refrigerant.density;
+ double const specific_heat_capacity = refrigerant.specific_heat_capacity;
+ double const rho_g = grout.rho_g;
+ double const porosity_g = grout.porosity_g;
+ double const heat_cap_g = grout.heat_cap_g;
return {{/*i*/ rho_r * specific_heat_capacity,
/*o*/ rho_r * specific_heat_capacity,
@@ -46,12 +46,12 @@ double BHECommonCoaxial::updateFlowRateAndTemperature(double const T_out,
std::array<double, BHECommonCoaxial::number_of_unknowns>
BHECommonCoaxial::pipeHeatConductions() const
{
- double const& lambda_r = refrigerant.thermal_conductivity;
- double const& rho_r = refrigerant.density;
- double const& Cp_r = refrigerant.specific_heat_capacity;
- double const& alpha_L = _pipes.longitudinal_dispersion_length;
- double const& porosity_g = grout.porosity_g;
- double const& lambda_g = grout.lambda_g;
+ double const lambda_r = refrigerant.thermal_conductivity;
+ double const rho_r = refrigerant.density;
+ double const Cp_r = refrigerant.specific_heat_capacity;
+ double const alpha_L = _pipes.longitudinal_dispersion_length;
+ double const porosity_g = grout.porosity_g;
+ double const lambda_g = grout.lambda_g;
auto v = velocities();
// Here we calculate the laplace coefficients in the governing
@@ -69,8 +69,8 @@ BHECommonCoaxial::pipeHeatConductions() const
std::array<Eigen::Vector3d, BHECommonCoaxial::number_of_unknowns>
BHECommonCoaxial::pipeAdvectionVectors() const
{
- double const& rho_r = refrigerant.density;
- double const& Cp_r = refrigerant.specific_heat_capacity;
+ double const rho_r = refrigerant.density;
+ double const Cp_r = refrigerant.specific_heat_capacity;
auto v = velocities();
return {{// pipe i, Eq. 26 and Eq. 23
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
index 67741a3c0f37cabfcda94e9f4a62b11b429a77a8..a6d0d2e92582f1387e63f33d11f7daf5d6be90e6 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
@@ -36,6 +36,8 @@ public:
_pipes.outer_pipe.area() - _pipes.inner_pipe.outsideArea();
cross_section_area_grout =
borehole_geometry.area() - _pipes.outer_pipe.outsideArea();
+
+ _thermal_resistances.fill(std::numeric_limits<double>::quiet_NaN());
}
static constexpr int number_of_unknowns = 3;
@@ -82,7 +84,8 @@ protected:
std::array<double, number_of_unknowns> _thermal_resistances;
/// Flow velocity inside the pipes and annulus. Depends on the flow_rate.
- double _flow_velocity_inner, _flow_velocity_annulus;
+ double _flow_velocity_inner = std::numeric_limits<double>::quiet_NaN(),
+ _flow_velocity_annulus = std::numeric_limits<double>::quiet_NaN();
};
} // namespace BHE
} // namespace HeatTransportBHE
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHECommonUType.h b/ProcessLib/HeatTransportBHE/BHE/BHECommonUType.h
new file mode 100644
index 0000000000000000000000000000000000000000..b4381a5016edc45633369bf03d5c0491b1b1980f
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/BHECommonUType.h
@@ -0,0 +1,44 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#pragma once
+
+#include <Eigen/Eigen>
+#include "BHECommon.h"
+#include "FlowAndTemperatureControl.h"
+#include "PipeConfigurationUType.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+class BHECommonUType : public BHECommon
+{
+public:
+ BHECommonUType(BoreholeGeometry const& borehole,
+ RefrigerantProperties const& refrigerant,
+ GroutParameters const& grout,
+ FlowAndTemperatureControl const& flowAndTemperatureControl,
+ PipeConfigurationUType const& pipes)
+ : BHECommon{borehole, refrigerant, grout, flowAndTemperatureControl},
+ _pipes(pipes)
+ {
+ }
+
+protected:
+ PipeConfigurationUType const _pipes;
+
+ /// Flow velocity inside the pipes. Depends on the flow_rate.
+ double _flow_velocity = std::numeric_limits<double>::quiet_NaN();
+};
+} // namespace BHE
+} // namespace HeatTransportBHE
+} // end of namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHETypes.h b/ProcessLib/HeatTransportBHE/BHE/BHETypes.h
index 2523cb5ce760cdfd47b1ce0ce32390265c4c2b12..76c7f9ba6858dc28a2c1e9b59aec2a65a74b9148 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHETypes.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHETypes.h
@@ -12,6 +12,7 @@
#include <variant>
#include "BHE_1U.h"
+#include "BHE_2U.h"
#include "BHE_CXA.h"
#include "BHE_CXC.h"
@@ -21,7 +22,7 @@ namespace HeatTransportBHE
{
namespace BHE
{
-using BHETypes = std::variant<BHE_1U, BHE_CXA, BHE_CXC>;
+using BHETypes = std::variant<BHE_1U, BHE_CXA, BHE_CXC, BHE_2U>;
} // end of namespace BHE
} // end of namespace HeatTransportBHE
} // end of namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
index af9806f5f67acb3d3cec8fa74b907f038eb55c32..a2341ea96ecc835a4dca904f367600ff9954a31f 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
@@ -21,32 +21,14 @@ namespace HeatTransportBHE
{
namespace BHE
{
-BHE_1U::BHE_1U(BoreholeGeometry const& borehole,
- RefrigerantProperties const& refrigerant,
- GroutParameters const& grout,
- FlowAndTemperatureControl const& flowAndTemperatureControl,
- PipeConfiguration1U const& pipes)
- : BHECommon{borehole, refrigerant, grout, flowAndTemperatureControl},
- _pipes(pipes)
-{
- // Initialize thermal resistances.
- auto values = visit(
- [&](auto const& control) {
- return control(refrigerant.reference_temperature,
- 0. /* initial time */);
- },
- flowAndTemperatureControl);
- updateHeatTransferCoefficients(values.flow_rate);
-}
-
std::array<double, BHE_1U::number_of_unknowns> BHE_1U::pipeHeatCapacities()
const
{
- double const& rho_r = refrigerant.density;
- double const& specific_heat_capacity = refrigerant.specific_heat_capacity;
- double const& rho_g = grout.rho_g;
- double const& porosity_g = grout.porosity_g;
- double const& heat_cap_g = grout.heat_cap_g;
+ double const rho_r = refrigerant.density;
+ double const specific_heat_capacity = refrigerant.specific_heat_capacity;
+ double const rho_g = grout.rho_g;
+ double const porosity_g = grout.porosity_g;
+ double const heat_cap_g = grout.heat_cap_g;
return {{/*i1*/ rho_r * specific_heat_capacity,
/*o1*/ rho_r * specific_heat_capacity,
@@ -57,12 +39,12 @@ std::array<double, BHE_1U::number_of_unknowns> BHE_1U::pipeHeatCapacities()
std::array<double, BHE_1U::number_of_unknowns> BHE_1U::pipeHeatConductions()
const
{
- double const& lambda_r = refrigerant.thermal_conductivity;
- double const& rho_r = refrigerant.density;
- double const& Cp_r = refrigerant.specific_heat_capacity;
- double const& alpha_L = _pipes.longitudinal_dispersion_length;
- double const& porosity_g = grout.porosity_g;
- double const& lambda_g = grout.lambda_g;
+ double const lambda_r = refrigerant.thermal_conductivity;
+ double const rho_r = refrigerant.density;
+ double const Cp_r = refrigerant.specific_heat_capacity;
+ double const alpha_L = _pipes.longitudinal_dispersion_length;
+ double const porosity_g = grout.porosity_g;
+ double const lambda_g = grout.lambda_g;
double const velocity_norm = std::abs(_flow_velocity) * std::sqrt(2);
@@ -133,34 +115,22 @@ std::pair<double, double> thermalResistancesGroutSoil(double chi,
return 1.0 / ((1.0 / R_gg) + (1.0 / (2.0 * R_gs)));
};
- int count = 0;
- while (constraint() < 0.0)
+ std::array<double, 3> const multiplier{chi * 0.66, chi * 0.5 * 0.66, 0.0};
+ for (double m_chi : multiplier)
{
- if (count == 0)
+ if (constraint() >= 0)
{
- chi *= 0.66;
- R_gs = compute_R_gs(chi, R_g);
- R_gg = compute_R_gg(chi, R_gs, R_ar, R_g);
- }
- if (count == 1)
- {
- chi *= 0.5;
- R_gs = compute_R_gs(chi, R_g);
- R_gg = compute_R_gg(chi, R_gs, R_ar, R_g);
- }
- if (count == 2)
- {
- chi = 0.0;
- R_gs = compute_R_gs(chi, R_g);
- R_gg = compute_R_gg(chi, R_gs, R_ar, R_g);
break;
}
DBUG(
"Warning! Correction procedure was applied due to negative thermal "
- "resistance! Correction step #%d.\n",
- count);
- count++;
+ "resistance! Chi = %f.\n",
+ m_chi);
+
+ R_gs = compute_R_gs(m_chi, R_g);
+ R_gg = compute_R_gg(m_chi, R_gs, R_ar, R_g);
}
+
return {R_gg, R_gs};
}
@@ -181,9 +151,9 @@ std::array<double, BHE_1U::number_of_unknowns> BHE_1U::calcThermalResistances(
{
constexpr double pi = boost::math::constants::pi<double>();
- double const& lambda_r = refrigerant.thermal_conductivity;
- double const& lambda_g = grout.lambda_g;
- double const& lambda_p = _pipes.inlet.wall_thermal_conductivity;
+ double const lambda_r = refrigerant.thermal_conductivity;
+ double const lambda_g = grout.lambda_g;
+ double const lambda_p = _pipes.inlet.wall_thermal_conductivity;
// thermal resistances due to advective flow of refrigerant in the _pipes
// Eq. 36 in Diersch_2011_CG
@@ -197,8 +167,8 @@ std::array<double, BHE_1U::number_of_unknowns> BHE_1U::calcThermalResistances(
(2.0 * pi * lambda_p);
// the average outer diameter of the _pipes
- double const& d0 = _pipes.outlet.diameter;
- double const& D = borehole_geometry.diameter;
+ double const d0 = _pipes.outlet.diameter;
+ double const D = borehole_geometry.diameter;
// Eq. 51
double const chi = std::log(std::sqrt(D * D + 2 * d0 * d0) / 2 / d0) /
std::log(D / std::sqrt(2) / d0);
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
index 5fbbdd4604553af8da5f420b6eef804b862b7182..c57ee7ac9d195ffceface7387e57c88468482658 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
@@ -15,8 +15,9 @@
#include "BaseLib/Error.h"
#include "BHECommon.h"
+#include "BHECommonUType.h"
#include "FlowAndTemperatureControl.h"
-#include "PipeConfiguration1U.h"
+#include "PipeConfigurationUType.h"
namespace ProcessLib
{
@@ -37,14 +38,28 @@ namespace BHE
* through the thermal resistance values, which are calculated specifically
* during the initialization of the class.
*/
-class BHE_1U final : public BHECommon
+class BHE_1U final : public BHECommonUType
{
public:
BHE_1U(BoreholeGeometry const& borehole,
RefrigerantProperties const& refrigerant,
GroutParameters const& grout,
FlowAndTemperatureControl const& flowAndTemperatureControl,
- PipeConfiguration1U const& pipes);
+ PipeConfigurationUType const& pipes)
+ : BHECommonUType{borehole, refrigerant, grout,
+ flowAndTemperatureControl, pipes}
+ {
+ _thermal_resistances.fill(std::numeric_limits<double>::quiet_NaN());
+
+ // Initialize thermal resistances.
+ auto values = visit(
+ [&](auto const& control) {
+ return control(refrigerant.reference_temperature,
+ 0. /* initial time */);
+ },
+ flowAndTemperatureControl);
+ updateHeatTransferCoefficients(values.flow_rate);
+ }
static constexpr int number_of_unknowns = 4;
static constexpr int number_of_grout_zones = 2;
@@ -136,10 +151,6 @@ public:
static constexpr std::pair<int, int> inflow_outflow_bc_component_ids[] = {
{0, 1}};
-private:
- // Placing it here before using it in the cross_section_areas.
- PipeConfiguration1U const _pipes;
-
public:
std::array<double, number_of_unknowns> crossSectionAreas() const
{
@@ -162,9 +173,6 @@ private:
/// 1) Diersch (2013) FEFLOW book on page 958, M.3, or
/// 2) Diersch (2011) Comp & Geosci 37:1122-1135, Eq. 90-97.
std::array<double, number_of_unknowns> _thermal_resistances;
-
- /// Flow velocity inside the pipes. Depends on the flow_rate.
- double _flow_velocity;
};
} // namespace BHE
} // namespace HeatTransportBHE
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..013a0fc731b8035dc8db098c664ce07c795d263d
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
@@ -0,0 +1,251 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "BHE_2U.h"
+
+#include <boost/math/constants/constants.hpp>
+#include "FlowAndTemperatureControl.h"
+#include "Physics.h"
+#include "ThermoMechanicalFlowProperties.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+std::array<double, BHE_2U::number_of_unknowns> BHE_2U::pipeHeatCapacities()
+ const
+{
+ double const rho_r = refrigerant.density;
+ double const specific_heat_capacity = refrigerant.specific_heat_capacity;
+ double const rho_g = grout.rho_g;
+ double const porosity_g = grout.porosity_g;
+ double const heat_cap_g = grout.heat_cap_g;
+
+ return {{/*i1*/ rho_r * specific_heat_capacity,
+ /*i2*/ rho_r * specific_heat_capacity,
+ /*o1*/ rho_r * specific_heat_capacity,
+ /*o2*/ rho_r * specific_heat_capacity,
+ /*g1*/ (1.0 - porosity_g) * rho_g * heat_cap_g,
+ /*g2*/ (1.0 - porosity_g) * rho_g * heat_cap_g,
+ /*g3*/ (1.0 - porosity_g) * rho_g * heat_cap_g,
+ /*g4*/ (1.0 - porosity_g) * rho_g * heat_cap_g}};
+}
+
+std::array<double, BHE_2U::number_of_unknowns> BHE_2U::pipeHeatConductions()
+ const
+{
+ double const lambda_r = refrigerant.thermal_conductivity;
+ double const rho_r = refrigerant.density;
+ double const Cp_r = refrigerant.specific_heat_capacity;
+ double const alpha_L = _pipes.longitudinal_dispersion_length;
+ double const porosity_g = grout.porosity_g;
+ double const lambda_g = grout.lambda_g;
+
+ // Here we calculate the laplace coefficients in the governing
+ // equations of BHE. These governing equations can be found in
+ // 1) Diersch (2013) FEFLOW book on page 952, M.120-122, or
+ // 2) Diersch (2011) Comp & Geosci 37:1122-1135, Eq. 19-22.
+ return {{// pipe i1
+ (lambda_r + rho_r * Cp_r * alpha_L * _flow_velocity),
+ // pipe i2
+ (lambda_r + rho_r * Cp_r * alpha_L * _flow_velocity),
+ // pipe o1
+ (lambda_r + rho_r * Cp_r * alpha_L * _flow_velocity),
+ // pipe o2
+ (lambda_r + rho_r * Cp_r * alpha_L * _flow_velocity),
+ // pipe g1
+ (1.0 - porosity_g) * lambda_g,
+ // pipe g2
+ (1.0 - porosity_g) * lambda_g,
+ // pipe g3
+ (1.0 - porosity_g) * lambda_g,
+ // pipe g4
+ (1.0 - porosity_g) * lambda_g}};
+}
+
+std::array<Eigen::Vector3d, BHE_2U::number_of_unknowns>
+BHE_2U::pipeAdvectionVectors() const
+{
+ double const rho_r = refrigerant.density;
+ double const Cp_r = refrigerant.specific_heat_capacity;
+
+ return {{// pipe i1
+ {0, 0, -rho_r * Cp_r * _flow_velocity},
+ // pipe i2
+ {0, 0, -rho_r * Cp_r * _flow_velocity},
+ // pipe o1
+ {0, 0, rho_r * Cp_r * _flow_velocity},
+ // pipe o2
+ {0, 0, rho_r * Cp_r * _flow_velocity},
+ // grout g1
+ {0, 0, 0},
+ // grout g2
+ {0, 0, 0},
+ // grout g3
+ {0, 0, 0},
+ // grout g4
+ {0, 0, 0}}};
+}
+
+double compute_R_gs_2U(double const chi, double const R_g)
+{
+ return (1 - chi) * R_g;
+}
+
+double compute_R_gg_2U(double const chi, double const R_gs, double const R_ar,
+ double const R_g)
+{
+ double const R_gg = 2.0 * R_gs * (R_ar - 2.0 * chi * R_g) /
+ (2.0 * R_gs - R_ar + 2.0 * chi * R_g);
+ if (!std::isfinite(R_gg))
+ {
+ OGS_FATAL(
+ "Error!!! Grout Thermal Resistance is an infinite number! The "
+ "simulation will be stopped!");
+ }
+
+ return R_gg;
+}
+
+/// Thermal resistances due to grout-soil exchange.
+///
+/// Check if constraints regarding negative thermal resistances are violated
+/// apply correction procedure.
+/// Section (1.5.5) in FEFLOW White Papers Vol V.
+std::vector<double> thermalResistancesGroutSoil2U(double chi,
+ double const R_ar_1,
+ double const R_ar_2,
+ double const R_g)
+{
+ double R_gs = compute_R_gs_2U(chi, R_g);
+ double R_gg_1 = compute_R_gg_2U(chi, R_gs, R_ar_1, R_g);
+ double R_gg_2 = compute_R_gg_2U(chi, R_gs, R_ar_2,
+ R_g); // Resulting thermal resistances.
+
+ auto constraint = [&]() {
+ return 1.0 / ((1.0 / R_gg_1) + (1.0 / (2.0 * R_gs)));
+ };
+
+ std::array<double, 3> const multiplier{chi * 0.66, chi * 0.5 * 0.66, 0.0};
+ for (double m_chi : multiplier)
+ {
+ if (constraint() >= 0)
+ {
+ break;
+ }
+ DBUG(
+ "Warning! Correction procedure was applied due to negative thermal "
+ "resistance! Chi = %f.\n",
+ m_chi);
+ R_gs = compute_R_gs_2U(m_chi, R_g);
+ R_gg_1 = compute_R_gg_2U(m_chi, R_gs, R_ar_1, R_g);
+ R_gg_2 = compute_R_gg_2U(m_chi, R_gs, R_ar_2, R_g);
+ }
+
+ return {R_gg_1, R_gg_2, R_gs};
+}
+
+void BHE_2U::updateHeatTransferCoefficients(double const flow_rate)
+
+{
+ auto const tm_flow_properties = calculateThermoMechanicalFlowPropertiesPipe(
+ _pipes.inlet, borehole_geometry.length, refrigerant, flow_rate);
+
+ _flow_velocity = tm_flow_properties.velocity;
+ _thermal_resistances =
+ calcThermalResistances(tm_flow_properties.nusselt_number);
+}
+
+/// Nu is the Nusselt number.
+std::array<double, BHE_2U::number_of_unknowns> BHE_2U::calcThermalResistances(
+ double const Nu)
+{
+ constexpr double pi = boost::math::constants::pi<double>();
+
+ double const lambda_r = refrigerant.thermal_conductivity;
+ double const lambda_g = grout.lambda_g;
+ double const lambda_p = _pipes.inlet.wall_thermal_conductivity;
+
+ // thermal resistances due to advective flow of refrigerant in the _pipes
+ // Eq. 36 in Diersch_2011_CG
+ double const R_adv_i = 1.0 / (Nu * lambda_r * pi);
+ double const R_adv_o = 1.0 / (Nu * lambda_r * pi);
+
+ // thermal resistance due to thermal conductivity of the pipe wall material
+ // Eq. 49
+ double const R_con_a =
+ std::log(_pipes.inlet.outsideDiameter() / _pipes.inlet.diameter) /
+ (2.0 * pi * lambda_p);
+
+ // the average outer diameter of the _pipes
+ double const d0 = _pipes.outlet.outsideDiameter();
+ double const D = borehole_geometry.diameter;
+ // Eq. 38
+ double const chi =
+ std::log(std::sqrt(D * D + 4 * d0 * d0) / 2 / std::sqrt(2) / d0) /
+ std::log(D / 2 / d0);
+ // Eq. 39
+ // thermal resistances of the grout
+ double const R_g =
+ std::acosh((D * D + d0 * d0 - 2 * _pipes.distance * _pipes.distance) /
+ (2 * D * d0)) /
+ (2 * pi * lambda_g) *
+ (3.098 - 4.432 * std::sqrt(2) * _pipes.distance / D +
+ 2.364 * 2 * _pipes.distance * _pipes.distance / D / D);
+ // thermal resistance due to the grout transition.
+ double const R_con_b = chi * R_g;
+
+ // Eq. 29 and 30
+ double const R_fig = 2 * R_adv_i + 2 * R_con_a + R_con_b;
+ double const R_fog = 2 * R_adv_o + 2 * R_con_a + R_con_b;
+
+ // thermal resistance due to inter-grout exchange
+ double const R_ar_1 =
+ std::acosh((2.0 * _pipes.distance * _pipes.distance - d0 * d0) / d0 /
+ d0) /
+ (2.0 * pi * lambda_g);
+
+ double const R_ar_2 =
+ std::acosh((2.0 * 2.0 * _pipes.distance * _pipes.distance - d0 * d0) /
+ d0 / d0) /
+ (2.0 * pi * lambda_g);
+
+ std::vector<double> _intergrout_thermal_exchange;
+ _intergrout_thermal_exchange =
+ thermalResistancesGroutSoil2U(chi, R_ar_1, R_ar_2, R_g);
+
+ return {{R_fig, R_fog, _intergrout_thermal_exchange[0],
+ _intergrout_thermal_exchange[1], _intergrout_thermal_exchange[2]}};
+
+ // keep the following lines------------------------------------------------
+ // when debugging the code, printing the R and phi values are needed--------
+ // std::cout << "Rfig =" << R_fig << " Rfog =" << R_fog << " Rgg =" <<
+ // R_gg << " Rgs =" << R_gs << "\n"; double phi_fig = 1.0 / (R_fig *
+ // S_i); double phi_fog = 1.0 / (R_fog * S_o); double phi_gg = 1.0 / (R_gg
+ // * S_g1); double phi_gs = 1.0 / (R_gs * S_gs); std::cout << "phi_fig ="
+ // << phi_fig << " phi_fog =" << phi_fog << " phi_gg =" << phi_gg << "
+ // phi_gs =" << phi_gs << "\n";
+ // -------------------------------------------------------------------------
+}
+
+double BHE_2U::updateFlowRateAndTemperature(double const T_out,
+ double const current_time)
+{
+ auto values =
+ visit([&](auto const& control) { return control(T_out, current_time); },
+ flowAndTemperatureControl);
+ updateHeatTransferCoefficients(values.flow_rate);
+ return values.temperature;
+}
+} // namespace BHE
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h
new file mode 100644
index 0000000000000000000000000000000000000000..b1c2c4a97c07bd4875f5c7ce9267f0af077b4091
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h
@@ -0,0 +1,248 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <Eigen/Eigen>
+
+#include "BaseLib/Error.h"
+
+#include "BHECommon.h"
+#include "BHECommonUType.h"
+#include "FlowAndTemperatureControl.h"
+#include "PipeConfigurationUType.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+/**
+ * The BHE_2U class is the realization of 2U type of Borehole Heate Exchanger.
+ * In this class, the pipe heat capacity, pipe heat conductiion, pipe advection
+ * vectors are intialized according to the geometry of 2U type of BHE.
+ * For 2U type of BHE, 8 primary unknowns are assigned on the 1D BHE elements.
+ * They are the temperature in inflow pipe T_i1 and T_i2, temperature in outflow
+ * pipe T_o1 and T_o2, temperature of the four grout zones surrounding the
+ * inflow and outflow pipe T_g1, T_g2, T_g3 and T_g4. These primary varaibles
+ * are solved according to heat convection and conduction equations on the pipes
+ * and also in the grout zones. The interaction of the 2U type of BHE and the
+ * surrounding soil is regulated through the thermal resistance values, which
+ * are calculated specifically during the initialization of the class.
+ */
+class BHE_2U final : public BHECommonUType
+{
+public:
+ BHE_2U(BoreholeGeometry const& borehole,
+ RefrigerantProperties const& refrigerant,
+ GroutParameters const& grout,
+ FlowAndTemperatureControl const& flowAndTemperatureControl,
+ PipeConfigurationUType const& pipes)
+ : BHECommonUType{borehole, refrigerant, grout,
+ flowAndTemperatureControl, pipes}
+ {
+ _thermal_resistances.fill(std::numeric_limits<double>::quiet_NaN());
+
+ // Initialize thermal resistances.
+ auto values = visit(
+ [&](auto const& control) {
+ return control(refrigerant.reference_temperature,
+ 0. /* initial time */);
+ },
+ flowAndTemperatureControl);
+ updateHeatTransferCoefficients(values.flow_rate);
+ }
+
+ static constexpr int number_of_unknowns = 8;
+ static constexpr int number_of_grout_zones = 4;
+
+ std::array<double, number_of_unknowns> pipeHeatCapacities() const;
+
+ std::array<double, number_of_unknowns> pipeHeatConductions() const;
+
+ std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors()
+ const;
+
+ template <int NPoints,
+ typename SingleUnknownMatrixType,
+ typename RMatrixType,
+ typename RPiSMatrixType,
+ typename RSMatrixType>
+ void assembleRMatrices(
+ int const idx_bhe_unknowns,
+ Eigen::MatrixBase<SingleUnknownMatrixType> const& matBHE_loc_R,
+ Eigen::MatrixBase<RMatrixType>& R_matrix,
+ Eigen::MatrixBase<RPiSMatrixType>& R_pi_s_matrix,
+ Eigen::MatrixBase<RSMatrixType>& R_s_matrix) const
+ {
+ switch (idx_bhe_unknowns)
+ {
+ case 0: // PHI_fig
+ R_matrix.block(0, 4 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(4 * NPoints, 0, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R; // R_i1
+
+ R_matrix.block(2, 5 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(5 * NPoints, 2, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R; // R_i2
+
+ R_matrix.block(0, 0, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R; // K_i1
+ R_matrix.block(NPoints, NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R; // K_i2
+ R_matrix.block(4 * NPoints, 4 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R;
+ R_matrix.block(5 * NPoints, 5 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R; // K_ig
+ return;
+ case 1: // PHI_fog
+ R_matrix.block(2 * NPoints, 6 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(6 * NPoints, 2 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R; // R_o1
+ R_matrix.block(3 * NPoints, 7 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(7 * NPoints, 3 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R; // R_o2
+
+ R_matrix.block(2 * NPoints, 2 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R; // K_o1
+ R_matrix.block(3 * NPoints, 3 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R; // K_o2
+ R_matrix.block(6 * NPoints, 6 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R;
+ R_matrix.block(7 * NPoints, 7 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R; // K_og
+ return;
+ case 2: // PHI_gg_1
+ R_matrix.block(4 * NPoints, 6 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(6 * NPoints, 4 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(4 * NPoints, 7 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(7 * NPoints, 4 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(5 * NPoints, 6 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(6 * NPoints, 5 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(5 * NPoints, 7 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(7 * NPoints, 5 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R; // R_g1
+
+ R_matrix.block(4 * NPoints, 4 * NPoints, NPoints, NPoints) +=
+ 2.0 * matBHE_loc_R;
+ R_matrix.block(5 * NPoints, 5 * NPoints, NPoints, NPoints) +=
+ 2.0 * matBHE_loc_R; // K_ig // notice we have two
+ // PHI_gg_1 term here.
+ R_matrix.block(6 * NPoints, 6 * NPoints, NPoints, NPoints) +=
+ 2.0 * matBHE_loc_R;
+ R_matrix.block(7 * NPoints, 7 * NPoints, NPoints, NPoints) +=
+ 2.0 * matBHE_loc_R; // K_og // see Diersch 2013 FEFLOW
+ // book page 954 Table M.2
+ return;
+ case 3: // PHI_gg_2
+ R_matrix.block(4 * NPoints, 5 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(5 * NPoints, 4 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(6 * NPoints, 7 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(7 * NPoints, 6 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R; // R_g2
+
+ R_matrix.block(4 * NPoints, 4 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R;
+ R_matrix.block(5 * NPoints, 5 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R; // K_ig // notice we only have
+ // 1 PHI_gg term here.
+ R_matrix.block(6 * NPoints, 6 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R;
+ R_matrix.block(7 * NPoints, 7 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R; // K_og // see Diersch 2013 FEFLOW
+ // book page 954 Table M.2
+ return;
+ case 4: // PHI_gs
+ R_s_matrix.template block<NPoints, NPoints>(0, 0).noalias() +=
+ 1.0 * matBHE_loc_R;
+
+ R_pi_s_matrix.block(4 * NPoints, 0, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_pi_s_matrix.block(5 * NPoints, 0, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_pi_s_matrix.block(6 * NPoints, 0, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_pi_s_matrix.block(7 * NPoints, 0, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(4 * NPoints, 4 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R;
+ R_matrix.block(5 * NPoints, 5 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R; // K_ig
+ R_matrix.block(6 * NPoints, 6 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R;
+ R_matrix.block(7 * NPoints, 7 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R; // K_og
+ return;
+ default:
+ OGS_FATAL(
+ "Error!!! In the function BHE_2U::assembleRMatrices, "
+ "the index of bhe unknowns is out of range! ");
+ }
+ }
+
+ /// Return the inflow temperature for the boundary condition.
+ double updateFlowRateAndTemperature(double T_out, double current_time);
+
+ double thermalResistance(int const unknown_index) const
+ {
+ return _thermal_resistances[unknown_index];
+ }
+
+ static constexpr std::pair<int, int> inflow_outflow_bc_component_ids[] = {
+ {0, 2}, {1, 3}};
+
+public:
+ std::array<double, number_of_unknowns> crossSectionAreas() const
+ {
+ return {{
+ _pipes.inlet.area(),
+ _pipes.inlet.area(),
+ _pipes.outlet.area(),
+ _pipes.outlet.area(),
+ borehole_geometry.area() / 4 - _pipes.inlet.outsideArea(),
+ borehole_geometry.area() / 4 - _pipes.inlet.outsideArea(),
+ borehole_geometry.area() / 4 - _pipes.outlet.outsideArea(),
+ borehole_geometry.area() / 4 - _pipes.outlet.outsideArea(),
+ }};
+ }
+
+ void updateHeatTransferCoefficients(double const flow_rate);
+
+private:
+ std::array<double, number_of_unknowns> calcThermalResistances(
+ double const Nu);
+
+private:
+ /// PHI_fig, PHI_fog, PHI_gg, PHI_gs;
+ /// Here we store the thermal resistances needed for computation of the heat
+ /// exchange coefficients in the governing equations of BHE.
+ /// These governing equations can be found in
+ /// 1) Diersch (2013) FEFLOW book on page 958, M.3, or
+ /// 2) Diersch (2011) Comp & Geosci 37:1122-1135, Eq. 90-97.
+ std::array<double, number_of_unknowns> _thermal_resistances;
+};
+} // namespace BHE
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.cpp b/ProcessLib/HeatTransportBHE/BHE/CreateBHEUType.cpp
similarity index 70%
rename from ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.cpp
rename to ProcessLib/HeatTransportBHE/BHE/CreateBHEUType.cpp
index 83f822dee9530c21ccd87c1215d4002cd8d0fa2d..71af576c6ab4786b561edbdfaec57a937eb7461c 100644
--- a/ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/CreateBHEUType.cpp
@@ -8,8 +8,11 @@
* http://www.opengeosys.org/project/license
*/
-#include "CreateBHE1U.h"
+#include "CreateBHEUType.h"
#include "BaseLib/ConfigTree.h"
+
+#include "BHE_1U.h"
+#include "BHE_2U.h"
#include "CreateFlowAndTemperatureControl.h"
namespace ProcessLib
{
@@ -17,7 +20,12 @@ namespace HeatTransportBHE
{
namespace BHE
{
-BHE::BHE_1U createBHE1U(
+static std::tuple<BoreholeGeometry,
+ RefrigerantProperties,
+ GroutParameters,
+ FlowAndTemperatureControl,
+ PipeConfigurationUType>
+parseBHEUTypeConfig(
BaseLib::ConfigTree const& config,
std::map<std::string,
std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
@@ -41,8 +49,8 @@ BHE::BHE_1U createBHE1U(
//! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__pipes__longitudinal_dispersion_length}
pipes_config.getConfigParameter<double>(
"longitudinal_dispersion_length");
- PipeConfiguration1U const pipes{inlet_pipe, outlet_pipe, pipe_distance,
- pipe_longitudinal_dispersion_length};
+ PipeConfigurationUType const pipes{inlet_pipe, outlet_pipe, pipe_distance,
+ pipe_longitudinal_dispersion_length};
//! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__grout}
auto const grout = createGroutParameters(config.getConfigSubtree("grout"));
@@ -60,6 +68,30 @@ BHE::BHE_1U createBHE1U(
return {borehole_geometry, refrigerant, grout, flowAndTemperatureControl,
pipes};
}
+
+template <typename T_BHE>
+T_BHE createBHEUType(
+ BaseLib::ConfigTree const& config,
+ std::map<std::string,
+ std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
+ curves)
+{
+ auto UType = parseBHEUTypeConfig(config, curves);
+ return {std::get<0>(UType), std::get<1>(UType), std::get<2>(UType),
+ std::get<3>(UType), std::get<4>(UType)};
+}
+
+template BHE_1U createBHEUType<BHE_1U>(
+ BaseLib::ConfigTree const& config,
+ std::map<std::string,
+ std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
+ curves);
+
+template BHE_2U createBHEUType<BHE_2U>(
+ BaseLib::ConfigTree const& config,
+ std::map<std::string,
+ std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
+ curves);
} // namespace BHE
} // namespace HeatTransportBHE
} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.h b/ProcessLib/HeatTransportBHE/BHE/CreateBHEUType.h
similarity index 79%
rename from ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.h
rename to ProcessLib/HeatTransportBHE/BHE/CreateBHEUType.h
index 304826268a29b259bcfa27be8e31545684f4abb3..df3af29a213a0c8dfa3611388939f9453a191bbb 100644
--- a/ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.h
+++ b/ProcessLib/HeatTransportBHE/BHE/CreateBHEUType.h
@@ -10,7 +10,10 @@
#pragma once
-#include "BHE_1U.h"
+#include <Eigen/Eigen>
+#include "BHECommon.h"
+#include "FlowAndTemperatureControl.h"
+#include "PipeConfigurationUType.h"
namespace BaseLib
{
@@ -22,7 +25,8 @@ namespace HeatTransportBHE
{
namespace BHE
{
-BHE::BHE_1U createBHE1U(
+template <typename T_BHE>
+T_BHE createBHEUType(
BaseLib::ConfigTree const& config,
std::map<std::string,
std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
diff --git a/ProcessLib/HeatTransportBHE/BHE/PipeConfiguration1U.h b/ProcessLib/HeatTransportBHE/BHE/PipeConfigurationUType.h
similarity index 90%
rename from ProcessLib/HeatTransportBHE/BHE/PipeConfiguration1U.h
rename to ProcessLib/HeatTransportBHE/BHE/PipeConfigurationUType.h
index dde87c8f6e33222c1b61ecf58854986be59a37ca..87f8016ee4ede37ba61c652ad3e471c6c624a624 100644
--- a/ProcessLib/HeatTransportBHE/BHE/PipeConfiguration1U.h
+++ b/ProcessLib/HeatTransportBHE/BHE/PipeConfigurationUType.h
@@ -9,7 +9,7 @@
*/
#pragma once
-
+#include "BaseLib/ConfigTree.h"
#include "Pipe.h"
namespace ProcessLib
@@ -18,7 +18,7 @@ namespace HeatTransportBHE
{
namespace BHE
{
-struct PipeConfiguration1U
+struct PipeConfigurationUType
{
Pipe const inlet;
Pipe const outlet;
diff --git a/ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.cpp b/ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.cpp
index 642103b0a1f1a7a540c46ad01ccc333f873bdc24..0072e00d280581e6304ee919f9d150323163e1b2 100644
--- a/ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.cpp
+++ b/ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.cpp
@@ -16,8 +16,8 @@
#include "ProcessLib/Output/CreateSecondaryVariables.h"
#include "BHE/BHETypes.h"
-#include "BHE/CreateBHE1U.h"
#include "BHE/CreateBHECoaxial.h"
+#include "BHE/CreateBHEUType.h"
#include "HeatTransportBHEProcess.h"
#include "HeatTransportBHEProcessData.h"
@@ -186,7 +186,8 @@ std::unique_ptr<Process> createHeatTransportBHEProcess(
if (bhe_type == "1U")
{
- bhes.emplace_back(BHE::createBHE1U(bhe_config, curves));
+ bhes.emplace_back(
+ BHE::createBHEUType<BHE::BHE_1U>(bhe_config, curves));
continue;
}
@@ -203,6 +204,13 @@ std::unique_ptr<Process> createHeatTransportBHEProcess(
BHE::createBHECoaxial<BHE::BHE_CXC>(bhe_config, curves));
continue;
}
+
+ if (bhe_type == "2U")
+ {
+ bhes.emplace_back(
+ BHE::createBHEUType<BHE::BHE_2U>(bhe_config, curves));
+ continue;
+ }
OGS_FATAL("Unknown BHE type '%s'.", bhe_type.c_str());
}
// end of reading BHE parameters -------------------------------------------
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
index 512549391c20ad24b348c942b73e444bdd4cc647..d6d5089a19cd4b724161f57d888199cbedf4964f 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
@@ -11,7 +11,6 @@
#pragma once
#include "HeatTransportBHELocalAssemblerBHE.h"
-
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
#include "ProcessLib/Utils/InitShapeMatrices.h"
@@ -61,6 +60,7 @@ HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod, BHEType>::
_R_matrix.setZero(bhe_unknowns_size, bhe_unknowns_size);
_R_pi_s_matrix.setZero(bhe_unknowns_size, temperature_size);
_R_s_matrix.setZero(temperature_size, temperature_size);
+ static constexpr int max_num_thermal_exchange_terms = 5;
// formulate the local BHE R matrix
for (int idx_bhe_unknowns = 0; idx_bhe_unknowns < bhe_unknowns;
idx_bhe_unknowns++)
@@ -84,13 +84,24 @@ HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod, BHEType>::
// The following assembly action is according to Diersch (2013) FEFLOW
// book please refer to M.127 and M.128 on page 955 and 956
- _bhe.template assembleRMatrices<ShapeFunction::NPOINTS>(
- idx_bhe_unknowns, matBHE_loc_R, _R_matrix, _R_pi_s_matrix,
- _R_s_matrix);
+ // The if check is absolutely necessary because
+ // (i) In the CXA and CXC case, there are 3 exchange terms,
+ // and it is the same as the number of unknowns;
+ // (ii) In the 1U case, there are 4 exchange terms,
+ // and it is again same as the number of unknowns;
+ // (iii) In the 2U case, there are 5 exchange terms,
+ // and it is less than the number of unknowns (8).
+ if (idx_bhe_unknowns < max_num_thermal_exchange_terms)
+ {
+ _bhe.template assembleRMatrices<ShapeFunction::NPOINTS>(
+ idx_bhe_unknowns, matBHE_loc_R, _R_matrix, _R_pi_s_matrix,
+ _R_s_matrix);
+ }
} // end of loop over BHE unknowns
// debugging
- // std::string sep = "\n----------------------------------------\n";
+ // std::string sep =
+ // "\n----------------------------------------\n";
// Eigen::IOFormat CleanFmt(4, 0, ", ", "\n", "[", "]");
// std::cout << "_R_matrix: \n" << sep;
// std::cout << _R_matrix.format(CleanFmt) << sep;
@@ -189,12 +200,10 @@ void HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod,
.noalias() += _bhe.number_of_grout_zones * _R_s_matrix;
// debugging
- // std::string sep =
- // "\n----------------------------------------\n";
+ // std::string sep = "\n----------------------------------------\n";
// Eigen::IOFormat CleanFmt(4, 0, ", ", "\n", "[", "]");
// std::cout << local_K.format(CleanFmt) << sep;
// std::cout << local_M.format(CleanFmt) << sep;
}
-
} // namespace HeatTransportBHE
} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h
index f0402efa58923cee28a564a841c51b161e50bfd2..f40923aadb0b1023b7db8c1dae893c9aafddd4d5 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h
@@ -276,6 +276,13 @@ private:
e, std::get<BHE::BHE_CXC>(bhe),
std::forward<ConstructorArgs>(args)...}};
}
+
+ if (std::holds_alternative<BHE::BHE_2U>(bhe))
+ {
+ return LADataIntfPtr{new LADataBHE<ShapeFunction, BHE::BHE_2U>{
+ e, std::get<BHE::BHE_2U>(bhe),
+ std::forward<ConstructorArgs>(args)...}};
+ }
OGS_FATAL(
"Trying to create local assembler for an unknown BHE type.");
};
diff --git a/ProcessLib/HeatTransportBHE/Tests.cmake b/ProcessLib/HeatTransportBHE/Tests.cmake
index 2dfdb9361c260858d8d68f96af4ee17cd5b19caf..93f9a91523a07d9708688b5d4239aa3ad41b420a 100644
--- a/ProcessLib/HeatTransportBHE/Tests.cmake
+++ b/ProcessLib/HeatTransportBHE/Tests.cmake
@@ -51,3 +51,16 @@ AddTest(
3D_deep_BHE_CXC_pcs_0_ts_10_t_600.000000.vtu 3D_deep_BHE_CXC_pcs_0_ts_10_t_600.000000.vtu temperature_BHE1 temperature_BHE1 0 5e-15
3D_deep_BHE_CXC_pcs_0_ts_10_t_600.000000.vtu 3D_deep_BHE_CXC_pcs_0_ts_10_t_600.000000.vtu temperature_soil temperature_soil 0 1e-13
)
+
+AddTest(
+ NAME HeatTransportBHE_3D_2U_BHE
+ PATH Parabolic/T/3D_2U_BHE
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS 3D_2U_BHE.prj
+ WRAPPER time
+ TESTER vtkdiff
+ REQUIREMENTS NOT OGS_USE_MPI
+ DIFF_DATA
+ 3D_2U_BHE_pcs_0_ts_10_t_600.000000.vtu 3D_2U_BHE_pcs_0_ts_10_t_600.000000.vtu temperature_BHE1 temperature_BHE1 1e-12 1e-14
+ 3D_2U_BHE_pcs_0_ts_10_t_600.000000.vtu 3D_2U_BHE_pcs_0_ts_10_t_600.000000.vtu temperature_soil temperature_soil 1e-12 1e-13
+)
diff --git a/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE.gml b/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE.gml
new file mode 100644
index 0000000000000000000000000000000000000000..af1d206400ee077409b0bb37dc3dc8b68ec9d4a3
--- /dev/null
+++ b/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE.gml
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:fbd0fbbdff9aab09706e7cf2a33a9f32b4b75c2136ebc28106825802ba686ccd
+size 641
diff --git a/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE.prj b/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE.prj
new file mode 100644
index 0000000000000000000000000000000000000000..6791bf41659ff7d55a41f85f94d7ddaab0b7f003
--- /dev/null
+++ b/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE.prj
@@ -0,0 +1,220 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <mesh>3D_2U_BHE.vtu</mesh>
+ <geometry>3D_2U_BHE.gml</geometry>
+ <processes>
+ <process>
+ <name>HeatTransportBHE</name>
+ <type>HEAT_TRANSPORT_BHE</type>
+ <integration_order>2</integration_order>
+ <process_variables>
+ <process_variable>temperature_soil</process_variable>
+ <process_variable>temperature_BHE1</process_variable>
+ </process_variables>
+ <thermal_conductivity_solid>K_s</thermal_conductivity_solid>
+ <heat_capacity_solid>Cp_s</heat_capacity_solid>
+ <density_solid>rho_s</density_solid>
+ <thermal_conductivity_fluid>K_f</thermal_conductivity_fluid>
+ <heat_capacity_fluid>Cp_f</heat_capacity_fluid>
+ <density_fluid>rho_f</density_fluid>
+ <thermal_conductivity_gas>K_g</thermal_conductivity_gas>
+ <heat_capacity_gas>Cp_g</heat_capacity_gas>
+ <density_gas>rho_g</density_gas>
+ <borehole_heat_exchangers>
+ <borehole_heat_exchanger>
+ <type>2U</type>
+ <flow_and_temperature_control>
+ <type>TemperatureCurveConstantFlow</type>
+ <flow_rate>2.0e-4</flow_rate>
+ <temperature_curve>inflow_temperature</temperature_curve>
+ </flow_and_temperature_control>
+ <borehole>
+ <length>18.0</length>
+ <diameter>0.13</diameter>
+ </borehole>
+ <grout>
+ <density>2190.0</density>
+ <porosity>0.0</porosity>
+ <heat_capacity>1141.55</heat_capacity>
+ <thermal_conductivity>1</thermal_conductivity>
+ </grout>
+ <pipes>
+ <inlet>
+ <diameter> 0.0378</diameter>
+ <wall_thickness>0.0029</wall_thickness>
+ <wall_thermal_conductivity>0.42</wall_thermal_conductivity>
+ </inlet>
+ <outlet>
+ <diameter>0.0378</diameter>
+ <wall_thickness>0.0029</wall_thickness>
+ <wall_thermal_conductivity>0.42</wall_thermal_conductivity>
+ </outlet>
+ <distance_between_pipes>0.053</distance_between_pipes>
+ <longitudinal_dispersion_length>0.001</longitudinal_dispersion_length>
+ </pipes>
+ <refrigerant>
+ <density>1052</density>
+ <viscosity>0.0003</viscosity>
+ <specific_heat_capacity>3802.28</specific_heat_capacity>
+ <thermal_conductivity>0.48</thermal_conductivity>
+ <reference_temperature>22</reference_temperature>
+ </refrigerant>
+ </borehole_heat_exchanger>
+ </borehole_heat_exchangers>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="HeatTransportBHE">
+ <nonlinear_solver>basic_picard</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <reltol>1e-6</reltol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <time_stepping>
+ <type>FixedTimeStepping</type>
+ <t_initial> 0.0 </t_initial>
+ <!-- use the following for full simulation
+ <t_end> 186420 </t_end>
+ -->
+ <t_end> 600 </t_end>
+ <timesteps>
+ <!-- use the following for full simulation
+ <pair><repeat>3107</repeat><delta_t>60</delta_t></pair>
+ -->
+ <pair>
+ <repeat>10</repeat>
+ <delta_t>60</delta_t>
+ </pair>
+ </timesteps>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>3D_2U_BHE</prefix>
+ <timesteps>
+ <!-- use the following for full simulation
+ <pair><repeat> 10</repeat><each_steps> 1 </each_steps></pair>
+ <pair><repeat> 310</repeat><each_steps> 10 </each_steps></pair>
+ -->
+ <pair>
+ <repeat> 100</repeat>
+ <each_steps> 1 </each_steps>
+ </pair>
+ </timesteps>
+ <variables>
+ <variable>temperature_soil</variable>
+ <variable>temperature_BHE1</variable>
+ </variables>
+ </output>
+ </time_loop>
+ <parameters>
+ <parameter>
+ <name>K_s</name>
+ <type>Constant</type>
+ <value>2.78018</value>
+ </parameter>
+ <parameter>
+ <name>Cp_s</name>
+ <type>Constant</type>
+ <value>1778</value>
+ </parameter>
+ <parameter>
+ <name>rho_s</name>
+ <type>Constant</type>
+ <value>1800</value>
+ </parameter>
+ <parameter>
+ <name>K_f</name>
+ <type>Constant</type>
+ <value>0.1</value>
+ </parameter>
+ <parameter>
+ <name>Cp_f</name>
+ <type>Constant</type>
+ <value>4068</value>
+ </parameter>
+ <parameter>
+ <name>rho_f</name>
+ <type>Constant</type>
+ <value>992.92</value>
+ </parameter>
+ <parameter>
+ <name>K_g</name>
+ <type>Constant</type>
+ <value>3.2</value>
+ </parameter>
+ <parameter>
+ <name>Cp_g</name>
+ <type>Constant</type>
+ <value>1000</value>
+ </parameter>
+ <parameter>
+ <name>rho_g</name>
+ <type>Constant</type>
+ <value>2500</value>
+ </parameter>
+ <parameter>
+ <name>T0</name>
+ <type>Constant</type>
+ <value>22</value>
+ </parameter>
+ <parameter>
+ <name>T0_BHE1</name>
+ <type>Constant</type>
+ <values>22 22 22 22 22 22 22 22</values>
+ </parameter>
+ </parameters>
+ <process_variables>
+ <process_variable>
+ <name>temperature_soil</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>T0</initial_condition>
+ <boundary_conditions>
+ </boundary_conditions>
+ </process_variable>
+ <process_variable>
+ <name>temperature_BHE1</name>
+ <components>8</components>
+ <order>1</order>
+ <initial_condition>T0_BHE1</initial_condition>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_picard</name>
+ <type>Picard</type>
+ <max_iter>100</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+ <eigen>
+ <solver_type>BiCGSTAB</solver_type>
+ <precon_type>ILUT</precon_type>
+ <max_iteration_step>1000</max_iteration_step>
+ <error_tolerance>1e-16</error_tolerance>
+ </eigen>
+ <petsc>
+ <prefix>gw</prefix>
+ <parameters>-gw_ksp_type cg -gw_pc_type bjacobi -gw_ksp_rtol 1e-16 -gw_ksp_max_it 10000</parameters>
+ </petsc>
+ </linear_solver>
+ </linear_solvers>
+ <curves>
+ <curve>
+ <name>inflow_temperature</name>
+ <coords>0 60 120 180 240 300 360 420 480 540 600 1200 1800 2400 3000 3600</coords>
+ <values>22 22.2 22.4 22.6 22.8 23 23.2 23.4 23.6 23.8 24 26 28 30 32 34</values>
+ </curve>
+ </curves>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE.vtu b/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..6be1b976a9c13e9459614724c796a27b83cc8470
--- /dev/null
+++ b/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ccf17b79ecd727f0f7019717b24ee8467ebd5cbc6fbe1d974b65769e1280a915
+size 683856
diff --git a/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE_pcs_0_ts_10_t_600.000000.vtu b/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE_pcs_0_ts_10_t_600.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..04de3f6c5ad2f08d030b34d285dc16707944e44a
--- /dev/null
+++ b/Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE_pcs_0_ts_10_t_600.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:0a1fcb70c3ab3609b505e86cff2f00e2e3b5fe53b08aa1d031b4e6b4fa52c03b
+size 263370
diff --git a/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE.pandoc b/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE.pandoc
new file mode 100644
index 0000000000000000000000000000000000000000..7b7a2483be0bce9779f0ba12f3c9398eeb437517
--- /dev/null
+++ b/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE.pandoc
@@ -0,0 +1,51 @@
++++
+author = "Chaofan Chen, Haibing Shao"
+date = "2018-09-25T13:44:00+01:00"
+title = "3D 2U BHE"
+weight = 123
+project = "Parabolic/T/3D_2U_BHE/3D_2U_BHE.prj"
+
+[menu]
+ [menu.benchmarks]
+ parent = "heat-transport-bhe"
+
++++
+
+{{< data-link >}}
+
+## Model Setup
+
+The numerical model was established based on dual continuum method developed by Diersch et al. (2011). In the 3D 2U type BHE benchmark, the depth of the whole domain is 18.5 m with a square cross section of 5 m per side. The BHE bottom is situated at -18 m and its total length is 18 m. Detailed parameters for the model configuration are listed in the follwoing table.
+
+| Parameter | Symbol | Value | Unit |
+| ---------------------------------- |:------------------ | -------------------:| --------------------------: |
+| Soil thermal conductivity | $\lambda_{s}$ | 2.78018 | $\mathrm{W m^{-1} K^{-1}}$ |
+| Soil heat capacity | $(\rho c)_{soil}$ | $3.2\times10^{6}$ | $\mathrm{Jm^{-3}K^{-1}}$ |
+| Length of the BHE | $L$ | 18 | $\mathrm{m}$ |
+| Diameter of the BHE | $D$ | 0.13 | $\mathrm{m}$ |
+| Diameter of the pipeline | $d$ | 3.78 | $\mathrm{cm}$ |
+| Wall thickness of the pipeline | $b$ | 0.29 | $\mathrm{cm}$ |
+| Distance between pipelines | $w$ | 5.3 | $\mathrm{cm}$ |
+| Pipeline wall thermal conductivity | $\lambda_{p}$ | 0.42 | $\mathrm{W m^{-1} K^{-1}}$ |
+| Grout thermal conductivity | $\lambda_{g}$ | 1 | $\mathrm{W m^{-1} K^{-1}}$ |
+| Grout heat capacity | $(\rho c)_{grout}$ | $2.5\times10^{6}$ | $\mathrm{Jm^{-3}K^{-1}}$ |
+
+## OGS Input Files
+
+The detailed input parameters can be seen from the 3D_2U_BHE.prj file. The inflow temperature of the BHE, which was imposed as boundary condition of the BHE is shown in Figure 1. All the initial temperatures are set as 22 $^{\circ}$C.
+
+{{< img src="../3D_2U_BHE_figures/In_out_temperature_comparison.png" width="200">}}
+
+Figure 1: Inflow temperature curve and outflow temperature comparison
+
+## Results
+
+The OGS numerical outflow temperature over time was compared against results of the FEFLOW software as shown in the Figure 1. And the vertical distributed temperature of circulating water was presented in Figure 2 after operation for 3000 s. The comparison figures demonstrate that the OGS numerical results and FEFLOW results can match very well and the biggest absolute error of outflow temperature is 0.19 $^{\circ}$C at the starting up stage, while such error decreases to 0.05 $^{\circ}$C after 3600 s' operation. The maximum relative error of vertical temperature is 0.015 \% after operation for 3000 s.
+
+{{< img src="../3D_2U_BHE_figures/vertical_temperature_distribution.png" width="200">}}
+
+Figure 2: Comparison of vertical temperature distribution
+
+## References
+
+[1] Diersch, H. J., Bauer, D., Heidemann, W., Rühaak, W., & Schätzl, P. (2011). Finite element modeling of borehole heat exchanger systems: Part 1. Fundamentals. Computers & Geosciences, 37(8), 1122-1135.
diff --git a/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/In_out_temperature_comparison.png b/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/In_out_temperature_comparison.png
new file mode 100644
index 0000000000000000000000000000000000000000..11d9e029025ff5516e75ce141ce33e14b70e41d0
--- /dev/null
+++ b/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/In_out_temperature_comparison.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:20eb7531c6e2e0f870909087e6fc6738030d252b451a5557ef464e4de869c7b6
+size 73683
diff --git a/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution.png b/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution.png
new file mode 100644
index 0000000000000000000000000000000000000000..5e176dd82f6980de7190341b65ca360ee98f5102
--- /dev/null
+++ b/web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f4ced51d7c7f7cbfa480bb3c6b28107f26fb386d9f19704bd2f9f21d5e632e49
+size 71358