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.01835 │ 1.98198 │ 220 (0) │
│ Linearization │ 3.0367 │ 2.98198 │ 331 (0) │
│ Linear solver │ 2.01835 │ 1.98198 │ 220 (0) │
│ Precond apply │ 4.0367 │ 3.96396 │ 440 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬─────────┬────────────┬────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ ms │ Percentage │ s │
├───────────────┼─────────┼────────────┼────────┤
│ Properties │ 0.1079 │ 0.48 % │ 0.0237 │
│ Equations │ 3.9437 │ 26.20 % │ 1.3054 │
│ Assembly │ 0.7716 │ 5.12 % │ 0.2554 │
│ Linear solve │ 0.6380 │ 2.82 % │ 0.1404 │
│ Linear setup │ 0.9331 │ 4.12 % │ 0.2053 │
│ Precond apply │ 0.0393 │ 0.35 % │ 0.0173 │
│ Update │ 2.2546 │ 9.95 % │ 0.4960 │
│ Convergence │ 2.9758 │ 19.77 % │ 0.9850 │
│ Input/Output │ 0.5610 │ 1.25 % │ 0.0623 │
│ Other │ 6.7841 │ 29.95 % │ 1.4925 │
├───────────────┼─────────┼────────────┼────────┤
│ Total │ 22.6510 │ 100.00 % │ 4.9832 │
╰───────────────┴─────────┴────────────┴────────╯
[ 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.01176 │ 1.96552 │ 171 (0) │
│ Linearization │ 3.03529 │ 2.96552 │ 258 (0) │
│ Linear solver │ 2.01176 │ 1.96552 │ 171 (0) │
│ Precond apply │ 4.02353 │ 3.93103 │ 342 (0) │
╰────────────────┴──────────┴──────────────┴──────────╯
╭───────────────┬──────────┬────────────┬──────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ μs │ Percentage │ ms │
├───────────────┼──────────┼────────────┼──────────┤
│ Properties │ 91.8709 │ 13.36 % │ 15.7099 │
│ Equations │ 56.4772 │ 12.39 % │ 14.5711 │
│ Assembly │ 13.9095 │ 3.05 % │ 3.5886 │
│ Linear solve │ 31.7454 │ 4.62 % │ 5.4285 │
│ Linear setup │ 246.3743 │ 35.82 % │ 42.1300 │
│ Precond apply │ 39.2246 │ 11.41 % │ 13.4148 │
│ Update │ 15.5376 │ 2.26 % │ 2.6569 │
│ Convergence │ 36.2319 │ 7.95 % │ 9.3478 │
│ Input/Output │ 22.1160 │ 1.64 % │ 1.9241 │
│ Other │ 51.6937 │ 7.52 % │ 8.8396 │
├───────────────┼──────────┼────────────┼──────────┤
│ Total │ 687.7861 │ 100.00 % │ 117.6114 │
╰───────────────┴──────────┴────────────┴──────────╯
[ 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.29358 │ 2.27273 │ 500 (0) │
│ Linearization │ 3.30275 │ 3.27273 │ 720 (0) │
│ Linear solver │ 2.29358 │ 2.27273 │ 500 (0) │
│ Precond apply │ 4.58716 │ 4.54545 │ 1000 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬──────────┬────────────┬──────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ μs │ Percentage │ ms │
├───────────────┼──────────┼────────────┼──────────┤
│ Properties │ 85.0163 │ 9.56 % │ 42.5082 │
│ Equations │ 54.2583 │ 8.78 % │ 39.0660 │
│ Assembly │ 184.4402 │ 29.86 % │ 132.7969 │
│ Linear solve │ 30.8842 │ 3.47 % │ 15.4421 │
│ Linear setup │ 245.3428 │ 27.58 % │ 122.6714 │
│ Precond apply │ 38.1896 │ 8.59 % │ 38.1896 │
│ Update │ 11.9179 │ 1.34 % │ 5.9590 │
│ Convergence │ 31.8347 │ 5.15 % │ 22.9210 │
│ Input/Output │ 16.3494 │ 0.81 % │ 3.5969 │
│ Other │ 43.3074 │ 4.87 % │ 21.6537 │
├───────────────┼──────────┼────────────┼──────────┤
│ Total │ 889.6093 │ 100.00 % │ 444.8047 │
╰───────────────┴──────────┴────────────┴──────────╯
[ 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.45946 │ 2.42667 │ 364 (0) │
│ Linearization │ 3.47297 │ 3.42667 │ 514 (0) │
│ Linear solver │ 2.45946 │ 2.42667 │ 364 (0) │
│ Precond apply │ 4.91892 │ 4.85333 │ 728 (0) │
╰────────────────┴───────────┴───────────────┴──────────╯
╭───────────────┬──────────┬────────────┬──────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ μs │ Percentage │ ms │
├───────────────┼──────────┼────────────┼──────────┤
│ Properties │ 87.3546 │ 13.73 % │ 31.7971 │
│ Equations │ 54.1595 │ 12.02 % │ 27.8380 │
│ Assembly │ 13.0481 │ 2.90 % │ 6.7067 │
│ Linear solve │ 30.5963 │ 4.81 % │ 11.1371 │
│ Linear setup │ 238.0665 │ 37.41 % │ 86.6562 │
│ Precond apply │ 38.6790 │ 12.16 % │ 28.1583 │
│ Update │ 11.9514 │ 1.88 % │ 4.3503 │
│ Convergence │ 32.9937 │ 7.32 % │ 16.9588 │
│ Input/Output │ 16.9006 │ 1.09 % │ 2.5351 │
│ Other │ 42.5019 │ 6.68 % │ 15.4707 │
├───────────────┼──────────┼────────────┼──────────┤
│ Total │ 636.2864 │ 100.00 % │ 231.6083 │
╰───────────────┴──────────┴────────────┴──────────╯
[ 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 │ 85.9678 │ 13.39 % │ 8.2529 │
│ Equations │ 55.1410 │ 11.90 % │ 7.3338 │
│ Assembly │ 13.2376 │ 2.86 % │ 1.7606 │
│ Linear solve │ 31.4250 │ 4.90 % │ 3.0168 │
│ Linear setup │ 244.3666 │ 38.07 % │ 23.4592 │
│ Precond apply │ 38.7342 │ 12.07 % │ 7.4370 │
│ Update │ 11.9990 │ 1.87 % │ 1.1519 │
│ Convergence │ 33.8983 │ 7.32 % │ 4.5085 │
│ Input/Output │ 17.1036 │ 1.03 % │ 0.6328 │
│ Other │ 42.4027 │ 6.61 % │ 4.0707 │
├───────────────┼──────────┼────────────┼─────────┤
│ Total │ 641.9178 │ 100.00 % │ 61.6241 │
╰───────────────┴──────────┴────────────┴─────────╯
[ 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 │ 85.9666 │ 13.44 % │ 8.0809 │
│ Equations │ 54.4305 │ 11.86 % │ 7.1304 │
│ Assembly │ 13.2943 │ 2.90 % │ 1.7415 │
│ Linear solve │ 31.1010 │ 4.86 % │ 2.9235 │
│ Linear setup │ 242.6521 │ 37.92 % │ 22.8093 │
│ Precond apply │ 39.0524 │ 12.21 % │ 7.3418 │
│ Update │ 12.2606 │ 1.92 % │ 1.1525 │
│ Convergence │ 33.5652 │ 7.31 % │ 4.3970 │
│ Input/Output │ 16.3181 │ 1.00 % │ 0.6038 │
│ Other │ 42.1716 │ 6.59 % │ 3.9641 │
├───────────────┼──────────┼────────────┼─────────┤
│ Total │ 639.8391 │ 100.00 % │ 60.1449 │
╰───────────────┴──────────┴────────────┴─────────╯
[ 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.97344 │ 1.97216 │ 6092 (0) │
│ Linearization │ 2.97408 │ 2.97216 │ 9181 (0) │
│ Linear solver │ 1.97344 │ 1.97216 │ 6092 (0) │
│ Precond apply │ 3.94687 │ 3.94432 │ 12184 (0) │
╰────────────────┴────────────┴────────────────┴───────────╯
╭───────────────┬──────────┬────────────┬────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ μs │ Percentage │ s │
├───────────────┼──────────┼────────────┼────────┤
│ Properties │ 106.3961 │ 15.17 % │ 0.6482 │
│ Equations │ 55.6798 │ 11.97 % │ 0.5112 │
│ Assembly │ 13.6502 │ 2.93 % │ 0.1253 │
│ Linear solve │ 49.8776 │ 7.11 % │ 0.3039 │
│ Linear setup │ 240.7145 │ 34.33 % │ 1.4664 │
│ Precond apply │ 38.9140 │ 11.10 % │ 0.4741 │
│ Update │ 13.1582 │ 1.88 % │ 0.0802 │
│ Convergence │ 33.8899 │ 7.28 % │ 0.3111 │
│ Input/Output │ 18.4035 │ 1.33 % │ 0.0568 │
│ Other │ 48.2768 │ 6.89 % │ 0.2941 │
├───────────────┼──────────┼────────────┼────────┤
│ Total │ 701.1412 │ 100.00 % │ 4.2714 │
╰───────────────┴──────────┴────────────┴────────╯
[ 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.07702 │ 2.07446 │ 3371 (0) │
│ Linearization │ 3.07825 │ 3.07446 │ 4996 (0) │
│ Linear solver │ 2.07702 │ 2.07446 │ 3371 (0) │
│ Precond apply │ 4.15404 │ 4.14892 │ 6742 (0) │
╰────────────────┴────────────┴────────────────┴──────────╯
╭───────────────┬──────────┬────────────┬────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ μs │ Percentage │ s │
├───────────────┼──────────┼────────────┼────────┤
│ Properties │ 92.3898 │ 13.90 % │ 0.3114 │
│ Equations │ 55.6053 │ 12.40 % │ 0.2778 │
│ Assembly │ 14.2523 │ 3.18 % │ 0.0712 │
│ Linear solve │ 31.2636 │ 4.71 % │ 0.1054 │
│ Linear setup │ 240.4617 │ 36.19 % │ 0.8106 │
│ Precond apply │ 39.0510 │ 11.75 % │ 0.2633 │
│ Update │ 12.9090 │ 1.94 % │ 0.0435 │
│ Convergence │ 33.8224 │ 7.54 % │ 0.1690 │
│ Input/Output │ 18.0132 │ 1.31 % │ 0.0293 │
│ Other │ 46.9860 │ 7.07 % │ 0.1584 │
├───────────────┼──────────┼────────────┼────────┤
│ Total │ 664.4547 │ 100.00 % │ 2.2399 │
╰───────────────┴──────────┴────────────┴────────╯
[ 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.25161 │ 2.24861 │ 10121 (30) │
│ Linearization │ 3.25295 │ 3.24861 │ 14622 (32) │
│ Linear solver │ 2.25161 │ 2.24861 │ 10121 (30) │
│ Precond apply │ 4.50323 │ 4.49722 │ 20242 (60) │
╰────────────────┴────────────┴────────────────┴────────────╯
╭───────────────┬──────────┬────────────┬────────╮
│ Timing type │ Each │ Relative │ Total │
│ │ μs │ Percentage │ s │
├───────────────┼──────────┼────────────┼────────┤
│ Properties │ 94.0778 │ 11.68 % │ 0.9522 │
│ Equations │ 56.2949 │ 10.10 % │ 0.8231 │
│ Assembly │ 13.6921 │ 2.46 % │ 0.2002 │
│ Linear solve │ 31.6006 │ 3.92 % │ 0.3198 │
│ Linear setup │ 322.4148 │ 40.04 % │ 3.2632 │
│ Precond apply │ 38.0789 │ 9.46 % │ 0.7708 │
│ Update │ 13.4604 │ 1.67 % │ 0.1362 │
│ Convergence │ 59.3059 │ 10.64 % │ 0.8672 │
│ Input/Output │ 59.0880 │ 3.26 % │ 0.2660 │
│ Other │ 54.3817 │ 6.75 % │ 0.5504 │
├───────────────┼──────────┼────────────┼────────┤
│ Total │ 805.1626 │ 100.00 % │ 8.1491 │
╰───────────────┴──────────┴────────────┴────────╯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 55.413959634 seconds to complete.This page was generated using Literate.jl.