From 4b06176100779052005d29461c869b3a725b96bd Mon Sep 17 00:00:00 2001 From: Marc Basquens Date: Fri, 19 Jul 2024 12:40:26 +0200 Subject: [PATCH] Add boundary conditions --- Project.toml | 2 +- README.md | 8 ++-- examples/modular/test.jl | 32 +++++++++----- src/BoreholeNetworksSimulator.jl | 1 + .../DirichletBoundaryCondition.jl | 1 + src/modular/boundary_conditions/NoBoundary.jl | 1 + src/modular/core/core.jl | 2 +- src/modular/interfaces/BoundaryCondition.jl | 1 + src/modular/methods/ConvolutionMethod.jl | 6 +-- src/modular/methods/NonHistoryMethod.jl | 43 ++++++++++++++----- src/modular/methods/ground_response.jl | 17 ++++++-- src/modular/methods/ground_water_response.jl | 2 +- src/modular/simulate.jl | 4 +- 13 files changed, 83 insertions(+), 37 deletions(-) create mode 100644 src/modular/boundary_conditions/DirichletBoundaryCondition.jl create mode 100644 src/modular/boundary_conditions/NoBoundary.jl create mode 100644 src/modular/interfaces/BoundaryCondition.jl diff --git a/Project.toml b/Project.toml index 8b26b65..f9ec9ba 100644 --- a/Project.toml +++ b/Project.toml @@ -32,4 +32,4 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "Aqua"] \ No newline at end of file +test = ["Test", "Aqua"] diff --git a/README.md b/README.md index a670456..76539af 100644 --- a/README.md +++ b/README.md @@ -15,7 +15,7 @@ Simply clone the project using e.g. git clone, cd to the project directory, ente ## Example The example considers a borehole field of 48 boreholes connected according to scheme utilized in the installation in Braedsturp, Denmark. - + ## Run the example ``` @@ -45,13 +45,13 @@ julia > include("examples/plots/sim1.jl") Inlet borehole temperatures and heat flows for boreholes along two branches in the borehole field. The time series are color coded according to the previous configuration plot above. In addition to the inlet temperature, the output temperature from the branch (grey dot), and the mean output temperature from the field (black dot) are displayed. - - + + Finally we can display the heatmap of the temperature field in the borehole region during the 10th year of operation - + ## Running the code in Python diff --git a/examples/modular/test.jl b/examples/modular/test.jl index 60dfed3..8a70b99 100644 --- a/examples/modular/test.jl +++ b/examples/modular/test.jl @@ -2,27 +2,29 @@ using BoreholeNetworksSimulator using FiniteLineSource Δt = 8760*3600/12. -Nt = 10 -network = [ +Nt = 1 +configurations = [ BoreholeNetwork([[1]]) ] -Nbr = length(network) α = 1e-6 λ = 3. rb = 0.1 +D = 0. +H = 1. + positions = [(0., i) for i in 0:0] T0 = 0. -borehole = SingleUPipeBorehole(H=1., D=0., rb=rb) +borehole = SingleUPipeBorehole(H=H, D=D, rb=rb) medium = GroundMedium(α=α, λ=λ) borefield = EqualBoreholesBorefield(borehole_prototype=borehole, positions=positions, medium=medium, T0 = T0) -constraint = InletTempConstraint(10*ones(1)) -#constraint = HeatLoadConstraint([1.]) +#constraint = InletTempConstraint(10*ones(1)) +constraint = HeatLoadConstraint([1.]) fluid = Fluid(cpf = 4182., name = "INCOMP::MEA-20%") function operator(i, Tin, Tout, Tb, q) - BoreholeOperation(network[1], 0.1 .* ones(1)) + BoreholeOperation(configurations[1], 0.1 .* ones(1)) end @@ -31,6 +33,7 @@ options_c = SimulationOptions( constraint = constraint, borefield = borefield, fluid = fluid, + #boundary_condition=NoBoundary(), Δt = Δt, Nt = Nt ) @@ -45,6 +48,7 @@ options_nh = SimulationOptions( constraint = constraint, borefield = borefield, fluid = fluid, + #boundary_condition=NoBoundary(), Δt = Δt, Nt = Nt ) @@ -57,12 +61,18 @@ X1-X2 ### FLS -q = [1 for i=1:10] +q = [1 for i=1:Nt] I = zeros(length(q)) - -setup = FiniteLineSource.SegmentToSegment(D1=0., H1=1., D2=0., H2=1., σ=.1) params = FiniteLineSource.Constants(Δt=Δt, rb=rb, b=2.) + +setup = FiniteLineSource.SegmentToSegment(D1=D, H1=H, D2=D, H2=H, σ=rb) precomp = FiniteLineSource.precompute_parameters(setup, params=params) compute_integral_throught_history!(setup, I=I, q=q, precomp=precomp, params=params) -BoreholeNetworksSimulator.sts(setup, params, tstep) +I2 = zeros(length(q)) +setup2 = FiniteLineSource.SegmentToSegment(D1=-D, H1=-H, D2=D, H2=H, σ=rb) +precomp2 = FiniteLineSource.precompute_parameters(setup2, params=params) +compute_integral_throught_history!(setup2, I=I2, q=q, precomp=precomp2, params=params) + +I-I2 + diff --git a/src/BoreholeNetworksSimulator.jl b/src/BoreholeNetworksSimulator.jl index e6e6162..cee6b92 100644 --- a/src/BoreholeNetworksSimulator.jl +++ b/src/BoreholeNetworksSimulator.jl @@ -30,6 +30,7 @@ export Borefield, EqualBoreholesBorefield export Medium, GroundMedium, FlowInPorousMedium export Constraint, HeatLoadConstraint, InletTempConstraint export TimeSuperpositionMethod, ConvolutionMethod, NonHistoryMethod +export BoundaryCondition, NoBoundary, DirichletBoundaryCondition export SimulationOptions, SimulationContainers, BoreholeOperation, BoreholeNetwork, Fluid export simulate, compute_parameters, load_cache!, save_cache, initialize diff --git a/src/modular/boundary_conditions/DirichletBoundaryCondition.jl b/src/modular/boundary_conditions/DirichletBoundaryCondition.jl new file mode 100644 index 0000000..2d27d72 --- /dev/null +++ b/src/modular/boundary_conditions/DirichletBoundaryCondition.jl @@ -0,0 +1 @@ +struct DirichletBoundaryCondition <: BoundaryCondition end \ No newline at end of file diff --git a/src/modular/boundary_conditions/NoBoundary.jl b/src/modular/boundary_conditions/NoBoundary.jl new file mode 100644 index 0000000..b9928a6 --- /dev/null +++ b/src/modular/boundary_conditions/NoBoundary.jl @@ -0,0 +1 @@ +struct NoBoundary <: BoundaryCondition end \ No newline at end of file diff --git a/src/modular/core/core.jl b/src/modular/core/core.jl index ea42d30..4839e2b 100644 --- a/src/modular/core/core.jl +++ b/src/modular/core/core.jl @@ -2,7 +2,6 @@ @with_kw struct BoreholeOperation network mass_flows:: Vector # Mass flow for each branch in network - #cpf::T = 4182. # specific heat capacity end BoreholeOperation(::Nothing) = BoreholeOperation([[]], [0.]) @@ -23,6 +22,7 @@ first_boreholes(network::BoreholeNetwork) = map(first, network.branches) constraint::Constraint borefield::Borefield fluid::Fluid + boundary_condition::BoundaryCondition = DirichletBoundaryCondition() Δt Nt Nb = n_boreholes(borefield) diff --git a/src/modular/interfaces/BoundaryCondition.jl b/src/modular/interfaces/BoundaryCondition.jl new file mode 100644 index 0000000..edc8cd9 --- /dev/null +++ b/src/modular/interfaces/BoundaryCondition.jl @@ -0,0 +1 @@ +abstract type BoundaryCondition end diff --git a/src/modular/methods/ConvolutionMethod.jl b/src/modular/methods/ConvolutionMethod.jl index ee6005c..ce8e81d 100644 --- a/src/modular/methods/ConvolutionMethod.jl +++ b/src/modular/methods/ConvolutionMethod.jl @@ -6,10 +6,10 @@ end ConvolutionMethod() = ConvolutionMethod(zeros(0,0,0), zeros(0,0)) function precompute_auxiliaries!(model::ConvolutionMethod; options::SimulationOptions) - @unpack Nb, Nt, t, borefield = options + @unpack Nb, Nt, t, borefield, boundary_condition = options model.g = zeros(Nb, Nb, Nt) model.q = zeros(Nb, Nt) - compute_response!(model.g, borefield.medium, borefield, t) + compute_response!(model.g, borefield.medium, borefield, boundary_condition, t) return model end @@ -18,7 +18,7 @@ function update_auxiliaries!(method::ConvolutionMethod, X, borefield::Borefield, method.q[:, step] = @view X[3Nb+1:end, step] end -function method_coeffs!(M, method::ConvolutionMethod, borefield::Borefield) +function method_coeffs!(M, method::ConvolutionMethod, borefield::Borefield, ::BoundaryCondition) Nb = n_boreholes(borefield) Ns = n_segments(borefield) M[1:Ns, 3Nb+1:3Nb+Ns] = @view method.g[:,:,1] diff --git a/src/modular/methods/NonHistoryMethod.jl b/src/modular/methods/NonHistoryMethod.jl index e26aedd..e499be1 100644 --- a/src/modular/methods/NonHistoryMethod.jl +++ b/src/modular/methods/NonHistoryMethod.jl @@ -1,3 +1,5 @@ +using .FiniteLineSource: SegmentToSegment, Constants, adaptive_gk_segments, discretization_parameters, f_guess, precompute_coefficients + mutable struct NonHistoryMethod{T} <: TimeSuperpositionMethod F::Matrix{T} ζ::Vector{T} @@ -6,10 +8,11 @@ mutable struct NonHistoryMethod{T} <: TimeSuperpositionMethod end NonHistoryMethod() = NonHistoryMethod(zeros(0, 0), zeros(0), zeros(0, 0), zeros(0)) -FiniteLineSource.SegmentToSegment(s::MeanSegToSegEvParams) = FiniteLineSource.SegmentToSegment(D1=s.D1, H1=s.H1, D2=s.D2, H2=s.H2, σ=s.σ) +SegmentToSegment(s::MeanSegToSegEvParams) = SegmentToSegment(D1=s.D1, H1=s.H1, D2=s.D2, H2=s.H2, σ=s.σ) +image(s::SegmentToSegment) = SegmentToSegment(D1=-s.D1, H1=-s.H1, D2=s.D2, H2=s.H2, σ=s.σ) function precompute_auxiliaries!(model::NonHistoryMethod; options::SimulationOptions) - @unpack Nb, Nt, Ns, Δt, borefield = options + @unpack Nb, Nt, Ns, Δt, borefield, boundary_condition = options b = 2. α = get_α(borefield.medium) rb = get_rb(borefield, 1) @@ -18,9 +21,9 @@ function precompute_auxiliaries!(model::NonHistoryMethod; options::SimulationOpt n_segment = 20 - constants = FiniteLineSource.Constants(Δt=Δt, α=α, rb=rb, kg=kg, b=b) - segments = FiniteLineSource.adaptive_gk_segments(FiniteLineSource.f_guess(FiniteLineSource.SegmentToSegment(get_sts(borefield, 1, 1)), constants), 0., b) - dps = @views [FiniteLineSource.discretization_parameters(s.a, s.b, n_segment) for s in segments] + constants = Constants(Δt=Δt, α=α, rb=rb, kg=kg, b=b) + segments = adaptive_gk_segments(f_guess(SegmentToSegment(get_sts(borefield, 1, 1)), constants), 0., b) + dps = @views [discretization_parameters(s.a, s.b, n_segment) for s in segments] ζ = reduce(vcat, (dp.x for dp in dps)) expΔt = @. exp(-ζ^2 * Δt̃) @@ -30,8 +33,8 @@ function precompute_auxiliaries!(model::NonHistoryMethod; options::SimulationOpt for i in 1:Ns for j in 1:Ns rb = get_rb(borefield, (i-1)*Ns+j) # Get the right rb - sts = get_sts(borefield, i, j) - w[:, (i-1)*Ns+j] = reduce(vcat, [FiniteLineSource.precompute_coefficients(FiniteLineSource.SegmentToSegment(sts), params=constants, dp=dp) for dp in dps]) + sts = SegmentToSegment(get_sts(borefield, i, j)) + w[:, (i-1)*Ns+j] = reduce(vcat, [coefficients_sts(boundary_condition, sts, constants, dp) for dp in dps]) end end @@ -43,6 +46,17 @@ function precompute_auxiliaries!(model::NonHistoryMethod; options::SimulationOpt model.expΔt = expΔt[perm] end +function coefficients_sts(::NoBoundary, sts::SegmentToSegment, params::Constants, dp) + precompute_coefficients(sts, params=params, dp=dp) +end + +function coefficients_sts(::DirichletBoundaryCondition, sts::SegmentToSegment, params::Constants, dp) + image_sts = image(sts) + w1 = precompute_coefficients(sts, params=params, dp=dp) + w2 = precompute_coefficients(image_sts, params=params, dp=dp) + w1 - w2 +end + function get_sts(borefield::Borefield, i, j) xi, yi, Di, Hi = segment_coordinates(borefield, i) xj, yj, Dj, Hj = segment_coordinates(borefield, j) @@ -59,16 +73,15 @@ function update_auxiliaries!(method::NonHistoryMethod, X, borefield::Borefield, end end -function method_coeffs!(M, method::NonHistoryMethod, borefield::Borefield) +function method_coeffs!(M, method::NonHistoryMethod, borefield::Borefield, boundary_condition::BoundaryCondition) Nb = n_boreholes(borefield) Ns = n_segments(borefield) λ = get_λ(borefield.medium) for i in 1:Ns for j in 1:Ns - sts = get_sts(borefield, i, j) - rb = get_rb(borefield, i) - M[i, 3Nb+j] = - q_coef(borefield.medium, method, sts, λ, (i-1)*Ns+j) + sts = SegmentToSegment(get_sts(borefield, i, j)) + M[i, 3Nb+j] = - q_coef(boundary_condition, borefield.medium, method, sts, λ, (i-1)*Ns+j) end end @@ -88,6 +101,14 @@ function method_b!(b, method::NonHistoryMethod, borefield::Borefield, step) end end +function q_coef(::NoBoundary, medium, method, sts, λ, i) + q_coef(medium, method, sts, λ, i) +end + +function q_coef(::DirichletBoundaryCondition, medium, method, sts, λ, i) + q_coef(medium, method, sts, λ, i) - q_coef(medium, method, image(sts), λ, i) +end + function q_coef(::GroundMedium, method, sts, λ, i) @unpack D1, H1, D2, H2, σ = sts @unpack expΔt, w, ζ = method diff --git a/src/modular/methods/ground_response.jl b/src/modular/methods/ground_response.jl index f729fa6..e19aff3 100644 --- a/src/modular/methods/ground_response.jl +++ b/src/modular/methods/ground_response.jl @@ -1,5 +1,5 @@ -function compute_response!(g, medium::GroundMedium, borefield::Borefield, t) +function compute_response!(g, medium::GroundMedium, borefield::Borefield, boundary_condition::BoundaryCondition, t) @unpack λ, α = medium Ns = n_segments(borefield) @@ -13,7 +13,7 @@ function compute_response!(g, medium::GroundMedium, borefield::Borefield, t) rb = get_rb(borefield, i) σ = i == j ? rb : sqrt((x2-x1)^2 + (y2-y1)^2) s = SegmentToSegment(D1=D1, H1=H1, D2=D2, H2=H2, σ=σ) - g[i, j, k] = sts(s, Constants(α=α, kg=λ, rb=rb), tt) + g[i, j, k] = sts(boundary_condition, s, Constants(α=α, kg=λ, rb=rb), tt) end end end @@ -41,7 +41,7 @@ function fls(t, x, y, z, H, D, α, kg, atol = 1e-8) return fls_step_response(t, σ, z, D, D+H, α, kg, atol=atol) end -function sts(s:: SegmentToSegment, params::Constants, t) +function sts_response(s::SegmentToSegment, params::Constants, t) @unpack α, kg = params params = FiniteLineSource.MeanSegToSegEvParams(s) h_mean_sts, r_min, r_max = FiniteLineSource.mean_sts_evaluation(params) @@ -49,3 +49,14 @@ function sts(s:: SegmentToSegment, params::Constants, t) x, w = FiniteLineSource.adaptive_gk(f, r_min, r_max) dot(f.(x), w) end + +function sts(::NoBoundary, s::SegmentToSegment, params::Constants, t) + sts_response(s, params, t) +end + +function sts(::DirichletBoundaryCondition, s::SegmentToSegment, params::Constants, t) + s_image = SegmentToSegment(D1=-s.D1, H1=-s.H1, D2=s.D2, H2=s.H2, σ=s.σ) + Ip = sts_response(s, params, t) + In = sts_response(s_image, params, t) + Ip - In +end diff --git a/src/modular/methods/ground_water_response.jl b/src/modular/methods/ground_water_response.jl index 76e14a1..c1f5251 100644 --- a/src/modular/methods/ground_water_response.jl +++ b/src/modular/methods/ground_water_response.jl @@ -1,5 +1,5 @@ -function compute_response!(g, medium::FlowInPorousMedium, borefield::Borefield, t) +function compute_response!(g, medium::FlowInPorousMedium, borefield::Borefield, boundary_condition::BoundaryCondition, t) @unpack θ, λ, vt, α = medium Ns = segment_amount(borefield) diff --git a/src/modular/simulate.jl b/src/modular/simulate.jl index 8281ae4..4c9c5a8 100644 --- a/src/modular/simulate.jl +++ b/src/modular/simulate.jl @@ -11,7 +11,7 @@ # (4) Heat balance equations function simulate(;operator, options::SimulationOptions, containers::SimulationContainers) - @unpack method, constraint, borefield, fluid = options + @unpack method, constraint, borefield, fluid, boundary_condition = options @unpack Nb, Ns, Nt, Ts = options @unpack M, b, X = containers @@ -45,7 +45,7 @@ function simulate(;operator, options::SimulationOptions, containers::SimulationC end @views constraints_coeffs!(M[constraints_eqs, :], constraint, operation) if i == Ts - @views method_coeffs!(M[method_eqs, :], method, borefield) + @views method_coeffs!(M[method_eqs, :], method, borefield, boundary_condition) end if last_operation.mass_flows != operation.mass_flows @views heat_balance_coeffs!(M[balance_eqs, :], borefield, operation, fluid)