Skip to content
Snippets Groups Projects
Commit f635a679 authored by Dmitri Naumov's avatar Dmitri Naumov
Browse files

Merge branch 'BHE' into 'master' 

[BHE] fixed the thermal resistance errors in 1U and 2U type BHE

See merge request ogs/ogs!3426
parents eb4f1f93 ec8181a1
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 339 additions and 66 deletions
ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
+ 15
− 13
View file @ f635a679
......@@ -124,19 +124,20 @@ double compute_R_gg(double const chi, double const R_gs, double const R_ar,
/// Check if constraints regarding negative thermal resistances are violated
/// apply correction procedure.
/// Section (1.5.5) in FEFLOW White Papers Vol V.
std::pair<double, double> thermalResistancesGroutSoil(double chi,
double const R_ar,
double const R_g)
std::tuple<double, double, double> thermalResistancesGroutSoil(
double const chi, double const R_ar, double const R_g)
{
double R_gs = compute_R_gs(chi, R_g);
double R_gg =
compute_R_gg(chi, R_gs, R_ar, R_g); // Resulting thermal resistances.
double new_chi = chi;
auto constraint = [&]() {
return 1.0 / ((1.0 / R_gg) + (1.0 / (2.0 * R_gs)));
};
std::array<double, 3> const multiplier{chi * 0.66, chi * 0.5 * 0.66, 0.0};
std::array<double, 3> const multiplier{chi * 2.0 / 3.0, chi * 1.0 / 3.0,
0.0};
for (double m_chi : multiplier)
{
if (constraint() >= 0)
......@@ -150,9 +151,10 @@ std::pair<double, double> thermalResistancesGroutSoil(double chi,
R_gs = compute_R_gs(m_chi, R_g);
R_gg = compute_R_gg(m_chi, R_gs, R_ar, R_g);
new_chi = m_chi;
}
return {R_gg, R_gs};
return {new_chi, R_gg, R_gs};
}
void BHE_1U::updateHeatTransferCoefficients(double const flow_rate)
......@@ -199,11 +201,6 @@ std::array<double, BHE_1U::number_of_unknowns> BHE_1U::calcThermalResistances(
std::acosh((D * D + d0 * d0 - _pipes.distance * _pipes.distance) /
(2 * D * d0)) /
(2 * pi * lambda_g) * (1.601 - 0.888 * _pipes.distance / D);
// thermal resistance due to the grout transition.
double const R_con_b = chi * R_g;
// Eq. 29 and 30
double const R_fig = R_adv_i1 + R_con_a + R_con_b;
double const R_fog = R_adv_o1 + R_con_a + R_con_b;
// thermal resistance due to inter-grout exchange
double const R_ar =
......@@ -211,9 +208,14 @@ std::array<double, BHE_1U::number_of_unknowns> BHE_1U::calcThermalResistances(
d0) /
(2.0 * pi * lambda_g);
double R_gg;
double R_gs;
std::tie(R_gg, R_gs) = thermalResistancesGroutSoil(chi, R_ar, R_g);
auto const [chi_new, R_gg, R_gs] =
thermalResistancesGroutSoil(chi, R_ar, R_g);
// thermal resistance due to the grout transition.
double const R_con_b = chi_new * R_g;
// Eq. 29 and 30
double const R_fig = R_adv_i1 + R_con_a + R_con_b;
double const R_fog = R_adv_o1 + R_con_a + R_con_b;
return {{R_fig, R_fog, R_gg, R_gs}};
......
ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
+ 19
− 16
View file @ f635a679
......@@ -142,21 +142,24 @@ double compute_R_gg_2U(double const chi, double const R_gs, double const R_ar,
/// 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)
std::tuple<double, double, double, double> thermalResistancesGroutSoil2U(
double const 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.
double chi_new = chi;
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};
std::array<double, 3> const multiplier{chi * 2.0 / 3.0, chi * 1.0 / 3.0,
0.0};
for (double m_chi : multiplier)
{
if (constraint() >= 0)
......@@ -170,9 +173,10 @@ std::vector<double> thermalResistancesGroutSoil2U(double 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);
chi_new = m_chi;
}
return {R_gg_1, R_gg_2, R_gs};
return {chi_new, R_gg_1, R_gg_2, R_gs};
}
void BHE_2U::updateHeatTransferCoefficients(double const flow_rate)
......@@ -222,12 +226,6 @@ std::array<double, BHE_2U::number_of_unknowns> BHE_2U::calcThermalResistances(
(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 =
......@@ -240,12 +238,17 @@ std::array<double, BHE_2U::number_of_unknowns> BHE_2U::calcThermalResistances(
d0 / d0) /
(2.0 * pi * lambda_g);
std::vector<double> _intergrout_thermal_exchange;
_intergrout_thermal_exchange =
auto const [chi_new, R_gg_1, R_gg_2, R_gs] =
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]}};
// thermal resistance due to the grout transition.
double const R_con_b = chi_new * R_g;
// Eq. 29 and 30
double const R_fig = R_adv_i + R_con_a + R_con_b;
double const R_fog = R_adv_o + R_con_a + R_con_b;
return {{R_fig, R_fog, R_gg_1, R_gg_2, R_gs}};
// keep the following lines------------------------------------------------
// when debugging the code, printing the R and phi values are needed--------
......
Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE.prj
+ 3
− 0
View file @ f635a679
<?xml version="1.0" encoding="ISO-8859-1"?>
<OpenGeoSysProject>
<mesh>3D_2U_BHE.vtu</mesh>
<!-- use the following mesh for benchmark with FEFLOW model
<mesh>3D_2U_BHE_benchmark.vtu</mesh>
-->
<geometry>3D_2U_BHE.gml</geometry>
<processes>
<process>
......
Tests/Data/Parabolic/T/3D_2U_BHE/3D_2U_BHE_powerBC.prj 0 → 100644
+ 231
− 0
View file @ f635a679
<?xml version="1.0" encoding="ISO-8859-1"?>
<OpenGeoSysProject>
<mesh>3D_2U_BHE_benchmark.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>
<borehole_heat_exchangers>
<borehole_heat_exchanger>
<type>2U</type>
<flow_and_temperature_control>
<!-- if power type boundary conditions are selected, the related
power value should be set with half of the presumed BHE power-->
<type>PowerCurveConstantFlow</type>
<power_curve>heating_load</power_curve>
<flow_rate>2.0e-4</flow_rate>
</flow_and_temperature_control>
<borehole>
<length>18.0</length>
<diameter>0.13</diameter>
</borehole>
<grout>
<density>2190.0</density>
<porosity>0.0</porosity>
<specific_heat_capacity>1141.55</specific_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>
<media>
<medium id="0">
<phases>
<phase>
<type>AqueousLiquid</type>
<properties>
<property>
<name>phase_velocity</name>
<type>Constant</type>
<value>0 0 0</value>
</property>
<property>
<name>specific_heat_capacity</name>
<type>Constant</type>
<value>4068</value>
</property>
<property>
<name>density</name>
<type>Constant</type>
<value>992.92</value>
</property>
</properties>
</phase>
<phase>
<type>Solid</type>
<properties>
<property>
<name>specific_heat_capacity</name>
<type>Constant</type>
<value>1778</value>
</property>
<property>
<name>density</name>
<type>Constant</type>
<value>1800</value>
</property>
</properties>
</phase>
<phase>
<type>Gas</type>
<properties>
<property>
<name>specific_heat_capacity</name>
<type>Constant</type>
<value>1000</value>
</property>
<property>
<name>density</name>
<type>Constant</type>
<value>2500</value>
</property>
</properties>
</phase>
</phases>
<properties>
<property>
<name>porosity</name>
<type>Constant</type>
<value>0</value>
</property>
<property>
<name>thermal_conductivity</name>
<type>Constant</type>
<value>2.78018</value>
</property>
<property>
<name>thermal_longitudinal_dispersivity</name>
<type>Constant</type>
<value>0</value>
</property>
<property>
<name>thermal_transversal_dispersivity</name>
<type>Constant</type>
<value>0</value>
</property>
</properties>
</medium>
</media>
<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>
<t_end> 3600 </t_end>
<timesteps>
<pair><repeat>60</repeat><delta_t>60</delta_t></pair>
</timesteps>
</time_stepping>
</process>
</processes>
<output>
<type>VTK</type>
<prefix>3D_2U_BHE_powerBC</prefix>
<timesteps>
<pair><repeat> 1</repeat><each_steps> 1 </each_steps></pair>
</timesteps>
<variables>
<variable>temperature_soil</variable>
<variable>temperature_BHE1</variable>
</variables>
<suffix>_ts_{:timestep}_t_{:time}</suffix>
</output>
</time_loop>
<parameters>
<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>heating_load</name>
<coords>0 3600</coords>
<values>-157.86 -157.86</values>
</curve>
</curves>
</OpenGeoSysProject>
web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE.md
+ 37
− 3
View file @ f635a679
+++
author = "Chaofan Chen, Haibing Shao"
author = "Chaofan Chen, Shuang Chen, Haibing Shao"
date = "2018-09-25T13:44:00+01:00"
title = "3D 2U BHE"
weight = 123
......@@ -32,20 +32,54 @@ The numerical model was established based on dual continuum method developed by
## OGS Input Files
For this benchmark, Two different scenarios were carried out by applying two different boundary conditions imposed on the BHE.
* Fixed inflow boundary condition
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. The flow rate within each U-pipe is set to $2.0\times10^{-4}$ $\mathrm{m^{3} s^{-1}}$ during the whole simulation time.
{{< img src="../3D_2U_BHE_figures/In_out_temperature_comparison.png" width="200">}}
Figure 1: Inflow temperature curve and outflow temperature comparison
* Fixed power boundary condition
The detailed input parameters can be seen from the 3D_2U_BHE_powerBC.prj file.
A -315.72 $W$ thermal load is imposed on the BHE to extract the heat from the subsurface during the entire simulation.
All the other parameters adopted in the model is same as the ones used in the scenario with fixed inflow boundary condition.
## FEFLOW Input Files
For the benchmark a FEFLOW model is presented.
The mesh used in the OGS model is directly coverted from the FEFLOW model mesh, to ensure that there is no influence to the comparison results from the mesh side.
Both the FEFLOW and ogs model mesh can be found in the ogs GitLab (<https://gitlab.opengeosys.org/ogs/ogs/-/merge_requests/3426>).
## FEFLOW Input Files
## 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.
The computed resutls from scenario by adopting the fixed inflow boundary condition are illustracted in Figure 1 and Figure 2.
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 3300 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.20 $^{\circ}$C after 360 s' operation, while such error decreases to 0.037 $^{\circ}$C after 3600 s' operation. The maximum relative error of vertical temperature is 0.019 \% after operation for 3300 s.
{{< img src="../3D_2U_BHE_figures/vertical_temperature_distribution.png" width="200">}}
Figure 2: Comparison of vertical temperature distribution
Figure 2: Comparison of vertical temperature distribution from scenario by adopting the fixed inflow boundary condition
Figure 3 shows the the vertical distributed temperature of circulating fluid after operation for 3300 s by adopting the power boundary condition in OGS and FEFLOW models.
A 0.18 $^{\circ}$C difference is found between the averaged vertical temperature from the two models.
The reason to the results difference is due to the different power boundary condition type adopted in the two software.
In FEFLOW the power boundary condtion is based on the outlet temperature calculated from the last time step (non-iterative).
Compared to it, the default power boundary condition adopted in the OGS `Heat_Transport_BHE` process is based on the outlet temperature calculated from the current time step (with-iterative).
Besides, by setting python bindings, the current OGS `Heat_Transport_BHE` process is capable to adopt the power boundary condition type used in the FEFLOW software.
In this way, the computed vertical distributed circulating fluid temperature difference between the OGS and FEFLOW models is becoming much closer to each other.
{{< img src="../3D_2U_BHE_figures/vertical_temperature_distribution_powerBC.png" width="200">}}
Figure 3: Comparison of vertical temperature distribution from scenario by adopting the power boundary condition
## 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.
[2] FEFLOW online documentation. URL: <http://www.feflow.info/html/help72/feflow/14_References/GUI/Dialogs/bhe_editor.html>.
\ No newline at end of file
web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/In_out_temperature_comparison.png
+ 0
− 0
View replaced file @ eb4f1f93 View file @ f635a679
web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/In_out_temperature_comparison.png

72 KiB | W: | H:

web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/In_out_temperature_comparison.png

138 KiB | W: | H:

web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/In_out_temperature_comparison.png
web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/In_out_temperature_comparison.png
web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/In_out_temperature_comparison.png
web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/In_out_temperature_comparison.png
  • 2-up
  • Swipe
  • Onion skin
web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution.png
+ 0
− 0
View replaced file @ eb4f1f93 View file @ f635a679
web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution.png

69.7 KiB | W: | H:

web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution.png

129 KiB | W: | H:

web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution.png
web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution.png
web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution.png
web/content/docs/benchmarks/heat-transport-bhe/3D_2U_BHE_figures/vertical_temperature_distribution.png
  • 2-up
  • Swipe
  • Onion skin
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment