Basic tutorial

In this tutorial we will learn how to simulate the expected temperature and heat extraction in a borehole field.

First, we need to specify the parameters and constraints of our system, as well as some options for the simulation. This is done through a SimulationOptions object. Its variables are modular components with which we can describe many scenarios by leveraging Julia's multiple dispatch. Let us see the components one by one with an example.

We start by specifying the simulation time step and the simulation duration. For our example, we will take monthly time steps during 10 years:

using BoreholeNetworksSimulator

Δt = 8760*3600/12.
Nt = 10*12
120

Suppose the ground we are interested in simulating is made of solid rock. This means the heat transfer will occur by pure conduction. Assume that the thermal diffusivity of the rock is $α = 10^{-6} \frac{m^2}{s}$ and its thermal conductivity is $λ = 3 \frac{W}{m \ K}$. The undisturbed temperature of the ground is $T_0=10 \ ^{\circ}C$. We model the ground with a subtype of Medium, in our case, as per our assumptions, we are particularly interested in GroundMedium:

α = 1e-6
λ = 3.
T0 = 10.
medium = GroundMedium(α=α, λ=λ, T0=T0)


D = 10.
H = 100.
borehole = SingleUPipeBorehole(H=H, D=D)
SingleUPipeBorehole{Float64}
  λg: Float64 2.5
  Cg: Float64 3.1e6
  αg: Float64 8.064516129032258e-7
  rp: Float64 0.02
  λp: Float64 0.42
  dpw: Float64 0.0023
  rpo: Float64 0.0223
  hp: Float64 725.0
  pipe_position: Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}
  rb: Float64 0.0575
  H: Float64 100.0
  D: Float64 10.0
  n_segments: Int64 1
  R_cache: StaticArraysCore.SMatrix{2, 2, Float64, 4}
  A: StaticArraysCore.MMatrix{2, 2, Float64, 4}
  method: ExponentialUtilities.ExpMethodHigham2005
  exp_cache: Tuple{Vector{StaticArraysCore.MMatrix{2, 2, Float64, 4}}, Vector{Float64}}

Next, we need to know how and where the boreholes in our system are. Suppose we are interested in simulating the evolution of two identical vertical boreholes of length $H=100m$ and buried at depth $D=10m$, separated by a distance of $σ=5m$. The boreholes are of the U-type. We need to enconde this information in a subtype of Borefield. Since the boreholes in our example have the same properties, we can use EqualBoreholesBorefield, which instantiates several identical boreholes from a prototype. The prototype is specified by a subtype of Borehole. In our case, we can use SingleUPipeBorehole to model a borehole with a single U-pipe.

D = 10.
H = 100.

borehole = SingleUPipeBorehole(H=H, D=D)
SingleUPipeBorehole{Float64}
  λg: Float64 2.5
  Cg: Float64 3.1e6
  αg: Float64 8.064516129032258e-7
  rp: Float64 0.02
  λp: Float64 0.42
  dpw: Float64 0.0023
  rpo: Float64 0.0223
  hp: Float64 725.0
  pipe_position: Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}
  rb: Float64 0.0575
  H: Float64 100.0
  D: Float64 10.0
  n_segments: Int64 1
  R_cache: StaticArraysCore.SMatrix{2, 2, Float64, 4}
  A: StaticArraysCore.MMatrix{2, 2, Float64, 4}
  method: ExponentialUtilities.ExpMethodHigham2005
  exp_cache: Tuple{Vector{StaticArraysCore.MMatrix{2, 2, Float64, 4}}, Vector{Float64}}

There are more parameters that are also relevant, such as borehole radius, grout properties, pipe resistance, etc., but for the moment we will use their default values.

Next, we need to specify where the borehole are located.

σ = 5.
positions = [(0., 0.), (0., σ)]
borefield = EqualBoreholesBorefield(borehole_prototype=borehole, positions=positions)
EqualBoreholesBorefield{SingleUPipeBorehole{Float64}, Float64}
  borehole_prototype: SingleUPipeBorehole{Float64}
  positions: Array{Tuple{Float64, Float64}}((2,))

Note that we didn't specify yet how the boreholes are connected. This is because the simulation allows for dynamical changes in the network hydraulic configuration. Using different configurations, we can simulate the reverse flow, or even simulate valves opening and closing. The reason behind this design is to allow for decision making on the operation of the borefield during the simulation, based on the inputs that we are interested in.

For now, let us state the possible hydraulic configurations that we will allows our borefield to be in. This is done with an vector of BoreholeNetwork. Borefields usually have several branches of boreholes. Each borehole in a branch is connected in series, while branches may o may not be connected in parallel between themselves. To create a BoreholeNetwork, we need to provide all of its branches in a vector. Each branch is, in turn, a vector containing the identifiers, in order, of the boreholes present in that branch. Each identifier i::Int refers to the borehole located at position positions[i].

In our example, we want to simulate two independent boreholes, so each of them must be in a separate branch. Also, for the moment, we are only interested in this configuration, so we define:

configurations = [BoreholeNetwork([[1], [2]])]
1-element Vector{BoreholeNetwork}:
 BoreholeNetwork([[1], [2]], 2)

Even with all these specifications, the evolution of the system is still not fully determined. The missing conditions are referred to as constraints, and are modeled by subtypes of Constraint. For instance, if we would like the two boreholes to be connected in parallel, we would still need to impose that their inlet temperatures be equal. In our example, since we want out boreholes to be independent, we will impose the total amount of heat that we want to extract from each borehole. We will impose a constant load, equal for both boreholes. This is specified by

q1 = 5.
q2 = 5.
constraint = constant_HeatLoadConstraint([q1, q2], Nt)
HeatLoadConstraint{Float64}([5.0 5.0 … 5.0 5.0; 5.0 5.0 … 5.0 5.0])

We can finally create the object with all the options:

options = SimulationOptions(
    method = ConvolutionMethod(),
    constraint = constraint,
    borefield = borefield,
    medium = medium,
    fluid = Water(),
    Δt = Δt,
    Nt = Nt,
    configurations = configurations
)
SimulationOptions{ConvolutionMethod{Float64}, HeatLoadConstraint{Float64}, EqualBoreholesBorefield{SingleUPipeBorehole{Float64}, Float64}, GroundMedium{Float64}, DirichletBoundaryCondition, Water}
  method: ConvolutionMethod{Float64}
  constraint: HeatLoadConstraint{Float64}
  borefield: EqualBoreholesBorefield{SingleUPipeBorehole{Float64}, Float64}
  medium: GroundMedium{Float64}
  fluid: Water
  boundary_condition: DirichletBoundaryCondition DirichletBoundaryCondition()
  Δt: Float64 2.628e6
  Nt: Int64 120
  Nb: Int64 2
  Ns: Int64 2
  Ts: Int64 1
  Tmax: Float64 3.1536e8
  t: StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}
  configurations: Array{BoreholeNetwork}((1,))

As we have mentioned, the simulation is designed to allow for a controllable opeartion during its duration. We do this by creating a subtype of Operator containing the operation strategy. Then, we must implement a method of the function operate with our Operator subtype, that takes as an input the current state of the borefield and outputs a BoreholeOperation object. BoreholeOperation has two variables: the first specifies which configuration will be used for the next time step. In our case, we only want a static, simple configuration. The second specifies the mass flow rate through each branch of the selected configuration, provided as a vector. In our example, we will keep this constant through the simulation. For this purpose, there exists the type SimpleOperator, that implements precisely this strategy.

operator = SimpleOperator(mass_flow = 2., branches = 2)
SimpleOperator{Float64}([2.0, 2.0])

Before simulating, we first need to call initialize to run some precomputations that will be used throught the simulation and to instantiate containers where the result will be written.

containers = initialize(options)
SimulationContainers{Float64, Matrix{Float64}}
  M: Array{Float64}((8, 8)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  X: Array{Float64}((8, 120)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  b: Array{Float64}((8,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  mf: Array{Float64}((2, 120)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

And finally, we can start the simulation.

@time simulate!(operator=operator, options=options, containers=containers)
  0.000778 seconds (607 allocations: 117.375 KiB)

The simulation is over! Note that the bulk of the time is spent in the precoputation, while the simulation itself is quite fast. This is especially good if we want to test different operation strategies.

The result is saved in

containers.X
8×120 Matrix{Float64}:
 10.0128  10.0139  10.0146  10.0152  …  10.0226  10.0226  10.0226  10.0226
 10.0122  10.0133  10.014   10.0146     10.022   10.022   10.022   10.022
 10.0128  10.0139  10.0146  10.0152     10.0226  10.0226  10.0226  10.0226
 10.0122  10.0133  10.014   10.0146     10.022   10.022   10.022   10.022
 10.0099  10.011   10.0117  10.0123     10.0197  10.0197  10.0197  10.0198
 10.0099  10.011   10.0117  10.0123  …  10.0197  10.0197  10.0197  10.0198
  0.05     0.05     0.05     0.05        0.05     0.05     0.05     0.05
  0.05     0.05     0.05     0.05        0.05     0.05     0.05     0.05

Each column in containers.X contains the results for a time step. In each column, the first $2 N_b$ correspond to the fluid temperatures. In particular, the odd entries containers.X[1:2:2Nb, :] are the inlet temperatures, while the even entries containers.X[2:2:2Nb, :] are the outlet temperatures of each borehole. The next $Nb$ values containers.X[2Nb+1:3Nb, :] are the borehole wall temperatures. The last $Nb$ values containers.X[3Nb+1:4Nb, :] are the heat extraction of each borehole.


This page was generated using Literate.jl.