From 78aa156fc5c7cc89051748385c93143c5fd1e1e4 Mon Sep 17 00:00:00 2001
From: "Dmitry Yu. Naumov" <github@naumov.de>
Date: Tue, 30 Oct 2018 16:14:50 +0100
Subject: [PATCH] [PL] Rewrite BHE process.
---
.../BHEBottomDirichletBoundaryCondition.cpp | 118 +--
.../BHEBottomDirichletBoundaryCondition.h | 26 +-
.../BHEInflowDirichletBoundaryCondition.cpp | 113 ---
.../BHEInflowDirichletBoundaryCondition.h | 73 +-
ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h | 458 ------------
ProcessLib/HeatTransportBHE/BHE/BHECommon.h | 48 ++
ProcessLib/HeatTransportBHE/BHE/BHETypes.h | 26 +
ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp | 696 +++---------------
ProcessLib/HeatTransportBHE/BHE/BHE_1U.h | 426 +++++------
.../HeatTransportBHE/BHE/BoreholeGeometry.h | 8 +
.../HeatTransportBHE/BHE/CreateBHE1U.cpp | 277 ++-----
ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.h | 21 +-
.../BHE/CreateFlowAndTemperatureControl.h | 81 ++
.../BHE/FlowAndTemperatureControl.h | 73 ++
ProcessLib/HeatTransportBHE/BHE/MeshUtils.cpp | 104 +--
ProcessLib/HeatTransportBHE/BHE/MeshUtils.h | 13 +-
ProcessLib/HeatTransportBHE/BHE/Physics.h | 64 ++
ProcessLib/HeatTransportBHE/BHE/Pipe.cpp | 32 +
ProcessLib/HeatTransportBHE/BHE/Pipe.h | 43 ++
.../BHE/PipeConfiguration1U.h | 37 +
.../HeatTransportBHE/BHE/PipeParameters.h | 67 --
...ntParameters.h => RefrigerantProperties.h} | 20 +-
.../BHE/ThermoMechanicalFlowProperties.h | 71 ++
.../CreateHeatTransportBHEProcess.cpp | 124 +---
.../HeatTransportBHEProcess.cpp | 185 +++--
.../HeatTransportBHEProcess.h | 9 +-
.../HeatTransportBHEProcessData.h | 15 +-
.../LocalAssemblers/CreateLocalAssemblers.h | 36 +-
.../HeatTransportBHELocalAssemblerBHE-impl.h | 199 +++++
.../HeatTransportBHELocalAssemblerBHE.h | 45 +-
.../HeatTransportBHELocalAssemblerBHE_impl.h | 223 ------
...HeatTransportBHELocalAssemblerSoil-impl.h} | 47 +-
.../HeatTransportBHELocalAssemblerSoil.h | 10 +-
...eatTransportBHEProcessAssemblerInterface.h | 21 -
.../LocalAssemblers/IntegrationPointDataBHE.h | 54 +-
.../IntegrationPointDataSoil.h | 7 +-
.../LocalAssemblers/LocalDataInitializer.h | 206 ++----
37 files changed, 1402 insertions(+), 2674 deletions(-)
delete mode 100644 ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.cpp
delete mode 100644 ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h
create mode 100644 ProcessLib/HeatTransportBHE/BHE/BHECommon.h
create mode 100644 ProcessLib/HeatTransportBHE/BHE/BHETypes.h
create mode 100644 ProcessLib/HeatTransportBHE/BHE/CreateFlowAndTemperatureControl.h
create mode 100644 ProcessLib/HeatTransportBHE/BHE/FlowAndTemperatureControl.h
create mode 100644 ProcessLib/HeatTransportBHE/BHE/Physics.h
create mode 100644 ProcessLib/HeatTransportBHE/BHE/Pipe.cpp
create mode 100644 ProcessLib/HeatTransportBHE/BHE/Pipe.h
create mode 100644 ProcessLib/HeatTransportBHE/BHE/PipeConfiguration1U.h
delete mode 100644 ProcessLib/HeatTransportBHE/BHE/PipeParameters.h
rename ProcessLib/HeatTransportBHE/BHE/{RefrigerantParameters.h => RefrigerantProperties.h} (59%)
create mode 100644 ProcessLib/HeatTransportBHE/BHE/ThermoMechanicalFlowProperties.h
create mode 100644 ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
delete mode 100644 ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE_impl.h
rename ProcessLib/HeatTransportBHE/LocalAssemblers/{HeatTransportBHELocalAssemblerSoil_impl.h => HeatTransportBHELocalAssemblerSoil-impl.h} (74%)
diff --git a/ProcessLib/BoundaryCondition/BHEBottomDirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BHEBottomDirichletBoundaryCondition.cpp
index b87ec86a71c..5e7ef9c5820 100644
--- a/ProcessLib/BoundaryCondition/BHEBottomDirichletBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryCondition/BHEBottomDirichletBoundaryCondition.cpp
@@ -8,106 +8,50 @@
*/
#include "BHEBottomDirichletBoundaryCondition.h"
-
-#include <algorithm>
-#include <logog/include/logog.hpp>
-#include <vector>
-#include "ProcessLib/Utils/ProcessUtils.h"
+#include "BaseLib/Error.h"
namespace ProcessLib
{
-BHEBottomDirichletBoundaryCondition::BHEBottomDirichletBoundaryCondition(
- GlobalIndexType const global_idx_T_in_bottom,
- GlobalIndexType const global_idx_T_out_bottom,
- MeshLib::Mesh const& bulk_mesh,
- std::vector<MeshLib::Node*> const& vec_outflow_bc_nodes,
- int const variable_id,
- int const component_id)
- : _bulk_mesh(bulk_mesh)
-{
- DBUG(
- "Found %d nodes for BHE bottom Dirichlet BCs for the variable %d "
- "and "
- "component %d",
- vec_outflow_bc_nodes.size(), variable_id, component_id);
-
- MeshLib::MeshSubset bc_mesh_subset{_bulk_mesh, vec_outflow_bc_nodes};
-
- // create memory to store Tout values
- _T_in_values.ids.clear();
- _T_in_values.values.clear();
-
- _bc_values.ids.clear();
- _bc_values.values.clear();
-
- // convert mesh node ids to global index for the given component
- _bc_values.ids.reserve(bc_mesh_subset.getNumberOfNodes());
- _bc_values.values.reserve(bc_mesh_subset.getNumberOfNodes());
-
- const auto g_T_out_idx = global_idx_T_out_bottom;
- if (g_T_out_idx >= 0)
- {
- _bc_values.ids.emplace_back(g_T_out_idx);
- _bc_values.values.emplace_back(298.15);
- }
-
- const auto g_T_in_idx = global_idx_T_in_bottom;
-
- if (g_T_in_idx >= 0)
- {
- _T_in_values.ids.emplace_back(g_T_in_idx);
- _T_in_values.values.emplace_back(298.15);
- }
-}
-
void BHEBottomDirichletBoundaryCondition::getEssentialBCValues(
const double /*t*/, GlobalVector const& x,
NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
{
- bc_values.ids.clear();
- bc_values.values.clear();
-
- bc_values.ids.resize(_bc_values.ids.size());
- bc_values.values.resize(_bc_values.values.size());
-
- const std::size_t n_nodes = _T_in_values.ids.size();
- for (std::size_t i = 0; i < n_nodes; i++)
- {
- bc_values.ids[i] = _bc_values.ids[i];
- // here, the outflow temperature is always
- // the same as the inflow temperature
- // get the inflow temperature from here.
- bc_values.values[i] = x[_T_in_values.ids[i]];
- }
-}
-
-// update new values and corresponding indices.
-void BHEBottomDirichletBoundaryCondition::preTimestep(const double /*t*/,
- const GlobalVector& x)
-{
- // At the bottom of each BHE, the outflow temperature
- // is the same as the inflow temperature.
- // Here the task is to get the inflow temperature and
- // save it locally
- auto const n_nodes = _bc_values.ids.size();
- for (std::size_t i = 0; i < n_nodes; i++)
- {
- // read the T_out
- _T_in_values.values[i] = x[_T_in_values.ids[i]];
- }
+ bc_values.ids.resize(1);
+ bc_values.values.resize(1);
+
+ bc_values.ids[0] = _in_out_global_indices.second;
+ // here, the outflow temperature is always
+ // the same as the inflow temperature
+ // get the inflow temperature from here.
+ bc_values.values[0] = x[_in_out_global_indices.first];
}
std::unique_ptr<BHEBottomDirichletBoundaryCondition>
createBHEBottomDirichletBoundaryCondition(
- GlobalIndexType global_idx_T_in_bottom,
- GlobalIndexType global_idx_T_out_bottom, MeshLib::Mesh const& bulk_mesh,
- std::vector<MeshLib::Node*> const& vec_outflow_bc_nodes,
- int const variable_id, int const component_id)
+ std::pair<GlobalIndexType, GlobalIndexType>&& in_out_global_indices)
{
- DBUG("Constructing BHEBottomDirichletBoundaryCondition from config.");
+ DBUG("Constructing BHEBottomDirichletBoundaryCondition.");
+
+ // In case of partitioned mesh the boundary could be empty, i.e. there is no
+ // boundary condition.
+#ifdef USE_PETSC
+ // For this special boundary condition the boundary condition is not empty
+ // if the global indices are non-negative.
+ if (in_out_global_indices.first < 0 && in_out_global_indices.second < 0)
+ {
+ return nullptr;
+ }
+ // If only one of the global indices (in or out) is negative the
+ // implementation is not valid.
+ if (in_out_global_indices.first < 0 || in_out_global_indices.second < 0)
+ {
+ OGS_FATAL(
+ "The partition cuts the BHE into two independent parts. This "
+ "behaviour is not implemented.");
+ }
+#endif // USE_PETSC
return std::make_unique<BHEBottomDirichletBoundaryCondition>(
- global_idx_T_in_bottom, global_idx_T_out_bottom, bulk_mesh,
- vec_outflow_bc_nodes, variable_id, component_id);
+ std::move(in_out_global_indices));
}
} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/BHEBottomDirichletBoundaryCondition.h b/ProcessLib/BoundaryCondition/BHEBottomDirichletBoundaryCondition.h
index 063eac02a09..df3e2b28039 100644
--- a/ProcessLib/BoundaryCondition/BHEBottomDirichletBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/BHEBottomDirichletBoundaryCondition.h
@@ -10,10 +10,7 @@
#pragma once
#include "BoundaryCondition.h"
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
#include "NumLib/IndexValueVector.h"
-#include "ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h"
-#include "ProcessLib/Parameter/Parameter.h"
namespace ProcessLib
{
@@ -21,31 +18,20 @@ class BHEBottomDirichletBoundaryCondition final : public BoundaryCondition
{
public:
BHEBottomDirichletBoundaryCondition(
- GlobalIndexType global_idx_T_in_bottom,
- GlobalIndexType global_idx_T_out_bottom,
- MeshLib::Mesh const& bulk_mesh,
- std::vector<MeshLib::Node*> const& vec_outflow_bc_nodes,
- int const variable_id,
- int const component_id);
+ std::pair<GlobalIndexType, GlobalIndexType>&& in_out_global_indices)
+ : _in_out_global_indices(std::move(in_out_global_indices))
+ {
+ }
void getEssentialBCValues(
const double t, GlobalVector const& x,
NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override;
- void preTimestep(const double t, const GlobalVector& x) override;
-
private:
- MeshLib::Mesh const& _bulk_mesh;
-
- NumLib::IndexValueVector<GlobalIndexType> _bc_values;
-
- NumLib::IndexValueVector<GlobalIndexType> _T_in_values;
+ std::pair<GlobalIndexType, GlobalIndexType> const _in_out_global_indices;
};
std::unique_ptr<BHEBottomDirichletBoundaryCondition>
createBHEBottomDirichletBoundaryCondition(
- GlobalIndexType global_idx_T_in_bottom,
- GlobalIndexType global_idx_T_out_bottom, MeshLib::Mesh const& bulk_mesh,
- std::vector<MeshLib::Node*> const& vec_outflow_bc_nodes,
- int const variable_id, int const component_id);
+ std::pair<GlobalIndexType, GlobalIndexType>&& in_out_global_indices);
} // namespace ProcessLib
diff --git a/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.cpp b/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.cpp
deleted file mode 100644
index f26e9b6623c..00000000000
--- a/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.cpp
+++ /dev/null
@@ -1,113 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2018, 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 "BHEInflowDirichletBoundaryCondition.h"
-
-#include <algorithm>
-#include <logog/include/logog.hpp>
-#include <vector>
-#include "ProcessLib/Utils/ProcessUtils.h"
-
-namespace ProcessLib
-{
-BHEInflowDirichletBoundaryCondition::BHEInflowDirichletBoundaryCondition(
- std::pair<GlobalIndexType, GlobalIndexType>&& in_out_global_indices,
- MeshLib::Mesh const& bc_mesh,
- std::vector<MeshLib::Node*> const& vec_inflow_bc_nodes,
- int const variable_id,
- int const component_id,
- std::unique_ptr<ProcessLib::HeatTransportBHE::BHE::BHEAbstract> const&
- pt_bhe)
- : _bc_mesh(bc_mesh), _pt_bhe(pt_bhe)
-{
- DBUG(
- "Found %d nodes for BHE Inflow Dirichlet BCs for the variable %d and "
- "component %d",
- vec_inflow_bc_nodes.size(), variable_id, component_id);
-
- MeshLib::MeshSubset bc_mesh_subset{_bc_mesh, vec_inflow_bc_nodes};
-
- // create memory to store Tout values
- _T_out_values.clear();
- _T_out_indices.clear();
-
- _bc_values.ids.clear();
- _bc_values.values.clear();
-
- // convert mesh node ids to global index for the given component
- assert(bc_mesh_subset.getNumberOfNodes() == 1);
- _bc_values.ids.reserve(bc_mesh_subset.getNumberOfNodes());
- _bc_values.values.reserve(bc_mesh_subset.getNumberOfNodes());
-
- // that might be slow, but only done once
- const auto g_idx_T_in = in_out_global_indices.first;
- const auto g_idx_T_out = in_out_global_indices.second;
-
- if (g_idx_T_in >= 0 && g_idx_T_out >= 0)
- {
- _T_out_indices.emplace_back(g_idx_T_out);
- _T_out_values.emplace_back(320.0 /*using initial value*/);
- _bc_values.ids.emplace_back(g_idx_T_in);
- _bc_values.values.emplace_back(320.0 /*using initial value*/);
- }
-}
-
-void BHEInflowDirichletBoundaryCondition::getEssentialBCValues(
- const double t, GlobalVector const& x,
- NumLib::IndexValueVector<GlobalIndexType>& bc_values) const
-{
- bc_values.ids.clear();
- bc_values.values.clear();
-
- bc_values.ids.resize(_bc_values.ids.size());
- bc_values.values.resize(_bc_values.values.size());
-
- const std::size_t n_nodes = _T_out_values.size();
-
- for (std::size_t i = 0; i < n_nodes; i++)
- {
- bc_values.ids[i] = _bc_values.ids[i];
- // here call the corresponding BHE functions
- auto const tmp_T_out = x[_T_out_indices[i]];
- bc_values.values[i] = _pt_bhe->getTinByTout(tmp_T_out, t);
- }
-}
-
-// update new values and corresponding indices.
-void BHEInflowDirichletBoundaryCondition::preTimestep(const double /*t*/,
- const GlobalVector& x)
-{
- // for each BHE, the inflow temperature is dependent on
- // the ouflow temperature of the BHE.
- // Here the task is to get the outflow temperature and
- // save it locally
- auto const n_nodes = _bc_values.ids.size();
- for (std::size_t i = 0; i < n_nodes; i++)
- {
- // read the T_out
- _T_out_values[i] = x[_T_out_indices[i]];
- }
-}
-
-std::unique_ptr<BHEInflowDirichletBoundaryCondition>
-createBHEInflowDirichletBoundaryCondition(
- std::pair<GlobalIndexType, GlobalIndexType>&& in_out_global_indices,
- MeshLib::Mesh const& bc_mesh,
- std::vector<MeshLib::Node*> const& vec_inflow_bc_nodes,
- int const variable_id, int const component_id,
- std::unique_ptr<ProcessLib::HeatTransportBHE::BHE::BHEAbstract> const&
- pt_bhe)
-{
- DBUG("Constructing BHEInflowDirichletBoundaryCondition.");
-
- return std::make_unique<BHEInflowDirichletBoundaryCondition>(
- std::move(in_out_global_indices), bc_mesh, vec_inflow_bc_nodes,
- variable_id, component_id, pt_bhe);
-}
-} // namespace ProcessLib
\ No newline at end of file
diff --git a/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.h b/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.h
index 442473b1492..b955493cda3 100644
--- a/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.h
+++ b/ProcessLib/BoundaryCondition/BHEInflowDirichletBoundaryCondition.h
@@ -10,51 +10,68 @@
#pragma once
#include "BoundaryCondition.h"
-#include "NumLib/DOF/LocalToGlobalIndexMap.h"
#include "NumLib/IndexValueVector.h"
-#include "ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h"
-#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/HeatTransportBHE/BHE/BHETypes.h"
namespace ProcessLib
{
+template <typename BHEType>
class BHEInflowDirichletBoundaryCondition final : public BoundaryCondition
{
public:
BHEInflowDirichletBoundaryCondition(
std::pair<GlobalIndexType, GlobalIndexType>&& in_out_global_indices,
- MeshLib::Mesh const& bc_mesh,
- std::vector<MeshLib::Node*> const& vec_inflow_bc_nodes,
- int const variable_id,
- int const component_id,
- std::unique_ptr<ProcessLib::HeatTransportBHE::BHE::BHEAbstract> const&
- pt_bhe);
+ BHEType& bhe)
+ : _in_out_global_indices(std::move(in_out_global_indices)), _bhe(bhe)
+ {
+ }
void getEssentialBCValues(
const double t, GlobalVector const& x,
- NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override;
+ NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override
+ {
+ bc_values.ids.resize(1);
+ bc_values.values.resize(1);
- void preTimestep(const double t, const GlobalVector& x) override;
+ bc_values.ids[0] = _in_out_global_indices.first;
+ // here call the corresponding BHE functions
+ auto const T_out = x[_in_out_global_indices.second];
+ bc_values.values[0] = _bhe.getTinByTout(T_out, t);
+ }
private:
- // TODO (haibing) re-organize as the bottom BC data structure
- MeshLib::Mesh const& _bc_mesh;
-
- /// Stores the results of the outflow temperatures per boundary node.
- std::vector<double> _T_out_values;
- std::vector<GlobalIndexType> _T_out_indices;
-
- NumLib::IndexValueVector<GlobalIndexType> _bc_values;
-
- std::unique_ptr<ProcessLib::HeatTransportBHE::BHE::BHEAbstract> const&
- _pt_bhe;
+ std::pair<GlobalIndexType, GlobalIndexType> const _in_out_global_indices;
+ BHEType& _bhe;
};
-std::unique_ptr<BHEInflowDirichletBoundaryCondition>
+template <typename BHEType>
+std::unique_ptr<BHEInflowDirichletBoundaryCondition<BHEType>>
createBHEInflowDirichletBoundaryCondition(
std::pair<GlobalIndexType, GlobalIndexType>&& in_out_global_indices,
- MeshLib::Mesh const& bc_mesh,
- std::vector<MeshLib::Node*> const& vec_outflow_bc_nodes,
- int const variable_id, int const component_id,
- std::unique_ptr<ProcessLib::HeatTransportBHE::BHE::BHEAbstract> const&
- pt_bhe);
+ BHEType& bhe)
+{
+ DBUG("Constructing BHEInflowDirichletBoundaryCondition.");
+
+ // In case of partitioned mesh the boundary could be empty, i.e. there is no
+ // boundary condition.
+#ifdef USE_PETSC
+ // For this special boundary condition the boundary condition is not empty
+ // if the global indices are non-negative.
+ if (in_out_global_indices.first < 0 && in_out_global_indices.second < 0)
+ {
+ return nullptr;
+ }
+ // If only one of the global indices (in or out) is negative the
+ // implementation is not valid.
+ if (in_out_global_indices.first < 0 || in_out_global_indices.second < 0)
+ {
+ OGS_FATAL(
+ "The partition cuts the BHE into two independent parts. This "
+ "behaviour is not implemented.");
+ }
+#endif // USE_PETSC
+
+ return std::make_unique<BHEInflowDirichletBoundaryCondition<BHEType>>(
+ std::move(in_out_global_indices), bhe);
+}
} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h b/ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h
deleted file mode 100644
index 09126e877d2..00000000000
--- a/ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h
+++ /dev/null
@@ -1,458 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
- * Distributed under a Modified BSD License.
- * See accompanying file LICENSE.txt or
- * http://www.opengeosys.org/project/license
- *
- */
-
-/**
- * \file BHEAbstract.h
- * 2014/06/04 HS inital implementation
- * borehole heat exchanger abstract class
- *
- * 1) Diersch_2011_CG
- * Two very important references to understand this class implementations are:
- * H.-J.G. Diersch, D. Bauer, W. Heidemann, W. R黨aak, P. Sch鋞zl,
- * Finite element modeling of borehole heat exchanger systems:
- * Part 1. Fundamentals, Computers & Geosciences,
- * Volume 37, Issue 8, August 2011, Pages 1122-1135, ISSN 0098-3004,
- * http://dx.doi.org/10.1016/j.cageo.2010.08.003.
- *
- * 2) FEFLOW_2014_Springer
- * FEFLOW: Finite Element Modeling of Flow, Mass and Heat Transport in Porous
- * and Fractured Media Diersch, Hans-Joerg, 2014, XXXV, 996 p, Springer.
- *
- */
-
-#pragma once
-
-#include <Eigen/Eigen>
-#include <map>
-#include "GeoLib/Polyline.h"
-#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
-#include "ProcessLib/Utils/ProcessUtils.h"
-#include "boost/math/constants/constants.hpp"
-
-#include "BoreholeGeometry.h"
-#include "GroutParameters.h"
-#include "PipeParameters.h"
-#include "RefrigerantParameters.h"
-
-namespace ProcessLib
-{
-namespace HeatTransportBHE
-{
-namespace BHE // namespace of borehole heat exchanger
-{
-/**
- * list the types of borehole heat exchanger
- */
-enum class BHE_TYPE
-{
- TYPE_2U, // two u-tube borehole heat exchanger
- TYPE_1U, // one u-tube borehole heat exchanger
- TYPE_CXC, // coaxial pipe with annualar inlet
- TYPE_CXA // coaxial pipe with centreed inlet
-};
-
-enum class BHE_BOUNDARY_TYPE
-{
- FIXED_INFLOW_TEMP_BOUNDARY,
- FIXED_INFLOW_TEMP_CURVE_BOUNDARY,
- POWER_IN_WATT_BOUNDARY,
- POWER_IN_WATT_CURVE_FIXED_DT_BOUNDARY,
- BUILDING_POWER_IN_WATT_CURVE_FIXED_DT_BOUNDARY,
- BUILDING_POWER_IN_WATT_CURVE_FIXED_FLOW_RATE_BOUNDARY,
- POWER_IN_WATT_CURVE_FIXED_FLOW_RATE_BOUNDARY,
- FIXED_TEMP_DIFF_BOUNDARY
-};
-
-/**
- * discharge type of the 2U BHE
- */
-enum class BHE_DISCHARGE_TYPE
-{
- BHE_DISCHARGE_TYPE_PARALLEL, // parallel discharge
- BHE_DISCHARGE_TYPE_SERIAL // serial discharge
-};
-
-class BHEAbstract
-{
-public:
- struct ExternallyDefinedRaRb
- {
- /**
- * whether or not using external given borehole thermal resistance
- * values Ra, Rb
- */
- bool use_extern_Ra_Rb;
-
- /**
- * external given borehole internal thermal resistance value
- */
- double ext_Ra;
-
- /**
- * external given borehole thermal resistance value
- */
- double ext_Rb;
- };
-
- struct ExternallyDefinedThermalResistances
- {
- /**
- * whether or not using user defined borehole thermal resistance Rfig,
- * Rfog, Rgg, Rgs
- */
- bool if_use_defined_therm_resis;
-
- /**
- * external given borehole internal thermal resistance value
- */
- double ext_Rfig;
-
- /**
- * external given borehole internal thermal resistance value
- */
- double ext_Rfog;
-
- /**
- * external given borehole internal thermal resistance value
- */
- double ext_Rgg1;
-
- /**
- * external given borehole internal thermal resistance value
- */
- double ext_Rgg2;
-
- /**
- * external given borehole internal thermal resistance value
- */
- double ext_Rgs;
- };
-
- /**
- * constructor
- */
- BHEAbstract(
- const std::string name_,
- const BHE_TYPE bhe_type_,
- BoreholeGeometry borehole_geometry_,
- PipeParameters pipe_param_,
- RefrigerantParameters refrigerant_param_,
- GroutParameters grout_param_,
- ExternallyDefinedRaRb extern_Ra_Rb_,
- ExternallyDefinedThermalResistances extern_def_thermal_resistances_,
- std::map<std::string,
- std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
- bhe_curves_,
- BHE_BOUNDARY_TYPE my_bound_type =
- BHE_BOUNDARY_TYPE::FIXED_INFLOW_TEMP_BOUNDARY,
- bool if_use_ext_Ra_Rb = false,
- bool user_defined_R_vals = false,
- bool if_flowrate_curve = false)
- : name(name_),
- bhe_type(bhe_type_),
- boundary_type(my_bound_type),
- borehole_geometry(borehole_geometry_),
- pipe_param(pipe_param_),
- refrigerant_param(refrigerant_param_),
- grout_param(grout_param_),
- extern_Ra_Rb(extern_Ra_Rb_),
- extern_def_thermal_resistances(extern_def_thermal_resistances_),
- bhe_curves(bhe_curves_),
- use_flowrate_curve(if_flowrate_curve),
- if_use_ext_Ra_Rb(if_use_ext_Ra_Rb),
- user_defined_R_vals(user_defined_R_vals),
- PI(boost::math::constants::pi<double>()){};
-
- /**
- * destructor
- */
- virtual ~BHEAbstract() = default;
-
- /**
- * return the number of unknowns needed for this BHE
- * abstract function, need to be realized.
- */
- virtual std::size_t getNumUnknowns() const = 0;
-
- /**
- * initialization calcultion,
- * need to be overwritten.
- */
- virtual void initialize() = 0;
-
- /**
- * update all parameters based on the new flow rate
- * not necessarily needs to be overwritten.
- */
- virtual void updateFlowRate(double new_flow_rate)
- {
- Q_r = new_flow_rate;
- initialize();
- };
-
- virtual void updateFlowRateFromCurve(double current_time) = 0;
-
- /**
- * thermal resistance calculation,
- * need to be overwritten.
- */
- virtual void calcThermalResistances() = 0;
-
- /**
- * flow velocity inside the pipeline
- */
- double calcPipeFlowVelocity(double const& flow_rate,
- double const& pipe_diameter)
- {
- return 4.0 * flow_rate / (PI * pipe_diameter * pipe_diameter);
- }
-
- double calcPrandtlNumber(double const& viscosity,
- double const& heat_capacity,
- double const& heat_conductivity)
- {
- return viscosity * heat_capacity / heat_conductivity;
- }
-
- double calcRenoldsNumber(double const velocity_norm,
- double const pipe_diameter,
- double const viscosity,
- double const density)
- {
- return velocity_norm * pipe_diameter / (viscosity / density);
- };
- double calcNusseltNumber(double const pipe_diameter,
- double const pipe_length)
- {
- double tmp_Nu(0.0);
- double gamma(0.0), xi(0.0);
-
- if (Re < 2300.0)
- {
- tmp_Nu = 4.364;
- }
- else if (Re >= 2300.0 && Re < 10000.0)
- {
- gamma = (Re - 2300) / (10000 - 2300);
-
- tmp_Nu = (1.0 - gamma) * 4.364;
- tmp_Nu +=
- gamma *
- ((0.0308 / 8.0 * 1.0e4 * Pr) /
- (1.0 + 12.7 * std::sqrt(0.0308 / 8.0) *
- (std::pow(Pr, 2.0 / 3.0) - 1.0)) *
- (1.0 + std::pow(pipe_diameter / pipe_length, 2.0 / 3.0)));
- }
- else if (Re > 10000.0)
- {
- xi = pow(1.8 * std::log10(Re) - 1.5, -2.0);
- tmp_Nu = (xi / 8.0 * Re * Re) /
- (1.0 + 12.7 * std::sqrt(xi / 8.0) *
- (std::pow(Pr, 2.0 / 3.0) - 1.0)) *
- (1.0 + std::pow(pipe_diameter / pipe_length, 2.0 / 3.0));
- }
-
- return tmp_Nu;
- };
-
- /**
- * heat transfer coefficient,
- * need to be overwritten.
- */
- virtual void calcHeatTransferCoefficients() = 0;
-
- /**
- * return the coeff of mass matrix,
- * depending on the index of unknown.
- */
- virtual double getMassCoeff(std::size_t idx_unknown) const = 0;
-
- /**
- * return the coeff of laplace matrix,
- * depending on the index of unknown.
- */
- virtual void getLaplaceMatrix(std::size_t idx_unknown,
- Eigen::MatrixXd& mat_laplace) const = 0;
-
- /**
- * return the coeff of advection matrix,
- * depending on the index of unknown.
- */
- virtual void getAdvectionVector(std::size_t idx_unknown,
- Eigen::VectorXd& vec_advection) const = 0;
-
- /**
- * return the coeff of boundary heat exchange matrix,
- * depending on the index of unknown.
- */
- virtual double getBoundaryHeatExchangeCoeff(
- std::size_t idx_unknown) const = 0;
-
- /**
- * return the inflow temperature based on outflow temperature and fixed
- * power.
- */
- virtual double getTinByTout(double T_in, double current_time) = 0;
-
- /**
- * return the _R_matrix, _R_pi_s_matrix, _R_s_matrix
- */
- virtual void setRMatrices(
- const int idx_bhe_unknowns, const int NumNodes,
- Eigen::MatrixXd& matBHE_loc_R,
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
- R_matrix,
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
- R_pi_s_matrix,
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
- R_s_matrix) const = 0;
-
- virtual std::vector<std::pair<int, int>> const&
- inflowOutflowBcComponentIds() const = 0;
-
-public:
- /**
- * name of the borehole heat exchanger
- */
- const std::string name;
-
- /**
- * the type of this BHE
- */
- const BHE_TYPE bhe_type;
-
- /**
- * the type of the boundary condition on this BHE
- */
- const BHE_BOUNDARY_TYPE boundary_type;
-
- /**
- * geometry of the borehole
- */
- BoreholeGeometry const borehole_geometry;
-
- /**
- * geometry of the pipes in the borehole
- */
- PipeParameters const pipe_param;
-
- /**
- * parameters of the refrigerant
- */
- RefrigerantParameters const refrigerant_param;
-
- /**
- * parameters of the grout
- */
- GroutParameters const grout_param;
-
- /**
- * Ra Rb values defined by the user externally
- */
- ExternallyDefinedRaRb const extern_Ra_Rb;
-
- /**
- * thermal resistance values defined by the user externally
- */
- ExternallyDefinedThermalResistances const extern_def_thermal_resistances;
-
- /**
- * power extracted from or injected into the BHE
- * unit is Watt
- * if value positive, then injecting power
- * if value negative, then extracting power
- */
- double power_in_watt_val;
-
- /**
- * temperature difference between inflow and
- * outflow pipelines
- */
- double delta_T_val;
-
- /**
- * threshold Q value for switching off the BHE
- * when using the Q_curve_fixed_dT B.C.
- */
- double threshold;
-
- /**
- * map strucutre that contains all the curves related to this BHE
- */
- std::map<std::string,
- std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
- bhe_curves;
-
- /**
- * power in watt curve
- */
- MathLib::PiecewiseLinearInterpolation* power_in_watt_curve;
-
- /**
- * heating COP curve
- */
- MathLib::PiecewiseLinearInterpolation* heating_cop_curve;
-
- /**
- * cooling COP curve
- */
- MathLib::PiecewiseLinearInterpolation* cooling_cop_curve;
-
- /**
- * refrigerant flow rate curve
- */
- MathLib::PiecewiseLinearInterpolation* flowrate_curve;
-
- /**
- * inflow temperature curve
- */
- MathLib::PiecewiseLinearInterpolation* inflow_temperature_curve;
-
- /**
- * use refrigerant flow rate curve
- */
- bool use_flowrate_curve;
-
- /**
- * whether external Ra an Rb values are supplied by the user
- */
- const bool if_use_ext_Ra_Rb;
-
- /**
- * whether external R values are supplied by the user
- */
- const bool user_defined_R_vals;
-
-protected:
- /**
- * total refrigerant flow discharge of BHE
- * unit is m^3/sec
- */
- double Q_r;
-
- const double PI;
-
- /**
- * Reynolds number
- */
- double Re;
-
- /**
- * Prandtl number
- */
- double Pr;
-
- /**
- * pipe distance
- */
- double omega;
-};
-} // end of namespace BHE
-} // end of namespace HeatTransportBHE
-} // end of namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHECommon.h b/ProcessLib/HeatTransportBHE/BHE/BHECommon.h
new file mode 100644
index 00000000000..188f0a9a067
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/BHECommon.h
@@ -0,0 +1,48 @@
+/**
+ * \file
+ *
+ * 2014/06/04 HS inital implementation
+ * borehole heat exchanger abstract class
+ *
+ * 1) Diersch_2011_CG
+ * Two very important references to understand this class implementations are:
+ * Diersch, H-JG, D. Bauer, W. Heidemann, Wolfram Rühaak, and Peter Schätzl.
+ * Finite element modeling of borehole heat exchanger systems:
+ * Part 1. Fundamentals, Computers & Geosciences,
+ * Volume 37, Issue 8, August 2011, Pages 1122-1135, ISSN 0098-3004,
+ * http://dx.doi.org/10.1016/j.cageo.2010.08.003.
+ *
+ * 2) FEFLOW_2014_Springer
+ * FEFLOW: Finite Element Modeling of Flow, Mass and Heat Transport in Porous
+ * and Fractured Media Diersch, Hans-Joerg, 2014, XXXV, 996 p, Springer.
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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 "BoreholeGeometry.h"
+#include "FlowAndTemperatureControl.h"
+#include "GroutParameters.h"
+#include "RefrigerantProperties.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+struct BHECommon
+{
+ BoreholeGeometry const borehole_geometry;
+ RefrigerantProperties const refrigerant;
+ GroutParameters const grout;
+ FlowAndTemperatureControl const flowAndTemperatureControl;
+};
+} // end of namespace BHE
+} // end of namespace HeatTransportBHE
+} // end of namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHETypes.h b/ProcessLib/HeatTransportBHE/BHE/BHETypes.h
new file mode 100644
index 00000000000..4c6a45a6ffb
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/BHETypes.h
@@ -0,0 +1,26 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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 <boost/variant.hpp>
+#include "BHE_1U.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+using BHETypes = boost::variant<BHE_1U>;
+} // 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 85bbd24c3d1..1fd83dc95d1 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
@@ -9,644 +9,170 @@
#include "BHE_1U.h"
-using namespace ProcessLib::HeatTransportBHE::BHE;
+#include <boost/math/constants/constants.hpp>
+#include "FlowAndTemperatureControl.h"
+#include "Physics.h"
-void BHE_1U::initialize()
+using namespace ProcessLib::HeatTransportBHE::BHE;
+namespace
{
- double tmp_u = calcPipeFlowVelocity(Q_r, 2.0 * pipe_param.r_inner);
- _u(0) = tmp_u;
- _u(1) = tmp_u;
-
- // calculate Renolds number
- Re = calcRenoldsNumber(std::abs(_u(0)),
- 2.0 * pipe_param.r_inner,
- refrigerant_param.mu_r,
- refrigerant_param.rho_r);
- // calculate Prandtl number
- Pr = calcPrandtlNumber(refrigerant_param.mu_r,
- refrigerant_param.heat_cap_r,
- refrigerant_param.lambda_r);
-
- // calculate Nusselt number
- double tmp_Nu =
- calcNusseltNumber(2.0 * pipe_param.r_inner, borehole_geometry.length);
- _Nu(0) = tmp_Nu;
- _Nu(1) = tmp_Nu;
-
- calcThermalResistances();
- calcHeatTransferCoefficients();
+double compute_R_gs(double const chi, double const R_g)
+{
+ return (1 - chi) * R_g;
}
-/**
- * calculate thermal resistance
- */
-void BHE_1U::calcThermalResistances()
+double compute_R_gg(double const chi, double const R_gs, double const R_ar,
+ double const R_g)
{
- // thermal resistance due to the grout transition
- double chi;
- double d0; // the average outer diameter of the pipes
- // double s; // diagonal distances of pipes
- double R_adv, R_con;
- double const& D = borehole_geometry.diameter;
- // double const& L = borehole_geometry.L;
- double const& lambda_r = refrigerant_param.lambda_r;
- double const& lambda_g = grout_param.lambda_g;
- double const& lambda_p = pipe_param.lambda_p;
-
- // thermal resistance due to thermal conductivity of the pip wall material
- // Eq. 36 in Diersch_2011_CG
- _R_adv_i1 = 1.0 / (_Nu(0) * lambda_r * PI);
- _R_adv_o1 = 1.0 / (_Nu(1) * lambda_r * PI);
-
- d0 = 2.0 * pipe_param.r_outer;
- // s = omega * std::sqrt(2);
- // Eq. 49
- _R_con_a_i1 = _R_con_a_o1 =
- std::log(pipe_param.r_outer / pipe_param.r_inner) /
- (2.0 * PI * lambda_p);
- // Eq. 51
- chi = std::log(std::sqrt(D * D + 2 * d0 * d0) / 2 / d0) /
- std::log(D / std::sqrt(2) / d0);
- if (extern_Ra_Rb.use_extern_Ra_Rb)
- {
- R_adv = 0.5 * (_R_adv_i1 + _R_adv_o1);
- R_con = 0.5 * (_R_con_a_i1 + _R_con_a_o1);
- _R_g = 2 * extern_Ra_Rb.ext_Rb - R_adv - R_con;
- }
- else
- {
- // Eq. 52
- _R_g = acosh((D * D + d0 * d0 - omega * omega) / (2 * D * d0)) /
- (2 * PI * lambda_g) * (1.601 - 0.888 * omega / D);
- }
- _R_con_b = chi * _R_g;
- // Eq. 29 and 30
- if (extern_def_thermal_resistances.if_use_defined_therm_resis)
- {
- _R_fig = extern_def_thermal_resistances.ext_Rfig;
- _R_fog = extern_def_thermal_resistances.ext_Rfog;
- }
- else
- {
- _R_fig = _R_adv_i1 + _R_con_a_i1 + _R_con_b;
- _R_fog = _R_adv_o1 + _R_con_a_o1 + _R_con_b;
- }
-
- // thermal resistance due to grout-soil exchange
- if (extern_def_thermal_resistances.if_use_defined_therm_resis)
- _R_gs = extern_def_thermal_resistances.ext_Rgs;
- else
- _R_gs = (1 - chi) * _R_g;
-
- // thermal resistance due to inter-grout exchange
- double R_ar;
- if (extern_Ra_Rb.use_extern_Ra_Rb)
- {
- R_ar = extern_Ra_Rb.ext_Ra - 2 * (R_adv + R_con);
- }
- else
- {
- R_ar = acosh((2.0 * omega * omega - d0 * d0) / d0 / d0) /
- (2.0 * PI * lambda_g);
- }
-
- if (extern_def_thermal_resistances.if_use_defined_therm_resis)
- _R_gg = extern_def_thermal_resistances.ext_Rgg1;
- else
- _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))
+ 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!");
}
- // check if constraints regarding negative thermal resistances are violated
- // apply correction procedure
- // Section (1.5.5) in FEFLOW White Papers Vol V.
- double constraint = 1.0 / ((1.0 / _R_gg) + (1.0 / (2.0 * _R_gs)));
+ 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::pair<double, double> thermalResistancesGroutSoil(double 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.
+
+ auto constraint = [&]() {
+ return 1.0 / ((1.0 / R_gg) + (1.0 / (2.0 * R_gs)));
+ };
+
int count = 0;
- while (constraint < 0.0)
+ while (constraint() < 0.0)
{
- if (extern_def_thermal_resistances.if_use_defined_therm_resis ||
- extern_Ra_Rb.use_extern_Ra_Rb)
- {
- OGS_FATAL(
- "Error!!! Constraints on thermal resistances are violated! "
- "Correction procedure can't be applied due to user defined "
- "thermal resistances! The simulation will be stopped!");
- }
if (count == 0)
{
chi *= 0.66;
- _R_gs = (1 - chi) * _R_g;
- _R_gg = 2.0 * _R_gs * (R_ar - 2.0 * chi * _R_g) /
- (2.0 * _R_gs - R_ar + 2.0 * chi * _R_g);
+ 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 = (1 - chi) * _R_g;
- _R_gg = 2.0 * _R_gs * (R_ar - 2.0 * chi * _R_g) /
- (2.0 * _R_gs - R_ar + 2.0 * chi * _R_g);
+ 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 = (1 - chi) * _R_g;
- _R_gg = 2.0 * _R_gs * (R_ar - 2.0 * chi * _R_g) /
- (2.0 * _R_gs - R_ar + 2.0 * chi * _R_g);
+ 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);
- constraint = 1.0 / ((1.0 / _R_gg) + (1.0 / (2.0 * _R_gs)));
count++;
}
-
- // kleep 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";
- // -------------------------------------------------------------------------
+ return {R_gg, R_gs};
}
+} // namespace
-/**
- * calculate heat transfer coefficient
- */
-void BHE_1U::calcHeatTransferCoefficients()
-{
- _PHI_fig = 1.0 / _R_fig;
- _PHI_fog = 1.0 / _R_fog;
- _PHI_gg = 1.0 / _R_gg;
- _PHI_gs = 1.0 / _R_gs;
-}
-
-double BHE_1U::getMassCoeff(std::size_t idx_unknown) const
-{
- double const& rho_r = refrigerant_param.rho_r;
- double const& heat_cap_r = refrigerant_param.heat_cap_r;
- double const& rho_g = grout_param.rho_g;
- double const& porosity_g = grout_param.porosity_g;
- double const& heat_cap_g = grout_param.heat_cap_g;
+constexpr std::pair<int, int> BHE_1U::inflow_outflow_bc_component_ids[];
- double mass_coeff = 0.0;
+void BHE_1U::updateHeatTransferCoefficients(double const flow_rate)
- switch (idx_unknown)
- {
- case 0: // i1
- mass_coeff = rho_r * heat_cap_r * CSA_i;
- break;
- case 1: // o1
- mass_coeff = rho_r * heat_cap_r * CSA_o;
- break;
- case 2: // g1
- mass_coeff = (1.0 - porosity_g) * rho_g * heat_cap_g * CSA_g1;
- break;
- case 3: // g2
- mass_coeff = (1.0 - porosity_g) * rho_g * heat_cap_g * CSA_g2;
- break;
- default:
- break;
- }
-
- return mass_coeff;
-}
-
-void BHE_1U::getLaplaceMatrix(std::size_t idx_unknown,
- Eigen::MatrixXd& mat_laplace) const
{
- double const& lambda_r = refrigerant_param.lambda_r;
- double const& rho_r = refrigerant_param.rho_r;
- double const& heat_cap_r = refrigerant_param.heat_cap_r;
- double const& alpha_L = refrigerant_param.alpha_L;
- double const& porosity_g = grout_param.porosity_g;
- double const& lambda_g = grout_param.lambda_g;
+ _flow_velocity = flow_rate / _pipes.inlet.area();
- // Here we calculates 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.
- double laplace_coeff(0.0);
- mat_laplace.setZero();
+ double const Re = reynoldsNumber(std::abs(_flow_velocity),
+ _pipes.inlet.diameter,
+ refrigerant.dynamic_viscosity,
+ refrigerant.density);
+ double const Pr = prandtlNumber(refrigerant.dynamic_viscosity,
+ refrigerant.specific_heat_capacity,
+ refrigerant.thermal_conductivity);
- switch (idx_unknown)
- {
- case 0:
- // pipe i1, Eq. 19
- laplace_coeff =
- (lambda_r + rho_r * heat_cap_r * alpha_L * _u.norm()) * CSA_i;
- break;
- case 1:
- // pipe o1, Eq. 20
- laplace_coeff =
- (lambda_r + rho_r * heat_cap_r * alpha_L * _u.norm()) * CSA_o;
- break;
- case 2:
- // pipe g1, Eq. 21
- laplace_coeff = (1.0 - porosity_g) * lambda_g * CSA_g1;
- break;
- case 3:
- // pipe g2, Eq. 22
- laplace_coeff = (1.0 - porosity_g) * lambda_g * CSA_g2;
- break;
- default:
- OGS_FATAL(
- "Error !!! The index passed to get_laplace_coeff for BHE is "
- "not correct. ");
- break;
- }
+ double const Nu =
+ nusseltNumber(Re, Pr, _pipes.inlet.diameter, borehole_geometry.length);
- mat_laplace(0, 0) = laplace_coeff;
- mat_laplace(1, 1) = laplace_coeff;
- mat_laplace(2, 2) = laplace_coeff;
+ _thermal_resistances = calcThermalResistances(Nu);
}
-void BHE_1U::getAdvectionVector(std::size_t idx_unknown,
- Eigen::VectorXd& vec_advection) const
+/// Nu is the Nusselt number.
+std::array<double, BHE_1U::number_of_unknowns> BHE_1U::calcThermalResistances(
+ double const Nu)
{
- double const& rho_r = refrigerant_param.rho_r;
- double const& heat_cap_r = refrigerant_param.heat_cap_r;
- double advection_coeff(0);
- vec_advection.setZero();
-
- switch (idx_unknown)
- {
- case 0:
- // pipe i1, Eq. 19
- advection_coeff = rho_r * heat_cap_r * _u(0) * CSA_i;
- // z direction
- vec_advection(2) = -1.0 * advection_coeff;
- break;
- case 1:
- // pipe o1, Eq. 20
- advection_coeff = rho_r * heat_cap_r * _u(0) * CSA_o;
- // z direction
- vec_advection(2) = 1.0 * advection_coeff;
- break;
- case 2:
- // grout g1, Eq. 21
- advection_coeff = 0.0;
- break;
- case 3:
- // grout g2, Eq. 22
- advection_coeff = 0.0;
- break;
- default:
- OGS_FATAL(
- "Error !!! The index passed to get_advection_coeff for BHE is "
- "not correct. ");
- break;
- }
-}
+ static constexpr double pi = boost::math::constants::pi<double>();
-void ProcessLib::HeatTransportBHE::BHE::BHE_1U::setRMatrices(
- const int idx_bhe_unknowns, const int NumNodes,
- Eigen::MatrixXd& matBHE_loc_R,
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
- R_matrix,
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
- R_pi_s_matrix,
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
- R_s_matrix) const
-{
- switch (idx_bhe_unknowns)
- {
- case 0: // PHI_fig
- R_matrix.block(0, 2 * NumNodes, NumNodes, NumNodes) +=
- -1.0 * matBHE_loc_R;
- R_matrix.block(2 * NumNodes, 0, NumNodes, NumNodes) +=
- -1.0 * matBHE_loc_R;
+ double const& lambda_r = refrigerant.thermal_conductivity;
+ double const& lambda_g = grout.lambda_g;
+ double const& lambda_p = _pipes.inlet.wall_thermal_conductivity;
- R_matrix.block(0, 0, NumNodes, NumNodes) +=
- 1.0 * matBHE_loc_R; // K_i1
- R_matrix.block(2 * NumNodes, 2 * NumNodes, NumNodes, NumNodes) +=
- 1.0 * matBHE_loc_R; // K_ig
- break;
- case 1: // PHI_fog
- R_matrix.block(NumNodes, 3 * NumNodes, NumNodes, NumNodes) +=
- -1.0 * matBHE_loc_R;
- R_matrix.block(3 * NumNodes, NumNodes, NumNodes, NumNodes) +=
- -1.0 * matBHE_loc_R;
+ // thermal resistances due to advective flow of refrigerant in the _pipes
+ // Eq. 36 in Diersch_2011_CG
+ double const R_adv_i1 = 1.0 / (Nu * lambda_r * pi);
+ double const R_adv_o1 = 1.0 / (Nu * lambda_r * pi);
- R_matrix.block(NumNodes, NumNodes, NumNodes, NumNodes) +=
- 1.0 * matBHE_loc_R; // K_o1
- R_matrix.block(3 * NumNodes, 3 * NumNodes, NumNodes, NumNodes) +=
- 1.0 * matBHE_loc_R; // K_og
- break;
- case 2: // PHI_gg
- R_matrix.block(2 * NumNodes, 3 * NumNodes, NumNodes, NumNodes) +=
- -1.0 * matBHE_loc_R;
- R_matrix.block(3 * NumNodes, 2 * NumNodes, NumNodes, NumNodes) +=
- -1.0 * matBHE_loc_R;
+ // thermal resistance due to thermal conductivity of the pipe wall material
+ // Eq. 49
+ double const R_con_a =
+ std::log(_pipes.outlet.diameter / _pipes.inlet.diameter) /
+ (2.0 * pi * lambda_p);
- R_matrix.block(2 * NumNodes, 2 * NumNodes, NumNodes, NumNodes) +=
- 1.0 * matBHE_loc_R; // K_ig // notice we only have
- // 1 PHI_gg term here.
- R_matrix.block(3 * NumNodes, 3 * NumNodes, NumNodes, NumNodes) +=
- 1.0 * matBHE_loc_R; // K_og // see Diersch 2013 FEFLOW
- // book page 954 Table M.2
- break;
- case 3: // PHI_gs
- R_s_matrix += 1.0 * matBHE_loc_R;
+ // the average outer diameter of the _pipes
+ 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);
+ // Eq. 52
+ // thermal resistances of the grout
+ double const R_g =
+ 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;
- R_pi_s_matrix.block(2 * NumNodes, 0, NumNodes, NumNodes) +=
- -1.0 * matBHE_loc_R;
- R_pi_s_matrix.block(3 * NumNodes, 0, NumNodes, NumNodes) +=
- -1.0 * matBHE_loc_R;
- R_matrix.block(2 * NumNodes, 2 * NumNodes, NumNodes, NumNodes) +=
- 1.0 * matBHE_loc_R; // K_ig
- R_matrix.block(3 * NumNodes, 3 * NumNodes, NumNodes, NumNodes) +=
- 1.0 * matBHE_loc_R; // K_og
- break;
- }
-}
+ // thermal resistance due to inter-grout exchange
+ double const R_ar =
+ std::acosh((2.0 * _pipes.distance * _pipes.distance - d0 * d0) / d0 /
+ d0) /
+ (2.0 * pi * lambda_g);
-double BHE_1U::getBoundaryHeatExchangeCoeff(std::size_t idx_unknown) const
-{
- // Here we calculates the boundary 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.
+ double R_gg, R_gs;
+ std::tie(R_gg, R_gs) = thermalResistancesGroutSoil(chi, R_ar, R_g);
- double exchange_coeff(0);
+ return {R_fig, R_fog, R_gg, R_gs};
- switch (idx_unknown)
- {
- case 0:
- // PHI_fig
- exchange_coeff = _PHI_fig;
- break;
- case 1:
- // PHI_fog
- exchange_coeff = _PHI_fog;
- break;
- case 2:
- // PHI_gg
- exchange_coeff = _PHI_gg;
- break;
- case 3:
- // PHI_gs
- exchange_coeff = _PHI_gs;
- break;
- default:
- OGS_FATAL(
- "Error !!! The index passed to "
- "getBoundaryHeatExchangeCoeff for BHE is not correct. ");
- break;
- }
- return exchange_coeff;
+ // 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_1U::getTinByTout(double T_out, double current_time = -1.0)
+double BHE_1U::getTinByTout(double const T_out, double const current_time)
{
- double T_in(0.0);
- double power_tmp(0.0);
- double building_power_tmp(0.0);
- // double power_elect_tmp(0.0);
- double Q_r_tmp(0.0);
- double COP_tmp(0.0);
- double fac_dT = 1.0;
- double const& rho_r = refrigerant_param.rho_r;
- double const& heat_cap_r = refrigerant_param.heat_cap_r;
-
- switch (this->boundary_type)
- {
- case BHE_BOUNDARY_TYPE::FIXED_INFLOW_TEMP_CURVE_BOUNDARY:
- {
- T_in = inflow_temperature_curve->getValue(current_time);
- }
- break;
- case BHE_BOUNDARY_TYPE::POWER_IN_WATT_BOUNDARY:
- if (use_flowrate_curve)
- {
- Q_r_tmp = flowrate_curve->getValue(current_time);
-
- updateFlowRate(Q_r_tmp);
- }
- else
- Q_r_tmp = Q_r;
- T_in = power_in_watt_val / Q_r_tmp / heat_cap_r / rho_r + T_out;
- break;
- case BHE_BOUNDARY_TYPE::FIXED_TEMP_DIFF_BOUNDARY:
- if (use_flowrate_curve)
- {
- Q_r_tmp = flowrate_curve->getValue(current_time);
-
- updateFlowRate(Q_r_tmp);
- }
- T_in = T_out + delta_T_val;
- break;
- case BHE_BOUNDARY_TYPE::POWER_IN_WATT_CURVE_FIXED_DT_BOUNDARY:
- // get the power value in the curve
- // power_tmp = GetCurveValue(power_in_watt_curve_idx, 0,
- // current_time, &flag_valid);
- power_tmp = power_in_watt_curve->getValue(current_time);
-
- if (power_tmp < 0)
- fac_dT = -1.0;
- else
- fac_dT = 1.0;
- // if power value exceeds threshold, calculate new values
- if (fabs(power_tmp) > threshold)
- {
- // calculate the corresponding flow rate needed
- // using the defined delta_T value
- Q_r_tmp =
- power_tmp / (fac_dT * delta_T_val) / heat_cap_r / rho_r;
- // update all values dependent on the flow rate
- updateFlowRate(Q_r_tmp);
- // calculate the new T_in
- T_in = T_out + (fac_dT * delta_T_val);
- // print out updated flow rate
- // std::cout << "Qr: " << Q_r_tmp << std::endl;
- }
- else
- {
- Q_r_tmp = 1.0e-12; // this has to be a small value to avoid
- // division by zero
- // update all values dependent on the flow rate
- updateFlowRate(Q_r_tmp);
- // calculate the new T_in
- T_in = T_out;
- // print out updated flow rate
- // std::cout << "Qr: " << Q_r_tmp << std::endl;
- }
- break;
- case BHE_BOUNDARY_TYPE::BUILDING_POWER_IN_WATT_CURVE_FIXED_DT_BOUNDARY:
- // get the building power value in the curve
- // building_power_tmp = GetCurveValue(power_in_watt_curve_idx, 0,
- // current_time, &flag_valid);
- building_power_tmp = power_in_watt_curve->getValue(current_time);
-
- if (building_power_tmp <= 0.0)
- {
- // get COP value based on T_out in the curve
- COP_tmp = heating_cop_curve->getValue(T_out);
-
- // now calculate how much power needed from BHE
- power_tmp = building_power_tmp * (COP_tmp - 1.0) / COP_tmp;
- // also how much power from electricity
- // power_elect_tmp = building_power_tmp - power_tmp;
- // print the amount of power needed
- // std::cout << "COP: " << COP_tmp << ", Q_bhe: " << power_tmp
- // << ", Q_elect: " << power_elect_tmp << std::endl;
- fac_dT = -1.0;
- }
- else
- {
- // get COP value based on T_out in the curve
- COP_tmp = cooling_cop_curve->getValue(T_out);
-
- // now calculate how much power needed from BHE
- power_tmp = building_power_tmp * (COP_tmp + 1.0) / COP_tmp;
- // also how much power from electricity
- // power_elect_tmp = -building_power_tmp + power_tmp;
- // print the amount of power needed
- // std::cout << "COP: " << COP_tmp << ", Q_bhe: " << power_tmp
- // << ", Q_elect: " << power_elect_tmp << std::endl;
- fac_dT = 1.0;
- }
- // if power value exceeds threshold, calculate new values
- if (fabs(power_tmp) > threshold)
- {
- // calculate the corresponding flow rate needed
- // using the defined delta_T value
- Q_r_tmp =
- power_tmp / (fac_dT * delta_T_val) / heat_cap_r / rho_r;
- // update all values dependent on the flow rate
- updateFlowRate(Q_r_tmp);
- // calculate the new T_in
- T_in = T_out + (fac_dT * delta_T_val);
- // print out updated flow rate
- // std::cout << "Qr: " << Q_r_tmp << std::endl;
- }
- else
- {
- Q_r_tmp = 1.0e-12; // this has to be a small value to avoid
- // division by zero
- // update all values dependent on the flow rate
- updateFlowRate(Q_r_tmp);
- // calculate the new T_in
- T_in = T_out;
- // print out updated flow rate
- // std::cout << "Qr: " << Q_r_tmp << std::endl;
- }
- break;
- case BHE_BOUNDARY_TYPE::
- BUILDING_POWER_IN_WATT_CURVE_FIXED_FLOW_RATE_BOUNDARY:
- // get the building power value in the curve
- // building_power_tmp = GetCurveValue(power_in_watt_curve_idx, 0,
- // current_time, &flag_valid);
- building_power_tmp = power_in_watt_curve->getValue(current_time);
-
- if (building_power_tmp <= 0)
- {
- // get COP value based on T_out in the curve
- COP_tmp = heating_cop_curve->getValue(T_out);
- // now calculate how much power needed from BHE
- power_tmp = building_power_tmp * (COP_tmp - 1.0) / COP_tmp;
- // also how much power from electricity
- // power_elect_tmp = building_power_tmp - power_tmp;
- // print the amount of power needed
- // std::cout << "COP: " << COP_tmp << ", Q_bhe: " << power_tmp
- // << ", Q_elect: " << power_elect_tmp << std::endl;
- }
- else
- {
- // get COP value based on T_out in the curve
- COP_tmp = cooling_cop_curve->getValue(T_out);
- // now calculate how much power needed from BHE
- power_tmp = building_power_tmp * (COP_tmp + 1.0) / COP_tmp;
- // also how much power from electricity
- // power_elect_tmp = -building_power_tmp + power_tmp;
- // print the amount of power needed
- // std::cout << "COP: " << COP_tmp << ", Q_bhe: " << power_tmp
- // << ", Q_elect: " << power_elect_tmp << std::endl;
- }
- // Assign Qr whether from curve or fixed value
- if (use_flowrate_curve)
- {
- // Q_r_tmp = GetCurveValue(flowrate_curve_idx, 0, current_time,
- // &flag_valid);
- Q_r_tmp = flowrate_curve->getValue(current_time);
- updateFlowRate(Q_r_tmp);
- }
- else
- Q_r_tmp = Q_r;
- if (fabs(power_tmp) < threshold)
- {
- Q_r_tmp = 1.0e-12; // this has to be a small value to avoid
- // division by zero update all values
- // dependent on the flow rate
- updateFlowRate(Q_r_tmp);
- // calculate the new T_in
- T_in = T_out;
- // print out updated flow rate
- // std::cout << "Qr: " << Q_r_tmp << std::endl;
- }
- else
- {
- Q_r_tmp = Q_r;
- updateFlowRate(Q_r_tmp);
- // calculate the dT value based on fixed flow rate
- delta_T_val = power_tmp / Q_r_tmp / heat_cap_r / rho_r;
- // calcuate the new T_in
- T_in = T_out + delta_T_val;
- }
- break;
- case BHE_BOUNDARY_TYPE::POWER_IN_WATT_CURVE_FIXED_FLOW_RATE_BOUNDARY:
- // get the power value in the curve
- // power_tmp = GetCurveValue(power_in_watt_curve_idx, 0,
- // current_time, &flag_valid);
- power_tmp = power_in_watt_curve->getValue(current_time);
-
- // Assign Qr whether from curve or fixed value
- if (use_flowrate_curve)
- {
- // Q_r_tmp = GetCurveValue(flowrate_curve_idx, 0, current_time,
- // &flag_valid);
- Q_r_tmp = flowrate_curve->getValue(current_time);
- updateFlowRate(Q_r_tmp);
- }
- else
- Q_r_tmp = Q_r;
- // calculate the dT value based on fixed flow rate
- if (fabs(power_tmp) < threshold)
- {
- Q_r_tmp = 1.0e-12; // this has to be a small value to avoid
- // division by zero update all values
- // dependent on the flow rate
- updateFlowRate(Q_r_tmp);
- // calculate the new T_in
- T_in = T_out;
- // print out updated flow rate
- // std::cout << "Qr: " << Q_r_tmp << std::endl;
- }
- else
- {
- Q_r_tmp = Q_r;
- updateFlowRate(Q_r_tmp);
- // calculate the dT value based on fixed flow rate
- delta_T_val = power_tmp / Q_r_tmp / heat_cap_r / rho_r;
- // calcuate the new T_in
- T_in = T_out + delta_T_val;
- }
- break;
- default:
- T_in = T_out;
- break;
- }
-
- return T_in;
+ auto values = apply_visitor(
+ [&](auto const& control) { return control(T_out, current_time); },
+ flowAndTemperatureControl);
+ updateHeatTransferCoefficients(values.flow_rate);
+ return values.temperature;
}
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
index ec9620a6b98..1aaa0c91eb3 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
@@ -9,220 +9,162 @@
#pragma once
-#include "BHEAbstract.h"
+#include <Eigen/Eigen>
+
+#include "BaseLib/Error.h"
+
+#include "BHECommon.h"
+#include "FlowAndTemperatureControl.h"
+#include "PipeConfiguration1U.h"
namespace ProcessLib
{
namespace HeatTransportBHE
{
-namespace BHE // namespace of borehole heat exchanger
+namespace BHE
{
-class BHE_1U final : public BHEAbstract
+class BHE_1U final : public BHECommon
{
public:
- /**
- * constructor
- */
- BHE_1U(
- const std::string name /* name of the BHE */,
- BHE::BHE_BOUNDARY_TYPE bound_type /* type of BHE boundary */,
- std::map<std::string,
- std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
- bhe_curves /* bhe related curves */,
- BoreholeGeometry borehole_geometry = {100, 0.013},
- PipeParameters pipe_geometry =
- {0.016 /* inner radius of the pipline */,
- 0.016 /* outer radius of the pipline */,
- 0.0029 /* pipe-in wall thickness*/,
- 0.0029 /* pipe-out wall thickness*/,
- 0.38 /* thermal conductivity of the pipe wall */,
- 0.38 /* thermal conductivity of the inner pipe wall */,
- 0.38 /* thermal conductivity of the outer pipe wall */},
- RefrigerantParameters refrigerant_param =
- {
- 0.00054741 /* dynamic viscosity of the refrigerant */,
- 988.1 /* density of the refrigerant */,
- 0.6405 /* thermal conductivity of the refrigerant */,
- 4180 /* specific heat capacity of the refrigerant */, 1.0e-4 /* longitudinal dispersivity of the refrigerant in the pipeline */},
- GroutParameters grout_param =
- {2190 /* density of the grout */, 0.5 /* porosity of the grout */,
- 1000 /* specific heat capacity of the grout */,
- 2.3 /* thermal conductivity of the grout */},
- ExternallyDefinedRaRb extern_Ra_Rb =
- {false /* whether Ra and Rb values are used */,
- 0.0 /* external defined borehole internal thermal resistance */,
- 0.0 /* external defined borehole thermal resistance */},
- ExternallyDefinedThermalResistances extern_def_thermal_resistances =
- {false /* whether user defined R values are used */,
- 0.0 /* external defined borehole thermal resistance */,
- 0.0 /* external defined borehole thermal resistance */,
- 0.0 /* external defined borehole thermal resistance */,
- 0.0 /* external defined borehole thermal resistance */,
- 0.0 /* external defined borehole thermal resistance */},
- double my_Qr = 21.86 /
- 86400 /* total refrigerant flow discharge of BHE */,
- double my_omega = 0.06 /* pipe distance */,
- double my_power_in_watt = 0.0 /* injected or extracted power */,
- double my_delta_T_val =
- 0.0 /* Temperature difference btw inflow and outflow temperature */,
-
- bool if_flowrate_curve = false /* whether flowrate curve is used*/,
- double my_threshold = 0.1) /* Threshold Q value for switching off the
- BHE when using Q_Curve_fixed_dT B.C.*/
- : BHEAbstract(name, BHE_TYPE::TYPE_1U, borehole_geometry, pipe_geometry,
- refrigerant_param, grout_param, extern_Ra_Rb,
- extern_def_thermal_resistances, std::move(bhe_curves),
- bound_type, if_flowrate_curve)
+ BHE_1U(BoreholeGeometry const& borehole,
+ RefrigerantProperties const& refrigerant,
+ GroutParameters const& grout,
+ FlowAndTemperatureControl const& flowAndTemperatureControl,
+ PipeConfiguration1U const& pipes)
+ : BHECommon{borehole, refrigerant, grout, flowAndTemperatureControl},
+ _pipes(pipes)
{
- _u = Eigen::Vector2d::Zero();
- _Nu = Eigen::Vector2d::Zero();
- // 1U type of BHE has 2 pipes
-
- Q_r = my_Qr;
-
- omega = my_omega;
- power_in_watt_val = my_power_in_watt;
- delta_T_val = my_delta_T_val;
- threshold = my_threshold;
-
- // get the corresponding curve
- std::map<std::string,
- std::unique_ptr<MathLib::PiecewiseLinearInterpolation>>::
- const_iterator it;
- if (bound_type ==
- BHE_BOUNDARY_TYPE::POWER_IN_WATT_CURVE_FIXED_DT_BOUNDARY ||
- bound_type == BHE_BOUNDARY_TYPE::
- POWER_IN_WATT_CURVE_FIXED_FLOW_RATE_BOUNDARY ||
- bound_type ==
- BHE_BOUNDARY_TYPE::
- BUILDING_POWER_IN_WATT_CURVE_FIXED_FLOW_RATE_BOUNDARY)
- {
- it = bhe_curves.find("power_in_watt_curve");
- if (it == bhe_curves.end())
- {
- // curve not found, fatal error
- OGS_FATAL(
- "Required pow_in_watt_curve cannot be found in the BHE "
- "parameters!");
- }
-
- // curve successfully found
- power_in_watt_curve = it->second.get();
- }
-
- if (if_flowrate_curve)
- {
- use_flowrate_curve = true;
-
- it = bhe_curves.find("flow_rate_curve");
- if (it == bhe_curves.end())
- {
- OGS_FATAL(
- "Required flow_rate_curve annot be found in the BHE "
- "parameters!");
- }
-
- // curve successfully found
- flowrate_curve = it->second.get();
- }
- if (bound_type == BHE_BOUNDARY_TYPE::FIXED_INFLOW_TEMP_CURVE_BOUNDARY)
- {
- it = bhe_curves.find("inflow_temp_curve");
- if (it == bhe_curves.end())
- {
- OGS_FATAL(
- "Required inflow_temp_curve annot be found in the BHE "
- "parameters!");
- }
-
- // curve successfully found
- inflow_temperature_curve = it->second.get();
- }
-
- // Table 1 in Diersch_2011_CG
- S_i = PI * 2.0 * pipe_geometry.r_inner;
- S_o = PI * 2.0 * pipe_geometry.r_inner;
- S_g1 = borehole_geometry.diameter;
- S_gs =
- 0.5 * PI *
- borehole_geometry.diameter; // Corrected value according to FEFLOW
- // White Papers Vol V, Table 1-71
+ // Initialize thermal resistances.
+ auto values = apply_visitor(
+ [&](auto const& control) {
+ return control(refrigerant.reference_temperature,
+ 0. /* initial time */);
+ },
+ flowAndTemperatureControl);
+ updateHeatTransferCoefficients(values.flow_rate);
+ }
- // cross section area calculation
- CSA_i = CSA_o = PI * pipe_geometry.r_inner * pipe_geometry.r_inner;
- CSA_g1 = CSA_g2 = PI * (0.125 * borehole_geometry.diameter *
- borehole_geometry.diameter -
- pipe_geometry.r_outer * pipe_geometry.r_outer);
+ static constexpr int number_of_unknowns = 4;
- // initialization calculation
- initialize();
- };
+ std::array<double, number_of_unknowns> 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,
+ /*o1*/ 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}};
+ }
- /**
- * return the number of unknowns needed for 1U BHE
- */
- std::size_t getNumUnknowns() const { return 4; }
+ std::array<double, number_of_unknowns> 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_despersion_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);
+
+ // 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, Eq. 19
+ (lambda_r + rho_r * Cp_r * alpha_L * velocity_norm),
+ // pipe o1, Eq. 20
+ (lambda_r + rho_r * Cp_r * alpha_L * velocity_norm),
+ // pipe g1, Eq. 21
+ (1.0 - porosity_g) * lambda_g,
+ // pipe g2, Eq. 22
+ (1.0 - porosity_g) * lambda_g}};
+ }
- void initialize();
+ std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors() const
+ {
+ double const& rho_r = refrigerant.density;
+ double const& Cp_r = refrigerant.specific_heat_capacity;
+
+ return {{// pipe i1, Eq. 19
+ {0, 0, -rho_r * Cp_r * _flow_velocity},
+ // pipe o1, Eq. 20
+ {0, 0, rho_r * Cp_r * _flow_velocity},
+ // grout g1, Eq. 21
+ {0, 0, 0},
+ // grout g2, Eq. 22
+ {0, 0, 0}}};
+ }
- void updateFlowRateFromCurve(double current_time)
+ 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
{
- if (use_flowrate_curve)
+ switch (idx_bhe_unknowns)
{
- double Q_r_tmp(0.0);
- Q_r_tmp = flowrate_curve->getValue(current_time);
- updateFlowRate(Q_r_tmp);
+ case 0: // PHI_fig
+ R_matrix.block(0, 2 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(2 * NPoints, 0, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+
+ R_matrix.block(0, 0, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R; // K_i1
+ R_matrix.block(2 * NPoints, 2 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R; // K_ig
+ break;
+ case 1: // PHI_fog
+ R_matrix.block(NPoints, 3 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(3 * NPoints, NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+
+ R_matrix.block(NPoints, 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_og
+ break;
+ case 2: // PHI_gg
+ R_matrix.block(2 * NPoints, 3 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(3 * NPoints, 2 * NPoints, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+
+ R_matrix.block(2 * NPoints, 2 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R; // K_ig // notice we only have
+ // 1 PHI_gg term here.
+ R_matrix.block(3 * NPoints, 3 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R; // K_og // see Diersch 2013 FEFLOW
+ // book page 954 Table M.2
+ break;
+ case 3: // PHI_gs
+ R_s_matrix.template block<NPoints, NPoints>(0, 0).noalias() +=
+ 1.0 * matBHE_loc_R;
+
+ R_pi_s_matrix.block(2 * NPoints, 0, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_pi_s_matrix.block(3 * NPoints, 0, NPoints, NPoints) +=
+ -1.0 * matBHE_loc_R;
+ R_matrix.block(2 * NPoints, 2 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R; // K_ig
+ R_matrix.block(3 * NPoints, 3 * NPoints, NPoints, NPoints) +=
+ 1.0 * matBHE_loc_R; // K_og
+ break;
}
- };
-
- /**
- * calculate thermal resistance
- */
- void calcThermalResistances();
-
- /**
- * calculate heat transfer coefficient
- */
- void calcHeatTransferCoefficients();
-
- /**
- * return the coeff of mass matrix,
- * depending on the index of unknown.
- */
- double getMassCoeff(std::size_t idx_unknown) const;
-
- /**
- * return the coeff of laplace matrix,
- * depending on the index of unknown.
- */
- void getLaplaceMatrix(std::size_t idx_unknown,
- Eigen::MatrixXd& mat_laplace) const;
-
- /**
- * return the coeff of advection matrix,
- * depending on the index of unknown.
- */
- void getAdvectionVector(std::size_t idx_unknown,
- Eigen::VectorXd& vec_advection) const;
-
- /**
- * return the _R_matrix, _R_pi_s_matrix, _R_s_matrix
- */
- void setRMatrices(
- const int idx_bhe_unknowns, const int NumNodes,
- Eigen::MatrixXd& matBHE_loc_R,
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
- R_matrix,
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
- R_pi_s_matrix,
- Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
- R_s_matrix) const;
-
- /**
- * return the coeff of boundary heat exchange matrix,
- * depending on the index of unknown.
- */
- double getBoundaryHeatExchangeCoeff(std::size_t idx_unknown) const;
+ }
/**
* return the inflow temperature based on outflow temperature and fixed
@@ -230,79 +172,41 @@ public:
*/
double getTinByTout(double T_out, double current_time);
- std::vector<std::pair<int, int>> const& inflowOutflowBcComponentIds()
- const override
+ double thermalResistance(int const unknown_index) const
{
- return _inflow_outflow_bc_component_ids;
+ return _thermal_resistances[unknown_index];
}
- /**
- * required by eigen library,
- * to make sure the dynamically allocated class has
- * aligned "operator new"
- */
- EIGEN_MAKE_ALIGNED_OPERATOR_NEW
-
-private:
- std::vector<std::pair<int, int>> const _inflow_outflow_bc_component_ids = {
+ static constexpr std::pair<int, int> inflow_outflow_bc_component_ids[] = {
{0, 1}};
- /**
- * thermal resistances
- */
- double _R_fig, _R_fog;
-
- /**
- * thermal resistances due to advective flow of refrigerant in the pipes
- */
- double _R_adv_i1, _R_adv_o1;
-
- /**
- * thermal resistances due to the pipe wall material
- */
- double _R_con_a_i1, _R_con_a_o1;
-
- /**
- * thermal resistances due to the grout transition
- */
- double _R_con_b;
+private:
+ // Placing it here before using it in the cross_section_areas.
+ PipeConfiguration1U const _pipes;
- /**
- * thermal resistances of the grout
- */
- double _R_g;
+public:
+ std::array<double, number_of_unknowns> const cross_section_areas = {
+ {_pipes.inlet.area(), _pipes.inlet.area(),
+ borehole_geometry.area() / 2 - _pipes.outlet.area(),
+ borehole_geometry.area() / 2 - _pipes.outlet.area()}};
- /**
- * thermal resistances of the grout soil exchange
- */
- double _R_gs;
+private:
+ void updateHeatTransferCoefficients(double const flow_rate);
- /**
- * thermal resistances due to inter-grout exchange
- */
- double _R_gg;
+ std::array<double, number_of_unknowns> calcThermalResistances(
+ double const Nu);
- /**
- * heat transfer coefficients
- */
- double _PHI_fig, _PHI_fog, _PHI_gg, _PHI_gs;
-
- /**
- * specific exchange surfaces S
- */
- double S_i, S_o, S_g1, S_gs;
- /**
- * cross section area
- */
- double CSA_i, CSA_o, CSA_g1, CSA_g2;
- /**
- * Nusselt number
- */
- Eigen::Vector2d _Nu;
- /**
- * flow velocity inside the pipeline
- */
- Eigen::Vector2d _u;
+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;
+
+ /// Flow velocity inside the pipes. Depends on the flow_rate.
+ double _flow_velocity;
};
} // namespace BHE
} // namespace HeatTransportBHE
diff --git a/ProcessLib/HeatTransportBHE/BHE/BoreholeGeometry.h b/ProcessLib/HeatTransportBHE/BHE/BoreholeGeometry.h
index 658aa9816a3..520aec45472 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BoreholeGeometry.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BoreholeGeometry.h
@@ -11,6 +11,8 @@
#pragma once
+#include <boost/math/constants/constants.hpp>
+
namespace ProcessLib
{
namespace HeatTransportBHE
@@ -30,6 +32,12 @@ struct BoreholeGeometry
* unit is m
*/
double const diameter;
+
+ double area() const
+ {
+ constexpr double pi = boost::math::constants::pi<double>();
+ return pi * diameter * diameter / 4;
+ }
};
} // namespace BHE
} // namespace HeatTransportBHE
diff --git a/ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.cpp b/ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.cpp
index 1ec84499d05..dcdb87183c7 100644
--- a/ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.cpp
@@ -1,4 +1,6 @@
/**
+ * \file
+ *
* \copyright
* Copyright (c) 2012-2018, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
@@ -8,242 +10,75 @@
*/
#include "CreateBHE1U.h"
-#include "MaterialLib/Fluid/Density/CreateFluidDensityModel.h"
-#include "MaterialLib/Fluid/FluidProperty.h"
-#include "MaterialLib/Fluid/SpecificHeatCapacity/CreateSpecificFluidHeatCapacityModel.h"
-#include "MaterialLib/Fluid/ThermalConductivity/CreateFluidThermalConductivityModel.h"
-#include "MaterialLib/Fluid/Viscosity/CreateViscosityModel.h"
+#include "BaseLib/ConfigTree.h"
+#include "CreateFlowAndTemperatureControl.h"
namespace ProcessLib
{
namespace HeatTransportBHE
{
namespace BHE
{
-BHE::BHE_1U* CreateBHE1U(
- BaseLib::ConfigTree const& bhe_conf,
+BHE::BHE_1U createBHE1U(
+ BaseLib::ConfigTree const& config,
std::map<std::string,
std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
- curves,
- std::unique_ptr<MaterialLib::Fluid::FluidProperty> const&
- bhe_refrigerant_density,
- std::unique_ptr<MaterialLib::Fluid::FluidProperty> const&
- bhe_refrigerant_viscosity,
- std::unique_ptr<MaterialLib::Fluid::FluidProperty> const&
- bhe_refrigerant_heat_capacity,
- std::unique_ptr<MaterialLib::Fluid::FluidProperty> const&
- bhe_regrigerant_heat_conductivity)
+ curves)
{
- // read in the parameters
- const std::string bhe_ply_name =
- bhe_conf.getConfigParameter<std::string>("bhe_polyline");
- const std::string bhe_bound_type_str =
- bhe_conf.getConfigParameter<std::string>("bhe_bound_type");
- const bool bhe_if_use_flow_rate_curve = false;
- const double bhe_length = bhe_conf.getConfigParameter<double>("bhe_length");
- const double bhe_diameter =
- bhe_conf.getConfigParameter<double>("bhe_diameter");
- const double bhe_refrigerant_flow_rate =
- bhe_conf.getConfigParameter<double>("bhe_refrigerant_flow_rate");
- const double bhe_pipe_inner_radius =
- bhe_conf.getConfigParameter<double>("bhe_pipe_inner_radius");
- const double bhe_pipe_outer_radius =
- bhe_conf.getConfigParameter<double>("bhe_pipe_outer_radius");
- const double bhe_pipe_in_wall_thickness =
- bhe_conf.getConfigParameter<double>("bhe_pipe_in_wall_thickness");
- const double bhe_pipe_out_wall_thickness =
- bhe_conf.getConfigParameter<double>("bhe_pipe_out_wall_thickness");
- const double bhe_fluid_longitudinal_dispsion_length =
- bhe_conf.getConfigParameter<double>(
- "bhe_fluid_longitudinal_dispsion_length");
- const double bhe_grout_density =
- bhe_conf.getConfigParameter<double>("bhe_grout_density");
- const double bhe_grout_porosity =
- bhe_conf.getConfigParameter<double>("bhe_grout_porosity");
- const double bhe_grout_heat_capacity =
- bhe_conf.getConfigParameter<double>("bhe_grout_heat_capacity");
- const double bhe_pipe_wall_thermal_conductivity =
- bhe_conf.getConfigParameter<double>(
- "bhe_pipe_wall_thermal_conductivity");
- const double bhe_grout_thermal_conductivity =
- bhe_conf.getConfigParameter<double>("bhe_grout_thermal_conductivity");
- const double bhe_pipe_distance =
- bhe_conf.getConfigParameter<double>("bhe_pipe_distance");
-
- auto const bhe_use_ext_therm_resis_conf =
- bhe_conf.getConfigParameterOptional<bool>(
- "bhe_use_external_therm_resis");
-
- // optional parameters
- double bhe_power_in_watt_val = 0.0;
- double bhe_intern_resistance = 1.0;
- double bhe_therm_resistance = 1.0;
- double bhe_R_fig = 1.0;
- double bhe_R_fog = 1.0;
- double bhe_R_gg1 = 1.0;
- double bhe_R_gg2 = 1.0;
- double bhe_R_gs = 1.0;
- double bhe_delta_T_val = 0.0;
- double bhe_switch_off_threshold = 1.0e-6;
-
- // give default values to optional parameters
- // if the BHE is using external given thermal resistance values
- bool bhe_use_ext_therm_resis = false;
- if (*bhe_use_ext_therm_resis_conf == true)
- {
- DBUG("If using external given thermal resistance values : %s",
- (*bhe_use_ext_therm_resis_conf) ? "true" : "false");
- bhe_use_ext_therm_resis = *bhe_use_ext_therm_resis_conf;
-
- // only when using it, read the two resistance values
- if (bhe_use_ext_therm_resis)
- {
- bhe_intern_resistance =
- *bhe_conf.getConfigParameterOptional<double>(
- "bhe_internal_resistance");
-
- bhe_therm_resistance =
- bhe_conf
- .getConfigParameterOptional<double>("bhe_therm_resistance")
- .get();
- }
- }
-
- // if the BHE is using user defined thermal resistance values
- bool bhe_user_defined_therm_resis = false;
- if (auto const bhe_user_defined_therm_resis_conf =
- bhe_conf.getConfigParameterOptional<bool>(
- "bhe_user_defined_therm_resis"))
- {
- DBUG("If appplying user defined thermal resistance values : %s",
- (*bhe_user_defined_therm_resis_conf) ? "true" : "false");
- bhe_user_defined_therm_resis = *bhe_user_defined_therm_resis_conf;
+ auto const& borehole_config = config.getConfigSubtree("borehole");
+ const double borehole_length =
+ borehole_config.getConfigParameter<double>("length");
+ const double borehole_diameter =
+ borehole_config.getConfigParameter<double>("diameter");
+ BoreholeGeometry const borehole{borehole_length, borehole_diameter};
- // only when using it, read the two resistance values
- if (bhe_user_defined_therm_resis)
- {
- bhe_R_fig =
- bhe_conf.getConfigParameterOptional<double>("bhe_R_fig").get();
- bhe_R_fog =
- bhe_conf.getConfigParameterOptional<double>("bhe_R_fog").get();
- bhe_R_gg1 =
- bhe_conf.getConfigParameterOptional<double>("bhe_R_gg1").get();
- bhe_R_gg2 =
- bhe_conf.getConfigParameterOptional<double>("bhe_R_gg2").get();
- bhe_R_gs =
- bhe_conf.getConfigParameterOptional<double>("bhe_R_gs").get();
- }
- }
+ auto const& pipes_config = config.getConfigSubtree("pipes");
+ Pipe const inlet_pipe = createPipe(pipes_config.getConfigSubtree("inlet"));
+ Pipe const outlet_pipe =
+ createPipe(pipes_config.getConfigSubtree("outlet"));
+ const double pipe_distance =
+ pipes_config.getConfigParameter<double>("distance_between_pipes");
+ const double pipe_longitudinal_dispersion_length =
+ pipes_config.getConfigParameter<double>(
+ "longitudinal_dispersion_length");
+ PipeConfiguration1U const pipes{inlet_pipe, outlet_pipe, pipe_distance,
+ pipe_longitudinal_dispersion_length};
- // convert BHE boundary type
- BHE::BHE_BOUNDARY_TYPE bhe_bound_type;
- if (bhe_bound_type_str.compare("FIXED_INFLOW_TEMP") == 0)
- bhe_bound_type = BHE_BOUNDARY_TYPE::FIXED_INFLOW_TEMP_BOUNDARY;
- else if (bhe_bound_type_str.compare("FIXED_INFLOW_TEMP_CURVE") == 0)
- bhe_bound_type = BHE_BOUNDARY_TYPE::FIXED_INFLOW_TEMP_CURVE_BOUNDARY;
- else if (bhe_bound_type_str.compare("POWER_IN_WATT") == 0)
- {
- bhe_bound_type = BHE_BOUNDARY_TYPE::POWER_IN_WATT_BOUNDARY;
- bhe_power_in_watt_val =
- bhe_conf
- .getConfigParameterOptional<double>("bhe_power_in_watt_value")
- .get();
- bhe_switch_off_threshold =
- bhe_conf
- .getConfigParameterOptional<double>("bhe_switch_off_threshold")
- .get();
- }
- else if (bhe_bound_type_str.compare("POWER_IN_WATT_CURVE_FIXED_DT") == 0)
- {
- bhe_bound_type =
- BHE_BOUNDARY_TYPE::POWER_IN_WATT_CURVE_FIXED_DT_BOUNDARY;
- bhe_switch_off_threshold =
- bhe_conf
- .getConfigParameterOptional<double>("bhe_switch_off_threshold")
- .get();
- }
- else if (bhe_bound_type_str.compare(
- "BHE_BOUND_BUILDING_POWER_IN_WATT_CURVE_FIXED_DT") == 0)
- {
- bhe_bound_type =
- BHE_BOUNDARY_TYPE::BUILDING_POWER_IN_WATT_CURVE_FIXED_DT_BOUNDARY;
- bhe_switch_off_threshold =
- bhe_conf
- .getConfigParameterOptional<double>("bhe_switch_off_threshold")
- .get();
- }
- else if (bhe_bound_type_str.compare(
- "BHE_BOUND_BUILDING_POWER_IN_WATT_CURVE_FIXED_FLOW_RATE") == 0)
- {
- bhe_bound_type = BHE_BOUNDARY_TYPE::
- BUILDING_POWER_IN_WATT_CURVE_FIXED_FLOW_RATE_BOUNDARY;
- bhe_switch_off_threshold =
- bhe_conf
- .getConfigParameterOptional<double>("bhe_switch_off_threshold")
- .get();
- }
- else if (bhe_bound_type_str.compare(
- "POWER_IN_WATT_CURVE_FIXED_FLOW_RATE") == 0)
- {
- bhe_bound_type =
- BHE_BOUNDARY_TYPE::POWER_IN_WATT_CURVE_FIXED_FLOW_RATE_BOUNDARY;
- bhe_switch_off_threshold =
- bhe_conf
- .getConfigParameterOptional<double>("bhe_switch_off_threshold")
- .get();
- }
- else if (bhe_bound_type_str.compare("FIXED_TEMP_DIFF") == 0)
- {
- bhe_bound_type = BHE_BOUNDARY_TYPE::FIXED_TEMP_DIFF_BOUNDARY;
- bhe_delta_T_val = bhe_conf
- .getConfigParameterOptional<double>(
- "bhe_inout_delta_temperature")
- .get();
- }
- else
- {
- bhe_bound_type = BHE_BOUNDARY_TYPE::FIXED_INFLOW_TEMP_BOUNDARY;
- ERR("The BHE boundary condition is not correctly set! ");
- }
+ auto const& grout_config = config.getConfigSubtree("grout");
+ const double grout_density =
+ grout_config.getConfigParameter<double>("density");
+ const double grout_porosity =
+ grout_config.getConfigParameter<double>("porosity");
+ const double grout_heat_capacity =
+ grout_config.getConfigParameter<double>("heat_capacity");
+ const double grout_thermal_conductivity =
+ grout_config.getConfigParameter<double>("thermal_conductivity");
+ GroutParameters const grout{grout_density, grout_porosity,
+ grout_heat_capacity,
+ grout_thermal_conductivity};
- MaterialLib::Fluid::FluidProperty::ArrayType vars;
- vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] =
- 298.15;
- vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::p)] =
- 101325.0;
+ auto const& refrigerant_config = config.getConfigSubtree("refrigerant");
+ double const refrigerant_density =
+ refrigerant_config.getConfigParameter<double>("density");
+ double const refrigerant_viscosity =
+ refrigerant_config.getConfigParameter<double>("viscosity");
+ double const refrigerant_heat_capacity =
+ refrigerant_config.getConfigParameter<double>("specific_heat_capacity");
+ double const refrigerant_thermal_conductivity =
+ refrigerant_config.getConfigParameter<double>("thermal_conductivity");
+ double const refrigerant_reference_temperature =
+ refrigerant_config.getConfigParameter<double>("reference_temperature");
+ RefrigerantProperties const refrigerant{
+ refrigerant_viscosity, refrigerant_density,
+ refrigerant_thermal_conductivity, refrigerant_heat_capacity,
+ refrigerant_reference_temperature};
- BHE::BHE_1U* m_bhe_1u = new BHE::BHE_1U(
- bhe_ply_name,
- bhe_bound_type,
+ auto const flowAndTemperatureControl = createFlowAndTemperatureControl(
+ config.getConfigSubtree("flow_and_temperature_control"),
curves,
- {bhe_length, bhe_diameter} /* Borehole Geometry */,
- {bhe_pipe_inner_radius, bhe_pipe_outer_radius,
- bhe_pipe_in_wall_thickness, bhe_pipe_out_wall_thickness,
- bhe_pipe_wall_thermal_conductivity, 0.38,
- 0.38 /*lambda_p_i and lambda_p_o will not be used for 1U BHE*/}
- /* Pipe
- Parameters
- */
- ,
- {bhe_refrigerant_viscosity->getValue(vars),
- bhe_refrigerant_density->getValue(vars),
- bhe_regrigerant_heat_conductivity->getValue(vars),
- bhe_refrigerant_heat_capacity->getValue(vars),
- bhe_fluid_longitudinal_dispsion_length} /* Refrigerant Parameters */,
- {bhe_grout_density, bhe_grout_porosity, bhe_grout_heat_capacity,
- bhe_grout_thermal_conductivity} /* Grout Parameters */,
- {bhe_use_ext_therm_resis, bhe_intern_resistance,
- bhe_therm_resistance} /* If using given Ra Rb values*/,
- {bhe_user_defined_therm_resis, bhe_R_fig, bhe_R_fog, bhe_R_gg1,
- bhe_R_gg2, bhe_R_gs} /* If using customed thermal resistance values*/,
- bhe_refrigerant_flow_rate,
- bhe_pipe_distance,
- bhe_power_in_watt_val,
- bhe_delta_T_val,
- bhe_if_use_flow_rate_curve,
- bhe_switch_off_threshold);
+ refrigerant);
- return m_bhe_1u;
+ return {borehole, refrigerant, grout, flowAndTemperatureControl, pipes};
}
} // namespace BHE
} // namespace HeatTransportBHE
-} // namespace ProcessLib
\ No newline at end of file
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.h b/ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.h
index cfbae1c1efd..18dd0a28bcf 100644
--- a/ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.h
+++ b/ProcessLib/HeatTransportBHE/BHE/CreateBHE1U.h
@@ -10,27 +10,22 @@
#pragma once
#include "BHE_1U.h"
-#include "MaterialLib/Fluid/FluidProperty.h"
-#include "ProcessLib/Process.h"
+
+namespace BaseLib
+{
+class ConfigTree;
+}
namespace ProcessLib
{
namespace HeatTransportBHE
{
-namespace BHE // namespace of borehole heat exchanger
+namespace BHE
{
-BHE::BHE_1U* CreateBHE1U(
+BHE::BHE_1U createBHE1U(
BaseLib::ConfigTree const& bhe_conf,
std::map<std::string,
std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
- curves,
- std::unique_ptr<MaterialLib::Fluid::FluidProperty> const&
- bhe_refrigerant_density,
- std::unique_ptr<MaterialLib::Fluid::FluidProperty> const&
- bhe_refrigerant_viscosity,
- std::unique_ptr<MaterialLib::Fluid::FluidProperty> const&
- bhe_refrigerant_heat_capacity,
- std::unique_ptr<MaterialLib::Fluid::FluidProperty> const&
- bhe_regrigerant_heat_conductivity);
+ curves);
} // namespace BHE
} // namespace HeatTransportBHE
} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/CreateFlowAndTemperatureControl.h b/ProcessLib/HeatTransportBHE/BHE/CreateFlowAndTemperatureControl.h
new file mode 100644
index 00000000000..ffaf74f7828
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/CreateFlowAndTemperatureControl.h
@@ -0,0 +1,81 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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 "RefrigerantProperties.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+FlowAndTemperatureControl createFlowAndTemperatureControl(
+ BaseLib::ConfigTree const& config,
+ std::map<std::string,
+ std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
+ curves,
+ RefrigerantProperties const& refrigerant)
+{
+ auto find_curve_or_error = [&](std::string const& name,
+ std::string const& error_message)
+ -> MathLib::PiecewiseLinearInterpolation const& {
+ auto const it = curves.find(name);
+ if (it == curves.end())
+ {
+ ERR(error_message.c_str());
+ OGS_FATAL(
+ "Curve with name '%s' could not be found in the curves list.",
+ name.c_str());
+ }
+ return *it->second;
+ };
+
+ auto const type = config.getConfigParameter<std::string>("type");
+ if (type == "TemperatureCurveConstantFlow")
+ {
+ double const flow_rate = config.getConfigParameter<double>("flow_rate");
+
+ auto const& temperature_curve = find_curve_or_error(
+ config.getConfigParameter<std::string>("temperature_curve"),
+ "Required temperature curve not found.");
+
+ return TemperatureCurveConstantFlow{flow_rate, temperature_curve};
+ }
+ if (type == "FixedPowerConstantFlow")
+ {
+ double const power = config.getConfigParameter<double>("power");
+
+ double const flow_rate = config.getConfigParameter<double>("flow_rate");
+ return FixedPowerConstantFlow{flow_rate, power,
+ refrigerant.specific_heat_capacity,
+ refrigerant.density};
+ }
+
+ if (type == "FixedPowerFlowCurve")
+ {
+ auto const& flow_rate_curve = find_curve_or_error(
+ config.getConfigParameter<std::string>("flow_rate_curve"),
+ "Required flow rate curve not found.");
+
+ double const power = config.getConfigParameter<double>("power");
+
+ return FixedPowerFlowCurve{flow_rate_curve, power,
+ refrigerant.specific_heat_capacity,
+ refrigerant.density};
+ }
+ OGS_FATAL("FlowAndTemperatureControl type '%s' is not implemented.",
+ type.c_str());
+}
+} // namespace BHE
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/FlowAndTemperatureControl.h b/ProcessLib/HeatTransportBHE/BHE/FlowAndTemperatureControl.h
new file mode 100644
index 00000000000..9897b4d13ca
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/FlowAndTemperatureControl.h
@@ -0,0 +1,73 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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 <boost/variant.hpp>
+#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+struct FlowAndTemperature
+{
+ double const flow_rate;
+ double const temperature;
+};
+
+struct TemperatureCurveConstantFlow
+{
+ FlowAndTemperature operator()(double const /*T_out*/,
+ double const time) const
+ {
+ return {flow_rate, temperature_curve.getValue(time)};
+ }
+ double flow_rate;
+ MathLib::PiecewiseLinearInterpolation const& temperature_curve;
+};
+
+struct FixedPowerConstantFlow
+{
+ FlowAndTemperature operator()(double const T_out,
+ double const /*time*/) const
+ {
+ return {flow_rate, power / flow_rate / heat_capacity / density + T_out};
+ }
+ double flow_rate;
+ double power; // Value is expected to be in Watt.
+ double heat_capacity;
+ double density;
+};
+
+struct FixedPowerFlowCurve
+{
+ FlowAndTemperature operator()(double const T_out, double const time) const
+ {
+ double const flow_rate = flow_curve.getValue(time);
+ return {flow_rate, power / flow_rate / heat_capacity / density + T_out};
+ }
+ MathLib::PiecewiseLinearInterpolation const& flow_curve;
+
+ double power; // Value is expected to be in Watt.
+ double heat_capacity;
+ double density;
+};
+
+using FlowAndTemperatureControl = boost::variant<TemperatureCurveConstantFlow,
+ FixedPowerConstantFlow,
+ FixedPowerFlowCurve>;
+
+} // namespace BHE
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/MeshUtils.cpp b/ProcessLib/HeatTransportBHE/BHE/MeshUtils.cpp
index 844ce563fa9..54059b5f75b 100644
--- a/ProcessLib/HeatTransportBHE/BHE/MeshUtils.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/MeshUtils.cpp
@@ -14,58 +14,82 @@
#include "MeshLib/MeshSearch/NodeSearch.h"
#include "MeshLib/Node.h"
+namespace
+{
+std::vector<MeshLib::Element*> extractOneDimensionalElements(
+ std::vector<MeshLib::Element*> const& elements)
+{
+ std::vector<MeshLib::Element*> one_dimensional_elements;
+
+ copy_if(begin(elements), end(elements),
+ back_inserter(one_dimensional_elements),
+ [](MeshLib::Element* e) { return e->getDimension() == 1; });
+
+ return one_dimensional_elements;
+}
+
+std::vector<int> getUniqueMaterialIds(
+ std::vector<int> const& material_ids,
+ std::vector<MeshLib::Element*> const& elements)
+{
+ std::set<int> unique_material_ids;
+ std::transform(begin(elements), end(elements),
+ inserter(unique_material_ids, end(unique_material_ids)),
+ [&material_ids](MeshLib::Element const* const e) {
+ return material_ids[e->getID()];
+ });
+ return {begin(unique_material_ids), end(unique_material_ids)};
+}
+} // namespace
+
namespace ProcessLib
{
namespace HeatTransportBHE
{
BHEMeshData getBHEDataInMesh(MeshLib::Mesh const& mesh)
{
- BHEMeshData bheMeshData;
- // partition all the mesh elements, and copy them into
- // two seperate vectors, one with only matrix elements
- // and the other only BHE elements
- bheMeshData.soil_elements.reserve(mesh.getNumberOfElements());
- std::vector<MeshLib::Element*> all_BHE_elements;
- auto& all_mesh_elements = mesh.getElements();
- std::partition_copy(
- std::begin(all_mesh_elements),
- std::end(all_mesh_elements),
- std::back_inserter(all_BHE_elements),
- std::back_inserter(bheMeshData.soil_elements),
- [](MeshLib::Element* e) { return e->getDimension() == 1; });
+ std::vector<MeshLib::Element*> const all_bhe_elements =
+ extractOneDimensionalElements(mesh.getElements());
+
+ // finally counting two types of elements
+ // They are (i) soil, and (ii) BHE type of elements
+ DBUG("-> found total %d soil elements and %d BHE elements",
+ mesh.getNumberOfElements() - all_bhe_elements.size(),
+ all_bhe_elements.size());
// get BHE material IDs
- auto opt_material_ids = MeshLib::materialIDs(mesh);
+ auto const* const opt_material_ids = MeshLib::materialIDs(mesh);
if (opt_material_ids == nullptr)
{
OGS_FATAL("Not able to get material IDs! ");
}
- for (MeshLib::Element* e : all_BHE_elements)
- bheMeshData.BHE_mat_IDs.push_back((*opt_material_ids)[e->getID()]);
- BaseLib::makeVectorUnique(bheMeshData.BHE_mat_IDs);
- DBUG("-> found %d BHE material groups", bheMeshData.BHE_mat_IDs.size());
+ auto const& material_ids = *opt_material_ids;
+
+ auto const& bhe_material_ids =
+ getUniqueMaterialIds(material_ids, all_bhe_elements);
+ DBUG("-> found %d BHE material groups", bhe_material_ids.size());
// create a vector of BHE elements for each group
- bheMeshData.BHE_elements.resize(bheMeshData.BHE_mat_IDs.size());
- for (unsigned bhe_id = 0; bhe_id < bheMeshData.BHE_mat_IDs.size(); bhe_id++)
+ std::vector<std::vector<MeshLib::Element*>> bhe_elements;
+ bhe_elements.resize(bhe_material_ids.size());
+ for (unsigned bhe_id = 0; bhe_id < bhe_material_ids.size(); bhe_id++)
{
- const auto bhe_mat_id = bheMeshData.BHE_mat_IDs[bhe_id];
- std::vector<MeshLib::Element*>& vec_elements =
- bheMeshData.BHE_elements[bhe_id];
- std::copy_if(all_BHE_elements.begin(), all_BHE_elements.end(),
- std::back_inserter(vec_elements),
- [&](MeshLib::Element* e) {
- return (*opt_material_ids)[e->getID()] == bhe_mat_id;
- });
+ const auto bhe_mat_id = bhe_material_ids[bhe_id];
+ std::vector<MeshLib::Element*>& vec_elements = bhe_elements[bhe_id];
+ copy_if(begin(all_bhe_elements), end(all_bhe_elements),
+ back_inserter(vec_elements), [&](MeshLib::Element* e) {
+ return material_ids[e->getID()] == bhe_mat_id;
+ });
DBUG("-> found %d elements on the BHE_%d", vec_elements.size(), bhe_id);
}
// get a vector of BHE nodes
- bheMeshData.BHE_nodes.resize(bheMeshData.BHE_mat_IDs.size());
- for (unsigned bhe_id = 0; bhe_id < bheMeshData.BHE_mat_IDs.size(); bhe_id++)
+ std::vector<std::vector<MeshLib::Node*>> bhe_nodes;
+ bhe_nodes.resize(bhe_material_ids.size());
+ for (unsigned bhe_id = 0; bhe_id < bhe_material_ids.size(); bhe_id++)
{
- std::vector<MeshLib::Node*>& vec_nodes = bheMeshData.BHE_nodes[bhe_id];
- for (MeshLib::Element* e : bheMeshData.BHE_elements[bhe_id])
+ std::vector<MeshLib::Node*>& vec_nodes = bhe_nodes[bhe_id];
+ for (MeshLib::Element* e : bhe_elements[bhe_id])
{
for (unsigned i = 0; i < e->getNumberOfNodes(); i++)
{
@@ -79,21 +103,7 @@ BHEMeshData getBHEDataInMesh(MeshLib::Mesh const& mesh)
DBUG("-> found %d nodes on the BHE_%d", vec_nodes.size(), bhe_id);
}
- // get all the nodes into the pure soil nodes vector
- bheMeshData.soil_nodes.reserve(mesh.getNumberOfNodes());
- for (MeshLib::Node* n : mesh.getNodes())
- {
- // All elements are counted as a soil element
- bheMeshData.soil_nodes.push_back(n);
- }
-
- // finalLy counting two types of elements
- // They are (i) soil, and (ii) BHE type of elements
- DBUG("-> found total %d soil elements and %d BHE elements",
- bheMeshData.soil_elements.size(),
- bheMeshData.BHE_elements.size());
-
- return bheMeshData;
+ return {bhe_material_ids, bhe_elements, bhe_nodes};
}
} // end of namespace HeatTransportBHE
} // namespace ProcessLib
\ No newline at end of file
diff --git a/ProcessLib/HeatTransportBHE/BHE/MeshUtils.h b/ProcessLib/HeatTransportBHE/BHE/MeshUtils.h
index d275b896ec2..5b5a7978cf3 100644
--- a/ProcessLib/HeatTransportBHE/BHE/MeshUtils.h
+++ b/ProcessLib/HeatTransportBHE/BHE/MeshUtils.h
@@ -21,13 +21,22 @@ namespace ProcessLib
{
namespace HeatTransportBHE
{
+/* TODO (naumov) Just an idea
+struct BheMeshSubset
+{
+ int material_id;
+ std::vector<MeshLib::Element*> elements;
+ std::vector<MeshLib::Node*> nodes;
+};
+*/
+
struct BHEMeshData
{
- std::vector<MeshLib::Element*> soil_elements;
std::vector<int> BHE_mat_IDs;
std::vector<std::vector<MeshLib::Element*>> BHE_elements;
- std::vector<MeshLib::Node*> soil_nodes;
std::vector<std::vector<MeshLib::Node*>> BHE_nodes;
+
+ // TODO (naumov) Just an idea: std::vector<BheMeshSubset> mesh_subsets;
};
/**
diff --git a/ProcessLib/HeatTransportBHE/BHE/Physics.h b/ProcessLib/HeatTransportBHE/BHE/Physics.h
new file mode 100644
index 00000000000..446543da84e
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/Physics.h
@@ -0,0 +1,64 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE // namespace of borehole heat exchanger
+{
+inline double prandtlNumber(double const& viscosity,
+ double const& heat_capacity,
+ double const& heat_conductivity)
+{
+ return viscosity * heat_capacity / heat_conductivity;
+}
+
+inline double reynoldsNumber(double const velocity_norm,
+ double const pipe_diameter,
+ double const viscosity,
+ double const density)
+{
+ return velocity_norm * pipe_diameter / (viscosity / density);
+}
+
+inline double nusseltNumber(double const reynolds_number,
+ double const prandtl_number,
+ double const pipe_diameter,
+ double const pipe_length)
+{
+ if (reynolds_number < 2300.0)
+ {
+ return 4.364;
+ }
+ if (reynolds_number < 10000.0)
+ {
+ double const gamma = (reynolds_number - 2300) / (10000 - 2300);
+
+ return (1.0 - gamma) * 4.364 +
+ gamma *
+ ((0.0308 / 8.0 * 1.0e4 * prandtl_number) /
+ (1.0 + 12.7 * std::sqrt(0.0308 / 8.0) *
+ (std::pow(prandtl_number, 2.0 / 3.0) - 1.0)) *
+ (1.0 + std::pow(pipe_diameter / pipe_length, 2.0 / 3.0)));
+ }
+
+ double const xi = std::pow(1.8 * std::log10(reynolds_number) - 1.5, -2.0);
+ return (xi / 8.0 * reynolds_number * prandtl_number) /
+ (1.0 + 12.7 * std::sqrt(xi / 8.0) *
+ (std::pow(prandtl_number, 2.0 / 3.0) - 1.0)) *
+ (1.0 + std::pow(pipe_diameter / pipe_length, 2.0 / 3.0));
+}
+} // end of namespace BHE
+} // end of namespace HeatTransportBHE
+} // end of namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/Pipe.cpp b/ProcessLib/HeatTransportBHE/BHE/Pipe.cpp
new file mode 100644
index 00000000000..08356225745
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/Pipe.cpp
@@ -0,0 +1,32 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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 "Pipe.h"
+#include "BaseLib/ConfigTree.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+Pipe createPipe(BaseLib::ConfigTree const& config)
+{
+ const double radius = config.getConfigParameter<double>("radius");
+ const double wall_thickness =
+ config.getConfigParameter<double>("wall_thickness");
+ const double wall_thermal_conductivity =
+ config.getConfigParameter<double>("wall_thermal_conductivity");
+ return {radius, wall_thickness, wall_thermal_conductivity};
+}
+} // namespace BHE
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/Pipe.h b/ProcessLib/HeatTransportBHE/BHE/Pipe.h
new file mode 100644
index 00000000000..1eed3b50187
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/Pipe.h
@@ -0,0 +1,43 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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 <boost/math/constants/constants.hpp>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+struct Pipe
+{
+ double const diameter;
+ double const wall_thickness;
+ double const wall_thermal_conductivity;
+
+ double area() const
+ {
+ constexpr double pi = boost::math::constants::pi<double>();
+ return pi * diameter * diameter / 4;
+ }
+};
+
+Pipe createPipe(BaseLib::ConfigTree const& config);
+
+} // namespace BHE
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/PipeConfiguration1U.h b/ProcessLib/HeatTransportBHE/BHE/PipeConfiguration1U.h
new file mode 100644
index 00000000000..72f4d8b4383
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/PipeConfiguration1U.h
@@ -0,0 +1,37 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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 <boost/math/constants/constants.hpp>
+
+#include "BaseLib/ConfigTree.h"
+#include "Pipe.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+struct PipeConfiguration1U
+{
+ Pipe const inlet;
+ Pipe const outlet;
+
+ /// Distance between pipes.
+ double const distance;
+
+ double const longitudinal_despersion_length;
+};
+} // namespace BHE
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/PipeParameters.h b/ProcessLib/HeatTransportBHE/BHE/PipeParameters.h
deleted file mode 100644
index 1c016d74911..00000000000
--- a/ProcessLib/HeatTransportBHE/BHE/PipeParameters.h
+++ /dev/null
@@ -1,67 +0,0 @@
-/**
- * \file
- *
- * \copyright
- * Copyright (c) 2012-2018, 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
-
-namespace ProcessLib
-{
-namespace HeatTransportBHE
-{
-namespace BHE
-{
-struct PipeParameters
-{
- /**
- * radius of the pipline inner side
- * unit is m
- */
- double const r_inner;
-
- /**
- * radius of the pipline outer side
- * unit is m
- */
- double const r_outer;
-
- /**
- * pipe-in wall thickness
- * unit is m
- */
- double const b_in;
-
- /**
- * pipe-out wall thickness
- * unit is m
- */
- double const b_out;
-
- /**
- * thermal conductivity of the pipe wall
- * unit is kg m sec^-3 K^-1
- */
-
- double const lambda_p;
-
- /**
- * thermal conductivity of the inner pipe wall
- * unit is kg m sec^-3 K^-1
- */
- double const lambda_p_i;
-
- /**
- * thermal conductivity of the outer pipe wall
- * unit is kg m sec^-3 K^-1
- */
- double const lambda_p_o;
-};
-} // namespace BHE
-} // namespace HeatTransportBHE
-} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/BHE/RefrigerantParameters.h b/ProcessLib/HeatTransportBHE/BHE/RefrigerantProperties.h
similarity index 59%
rename from ProcessLib/HeatTransportBHE/BHE/RefrigerantParameters.h
rename to ProcessLib/HeatTransportBHE/BHE/RefrigerantProperties.h
index 8b46b31a351..ea54310ab99 100644
--- a/ProcessLib/HeatTransportBHE/BHE/RefrigerantParameters.h
+++ b/ProcessLib/HeatTransportBHE/BHE/RefrigerantProperties.h
@@ -17,37 +17,29 @@ namespace HeatTransportBHE
{
namespace BHE
{
-struct RefrigerantParameters
+struct RefrigerantProperties
{
/**
- * dynamics viscosity of the refrigerant
* unit is kg m-1 sec-1
*/
- double const mu_r;
+ double const dynamic_viscosity;
/**
- * density of the refrigerant
* unit is kg m-3
*/
- double const rho_r;
+ double const density;
/**
- * thermal conductivity of the refrigerant
* unit is kg m sec^-3 K^-1
*/
- double const lambda_r;
+ double const thermal_conductivity;
/**
- * specific heat capacity of the refrigerant
* unit is m^2 sec^-2 K^-1
*/
- double const heat_cap_r;
+ double const specific_heat_capacity;
- /**
- * longitudinal dispersivity of the
- * referigerant flow in the pipeline
- */
- double const alpha_L;
+ double const reference_temperature;
};
} // namespace BHE
} // namespace HeatTransportBHE
diff --git a/ProcessLib/HeatTransportBHE/BHE/ThermoMechanicalFlowProperties.h b/ProcessLib/HeatTransportBHE/BHE/ThermoMechanicalFlowProperties.h
new file mode 100644
index 00000000000..1d2b2e0128e
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/BHE/ThermoMechanicalFlowProperties.h
@@ -0,0 +1,71 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2018, 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 "Physics.h"
+#include "PipeParameters.h"
+#include "RefrigerantProperties.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+namespace BHE
+{
+struct ThermoMechanicalFlowProperties
+{
+ double velocity;
+ double nusselt_number;
+};
+
+inline ThermoMechanicalFlowProperties
+calculateThermoMechanicalFlowPropertiesPipe(PipeParameters const& pipe,
+ double const length,
+ RefrigerantProperties const& fluid,
+ double const flow_rate)
+{
+ double const Pr =
+ prandtlNumber(fluid.mu_r, fluid.specific_heat_capacity, fluid.lambda_r);
+
+ double const velocity = pipeFlowVelocity(flow_rate, pipe.r_inner);
+ double const Re =
+ reynoldsNumber(velocity, 2.0 * pipe.r_inner, fluid.mu_r, fluid.density);
+ double const nusselt_number =
+ nusseltNumber(Re, Pr, 2 * pipe.r_inner, length);
+ return {velocity, nusselt_number};
+}
+
+inline ThermoMechanicalFlowProperties
+calculateThermoMechanicalFlowPropertiesAnnulus(
+ PipeParameters const& pipe, double const length,
+ RefrigerantProperties const& fluid, double const flow_rate)
+{
+ double const Pr =
+ prandtlNumber(fluid.mu_r, fluid.specific_heat_capacity, fluid.lambda_r);
+
+ double const velocity =
+ annulusFlowVelocity(flow_rate, pipe.r_outer, pipe.r_inner + pipe.b_in);
+
+ double const Re = reynoldsNumber(
+ velocity, 2.0 * (pipe.r_outer - (pipe.r_inner + pipe.b_in)), fluid.mu_r,
+ fluid.density);
+
+ double const diameter_ratio = (pipe.r_inner + pipe.b_in) / pipe.r_outer;
+ double const pipe_aspect_ratio =
+ (2 * pipe.r_outer - 2 * (pipe.r_inner + pipe.b_in)) / length;
+ double const nusselt_number =
+ nusseltNumberAnnulus(Re, Pr, diameter_ratio, pipe_aspect_ratio);
+ return {velocity, nusselt_number};
+}
+} // namespace BHE
+} // namespace HeatTransportBHE
+} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.cpp b/ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.cpp
index ece955939dc..159eef896ae 100644
--- a/ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.cpp
+++ b/ProcessLib/HeatTransportBHE/CreateHeatTransportBHEProcess.cpp
@@ -9,25 +9,15 @@
#include "CreateHeatTransportBHEProcess.h"
-#include "BHE/BHEAbstract.h"
-#include "BHE/BHE_1U.h"
-#include "BHE/BHE_2U.h"
-#include "BHE/BHE_CXA.h"
-#include "BHE/BHE_CXC.h"
+#include <vector>
+
+#include "ProcessLib/Output/CreateSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+#include "BHE/BHETypes.h"
#include "BHE/CreateBHE1U.h"
-#include "BHE/CreateBHE2U.h"
-#include "BHE/CreateBHECXA.h"
-#include "BHE/CreateBHECXC.h"
-#include "BaseLib/Algorithm.h"
#include "HeatTransportBHEProcess.h"
#include "HeatTransportBHEProcessData.h"
-#include "MaterialLib/Fluid/Density/CreateFluidDensityModel.h"
-#include "MaterialLib/Fluid/FluidProperty.h"
-#include "MaterialLib/Fluid/SpecificHeatCapacity/CreateSpecificFluidHeatCapacityModel.h"
-#include "MaterialLib/Fluid/ThermalConductivity/CreateFluidThermalConductivityModel.h"
-#include "MaterialLib/Fluid/Viscosity/CreateViscosityModel.h"
-#include "ProcessLib/Output/CreateSecondaryVariables.h"
-#include "ProcessLib/Utils/ProcessUtils.h"
namespace ProcessLib
{
@@ -177,106 +167,30 @@ std::unique_ptr<Process> createHeatTransportBHEProcess(
DBUG("Use \'%s\' as gas phase density parameter.",
density_gas.name.c_str());
- // get the refrigerant properties from fluid property class
- //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__material_property__fluid}
- auto const& fluid_config = config.getConfigSubtree("fluid");
- //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__material_property__refrigerant_density}
- auto const& rho_conf = fluid_config.getConfigSubtree("refrigerant_density");
- auto bhe_refrigerant_density =
- MaterialLib::Fluid::createFluidDensityModel(rho_conf);
- //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__material_property__refrigerant_viscosity}
- auto const& mu_conf =
- fluid_config.getConfigSubtree("refrigerant_viscosity");
- auto bhe_refrigerant_viscosity =
- MaterialLib::Fluid::createViscosityModel(mu_conf);
- //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__material_property__refrigerant_specific_heat_capacity}
- auto const& cp_conf =
- fluid_config.getConfigSubtree("refrigerant_specific_heat_capacity");
- auto bhe_refrigerant_heat_capacity =
- MaterialLib::Fluid::createSpecificFluidHeatCapacityModel(cp_conf);
- //! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__material_property__refrigerant_thermal_conductivity}
- auto const& lambda_conf =
- fluid_config.getConfigSubtree("refrigerant_thermal_conductivity");
- auto bhe_regrigerant_heat_conductivity =
- MaterialLib::Fluid::createFluidThermalConductivityModel(lambda_conf);
-
- // reading BHE
- // parameters--------------------------------------------------------------
- std::vector<std::unique_ptr<BHE::BHEAbstract>> vec_BHEs;
- // BHE::BHE_Net BHE_network;
-
- // now read the BHE configurations
+ // reading BHE parameters --------------------------------------------------
+ std::vector<BHE::BHETypes> bhes;
auto const& bhe_configs =
//! \ogs_file_param{prj__processes__process__HEAT_TRANSPORT_BHE__borehole_heat_exchangers}
config.getConfigSubtree("borehole_heat_exchangers");
for (
- auto const& bhe_conf :
- //! \ogs_file_param
- //! prj__processes__process___HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger}
+ auto const& bhe_config :
+ //! \ogs_file_param{prj__processes__process___HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger}
bhe_configs.getConfigSubtreeList("borehole_heat_exchanger"))
{
// read in the parameters
- using namespace BHE;
- const std::string bhe_type_str =
- bhe_conf.getConfigParameter<std::string>("bhe_type");
+ const std::string bhe_type =
+ //! \ogs_file_param{prj__processes__process___HEAT_TRANSPORT_BHE__borehole_heat_exchangers__borehole_heat_exchanger__bhe_type}
+ bhe_config.getConfigParameter<std::string>("type");
- // convert BHE type
- if (bhe_type_str == "BHE_TYPE_1U")
- {
- // initialize the 1U type BHE
- BHE::BHE_1U* m_bhe_1u =
- BHE::CreateBHE1U(bhe_conf,
- curves,
- bhe_refrigerant_density,
- bhe_refrigerant_viscosity,
- bhe_refrigerant_heat_capacity,
- bhe_regrigerant_heat_conductivity);
-
- vec_BHEs.emplace_back(std::make_unique<BHE_1U>(*m_bhe_1u));
- }
- else if (bhe_type_str == "BHE_TYPE_2U")
- {
- // initialize the 2U type BHE
- BHE::BHE_2U* m_bhe_2u =
- BHE::CreateBHE2U(bhe_conf,
- curves,
- bhe_refrigerant_density,
- bhe_refrigerant_viscosity,
- bhe_refrigerant_heat_capacity,
- bhe_regrigerant_heat_conductivity);
-
- vec_BHEs.emplace_back(std::make_unique<BHE_2U>(*m_bhe_2u));
- }
- else if (bhe_type_str == "BHE_TYPE_CXC")
- {
- // initialize the CXC type BHE
- BHE::BHE_CXC* m_bhe_CXC =
- BHE::CreateBHECXC(bhe_conf,
- curves,
- bhe_refrigerant_density,
- bhe_refrigerant_viscosity,
- bhe_refrigerant_heat_capacity,
- bhe_regrigerant_heat_conductivity);
-
- vec_BHEs.emplace_back(std::make_unique<BHE_CXC>(*m_bhe_CXC));
- }
- else if (bhe_type_str == "BHE_TYPE_CXA")
+ if (bhe_type == "1U")
{
- // initialize the CXA type BHE
- BHE::BHE_CXA* m_bhe_CXA =
- BHE::CreateBHECXA(bhe_conf,
- curves,
- bhe_refrigerant_density,
- bhe_refrigerant_viscosity,
- bhe_refrigerant_heat_capacity,
- bhe_regrigerant_heat_conductivity);
-
- vec_BHEs.emplace_back(std::make_unique<BHE_CXA>(*m_bhe_CXA));
+ bhes.push_back(BHE::createBHE1U(bhe_config, curves));
+ continue;
}
+ OGS_FATAL("Unknown BHE type '%s'.", bhe_type.c_str());
}
- // end of reading BHE
- // parameters-------------------------------------------------------
+ // end of reading BHE parameters -------------------------------------------
HeatTransportBHEProcessData process_data{thermal_conductivity_solid,
thermal_conductivity_fluid,
@@ -287,7 +201,7 @@ std::unique_ptr<Process> createHeatTransportBHEProcess(
density_solid,
density_fluid,
density_gas,
- std::move(vec_BHEs)};
+ std::move(bhes)};
SecondaryVariableCollection secondary_variables;
diff --git a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
index 5a9da59c74d..2d609b204ad 100644
--- a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
+++ b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.cpp
@@ -10,7 +10,6 @@
#include "HeatTransportBHEProcess.h"
#include <cassert>
-#include <iostream>
#include "ProcessLib/HeatTransportBHE/BHE/MeshUtils.h"
#include "ProcessLib/HeatTransportBHE/LocalAssemblers/CreateLocalAssemblers.h"
@@ -44,9 +43,8 @@ HeatTransportBHEProcess::HeatTransportBHEProcess(
_process_data._vec_BHE_property.size())
{
OGS_FATAL(
- "The number of the given BHE properties (%d) are not "
- "consistent"
- " with the number of BHE groups in a mesh (%d).",
+ "The number of the given BHE properties (%d) are not consistent "
+ "with the number of BHE groups in a mesh (%d).",
_process_data._vec_BHE_property.size(),
_bheMeshData.BHE_mat_IDs.size());
}
@@ -57,10 +55,7 @@ HeatTransportBHEProcess::HeatTransportBHEProcess(
OGS_FATAL("Not able to get material IDs! ");
}
// create a map from a material ID to a BHE ID
- auto max_BHE_mat_id =
- std::max_element(material_ids->begin(), material_ids->end());
- _process_data._map_materialID_to_BHE_ID.resize(*max_BHE_mat_id + 1);
- for (std::size_t i = 0; i < _bheMeshData.BHE_mat_IDs.size(); i++)
+ for (int i = 0; i < static_cast<int>(_bheMeshData.BHE_mat_IDs.size()); i++)
{
// fill in the map structure
_process_data._map_materialID_to_BHE_ID[_bheMeshData.BHE_mat_IDs[i]] =
@@ -76,62 +71,49 @@ void HeatTransportBHEProcess::constructDofTable()
_mesh_subset_all_nodes =
std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
- //------------------------------------------------------------
- // prepare mesh subsets to define DoFs
- //------------------------------------------------------------
- // first all the soil nodes
- _mesh_subset_pure_soil_nodes =
- std::make_unique<MeshLib::MeshSubset>(_mesh, _bheMeshData.soil_nodes);
+ //
+ // Soil temperature variable defined on the whole mesh.
+ //
+ _mesh_subset_soil_nodes =
+ std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
+ std::vector<MeshLib::MeshSubset> all_mesh_subsets{*_mesh_subset_soil_nodes};
- std::vector<int> vec_n_BHE_unknowns;
- vec_n_BHE_unknowns.reserve(_bheMeshData.BHE_nodes.size());
- // the BHE nodes need to be cherry-picked from the vector
- for (unsigned i = 0; i < _bheMeshData.BHE_nodes.size(); i++)
- {
- _mesh_subset_BHE_nodes.push_back(
- std::make_unique<MeshLib::MeshSubset const>(
- _mesh, _bheMeshData.BHE_nodes[i]));
- int n_BHE_unknowns =
- _process_data._vec_BHE_property[i]->getNumUnknowns();
- vec_n_BHE_unknowns.emplace_back(n_BHE_unknowns);
- }
+ std::vector<std::vector<MeshLib::Element*> const*> vec_var_elements;
+ vec_var_elements.push_back(&(_mesh.getElements()));
- // Collect the mesh subsets in a vector.
- // All the soil nodes has 1 temperature variable
- std::vector<MeshLib::MeshSubset> all_mesh_subsets{
- *_mesh_subset_pure_soil_nodes};
+ std::vector<int> vec_n_components{
+ 1}; // one component for the soil temperature variable.
- // All the BHE nodes have additinal variables
- std::size_t count = 0;
- for (auto& ms : _mesh_subset_BHE_nodes)
- {
- std::generate_n(std::back_inserter(all_mesh_subsets),
- // Here the number of components equals to
- // the number of unknowns on the BHE
- vec_n_BHE_unknowns[count],
- [&]() { return *ms; });
- count++;
- }
+ //
+ // BHE nodes with BHE type dependend number of variables.
+ //
+ int const n_BHEs = _process_data._vec_BHE_property.size();
+ assert(n_BHEs == static_cast<int>(_bheMeshData.BHE_mat_IDs.size()));
+ assert(n_BHEs == static_cast<int>(_bheMeshData.BHE_nodes.size()));
+ assert(n_BHEs == static_cast<int>(_bheMeshData.BHE_elements.size()));
- std::vector<int> vec_n_components;
- // This is the soil temperature for first mesh subset.
- // 1, because for the soil part there is just one variable which is the soil
- // temperature.
- vec_n_components.push_back(1);
- // now the BHE subsets
- for (std::size_t i = 0; i < _bheMeshData.BHE_mat_IDs.size(); i++)
+ // the BHE nodes need to be cherry-picked from the vector
+ for (int i = 0; i < n_BHEs; i++)
{
- // Here the number of components equals to
- // the number of unknowns on the BHE
- vec_n_components.push_back(vec_n_BHE_unknowns[i]);
- }
+ auto const number_of_unknowns = apply_visitor(
+ [](auto const& bhe) { return bhe.number_of_unknowns; },
+ _process_data._vec_BHE_property[i]);
+ auto const& bhe_nodes = _bheMeshData.BHE_nodes[i];
+ auto const& bhe_elements = _bheMeshData.BHE_elements[i];
- std::vector<std::vector<MeshLib::Element*> const*> vec_var_elements;
- // vec_var_elements.push_back(&_vec_pure_soil_elements);
- vec_var_elements.push_back(&(_mesh.getElements()));
- for (std::size_t i = 0; i < _bheMeshData.BHE_elements.size(); i++)
- {
- vec_var_elements.push_back(&_bheMeshData.BHE_elements[i]);
+ // All the BHE nodes have additional variables.
+ _mesh_subset_BHE_nodes.push_back(
+ std::make_unique<MeshLib::MeshSubset const>(_mesh, bhe_nodes));
+
+ std::generate_n(
+ std::back_inserter(all_mesh_subsets),
+ // Here the number of components equals to the
+ // number of unknowns on the BHE
+ number_of_unknowns,
+ [& ms = _mesh_subset_BHE_nodes.back()]() { return *ms; });
+
+ vec_n_components.push_back(number_of_unknowns);
+ vec_var_elements.push_back(&bhe_elements);
}
_local_to_global_index_map =
@@ -150,27 +132,27 @@ void HeatTransportBHEProcess::initializeConcreteProcess(
MeshLib::Mesh const& mesh,
unsigned const integration_order)
{
+ // Quick access map to BHE's through element ids.
+ std::unordered_map<std::size_t, BHE::BHETypes*> element_to_bhe_map;
+ int const n_BHEs = _process_data._vec_BHE_property.size();
+ for (int i = 0; i < n_BHEs; i++)
+ {
+ auto const& bhe_elements = _bheMeshData.BHE_elements[i];
+ for (auto const& e : bhe_elements)
+ {
+ element_to_bhe_map[e->getID()] =
+ &_process_data._vec_BHE_property[i];
+ }
+ }
+
assert(mesh.getDimension() == 3);
ProcessLib::HeatTransportBHE::createLocalAssemblers<
- 3 /* mesh dimension */, HeatTransportBHELocalAssemblerSoil,
- HeatTransportBHELocalAssemblerBHE>(
- mesh.getElements(), dof_table, _local_assemblers,
- _process_data._vec_ele_connected_BHE_IDs,
- _process_data._vec_BHE_property, mesh.isAxiallySymmetric(),
- integration_order, _process_data);
-
- // create BHE boundary conditions
- // for each BHE, one BC on the top
- // and one BC at the bottom.
- std::vector<std::unique_ptr<BoundaryCondition>> bc_collections =
- createBHEBoundaryConditionTopBottom();
+ HeatTransportBHELocalAssemblerSoil, HeatTransportBHELocalAssemblerBHE>(
+ mesh.getElements(), dof_table, _local_assemblers, element_to_bhe_map,
+ mesh.isAxiallySymmetric(), integration_order, _process_data);
- const int process_id = 0;
- auto& current_process_BCs = _boundary_conditions[process_id];
- for (auto& bc_collection : bc_collections)
- {
- current_process_BCs.addCreatedBC(std::move(bc_collection));
- }
+ // Create BHE boundary conditions for each of the BHEs
+ createBHEBoundaryConditionTopBottom(_bheMeshData.BHE_nodes);
}
void HeatTransportBHEProcess::assembleConcreteProcess(const double t,
@@ -208,20 +190,18 @@ void HeatTransportBHEProcess::computeSecondaryVariableConcrete(
_coupled_solutions);
}
-std::vector<std::unique_ptr<BoundaryCondition>>
-HeatTransportBHEProcess::createBHEBoundaryConditionTopBottom()
+void HeatTransportBHEProcess::createBHEBoundaryConditionTopBottom(
+ std::vector<std::vector<MeshLib::Node*>> const& all_bhe_nodes)
{
- /**
- * boundary conditions
- */
- std::vector<std::unique_ptr<BoundaryCondition>> bcs;
+ const int process_id = 0;
+ auto& bcs = _boundary_conditions[process_id];
int const n_BHEs = static_cast<int>(_process_data._vec_BHE_property.size());
// for each BHE
for (int bhe_i = 0; bhe_i < n_BHEs; bhe_i++)
{
- auto const& bhe_nodes = _bheMeshData.BHE_nodes[bhe_i];
+ auto const& bhe_nodes = all_bhe_nodes[bhe_i];
// find the variable ID
// the soil temperature is 0-th variable
// the BHE temperature is therefore bhe_i + 1
@@ -234,7 +214,7 @@ HeatTransportBHEProcess::createBHEBoundaryConditionTopBottom()
return a->getCoords()[2] < b->getCoords()[2];
});
auto const bc_bottom_node_id = (*bottom_top_nodes.first)->getID();
- MeshLib::Node* const bc_top_node = *bottom_top_nodes.second;
+ auto const bc_top_node_id = (*bottom_top_nodes.second)->getID();
auto get_global_bhe_bc_indices =
[&](std::size_t const node_id,
@@ -248,25 +228,26 @@ HeatTransportBHEProcess::createBHEBoundaryConditionTopBottom()
variable_id, in_out_component_id.second));
};
- for (auto const& in_out_component_id :
- _process_data._vec_BHE_property[bhe_i]
- ->inflowOutflowBcComponentIds())
- {
- // Top, inflow.
- bcs.push_back(ProcessLib::createBHEInflowDirichletBoundaryCondition(
- get_global_bhe_bc_indices(bc_top_node->getID(),
- in_out_component_id),
- _mesh, {bc_top_node}, variable_id, in_out_component_id.first,
- _process_data._vec_BHE_property[bhe_i]));
-
- // Bottom, outflow.
- bcs.push_back(ProcessLib::createBHEBottomDirichletBoundaryCondition(
- get_global_bhe_bc_indices(bc_bottom_node_id,
- in_out_component_id)));
- }
+ auto createBCs = [&](auto& bhe) {
+ for (auto const& in_out_component_id :
+ bhe.inflow_outflow_bc_component_ids)
+ {
+ // Top, inflow.
+ bcs.addBoundaryCondition(
+ ProcessLib::createBHEInflowDirichletBoundaryCondition(
+ get_global_bhe_bc_indices(bc_top_node_id,
+ in_out_component_id),
+ bhe));
+
+ // Bottom, outflow.
+ bcs.addBoundaryCondition(
+ ProcessLib::createBHEBottomDirichletBoundaryCondition(
+ get_global_bhe_bc_indices(bc_bottom_node_id,
+ in_out_component_id)));
+ }
+ };
+ apply_visitor(createBCs, _process_data._vec_BHE_property[bhe_i]);
}
-
- return bcs;
}
} // namespace HeatTransportBHE
} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.h b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.h
index 6af02ecd02a..1cdcba286a8 100644
--- a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.h
+++ b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcess.h
@@ -59,8 +59,8 @@ private:
const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
- std::vector<std::unique_ptr<BoundaryCondition>>
- createBHEBoundaryConditionTopBottom();
+ void createBHEBoundaryConditionTopBottom(
+ std::vector<std::vector<MeshLib::Node*>> const& all_bhe_nodes);
HeatTransportBHEProcessData _process_data;
@@ -73,10 +73,7 @@ private:
std::vector<std::unique_ptr<MeshLib::MeshSubset const>>
_mesh_subset_BHE_soil_nodes;
- std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_pure_soil_nodes;
-
- std::unique_ptr<MeshLib::MeshSubset const>
- _mesh_subset_soil_nodes_connected_with_BHE;
+ std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_soil_nodes;
const BHEMeshData _bheMeshData;
};
diff --git a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcessData.h b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcessData.h
index 28224822f0e..847dce20919 100644
--- a/ProcessLib/HeatTransportBHE/HeatTransportBHEProcessData.h
+++ b/ProcessLib/HeatTransportBHE/HeatTransportBHEProcessData.h
@@ -9,8 +9,10 @@
#pragma once
+#include <unordered_map>
+
#include "MeshLib/PropertyVector.h"
-#include "ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h"
+#include "ProcessLib/HeatTransportBHE/BHE/BHETypes.h"
namespace MeshLib
{
@@ -24,8 +26,6 @@ struct Parameter;
namespace HeatTransportBHE
{
-using namespace BHE;
-
struct HeatTransportBHEProcessData
{
HeatTransportBHEProcessData(
@@ -38,7 +38,7 @@ struct HeatTransportBHEProcessData
Parameter<double> const& density_solid_,
Parameter<double> const& density_fluid_,
Parameter<double> const& density_gas_,
- std::vector<std::unique_ptr<BHE::BHEAbstract>>&& vec_BHEs_)
+ std::vector<BHE::BHETypes>&& vec_BHEs_)
: thermal_conductivity_solid(thermal_conductivity_solid_),
thermal_conductivity_fluid(thermal_conductivity_fluid_),
thermal_conductivity_gas(thermal_conductivity_gas_),
@@ -79,12 +79,9 @@ struct HeatTransportBHEProcessData
Parameter<double> const& density_gas;
MeshLib::PropertyVector<int> const* _mesh_prop_materialIDs = nullptr;
- std::vector<std::size_t> _map_materialID_to_BHE_ID;
-
- std::vector<std::unique_ptr<BHE::BHEAbstract>> _vec_BHE_property;
+ std::unordered_map<int, int> _map_materialID_to_BHE_ID;
- // a table of connected BHE IDs for each element
- std::vector<std::vector<std::size_t>> _vec_ele_connected_BHE_IDs;
+ std::vector<BHE::BHETypes> _vec_BHE_property;
};
} // namespace HeatTransportBHE
} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/CreateLocalAssemblers.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/CreateLocalAssemblers.h
index 03f6444c70d..93ac0d5ce23 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/CreateLocalAssemblers.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/CreateLocalAssemblers.h
@@ -14,7 +14,7 @@
#include "LocalDataInitializer.h"
#include "NumLib/DOF/LocalToGlobalIndexMap.h"
-#include "ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h"
+#include "ProcessLib/HeatTransportBHE/BHE/BHETypes.h"
namespace ProcessLib
{
@@ -22,23 +22,22 @@ namespace HeatTransportBHE
{
namespace detail
{
-template <
- int GlobalDim,
- template <typename, typename, int> class LocalAssemblerSoilImplementation,
- template <typename, typename, int> class LocalAssemblerBHEImplementation,
- typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+template <template <typename, typename> class LocalAssemblerSoilImplementation,
+ template <typename, typename, typename>
+ class LocalAssemblerBHEImplementation,
+ typename LocalAssemblerInterface, typename... ExtraCtorArgs>
void createLocalAssemblers(
NumLib::LocalToGlobalIndexMap const& dof_table,
std::vector<MeshLib::Element*> const& mesh_elements,
std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
- const std::vector<std::vector<std::size_t>>& vec_ele_connected_BHE_IDs,
- const std::vector<std::unique_ptr<BHE::BHEAbstract>>& vec_BHE_property,
+ std::unordered_map<std::size_t, BHE::BHETypes*> const& element_to_bhe_map,
ExtraCtorArgs&&... extra_ctor_args)
{
// Shape matrices initializer
- using LocalDataInitializer = LocalDataInitializer<
- LocalAssemblerInterface, LocalAssemblerSoilImplementation,
- LocalAssemblerBHEImplementation, 3, ExtraCtorArgs...>;
+ using LocalDataInitializer =
+ LocalDataInitializer<LocalAssemblerInterface,
+ LocalAssemblerSoilImplementation,
+ LocalAssemblerBHEImplementation, ExtraCtorArgs...>;
DBUG("Create local assemblers for the HeatTransportBHE process.");
// Populate the vector of local assemblers.
@@ -48,8 +47,8 @@ void createLocalAssemblers(
DBUG("Calling local assembler builder for all mesh elements.");
GlobalExecutor::transformDereferenced(
- initializer, mesh_elements, local_assemblers, vec_ele_connected_BHE_IDs,
- vec_BHE_property, std::forward<ExtraCtorArgs>(extra_ctor_args)...);
+ initializer, mesh_elements, local_assemblers, element_to_bhe_map,
+ std::forward<ExtraCtorArgs>(extra_ctor_args)...);
}
} // namespace detail
@@ -64,11 +63,10 @@ void createLocalAssemblers(
* The first two template parameters cannot be deduced from the arguments.
* Therefore they always have to be provided manually.
*/
-template <
- int GlobalDim,
- template <typename, typename, int> class LocalAssemblerSoilImplementation,
- template <typename, typename, int> class LocalAssemblerBHEImplementation,
- typename LocalAssemblerInterface, typename... ExtraCtorArgs>
+template <template <typename, typename> class LocalAssemblerSoilImplementation,
+ template <typename, typename, typename>
+ class LocalAssemblerBHEImplementation,
+ typename LocalAssemblerInterface, typename... ExtraCtorArgs>
void createLocalAssemblers(
std::vector<MeshLib::Element*> const& mesh_elements,
NumLib::LocalToGlobalIndexMap const& dof_table,
@@ -77,7 +75,7 @@ void createLocalAssemblers(
{
DBUG("Create local assemblers for the HeatTransportBHE process.");
- detail::createLocalAssemblers<GlobalDim, LocalAssemblerSoilImplementation,
+ detail::createLocalAssemblers<LocalAssemblerSoilImplementation,
LocalAssemblerBHEImplementation>(
dof_table, mesh_elements, local_assemblers,
std::forward<ExtraCtorArgs>(extra_ctor_args)...);
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
new file mode 100644
index 00000000000..3fb23405916
--- /dev/null
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
@@ -0,0 +1,199 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, 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 "HeatTransportBHELocalAssemblerBHE.h"
+
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+namespace ProcessLib
+{
+namespace HeatTransportBHE
+{
+template <typename ShapeFunction, typename IntegrationMethod, typename BHEType>
+HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod, BHEType>::
+ HeatTransportBHELocalAssemblerBHE(MeshLib::Element const& e,
+ BHEType const& bhe,
+ bool const is_axially_symmetric,
+ unsigned const integration_order,
+ HeatTransportBHEProcessData& process_data)
+ : _process_data(process_data),
+ _integration_method(integration_order),
+ _bhe(bhe),
+ element_id(e.getID())
+{
+ // need to make sure that the BHE elements are one-dimensional
+ assert(e.getDimension() == 1);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ _ip_data.reserve(n_integration_points);
+ _secondary_data.N.resize(n_integration_points);
+
+ auto const shape_matrices =
+ initShapeMatrices<ShapeFunction, ShapeMatricesType, IntegrationMethod,
+ 3 /* GlobalDim */>(e, is_axially_symmetric,
+ _integration_method);
+
+ // ip data initialization
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto const& sm = shape_matrices[ip];
+ // create the class IntegrationPointDataBHE in place
+ _ip_data.push_back(
+ {sm.N, sm.dNdx,
+ _integration_method.getWeightedPoint(ip).getWeight() *
+ sm.integralMeasure * sm.detJ});
+
+ _secondary_data.N[ip] = sm.N;
+ }
+
+ _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);
+ // formulate the local BHE R matrix
+ for (int idx_bhe_unknowns = 0; idx_bhe_unknowns < bhe_unknowns;
+ idx_bhe_unknowns++)
+ {
+ typename ShapeMatricesType::template MatrixType<
+ single_bhe_unknowns_size, single_bhe_unknowns_size>
+ matBHE_loc_R = ShapeMatricesType::template MatrixType<
+ single_bhe_unknowns_size,
+ single_bhe_unknowns_size>::Zero(single_bhe_unknowns_size,
+ single_bhe_unknowns_size);
+ // Loop over Gauss points
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto const& N = _ip_data[ip].N;
+ auto const& w = _ip_data[ip].integration_weight;
+
+ auto const& R = _bhe.thermalResistance(idx_bhe_unknowns);
+ // calculate mass matrix for current unknown
+ matBHE_loc_R += N.transpose() * N * (1 / R) * w;
+ } // end of loop over integration point
+
+ // 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);
+ } // end of loop over BHE unknowns
+
+ // debugging
+ // std::string sep = "\n----------------------------------------\n";
+ // Eigen::IOFormat CleanFmt(4, 0, ", ", "\n", "[", "]");
+ // std::cout << "_R_matrix: \n" << sep;
+ // std::cout << _R_matrix.format(CleanFmt) << sep;
+ // std::cout << "_R_s_matrix: \n" << sep;
+ // std::cout << _R_s_matrix.format(CleanFmt) << sep;
+ // std::cout << "_R_pi_s_matrix: \n" << sep;
+ // std::cout << _R_pi_s_matrix.format(CleanFmt) << sep;
+}
+
+template <typename ShapeFunction, typename IntegrationMethod, typename BHEType>
+void HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod,
+ BHEType>::
+ assemble(
+ double const /*t*/, std::vector<double> const& /*local_x*/,
+ std::vector<double>& local_M_data, std::vector<double>& local_K_data,
+ std::vector<double>& /*local_b_data*/) // local b vector is not touched
+{
+ auto local_M = MathLib::createZeroedMatrix<BheLocalMatrixType>(
+ local_M_data, local_matrix_size, local_matrix_size);
+ auto local_K = MathLib::createZeroedMatrix<BheLocalMatrixType>(
+ local_K_data, local_matrix_size, local_matrix_size);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ auto const& pipe_heat_capacities = _bhe.pipeHeatCapacities();
+ auto const& pipe_heat_conductions = _bhe.pipeHeatConductions();
+ auto const& pipe_advection_vectors = _bhe.pipeAdvectionVectors();
+ auto const& cross_section_areas = _bhe.cross_section_areas;
+
+ // the mass and conductance matrix terms
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto& ip_data = _ip_data[ip];
+
+ auto const& w = ip_data.integration_weight;
+ auto const& N = ip_data.N;
+ auto const& dNdx = ip_data.dNdx;
+
+ // looping over all unknowns.
+ for (int idx_bhe_unknowns = 0; idx_bhe_unknowns < bhe_unknowns;
+ idx_bhe_unknowns++)
+ {
+ // get coefficient of mass from corresponding BHE.
+ auto const& mass_coeff = pipe_heat_capacities[idx_bhe_unknowns];
+ auto const& lambda = pipe_heat_conductions[idx_bhe_unknowns];
+ auto const& advection_vector =
+ pipe_advection_vectors[idx_bhe_unknowns];
+ auto const& A = cross_section_areas[idx_bhe_unknowns];
+
+ int const single_bhe_unknowns_index =
+ bhe_unknowns_index +
+ single_bhe_unknowns_size * idx_bhe_unknowns;
+ // local M
+ local_M
+ .template block<single_bhe_unknowns_size,
+ single_bhe_unknowns_size>(
+ single_bhe_unknowns_index, single_bhe_unknowns_index)
+ .noalias() += N.transpose() * N * mass_coeff * A * w;
+
+ // local K
+ // laplace part
+ local_K
+ .template block<single_bhe_unknowns_size,
+ single_bhe_unknowns_size>(
+ single_bhe_unknowns_index, single_bhe_unknowns_index)
+ .noalias() += dNdx.transpose() * dNdx * lambda * A * w;
+ // advection part
+ local_K
+ .template block<single_bhe_unknowns_size,
+ single_bhe_unknowns_size>(
+ single_bhe_unknowns_index, single_bhe_unknowns_index)
+ .noalias() +=
+ N.transpose() * advection_vector.transpose() * dNdx * A * w;
+ }
+ }
+
+ // add the R matrix to local_K
+ local_K.template block<bhe_unknowns_size, bhe_unknowns_size>(
+ bhe_unknowns_index, bhe_unknowns_index) += _R_matrix;
+
+ // add the R_pi_s matrix to local_K
+ local_K
+ .template block<bhe_unknowns_size, temperature_size>(bhe_unknowns_index,
+ temperature_index)
+ .noalias() += _R_pi_s_matrix;
+ local_K
+ .template block<temperature_size, bhe_unknowns_size>(temperature_index,
+ bhe_unknowns_index)
+ .noalias() += _R_pi_s_matrix.transpose();
+
+ // add the R_s matrix to local_K
+ local_K
+ .template block<temperature_size, temperature_size>(temperature_index,
+ temperature_index)
+ .noalias() += 2.0 * _R_s_matrix;
+
+ // debugging
+ // 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/HeatTransportBHELocalAssemblerBHE.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h
index 0a904d254e2..ed74e2601ef 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE.h
@@ -23,18 +23,29 @@ namespace ProcessLib
{
namespace HeatTransportBHE
{
-template <typename ShapeFunction, typename IntegrationMethod, int BHE_Dim>
+template <typename ShapeFunction, typename IntegrationMethod, typename BHEType>
class HeatTransportBHELocalAssemblerBHE
: public HeatTransportBHELocalAssemblerInterface
{
+ static constexpr int bhe_unknowns = BHEType::number_of_unknowns;
+ static constexpr int single_bhe_unknowns_size = ShapeFunction::NPOINTS;
+ static constexpr int temperature_size = ShapeFunction::NPOINTS;
+ static constexpr int temperature_index = 0;
+ static constexpr int bhe_unknowns_size =
+ single_bhe_unknowns_size * bhe_unknowns;
+ static constexpr int bhe_unknowns_index = ShapeFunction::NPOINTS;
+ static constexpr int local_matrix_size =
+ temperature_size + bhe_unknowns_size;
+
public:
- using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, BHE_Dim>;
+ using ShapeMatricesType =
+ ShapeMatrixPolicyType<ShapeFunction, 3 /* GlobalDim */>;
using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
// Using dynamic size, because the number of unknowns in the BHE is runtime
// value.
using BheLocalMatrixType =
- typename ShapeMatricesType::template MatrixType<Eigen::Dynamic,
- Eigen::Dynamic>;
+ typename ShapeMatricesType::template MatrixType<local_matrix_size,
+ local_matrix_size>;
HeatTransportBHELocalAssemblerBHE(
HeatTransportBHELocalAssemblerBHE const&) = delete;
HeatTransportBHELocalAssemblerBHE(HeatTransportBHELocalAssemblerBHE&&) =
@@ -42,7 +53,7 @@ public:
HeatTransportBHELocalAssemblerBHE(
MeshLib::Element const& e,
- std::vector<unsigned> const& dofIndex_to_localIndex,
+ BHEType const& bhe,
bool const is_axially_symmetric,
unsigned const integration_order,
HeatTransportBHEProcessData& process_data);
@@ -52,14 +63,6 @@ public:
std::vector<double>& /*local_K_data*/,
std::vector<double>& /*local_b_data*/) override;
- void preTimestepConcrete(std::vector<double> const& /*local_x*/,
- double const /*t*/,
- double const /*delta_t*/) override
- {
- }
-
- void postTimestepConcrete(std::vector<double> const& /*local_x*/) override;
-
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
@@ -79,17 +82,25 @@ private:
IntegrationMethod _integration_method;
+ BHEType const& _bhe;
+
std::size_t const element_id;
SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
- BheLocalMatrixType _R_matrix;
+ typename ShapeMatricesType::template MatrixType<bhe_unknowns_size,
+ bhe_unknowns_size>
+ _R_matrix;
- BheLocalMatrixType _R_s_matrix;
+ typename ShapeMatricesType::template MatrixType<temperature_size,
+ temperature_size>
+ _R_s_matrix;
- BheLocalMatrixType _R_pi_s_matrix;
+ typename ShapeMatricesType::template MatrixType<bhe_unknowns_size,
+ temperature_size>
+ _R_pi_s_matrix;
};
} // namespace HeatTransportBHE
} // namespace ProcessLib
-#include "HeatTransportBHELocalAssemblerBHE_impl.h"
+#include "HeatTransportBHELocalAssemblerBHE-impl.h"
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE_impl.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE_impl.h
deleted file mode 100644
index 5d5d8c79caa..00000000000
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE_impl.h
+++ /dev/null
@@ -1,223 +0,0 @@
-/**
- * \copyright
- * Copyright (c) 2012-2018, 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 "HeatTransportBHELocalAssemblerBHE.h"
-#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
-#include "ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataBHE.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
-namespace ProcessLib
-{
-namespace HeatTransportBHE
-{
-using namespace BHE;
-
-template <typename ShapeFunction, typename IntegrationMethod, int BHE_Dim>
-HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod, BHE_Dim>::
- HeatTransportBHELocalAssemblerBHE(
- MeshLib::Element const& e,
- std::vector<unsigned> const& dofIndex_to_localIndex,
- bool const is_axially_symmetric,
- unsigned const integration_order,
- HeatTransportBHEProcessData& process_data)
- : HeatTransportBHELocalAssemblerInterface(
- ShapeFunction::NPOINTS * BHE_Dim, // no intersection
- dofIndex_to_localIndex),
- _process_data(process_data),
- _integration_method(integration_order),
- element_id(e.getID())
-{
- // need to make sure that the BHE elements are one-dimensional
- assert(e.getDimension() == 1);
-
- unsigned const n_integration_points =
- _integration_method.getNumberOfPoints();
-
- _ip_data.reserve(n_integration_points);
- _secondary_data.N.resize(n_integration_points);
-
- auto const shape_matrices =
- initShapeMatrices<ShapeFunction, ShapeMatricesType, IntegrationMethod,
- BHE_Dim>(e, is_axially_symmetric,
- _integration_method);
-
- auto mat_id = (*_process_data._mesh_prop_materialIDs)[e.getID()];
- auto BHE_id = _process_data._map_materialID_to_BHE_ID[mat_id];
-
- SpatialPosition x_position;
- x_position.setElementID(element_id);
-
- // ip data initialization
- for (unsigned ip = 0; ip < n_integration_points; ip++)
- {
- x_position.setIntegrationPoint(ip);
-
- // create the class IntegrationPointDataBHE in place
- _ip_data.emplace_back(*(_process_data._vec_BHE_property[BHE_id]));
- auto const& sm = shape_matrices[ip];
- auto& ip_data = _ip_data[ip];
- ip_data.integration_weight =
- _integration_method.getWeightedPoint(ip).getWeight() *
- sm.integralMeasure * sm.detJ;
- ip_data.N = sm.N;
- ip_data.dNdx = sm.dNdx;
-
- _secondary_data.N[ip] = sm.N;
- }
-
- const int BHE_n_unknowns = _ip_data[0]._bhe_instance.getNumUnknowns();
- const int NumRMatrixRows = ShapeFunction::NPOINTS * BHE_n_unknowns;
- _R_matrix.setZero(NumRMatrixRows, NumRMatrixRows);
- _R_pi_s_matrix.setZero(NumRMatrixRows, ShapeFunction::NPOINTS);
- _R_s_matrix.setZero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
- // formulate the local BHE R matrix
- Eigen::MatrixXd matBHE_loc_R =
- Eigen::MatrixXd::Zero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
- for (int idx_bhe_unknowns = 0; idx_bhe_unknowns < BHE_n_unknowns;
- idx_bhe_unknowns++)
- {
- matBHE_loc_R.setZero();
- // Loop over Gauss points
- for (unsigned ip = 0; ip < n_integration_points; ip++)
- {
- x_position.setIntegrationPoint(ip);
-
- auto const& N = _ip_data[ip].N;
- auto const& w = _ip_data[ip].integration_weight;
-
- // get coefficient of R matrix for corresponding BHE.
- auto R_coeff = _process_data._vec_BHE_property[BHE_id]
- ->getBoundaryHeatExchangeCoeff(idx_bhe_unknowns);
-
- // calculate mass matrix for current unknown
- matBHE_loc_R += N.transpose() * R_coeff * N * w;
- } // end of loop over integration point
-
- // The following assembly action is according to Diersch (2013) FEFLOW
- // book please refer to M.127 and M.128 on page 955 and 956
- _process_data._vec_BHE_property[BHE_id]->setRMatrices(
- idx_bhe_unknowns, ShapeFunction::NPOINTS, matBHE_loc_R, _R_matrix,
- _R_pi_s_matrix, _R_s_matrix);
- } // end of loop over BHE unknowns
-
- // debugging
- // std::string sep = "\n----------------------------------------\n";
- // Eigen::IOFormat CleanFmt(4, 0, ", ", "\n", "[", "]");
- // std::cout << "_R_matrix: \n" << sep;
- // std::cout << _R_matrix.format(CleanFmt) << sep;
- // std::cout << "_R_s_matrix: \n" << sep;
- // std::cout << _R_s_matrix.format(CleanFmt) << sep;
- // std::cout << "_R_pi_s_matrix: \n" << sep;
- // std::cout << _R_pi_s_matrix.format(CleanFmt) << sep;
-}
-
-template <typename ShapeFunction, typename IntegrationMethod, int BHE_Dim>
-void HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod,
- BHE_Dim>::
- assemble(
- double const /*t*/, std::vector<double> const& local_x,
- std::vector<double>& local_M_data, std::vector<double>& local_K_data,
- std::vector<double>& /*local_b_data*/) // local b vector is not touched
-{
- auto const local_matrix_size = local_x.size();
- const int BHE_n_unknowns = _ip_data[0]._bhe_instance.getNumUnknowns();
- // plus one because the soil temperature is included in local_x
- assert(local_matrix_size == ShapeFunction::NPOINTS * (BHE_n_unknowns + 1));
-
- auto local_M = MathLib::createZeroedMatrix<BheLocalMatrixType>(
- local_M_data, local_matrix_size, local_matrix_size);
- auto local_K = MathLib::createZeroedMatrix<BheLocalMatrixType>(
- local_K_data, local_matrix_size, local_matrix_size);
-
- unsigned const n_integration_points =
- _integration_method.getNumberOfPoints();
-
- SpatialPosition x_position;
- x_position.setElementID(element_id);
-
- int shift = 0;
- static const int local_BHE_matrix_size =
- ShapeFunction::NPOINTS * BHE_n_unknowns;
- static const int shift_start = ShapeFunction::NPOINTS;
-
- // the mass and conductance matrix terms
- for (unsigned ip = 0; ip < n_integration_points; ip++)
- {
- x_position.setIntegrationPoint(ip);
- auto& ip_data = _ip_data[ip];
-
- auto const& w = ip_data.integration_weight;
- auto const& N = ip_data.N;
- auto const& dNdx = ip_data.dNdx;
-
- // looping over all unknowns.
- for (int idx_bhe_unknowns = 0; idx_bhe_unknowns < BHE_n_unknowns;
- idx_bhe_unknowns++)
- {
- // get coefficient of mass from corresponding BHE.
- auto& mass_coeff = ip_data._vec_mass_coefficients[idx_bhe_unknowns];
- auto& laplace_mat = ip_data._vec_mat_Laplace[idx_bhe_unknowns];
- auto& advection_vec =
- ip_data._vec_Advection_vectors[idx_bhe_unknowns];
-
- // calculate shift.
- shift = shift_start + ShapeFunction::NPOINTS * idx_bhe_unknowns;
- // local M
- local_M
- .template block<ShapeFunction::NPOINTS, ShapeFunction::NPOINTS>(
- shift, shift)
- .noalias() += N.transpose() * mass_coeff * N * w;
-
- // local K
- // laplace part
- local_K
- .block(shift, shift, ShapeFunction::NPOINTS,
- ShapeFunction::NPOINTS)
- .noalias() += dNdx.transpose() * laplace_mat * dNdx * w;
- // advection part
- local_K
- .block(shift, shift, ShapeFunction::NPOINTS,
- ShapeFunction::NPOINTS)
- .noalias() +=
- N.transpose() * advection_vec.transpose() * dNdx * w;
- }
- }
-
- // add the R matrix to local_K
- local_K.block(shift_start, shift_start, local_BHE_matrix_size,
- local_BHE_matrix_size) += _R_matrix;
-
- // add the R_pi_s matrix to local_K
- local_K.block(shift_start, 0, local_BHE_matrix_size, shift_start) +=
- _R_pi_s_matrix;
- local_K.block(0, shift_start, shift_start, local_BHE_matrix_size) +=
- _R_pi_s_matrix.transpose();
-
- // add the R_s matrix to local_K
- local_K.block(0, 0, shift_start, shift_start) += 2.0 * _R_s_matrix;
-
- // debugging
- // 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;
-}
-
-template <typename ShapeFunction, typename IntegrationMethod, int BHE_Dim>
-void HeatTransportBHELocalAssemblerBHE<
- ShapeFunction, IntegrationMethod,
- BHE_Dim>::postTimestepConcrete(std::vector<double> const& /*local_x*/)
-{
-}
-} // namespace HeatTransportBHE
-} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil_impl.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil-impl.h
similarity index 74%
rename from ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil_impl.h
rename to ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil-impl.h
index ffe4d351484..058044c0dac 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil_impl.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil-impl.h
@@ -9,21 +9,14 @@
#pragma once
-#include "HeatTransportBHELocalAssemblerSoil.h"
-
#include <valarray>
#include <vector>
-#include <Eigen/Eigen>
-
#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
-
#include "NumLib/Fem/ShapeMatrixPolicy.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
#include "ProcessLib/HeatTransportBHE/HeatTransportBHEProcessData.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
-// #include "IntegrationPointDataMatrix.h"
#include "HeatTransportBHEProcessAssemblerInterface.h"
#include "SecondaryData.h"
@@ -31,22 +24,16 @@ namespace ProcessLib
{
namespace HeatTransportBHE
{
-template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
-HeatTransportBHELocalAssemblerSoil<ShapeFunction,
- IntegrationMethod,
- GlobalDim>::
+template <typename ShapeFunction, typename IntegrationMethod>
+HeatTransportBHELocalAssemblerSoil<ShapeFunction, IntegrationMethod>::
HeatTransportBHELocalAssemblerSoil(
MeshLib::Element const& e,
- std::vector<unsigned> const& dofIndex_to_localIndex,
bool const is_axially_symmetric,
unsigned const integration_order,
HeatTransportBHEProcessData& process_data)
- : HeatTransportBHELocalAssemblerInterface(
- ShapeFunction::NPOINTS * GlobalDim, dofIndex_to_localIndex),
- _process_data(process_data),
+ : _process_data(process_data),
_integration_method(integration_order),
- element_id(e.getID()),
- _is_axially_symmetric(is_axially_symmetric)
+ element_id(e.getID())
{
unsigned const n_integration_points =
_integration_method.getNumberOfPoints();
@@ -57,7 +44,7 @@ HeatTransportBHELocalAssemblerSoil<ShapeFunction,
_shape_matrices = initShapeMatrices<ShapeFunction,
ShapeMatricesType,
IntegrationMethod,
- GlobalDim>(
+ 3 /* GlobalDim */>(
e, is_axially_symmetric, _integration_method);
SpatialPosition x_position;
@@ -69,27 +56,23 @@ HeatTransportBHELocalAssemblerSoil<ShapeFunction,
x_position.setIntegrationPoint(ip);
// create the class IntegrationPointDataBHE in place
- _ip_data.emplace_back();
auto const& sm = _shape_matrices[ip];
- auto& ip_data = _ip_data[ip];
double const w = _integration_method.getWeightedPoint(ip).getWeight() *
sm.integralMeasure * sm.detJ;
- ip_data.NTN_product_times_w = sm.N.transpose() * sm.N * w;
- ip_data.dNdxTdNdx_product_times_w = sm.dNdx.transpose() * sm.dNdx * w;
+ _ip_data.push_back(
+ {sm.N.transpose() * sm.N * w, sm.dNdx.transpose() * sm.dNdx * w});
_secondary_data.N[ip] = sm.N;
}
}
-template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
-void HeatTransportBHELocalAssemblerSoil<
- ShapeFunction,
- IntegrationMethod,
- GlobalDim>::assemble(double const t,
- std::vector<double> const& local_x,
- std::vector<double>& local_M_data,
- std::vector<double>& local_K_data,
- std::vector<double>& /*local_b_data*/)
+template <typename ShapeFunction, typename IntegrationMethod>
+void HeatTransportBHELocalAssemblerSoil<ShapeFunction, IntegrationMethod>::
+ assemble(double const t,
+ std::vector<double> const& local_x,
+ std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data,
+ std::vector<double>& /*local_b_data*/)
{
assert(local_x.size() == ShapeFunction::NPOINTS);
(void)local_x; // Avoid unused arg warning.
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil.h
index f556f807109..51088a08481 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil.h
@@ -26,12 +26,13 @@ namespace ProcessLib
{
namespace HeatTransportBHE
{
-template <typename ShapeFunction, typename IntegrationMethod, int GlobalDim>
+template <typename ShapeFunction, typename IntegrationMethod>
class HeatTransportBHELocalAssemblerSoil
: public HeatTransportBHELocalAssemblerInterface
{
public:
- using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+ using ShapeMatricesType =
+ ShapeMatrixPolicyType<ShapeFunction, 3 /* GlobalDim */>;
using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
@@ -43,7 +44,6 @@ public:
HeatTransportBHELocalAssemblerSoil(
MeshLib::Element const& e,
- std::vector<unsigned> const& dofIndex_to_localIndex,
bool is_axially_symmetric,
unsigned const integration_order,
HeatTransportBHEProcessData& process_data);
@@ -77,11 +77,9 @@ private:
std::size_t const element_id;
- bool const _is_axially_symmetric;
-
SecondaryData<typename ShapeMatrices::ShapeType> _secondary_data;
};
} // namespace HeatTransportBHE
} // namespace ProcessLib
-#include "HeatTransportBHELocalAssemblerSoil_impl.h"
+#include "HeatTransportBHELocalAssemblerSoil-impl.h"
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHEProcessAssemblerInterface.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHEProcessAssemblerInterface.h
index d0390c17999..34f0964fd2d 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHEProcessAssemblerInterface.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHEProcessAssemblerInterface.h
@@ -9,11 +9,6 @@
#pragma once
-#include <utility>
-#include <vector>
-
-#include "BaseLib/Error.h"
-
#include "NumLib/Extrapolation/ExtrapolatableElement.h"
#include "ProcessLib/LocalAssemblerInterface.h"
@@ -25,22 +20,6 @@ class HeatTransportBHELocalAssemblerInterface
: public ProcessLib::LocalAssemblerInterface,
public NumLib::ExtrapolatableElement
{
-public:
- // The following function is necessary, because the BHE or Soil assemblers
- // create their local_x memory based on the passed on dofIndex_to_localIndex
- // Please keep it unchanged.
- HeatTransportBHELocalAssemblerInterface(std::size_t /*n_local_size*/,
- std::vector<unsigned>
- dofIndex_to_localIndex)
- : _dofIndex_to_localIndex(std::move(dofIndex_to_localIndex))
- {
- }
-
-private:
- // this dofTalbe must be kept here.
- // When initializing the assembler, this will be needed to
- // initialize the memory for local assemblers.
- std::vector<unsigned> const _dofIndex_to_localIndex;
};
} // namespace HeatTransportBHE
} // namespace ProcessLib
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataBHE.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataBHE.h
index 3f09811c486..983abb0f92a 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataBHE.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataBHE.h
@@ -9,64 +9,16 @@
#pragma once
-#include <Eigen/Eigen>
-
-#include "ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h"
-
namespace ProcessLib
{
namespace HeatTransportBHE
{
-using namespace BHE;
-
template <typename ShapeMatrixType>
struct IntegrationPointDataBHE final
{
- explicit IntegrationPointDataBHE(BHEAbstract& bhe_instance)
- : _bhe_instance(bhe_instance)
- {
- // depending on the type of BHE
- const int unknown_size = _bhe_instance.getNumUnknowns();
- // initialization
- _vec_mass_coefficients.resize(unknown_size);
- for (int i = 0; i < unknown_size; i++)
- {
- Eigen::MatrixXd mat_laplace(3, 3);
- mat_laplace.setZero();
- _vec_mat_Laplace.push_back(mat_laplace);
- Eigen::VectorXd vec_advection(3);
- vec_advection.setZero();
- _vec_Advection_vectors.push_back(vec_advection);
- }
-
- // set parameter values
- for (int j = 0; j < unknown_size; j++)
- {
- // mass matrix coefficients
- _vec_mass_coefficients[j] = _bhe_instance.getMassCoeff(j);
- // laplace matrix
- _bhe_instance.getLaplaceMatrix(j, _vec_mat_Laplace[j]);
- // advection vector
- _bhe_instance.getAdvectionVector(j, _vec_Advection_vectors[j]);
- }
- }
-
- BHEAbstract& _bhe_instance;
-
- typename ShapeMatrixType::NodalRowVectorType N;
- typename ShapeMatrixType::GlobalDimNodalMatrixType dNdx;
- double integration_weight;
-
- // mass coefficients, length depending on the type of BHE
- std::vector<double> _vec_mass_coefficients;
-
- // Laplace matrices
- std::vector<Eigen::MatrixXd> _vec_mat_Laplace;
-
- // Advection vectors
- std::vector<Eigen::VectorXd> _vec_Advection_vectors;
-
- void pushBackState() {}
+ typename ShapeMatrixType::NodalRowVectorType const N;
+ typename ShapeMatrixType::GlobalDimNodalMatrixType const dNdx;
+ double const integration_weight;
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataSoil.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataSoil.h
index 04e5c338af9..f67e762730d 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataSoil.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/IntegrationPointDataSoil.h
@@ -9,9 +9,6 @@
#pragma once
-#include <memory>
-#include <vector>
-
namespace ProcessLib
{
namespace HeatTransportBHE
@@ -19,8 +16,8 @@ namespace HeatTransportBHE
template <typename NodalMatrixType>
struct IntegrationPointDataSoil final
{
- NodalMatrixType NTN_product_times_w;
- NodalMatrixType dNdxTdNdx_product_times_w;
+ NodalMatrixType const NTN_product_times_w;
+ NodalMatrixType const dNdxTdNdx_product_times_w;
EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
};
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h
index ac08d79591b..9332acf884f 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/LocalDataInitializer.h
@@ -19,7 +19,7 @@
#include "MeshLib/Elements/Elements.h"
#include "NumLib/DOF/LocalToGlobalIndexMap.h"
#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
-#include "ProcessLib/HeatTransportBHE/BHE/BHEAbstract.h"
+#include "ProcessLib/HeatTransportBHE/BHE/BHETypes.h"
#ifndef OGS_MAX_ELEMENT_DIM
static_assert(false, "The macro OGS_MAX_ELEMENT_DIM is undefined.");
@@ -57,15 +57,6 @@ static_assert(false, "The macro OGS_MAX_ELEMENT_ORDER is undefined.");
#endif
// Dependent element types.
-// Faces of tets, pyramids and prisms are triangles
-#define ENABLED_ELEMENT_TYPE_TRI \
- ((ENABLED_ELEMENT_TYPE_SIMPLEX) | (ENABLED_ELEMENT_TYPE_PYRAMID) | \
- (ENABLED_ELEMENT_TYPE_PRISM))
-// Faces of hexes, pyramids and prisms are quads
-#define ENABLED_ELEMENT_TYPE_QUAD \
- ((ENABLED_ELEMENT_TYPE_CUBOID) | (ENABLED_ELEMENT_TYPE_PYRAMID) | \
- (ENABLED_ELEMENT_TYPE_PRISM))
-
// All enabled element types
#define OGS_ENABLED_ELEMENTS \
((ENABLED_ELEMENT_TYPE_SIMPLEX) | (ENABLED_ELEMENT_TYPE_CUBOID) | \
@@ -82,22 +73,11 @@ static_assert(false, "The macro OGS_MAX_ELEMENT_ORDER is undefined.");
#include "NumLib/Fem/ShapeFunction/ShapeTet4.h"
#endif
-#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0
-#include "NumLib/Fem/ShapeFunction/ShapeTri3.h"
-#include "NumLib/Fem/ShapeFunction/ShapeTri6.h"
-#endif
-
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0
#include "NumLib/Fem/ShapeFunction/ShapeHex20.h"
#include "NumLib/Fem/ShapeFunction/ShapeHex8.h"
#endif
-#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0
-#include "NumLib/Fem/ShapeFunction/ShapeQuad4.h"
-#include "NumLib/Fem/ShapeFunction/ShapeQuad8.h"
-#include "NumLib/Fem/ShapeFunction/ShapeQuad9.h"
-#endif
-
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_PRISM) != 0
#include "NumLib/Fem/ShapeFunction/ShapePrism15.h"
#include "NumLib/Fem/ShapeFunction/ShapePrism6.h"
@@ -118,9 +98,11 @@ namespace HeatTransportBHE
/// For example for MeshLib::Quad a local assembler data with template argument
/// NumLib::ShapeQuad4 is created.
template <typename LocalAssemblerInterface,
- template <typename, typename, int> class LocalAssemblerDataSoil,
- template <typename, typename, int> class LocalAssemblerDataBHE,
- int GlobalDim, typename... ConstructorArgs>
+ template <typename, typename>
+ class LocalAssemblerDataSoil,
+ template <typename, typename, typename>
+ class LocalAssemblerDataBHE,
+ typename... ConstructorArgs>
class LocalDataInitializer final
{
public:
@@ -132,27 +114,12 @@ public:
{
// REMARKS: At the moment, only a 3D mesh (soil) with 1D elements (BHE)
// are supported.
-
-#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 && \
- OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 1
- _builder[std::type_index(typeid(MeshLib::Quad))] =
- makeLocalAssemblerBuilder<NumLib::ShapeQuad4>();
-#endif
-
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 && \
OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 1
_builder[std::type_index(typeid(MeshLib::Hex))] =
makeLocalAssemblerBuilder<NumLib::ShapeHex8>();
#endif
-#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_QUAD) != 0 && \
- OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
- _builder[std::type_index(typeid(MeshLib::Quad8))] =
- makeLocalAssemblerBuilder<NumLib::ShapeQuad8>();
- _builder[std::type_index(typeid(MeshLib::Quad9))] =
- makeLocalAssemblerBuilder<NumLib::ShapeQuad9>();
-#endif
-
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_CUBOID) != 0 && \
OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
_builder[std::type_index(typeid(MeshLib::Hex20))] =
@@ -160,25 +127,12 @@ public:
#endif
// /// Simplices ////////////////////////////////////////////////
-
-#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 && \
- OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 1
- _builder[std::type_index(typeid(MeshLib::Tri))] =
- makeLocalAssemblerBuilder<NumLib::ShapeTri3>();
-#endif
-
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 && \
OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 1
_builder[std::type_index(typeid(MeshLib::Tet))] =
makeLocalAssemblerBuilder<NumLib::ShapeTet4>();
#endif
-#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_TRI) != 0 && \
- OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 2
- _builder[std::type_index(typeid(MeshLib::Tri6))] =
- makeLocalAssemblerBuilder<NumLib::ShapeTri6>();
-#endif
-
#if (OGS_ENABLED_ELEMENTS & ENABLED_ELEMENT_TYPE_SIMPLEX) != 0 && \
OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
_builder[std::type_index(typeid(MeshLib::Tet10))] =
@@ -216,12 +170,12 @@ public:
#if OGS_MAX_ELEMENT_DIM >= 2 && OGS_MAX_ELEMENT_ORDER >= 1
_builder[std::type_index(typeid(MeshLib::Line))] =
- makeLocalAssemblerBuilder<NumLib::ShapeLine2>();
+ makeLocalAssemblerBuilderBHE<NumLib::ShapeLine2>();
#endif
#if OGS_MAX_ELEMENT_DIM >= 3 && OGS_MAX_ELEMENT_ORDER >= 2
_builder[std::type_index(typeid(MeshLib::Line3))] =
- makeLocalAssemblerBuilder<NumLib::ShapeLine3>();
+ makeLocalAssemblerBuilderBHE<NumLib::ShapeLine3>();
#endif
}
@@ -230,15 +184,12 @@ public:
/// \attention
/// The index \c id is not necessarily the mesh item's id. Especially when
/// having multiple meshes it will differ from the latter.
- void operator()(
- std::size_t const id,
- MeshLib::Element const& mesh_item,
- LADataIntfPtr& data_ptr,
- const std::vector<
- std::vector<std::size_t>>& /*vec_ele_connected_BHE_IDs*/,
- const std::vector<
- std::unique_ptr<BHE::BHEAbstract>>& /*vec_BHE_property*/,
- ConstructorArgs&&... args) const
+ void operator()(std::size_t const /*id*/,
+ MeshLib::Element const& mesh_item,
+ LADataIntfPtr& data_ptr,
+ std::unordered_map<std::size_t, BHE::BHETypes*> const&
+ element_to_bhe_map,
+ ConstructorArgs&&... args) const
{
auto const type_idx = std::type_index(typeid(mesh_item));
auto const it = _builder.find(type_idx);
@@ -253,49 +204,16 @@ public:
type_idx.name());
}
- std::vector<unsigned> dofIndex_to_localIndex;
-
- // Create customed dof table when it is a BHE element.
- if (mesh_item.getDimension() == 1)
- {
- // this is a BHE element
- auto const n_local_dof = _dof_table.getNumberOfElementDOF(id);
- dofIndex_to_localIndex.resize(n_local_dof);
- unsigned dof_id = 0;
- unsigned local_id = 0;
- for (auto const i : _dof_table.getElementVariableIDs(id))
- {
- auto const n_global_components =
- _dof_table.getNumberOfElementComponents(i);
- for (unsigned j = 0; j < n_global_components; j++)
- {
- auto const& ms = _dof_table.getMeshSubset(i, j);
- auto const mesh_id = ms.getMeshID();
- for (unsigned k = 0; k < mesh_item.getNumberOfNodes(); k++)
- {
- MeshLib::Location l(mesh_id,
- MeshLib::MeshItemType::Node,
- mesh_item.getNodeIndex(k));
- auto const global_index =
- _dof_table.getGlobalIndex(l, i, j);
- if (global_index != NumLib::MeshComponentMap::nop)
- {
- dofIndex_to_localIndex[dof_id++] = local_id;
- }
- local_id++;
- }
- }
- }
- }
-
- data_ptr = it->second(mesh_item, dofIndex_to_localIndex,
+ data_ptr = it->second(mesh_item,
+ element_to_bhe_map,
std::forward<ConstructorArgs>(args)...);
}
private:
using LADataBuilder = std::function<LADataIntfPtr(
MeshLib::Element const& e,
- std::vector<unsigned> const& dofIndex_to_localIndex,
+ std::unordered_map<std::size_t, BHE::BHETypes*> const&
+ element_to_bhe_map,
ConstructorArgs&&...)>;
template <typename ShapeFunction>
@@ -305,87 +223,53 @@ private:
// local assembler builder implementations.
template <typename ShapeFunction>
using LADataSoil =
- LocalAssemblerDataSoil<ShapeFunction, IntegrationMethod<ShapeFunction>,
- GlobalDim>;
+ LocalAssemblerDataSoil<ShapeFunction, IntegrationMethod<ShapeFunction>>;
- template <typename ShapeFunction>
- using LADataBHE =
- LocalAssemblerDataBHE<ShapeFunction, IntegrationMethod<ShapeFunction>,
- GlobalDim>;
- /// A helper forwarding to the correct version of makeLocalAssemblerBuilder
- /// depending whether the global dimension is less than the shape function's
- /// dimension or not.
template <typename ShapeFunction>
static LADataBuilder makeLocalAssemblerBuilder()
- {
- return makeLocalAssemblerBuilder<ShapeFunction>(
- static_cast<std::integral_constant<
- bool, (GlobalDim >= ShapeFunction::DIM)>*>(nullptr));
- }
-
- /// Mapping of element types to local assembler constructors.
- std::unordered_map<std::type_index, LADataBuilder> _builder;
-
- NumLib::LocalToGlobalIndexMap const& _dof_table;
-
-private:
- /// Generates a function that creates a new LocalAssembler of type
- /// LAData<ShapeFunction>. Only functions with shape function's dimension
- /// less or equal to the global dimension are instantiated, e.g. following
- /// combinations of shape functions and global dimensions: (Line2, 1),
- /// (Line2, 2), (Line2, 3), (Hex20, 3) but not (Hex20, 2) or (Hex20, 1).
- template <typename ShapeFunction>
- static LADataBuilder makeLocalAssemblerBuilder(std::true_type*)
{
return [](MeshLib::Element const& e,
- std::vector<unsigned> const& dofIndex_to_localIndex,
+ std::unordered_map<std::size_t, BHE::BHETypes*> const&
+ /* unused */,
ConstructorArgs&&... args) -> LADataIntfPtr {
- if (e.getDimension() == GlobalDim) // soil elements
+ if (e.getDimension() == 3) // soil elements
{
return LADataIntfPtr{new LADataSoil<ShapeFunction>{
- e, dofIndex_to_localIndex,
- std::forward<ConstructorArgs>(args)...}};
+ e, std::forward<ConstructorArgs>(args)...}};
}
return nullptr;
};
}
- template <>
- static LADataBuilder makeLocalAssemblerBuilder<NumLib::ShapeLine2>(
- std::true_type*)
+ template <typename ShapeFunction, typename BHEType>
+ using LADataBHE = LocalAssemblerDataBHE<ShapeFunction,
+ IntegrationMethod<ShapeFunction>,
+ BHEType>;
+ template <typename ShapeFunction>
+ static LADataBuilder makeLocalAssemblerBuilderBHE()
{
return [](MeshLib::Element const& e,
- std::vector<unsigned> const& dofIndex_to_localIndex,
+ std::unordered_map<std::size_t, BHE::BHETypes*> const&
+ element_to_bhe_map,
ConstructorArgs&&... args) -> LADataIntfPtr {
- return LADataIntfPtr{new LADataBHE<NumLib::ShapeLine2>{
- // BHE elements
- e, dofIndex_to_localIndex,
- std::forward<ConstructorArgs>(args)...}};
- };
- }
+ auto& bhe = *element_to_bhe_map.at(e.getID());
- template <>
- static LADataBuilder makeLocalAssemblerBuilder<NumLib::ShapeLine3>(
- std::true_type*)
- {
- return [](MeshLib::Element const& e,
- std::vector<unsigned> const& dofIndex_to_localIndex,
- ConstructorArgs&&... args) -> LADataIntfPtr {
- return LADataIntfPtr{new LADataBHE<NumLib::ShapeLine3>{
- // BHE elements
- e, dofIndex_to_localIndex,
- std::forward<ConstructorArgs>(args)...}};
+ if (bhe.type() == typeid(BHE::BHE_1U))
+ {
+ return LADataIntfPtr{new LADataBHE<ShapeFunction, BHE::BHE_1U>{
+ e, boost::get<BHE::BHE_1U>(bhe),
+ std::forward<ConstructorArgs>(args)...}};
+ }
+ OGS_FATAL(
+ "Trying to create local assembler for an unknown BHE type.");
};
}
- /// Returns nullptr for shape functions whose dimensions are less than the
- /// global dimension.
- template <typename ShapeFunction>
- static LADataBuilder makeLocalAssemblerBuilder(std::false_type*)
- {
- return nullptr;
- }
+ /// Mapping of element types to local assembler constructors.
+ std::unordered_map<std::type_index, LADataBuilder> _builder;
+
+ NumLib::LocalToGlobalIndexMap const& _dof_table;
}; // namespace HeatTransportBHE
} // namespace HeatTransportBHE
} // namespace ProcessLib
@@ -394,6 +278,4 @@ private:
#undef ENABLED_ELEMENT_TYPE_CUBOID
#undef ENABLED_ELEMENT_TYPE_PYRAMID
#undef ENABLED_ELEMENT_TYPE_PRISM
-#undef ENABLED_ELEMENT_TYPE_TRI
-#undef ENABLED_ELEMENT_TYPE_QUAD
#undef OGS_ENABLED_ELEMENTS
--
GitLab