Skip to content
Snippets Groups Projects
Unverified Commit 427777a3 authored by Dmitri Naumov's avatar Dmitri Naumov Committed by GitHub
Browse files

Merge pull request #2711 from boyanmeng/BHE_GWF

Groundwater Advection Feature for Heat_Transport_BHE Process
parents d9f6b97a 4ce38960
No related branches found
No related tags found
No related merge requests found
Showing
with 576 additions and 94 deletions
...@@ -14,8 +14,10 @@ ...@@ -14,8 +14,10 @@
#include <vector> #include <vector>
#include "MaterialLib/MPL/Medium.h" #include "MaterialLib/MPL/Medium.h"
#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
#include "MathLib/LinAlg/Eigen/EigenMapTools.h" #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h" #include "NumLib/Fem/ShapeMatrixPolicy.h"
#include "NumLib/Function/Interpolation.h"
#include "ProcessLib/HeatTransportBHE/HeatTransportBHEProcessData.h" #include "ProcessLib/HeatTransportBHE/HeatTransportBHEProcessData.h"
#include "ProcessLib/Utils/InitShapeMatrices.h" #include "ProcessLib/Utils/InitShapeMatrices.h"
...@@ -43,10 +45,8 @@ HeatTransportBHELocalAssemblerSoil<ShapeFunction, IntegrationMethod>:: ...@@ -43,10 +45,8 @@ HeatTransportBHELocalAssemblerSoil<ShapeFunction, IntegrationMethod>::
_ip_data.reserve(n_integration_points); _ip_data.reserve(n_integration_points);
_secondary_data.N.resize(n_integration_points); _secondary_data.N.resize(n_integration_points);
_shape_matrices = initShapeMatrices<ShapeFunction, _shape_matrices = initShapeMatrices<ShapeFunction, ShapeMatricesType,
ShapeMatricesType, IntegrationMethod, 3 /* GlobalDim */>(
IntegrationMethod,
3 /* GlobalDim */>(
e, is_axially_symmetric, _integration_method); e, is_axially_symmetric, _integration_method);
ParameterLib::SpatialPosition x_position; ParameterLib::SpatialPosition x_position;
...@@ -61,8 +61,7 @@ HeatTransportBHELocalAssemblerSoil<ShapeFunction, IntegrationMethod>:: ...@@ -61,8 +61,7 @@ HeatTransportBHELocalAssemblerSoil<ShapeFunction, IntegrationMethod>::
auto const& sm = _shape_matrices[ip]; auto const& sm = _shape_matrices[ip];
double const w = _integration_method.getWeightedPoint(ip).getWeight() * double const w = _integration_method.getWeightedPoint(ip).getWeight() *
sm.integralMeasure * sm.detJ; sm.integralMeasure * sm.detJ;
_ip_data.push_back( _ip_data.push_back({sm.N, sm.dNdx, w});
{sm.N.transpose() * sm.N * w, sm.dNdx.transpose() * sm.dNdx * w});
_secondary_data.N[ip] = sm.N; _secondary_data.N[ip] = sm.N;
} }
...@@ -89,6 +88,7 @@ void HeatTransportBHELocalAssemblerSoil<ShapeFunction, IntegrationMethod>:: ...@@ -89,6 +88,7 @@ void HeatTransportBHELocalAssemblerSoil<ShapeFunction, IntegrationMethod>::
auto const& medium = *_process_data.media_map->getMedium(_element_id); auto const& medium = *_process_data.media_map->getMedium(_element_id);
auto const& solid_phase = medium.phase("Solid"); auto const& solid_phase = medium.phase("Solid");
auto const& liquid_phase = medium.phase("AqueousLiquid");
MaterialPropertyLib::VariableArray vars; MaterialPropertyLib::VariableArray vars;
...@@ -99,13 +99,19 @@ void HeatTransportBHELocalAssemblerSoil<ShapeFunction, IntegrationMethod>:: ...@@ -99,13 +99,19 @@ void HeatTransportBHELocalAssemblerSoil<ShapeFunction, IntegrationMethod>::
{ {
pos.setIntegrationPoint(ip); pos.setIntegrationPoint(ip);
auto& ip_data = _ip_data[ip]; auto& ip_data = _ip_data[ip];
auto const& M_ip = ip_data.NTN_product_times_w; auto const& N = ip_data.N;
auto const& K_ip = ip_data.dNdxTdNdx_product_times_w; auto const& dNdx = ip_data.dNdx;
auto const& w = ip_data.integration_weight;
auto const k_s = double T_int_pt = 0.0;
solid_phase NumLib::shapeFunctionInterpolate(local_x, N, T_int_pt);
.property(
MaterialPropertyLib::PropertyType::thermal_conductivity) vars[static_cast<int>(MaterialPropertyLib::Variable::temperature)] =
T_int_pt;
// for now only using the solid and liquid phase parameters
auto const density_s =
solid_phase.property(MaterialPropertyLib::PropertyType::density)
.template value<double>(vars, pos, t); .template value<double>(vars, pos, t);
auto const heat_capacity_s = auto const heat_capacity_s =
...@@ -114,16 +120,73 @@ void HeatTransportBHELocalAssemblerSoil<ShapeFunction, IntegrationMethod>:: ...@@ -114,16 +120,73 @@ void HeatTransportBHELocalAssemblerSoil<ShapeFunction, IntegrationMethod>::
MaterialPropertyLib::PropertyType::specific_heat_capacity) MaterialPropertyLib::PropertyType::specific_heat_capacity)
.template value<double>(vars, pos, t); .template value<double>(vars, pos, t);
auto const density_s = auto const density_f =
solid_phase.property(MaterialPropertyLib::PropertyType::density) liquid_phase.property(MaterialPropertyLib::PropertyType::density)
.template value<double>(vars, pos, t); .template value<double>(vars, pos, t);
// for now only using the solid phase parameters
auto const heat_capacity_f =
liquid_phase
.property(
MaterialPropertyLib::PropertyType::specific_heat_capacity)
.template value<double>(vars, pos, t);
auto const porosity =
medium.property(MaterialPropertyLib::PropertyType::porosity)
.template value<double>(vars, pos, t);
auto const velocity_arr =
liquid_phase
.property(MaterialPropertyLib::PropertyType::phase_velocity)
.value(vars, pos, t);
auto const velocity = Eigen::Map<Eigen::Vector3d const>(
std::get<MaterialPropertyLib::Vector>(velocity_arr).data(), 1, 3);
// calculate the hydrodynamic thermodispersion tensor
auto const thermal_conductivity =
MaterialPropertyLib::formEigenTensor<3>(
medium
.property(
MaterialPropertyLib::PropertyType::thermal_conductivity)
.value(vars, pos, t));
auto thermal_conductivity_dispersivity = thermal_conductivity;
double const velocity_magnitude = velocity.norm();
if (velocity_magnitude >= std::numeric_limits<double>::epsilon())
{
auto const thermal_dispersivity_longitudinal =
medium
.property(MaterialPropertyLib::PropertyType::
thermal_longitudinal_dispersivity)
.template value<double>();
auto const thermal_dispersivity_transversal =
medium
.property(MaterialPropertyLib::PropertyType::
thermal_transversal_dispersivity)
.template value<double>();
auto const thermal_dispersivity =
density_f * heat_capacity_f *
(thermal_dispersivity_transversal * velocity_magnitude *
Eigen::Matrix3d::Identity() +
(thermal_dispersivity_longitudinal -
thermal_dispersivity_transversal) /
velocity_magnitude * velocity * velocity.transpose());
thermal_conductivity_dispersivity += thermal_dispersivity;
}
// assemble Conductance matrix // assemble Conductance matrix
local_K.noalias() += K_ip * k_s; local_K.noalias() +=
(dNdx.transpose() * thermal_conductivity_dispersivity * dNdx +
N.transpose() * velocity.transpose() * dNdx * density_f *
heat_capacity_f) *
w;
// assemble Mass matrix // assemble Mass matrix
local_M.noalias() += M_ip * density_s * heat_capacity_s; local_M.noalias() += N.transpose() * N * w *
(density_s * heat_capacity_s * (1 - porosity) +
density_f * heat_capacity_f * porosity);
} }
// debugging // debugging
......
...@@ -36,7 +36,11 @@ public: ...@@ -36,7 +36,11 @@ public:
ShapeMatrixPolicyType<ShapeFunction, 3 /* GlobalDim */>; ShapeMatrixPolicyType<ShapeFunction, 3 /* GlobalDim */>;
using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType; using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
using NodalVectorType = typename ShapeMatricesType::NodalVectorType; using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
using NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType;
using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices; using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
using GlobalDimNodalMatrixType =
typename ShapeMatricesType::GlobalDimNodalMatrixType;
HeatTransportBHELocalAssemblerSoil( HeatTransportBHELocalAssemblerSoil(
HeatTransportBHELocalAssemblerSoil const&) = delete; HeatTransportBHELocalAssemblerSoil const&) = delete;
...@@ -68,8 +72,9 @@ private: ...@@ -68,8 +72,9 @@ private:
HeatTransportBHEProcessData& _process_data; HeatTransportBHEProcessData& _process_data;
std::vector< std::vector<
IntegrationPointDataSoil<NodalMatrixType>, IntegrationPointDataSoil<NodalRowVectorType, GlobalDimNodalMatrixType>,
Eigen::aligned_allocator<IntegrationPointDataSoil<NodalMatrixType>>> Eigen::aligned_allocator<IntegrationPointDataSoil<
NodalRowVectorType, GlobalDimNodalMatrixType>>>
_ip_data; _ip_data;
IntegrationMethod const _integration_method; IntegrationMethod const _integration_method;
......
...@@ -14,11 +14,12 @@ namespace ProcessLib ...@@ -14,11 +14,12 @@ namespace ProcessLib
{ {
namespace HeatTransportBHE namespace HeatTransportBHE
{ {
template <typename NodalMatrixType> template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType>
struct IntegrationPointDataSoil final struct IntegrationPointDataSoil final
{ {
NodalMatrixType const NTN_product_times_w; NodalRowVectorType const N;
NodalMatrixType const dNdxTdNdx_product_times_w; GlobalDimNodalMatrixType const dNdx;
double const integration_weight;
EIGEN_MAKE_ALIGNED_OPERATOR_NEW; EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
}; };
......
...@@ -64,3 +64,17 @@ AddTest( ...@@ -64,3 +64,17 @@ AddTest(
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_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 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
) )
AddTest(
NAME HeatTransportBHE_3D_BHE_groundwater_advection
PATH Parabolic/T/3D_BHE_GW_advection
RUNTIME 30
EXECUTABLE ogs
EXECUTABLE_ARGS BHE_GW_advection.prj
WRAPPER time
TESTER vtkdiff
REQUIREMENTS NOT OGS_USE_MPI
DIFF_DATA
BHE_GW_advection_pcs_0_ts_10_t_500.000000.vtu BHE_GW_advection_pcs_0_ts_10_t_500.000000.vtu temperature_BHE1 temperature_BHE1 1e-12 1e-14
BHE_GW_advection_pcs_0_ts_10_t_500.000000.vtu BHE_GW_advection_pcs_0_ts_10_t_500.000000.vtu temperature_soil temperature_soil 1e-12 1e-13
)
...@@ -61,14 +61,14 @@ ...@@ -61,14 +61,14 @@
<type>AqueousLiquid</type> <type>AqueousLiquid</type>
<properties> <properties>
<property> <property>
<name>specific_heat_capacity</name> <name>phase_velocity</name>
<type>Constant</type> <type>Constant</type>
<value>4068</value> <value>0 0 0</value>
</property> </property>
<property> <property>
<name>thermal_conductivity</name> <name>specific_heat_capacity</name>
<type>Constant</type> <type>Constant</type>
<value>0.1</value> <value>4068</value>
</property> </property>
<property> <property>
<name>density</name> <name>density</name>
...@@ -85,11 +85,6 @@ ...@@ -85,11 +85,6 @@
<type>Constant</type> <type>Constant</type>
<value>1778</value> <value>1778</value>
</property> </property>
<property>
<name>thermal_conductivity</name>
<type>Constant</type>
<value>2.78018</value>
</property>
<property> <property>
<name>density</name> <name>density</name>
<type>Constant</type> <type>Constant</type>
...@@ -105,11 +100,6 @@ ...@@ -105,11 +100,6 @@
<type>Constant</type> <type>Constant</type>
<value>1000</value> <value>1000</value>
</property> </property>
<property>
<name>thermal_conductivity</name>
<type>Constant</type>
<value>3.2</value>
</property>
<property> <property>
<name>density</name> <name>density</name>
<type>Constant</type> <type>Constant</type>
...@@ -119,6 +109,16 @@ ...@@ -119,6 +109,16 @@
</phase> </phase>
</phases> </phases>
<properties> <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> <property>
<name>thermal_longitudinal_dispersivity</name> <name>thermal_longitudinal_dispersivity</name>
<type>Constant</type> <type>Constant</type>
......
<?xml version="1.0" encoding="ISO-8859-1"?>
<OpenGeoSysProject>
<meshes>
<mesh>BHE_GW_advection.vtu</mesh>
<mesh>front.vtu</mesh>
<mesh>back.vtu</mesh>
<mesh>left.vtu</mesh>
<mesh>right.vtu</mesh>
</meshes>
<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>1U</type>
<flow_and_temperature_control>
<type>FixedPowerConstantFlow</type>
<flow_rate>2.0e-4</flow_rate>
<power>300</power>
</flow_and_temperature_control>
<borehole>
<length>15.0</length>
<diameter>0.13</diameter>
</borehole>
<grout>
<density>2190.0</density>
<porosity>0.0</porosity>
<heat_capacity>1735.160</heat_capacity>
<thermal_conductivity>0.73</thermal_conductivity>
</grout>
<pipes>
<inlet>
<diameter> 0.013665</diameter>
<wall_thickness>0.003035</wall_thickness>
<wall_thermal_conductivity>0.39</wall_thermal_conductivity>
</inlet>
<outlet>
<diameter>0.013665</diameter>
<wall_thickness>0.003035</wall_thickness>
<wall_thermal_conductivity>0.39</wall_thermal_conductivity>
</outlet>
<distance_between_pipes>0.053</distance_between_pipes>
<longitudinal_dispersion_length>0.001</longitudinal_dispersion_length>
</pipes>
<refrigerant>
<density>992.92</density>
<viscosity>0.00067418</viscosity>
<specific_heat_capacity>4068</specific_heat_capacity>
<thermal_conductivity>0.62863</thermal_conductivity>
<reference_temperature>25</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 1e-7 0</value>
</property>
<property>
<name>specific_heat_capacity</name>
<type>Constant</type>
<value>4200</value>
</property>
<property>
<name>density</name>
<type>Constant</type>
<value>1000</value>
</property>
</properties>
</phase>
<phase>
<type>Solid</type>
<properties>
<property>
<name>specific_heat_capacity</name>
<type>Constant</type>
<value>880</value>
</property>
<property>
<name>density</name>
<type>Constant</type>
<value>2650</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.26</value>
</property>
<property>
<name>thermal_conductivity</name>
<type>Constant</type>
<value>2.5</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> 500 </t_end>
<timesteps>
<pair>
<repeat>10</repeat>
<delta_t>50</delta_t>
</pair>
</timesteps>
</time_stepping>
</process>
</processes>
<output>
<type>VTK</type>
<prefix>BHE_GW_advection</prefix>
<timesteps>
<pair>
<repeat> 1</repeat>
<each_steps> 10 </each_steps>
</pair>
</timesteps>
<variables>
<variable>temperature_soil</variable>
<variable>temperature_BHE1</variable>
</variables>
</output>
</time_loop>
<parameters>
<parameter>
<name>T0</name>
<type>Constant</type>
<value>25</value>
</parameter>
<parameter>
<name>T_DBC</name>
<type>Constant</type>
<value>25</value>
</parameter>
<parameter>
<name>T0_BHE1</name>
<type>Constant</type>
<values>25.36 25.13 25.23 25.115</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_condition>
<mesh>front</mesh>
<type>Dirichlet</type>
<parameter>T_DBC</parameter>
</boundary_condition>
<boundary_condition>
<mesh>back</mesh>
<type>Dirichlet</type>
<parameter>T_DBC</parameter>
</boundary_condition>
<boundary_condition>
<mesh>left</mesh>
<type>Dirichlet</type>
<parameter>T_DBC</parameter>
</boundary_condition>
<boundary_condition>
<mesh>right</mesh>
<type>Dirichlet</type>
<parameter>T_DBC</parameter>
</boundary_condition>
</boundary_conditions>
</process_variable>
<process_variable>
<name>temperature_BHE1</name>
<components>4</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>50</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>
</OpenGeoSysProject>
File added
File added
File added
File added
File added
File added
...@@ -61,14 +61,14 @@ ...@@ -61,14 +61,14 @@
<type>AqueousLiquid</type> <type>AqueousLiquid</type>
<properties> <properties>
<property> <property>
<name>specific_heat_capacity</name> <name>phase_velocity</name>
<type>Constant</type> <type>Constant</type>
<value>4068</value> <value>0 0 0</value>
</property> </property>
<property> <property>
<name>thermal_conductivity</name> <name>specific_heat_capacity</name>
<type>Constant</type> <type>Constant</type>
<value>0.1</value> <value>4068</value>
</property> </property>
<property> <property>
<name>density</name> <name>density</name>
...@@ -85,11 +85,6 @@ ...@@ -85,11 +85,6 @@
<type>Constant</type> <type>Constant</type>
<value>1778</value> <value>1778</value>
</property> </property>
<property>
<name>thermal_conductivity</name>
<type>Constant</type>
<value>2.78018</value>
</property>
<property> <property>
<name>density</name> <name>density</name>
<type>Constant</type> <type>Constant</type>
...@@ -105,11 +100,6 @@ ...@@ -105,11 +100,6 @@
<type>Constant</type> <type>Constant</type>
<value>1000</value> <value>1000</value>
</property> </property>
<property>
<name>thermal_conductivity</name>
<type>Constant</type>
<value>3.2</value>
</property>
<property> <property>
<name>density</name> <name>density</name>
<type>Constant</type> <type>Constant</type>
...@@ -119,6 +109,16 @@ ...@@ -119,6 +109,16 @@
</phase> </phase>
</phases> </phases>
<properties> <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> <property>
<name>thermal_longitudinal_dispersivity</name> <name>thermal_longitudinal_dispersivity</name>
<type>Constant</type> <type>Constant</type>
......
...@@ -54,21 +54,21 @@ ...@@ -54,21 +54,21 @@
</borehole_heat_exchangers> </borehole_heat_exchangers>
</process> </process>
</processes> </processes>
<media> <media>
<medium id="0"> <medium id="0">
<phases> <phases>
<phase> <phase>
<type>AqueousLiquid</type> <type>AqueousLiquid</type>
<properties> <properties>
<property> <property>
<name>specific_heat_capacity</name> <name>phase_velocity</name>
<type>Constant</type> <type>Constant</type>
<value>4068</value> <value>0 0 0</value>
</property> </property>
<property> <property>
<name>thermal_conductivity</name> <name>specific_heat_capacity</name>
<type>Constant</type> <type>Constant</type>
<value>0.1</value> <value>4068</value>
</property> </property>
<property> <property>
<name>density</name> <name>density</name>
...@@ -85,11 +85,6 @@ ...@@ -85,11 +85,6 @@
<type>Constant</type> <type>Constant</type>
<value>1778</value> <value>1778</value>
</property> </property>
<property>
<name>thermal_conductivity</name>
<type>Constant</type>
<value>2.78018</value>
</property>
<property> <property>
<name>density</name> <name>density</name>
<type>Constant</type> <type>Constant</type>
...@@ -97,7 +92,7 @@ ...@@ -97,7 +92,7 @@
</property> </property>
</properties> </properties>
</phase> </phase>
<phase> <phase>
<type>Gas</type> <type>Gas</type>
<properties> <properties>
<property> <property>
...@@ -105,11 +100,6 @@ ...@@ -105,11 +100,6 @@
<type>Constant</type> <type>Constant</type>
<value>1000</value> <value>1000</value>
</property> </property>
<property>
<name>thermal_conductivity</name>
<type>Constant</type>
<value>3.2</value>
</property>
<property> <property>
<name>density</name> <name>density</name>
<type>Constant</type> <type>Constant</type>
...@@ -119,6 +109,16 @@ ...@@ -119,6 +109,16 @@
</phase> </phase>
</phases> </phases>
<properties> <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> <property>
<name>thermal_longitudinal_dispersivity</name> <name>thermal_longitudinal_dispersivity</name>
<type>Constant</type> <type>Constant</type>
......
...@@ -60,14 +60,14 @@ ...@@ -60,14 +60,14 @@
<type>AqueousLiquid</type> <type>AqueousLiquid</type>
<properties> <properties>
<property> <property>
<name>specific_heat_capacity</name> <name>phase_velocity</name>
<type>Constant</type> <type>Constant</type>
<value>4068</value> <value>0 0 0</value>
</property> </property>
<property> <property>
<name>thermal_conductivity</name> <name>specific_heat_capacity</name>
<type>Constant</type> <type>Constant</type>
<value>0.1</value> <value>4068</value>
</property> </property>
<property> <property>
<name>density</name> <name>density</name>
...@@ -84,11 +84,6 @@ ...@@ -84,11 +84,6 @@
<type>Constant</type> <type>Constant</type>
<value>2000</value> <value>2000</value>
</property> </property>
<property>
<name>thermal_conductivity</name>
<type>Constant</type>
<value>2.5</value>
</property>
<property> <property>
<name>density</name> <name>density</name>
<type>Constant</type> <type>Constant</type>
...@@ -104,11 +99,6 @@ ...@@ -104,11 +99,6 @@
<type>Constant</type> <type>Constant</type>
<value>1000</value> <value>1000</value>
</property> </property>
<property>
<name>thermal_conductivity</name>
<type>Constant</type>
<value>3.2</value>
</property>
<property> <property>
<name>density</name> <name>density</name>
<type>Constant</type> <type>Constant</type>
...@@ -118,6 +108,16 @@ ...@@ -118,6 +108,16 @@
</phase> </phase>
</phases> </phases>
<properties> <properties>
<property>
<name>porosity</name>
<type>Constant</type>
<value>0</value>
</property>
<property>
<name>thermal_conductivity</name>
<type>Constant</type>
<value>2.5</value>
</property>
<property> <property>
<name>thermal_longitudinal_dispersivity</name> <name>thermal_longitudinal_dispersivity</name>
<type>Constant</type> <type>Constant</type>
......
...@@ -60,14 +60,14 @@ ...@@ -60,14 +60,14 @@
<type>AqueousLiquid</type> <type>AqueousLiquid</type>
<properties> <properties>
<property> <property>
<name>specific_heat_capacity</name> <name>phase_velocity</name>
<type>Constant</type> <type>Constant</type>
<value>4068</value> <value>0 0 0</value>
</property> </property>
<property> <property>
<name>thermal_conductivity</name> <name>specific_heat_capacity</name>
<type>Constant</type> <type>Constant</type>
<value>0.1</value> <value>4068</value>
</property> </property>
<property> <property>
<name>density</name> <name>density</name>
...@@ -84,11 +84,6 @@ ...@@ -84,11 +84,6 @@
<type>Constant</type> <type>Constant</type>
<value>2000</value> <value>2000</value>
</property> </property>
<property>
<name>thermal_conductivity</name>
<type>Constant</type>
<value>2.5</value>
</property>
<property> <property>
<name>density</name> <name>density</name>
<type>Constant</type> <type>Constant</type>
...@@ -104,11 +99,6 @@ ...@@ -104,11 +99,6 @@
<type>Constant</type> <type>Constant</type>
<value>1000</value> <value>1000</value>
</property> </property>
<property>
<name>thermal_conductivity</name>
<type>Constant</type>
<value>3.2</value>
</property>
<property> <property>
<name>density</name> <name>density</name>
<type>Constant</type> <type>Constant</type>
...@@ -118,6 +108,16 @@ ...@@ -118,6 +108,16 @@
</phase> </phase>
</phases> </phases>
<properties> <properties>
<property>
<name>porosity</name>
<type>Constant</type>
<value>0</value>
</property>
<property>
<name>thermal_conductivity</name>
<type>Constant</type>
<value>2.5</value>
</property>
<property> <property>
<name>thermal_longitudinal_dispersivity</name> <name>thermal_longitudinal_dispersivity</name>
<type>Constant</type> <type>Constant</type>
......
+++
author = "Boyan Meng, Chaofan Chen, Haibing Shao"
date = "2019-11-11T11:52:00+01:00"
title = "BHE with groundwater advection"
weight = 123
project = "Parabolic/T/3D_BHE_GW_advection/BHE_GW_advection.prj"
[menu]
[menu.benchmarks]
parent = "heat-transport-bhe"
+++
{{< data-link >}}
## Problem description
When groundwater flow is present, advective heat transport in the soil matrix
has to be considered. Hence, the governing equation for advective and
conductive heat transport in an isotropic porous media can be expressed (in
2D form) as follows:
$$
\begin{equation}
\rho c \frac{\partial T}{\partial t} + \rho_wc_w \left(u_x\frac{\partial
T}{\partial x} + u_y\frac{\partial T}{\partial y}\right) - \lambda
\left(\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial
y^2}\right) = 0
\end{equation}
$$
where **$u$**$=(u_x,u_y)$ denotes the Darcy velocity. For simplification, we
assume that **$u$** is uniform and in the $x$ direction, i.e. $u_y=0$. In
addition, the porous medium is assumed to be infinite with homogeneous
initial temperature as well as hydraulic/thermal parameters. Under these
assumptions, the analytical solution for the ground temperature response to a
constant and uniform line source located at (0, 0) with infinite length along
the $z$ direction is expressed as (Diao et al. (2004)):
$$
\begin{equation}
\Delta T(x,y,t)=\frac{q_L}{4\pi\lambda}{\rm
exp}\left[\frac{v_Tx}{2a}\right]\int_{0}^{v_T^2t/4a} \frac{1}{\psi}{\rm
exp}\left[-\psi-\frac{v_T^2(x^2+y^2)}{16a^2\psi}\right] {\rm d}\psi
\end{equation}
$$
in which $v_T=u_x\rho_w c_w/\rho c$ is the effective heat transport velocity
(Molina-Giraldo et al. (2011)) and $q_L$ is the continuous heat exchange rate
per unit length. As a common practice, the BHE can be approximated by a line
source (e.g. Eskilson (1987) and Diao et al. (2004)). In cases where the BHE
penetrates the entire depth of the 3D porous medium, the above analytical
solution can be applied to solve the spatio-temporal distribution of the
induced ground temperature.
## Model Setup
The input files for the full simulation including the analytical solution for
the soil temperature can be found [here](../BHE_GW_advection_2years.zip). The
geometry of the model is
illustrated in Figure 1. The depth of the model domain is 15 m with an areal
extent of 80 m x 80 m. The BHE is 1U-type and is
represented by a straight line located at $x=0$ m and $y=30$ m. In this
benchmark the groundwater flow is set in the $y$ direction. Accordingly, the
mesh was intentionally extended downstream of the BHE, so that the boundary
effects on the ground temperature distribution can be neglected even for the
long-term simulation. Detailed parameters for the soil heat transport model
can be found in the following table.
| Parameter | Symbol | Value | Unit |
| -------------------------------------------------- |:------------------ | -------------------:| --------------------------: |
| Thermal conductivity of the porous medium | $\lambda$ | 2.5 | $\mathrm{W m^{-1} K^{-1}}$ |
| Volumetric heat capacity of the porous medium | $\rho c$ | $2.818\times10^{6}$ | $\mathrm{Jm^{-3}K^{-1}}$ |
| Length of the BHE | $L$ | 15 | $\mathrm{m}$ |
| Darcy velocity | $u_y$ | $1\times10^{-7}$ | $\mathrm{m s^{-1}}$ |
| Specific heat exchange rate of the BHE | $q_L$ | 20 | $\mathrm{W m^{-1}}$ |
| Initial ground temperature | $T_{ini}$ | 25 | $^{\circ}$C |
The BHE parameters are only relevant for the numerical model and are adopted
from the [3D Beier sandbox
benchmark](https://www.opengeosys.org/docs/benchmarks/heat-transport-bhe/3d_bei
er_sandbox/).
{{< img src="../3D_BHE_GW_advection_figures/mesh.png" width="150">}}
Figure 1: Geometry and mesh of the BHE model
## Results
In Figure 2, the numerically simulated ground temperature distribution from
OGS-6 is shown for the $z=-7$ m plane after $t=2$ years. Also, the result is
compared with the moving line source analytical solution (evaluated using
MATLAB®) in Figure 3. The comparison demonstrates that the numerical results
and analytical solution match very well as the maximum relative error of
ground temperature is less than 0.2 \%. The largest difference is found near
the BHE node towards which the analytical solution approaches infinity.
{{< img src="../3D_BHE_GW_advection_figures/temperature_soil_2years.png"
width="150">}}
Figure 2: Ground temperature distribution after two years at $z=-7$ m.
{{< img src="../3D_BHE_GW_advection_figures/rel_err.png" width="150">}}
Figure 3: Comparison of OGS-6 results and analytical solution. Note the
singularity of the analytical solution at the BHE node.
## References
[1] Diao, N., Li, Q., & Fang, Z. (2004). Heat transfer in ground heat
exchangers with groundwater advection. International Journal of Thermal
Sciences, 43(12), 1203-1211.
[2] Molina-Giraldo, N., Blum, P., Zhu, K., Bayer, P., & Fang, Z. (2011). A
moving finite line source model to simulate borehole heat exchangers with
groundwater advection. International Journal of Thermal Sciences, 50(12),
2506-2513.
[3] P. Eskilson, Thermal analysis of heat extraction boreholes, Ph.D. Thesis,
University of Lund, Lund, Sweden, 1987.
web/content/docs/benchmarks/heat-transport-bhe/3D_BHE_GW_advection_figures/mesh.PNG

131 B

web/content/docs/benchmarks/heat-transport-bhe/3D_BHE_GW_advection_figures/rel_err.PNG

131 B

0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment