diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index c1e023a..4e3ba82 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -5,7 +5,6 @@ on: - main tags: ['*'] pull_request: - workflow_dispatch: concurrency: # Skip intermediate builds: always. # Cancel intermediate builds: only if it is a pull request build. @@ -47,14 +46,6 @@ jobs: permissions: contents: write steps: - - name: Increase swapfile - run: | - sudo swapoff -a - sudo fallocate -l 15G /swapfile - sudo chmod 600 /swapfile - sudo mkswap /swapfile - sudo swapon /swapfile - sudo swapon --show - name: Add GLMakie/XFVB dependencies run: sudo apt-get update && sudo apt-get install -y xorg-dev mesa-utils xvfb libgl1 freeglut3-dev libxrandr-dev libxinerama-dev libxcursor-dev libxi-dev libxext-dev libcairo2-dev libfreetype6-dev libffi-dev libjpeg-dev libpng-dev libz-dev - name: Checkout @@ -80,4 +71,5 @@ jobs: using Documenter: DocMeta, doctest using Fimbul DocMeta.setdocmeta!(Fimbul, :DocTestSetup, :(using Fimbul); recursive=true) - doctest(Fimbul)' \ No newline at end of file + doctest(Fimbul)' + \ No newline at end of file diff --git a/docs/src/refs.bib b/docs/src/refs.bib index e76867b..87e8a5d 100644 --- a/docs/src/refs.bib +++ b/docs/src/refs.bib @@ -9,4 +9,29 @@ @article{egg_model year = {2014}, month = nov, pages = {192–195} -} \ No newline at end of file +} + +@article{eskilson_1988, + title = {Simulation model for thermally interacting heat extraction boreholes}, + volume = {13}, + ISSN = {0149-5720}, + url = {http://dx.doi.org/10.1080/10407788808913609}, + DOI = {10.1080/10407788808913609}, + number = {2}, + journal = {Numerical Heat Transfer}, + publisher = {Informa UK Limited}, + author = {Eskilson, Per and Claesson, Johan}, + year = {1988}, + month = mar, + pages = {149–165} +} + +@techreport{feflow_btes, + author = {Hans-J{\"u}rgen G. Diersch and Daniel Bauer and Wolfgang Heidemann and Wolfgang R{\"u}haak and Peter Sch{\"a}tzl}, + title = {Finite-Element Modeling of Borehole Heat Exchangers in Geothermal Heating Systems}, + institution = {DHI-WASY GmbH and Institute of Thermodynamics and Thermal Engineering (ITW), University of Stuttgart}, + type = {FEFLOW White Papers, Vol. 5}, + year = {2011}, + address = {Berlin, Germany}, + keywords = {FEFLOW, Borehole Heat Exchanger, Geothermal, Numerical Simulation}, +} diff --git a/examples/analytical/analytical_closed_loop.jl b/examples/analytical/analytical_closed_loop.jl new file mode 100644 index 0000000..fad3be7 --- /dev/null +++ b/examples/analytical/analytical_closed_loop.jl @@ -0,0 +1,296 @@ +# # Analytical solutions for closed-loop geothermal wells +# This example validates the numerical simulation of a closed-loop geothermal +# well system against analytical solutions for both U-tube and coaxial +# closed-loop configurations assuming fixed ground temperature. The validation +# demonstrates the accuracy of Fimbul's closed-loop well model implementation. + +# Add modules to namespace +using Jutul, JutulDarcy, Fimbul +using HYPRE +using GLMakie + +# Useful SI units +meter = si_unit(:meter) +kilogram = si_unit(:kilogram) +atm = si_unit(:atm) +second, day = si_units(:second, :day) +Kelvin, Joule, Watt = si_units(:kelvin, :joule, :watt) +darcy = si_unit(:darcy); + +# ## Problem parameters +# We set up common parameters for both BTES configurations. +# These represent typical values for shallow geothermal systems. + +# ### Operational conditions +# We inject fluid at 20.0 m^3/day and 80 °C temperature into the well, and +# assume a constant ground temperature of 10 °C. +T_in = convert_to_si(80.0, :Celsius) # Inlet temperature (supply fluid) +T_rock = convert_to_si(10.0, :Celsius) # Ground temperature (constant far-field) +Q = 20.0*meter^3/day # Volumetric flow rate through well + +# ### Fluid properties (water-based heat transfer fluid) +# To enable comparison with analytical solutions, we use constant fluid +# properties. +ρf = 988.1*kilogram/meter^3 # Fluid density at operating temperature +Cpf = 4184.0*Joule/(kilogram*Kelvin) # Specific heat capacity of fluid +λf = 0.6405*Watt/(meter*Kelvin) # Thermal conductivity of fluid + +# ### Rock formation properties (typical for sedimentary rock) +ϕ = 0.01 # Porosity (very low for compact rock) +K = 1e-3*darcy # Permeability (practically impermeable) +ρr = 2650.0*kilogram/meter^3 # Rock density +Cpr = 900.0*Joule/(kilogram*Kelvin) # Rock specific heat capacity +λr = 2.5*Watt/(meter*Kelvin) # Rock thermal conductivity + +# ### Closed-loop installation materials +ρg = 2000.0*kilogram/meter^3 # Grout density (cement/sand mixture) +Cpg = 1500.0*Joule/(kilogram*Kelvin) # Grout specific heat capacity +λg = 2.3*Watt/(meter*Kelvin) # Grout thermal conductivity +λp = 0.38*Watt/(meter*Kelvin); # Pipe wall thermal conductivity (HDPE) + +# ### Closed-loop configurations +# We will consider two different closed-loop configurations: U-tube and coaxial. +# +# **U-tube configuration**: The fluid flows down one pipe, makes a U-turn at +# the bottom, and returns to the surface through a parallel pipe. Both pipes +# are surrounded by grout within a single borehole. +# +# **Coaxial configuration**: The fluid flows down an outer pipe and returns +# through an inner pipe concentrically located inside the outer pipe, or vice +# versa. This design allows for more compact installations. +# +# In both configurations, the pipes are surrounded by grout material, which +# provides thermal contact with the rock formation while maintaining structural +# integrity of the borehole. + +L = 100.0*meter # Closed-loop well length (vertical depth) +cfg_u1 = ( # U-tube parameters + closed_loop_type = :u1, + radius_grout = 65e-3*meter, # Borehole radius (grout outer boundary) + radius_pipe = 16e-3*meter, # Pipe outer radius (both supply/return) + wall_thickness_pipe = 2.9e-3*meter, # Pipe wall thickness + pipe_spacing = 60e-3*meter, # Center-to-center distance between pipes + pipe_thermal_conductivity = λp, # Thermal conductivity of pipe material +) + +cfg_coax = ( # Coaxial parameters + closed_loop_type = :coaxial, + radius_grout = 50e-3*meter, # Borehole radius (grout outer boundary) + radius_pipe_inner = 12e-3*meter, # Inner pipe outer radius + wall_thickness_pipe_inner = 3e-3*meter, # Inner pipe wall thickness + radius_pipe_outer = 25e-3*meter, # Outer pipe outer radius + wall_thickness_pipe_outer = 4e-3*meter, # Outer pipe wall thickness + inner_pipe_thermal_conductivity = λp, # Inner pipe material conductivity + outer_pipe_thermal_conductivity = λp # Outer pipe material conductivity +); + +# ## Utility functions +# We set up convenience functions to create closed-loop simulation cases and +# visualize results against analytical solutions. To mimic constant ground +# temperature, we set up a very large reservoir domain with a single cell in the +# horizontal directions. + +function setup_closed_loop_single( # Utility function to set up closed-loop simulation case + type; # Type of closed-loop geometry (:u1 or :coaxial) + nz = 100, # Number of cells in vertical direction + n_step = 1, # Number of time steps + inlet = :outer # Inlet pipe for coaxial geometry (:outer or :inner + ) + + ## Select well configuration based on geometry type + if type == :u1 + well_args = cfg_u1 + elseif type == :coaxial + well_args = cfg_coax + else + error("Unknown closed loop type: $type") + end + ## Create computational domain + dims = (1, 1, nz) + Δ = (100000.0*dims[3]/L, 100000.0*dims[3]/L, L) + msh = CartesianMesh(dims, Δ) + reservoir = reservoir_domain(msh; + rock_density = ρr, + rock_heat_capacity = Cpr, + rock_thermal_conductivity = λr, + fluid_thermal_conductivity = λf, + component_heat_capacity = Cpf, + porosity = ϕ, + permeability = K + ) + ## Set up closed-loop well + wells = Fimbul.setup_vertical_btes_well(reservoir, 1, 1; + name = :CL, + well_args..., + grouting_density = ρg, + grouting_heat_capacity = Cpg, + grouting_thermal_conductivity = λg + ) + ## Set up reservoir model + sys = SinglePhaseSystem(AqueousPhase(), reference_density = ρf) + model = setup_reservoir_model(reservoir, sys; wells = wells, thermal = true) + ## Set up controls and boundary conditions + ctrl_inj = InjectorControl( # Injection control with fixed temperature + TotalRateTarget(Q), [1.0], + density=ρf, temperature=T_in) + ctrl_prod = ProducerControl( # Production control with fixed BHP + BottomHolePressureTarget(1.0*atm)) + geo = tpfv_geometry(msh) + bc_cells = geo.boundary_neighbors + domain = reservoir_model(model).data_domain + bc = flow_boundary_condition( # Fixed pressure and temperature BCs + bc_cells, domain, 5.0*atm, T_rock) + ## Set up forces and initial state + if type == :coaxial && inlet == :outer # Coax with injection in outer pipe + ctrl_supply = ctrl_prod + ctrl_return = ctrl_inj + else # Coax with injection in inner pipe or U1 + ctrl_supply = ctrl_inj + ctrl_return = ctrl_prod + end + controls = Dict( + :CL_supply => ctrl_supply, + :CL_return => ctrl_return + ) + forces = setup_reservoir_forces(model; control=controls, bc = bc) + state0 = setup_reservoir_state( # Initialize with uniform pressure and temperature + model; Pressure = 5*atm, Temperature = T_rock) + ## Set time large enough to ensure steady-state + rg = well_args.radius_grout + time = 5/4*2*rg*(ϕ*ρf*Cpf + (1 - ϕ)*ρr*Cpr)/(ϕ*λf + (1 - ϕ)*λr)*10 + dt = fill(time/n_step, n_step) + ## Create Jutul case + case = JutulCase(model, dt, forces; state0 = state0) + ## Set up simulator with temperature-based timestep selector + sim, cfg = setup_reservoir_simulator( + case; initial_dt = 1.0); + sel = VariableChangeTimestepSelector(:Temperature, 2.5; + relative = false, model = :CL_supply) + push!(cfg[:timestep_selectors], sel) + + return case, sim, cfg +end + +function plot_closed_loop( # Utility function to plot closed-loop simulation results + case, simulated, analytical; title = "Closed-loop temperature profiles" + ) + + ## Create figure + fig = Figure(size = (800, 600)) + ax = Axis(fig[1, 1], + title = title, + xlabel = "Temperature (°C)", + ylabel = "Depth (m)", + yreversed = true # Surface at top, depth increases downward + ) + ## Extract numerical solution from final state + well = case.model.models[:CL_supply].data_domain + sim_temp = convert_from_si.( + simulated.result.states[end][:CL_supply][:Temperature], :Celsius) + za = collect(range(0, L, step=0.1)) # Analytical depth points + section_handles = [] + solution_handles = [] + section_names = String[] + solution_names = ["Analytical", "Fimbul"] + ## Plot results for each well section + sections = last.(well[:section]) |> unique |> collect + colors = cgrad(:BrBG_4, 4, categorical=true)[[1,2,4,3]] + err_inf = -Inf + for (i, section) in enumerate(sections) + ## Section name for printing + name = titlecase(replace(string(section), "_" => " ")) + ## Plot analytical solution + Ta = convert_from_si.(analytical[section].(za), :Celsius) + la = lines!(ax, Ta, za; color=colors[i], linewidth = 8, linestyle = :dash, label=name) + ## Plot numerical solution + cells = last.(well[:section]) .== section + Tn = sim_temp[cells] + zn = well[:cell_centroids][3, cells] + ln = lines!(ax, Tn, zn; color=colors[i], linewidth = 2) + ## Store handles and names for legend + push!(section_handles, la) + push!(section_names, name) + if i == 1 + push!(solution_handles, la, ln) + end + ## Compute maximum error in section + Ta_n = convert_from_si.(analytical[section].(zn), :Celsius) + e_inf = maximum(abs.(Tn .- Ta_n)) + err_inf = max(err_inf, e_inf) + end + ## Add legend + Legend(fig[1, 2], [section_handles, solution_handles], [section_names, solution_names], + ["Sections", "Solutions"]; orientation = :vertical) + + return fig, err_inf + +end; + +# ## Validate U-tube configuration +# We set up and simulate a closed-loop system with U-tube geometry, and compare +# the numerical results against the analytical solution. This validates the +# implementation for the most common closed-loop configuration. +case_u1, sim, cfg = setup_closed_loop_single(:u1; nz=125); +res_u1 = simulate_reservoir(case_u1; simulator = sim, config = cfg) + +# ### Analytical solution +# The function `analytical_closed_loop_u1` computes the analytical steady-state +# solution for a vertical U-tube closed-loop system using the method of +# [eskilson_1988](@cite) as described in [feflow_btes](@cite). The function can +# be called with explicit parameters, but also has a convenience version that +# extracts parameters directly from the well model. +analytical_u1 = Fimbul.analytical_closed_loop_u1(Q, T_in, T_rock, + ρf, Cpf, case_u1.model.models[:CL_supply].data_domain) + +# ### Plot comparison +# The numerical solution agrees very well with the analytical solution. +fig, err_inf = plot_closed_loop(case_u1, res_u1, analytical_u1; + title = "U-tube temperature profiles") +println("Maximum error in U-tube configuration: $(round(err_inf, digits=3)) °C") +fig + +# ## Validate coaxial configuration with outer pipe inlet +# Next, we set up and simulate a closed-loop system with coaxial geometry, +# injecting fluid through the outer (annular) pipe, and compare the numerical +# results against the analytical solution. +case_coax_outer, sim, cfg = setup_closed_loop_single(:coaxial; nz=125, inlet = :outer) +res_coax_outer = simulate_reservoir(case_coax_outer; simulator = sim, config = cfg, info_level = 0) + +# ### Analytical solution +# Following [feflow_btes](@cite), we can also compute the analytical +# steady-state solution for coaxial closed-loop systems. This is implemented in +# `analytical_closed_loop_coaxial`. As or U-tube closed loops, this can be +# called with all parameters, or extract parameters from the well model. The +# function also accepts an `inlet` keyword argument to specify whether the +# fluid is injected through the outer or inner pipe. +analytical_coax_outer = Fimbul.analytical_closed_loop_coaxial(Q, T_in, T_rock, + ρf, Cpf, case_coax_outer.model.models[:CL_supply].data_domain; inlet = :outer) + +# ### Plot comparison +# Coaxial systems typically show different thermal behavior than U-tubes due to +# the closer thermal coupling between supply and return flows. +fig, err_inf = plot_closed_loop(case_coax_outer, res_coax_outer, analytical_coax_outer; + title = "Coaxial temperature profiles (outer inlet)") +println("Maximum error in coaxial outer inlet configuration: $(round(err_inf, digits=3)) °C") +fig + +# ## Validate coaxial configuration with inner pipe inlet +# Finally, we validate the coaxial configuration with inner pipe inlet. This +# configuration is often preferred when storing heat, as the hot fluid is not +# exposed to the rock before it reaches the bottom of the well, ensuring a +# higher temperature difference along the entire wellbore, and consequently more +# heat storage. +case_coax_inner, sim, cfg = setup_closed_loop_single(:coaxial; nz=125, inlet = :inner) +res_coax_inner = simulate_reservoir(case_coax_inner; simulator = sim, config = cfg, info_level = 0) + +# ### Analytical solution +# We compute the analytical solution for the coaxial configuration with inner +# pipe inlet. +analytical_coax_inner = Fimbul.analytical_closed_loop_coaxial(Q, T_in, T_rock, + ρf, Cpf, case_coax_inner.model.models[:CL_supply].data_domain; inlet = :inner) + +# ### Plot comparison +fig, err_inf = plot_closed_loop(case_coax_inner, res_coax_inner, analytical_coax_inner; + title = "Coaxial temperature profiles (inner inlet)") +println("Maximum error in coaxial inner inlet configuration: $(round(err_inf, digits=3)) °C") +fig \ No newline at end of file diff --git a/examples/production/ags_demo.jl b/examples/production/ags_demo.jl index 6c8cc4c..eabd4a4 100644 --- a/examples/production/ags_demo.jl +++ b/examples/production/ags_demo.jl @@ -182,7 +182,7 @@ hidexdecorations!(ax_tmp, grid = false) ax_pwr = Axis( # Panel 2: Lateral power fig[2, 1:3]; title = "Lateral Power", - ylabel = "Power (W)", xlabel = "Time (years)", + ylabel = "Power (MW)", xlabel = "Time (years)", limits = (nothing, (-0.05, 0.5))) MW = si_unit(:mega)*si_unit(:watt) plot_lateral_data!(ax_pwr, time, section_data[:Power]./MW, stacked = true) diff --git a/examples/production/doublet_demo.jl b/examples/production/doublet_demo.jl index 4da588f..6446220 100644 --- a/examples/production/doublet_demo.jl +++ b/examples/production/doublet_demo.jl @@ -89,15 +89,7 @@ fig # compute the change in reservoir states from the initial state and plot the # results interactively. The change in temperature is particularly interesting # as it shows the evolution of the cold front in the aquifer -Δstates = [] -for state in results.states - Δstate = Dict{Symbol, Any}() - for (k, v) in state - v0 = case.state0[:Reservoir][k] - Δstate[k] = v .- v0 - end - push!(Δstates, Δstate) -end +Δstates = JutulDarcy.delta_state(results.states, case.state0[:Reservoir]) plot_reservoir(case.model, Δstates; colormap = :seaborn_icefire_gradient, key = :Temperature, aspect = :data) diff --git a/examples/production/egs_demo.jl b/examples/production/egs_demo.jl index d1908d5..99f6878 100644 --- a/examples/production/egs_demo.jl +++ b/examples/production/egs_demo.jl @@ -137,17 +137,7 @@ plot_reservoir(case.model, results.states; plot_res_args...) # It is often more informative to visualize the deviation from the inital # conditions to highlight thermal depletion zones. We therefore compute the # change in reservoir variables to the initial state for all timesteps. -Δstates = [] -for state in results.states - Δstate = Dict{Symbol, Any}() - for (k, v) in state - if haskey(case.state0[:Reservoir], k) - v0 = case.state0[:Reservoir][k] - Δstate[k] = v .- v0 - end - end - push!(Δstates, Δstate) -end +Δstates = JutulDarcy.delta_state(results.states, case.state0[:Reservoir]) plot_reservoir(case.model, Δstates; plot_res_args...) # ### Well Performance Analysis diff --git a/examples/storage/ates_demo.jl b/examples/storage/ates_demo.jl index 141b332..37ae519 100644 --- a/examples/storage/ates_demo.jl +++ b/examples/storage/ates_demo.jl @@ -233,17 +233,16 @@ color_stages!(ax_eta) fig_wells # The efficiency improves with each cycle as the thermal plume matures and -# losses decrease. Notably, efficiency exceeds 100% after the second year -# because the hot well produces at higher rates than the cold well injects due -# to lower viscosity of heated water. This rate imbalance extracts more energy -# than stored, causing progressive aquifer cooling as shown below. While -# beneficial for energy recovery in the short term, this trend requires -# monitoring for system sustainability. +# losses decrease, starting at approximately 81% in the first year and reaching +# around 89% by year five. # ### Aquifer temperature evolution # Examine the spatial and temporal temperature distribution in the aquifer to -# understand thermal plume propagation. We analyze temperature profiles along -# a horizontal transect between wells and track aquifer-wide temperature statistics. +# understand thermal plume propagation. We analyze temperature profiles along a +# horizontal transect between wells and track aquifer-wide temperature +# statistics. Notice how the hot and cold plume interactions can be seen as +# asymmetries in the temperature profiles around the cold well during charging +# and hot well during discharging. # Extract temperature along a horizontal line in the aquifer layer # This transect passes through the center of the aquifer between the two wells @@ -322,6 +321,7 @@ lines!(ax, times, T_mean, color = lcolor, linewidth = 3) fig # The proposed setup results in a maximum increase in the mean aquifer -# temperature of approximately 1°C due to charging. However, the setup also -# results in a higher rate during discharging compared to charging, leading to a -# net cooling effect. \ No newline at end of file +# temperature of approximately 3.8°C due to charging. Due to a system efficiency +# of less than 100%, we also see a net heating effect, but this will likely +# stabilize over longer operational periods as thermal losses balance with +# recovery. \ No newline at end of file diff --git a/examples/storage/btes_demo.jl b/examples/storage/btes_demo.jl index 2cc8e8c..4926b17 100644 --- a/examples/storage/btes_demo.jl +++ b/examples/storage/btes_demo.jl @@ -56,8 +56,11 @@ simulator, config = setup_reservoir_simulator(case; presolve_wells = true, tol_cnv = 1e-2, tol_mb = 1e-6, + tol_cnve_well = Inf, + inc_tol_dT = 1e-2, timesteps = :auto, initial_dt = 5.0, + target_its = 10, relaxation = true); # The transition from charging to discharging creates a thermal shock that is diff --git a/src/Fimbul.jl b/src/Fimbul.jl index b2b3a03..9499ac7 100644 --- a/src/Fimbul.jl +++ b/src/Fimbul.jl @@ -35,7 +35,10 @@ module Fimbul include("meshing/fractured.jl") include("meshing/utils.jl") # Wells - include("wells/btes.jl") + include("wells/closed_loop.jl") + include("wells/closed_loop_u1.jl") + include("wells/closed_loop_coaxial.jl") + include("wells/closed_loop_analytical.jl") # Cases include("cases/cases.jl") # Optimization diff --git a/src/cases/ags.jl b/src/cases/ags.jl index b82260a..975a37e 100644 --- a/src/cases/ags.jl +++ b/src/cases/ags.jl @@ -96,7 +96,7 @@ function ags(; rock_heat_capacity = process_prop( :rock_heat_capacity, rock_heat_capacity, layers) - domain = reservoir_domain(msh; + domain = reservoir_domain(msh; porosity = porosity, permeability = permeability, rock_thermal_conductivity = rock_thermal_conductivity, @@ -217,7 +217,7 @@ function setup_ags_wells(domain, well_coords, well_connectivity) for k = 1:length(well_coords) from = well_connectivity[k,1] if from != 0 - ix = findfirst(rcells[k][1] .== rcells[from]) + ix = findfirst(rcells[from] .== rcells[k][1]) if isnothing(ix) error("Could not find connection from well section $from to well $k") end @@ -226,7 +226,7 @@ function setup_ags_wells(domain, well_coords, well_connectivity) end to = well_connectivity[k,2] if to != 0 - ix = findfirst(rcells[k][end] .== rcells[to]) + ix = findfirst(rcells[to] .== rcells[k][end]) if isnothing(ix) error("Could not find connection from well section $k to well $to") end @@ -248,6 +248,7 @@ function setup_ags_wells(domain, well_coords, well_connectivity) common_args = ( simple_well = false, type = :closed_loop, + radius = 150e-3, WI = 0.0 ) @@ -256,12 +257,14 @@ function setup_ags_wells(domain, well_coords, well_connectivity) perforation_cells_well = wcells, dir = directions, well_cell_centers = geo.cell_centroids[:, rcells], + end_nodes = [length(wcells)], name = :AGS_supply, common_args... ) w_return = setup_well(domain, rcells[end]; name = :AGS_return, + WIth = 0.0, common_args... ) diff --git a/src/cases/ates.jl b/src/cases/ates.jl index 289747c..3b0b0f7 100644 --- a/src/cases/ates.jl +++ b/src/cases/ates.jl @@ -47,6 +47,7 @@ function ates(; permeability = [1.0, 5.0, 1000.0, 5.0, 1.0].*1e-3.*si_unit(:darcy), rock_thermal_conductivity = [2.5, 2.0, 1.5, 2.0, 2.5].*watt/(meter*Kelvin), rock_heat_capacity = 900.0*joule/(kilogram*Kelvin), + rock_density = fill(2600.0*kilogram/meter^3, length(depths)-1), aquifer_layer = 3, temperature_charge = convert_to_si(95, :Celsius), temperature_discharge = convert_to_si(25, :Celsius), @@ -71,20 +72,20 @@ function ates(; rock_thermal_conductivity = process_value(rock_thermal_conductivity) rock_heat_capacity = process_value(rock_heat_capacity) # Calculate thermal properties for the aquifer layer - Cf, Cr = 4184.0, rock_heat_capacity[aquifer_layer] + + Cf = 4184.0*joule/(kilogram*Kelvin)*1000.0*kilogram/meter^3 # Fluid heat capacity + Cr = rock_heat_capacity[aquifer_layer]*rock_density[aquifer_layer] ϕ = porosity[aquifer_layer] Caq = Cf*ϕ + Cr*(1 - ϕ) Haq = layer_thickness[aquifer_layer] # Default charging duration (6 months) - # TODO: Use charge duration from input - charge_duration = 0.5si_unit(:year) ch_start = Fimbul.process_time(2025, charge_period[1]) ch_end = Fimbul.process_time(2025, charge_period[2], true) charge_duration = (ch_end - ch_start).value*1e-3 thermal_radius = missing # Calculate rate based on a target 250 m thermal radius if not provided if ismissing(rate_charge) - thermal_radius = ismissing(well_distance) ? 250.0 : well_distance/2 + thermal_radius = ismissing(well_distance) ? 125 : well_distance/4 rate_charge = (thermal_radius^2*Caq*π*Haq/Cf)/charge_duration # Adjust for 2D simulation (per unit width) if use_2d @@ -102,12 +103,12 @@ function ates(; end # Set well distance based on thermal radius if not specified if ismissing(well_distance) - well_distance = 2*thermal_radius + well_distance = 4*thermal_radius end # ## Create reservoir domain # Create mesh - nearwell_radius = 0.5*min(thermal_radius, well_distance/2) + nearwell_radius = min(thermal_radius, well_distance/2) msh, layers = make_ates_cart_mesh(well_distance, depths, aquifer_layer; nearwell_radius = nearwell_radius, use_2d = use_2d, @@ -209,7 +210,7 @@ function ates(; end # Hot well: inject hot water ctrl_hot = InjectorControl( - rate_target, [1.0], density = rho, temperature = temperature_charge) + rate_target, [1.0], density = rho, temperature = temperature_charge; check = false) # Cold well: produce water rate_target = TotalRateTarget(-rate_charge) ctrl_cold = ProducerControl(rate_target) @@ -227,7 +228,7 @@ function ates(; rate_target = TotalRateTarget(rate_discharge) end ctrl_cold = InjectorControl( - rate_target, [1.0], density = rho, temperature = temperature_discharge) + rate_target, [1.0], density = rho, temperature = temperature_discharge; check = false) # Set up forces control = Dict(:Hot => ctrl_hot, :Cold => ctrl_cold) forces_discharge = setup_reservoir_forces(model; bc = bc, control = control) @@ -363,7 +364,7 @@ function make_ates_cart_mesh(well_distance, depths, aquifer_layer; hz_min = missing, hz_max = missing, nearwell_radius = missing, - offset_rel = 1.0, + offset_rel = 3.0, use_2d = false, ) @@ -391,14 +392,12 @@ function make_ates_cart_mesh(well_distance, depths, aquifer_layer; xy_hot[1]-offset, # Left boundary xy_hot[1]-nearwell_radius, # Left thermal boundary xy_hot[1]-hxy_min/2, # Near hot well - xy_hot[1]+nearwell_radius, # Between wells - xy_cold[1]-nearwell_radius,# Between wells xy_cold[1]-hxy_min/2, # Near cold well xy_cold[1]+nearwell_radius,# Right thermal boundary xy_cold[1]+offset # Right boundary ] # Grid spacing array: coarse-fine-fine-coarse-fine-fine-coarse - hx = [hxy_max, hxy_min, hxy_min, hxy_max, hxy_min, hxy_min, hxy_max]; + hx = [hxy_max, hxy_min, hxy_min, hxy_min, hxy_max]; dx = diff(x) hx[dx .< hxy_max] .= hxy_min nx = round.(dx./hx) diff --git a/src/cases/btes.jl b/src/cases/btes.jl index 7f0f514..3dae81f 100644 --- a/src/cases/btes.jl +++ b/src/cases/btes.jl @@ -88,10 +88,10 @@ function btes(; v = (d > 0) ? xw./d : (1.0, 0.0) # Shift coordiates a bit to avoid being exactly on the node xw = xw .+ (hxy_min/2) .* v - trajectory = [xw[1] xw[2] 0.0; xw[1] xw[2] depths[end]] + trajectory = [xw[1] xw[2] 0.0 + 1e-3; xw[1] xw[2] depths[end] - 1e-3] cells = Jutul.find_enclosing_cells(mesh, trajectory, n = 100) filter!(c -> layers[c] ∈ well_layers, cells) - w_sup, w_ret = setup_btes_well(domain, cells, name=name, btes_type=:simple) + w_sup, w_ret = setup_btes_well(domain, cells, name=name, closed_loop_type=:u1) push!(well_models, w_sup, w_ret) end # Make the model diff --git a/src/cases/doublet.jl b/src/cases/doublet.jl index ebed64e..cdaccd2 100644 --- a/src/cases/doublet.jl +++ b/src/cases/doublet.jl @@ -118,9 +118,8 @@ function geothermal_doublet(; rho = reservoir_model(model).system.rho_ref[1] inj_target = JutulDarcy.ReinjectionTarget(NaN, [:Producer]) - # inj_target = TotalRateTarget(rate) ctrl_inj = InjectorControl(inj_target, [1.0], - density=rho, temperature=temperature_inj, tracers = [1.0]) + density=rho, temperature=temperature_inj, tracers = [1.0], check=false) # BHP control for producer prod_target = TotalRateTarget(-rate) ctrl_prod = ProducerControl(prod_target); diff --git a/src/cases/egs.jl b/src/cases/egs.jl index 960baab..6847cfd 100644 --- a/src/cases/egs.jl +++ b/src/cases/egs.jl @@ -141,7 +141,7 @@ function egs(well_coords, fracture_radius, fracture_spacing; rate_target = TotalRateTarget(rate) end ctrl_inj = InjectorControl( - rate_target, [1.0], density = rho, temperature = temperature_inj) + rate_target, [1.0], density = rho, temperature = temperature_inj, check = false) # Producer # TODO: support multiple periods with different rates rate_target = TotalRateTarget(-rate) diff --git a/src/meshing/extruded.jl b/src/meshing/extruded.jl index 0c15eae..8be1df7 100644 --- a/src/meshing/extruded.jl +++ b/src/meshing/extruded.jl @@ -4,12 +4,15 @@ function extruded_mesh(cell_constraints, depths; interpolation = :default, dist_min_factor = 1.1, dist_max_factor = 0.75, recombine_to_quads = true, + info_level = -1 ) @assert depths[1] == 0.0 # Initialize Gmsh gmsh.initialize() + gmsh.option.setNumber("General.Terminal", info_level < 0 ? 0 : 1) + gmsh.option.setNumber("General.Verbosity", max(info_level, 0)) gmsh.clear() gmsh.model.add("extruded_mesh") diff --git a/src/wells/btes.jl b/src/wells/btes.jl deleted file mode 100644 index ad3e032..0000000 --- a/src/wells/btes.jl +++ /dev/null @@ -1,298 +0,0 @@ -using LinearAlgebra - -## Utility functions for setting up BTES wells -function setup_btes_well(D::DataDomain, reservoir_cells; - btes_type = :simple, - name = :BTES, kwarg...) - - if btes_type == :simple - # Simple BTES well with direct pipe/reservoir thermal communication - return setup_btes_well_simple(D, reservoir_cells; name = name, kwarg...) - elseif btes_type == :u1 - # U-type BTES well with grout annulus - return setup_btes_well_u1(D, reservoir_cells; name = name, kwarg...) - else - # Unknown BTES type - error("Unknown BTES type: $btes_type") - end - -end - -function setup_vertical_btes_well(D::DataDomain, i, j; - heel = 1, toe = missing, kwarg...) - - # Get reservoir cells from ijk-indices - g = physical_representation(D) - if ismissing(toe) - toe = grid_dims_ijk(g)[3] - end - @assert heel <= toe - @assert heel > 0 - @assert toe > 0 - k_range = heel:toe - n = length(k_range) - @assert n > 0 - reservoir_cells = zeros(Int64, n) - for (ix, k) in enumerate(k_range) - reservoir_cells[ix] = cell_index(g, (i, j, k)) - end - # Set up BTES well - return setup_btes_well(D, reservoir_cells; kwarg...) -end - -function setup_btes_well_simple(D::DataDomain, reservoir_cells; - name = :BTES, - radius_pipe = 20e-3, - pipe_thickness = 2.5e-3, - grouting_thickness = 50e-3, - thermal_conductivity_pipe = 0.35, - kwarg...) - - # Common properties - args = ( - WI = 0.0, - radius = radius_pipe, - casing_thickness = pipe_thickness, - grouting_thickness = grouting_thickness, - thermal_conductivity_casing = thermal_conductivity_pipe, - end_nodes = [length(reservoir_cells)], - simple_well = false, - type = :closed_loop - ) - - # Set up supply and return wells - supply_well = setup_well(D::DataDomain, reservoir_cells; - name = Symbol(name, "_supply"), - args..., kwarg...) - return_well = setup_well(D::DataDomain, reservoir_cells; - name = Symbol(name, "_return"), - args..., kwarg...) - - return supply_well, return_well - -end - - -function setup_btes_well_u1(D::DataDomain, reservoir_cells; - cell_centers = D[:cell_centroids], - name = :BTES, - radius_grout = 65e-3, - radius_pipe_outer = 15e-3, - pipe_thickness = 3e-3, - pipe_spacing = 60e-3, - thermal_conductivity_grout = 2.3, - thermal_conductivity_pipe = 0.38, - heat_capacity_grout = 420.0, - density_grout = 1500.0, - kwarg...) - - ## Set up connectivity - nc_pipe = length(reservoir_cells) - nc_grout = length(reservoir_cells) - - pipe_cells = (1:nc_pipe) - grout_cells = (1:nc_grout) .+ nc_pipe - - pipe_to_pipe = vcat(pipe_cells[1:end-1]', pipe_cells[2:end]') - pipe_to_grout = vcat(pipe_cells', grout_cells') - grout_to_grout = vcat(grout_cells[1:end-1]', grout_cells[2:end]') - - N = hcat(pipe_to_pipe, pipe_to_grout) - - ## Set up segment flow models - segment_models = Vector{Any}() - - # Set centers and depths - well_cell_centers = repeat(cell_centers[:, reservoir_cells], 1, 2) - - # Add top node - nseg = size(N,2) - - # Set material thermal conducivities - nr = length(reservoir_cells) - - for seg in 1:nseg - if seg < nc_pipe - # Pipe segments use standard wellbore friction model - seg_model = SegmentWellBoreFrictionHB() - else - seg_model = JutulDarcy.ClosedSegment() - end - push!(segment_models, seg_model) - end - - # nc = nc_pipe + nc_grout + 1 - - ## Set up supply and return wells - args = ( - type = :closed_loop, - simple_well = false, - neighborship = N, - perforation_cells_well = collect(grout_cells), - well_cell_centers = well_cell_centers, - radius = radius_pipe_outer, - casing_thickness = pipe_thickness, - WI = 0.0, - thermal_conductivity_casing = thermal_conductivity_pipe, - thermal_conductivity_grout = thermal_conductivity_grout, - segment_models = segment_models, - end_nodes = [nc_pipe], - ) - - supply_well = setup_well(D, reservoir_cells; - name = Symbol(name, "_supply"), - args..., kwarg...) - augment_btes_domain!(supply_well, - radius_grout, - pipe_spacing, - heat_capacity_grout, - density_grout - ) - - return_well = setup_well(D, reservoir_cells; - name = Symbol(name, "_return"), - args...) - augment_btes_domain!(return_well, - radius_grout, - pipe_spacing, - heat_capacity_grout, - density_grout - ) - - set_default_btes_thermal_indices!(supply_well) - set_default_btes_thermal_indices!(return_well) - - return [supply_well, return_well] - -end - -function augment_btes_domain!(well::DataDomain, - radius_grout, - pipe_spacing, - heat_capacity_grout, - density_grout; - pipe_grout_thermal_index = missing, - grout_grout_thermal_index = missing -) - - c = Cells() - p = Perforations() - f = Faces() - - well[:radius_grout, p] = radius_grout - well[:pipe_spacing, p] = pipe_spacing - well[:heat_capacity_grout, p] = heat_capacity_grout - well[:density_grout, p] = density_grout - - treat_defaulted(x) = x - treat_defaulted(::Missing) = NaN - treat_defaulted(::Nothing) = NaN - - λpg = treat_defaulted(pipe_grout_thermal_index) - λgg = treat_defaulted(grout_grout_thermal_index) - well[:pipe_grout_thermal_index, p] = λpg - well[:grout_grout_thermal_index, p] = λgg - -end - -function set_default_btes_thermal_indices!(well::DataDomain) - - cell = Cells() - face = Faces() - perf = Perforations() - - volumes = zeros(Float64, well.representation.num_nodes) - - rep = physical_representation(well) - for (pno, cno) in enumerate(rep.perforations.self) - r_grout = well[:radius_grout, perf][pno] - r_pipe_outer = well[:perforation_radius, perf][pno] - r_pipe_inner = r_pipe_outer - well[:casing_thickness, cell][cno] - length = well[:cell_length, cell][cno] - vol_p, vol_g = btes_volume( - BTESTypeU1(), length, r_grout, r_pipe_inner, r_pipe_outer - ) - - volumes[pno] = vol_g - seg = findall(well.representation.neighborship[2, :] .== cno) - pcno = well.representation.neighborship[1, seg[1]] - println("pcno = $pcno") - volumes[pcno[1]] = vol_p - - pipe_spacing = well[:pipe_spacing, perf][pno] - λg = well[:thermal_conductivity_grout, cell][cno] - λp = well[:thermal_conductivity_casing, cell][cno] - λpg, λgr, λgg = btes_thermal_conductivity( - BTESTypeU1(), - r_grout, r_pipe_inner, r_pipe_outer, pipe_spacing, length, λg, λp) - - if well[:material_thermal_conductivity, face][seg] == 0.0 - well[:material_thermal_conductivity, face][seg] = λpg - end - - if isnan(well[:thermal_well_index, perf][pno]) - well[:thermal_well_index, perf][pno] = λgr - end - - if isnan(well[:grout_grout_thermal_index, perf][pno]) - well[:grout_grout_thermal_index, perf][pno] = λgg - end - - println("Perforation $pno: vol_p = $vol_p, vol_g = $vol_g, L = $length "* - "λpg = $λpg, λgr = $λgr, λgg = $λgg") - end - -end - -# Conveience types for multiple dispatching -abstract type AbstractBTESType end -struct BTESTypeU1 <: AbstractBTESType end - -function btes_volume(type::BTESTypeU1, length, radius_grout, radius_pipe_inner, radius_pipe_outer) - - # Compute pipe and grout volume - L = length - rg, rp_in, rp_out = radius_grout, radius_pipe_inner, radius_pipe_outer - vol_p = π*rp_in^2*L - vol_g = π*rg^2*L/2 - vol_p - - return vol_p, vol_g, L - -end - -function btes_thermal_conductivity(type::BTESTypeU1, - radius_grout, radius_pipe_inner, radius_pipe_outer, pipe_spacing, length, - thermal_conductivity_grout, thermal_conductivity_pipe) - # Conveient short-hand notation - rg, rp_in, rp_out, w, L = - radius_grout, radius_pipe_inner, radius_pipe_outer, pipe_spacing, length - λg, λp = thermal_conductivity_grout, thermal_conductivity_pipe - - ## Compute thermal resistances - # Advection-dependent pipe thermal resistance - Ra = 0.0 # TODO: Implement this - # Conduction-dependent pipe thermal resistance - Rc_a = log(rp_out/rp_in)/(2*π*λp) - dg, dp_out = 2*rg, 2*rp_out - x = log(sqrt(dg^2 + 2*dp_out^2)/(2*dp_out))/log(dg/(sqrt(2)*dp_out)) - # Grout thermal resistance - Rg = acosh((dg^2 + dp_out^2 - w^2)/(2*dg*dp_out))/(2*π*λg)*(1.601 - 0.888*w/dg) - # Conduction-dependent grout thermal resistance - Rc_b = x*Rg - # Combined thermal resistance of pipe and grout - Rpg = Ra + Rc_a + Rc_b - Rgr = (1-x)*Rg - # Grout-to-grout thermal resistance - Rar = acosh((2*w^2 - dp_out^2)/dp_out^2)/(2*π*λg) - Rgg = 2*Rgr*(Rar - 2*x*Rg)/(2*Rgr - Rar + 2*x*Rg) - - println("Rpg = $Rpg, Rgr = $Rgr, Rgg = $Rgg") - - # Compute thermal conducivities - λpg = L*2*π*rp_in/Rpg - λgr = L*π*rg/Rgr - λgg = L*2*rg/Rgg - - return λpg, λgr, λgg - -end diff --git a/src/wells/closed_loop.jl b/src/wells/closed_loop.jl new file mode 100644 index 0000000..c43a67f --- /dev/null +++ b/src/wells/closed_loop.jl @@ -0,0 +1,140 @@ +using LinearAlgebra + +# Utility functions for setting up BTES wells +function setup_closed_loop_well(D::DataDomain, reservoir_cells; + closed_loop_type = :simple, + name = :CL, kwarg...) + + if closed_loop_type == :simple + # Simple closed loop well with direct pipe/reservoir thermal communication + return setup_closed_loop_well_simple(D, reservoir_cells; name = name, kwarg...) + elseif closed_loop_type == :u1 + # U-type closed loop well with grout annulus + return setup_closed_loop_well_u1(D, reservoir_cells; name = name, kwarg...) + elseif closed_loop_type == :coaxial + # Coaxial closed loop well with grout annulus + return setup_closed_loop_well_coaxial(D, reservoir_cells; name = name, kwarg...) + else + # Unknown closed loop type + error("Unknown closed loop type: $closed_loop_type") + end + +end + +function setup_vertical_closed_loop_well(D::DataDomain, i, j; + heel = 1, toe = missing, kwarg...) + + # Get reservoir cells from ijk-indices + g = physical_representation(D) + if ismissing(toe) + toe = grid_dims_ijk(g)[3] + end + @assert heel <= toe + @assert heel > 0 + @assert toe > 0 + k_range = heel:toe + n = length(k_range) + @assert n > 0 + reservoir_cells = zeros(Int64, n) + for (ix, k) in enumerate(k_range) + reservoir_cells[ix] = cell_index(g, (i, j, k)) + end + # Set up BTES well + return setup_closed_loop_well(D, reservoir_cells; kwarg...) +end + +function setup_btes_well(D::DataDomain, reservoir_cells; + name = :BTES, kwarg...) + + return setup_closed_loop_well(D, reservoir_cells; + name = name, kwarg...) + +end + +function setup_vertical_btes_well(D::DataDomain, i, j; + name = :BTES, heel = 1, toe = missing, kwarg...) + + return setup_vertical_closed_loop_well(D, i, j; + name = name, heel = heel, toe = toe, kwarg...) + +end + +function setup_closed_loop_well_simple(D::DataDomain, reservoir_cells; + return_reservoir_cell = missing, + cell_centers = D[:cell_centroids], + neighborship = missing, + well_cell_centers = missing, + segment_models = missing, + end_nodes = missing, + section = missing, + name = :CL, + radius_pipe = 20e-3, + wall_thickness = 2.5e-3, + grouting_thickness = 50e-3, + pipe_thermal_conductivity = 0.35, + grouting_thermal_conductivity = 2.3, + kwarg...) + + if ismissing(neighborship) + # Create pipe cells + reservoir_cells = vcat(reservoir_cells, reverse(reservoir_cells)) + nc_pipe = length(reservoir_cells) + pipe_cells = collect(1:nc_pipe) + # Set up connectivity + pipe_to_pipe = vcat(pipe_cells[1:end-1]', pipe_cells[2:end]') + neighborship = hcat(pipe_to_pipe) + # Set up well cell centers and end nodes + nc_mid = div(nc_pipe, 2) + left_ix = Int.(1:nc_mid) + right_ix = Int.(1:nc_mid) .+ nc_mid + pipe_spacing = 60e-3 + wc_left = cell_centers[:, reservoir_cells[1:nc_mid]] .- [pipe_spacing/2, 0.0, 0.0] + wc_right = cell_centers[:, reservoir_cells[nc_mid+1:end]] .+ [pipe_spacing/2, 0.0, 0.0] + well_cell_centers = hcat(wc_left, wc_right) + end_nodes = [pipe_cells[end]] + # Set section for easy lookup + if !ismissing(section) + @warn(["section argument is ignored when neighborship is not provided. ", + "Sections will be created automatically."]) + end + section = Vector{Any}(undef, nc_pipe) + section[pipe_cells[left_ix]] .= [(1, :pipe_left)] + section[pipe_cells[right_ix]] .= [(1, :pipe_right)] + else + if ismissing(well_cell_centers) || ismissing(end_nodes) + error("""If neighborship is provided, well_cell_centers, and + end_nodes must also be provided.""") + end + end + # Set return reservoir cell if missing + if ismissing(return_reservoir_cell) + return_reservoir_cell = reservoir_cells[end] + end + + # Common properties + args = ( + WI = 0.0, + radius = radius_pipe, + simple_well = false, + type = :closed_loop + ) + + # Set up supply and return wells + supply_well = setup_well(D::DataDomain, reservoir_cells; + name = Symbol(name, "_supply"), + casing_thickness = wall_thickness, + grouting_thickness = grouting_thickness, + casing_thermal_conductivity = pipe_thermal_conductivity, + grouting_thermal_conductivity = grouting_thermal_conductivity, + end_nodes = end_nodes, + args..., kwarg...) + + return_well = setup_well(D::DataDomain, return_reservoir_cell; + name = Symbol(name, "_return"), + WIth = 0.0, + casing_thickness = wall_thickness, + args..., kwarg...) + + return [supply_well, return_well] + +end \ No newline at end of file diff --git a/src/wells/closed_loop_analytical.jl b/src/wells/closed_loop_analytical.jl new file mode 100644 index 0000000..c3ba5f3 --- /dev/null +++ b/src/wells/closed_loop_analytical.jl @@ -0,0 +1,233 @@ +function temperature_closed_loop_pipe(u, T_in, T_rock, ρ, Cp, A, L, R1Δ, R2Δ, R12Δ) + + f1, f2, f3, f4, f5 = f_functions(u, ρ, Cp, A, R1Δ, R2Δ, R12Δ) + + domain = z -> (0.0, z) + + integrand_supp = z -> (ξ,p) -> T_rock*f4(z - ξ) + I_supp = z -> solve(IntegralProblem(integrand_supp(z), domain(z)), QuadGKJL()) + + integrand_ret = z -> (ξ,p) -> T_rock*f5(z - ξ) + I_ret = z -> solve(IntegralProblem(integrand_ret(z), domain(z)), QuadGKJL()) + + # Find T_out from T_supply(L) = T_return(L) + T_out = (T_in*(f1(L) .+ f2(L)) .+ I_supp(L) .+ I_ret(L))/(f3(L) .- f2(L)) + T_supply = z -> T_in*f1(z) .+ T_out*f2(z) .+ I_supp(z) + T_return = z -> -T_in*f2(z) .+ T_out*f3(z) .- I_ret(z) + + return T_supply, T_return + +end + +function f_functions(u, ρ, Cp, A, R1Δ, R2Δ, R12Δ) + + β1 = 1/(R1Δ*A*ρ*Cp*u) + β2 = 1/(R2Δ*A*ρ*Cp*u) + β12 = 1/(R12Δ*A*ρ*Cp*u) + β = (β2 - β1)/2 + γ = sqrt((β1 + β2)^2/4 + β12*(β1 + β2)) + δ = 1/γ*(β12 + (β1 + β2)/2) + + f1 = z -> exp(β*z)*(cosh(γ*z) - δ*sinh(γ*z)) + f2 = z -> exp(β*z)*β12/γ*sinh(γ*z) + f3 = z -> exp(β*z)*(cosh(γ*z) + δ*sinh(γ*z)) + f4 = z -> exp(β*z)*(β1*cosh(γ*z) - (δ*β1 + β2*β12/γ)*sinh(γ*z)) + f5 = z -> exp(β*z)*(β2*cosh(γ*z) + (δ*β2 + β1*β12/γ)*sinh(γ*z)) + return (f1, f2, f3, f4, f5) +end + +function analytical_closed_loop_u1(rate, temperature_in, temperature_rock, + density_fluid, heat_capacity_fluid, well::DataDomain) + + # Identify pipe and grout cells + section = well[:section] + is_pipe_cell = last.(section) .== :pipe_left + is_grout_cell = last.(section) .== :grout_left + # Get length + length = sum(well[:cell_length][is_pipe_cell]) + # Get pipe outer radius + rp_i = well[:radius][is_pipe_cell][1] + wall_thickness_pipe = well[:casing_thickness][is_pipe_cell][1] + radius_pipe = rp_i + wall_thickness_pipe + # Get grout radius + L = well[:cell_length][is_grout_cell][1] + vg = well[:volume_override_grouting][is_grout_cell][1] + vp = radius_pipe^2*π*L + radius_grout = sqrt((2*vg + 2*vp)/(π*L)) + # Get other parameters + pipe_spacing = maximum(well[:pipe_spacing]) + thermal_conductivity_grout = maximum(well[:grouting_thermal_conductivity]) + thermal_conductivity_pipe = maximum(well[:casing_thermal_conductivity]) + + return analytical_closed_loop_u1( + rate, temperature_in, temperature_rock, + density_fluid, heat_capacity_fluid, + length, radius_grout, radius_pipe, wall_thickness_pipe, pipe_spacing, + thermal_conductivity_grout, thermal_conductivity_pipe) + +end + +function analytical_closed_loop_u1(rate, temperature_in, temperature_rock, + density_fluid, heat_capacity_fluid, + length, radius_grout, radius_pipe, wall_thickness_pipe, pipe_spacing, + thermal_conductivity_grout, thermal_conductivity_pipe) + + Rpg, Rgr, Rgg = Fimbul.closed_loop_thermal_resistance_u1( + radius_grout, radius_pipe, wall_thickness_pipe, pipe_spacing, + thermal_conductivity_grout, thermal_conductivity_pipe) + A = π*(radius_pipe - wall_thickness_pipe)^2 + + u = rate/A + Tps, Tpr = temperature_pipe_u1( + u, temperature_in, temperature_rock, density_fluid, heat_capacity_fluid, A, length, Rpg, Rgr, Rgg) + + Tgs, Tgr = temperature_grout_u1( + Tps, Tpr, temperature_rock, Rpg, Rgr, Rgg) + + sol = Dict() + sol[:pipe_left] = Tps + sol[:pipe_right] = Tpr + sol[:grout_left] = Tgs + sol[:grout_right] = Tgr + + return sol + +end + +function temperature_pipe_u1(u, T_in, T_rock, ρ, Cp, A, L, Rpg, Rgr, Rgg) + + u1 = 1/Rpg + 1/Rgr + 1/Rgg + R1Δ = Rpg + Rgr + R2Δ = Rpg + Rgr + R12Δ = ((u1*Rpg*Rgg)^2 - Rpg^2)/Rgg + + return temperature_closed_loop_pipe(u, T_in, T_rock, ρ, Cp, A, L, R1Δ, R2Δ, R12Δ) + +end + +function temperature_grout_u1(T_supply, T_return, T_rock, Rpg, Rgr, Rgg) + + u1 = 1/Rpg + 1/Rgr + 1/Rgg + Tgs = z -> ( + T_rock/Rgr + + T_return(z)/Rpg + + (T_rock/Rgr + T_supply(z)/Rpg)*u1*Rgg + )*Rgg/ + ((Rgg*u1)^2-1) + + Tgr = z -> ( + Tgs(z)/Rgg + + T_return(z)/Rpg + + T_rock/Rgr + )*1/u1 + + return Tgs, Tgr + +end + +function analytical_closed_loop_coaxial(rate, temperature_in, temperature_rock, + density_fluid, heat_capacity_fluid, well::DataDomain; kwargs...) + + # Identify pipe and grout cells + section = well[:section] + is_inner_pipe_cell = last.(section) .== :pipe_inner + is_outer_pipe_cell = last.(section) .== :pipe_outer + is_grout_cell = last.(section) .== :grout + # Get length + length = sum(well[:cell_length][is_inner_pipe_cell]) + # Get pipe radii and wall thicknesses + radius_inner_pipe = well[:radius_inner][is_outer_pipe_cell][1] + wall_thickness_inner_pipe = well[:casing_thickness][is_inner_pipe_cell][1] + rpo_i = well[:radius][is_outer_pipe_cell][1] + wall_thickness_outer_pipe = well[:casing_thickness][is_outer_pipe_cell][1] + radius_outer_pipe = rpo_i + wall_thickness_outer_pipe + # Get grout radius + L = well[:cell_length][is_grout_cell][1] + vg = well[:volume_override_grouting][is_grout_cell][1] + radius_grout = sqrt((vg + radius_outer_pipe^2*pi*L)/(π*L)) + # Get other parameters + thermal_conductivity_grout = maximum(well[:grouting_thermal_conductivity]) + thermal_conductivity_inner_pipe = maximum(well[:casing_thermal_conductivity][is_inner_pipe_cell]) + thermal_conductivity_outer_pipe = maximum(well[:casing_thermal_conductivity][is_outer_pipe_cell]) + + return analytical_closed_loop_coaxial(rate, temperature_in, temperature_rock, + density_fluid, heat_capacity_fluid, + length, radius_grout, + radius_inner_pipe, wall_thickness_inner_pipe, + radius_outer_pipe, wall_thickness_outer_pipe, + thermal_conductivity_grout, + thermal_conductivity_inner_pipe, + thermal_conductivity_outer_pipe; + kwargs...) + +end + +function analytical_closed_loop_coaxial(rate, temperature_in, temperature_rock, + density_fluid, heat_capacity_fluid, + length, radius_grout, + radius_inner_pipe, wall_thickness_inner_pipe, + radius_outer_pipe, wall_thickness_outer_pipe, + thermal_conductivity_grout, + thermal_conductivity_inner_pipe, + thermal_conductivity_outer_pipe; + inlet = :outer) + + Rpp, Rpg, Rgr = closed_loop_thermal_resistance_coaxial( + radius_grout, + radius_inner_pipe, wall_thickness_inner_pipe, + radius_outer_pipe, wall_thickness_outer_pipe, + thermal_conductivity_grout, + thermal_conductivity_inner_pipe, + thermal_conductivity_outer_pipe + ) + if inlet == :outer + A = π*((radius_outer_pipe - wall_thickness_outer_pipe)^2 - radius_inner_pipe^2) + elseif inlet == :inner + A = π*(radius_inner_pipe - wall_thickness_inner_pipe)^2 + else + error("Inlet must be :outer or :inner") + end + u = rate/A + Tps, Tpr = temperature_pipe_coaxial( + u, temperature_in, temperature_rock, + density_fluid, heat_capacity_fluid, + A, length, Rpg, Rgr, Rpp, inlet + ) + + To = inlet == :inner ? Tpr : Tps + Ti = inlet == :inner ? Tps : Tpr + Tg = temperature_grout_coaxial( + To, temperature_rock, Rpg, Rgr) + + sol = Dict() + sol[:pipe_inner] = Ti + sol[:pipe_outer] = To + sol[:grout] = Tg + + return sol + +end + +function temperature_pipe_coaxial(u, T_in, T_rock, ρ, Cp, A, L, Rpg, Rgr, Rpp, inlet) + + if inlet == :outer + R1Δ = Rpg + Rgr + R2Δ = Inf + elseif inlet == :inner + R1Δ = Inf + R2Δ = Rpg + Rgr + else + error("Inlet must be :outer or :inner") + end + R12Δ = Rpp + + return temperature_closed_loop_pipe(u, T_in, T_rock, ρ, Cp, A, L, R1Δ, R2Δ, R12Δ) + +end + +function temperature_grout_coaxial(T_annulus, T_rock, Rpg, Rgr) + + T_grout = z -> Rpg/(Rpg + Rgr)*(T_rock - T_annulus(z)) + T_annulus(z) + return T_grout + +end \ No newline at end of file diff --git a/src/wells/closed_loop_coaxial.jl b/src/wells/closed_loop_coaxial.jl new file mode 100644 index 0000000..b1e9c1c --- /dev/null +++ b/src/wells/closed_loop_coaxial.jl @@ -0,0 +1,315 @@ +function setup_closed_loop_well_coaxial(D::DataDomain, reservoir_cells; + name = :CL, + return_reservoir_cell = reservoir_cells[1], + cell_centers = D[:cell_centroids], + neighborship = missing, + pipe_cells_inner = missing, + pipe_cells_outer = missing, + grout_cells = missing, + section = missing, + well_cell_centers = missing, + segment_models = missing, + end_nodes = missing, + radius_grout = 65e-3, + radius_pipe_inner = 15e-3, + radius_pipe_outer = 20e-3, + wall_thickness_pipe_inner = 3e-3, + wall_thickness_pipe_outer = 3e-3, + grouting_thermal_conductivity = 2.3, + inner_pipe_thermal_conductivity = 0.38, + outer_pipe_thermal_conductivity = 0.38, + grouting_heat_capacity = 1460, + grouting_density = 1500.0, + kwarg...) + + + if ismissing(neighborship) + if !ismissing(pipe_cells_inner) && !ismissing(pipe_cells_outer) && !ismissing(grout_cells) + error("""Provide either neighborship, pipe_cells_inner, + pipe_cells_outer, grout_cells, "well_cell_centers, end_nodes, or + none of them.""") + end + # Set inner/outer pipe and grout cells + nc_pipe = length(reservoir_cells)*2 + nc_grout = length(reservoir_cells) + nc_r = length(reservoir_cells) + pipe_cells_inner = collect(1:nc_r) + pipe_cells_outer = pipe_cells_inner .+ nc_r + grout_cells = pipe_cells_outer .+ nc_r + # Set connectivity + num_cells = nc_pipe + nc_grout + pi2pi = vcat(pipe_cells_inner[1:end-1]', pipe_cells_inner[2:end]') + po2po = vcat(pipe_cells_outer[1:end-1]', pipe_cells_outer[2:end]') + pi2po_bottom = [pipe_cells_inner[end]', pipe_cells_outer[end]'] + pi2po = vcat(pipe_cells_inner[1:end-1]', pipe_cells_outer[1:end-1]') + po2g = vcat(pipe_cells_outer', grout_cells') + neighborship = hcat(pi2pi, po2po, pi2po_bottom, pi2po, po2g) + # Set centers and end nodes + well_cell_centers = repeat(cell_centers[:, reservoir_cells], 1, 3) + end_nodes = [nc_r+1] + if !ismissing(section) + warn(["section argument is ignored when neighborship is not provided. ", + "Sections will be created automatically."]) + end + section = Vector{Any}(undef, nc_pipe + nc_grout) + section[pipe_cells_inner] .= [(1, :pipe_inner)] + section[pipe_cells_outer] .= [(1, :pipe_outer)] + section[grout_cells] .= [(1, :grout)] + end + + if ismissing(segment_models) + # Set up segment flow models + segment_models = Vector{Any}() + flow_segs = size(pi2pi, 2) + size(po2po, 2) + size(pi2po_bottom, 2) + nseg = size(neighborship,2) + for s in 1:nseg + if s <= flow_segs + # Pipe segments use standard wellbore friction model + seg_model = SegmentWellBoreFrictionHB() + else + seg_model = JutulDarcy.ClosedSegment() + end + push!(segment_models, seg_model) + end + end + + # Set cell radii + cell_radius = zeros(num_cells) + cell_radius_inner = zeros(num_cells) + cell_radius[pipe_cells_inner] .= radius_pipe_inner - wall_thickness_pipe_inner + cell_radius[pipe_cells_outer] .= radius_pipe_outer - wall_thickness_pipe_outer + cell_radius_inner[pipe_cells_outer] .= radius_pipe_inner + # Set casing thicknesses + wall_thickness = zeros(num_cells) + wall_thickness[pipe_cells_inner] .= wall_thickness_pipe_inner + wall_thickness[pipe_cells_outer] .= wall_thickness_pipe_outer + # Set pipe wall thermal conductivities + λ_pipe = zeros(num_cells) + λ_pipe[pipe_cells_inner] .= inner_pipe_thermal_conductivity + λ_pipe[pipe_cells_outer] .= outer_pipe_thermal_conductivity + # Set grouting thermal conductivities + λ_grout = zeros(num_cells) + λ_grout[grout_cells] .= grouting_thermal_conductivity + + # Common well arguments + args = ( + type = :closed_loop, + simple_well = false, + WI = 0.0 + ) + + # Setup supply well + supply_well = setup_well(D, reservoir_cells; + name = Symbol(name, "_supply"), + neighborship = neighborship, + perforation_cells_well = collect(grout_cells), + well_cell_centers = well_cell_centers, + cell_radius = cell_radius, + cell_radius_inner = cell_radius_inner, + casing_thickness = wall_thickness, + casing_thermal_conductivity = λ_pipe, + grouting_thermal_conductivity = λ_grout, + segment_models = segment_models, + end_nodes = end_nodes, + args..., kwarg...) + # Augment supply well with closed loop specific data + augment_closed_loop_domain_coaxial!(supply_well, + radius_grout, + grouting_heat_capacity, + grouting_density, + section = section + ) + # Set default thermal indices + set_default_closed_loop_thermal_indices_coaxial!(supply_well) + + # Setup return well + return_well = setup_well(D, return_reservoir_cell; + name = Symbol(name, "_return"), + WIth = 0.0, + args...) + + return [supply_well, return_well] + +end + +function augment_closed_loop_domain_coaxial!(well::DataDomain, + radius_grout, + grouting_heat_capacity, + grouting_density; + pipe_pipe_thermal_index = missing, + pipe_grout_thermal_index = missing, + section = missing +) + + c = Cells() + p = Perforations() + f = Faces() + + well[:radius_grout, c] = radius_grout + well[:material_heat_capacity, c] = grouting_heat_capacity + well[:material_density, c] = grouting_density + + treat_defaulted(x) = x + treat_defaulted(::Missing) = NaN + treat_defaulted(::Nothing) = NaN + + λpp = treat_defaulted(pipe_pipe_thermal_index) + λpg = treat_defaulted(pipe_grout_thermal_index) + well[:pipe_pipe_thermal_index, p] = λpp + well[:pipe_grout_thermal_index, p] = λpg + + if !ismissing(section) + well[:section, c] = section + end + +end + +function set_default_closed_loop_thermal_indices_coaxial!(well::DataDomain) + + cell = Cells() + face = Faces() + perf = Perforations() + + num_nodes = well.representation.num_nodes + + hole_volumes = zeros(Float64, num_nodes) + casing_volumes = zeros(Float64, num_nodes) + grout_volumes = zeros(Float64, num_nodes) + + nc = Int.(num_nodes/3) + pipe_cells_inner = Int.(1:nc) + pipe_cells_outer = Int.(1:nc) .+ nc + grout_cells = Int.(1:nc) .+ 2*nc + N = well.representation.neighborship + rep = physical_representation(well) + for (pno, gc) in enumerate(rep.perforations.self) + + @assert gc ∈ grout_cells + seg_pg = findall(N[2, :] .== gc) + @assert length(seg_pg) == 1 + seg_pg = seg_pg[1] + pc_out = N[1, seg_pg] + @assert pc_out ∈ pipe_cells_outer + seg_pp = findall(N[2, :] .== pc_out) + @assert length(seg_pp) <= 2 + seg_pp = seg_pp[end] + pc_in = N[1, seg_pp] + @assert pc_in ∈ pipe_cells_inner + + r_grout = well[:radius_grout, cell][gc] + r_inner_pipe = well[:radius, cell][pc_in] + wall_thickness_inner = well[:casing_thickness, cell][pc_in] + r_outer_pipe = well[:radius, cell][pc_out] + wall_thickness_outer = well[:casing_thickness, cell][pc_out] + L = well[:cell_length, cell][pc_in] + vol_ip, vol_iw, vol_op, vol_ow, vol_g, L = closed_loop_volume_coaxial( + L, r_grout, + r_inner_pipe + wall_thickness_inner, wall_thickness_inner, + r_outer_pipe + wall_thickness_outer, wall_thickness_outer + ) + + hole_volumes[pc_in] = vol_ip + hole_volumes[pc_out] = vol_op + hole_volumes[gc] = 0.0 + casing_volumes[pc_in] = vol_iw + casing_volumes[pc_out] = vol_ow + casing_volumes[gc] = 0.0 + grout_volumes[pc_in] = 0.0 + grout_volumes[pc_out] = 0.0 + grout_volumes[gc] = vol_g + + λg = well[:grouting_thermal_conductivity, cell][gc] + λpi = well[:casing_thermal_conductivity, cell][pc_in] + λpo = well[:casing_thermal_conductivity, cell][pc_out] + Rpp, Rpg, Rgr = closed_loop_thermal_resistance_coaxial( + r_grout, + r_inner_pipe + wall_thickness_inner, wall_thickness_inner, + r_outer_pipe + wall_thickness_outer, wall_thickness_outer, + λg, λpi, λpo) + λpp = L/Rpp + λpg = L/Rpg + λgr = L/Rgr + + if isnan(well[:thermal_well_index, perf][pno]) + well[:thermal_well_index, perf][pno] = λgr + end + + if well[:material_thermal_conductivity, face][seg_pg] == 0.0 + well[:material_thermal_conductivity, face][seg_pg] = λpg + end + + if well[:material_thermal_conductivity, face][seg_pp] == 0.0 + well[:material_thermal_conductivity, face][seg_pp] = λpp + end + end + + well[:volume_override_hole, cell] = hole_volumes + well[:volume_override_casing, cell] = casing_volumes + well[:volume_override_grouting, cell] = grout_volumes +end + +function closed_loop_volume_coaxial( + length, + radius_grout, + radius_inner_pipe, wall_thickness_inner_pipe, + radius_outer_pipe, wall_thickness_outer_pipe + ) + + # Compute pipe and grout volume + r_g, r_ipi, r_ipo, r_opi, r_opo, L = + radius_grout, + radius_inner_pipe - wall_thickness_inner_pipe, + radius_inner_pipe, + radius_outer_pipe - wall_thickness_outer_pipe, + radius_outer_pipe, + length + + vol_hole_inner = π*r_ipi^2*L + vol_wall_inner = π*(r_ipo^2 - r_ipi^2)*L + vol_hole_outer = π*(r_opi^2 - r_ipo^2)*L + vol_wall_outer = π*(r_opo^2 - r_opi^2)*L + vol_grout = π*(r_g^2 - r_opo^2)*L + + return vol_hole_inner, vol_wall_inner, vol_hole_outer, vol_wall_outer, vol_grout, L + +end + +function closed_loop_thermal_resistance_coaxial( + radius_grout, + radius_inner_pipe, wall_thickness_inner_pipe, + radius_outer_pipe, wall_thickness_outer_pipe, + thermal_conductivity_grout, + thermal_conductivity_inner_pipe, + thermal_conductivity_outer_pipe + ) + + r_g, r_ipi, r_ipo, r_opi, r_opo = + radius_grout, + radius_inner_pipe - wall_thickness_inner_pipe, + radius_inner_pipe, + radius_outer_pipe - wall_thickness_outer_pipe, + radius_outer_pipe + + λ_g, λ_pi, λ_po = + thermal_conductivity_grout, + thermal_conductivity_inner_pipe, + thermal_conductivity_outer_pipe + + R_ai = 0.0 # TODO: Implement this + R_ao = 0.0 + + R_ci = log(r_ipo/r_ipi)/(2*π*λ_pi) + R_co = log(r_opo/r_opi)/(2*π*λ_po) + + d_g, d_opo = 2*r_g, 2*r_opo + x = log(sqrt((d_g^2 + d_opo^2)/(2*d_opo^2)))/log(d_g/d_opo) + R_g = log(d_g/d_opo)/(2*π*λ_g) + R_cg = x*R_g + + Rpp = R_ai + R_ci + Rpg = R_ao + R_co + R_cg + Rgr = (1-x)*R_g + + return Rpp, Rpg, Rgr + +end \ No newline at end of file diff --git a/src/wells/closed_loop_u1.jl b/src/wells/closed_loop_u1.jl new file mode 100644 index 0000000..2ab2185 --- /dev/null +++ b/src/wells/closed_loop_u1.jl @@ -0,0 +1,307 @@ +function setup_closed_loop_well_u1(D::DataDomain, reservoir_cells; + name = :CL, + return_reservoir_cell = missing, + cell_centers = D[:cell_centroids], + neighborship = missing, + pipe_cells = missing, + grout_cells = missing, + section = missing, + well_cell_centers = missing, + segment_models = missing, + end_nodes = missing, + radius_grout = 65e-3, + radius_pipe = 15e-3, + wall_thickness_pipe = 3e-3, + pipe_spacing = 60e-3, + grouting_thermal_conductivity = 2.3, + pipe_thermal_conductivity = 0.38, + grouting_heat_capacity = 1500, + grouting_density = 2000.0, + kwarg...) + + if ismissing(neighborship) + if !ismissing(pipe_cells) && !ismissing(grout_cells) + error(["Provide either neighborship, pipe_cells, grout_cells, ", + "well_cell_centers, end_nodes, or none of them."]) + end + # Create grout and pipe cells + reservoir_cells = vcat(reservoir_cells, reverse(reservoir_cells)) + nc_pipe = length(reservoir_cells) + nc_grout = length(reservoir_cells) + nc_mid = div(nc_pipe, 2) + left_ix = Int.(1:nc_mid) + right_ix = Int.(1:nc_mid) .+ nc_mid + pipe_cells = collect(1:nc_pipe) + grout_cells = collect(1:nc_grout) .+ nc_pipe + # Set up connectivity + pipe_to_pipe = vcat(pipe_cells[1:end-1]', pipe_cells[2:end]') + pipe_to_grout = vcat(pipe_cells', grout_cells') + grout_to_grout = vcat(grout_cells[left_ix]', reverse(grout_cells[right_ix]')) + neighborship = hcat(pipe_to_pipe, pipe_to_grout, grout_to_grout) + # Set up well cell centers and end nodes + wc_supply = cell_centers[:, reservoir_cells[1:nc_mid]] .- [pipe_spacing/2, 0.0, 0.0] + wc_return = cell_centers[:, reservoir_cells[nc_mid+1:end]] .+ [pipe_spacing/2, 0.0, 0.0] + well_cell_centers = repeat(hcat(wc_supply, wc_return), 1, 2) + end_nodes = [pipe_cells[end]] + # Set section for easy lookup + if !ismissing(section) + @warn(["section argument is ignored when neighborship is not provided. ", + "Sections will be created automatically."]) + end + section = Vector{Any}(undef, nc_pipe + nc_grout) + section[pipe_cells[left_ix]] .= [(1, :pipe_left)] + section[pipe_cells[right_ix]] .= [(1, :pipe_right)] + section[grout_cells[left_ix]] .= [(1, :grout_left)] + section[grout_cells[right_ix]] .= [(1, :grout_right)] + else + if ismissing(pipe_cells) || ismissing(grout_cells) || + ismissing(well_cell_centers) || ismissing(end_nodes) + error(["If neighborship is provided, pipe_cells, grout_cells, ", + "well_cell_centers, and end_nodes must also be provided."]) + end + end + # Set return reservoir cell if missing + if ismissing(return_reservoir_cell) + return_reservoir_cell = reservoir_cells[end] + end + + # Set up segment flow models + if ismissing(segment_models) + segment_models = Vector{Any}() + for seg in eachcol(neighborship) + if all([c ∈ pipe_cells for c in seg]) + # Pipe segments use standard wellbore friction model + seg_model = SegmentWellBoreFrictionHB() + else + seg_model = JutulDarcy.ClosedSegment() + end + push!(segment_models, seg_model) + end + end + + # Setup pipe radius + nc = maximum(neighborship) + if radius_pipe isa Number + radius_pipe = fill(radius_pipe, nc) + radius_pipe[grout_cells] .= 0.0 + end + # Setup wall thickness + if wall_thickness_pipe isa Number + wall_thickness_pipe = fill(wall_thickness_pipe, nc) + wall_thickness_pipe[grout_cells] .= 0.0 + end + # Setup pipe wall thermal conductivities + if pipe_thermal_conductivity isa Number + pipe_thermal_conductivity = fill(pipe_thermal_conductivity, nc) + pipe_thermal_conductivity[grout_cells] .= 0.0 + end + # Setup grouting thermal conductivities + if grouting_thermal_conductivity isa Number + grouting_thermal_conductivity = fill(grouting_thermal_conductivity, nc) + grouting_thermal_conductivity[pipe_cells] .= 0.0 + end + + # Set up common well properties + args = ( + type = :closed_loop, + simple_well = false, + WI = 0.0 + ) + # Set up supply well + supply_well = setup_well(D, reservoir_cells; + name = Symbol(name, "_supply"), + neighborship = neighborship, + perforation_cells_well = grout_cells, + well_cell_centers = well_cell_centers, + cell_radius = radius_pipe - wall_thickness_pipe, + casing_thickness = wall_thickness_pipe, + casing_thermal_conductivity = pipe_thermal_conductivity, + grouting_thermal_conductivity = grouting_thermal_conductivity, + segment_models = segment_models, + end_nodes = end_nodes, + args..., kwarg...) + # Augment supply well with closed_loop-specific properties + augment_closed_loop_domain_u1!(supply_well, + radius_grout, + pipe_spacing, + grouting_heat_capacity, + grouting_density; + section = section + ) + # Set default thermal indices + set_default_closed_loop_thermal_indices_u1!( + supply_well, pipe_cells, grout_cells) + # Set up return well + return_well = setup_well(D, return_reservoir_cell; + name = Symbol(name, "_return"), + WIth = 0.0, + args...) + + return [supply_well, return_well] + +end + +function augment_closed_loop_domain_u1!(well::DataDomain, + radius_grout, + pipe_spacing, + grouting_heat_capacity, + grouting_density; + pipe_grout_thermal_index = missing, + grout_grout_thermal_index = missing, + section = missing +) + + c = Cells() + p = Perforations() + f = Faces() + + well[:radius_grout, c] = radius_grout + well[:pipe_spacing, p] = pipe_spacing + well[:material_heat_capacity, c] = grouting_heat_capacity + well[:material_density, c] = grouting_density + + treat_defaulted(x) = x + treat_defaulted(::Missing) = NaN + treat_defaulted(::Nothing) = NaN + + λpg = treat_defaulted(pipe_grout_thermal_index) + λgg = treat_defaulted(grout_grout_thermal_index) + well[:pipe_grout_thermal_index, p] = λpg + well[:grout_grout_thermal_index, p] = λgg + + if !ismissing(section) + well[:section, c] = section + end + +end + +function set_default_closed_loop_thermal_indices_u1!(well::DataDomain, pipe_cells, grout_cells) + + cell = Cells() + face = Faces() + perf = Perforations() + + num_nodes = well.representation.num_nodes + hole_volumes = zeros(Float64, num_nodes) + casing_volumes = zeros(Float64, num_nodes) + grout_volumes = zeros(Float64, num_nodes) + N = well.representation.neighborship + + rep = physical_representation(well) + for (pno, grout_cell) in enumerate(rep.perforations.self) + + @assert grout_cell ∈ grout_cells + segs = findall(N[2, :] .== grout_cell) + @assert length(segs) <= 2 + seg_pg = segs[1] + pipe_cell = N[1, seg_pg] + @assert pipe_cell ∈ pipe_cells + seg_gg = nothing + if length(segs) > 1 + seg_gg = segs[2] + other_grout_cell = N[1, seg_gg] + end + + r_grout = well[:radius_grout, cell][grout_cell] + r_pipe = well[:radius, cell][pipe_cell] + wall_thickness = well[:casing_thickness, cell][pipe_cell] + L = well[:cell_length, cell][pipe_cell] + vol_p, vol_w, vol_g = closed_loop_volume_u1( + L, r_grout, r_pipe + wall_thickness, wall_thickness + ) + + hole_volumes[pipe_cell] = vol_p + hole_volumes[grout_cell] = 0.0 + casing_volumes[pipe_cell] = 0.0 + casing_volumes[grout_cell] = 0.0 + grout_volumes[pipe_cell] = 0.0 + grout_volumes[grout_cell] = vol_g + + pipe_spacing = well[:pipe_spacing, perf][pno] + λg = well[:grouting_thermal_conductivity, cell][grout_cell] + λp = well[:casing_thermal_conductivity, cell][pipe_cell] + + Rpg, Rgr, Rgg = closed_loop_thermal_resistance_u1( + r_grout, r_pipe + wall_thickness, wall_thickness, pipe_spacing, λg, λp) + λpg = L/Rpg + λgr = L/Rgr + λgg = L/Rgg + + if isnan(well[:thermal_well_index, perf][pno]) + well[:thermal_well_index, perf][pno] = λgr + end + + if well[:material_thermal_conductivity, face][seg_pg] == 0.0 + well[:material_thermal_conductivity, face][seg_pg] = λpg + end + + if seg_gg !== nothing && well[:material_thermal_conductivity, face][seg_gg] == 0.0 + well[:material_thermal_conductivity, face][seg_gg] = λgg + end + end + + for pipe_cell in pipe_cells + if hole_volumes[pipe_cell] == 0.0 + # This pipe cell was not assigned a volume, likely because it is not + # connected to a grout cell. Assign a default volume based on geometry. + r_pipe = well[:radius, cell][pipe_cell] + wall_thickness = well[:casing_thickness, cell][pipe_cell] + L = well[:cell_length, cell][pipe_cell] + vol_p, vol_w, vol_g = closed_loop_volume_u1( + L, 0.0, r_pipe + wall_thickness, wall_thickness + ) + hole_volumes[pipe_cell] = vol_p + casing_volumes[pipe_cell] = 0.0 + end + end + + well[:volume_override_hole, cell] = hole_volumes + well[:volume_override_casing, cell] = casing_volumes + well[:volume_override_grouting, cell] = grout_volumes +end + +function closed_loop_volume_u1(length, radius_grout, radius_pipe, wall_thickness) + + # Compute pipe and grout volume + L = length + rg, rp_in, rp_out = radius_grout, radius_pipe - wall_thickness, radius_pipe + vol_hole = π*rp_in^2*L + vol_pipe = π*rp_out^2*L + vol_wall = vol_pipe - vol_hole + vol_grout = π*rg^2*L/2 - vol_pipe + + return vol_hole, vol_wall, vol_grout + +end + +function closed_loop_thermal_resistance_u1( + radius_grout, radius_pipe, wall_thickness_pipe, pipe_spacing, + thermal_conductivity_grout, thermal_conductivity_pipe) + # Conveient short-hand notation + radius_pipe_outer = radius_pipe + radius_pipe_inner = radius_pipe - wall_thickness_pipe + rg, rp_in, rp_out, w = + radius_grout, radius_pipe_inner, radius_pipe_outer, pipe_spacing + λg, λp = thermal_conductivity_grout, thermal_conductivity_pipe + + ## Compute thermal resistances + # Advection-dependent pipe thermal resistance + Ra = 0.0 # TODO: Implement this + # Conduction-dependent pipe thermal resistance + Rc_a = log(rp_out/rp_in)/(2*π*λp) + dg, dp_out = 2*rg, 2*rp_out + x = log(sqrt(dg^2 + 2*dp_out^2)/(2*dp_out))/log(dg/(sqrt(2)*dp_out)) + # Grout thermal resistance + Rg = acosh((dg^2 + dp_out^2 - w^2)/(2*dg*dp_out))/(2*π*λg)*(1.601 - 0.888*w/dg) + # Conduction-dependent grout thermal resistance + Rc_b = x*Rg + # Combined thermal resistance of pipe and grout + Rpg = Ra + Rc_a + Rc_b + Rgr = (1-x)*Rg + # Grout-to-grout thermal resistance + Rar = acosh((2*w^2 - dp_out^2)/dp_out^2)/(2*π*λg) + Rgg = 2*Rgr*(Rar - 2*x*Rg)/(2*Rgr - Rar + 2*x*Rg) + + return Rpg, Rgr, Rgg + +end \ No newline at end of file diff --git a/test/validation.jl b/test/validation.jl index 0c4e9b5..6a8ec7b 100644 --- a/test/validation.jl +++ b/test/validation.jl @@ -27,4 +27,182 @@ using Jutul, JutulDarcy, Fimbul, Test, LinearAlgebra ϵ = compute_error(res, sol, x, t) @test ϵ < 1e-2 +end + +@testset "Validation closed-loop" begin + + meter = si_unit(:meter) + atm = si_unit(:atm) + kilogram = si_unit(:kilogram) + Kelvin, Joule, Watt = si_units(:Kelvin, :joule, :watt) + darcy = si_unit(:darcy) + day = si_unit(:day) + # Operational conditions + T_in = convert_to_si(80.0, :Celsius) # Inlet temperature + T_rock = convert_to_si(10.0, :Celsius) # Ground temperature (assumed constant) + Q = 20.0*meter^3/day # Volumetric flow rate + # Fluid properties + ρf = 988.1*kilogram/meter^3 # Fluid density + Cpf = 4184.0*Joule/(kilogram*Kelvin) # Fluid heat capacity + λf = 0.6405*Watt/(meter*Kelvin) # Fluid thermal conductivity + # Pipe properties + λp = 0.38*Watt/(meter*Kelvin) # Pipe thermal conductivity + # Rock properties + ρr = 2650.0*kilogram/meter^3 # Rock + Cpr = 900.0*Joule/(kilogram*Kelvin) # Rock heat capacity + λr = 2.5*Watt/(meter*Kelvin) # Rock + ϕ = 0.01 # Porosity + K = 1e-3*darcy # Permeability + # BTES geometry + L = 100.0*meter # BTES length + geo_u1 = ( # U1 geometry + closed_loop_type = :u1, + radius_grout = 65e-3*meter, # Grout radius + radius_pipe = 16e-3*meter, # Pipe outer radius + wall_thickness_pipe = 2.9e-3*meter, # Pipe wall thickness + pipe_spacing = 60e-3*meter, # Pipe spacing + pipe_thermal_conductivity = λp, # Pipe thermal conductivity + ) + geo_coax = ( # Coaxial geometry + closed_loop_type = :coaxial, + radius_grout = 50e-3*meter, # Grout radius + radius_pipe_inner = 12e-3*meter, # Inner pipe outer radius + wall_thickness_pipe_inner = 3e-3*meter, # Inner pipe wall thickness + radius_pipe_outer = 25e-3*meter, # Outer pipe outer radius + wall_thickness_pipe_outer = 4e-3*meter, # Outer pipe wall thickness + inner_pipe_thermal_conductivity = λp, # Pipe thermal conductivity + outer_pipe_thermal_conductivity = λp # Pipe thermal conductivity + ); + + function setup_btes_single( # Utility function to set up a BTES simulation case + type; # Type of closed-loop geometry (:u1 or :coaxial) + nz = 110, # Number of cells in vertical direction + n_step = 1, # Number of time steps + inlet = :outer # Inlet pipe for coaxial geometry (:outer or :inner + ) + + ## Select well arguments based on geometry type + if type == :u1 + well_args = geo_u1 + elseif type == :coaxial + well_args = geo_coax + else + error("Unknown closed loop type: $type") + end + ## Create reservoir domain + dims = (1,1,nz) + Δ = (100000.0*dims[3]/L, 100000.0*dims[3]/L, L) + msh = CartesianMesh(dims, Δ) + reservoir = reservoir_domain(msh; + rock_density = ρr, + rock_heat_capacity = Cpr, + rock_thermal_conductivity = λr, + fluid_thermal_conductivity = λf, + component_heat_capacity = Cpf, + porosity = ϕ, + permeability = K + ) + ## Set up BTES well + wells = Fimbul.setup_vertical_btes_well(reservoir, 1, 1; + name = :BTES, + well_args..., + grouting_density = 2000.0, + grouting_heat_capacity = 1500.0, + grouting_thermal_conductivity = 2.3 + ) + ## Set up reservoir model + sys = SinglePhaseSystem(AqueousPhase(), reference_density = 1000.0) + model = setup_reservoir_model(reservoir, sys; wells = wells, thermal = true) + ## Set up controls and boundary conditions + ctrl_inj = InjectorControl( # Injection control with fixed temperature + TotalRateTarget(Q), [1.0], + density=1000.0, temperature=T_in) + ctrl_prod = ProducerControl( # Production control with fixed BHP + BottomHolePressureTarget(5.0*atm)) + geo = tpfv_geometry(msh) + bc_cells = geo.boundary_neighbors + domain = reservoir_model(model).data_domain + bc = flow_boundary_condition( # Fixed pressure and temperature BCs + bc_cells, domain, 5.0*atm, T_rock) + ## Set up forces and initial state + if type == :coaxial && inlet == :outer # Coax with injection in outer pipe + ctrl_supply = ctrl_prod + ctrl_return = ctrl_inj + else # Coax with injection in inner pipe or U1 + ctrl_supply = ctrl_inj + ctrl_return = ctrl_prod + end + controls = Dict( + :BTES_supply => ctrl_supply, + :BTES_return => ctrl_return + ) + forces = setup_reservoir_forces(model; control=controls, bc = bc) + state0 = setup_reservoir_state( # Initialize with uniform pressure and temperature + model; Pressure = 5*atm, Temperature = T_rock) + ## Set time large enough to ensure steady-state + rg = well_args.radius_grout + time = 5/4*2*rg*(ϕ*ρf*Cpf + (1 - ϕ)*ρr*Cpr)/(ϕ*λf + (1 - ϕ)*λr)*10 + dt = fill(time/n_step, n_step) + ## Create Jutul case + case = JutulCase(model, dt, forces; state0 = state0) + ## Set up simulator with temperature-based timestep selector + sim, cfg = setup_reservoir_simulator( + case; initial_dt = 1.0, info_level = -1); + sel = VariableChangeTimestepSelector(:Temperature, 2.5; + relative = false, model = :BTES_supply) + push!(cfg[:timestep_selectors], sel); + + return case, sim, cfg + + end + + function compute_error(case, simulated, analytical) + + well = case.model.models[:BTES_supply].data_domain + sections = last.(well[:section]) |> unique |> collect + sim_temp = simulated.result.states[end][:BTES_supply][:Temperature] + + ϵ = 0.0 + for section in sections + cells = findall(last.(well[:section]) .== section) + z = well[:cell_centroids][3, cells] + dz = well[:cell_length][cells] + L = sum(dz) + Tn = sim_temp[cells] + Ta = analytical[section].(z) + ϵ += sqrt(sum((Tn .- Ta).^2.0.*dz))/L + end + + return ϵ + + end + + tolerance = 1e-1 + + # U1 validation + u1, sim_u1, cfg_u1 = setup_btes_single(:u1); + res_u1 = simulate_reservoir(u1; simulator=sim_u1, config=cfg_u1) + analytical_u1 = Fimbul.analytical_closed_loop_u1(Q, T_in, T_rock, + ρf, Cpf, u1.model.models[:BTES_supply].data_domain) + ϵ_u1 = compute_error(u1, res_u1, analytical_u1) + @test ϵ_u1 < tolerance + + # Coaxial validation with outer inlet + coax_out, sim_coax_out, cfg_coax_out = setup_btes_single(:coaxial, inlet = :outer); + res_coax_out = simulate_reservoir(coax_out; simulator=sim_coax_out, config=cfg_coax_out) + analytical_coax_out = Fimbul.analytical_closed_loop_coaxial(Q, T_in, T_rock, + ρf, Cpf, coax_out.model.models[:BTES_supply].data_domain; + inlet = :outer) + ϵ_coax_out = compute_error(coax_out, res_coax_out, analytical_coax_out) + @test ϵ_coax_out < tolerance + + # Coaxial validation with inner inlet + coax_in, sim_coax_in, cfg_coax_in = setup_btes_single(:coaxial, inlet = :inner); + res_coax_in = simulate_reservoir(coax_in; simulator=sim_coax_in, config=cfg_coax_in) + analytical_coax_in = Fimbul.analytical_closed_loop_coaxial(Q, T_in, T_rock, + ρf, Cpf, coax_in.model.models[:BTES_supply].data_domain; + inlet = :inner) + ϵ_coax_in = compute_error(coax_in, res_coax_in, analytical_coax_in) + @test ϵ_coax_in < tolerance + end \ No newline at end of file