Public API

initialize(options::SimulationOptions)

Precompute the objects of each TimeSuperpositionMethod that can be computed ahead of time and return the SimulationContainers of the required size.

simulate!(;options::SimulationOptions, operator, containers::SimulationContainers)

Run the simulation defined by options. At the end of simulation, containers.X will contain the results. containers should be the output of initialize.

operator should be a function that returns a BoreholeOperation and with signature operator(step, options, X):

• step::Int is the time step.
• options is the SimulationOptions object passed to simulate!.
• X is the matrix containers.X containing Tin, Tout, Tb, and q for each borehole.

simulate_steps!(;n, initial_step = nothing, options::SimulationOptions, operator, containers::SimulationContainers)

Similar to simulate!, but only run n steps, starting at the step initial_step, of the simulation defined by options. If initial_step is not specified, the index of first zero column of containers.X is taken as the initial step.

BoreholeNetwork(
    graph::SimpleDiGraph = SimpleDiGraph()
)

Represents the network of boreholes with a directed graph graph, whose edges represent the flow of fluid. The graph is expected to have number of vertices Nb+2, where Nb is the number of boreholes in the borefield. Vertexs number Nb+1 and Nb+2 represent the mass flow source and sink, respectively. All branches are expected to be connected to the source at the beggining and to the sink at the end.

Convenience constructors

• BoreholeNetwork(n::Int)

Create a network of n boreholes (and the source and sink) without any connections.

• all_parallel_network(n::Int)

Create a network of n boreholes connected in parallel.

• all_series_network(n::Int)

Create a network of n boreholes connected in series.

boreholes_in_branch(n::BoreholeNetwork; first_bh::Int)

Return all the boreholes in the branch that starts by borehole first_bh.

first_bhs_in_branch(network::BoreholeNetwork)

Return the first borehole in each branch of the borefield.

source(network::BoreholeNetwork)

Return the vertex index of the source.

sink(network::BoreholeNetwork)

Return the vertex index of the sink.

connect!(network::BoreholeNetwork, i::Int, j::Int)

Create a flow from borehole i to borehole j.

connect_to_source!(network::BoreholeNetwork, i::Int)

Create a connection injecting mass flow to borehole i.

connect_to_source!(network::BoreholeNetwork, I::Vector{Int})

Vector version of connecttosource!.

connect_to_sink!(network::BoreholeNetwork, i::Int)

Create a connection extracting mass flow from borehole i.

connect_to_sink!(network::BoreholeNetwork, I::Vector{Int})

Vector version of connecttosink!.

connect_in_series!(network::BoreholeNetwork, I::Vector{Int})

Connect all the boreholes in I, in order.

connect_in_parallel!(network::BoreholeNetwork, j::Int, I::Vector{Int})

Make a connection from borehole j to each of the boreholes in I.

BoreholeOperation{T <: Number}(
    valves::Dict{Int, Valve{T}}
    source_mass_flow::T
    network::BoreholeNetwork
)

Representation of the a state of operation of the borefield.

Arguments

• valves: A dicitonary specifying the divergences of the network. Each key is the borehole where the divergence occurs, and its corresponding value is a Valve object, describing the behaviour of the divergence.
• source_mass_flow: The total amount of mass flow injected into all branches.
• network: The representation of the pipe layout of the borefield, including direction of flow.

Convenience constructors

BoreholeOperation(;network::BoreholeNetwork, mass_flows::Vector{T}) Builds the operation corresponding to network sending mass flow mass_flows[i] in each branch identical. This method takes care of the initial valve from the source to the start of each branch.

Valve{T <: Number}(
    split::Dict{Int, T}
)

Models the behaviour of the fluid at pipe divergences. split is a dicitonary representing the divergence. Its keys represent the boreholes where the flow splits and its values are the proportion of flow going to each borehole.

valve(nodes::Vector{Int}, weights::Vector{T})

Return a Valve that splits the flow to each borehole nodes[i] in the proportion weights[i].

equal_valve(nodes::Vector{Int})

Return a Valve that splits the flow equally to all the boreholes nodes in the divergence.

absolute_valve(nodes::Vector{Int}, mass_flows::Vector{T})

Return a Valve that splits the flow such that each borehole nodes[i] receives an absolute amount of mass flow equal to mass_flows[i].

abstract type Operator

Interface for operation strategies.

Required functions:

• operate(::Operator, step, options, X)

The implementation of operate must return an instance of BoreholeOperation, specifying the network topology and the mass flows per branch for the current time step. step is the current time step, options contains the provided simulation options and X contains the time series up to the time step step-1 of the inlet fluid temperature, the outlet fluid temperature, the wall temperature and the heat extraction rate, for each borehole in the borefield. To get Tfin, Tfout, Tb, and q from X, use the functions extract_Tfin, extract_Tfout, extract_Tb, extract_q, or get them manually at the indices [1:2:2Nb, :], [2:2:2Nb, :], [2Nb+1:3Nb, :], [3Nb+1:end, :], respectively.

ConstantOperator{T <: Number} <: Operator (
    valves::Dict{Int, Valve{T}} = Dict{Int, Valve{Float64}}()
    mass_flow::T
)

Constant operation strategy. Its implementation of operate is: operate(op::ConstantOperator, step, options, X) = BoreholeOperation(op.valves, op.mass_flow, options.configurations[1])

Convenience constructors

ConstantOperator(network::BoreholeNetwork; mass_flows)

Builds the initial valve and total source mass flow needed to have mass flow equal to mass_flows[i] in branch i, according to the network network.

struct SimulationOptions{
                N <: Number,
                Tol <: Number,
                TSM <: TimeSuperpositionMethod,
                C <: Constraint,
                B <: Borefield, 
                M <: Medium, 
                BC <: BoundaryCondition,
                A <: Approximation,
                F <: Fluid
            }(
    method::TSM
    constraint::C
    borefield::B
    medium::M
    boundary_condition::BoundaryCondition = DirichletBoundaryCondition()
    approximation::A = MeanApproximation()
    fluid::Fluid{N} = Fluid(cpf=4182., name="INCOMP::MEA-20%")
    configurations::Vector{BoreholeNetwork}
    Δt
    Nt::Int
    atol::Tol = 0.
    rtol::Tol = sqrt(eps())
)

Specifies all the options for the simulation.

• method: time superposition method used to compute the response. Available options: ConvolutionMethod, NonHistoryMethod.
• constraint: constraint that the system must satisfy. Can be variable with time. Available options: HeatLoadConstraint, InletTempConstraint, TotalHeatLoadConstraint.
• borefield: describes the geometrical properties and the boreholes of the borefield on which the simulation will be performed. Available options: EqualBoreholesBorefield.
• medium: properties of the ground where the borefield is places. Available options: GroundMedium, FlowInPorousMedium.
• boundary_condition: boundary condition of the domain where the simulation is performed. Available options: NoBoundary, DirichletBoundaryCondition, NeumannBoundaryCondition.
• approximation: determines how the approximate value for each segment is computed. Available options: MeanApproximation, MidPointApproximation.
• fluid: properties of the fluid flowing through the hydraulic system.
• configurations: possible hydraulic topologies possible in the system, including reverse flow.
• Δt: time step used in the simulation.
• Nt: total amount of time steps of the simulation.
• atol: absolute tolerance used for the adaptive integration methods.
• rtol: relative tolerance used for the adaptive integration methods.

abstract type Fluid end

Interface for fluids.

Required functions:

• cpf(::Fluid): Return the scpecific heat capacity of the fluid.
• thermophysical_properties(::Fluid, T): Return the μ, ρ, cp, and k of the fluid at temperature T.

Water <: Fluid (
    stored_properties::ThermophysicalProperties{Float64}
)

Models water.

To initialize, use the convenience method: function Water() that will automatically compute stored_properties.

EthanolMix <: Fluid (
    stored_properties::ThermophysicalProperties{Float64}
)

Models a 20% ethanol and water mix.

To initialize, use the convenience method: function EthanolMix() that will automatically compute stored_properties.

abstract type Medium

Interface for mediums.

Required functions:

• get_λ(::Medium): Return the thermal conductivity of the medium.
• get_α(::Medium): Return the thermal diffusivity of the medium.
• get_T0(::Medium): Return the initial temperature of the medium.

GroundMedium{T <: Real} <: Medium @deftype T

Model pure conduction in the ground.

Arguments

• λ = 3.: ground conductivity
• α = 1e-6: ground thermal diffusivity
• C = λ/α: ground medium capacity
• T0 = 0.: initial ground temperature

FlowInPorousMedium{T <: Real} <: Medium @deftype T

Model a porous ground with a water flow.

Arguments

• λw = 0.6: water thermal conductivity
• λs = 2.: ground thermal conductivity
• Cw = 4.18*1e6:water thermal capacity
• Cs = 1.7*1e6: ground thermal capacity
• θ = 0.: angle of Darcy velocity
• Φ = 0.2: porosity
• λ = λs * (1-Φ) + λw*Φ: porous medium conductivity
• C = Cs * (1-Φ) + Cw*Φ: porous medium capacity
• α = λ/C: porous medium thermal diffusivity
• ux_in_meterperday = 1e-2: groundwater speed along the flow coordinate
• ux = ux_in_meterperday/(3600*24): groundwater speed in m/s
• vt = ux * Cw/C: porous medium darcy velocity
• T0 = 0.: initial ground temperature

abstract type Borefield

Interface for borefields.

Required functions

• n_boreholes(::Borefield): Return the amount of boreholes present in the borefield.
• get_H(::Borefield, i): Return the length of borehole i.
• get_rb(::Borefield, i): Return the radius of borehole i.
• segment_coordinates(::Borefield): Return a vector with the coordinates of each segment.
• internal_model_coeffs!(M, ::Borefield, medium, operation, fluid): Compute inplace in M the coefficients corresponding to the internal model equations, given the medium, fluid and operation in use in the system. Note that M is only a slice of Nb (number of boreholes) rows, provided as a view.
• internal_model_b!(b, ::Borefield): Compute inplace in b the independent terms corresponding to the internal model equations. Note that b is only a vector of length Nb (number of boreholes) rows, provided as a view.

EqualBoreholesBorefield{T <: Borehole, S <: Number} <: Borefield
EqualBoreholesBorefield(borehole_prototype::T, positions::Vector{Point2{S}}))

Model a borefield with boreholes all identical to the prototype borehole_prototype, placed at positions. Note that the length of positions determines the amount of boreholes in the field.

abstract type Borehole

Interface for boreholes.

Required functions:

• get_H(::Borehole): Return the length of the borehole.
• get_D(::Borehole): Return the burial depth of the borehole.
• get_rb(::Borehole): Return the radius of the borehole.
• uniform_Tb_coeffs(::Borehole, λ, mass_flow, Tref, fluid): Return the internal model coefficients for the resistance network between the pipes and the wall.

SingleUPipeBorehole{T <: Real} <: Borehole @deftype T
SingleUPipeBorehole(H, D)

Model a borehole with a single U-pipe with burial depth D and length H.

Arguments

• λg = 2.5: grout conductivity
• Cg = 2000. * 1550.: grout capacity
• αg = λg/Cg: grout thermal diffusivity
• rp = 0.02: pipe radius
• λp = 0.42: pipe material conductivity
• dpw = 0.0023: pipe thickness
• rpo = rp + dpw: equivalent pipe radius
• hp = 725.: heat transfer coefficient fluid to pipe
• pipe_position::NTuple{2, Tuple{T, T}} = [(0.03, 0.0), (-0.03, 0.0)]: positions of the downward and upward branches of the pipe. (0, 0) represents the center of the borehole.
• rb = 0.115/2: borehole radius

abstract type Constraint

Interface for constraints.

Required functions:

• constraints_coeffs!(M, ::Constraint, operation::BoreholeOperation): Compute inplace in M the coefficients corresponding to the constraints equations, given the current operation.network. Note that M is only a slice of Nbr (number of branches) rows, provided as a view.
• constraints_b!(b, ::Constraint, ::BoreholeOperation, step): Compute inplace in b the independent term corresponding to the constraints equations, given the current operation.network, at the time step step. Note that b is only a vector of length Nbr (number of branches), provided as a view.

HeatLoadConstraint(Q_tot::Matrix{T}){T <: Number} <: Constraint

Constrain the heat extracted per branch.

The heat constraint Q_tot must be a Matrix, whose column i are the loads per branch at the time step i. The amount of rows of Q_tot must equal to the amount of branches specified in BoreholeNetwork.

constant_HeatLoadConstraint(Q_tot::Vector{T}, Nt) where {T <: Number}

Convenience initializer for HeatLoadConstraint. It creates a constant heat load constraint through all the Nt time steps, where Q_tot are the heat load for each branch.

uniform_HeatLoadConstraint(Q_tot::Vector{T}, Nbr) where {T <: Number}

Convenience initializer for HeatLoadConstraint. It creates a uniform heat load constraint along all branches, where T_in are the heat load for each time step.

InletTempConstraint(T_in::Matrix{T}){T <: Number} <: Constraint

Constrain the inlet temperature of the first borehole in each branch.

The inlet temperature T_in must be a $N_{br} imes N_t$ Matrix, whose column i are the inlet temperatures per branch at the time step i. The amount of rows of T_in must equal to the amount of branches specified in BoreholeNetwork.

constant_InletTempConstraint(T_in::Vector{N}, Nt) where {N <: Number}

Convenience initializer for InletTempConstraint. It creates a constant inlet temperature constraint through all the Nt time steps, where T_in are the inlet temperatures for each branch.

uniform_InletTempConstraint(T_in::Vector{N}, Nbr) where {N <: Number}

Convenience initializer for InletTempConstraint. It creates a uniform inlet temperature constraint along all branches, where T_in are the inlet temperatures for each time step.

abstract type TimeSuperpositionMethod

Interface for time superposition methods.

Required functions:

• method_coeffs!(M, ::TimeSuperpositionMethod, options): Compute inplace in M the coefficients corresponding to the heat transfer equations, given options. Note that M is only a slice of Nbr (number of branches) rows, provided as a view.
• method_b!(b, ::TimeSuperpositionMethod, borefield, medium, step): Compute inplace in b the independent terms corresponding to the heat transfer equations, given the medium, at the given time step step. Note that b is only a vector of length Nbr (number of branches) rows, provided as a view.
• precompute_auxiliaries!(method::TimeSuperpositionMethod; options): Compute inplace in method the auxiliary quantities used in the simulation that can be performed ahead of time.
• update_auxiliaries!(::TimeSuperpositionMethod, X, borefield, step): Update inplace in method the auxiliaries after each time step step.

ConvolutionMethod{T} <: TimeSuperpositionMethod 
ConvolutionMethod()

Use the naïve convolution to compute the thermal response between boreholes. It should be initialized without arguments, but it contains the variables:

• g stores the unit response between each pair of boreholes evaluated at each time of the simulation.

It should be precomputed with initialize.

• q stores the heat extraction in each borehole at each time step. It is filled as the simulation runs.

NonHistoryMethod{T} <: TimeSuperpositionMethod 
NonHistoryMethod()

Use the non-history method to compute the thermal response between boreholes. See A non-history dependent temporal superposition algorithm for the finite line source solution. It should be initialized without arguments, but it contains the variables:

• F::Matrix{T}: each column contains the F function (encoding the load history) for each borehole. It is initially 0.
• ζ::Vector{T}: discretization nodes of the integration interval. Shared for all boreholes. Precomputed in initialize.
• w::Matrix{T}: weights of the ζ integration for each pair of boreholes. Precomputed in initialize.
• expΔt::Vector{T}: exp(-ζ^2*Δt). Precomputed in initialize.

This feature is experimental and might not work as expected in some cases.

abstract type BoundaryCondition

Interface for boundary conditions.

NoBoundary <: BoundaryCondition

Option to solve the problem in an infinite domain

DirichletBoundaryCondition <: BoundaryCondition

Option to enforce that the surface plane z=0 remains at temperature T=0

NeumannBoundaryCondition <: BoundaryCondition

Option to enforce zero heat flow at the surface plane z=0.