BHE_2U.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <Eigen/Core>
#include <optional>
 
#include "BHECommon.h"
#include "BHECommonUType.h"
#include "BaseLib/Error.h"
#include "FlowAndTemperatureControl.h"
#include "PipeConfigurationUType.h"
 
namespace ParameterLib
{
class SpatialPosition;
}
 
namespace ProcessLib
{
namespace HeatTransportBHE
{
namespace BHE
{
class BHE_2U final : public BHECommonUType
{
public:
    BHE_2U(BoreholeGeometry const& borehole,
           RefrigerantProperties const& refrigerant,
           GroutParameters const& grout,
           FlowAndTemperatureControl const& flowAndTemperatureControl,
           PipeConfigurationUType const& pipes,
           bool const use_python_bcs);
 
    static constexpr int number_of_unknowns = 8;
    static constexpr int number_of_grout_zones = 4;
    static constexpr int number_of_flow_legs = 4;

    std::array<double, number_of_flow_legs> flowLegs() const
    {
        return {+flow_velocity, +flow_velocity, -flow_velocity, -flow_velocity};
    }
 
    std::array<double, number_of_unknowns> pipeHeatCapacities() const;
 
    std::array<double, number_of_unknowns> pipeHeatConductions() const;
 
    std::array<Eigen::Vector3d, number_of_unknowns> pipeAdvectionVectors(
        Eigen::Vector3d const& elem_direction) const;

    std::vector<double> thermalResistances(
        ParameterLib::SpatialPosition const& pos) const;
 
    template <int NPoints,
              typename SingleUnknownMatrixType,
              typename RMatrixType,
              typename RPiSMatrixType,
              typename RSMatrixType>
    static 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)
    {
        switch (idx_bhe_unknowns)
        {
            case 0:  // PHI_fig
                R_matrix.block(0, 4 * NPoints, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;
                R_matrix.block(4 * NPoints, 0, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;  // R_i1
 
                R_matrix.block(2, 5 * NPoints, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;
                R_matrix.block(5 * NPoints, 2, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;  // R_i2
 
                R_matrix.block(0, 0, NPoints, NPoints) +=
                    1.0 * matBHE_loc_R;  // K_i1
                R_matrix.block(NPoints, NPoints, NPoints, NPoints) +=
                    1.0 * matBHE_loc_R;  // K_i2
                R_matrix.block(4 * NPoints, 4 * NPoints, NPoints, NPoints) +=
                    1.0 * matBHE_loc_R;
                R_matrix.block(5 * NPoints, 5 * NPoints, NPoints, NPoints) +=
                    1.0 * matBHE_loc_R;  // K_ig
                return;
            case 1:  // PHI_fog
                R_matrix.block(2 * NPoints, 6 * NPoints, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;
                R_matrix.block(6 * NPoints, 2 * NPoints, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;  // R_o1
                R_matrix.block(3 * NPoints, 7 * NPoints, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;
                R_matrix.block(7 * NPoints, 3 * NPoints, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;  // R_o2
 
                R_matrix.block(2 * NPoints, 2 * NPoints, NPoints, NPoints) +=
                    1.0 * matBHE_loc_R;  // K_o1
                R_matrix.block(3 * NPoints, 3 * NPoints, NPoints, NPoints) +=
                    1.0 * matBHE_loc_R;  // K_o2
                R_matrix.block(6 * NPoints, 6 * NPoints, NPoints, NPoints) +=
                    1.0 * matBHE_loc_R;
                R_matrix.block(7 * NPoints, 7 * NPoints, NPoints, NPoints) +=
                    1.0 * matBHE_loc_R;  // K_og
                return;
            case 2:  // PHI_gg_1
                R_matrix.block(4 * NPoints, 6 * NPoints, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;
                R_matrix.block(6 * NPoints, 4 * NPoints, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;
                R_matrix.block(4 * NPoints, 7 * NPoints, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;
                R_matrix.block(7 * NPoints, 4 * NPoints, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;
                R_matrix.block(5 * NPoints, 6 * NPoints, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;
                R_matrix.block(6 * NPoints, 5 * NPoints, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;
                R_matrix.block(5 * NPoints, 7 * NPoints, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;
                R_matrix.block(7 * NPoints, 5 * NPoints, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;  // R_g1
 
                R_matrix.block(4 * NPoints, 4 * NPoints, NPoints, NPoints) +=
                    2.0 * matBHE_loc_R;
                R_matrix.block(5 * NPoints, 5 * NPoints, NPoints, NPoints) +=
                    2.0 * matBHE_loc_R;  // K_ig
                R_matrix.block(6 * NPoints, 6 * NPoints, NPoints, NPoints) +=
                    2.0 * matBHE_loc_R;
                R_matrix.block(7 * NPoints, 7 * NPoints, NPoints, NPoints) +=
                    2.0 * matBHE_loc_R;  // K_og
                return;
            case 3:  // PHI_gg_2
                R_matrix.block(4 * NPoints, 5 * NPoints, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;
                R_matrix.block(5 * NPoints, 4 * NPoints, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;
                R_matrix.block(6 * NPoints, 7 * NPoints, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;
                R_matrix.block(7 * NPoints, 6 * NPoints, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;  // R_g2
 
                R_matrix.block(4 * NPoints, 4 * NPoints, NPoints, NPoints) +=
                    1.0 * matBHE_loc_R;
                R_matrix.block(5 * NPoints, 5 * NPoints, NPoints, NPoints) +=
                    1.0 * matBHE_loc_R;  // K_ig
                R_matrix.block(6 * NPoints, 6 * NPoints, NPoints, NPoints) +=
                    1.0 * matBHE_loc_R;
                R_matrix.block(7 * NPoints, 7 * NPoints, NPoints, NPoints) +=
                    1.0 * matBHE_loc_R;  // K_og
                return;
            case 4:  // PHI_gs
                R_s_matrix.template block<NPoints, NPoints>(0, 0).noalias() +=
                    1.0 * matBHE_loc_R;
 
                R_pi_s_matrix.block(4 * NPoints, 0, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;
                R_pi_s_matrix.block(5 * NPoints, 0, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;
                R_pi_s_matrix.block(6 * NPoints, 0, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;
                R_pi_s_matrix.block(7 * NPoints, 0, NPoints, NPoints) +=
                    -1.0 * matBHE_loc_R;
                R_matrix.block(4 * NPoints, 4 * NPoints, NPoints, NPoints) +=
                    1.0 * matBHE_loc_R;
                R_matrix.block(5 * NPoints, 5 * NPoints, NPoints, NPoints) +=
                    1.0 * matBHE_loc_R;  // K_ig
                R_matrix.block(6 * NPoints, 6 * NPoints, NPoints, NPoints) +=
                    1.0 * matBHE_loc_R;
                R_matrix.block(7 * NPoints, 7 * NPoints, NPoints, NPoints) +=
                    1.0 * matBHE_loc_R;  // K_og
                return;
            default:
                OGS_FATAL(
                    "BHE_2U::assembleRMatrices: unknown index {:d} "
                    "out of range.",
                    idx_bhe_unknowns);
        }
    }

    double updateFlowRateAndTemperature(double T_out, double current_time);
 
    static constexpr std::pair<int, int> inflow_outflow_bc_component_ids[] = {
        {0, 2}, {1, 3}};
 
    static std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>
    getBHEInflowDirichletBCNodesAndComponents(
        std::size_t const top_node_id,
        std::size_t const /*bottom_node_id*/,
        int const in_component_id);
 
    static std::optional<
        std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>>
    getBHEBottomDirichletBCNodesAndComponents(std::size_t const bottom_node_id,
                                              int const in_component_id,
                                              int const out_component_id);
 
public:
    std::array<double, number_of_unknowns> crossSectionAreas(
        ParameterLib::SpatialPosition const& pos) const;
 
    void updateHeatTransferCoefficients(double const flow_rate);
 
private:
    double cached_nu_ = 0.0;
 
    std::vector<double> calcThermalResistances(
        double const Nu, ParameterLib::SpatialPosition const& pos) const;
};
}  // namespace BHE
}  // namespace HeatTransportBHE
}  // namespace ProcessLib