1D high-temperature geothermal benchmark
ValidationThis example reproduces the 1D high-temperature geothermal benchmark cases from [] using the pressure-enthalpy formulation in Fimbul, and validates the results against HYDROTHERM [] reference solutions.
julia
using Jutul, JutulDarcy, Fimbul, GLMakie
to_celsius(T) = convert_from_si.(T, :Celsius)
to_megapascal(p) = convert_from_si.(p, "megapascal")
to_kj_per_kg(h) = h ./ 1e3
const SINGLE_PHASE_CASES = (:a, :b, :c)
const TWO_PHASE_CASES = (:d, :e)
const CASE_COLORS = Dict(
:a => :blue,
:b => :red,
:c => :lightgreen,
:d => :cyan,
:e => :magenta,
)
const PROFILE_SPECS = (
(name = :Pressure, label = "Pressure [MPa]", transform = to_megapascal,
hydrotherm_column = "pressure_mpa"),
(name = :Temperature, label = "Temperature [°C]", transform = to_celsius,
hydrotherm_column = "temperature_c"),
(name = :Saturations, label = "Liquid saturation [-]", transform = x -> vec(x[1,:]),
hydrotherm_column = "liquid_saturation"),
)
nx = 200
cell_size = 10.0;Utilities
julia
function simulate_benchmark_case(case_symbol; vertical = false, nx = 100, cell_size = 10.0)
case = benchmark_ht_1d(
benchmark_case = case_symbol,
nx = nx,
cell_size = cell_size,
vertical = vertical,
)
case = Fimbul.replace_case_timesteps(
case,
Fimbul.load_hydrotherm_1d_timesteps(case_symbol, vertical);
check_sum = true,
)
simulator, config = setup_reservoir_simulator(
case;
tol_cnv = 1e-3,
tol_mb = 1e-7,
max_timestep = Inf,
timesteps = :none,
relaxation = true,
)
results = simulate_reservoir(case; simulator = simulator, config = config)
return (case = case, results = results)
end
function simulate_case_family(case_symbols; nx = 100, cell_size = 10.0)
outputs = Dict{Tuple{Symbol, Bool}, Any}()
for case_symbol in case_symbols
for vertical in Fimbul.available_vertical_modes_ht_1d(case_symbol)
@info "Simulating case $(case_symbol) with vertical = $(vertical)"
outputs[(case_symbol, vertical)] = simulate_benchmark_case(
case_symbol;
vertical = vertical,
nx = nx,
cell_size = cell_size,
)
end
end
return outputs
end
function reservoir_coordinate(out)
domain = out.case.model.models[:Reservoir].data_domain
centroids = tpfv_geometry(physical_representation(domain)).cell_centroids
if out.case.input_data[:vertical]
return vec(centroids[3, :])
else
return vec(centroids[1, :])
end
end
function ordered_values(values, vertical)
values = vec(values)
return vertical ? reverse(values) : values
end
function plot_property_maps(tables)
specs = (
(
variable = :temperature,
label = "Temperature [°C]",
),
(
variable = :density_mix,
label = "Density [kg/m³]",
),
(
variable = :saturation_vapor_ph,
label = "Vapor saturation [-]",
),
)
fig = Figure(size = (1200, 500))
for (i, spec) in enumerate(specs)
yaxisposition = ifelse(i < length(specs), :left, :right)
ax = Axis(fig[2, i];
xticksmirrored = true, yticksmirrored = true,
yaxisposition = yaxisposition, aspect = AxisAspect(1))
handles = Fimbul.plot_phase_diagram_contours!(
ax,
tables;
variable = spec.variable,
pressure_limits = (1e5, 52.5e6),
enthalpy_limits = (500e3, 3500e3),
n_pressure = 500,
n_enthalpy = 500,
levels = 18,
lines = true,
contourf_kwargs = (; colormap = :seaborn_icefire_gradient),
)
Colorbar(fig[1, i], handles.filled, vertical = false, flipaxis = true, width = ax.width,
label = spec.label, labelsize = 20)
if i > 1
hideydecorations!(ax, ticks = false)
end
end
is_ax = [f isa Axis for f in fig.content]
linkaxes!(fig.content[is_ax]...)
return fig
end
function plot_case_profiles(case_symbol, results)
vertical_modes = Fimbul.available_vertical_modes_ht_1d(case_symbol)
fig_height = 400 * length(vertical_modes)
fig = Figure(size = (800 + 400*(case_symbol ∈ TWO_PHASE_CASES), fig_height))
for (k, vertical) in enumerate(vertical_modes)
out = results[(case_symbol, vertical)]
x = reservoir_coordinate(out)
state = out.results.states[end]
row = 2*(k-1)
vertical_label = vertical ? "vertical" : "horizontal"
for (col, spec) in enumerate(PROFILE_SPECS)
if spec.name == :Saturations && case_symbol ∉ TWO_PHASE_CASES
continue
end
values = spec.transform(state[spec.name])
hydrotherm = Fimbul.load_hydrotherm_1d_property(case_symbol, vertical, spec.hydrotherm_column)
if vertical
ax = Axis(
fig[row+1, col];
xlabel = spec.label,
ylabel = "Depth [m]",
yreversed = true,
)
if hydrotherm !== nothing
lines!(ax, hydrotherm.values, hydrotherm.coordinate_m .- cell_size/2;
linewidth = 8, linestyle = :dash, color = :black, label = "HYDROTHERM")
end
lines!(ax, values, x, color = CASE_COLORS[case_symbol], linewidth = 3, label = "Fimbul")
else
ax = Axis(
fig[row+1, col];
xlabel = "Distance [m]",
ylabel = spec.label,
)
y = ordered_values(values, vertical)
if hydrotherm !== nothing
lines!(ax, hydrotherm.coordinate_m .- cell_size/2, hydrotherm.values;
linewidth = 8, linestyle = :dash, color = :black, label = "HYDROTHERM")
end
lines!(ax, x, y, color = CASE_COLORS[case_symbol], linewidth = 3, label = "Fimbul")
if col == 1
axislegend(ax; position = :rt)
end
end
end
fig[row, :] = Label(fig, "Case $(case_symbol) ($(vertical_label))"; fontsize = 20)
end
return fig
end;Simulate all cases
julia
all_results = simulate_case_family(
(SINGLE_PHASE_CASES..., TWO_PHASE_CASES...);
nx = nx, cell_size = cell_size);[ Info: Simulating case a with vertical = false
Jutul: Simulating 250 years, 3.993 weeks as 109 report steps
╭────────────────┬───────────┬───────────────┬──────────╮
│ Iteration type │ Avg/step │ Avg/ministep │ Total │
│ │ 109 steps │ 111 ministeps │ (wasted) │
├────────────────┼───────────┼───────────────┼──────────┤
│ Newton │ 2.02752 │ 1.99099 │ 221 (0) │
│ Linearization │ 3.04587 │ 2.99099 │ 332 (0) │
│ Linear solver │ 2.02752 │ 1.99099 │ 221 (0) │
│ Precond apply │ 4.05505 │ 3.98198 │ 442 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬────────┬────────────┬────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ ms │ Percentage │ s │
├───────────────┼────────┼────────────┼────────┤
│ Properties │ 0.0560 │ 0.64 % │ 0.0124 │
│ Equations │ 2.9159 │ 49.99 % │ 0.9681 │
│ Assembly │ 0.2136 │ 3.66 % │ 0.0709 │
│ Linear solve │ 0.0128 │ 0.15 % │ 0.0028 │
│ Linear setup │ 0.2801 │ 3.20 % │ 0.0619 │
│ Precond apply │ 0.0230 │ 0.53 % │ 0.0102 │
│ Update │ 0.5623 │ 6.42 % │ 0.1243 │
│ Convergence │ 1.6227 │ 27.82 % │ 0.5387 │
│ Input/Output │ 0.1392 │ 0.80 % │ 0.0154 │
│ Other │ 0.5963 │ 6.81 % │ 0.1318 │
├───────────────┼────────┼────────────┼────────┤
│ Total │ 8.7626 │ 100.00 % │ 1.9365 │
╰───────────────┴────────┴────────────┴────────╯
[ Info: Simulating case a with vertical = true
Jutul: Simulating 749 years, 35.04 weeks as 85 report steps
╭────────────────┬──────────┬──────────────┬──────────╮
│ Iteration type │ Avg/step │ Avg/ministep │ Total │
│ │ 85 steps │ 87 ministeps │ (wasted) │
├────────────────┼──────────┼──────────────┼──────────┤
│ Newton │ 2.03529 │ 1.98851 │ 173 (0) │
│ Linearization │ 3.05882 │ 2.98851 │ 260 (0) │
│ Linear solver │ 2.03529 │ 1.98851 │ 173 (0) │
│ Precond apply │ 4.07059 │ 3.97701 │ 346 (0) │
╰────────────────┴──────────┴──────────────┴──────────╯
╭───────────────┬──────────┬────────────┬─────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ μs │ Percentage │ ms │
├───────────────┼──────────┼────────────┼─────────┤
│ Properties │ 43.1580 │ 11.31 % │ 7.4663 │
│ Equations │ 25.0966 │ 9.88 % │ 6.5251 │
│ Assembly │ 7.2081 │ 2.84 % │ 1.8741 │
│ Linear solve │ 11.8793 │ 3.11 % │ 2.0551 │
│ Linear setup │ 121.2265 │ 31.76 % │ 20.9722 │
│ Precond apply │ 22.9675 │ 12.03 % │ 7.9467 │
│ Update │ 10.7252 │ 2.81 % │ 1.8555 │
│ Convergence │ 35.4443 │ 13.96 % │ 9.2155 │
│ Input/Output │ 8.2354 │ 1.09 % │ 0.7165 │
│ Other │ 42.8186 │ 11.22 % │ 7.4076 │
├───────────────┼──────────┼────────────┼─────────┤
│ Total │ 381.7033 │ 100.00 % │ 66.0347 │
╰───────────────┴──────────┴────────────┴─────────╯
[ Info: Simulating case b with vertical = false
Jutul: Simulating 120 years, 1.331 week as 218 report steps
╭────────────────┬───────────┬───────────────┬──────────╮
│ Iteration type │ Avg/step │ Avg/ministep │ Total │
│ │ 218 steps │ 220 ministeps │ (wasted) │
├────────────────┼───────────┼───────────────┼──────────┤
│ Newton │ 2.29817 │ 2.27727 │ 501 (0) │
│ Linearization │ 3.30734 │ 3.27727 │ 721 (0) │
│ Linear solver │ 2.29817 │ 2.27727 │ 501 (0) │
│ Precond apply │ 4.59633 │ 4.55455 │ 1002 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬──────────┬────────────┬──────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ μs │ Percentage │ ms │
├───────────────┼──────────┼────────────┼──────────┤
│ Properties │ 42.1621 │ 10.94 % │ 21.1232 │
│ Equations │ 26.0420 │ 9.72 % │ 18.7763 │
│ Assembly │ 7.3315 │ 2.74 % │ 5.2860 │
│ Linear solve │ 11.8147 │ 3.07 % │ 5.9191 │
│ Linear setup │ 119.2223 │ 30.93 % │ 59.7304 │
│ Precond apply │ 22.8194 │ 11.84 % │ 22.8650 │
│ Update │ 11.7020 │ 3.04 % │ 5.8627 │
│ Convergence │ 41.4214 │ 15.47 % │ 29.8648 │
│ Input/Output │ 10.1234 │ 1.15 % │ 2.2271 │
│ Other │ 42.8140 │ 11.11 % │ 21.4498 │
├───────────────┼──────────┼────────────┼──────────┤
│ Total │ 385.4382 │ 100.00 % │ 193.1045 │
╰───────────────┴──────────┴────────────┴──────────╯
[ Info: Simulating case b with vertical = true
Jutul: Simulating 349 years, 48.81 weeks as 148 report steps
╭────────────────┬───────────┬───────────────┬──────────╮
│ Iteration type │ Avg/step │ Avg/ministep │ Total │
│ │ 148 steps │ 150 ministeps │ (wasted) │
├────────────────┼───────────┼───────────────┼──────────┤
│ Newton │ 2.75 │ 2.71333 │ 407 (0) │
│ Linearization │ 3.76351 │ 3.71333 │ 557 (0) │
│ Linear solver │ 2.75 │ 2.71333 │ 407 (0) │
│ Precond apply │ 5.5 │ 5.42667 │ 814 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬──────────┬────────────┬──────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ μs │ Percentage │ ms │
├───────────────┼──────────┼────────────┼──────────┤
│ Properties │ 43.0446 │ 11.08 % │ 17.5192 │
│ Equations │ 26.8949 │ 9.48 % │ 14.9805 │
│ Assembly │ 7.7790 │ 2.74 % │ 4.3329 │
│ Linear solve │ 12.1466 │ 3.13 % │ 4.9437 │
│ Linear setup │ 123.5493 │ 31.82 % │ 50.2846 │
│ Precond apply │ 24.4134 │ 12.57 % │ 19.8725 │
│ Update │ 11.9264 │ 3.07 % │ 4.8540 │
│ Convergence │ 42.7761 │ 15.08 % │ 23.8263 │
│ Input/Output │ 10.6993 │ 1.02 % │ 1.6049 │
│ Other │ 38.8917 │ 10.02 % │ 15.8289 │
├───────────────┼──────────┼────────────┼──────────┤
│ Total │ 388.3229 │ 100.00 % │ 158.0474 │
╰───────────────┴──────────┴────────────┴──────────╯
[ Info: Simulating case c with vertical = false
Jutul: Simulating 1499 years, 42.08 weeks as 35 report steps
╭────────────────┬──────────┬──────────────┬──────────╮
│ Iteration type │ Avg/step │ Avg/ministep │ Total │
│ │ 35 steps │ 37 ministeps │ (wasted) │
├────────────────┼──────────┼──────────────┼──────────┤
│ Newton │ 2.74286 │ 2.59459 │ 96 (0) │
│ Linearization │ 3.8 │ 3.59459 │ 133 (0) │
│ Linear solver │ 2.74286 │ 2.59459 │ 96 (0) │
│ Precond apply │ 5.48571 │ 5.18919 │ 192 (0) │
╰────────────────┴──────────┴──────────────┴──────────╯
╭───────────────┬──────────┬────────────┬─────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ μs │ Percentage │ ms │
├───────────────┼──────────┼────────────┼─────────┤
│ Properties │ 45.0789 │ 11.29 % │ 4.3276 │
│ Equations │ 26.7189 │ 9.27 % │ 3.5536 │
│ Assembly │ 8.1845 │ 2.84 % │ 1.0885 │
│ Linear solve │ 12.7584 │ 3.20 % │ 1.2248 │
│ Linear setup │ 127.5283 │ 31.94 % │ 12.2427 │
│ Precond apply │ 24.5368 │ 12.29 % │ 4.7111 │
│ Update │ 12.2782 │ 3.07 % │ 1.1787 │
│ Convergence │ 42.6139 │ 14.78 % │ 5.6677 │
│ Input/Output │ 10.5697 │ 1.02 % │ 0.3911 │
│ Other │ 41.1263 │ 10.30 % │ 3.9481 │
├───────────────┼──────────┼────────────┼─────────┤
│ Total │ 399.3113 │ 100.00 % │ 38.3339 │
╰───────────────┴──────────┴────────────┴─────────╯
[ Info: Simulating case c with vertical = true
Jutul: Simulating 1499 years, 42.08 weeks as 35 report steps
╭────────────────┬──────────┬──────────────┬──────────╮
│ Iteration type │ Avg/step │ Avg/ministep │ Total │
│ │ 35 steps │ 37 ministeps │ (wasted) │
├────────────────┼──────────┼──────────────┼──────────┤
│ Newton │ 2.68571 │ 2.54054 │ 94 (0) │
│ Linearization │ 3.74286 │ 3.54054 │ 131 (0) │
│ Linear solver │ 2.68571 │ 2.54054 │ 94 (0) │
│ Precond apply │ 5.37143 │ 5.08108 │ 188 (0) │
╰────────────────┴──────────┴──────────────┴──────────╯
╭───────────────┬──────────┬────────────┬─────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ μs │ Percentage │ ms │
├───────────────┼──────────┼────────────┼─────────┤
│ Properties │ 43.3399 │ 11.04 % │ 4.0739 │
│ Equations │ 26.3044 │ 9.34 % │ 3.4459 │
│ Assembly │ 8.1884 │ 2.91 % │ 1.0727 │
│ Linear solve │ 12.1855 │ 3.10 % │ 1.1454 │
│ Linear setup │ 125.4968 │ 31.96 % │ 11.7967 │
│ Precond apply │ 24.1484 │ 12.30 % │ 4.5399 │
│ Update │ 11.7909 │ 3.00 % │ 1.1083 │
│ Convergence │ 41.9127 │ 14.87 % │ 5.4906 │
│ Input/Output │ 10.0931 │ 1.01 % │ 0.3734 │
│ Other │ 41.1262 │ 10.47 % │ 3.8659 │
├───────────────┼──────────┼────────────┼─────────┤
│ Total │ 392.6887 │ 100.00 % │ 36.9127 │
╰───────────────┴──────────┴────────────┴─────────╯
[ Info: Simulating case d with vertical = false
Jutul: Simulating 200 years, 2.002 days as 3087 report steps
╭────────────────┬────────────┬────────────────┬───────────╮
│ Iteration type │ Avg/step │ Avg/ministep │ Total │
│ │ 3087 steps │ 3089 ministeps │ (wasted) │
├────────────────┼────────────┼────────────────┼───────────┤
│ Newton │ 1.97376 │ 1.97248 │ 6093 (0) │
│ Linearization │ 2.97441 │ 2.97248 │ 9182 (0) │
│ Linear solver │ 1.97376 │ 1.97248 │ 6093 (0) │
│ Precond apply │ 3.94752 │ 3.94497 │ 12186 (0) │
╰────────────────┴────────────┴────────────────┴───────────╯
╭───────────────┬──────────┬────────────┬────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ μs │ Percentage │ s │
├───────────────┼──────────┼────────────┼────────┤
│ Properties │ 42.4053 │ 10.16 % │ 0.2584 │
│ Equations │ 25.4679 │ 9.20 % │ 0.2338 │
│ Assembly │ 7.4095 │ 2.68 % │ 0.0680 │
│ Linear solve │ 11.4642 │ 2.75 % │ 0.0699 │
│ Linear setup │ 119.4378 │ 28.63 % │ 0.7277 │
│ Precond apply │ 23.4078 │ 11.22 % │ 0.2852 │
│ Update │ 11.0872 │ 2.66 % │ 0.0676 │
│ Convergence │ 35.8314 │ 12.94 % │ 0.3290 │
│ Input/Output │ 8.7076 │ 1.06 % │ 0.0269 │
│ Other │ 78.0731 │ 18.71 % │ 0.4757 │
├───────────────┼──────────┼────────────┼────────┤
│ Total │ 417.2404 │ 100.00 % │ 2.5422 │
╰───────────────┴──────────┴────────────┴────────╯
[ Info: Simulating case d with vertical = true
Jutul: Simulating 1000 years, 2.127 weeks as 1623 report steps
╭────────────────┬────────────┬────────────────┬──────────╮
│ Iteration type │ Avg/step │ Avg/ministep │ Total │
│ │ 1623 steps │ 1625 ministeps │ (wasted) │
├────────────────┼────────────┼────────────────┼──────────┤
│ Newton │ 2.07086 │ 2.06831 │ 3361 (0) │
│ Linearization │ 3.07209 │ 3.06831 │ 4986 (0) │
│ Linear solver │ 2.07086 │ 2.06831 │ 3361 (0) │
│ Precond apply │ 4.14171 │ 4.13662 │ 6722 (0) │
╰────────────────┴────────────┴────────────────┴──────────╯
╭───────────────┬──────────┬────────────┬────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ μs │ Percentage │ s │
├───────────────┼──────────┼────────────┼────────┤
│ Properties │ 45.2845 │ 10.12 % │ 0.1522 │
│ Equations │ 26.8900 │ 8.91 % │ 0.1341 │
│ Assembly │ 7.8782 │ 2.61 % │ 0.0393 │
│ Linear solve │ 12.6962 │ 2.84 % │ 0.0427 │
│ Linear setup │ 124.6011 │ 27.84 % │ 0.4188 │
│ Precond apply │ 23.8485 │ 10.66 % │ 0.1603 │
│ Update │ 11.5685 │ 2.59 % │ 0.0389 │
│ Convergence │ 71.2241 │ 23.61 % │ 0.3551 │
│ Input/Output │ 9.2052 │ 0.99 % │ 0.0150 │
│ Other │ 43.9589 │ 9.82 % │ 0.1477 │
├───────────────┼──────────┼────────────┼────────┤
│ Total │ 447.4951 │ 100.00 % │ 1.5040 │
╰───────────────┴──────────┴────────────┴────────╯
[ Info: Simulating case e with vertical = false
Jutul: Simulating 1999 years, 51.39 weeks as 4495 report steps
╭────────────────┬────────────┬────────────────┬────────────╮
│ Iteration type │ Avg/step │ Avg/ministep │ Total │
│ │ 4495 steps │ 4501 ministeps │ (wasted) │
├────────────────┼────────────┼────────────────┼────────────┤
│ Newton │ 2.25184 │ 2.24883 │ 10122 (30) │
│ Linearization │ 3.25317 │ 3.24883 │ 14623 (32) │
│ Linear solver │ 2.25184 │ 2.24883 │ 10122 (30) │
│ Precond apply │ 4.50367 │ 4.49767 │ 20244 (60) │
╰────────────────┴────────────┴────────────────┴────────────╯
╭───────────────┬──────────┬────────────┬────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ μs │ Percentage │ s │
├───────────────┼──────────┼────────────┼────────┤
│ Properties │ 44.6031 │ 10.21 % │ 0.4515 │
│ Equations │ 25.4196 │ 8.40 % │ 0.3717 │
│ Assembly │ 7.3454 │ 2.43 % │ 0.1074 │
│ Linear solve │ 12.1071 │ 2.77 % │ 0.1225 │
│ Linear setup │ 146.5793 │ 33.54 % │ 1.4837 │
│ Precond apply │ 23.0783 │ 10.56 % │ 0.4672 │
│ Update │ 27.7725 │ 6.36 % │ 0.2811 │
│ Convergence │ 46.9419 │ 15.52 % │ 0.6864 │
│ Input/Output │ 8.7196 │ 0.89 % │ 0.0392 │
│ Other │ 40.7284 │ 9.32 % │ 0.4123 │
├───────────────┼──────────┼────────────┼────────┤
│ Total │ 436.9749 │ 100.00 % │ 4.4231 │
╰───────────────┴──────────┴────────────┴────────╯H2O properties in pressure-enthalpy space
The steam tables have been generated using the CoolProp library []–see also FimbulCoolPropExt``). They are available in the model's fluid system, but can also be loaded directly using Artifacts. To generate your own steam tables, see thebuild_steam_tables_h2ofunction inFimbulCoolPropExt`.
We first inspect the steam tables directly. The figure below shows three key properties in
julia
tables = Fimbul.steam_tables_h2o()
fig_properties = plot_property_maps(tables)
fig_properties
Single-phase benchmark cases
Cases :a to :c remain in the single-phase region for the plotted states. We therefore use them to compare Fimbul and HYDROTHERM pressure and temperature profiles along the 1D column, both with and without gravity.
Case a
Case :a spans 50 to 25 MPa and 350 to 150 °C. The pressure stays high enough that the fluid remains in the compressed-liquid region throughout the column, so the solution is single-phase liquid with smooth pressure and temperature variations.
julia
fig_case_a = plot_case_profiles(:a, all_results)
fig_case_a
Case b
Case :b spans 40 to 20 MPa and 450 to 300 °C. These conditions are hotter than case :a, but the pressure is still high enough to avoid flashing, so the response remains single-phase while showing stronger thermal contrasts.
julia
fig_case_b = plot_case_profiles(:b, all_results)
fig_case_b
Case c
Case :c spans 15 to 1 MPa and 500 to 350 °C. Here the fluid is much hotter and less compressed, placing it in the vapor-dominated single-phase region, so the profiles represent hot steam rather than liquid water.
julia
fig_case_c = plot_case_profiles(:c, all_results)
fig_case_c
Two-phase benchmark cases
Cases :d and :e traverse the two-phase region. Here we compare Fimbul and HYDROTHERM pressure, temperature, and liquid-saturation profiles, and then compare the no-gravity paths in pressure-enthalpy space.
Case d
Case :d spans 20 to 1 MPa and 400 to 150 °C. The inlet starts as hot, pressurized water, but the strong pressure and temperature drop drives the state path across the saturation envelope, producing a genuine two-phase liquid-vapor transition along the column.
julia
fig_case_d = plot_case_profiles(:d, all_results)
fig_case_d
Case e
Case :e spans 4 to 1 MPa and 300 to 150 °C. Because the pressures are low, boiling is easier to trigger than in case :d, so the system develops a broad two-phase region with more pronounced vapor formation.
julia
fig_case_e = plot_case_profiles(:e, all_results)
fig_case_e
Phase diagram comparison
We compare the horizontal cases directly in pressure-enthalpy space, overlaying Fimbul state paths and HYDROTHERM reference paths on temperature contours.
julia
fig = Figure(size = (700, 640))
Label(fig[0, 1:2], "Phase diagram comparison"; fontsize = 22)
ax = Axis(fig[1, 1]; xticksmirrored = true, yticksmirrored = true, aspect = AxisAspect(1))
handles = Fimbul.plot_phase_diagram_contours!(
ax,
tables;
variable = :temperature,
pressure_limits = (1e5, 52.5e6),
enthalpy_limits = (500e3, 3500e3),
levels = 20,
lines = true,
)
for case_symbol in (SINGLE_PHASE_CASES..., TWO_PHASE_CASES...)
hydrotherm = Fimbul.load_hydrotherm_1d_phase_path(case_symbol)
lines!(ax, hydrotherm.enthalpy_kj_per_kg, hydrotherm.pressure_mpa; linewidth = 6, linestyle = :dash, color = CASE_COLORS[case_symbol])
out = all_results[(case_symbol, false)]
state = out.results.states[end]
Fimbul.plot_reservoir_state_ph!(
ax,
state;
color = CASE_COLORS[case_symbol],
linewidth = 3,
label = "Case $(case_symbol)",
)
end
axislegend(ax; position = :rt)
Colorbar(fig[1, 2], handles.filled, label = "Temperature [°C]")
fig
Example on GitHub
If you would like to run this example yourself, it can be downloaded from the Fimbul.jl GitHub repository as a script.
This example took 21.122348788 seconds to complete.This page was generated using Literate.jl.