diff --git a/Project.toml b/Project.toml index c8d0a66..9a98f41 100644 --- a/Project.toml +++ b/Project.toml @@ -1,25 +1,27 @@ name = "MultiComponentFlash" uuid = "35e5bd01-9722-4017-9deb-64a5d32478ff" +version = "1.1.19" authors = ["Olav Møyner "] -version = "1.1.18" [deps] ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] ForwardDiff = "0.10, 1" +LinearAlgebra = "1" +Printf = "1" Roots = "1, 2" StaticArrays = "1" -LinearAlgebra = "1" julia = "1.6" [extras] -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] test = ["Test", "ForwardDiff", "StaticArrays"] diff --git a/docs/src/examples/advanced.md b/docs/src/examples/advanced.md index b0519e1..7b225cc 100644 --- a/docs/src/examples/advanced.md +++ b/docs/src/examples/advanced.md @@ -139,3 +139,11 @@ xlabel!("T [°Celsius]") ``` ![Phase diagram](../assets/phase_diagram_simple.png) + +### PVT table generation + +There is experimental support for generating simulator input blackoil tables (e.g. PVTG/PVDG and PVTO/PVDO). + +```@docs +generate_pvt_tables +``` diff --git a/src/MultiComponentFlash.jl b/src/MultiComponentFlash.jl index 0c9f821..6514b7e 100644 --- a/src/MultiComponentFlash.jl +++ b/src/MultiComponentFlash.jl @@ -12,7 +12,7 @@ module MultiComponentFlash include("flash_types.jl") include("flow_coupler_types.jl") # Types that define specific cubic equations of state - export SoaveRedlichKwong, RedlichKwong, PengRobinson, PengRobinsonCorrected, ZudkevitchJoffe + export SoaveRedlichKwong, RedlichKwong, PengRobinson, PengRobinsonCorrected, ZudkevitchJoffe, SoreideWhitson # The generic cubic form that supports the above export GenericCubicEOS # KValue "fake" EOS @@ -57,4 +57,8 @@ module MultiComponentFlash include("flow_coupler.jl") include("utils.jl") + + include("PVTExperiments/PVTExperiments.jl") + using .PVTExperiments + export generate_pvt_tables end diff --git a/src/PVTExperiments/PVTExperiments.jl b/src/PVTExperiments/PVTExperiments.jl new file mode 100644 index 0000000..83761bc --- /dev/null +++ b/src/PVTExperiments/PVTExperiments.jl @@ -0,0 +1,25 @@ +module PVTExperiments + # Notes: This module was in part generated with the help of Claude Opus. We + # export a single function, `generate_pvt_tables`, which is the main + # interface for generating all the tables at once. Additional experiments + # are available, but their interfaces are not exported and can be changed + # without breaking versions. + using MultiComponentFlash + using Printf + + import MultiComponentFlash: flashed_mixture_2ph, number_of_components + + const WATER_DENSITY_SC = 999.0 # kg/m³ at standard conditions + const R_GAS = MultiComponentFlash.IDEAL_GAS_CONSTANT # 8.3144598 J/(mol·K) + + include("types.jl") + include("utils.jl") + include("cce.jl") + include("dle.jl") + include("cvd.jl") + include("mss.jl") + include("tables.jl") + include("interface.jl") + + export generate_pvt_tables +end diff --git a/src/PVTExperiments/cce.jl b/src/PVTExperiments/cce.jl new file mode 100644 index 0000000..1a0680a --- /dev/null +++ b/src/PVTExperiments/cce.jl @@ -0,0 +1,63 @@ +""" + cce(eos, z, T; p_range, n_points) + +Perform a Constant Composition Expansion (CCE) experiment. + +The mixture with composition `z` is held at temperature `T` (K) and the +pressure is decreased from a high value through the saturation pressure. +The volume is measured relative to the volume at the saturation point. + +# Arguments +- `eos`: Equation of state +- `z`: Overall mole fractions +- `T`: Temperature (K) +- `p_range`: Tuple (p_max, p_min) in Pa. Default: (50e6, 1e5) +- `n_points`: Number of pressure steps. Default: 40 +""" +function cce(eos, z, T; + p_range = (50e6, 1e5), + n_points = 40 + ) + z = collect(Float64, z) + z ./= sum(z) + nc = number_of_components(eos) + + # Find saturation pressure + p_sat, is_bp = find_saturation_pressure(eos, z, T; + p_min = p_range[2], p_max = p_range[1]) + + # Generate pressure steps + pressures = collect(range(p_range[1], p_range[2], length = n_points)) + sort!(pressures, rev = true) + + # Compute properties at saturation point for reference volume + props_sat = flash_and_properties(eos, p_sat, T, z) + V_mol_sat = (1.0 - props_sat.V) * props_sat.V_mol_l + props_sat.V * props_sat.V_mol_v + + # Storage + Z_factors = zeros(n_points) + V_factors = zeros(n_points) + rel_vol = zeros(n_points) + ρ_oil = zeros(n_points) + ρ_gas = zeros(n_points) + μ_oil = zeros(n_points) + μ_gas = zeros(n_points) + + for (i, p) in enumerate(pressures) + props = flash_and_properties(eos, p, T, z) + V_factors[i] = props.V + Z_factors[i] = (1.0 - props.V) * props.Z_L + props.V * props.Z_V + + # Total molar volume at this pressure + V_mol = (1.0 - props.V) * props.V_mol_l + props.V * props.V_mol_v + rel_vol[i] = V_mol / V_mol_sat + + ρ_oil[i] = props.ρ_l + ρ_gas[i] = props.ρ_v + μ_oil[i] = props.μ_l + μ_gas[i] = props.μ_v + end + + return CCEResult(T, z, pressures, Z_factors, V_factors, rel_vol, + ρ_oil, ρ_gas, μ_oil, μ_gas, p_sat, is_bp) +end diff --git a/src/PVTExperiments/cvd.jl b/src/PVTExperiments/cvd.jl new file mode 100644 index 0000000..52ce854 --- /dev/null +++ b/src/PVTExperiments/cvd.jl @@ -0,0 +1,125 @@ +""" + cvd(eos, z, T; p_range, n_points, p_sc, T_sc) + +Perform a Constant Volume Depletion (CVD) experiment. + +Starting at the dew point, the pressure is reduced in steps. At each step, +gas is removed to restore the cell to its original volume. The liquid dropout +and produced gas properties are recorded. + +This experiment is typical for gas condensate systems. + +# Arguments +- `eos`: Equation of state +- `z`: Overall mole fractions +- `T`: Temperature (K) +- `p_range`: (p_max, p_min) for the experiment. Default: (50e6, 1e5) +- `n_points`: Number of pressure steps. Default: 20 +- `p_sc`: Standard condition pressure (Pa). Default: 101325.0 +- `T_sc`: Standard condition temperature (K). Default: 288.706 (60°F) +""" +function cvd(eos, z, T; + p_range = (50e6, 1e5), + n_points = 20, + p_sc = 101325.0, + T_sc = 288.706 + ) + z = collect(Float64, z) + z ./= sum(z) + nc = number_of_components(eos) + + # Find saturation pressure + p_sat, is_bp = find_saturation_pressure(eos, z, T; + p_min = p_range[2], p_max = p_range[1]) + + if is_bp + @warn "CVD is designed for dew point (gas condensate) fluids. Found bubble point fluid." + end + + # Pressure steps from p_sat down + pressures = collect(range(p_sat, p_range[2], length = n_points + 1)) + + # Storage + Z_factors = zeros(n_points + 1) + liq_dropout = zeros(n_points + 1) + gas_produced_cum = zeros(n_points + 1) + ρ_gas_arr = zeros(n_points + 1) + μ_gas_arr = zeros(n_points + 1) + Z_gas_arr = zeros(n_points + 1) + gas_comp = Vector{Vector{Float64}}(undef, n_points + 1) + + # Initial state at dew point: single phase gas + props_init = flash_and_properties(eos, p_sat, T, z) + V_cell = props_init.V_mol_v # Reference cell volume per mole + + n_total = 1.0 # Total moles in cell + z_cell = copy(z) # Current overall composition in cell + cum_gas_produced = 0.0 # Cumulative gas produced (moles) + + for (i, p) in enumerate(pressures) + if i == 1 + # At dew point - single phase gas + Z_factors[i] = props_init.Z_V + liq_dropout[i] = 0.0 + gas_produced_cum[i] = 0.0 + ρ_gas_arr[i] = props_init.ρ_v + μ_gas_arr[i] = props_init.μ_v + Z_gas_arr[i] = props_init.Z_V + gas_comp[i] = copy(z) + else + # Flash at lower pressure + props = flash_and_properties(eos, p, T, z_cell) + + Z_gas_arr[i] = props.Z_V + gas_comp[i] = copy(props.y) + ρ_gas_arr[i] = props.ρ_v + μ_gas_arr[i] = props.μ_v + + # Two-phase Z factor + Z_factors[i] = (1.0 - props.V) * props.Z_L + props.V * props.Z_V + + # Volumes + V_liq = n_total * (1.0 - props.V) * props.V_mol_l + V_gas = n_total * props.V * props.V_mol_v + V_total = V_liq + V_gas + + # Liquid dropout = liquid volume / cell volume + liq_dropout[i] = V_liq / V_cell + + # Gas to remove to restore cell volume + V_excess = V_total - V_cell + if V_excess > 0 && props.V_mol_v > 0 + n_gas_remove = V_excess / props.V_mol_v + else + n_gas_remove = 0.0 + end + + cum_gas_produced += n_gas_remove + gas_produced_cum[i] = cum_gas_produced + + # Update cell: remove gas, keep liquid + remaining gas + n_gas_remaining = n_total * props.V - n_gas_remove + n_liq = n_total * (1.0 - props.V) + + if n_gas_remaining < 0 + n_gas_remaining = 0.0 + end + + n_total_new = n_liq + n_gas_remaining + + # New composition + if n_total_new > 0 + z_cell = (n_liq .* props.x .+ n_gas_remaining .* props.y) ./ n_total_new + z_cell ./= sum(z_cell) + end + n_total = n_total_new + end + end + + # Normalize cumulative gas produced + gas_produced_cum ./= max(gas_produced_cum[end], 1e-30) + + return CVDResult(T, collect(Float64, z ./ sum(z)), pressures, + Z_factors, liq_dropout, gas_produced_cum, + ρ_gas_arr, μ_gas_arr, Z_gas_arr, gas_comp, p_sat) +end diff --git a/src/PVTExperiments/dle.jl b/src/PVTExperiments/dle.jl new file mode 100644 index 0000000..6cf1c8b --- /dev/null +++ b/src/PVTExperiments/dle.jl @@ -0,0 +1,154 @@ +""" + dle(eos, z, T; p_range, n_points, p_sc, T_sc) + +Perform a Differential Liberation Expansion (DLE) experiment. + +Starting at the bubble point, the pressure is reduced in steps. At each step, +the equilibrium gas is completely removed and the liquid is re-equilibrated at +the next lower pressure. Properties (Bo, Rs, etc.) are referenced to the +final stock-tank oil volume at standard conditions. + +# Arguments +- `eos`: Equation of state +- `z`: Overall mole fractions of initial fluid +- `T`: Reservoir temperature (K) +- `p_range`: (p_max, p_min) for the experiment. Default: (50e6, 1e5) +- `n_points`: Number of pressure steps below the saturation point. Default: 20 +- `p_sc`: Standard condition pressure (Pa). Default: 101325.0 +- `T_sc`: Standard condition temperature (K). Default: 288.706 (60°F) +""" +function dle(eos, z, T; + p_range = (50e6, 1e5), + n_points = 20, + p_sc = 101325.0, + T_sc = 288.706 + ) + z = collect(Float64, z) + z ./= sum(z) + nc = number_of_components(eos) + + # Find saturation pressure + p_sat, is_bp = find_saturation_pressure(eos, z, T; + p_min = p_range[2], p_max = p_range[1]) + + if !is_bp + @warn "DLE is designed for bubble point fluids. Found dew point fluid." + end + + # Generate pressure steps from p_sat down to p_min + pressures = collect(range(p_sat, p_range[2], length = n_points + 1)) + + # Storage + Bo_arr = zeros(n_points + 1) + Rs_arr = zeros(n_points + 1) + ρ_oil_arr = zeros(n_points + 1) + ρ_gas_arr = zeros(n_points + 1) + μ_oil_arr = zeros(n_points + 1) + μ_gas_arr = zeros(n_points + 1) + Bg_arr = zeros(n_points + 1) + Z_gas_arr = zeros(n_points + 1) + gas_comp = Vector{Vector{Float64}}(undef, n_points + 1) + oil_comp = Vector{Vector{Float64}}(undef, n_points + 1) + + # Track remaining oil and gas + z_oil = copy(z) # Current oil composition + n_oil = 1.0 # Moles of oil remaining (normalized) + + # Flash the residual oil to standard conditions to get reference volume + # First, step through the liberation + moles_gas_removed = Float64[] # Gas removed at each step + oil_molar_volumes = Float64[] # Oil molar volume at reservoir conditions + gas_molar_volumes_sc = Float64[] # Gas molar volume at standard conditions + + for (i, p) in enumerate(pressures) + if i == 1 + # At bubble point - single phase liquid + props = flash_and_properties(eos, p, T, z_oil) + gas_comp[i] = copy(z_oil) + oil_comp[i] = copy(z_oil) + ρ_oil_arr[i] = props.ρ_l + ρ_gas_arr[i] = props.ρ_v + μ_oil_arr[i] = props.μ_l + μ_gas_arr[i] = props.μ_v + Z_gas_arr[i] = props.Z_V + Bg_arr[i] = 0.0 + push!(oil_molar_volumes, props.V_mol_l) + push!(moles_gas_removed, 0.0) + else + # Flash at this pressure + props = flash_and_properties(eos, p, T, z_oil) + + gas_comp[i] = copy(props.y) + oil_comp[i] = copy(props.x) + ρ_oil_arr[i] = props.ρ_l + ρ_gas_arr[i] = props.ρ_v + μ_oil_arr[i] = props.μ_l + μ_gas_arr[i] = props.μ_v + Z_gas_arr[i] = props.Z_V + Bg_arr[i] = phase_molar_volume_from_Z(p, T, props.Z_V) / + phase_molar_volume_from_Z(p_sc, T_sc, 1.0) + push!(oil_molar_volumes, props.V_mol_l) + + # Gas removed at this step + n_gas = n_oil * props.V + push!(moles_gas_removed, n_gas) + + # Update oil: only liquid remains + n_oil = n_oil * (1.0 - props.V) + z_oil = copy(props.x) + end + end + + # Get stock tank oil properties: flash remaining oil at standard conditions + props_sc = flash_and_properties(eos, p_sc, T_sc, z_oil) + V_mol_oil_sc = props_sc.V_mol_l + mw_oil_sc = mixture_molar_mass(props_sc.x, eos) + + # Reference: 1 mol of original fluid at bubble point + # Volume of stock tank oil per mole of original fluid + V_oil_sc_per_mol_init = n_oil * (1.0 - props_sc.V) * V_mol_oil_sc + + # Compute Bo and Rs + # Bo = volume of oil at reservoir conditions / volume of stock tank oil + # Rs = volume of gas at standard conditions / volume of stock tank oil + V_gas_sc_per_mol = phase_molar_volume_from_Z(p_sc, T_sc, 1.0) # ideal gas at SC + + # Cumulative gas remaining in solution from step i downward + cum_gas_removed_below = zeros(n_points + 1) + for i in (n_points + 1):-1:2 + cum_gas_removed_below[i] = 0.0 + end + # Gas still in solution at step i = total gas from step i+1 to end + gas_in_solution = zeros(n_points + 1) + running_gas = 0.0 + for i in (n_points + 1):-1:1 + gas_in_solution[i] = running_gas + if i > 1 + running_gas += moles_gas_removed[i] + end + end + + for (i, p) in enumerate(pressures) + if V_oil_sc_per_mol_init > 0.0 + # Bo: volume of oil at (p,T) per unit stock tank oil volume + # At step i, we have n_oil_at_i moles of oil + # Compute n_oil_at_i + n_oil_i = 1.0 + for j in 2:i + n_oil_i *= (1.0 - flash_and_properties(eos, pressures[j], T, oil_comp[max(1, j - 1)]).V) + end + if i == 1 + Bo_arr[i] = oil_molar_volumes[i] / V_oil_sc_per_mol_init + else + Bo_arr[i] = n_oil_i * oil_molar_volumes[i] / V_oil_sc_per_mol_init + end + + # Rs: gas in solution at standard conditions / stock tank oil volume + Rs_arr[i] = gas_in_solution[i] * V_gas_sc_per_mol / V_oil_sc_per_mol_init + end + end + + return DLEResult(T, collect(Float64, z ./ sum(z)), pressures, + Bo_arr, Rs_arr, ρ_oil_arr, ρ_gas_arr, μ_oil_arr, μ_gas_arr, + Bg_arr, Z_gas_arr, gas_comp, oil_comp, p_sat) +end diff --git a/src/PVTExperiments/interface.jl b/src/PVTExperiments/interface.jl new file mode 100644 index 0000000..45b635b --- /dev/null +++ b/src/PVTExperiments/interface.jl @@ -0,0 +1,131 @@ +""" + generate_pvt_tables(eos, z, T_res; ) + +High-level interface that generates complete black oil PVT tables from a +compositional fluid description. + +Goes from a fluid sample, reservoir temperature and surface conditions to +PVTO/PVDG (or PVTG/PVDG) tables + surface densities. + +# Arguments +- `eos`: Equation of state (e.g., `GenericCubicEOS(mixture, PengRobinson())`) +- `z`: Overall mole fractions of reservoir fluid +- `T_res`: Reservoir temperature (K) + +# Keyword Arguments +- `disgas::Union{Bool,Nothing} = nothing`: Allow dissolved gas in oil phase. + If `nothing`, auto-detected from fluid type. +- `vapoil::Union{Bool,Nothing} = nothing`: Allow vaporized oil in gas phase. + If `nothing`, auto-detected from fluid type. +- `p_range`: Pressure range for tables. Default: (50e6, 1e5) +- `p_sc`: Standard condition pressure (Pa). Default: 101325.0 +- `T_sc`: Standard condition temperature (K). Default: 288.706 +- `separator_stages`: Separator stage definitions. If provided, MSS is used + for surface reference; otherwise DLE is used. +- `n_pvto`: Number of Rs levels for PVTO. Default: 15 +- `n_pvtg`: Number of pressure levels for PVTG. Default: 15 +- `n_pvdg`: Number of pressure points for PVDG. Default: 20 +- `n_undersaturated`: Undersaturated points per level. Default: 5 + +# Returns +A NamedTuple with fields: +- `pvto`: PVTOTable or `nothing` +- `pvtg`: PVTGTable or `nothing` +- `pvdg`: PVDGTable +- `surface_densities`: SurfaceDensities +- `saturation_pressure`: Saturation pressure (Pa) +- `is_bubblepoint`: Whether the fluid has a bubble point (oil) or dew point (gas) +""" +function generate_pvt_tables(eos, z, T_res; + disgas::Union{Bool,Nothing} = nothing, + vapoil::Union{Bool,Nothing} = nothing, + p_range = (50e6, 1e5), + p_sc = 101325.0, + T_sc = 288.706, + separator_stages = nothing, + n_pvto = 15, + n_pvtg = 15, + n_pvdg = 20, + n_undersaturated = 5 + ) + z = collect(Float64, z) + z ./= sum(z) + + # Find saturation pressure and fluid type + p_sat, is_bp = find_saturation_pressure(eos, z, T_res; + p_min = p_range[2], p_max = p_range[1]) + + # Auto-detect table types if not specified + if isnothing(disgas) + disgas = is_bp # Oil system → dissolved gas + end + if isnothing(vapoil) + vapoil = !is_bp # Gas system → vaporized oil + end + + # Surface stages + if isnothing(separator_stages) + surf_stages = [SeparatorStage(p_sc, T_sc)] + else + surf_stages = separator_stages + end + + # Generate tables + pvto_result = nothing + pvtg_result = nothing + + if disgas + pvto_result = pvto_table(eos, z, T_res; + p_range = p_range, + n_rs = n_pvto, + n_undersaturated = n_undersaturated, + p_sc = p_sc, + T_sc = T_sc, + separator_stages = isnothing(separator_stages) ? nothing : separator_stages + ) + end + + if vapoil + pvtg_result = pvtg_table(eos, z, T_res; + p_range = p_range, + n_rv = n_pvtg, + n_undersaturated = 3, + p_sc = p_sc, + T_sc = T_sc + ) + end + + # PVDG is always generated + # For oil systems, use the gas composition from the first liberation step + # For gas systems, use the overall composition + if is_bp + # Get gas composition at bubble point + props_bp = flash_and_properties(eos, p_sat * 0.95, T_res, z) + z_gas = props_bp.y + else + z_gas = z + end + + pvdg_result = pvdg_table(eos, z_gas, T_res; + p_range = p_range, + n_points = n_pvdg, + p_sc = p_sc, + T_sc = T_sc + ) + + # Surface densities + sd = surface_densities(eos, z, T_res; + p_sc = p_sc, + T_sc = T_sc, + separator_stages = separator_stages + ) + + return ( + pvto = pvto_result, + pvtg = pvtg_result, + pvdg = pvdg_result, + surface_densities = sd, + saturation_pressure = p_sat, + is_bubblepoint = is_bp + ) +end diff --git a/src/PVTExperiments/mss.jl b/src/PVTExperiments/mss.jl new file mode 100644 index 0000000..3cb6624 --- /dev/null +++ b/src/PVTExperiments/mss.jl @@ -0,0 +1,97 @@ +""" + mss(eos, z_feed, T_res, stages::Vector{SeparatorStage}; p_res) + +Perform a Multistage Separator Test (MSS). + +The feed fluid at reservoir conditions is flashed through a sequence of +separator stages at specified (p, T) conditions. Gas is removed at each stage, +and the liquid proceeds to the next stage. The final stage is the stock tank. + +Returns Bo, Rs (referenced to stock tank oil), and properties at each stage. + +# Arguments +- `eos`: Equation of state +- `z_feed`: Feed composition (overall mole fractions at reservoir conditions) +- `T_res`: Reservoir temperature (K) +- `stages`: Vector of `SeparatorStage(p, T)` (last stage is stock tank) +- `p_res`: Reservoir pressure (Pa). Default: will be determined from saturation pressure. +""" +function mss(eos, z_feed, T_res, stages::Vector{SeparatorStage}; + p_res = nothing + ) + z_feed = collect(Float64, z_feed) + z_feed ./= sum(z_feed) + nc = number_of_components(eos) + n_stages = length(stages) + + # Find reservoir pressure if not given + if isnothing(p_res) + p_res_val, _ = find_saturation_pressure(eos, z_feed, T_res) + else + p_res_val = p_res + end + + # Flash feed through separator stages + sep_result = separator_flash(eos, z_feed, stages) + + # Compute Bo: oil volume at reservoir conditions / stock tank oil volume + # Oil at reservoir conditions (1 mole basis) + props_res = flash_and_properties(eos, p_res_val, T_res, z_feed) + V_mol_oil_res = props_res.V_mol_l + + # Track liquid mole fraction through stages + n_oil = 1.0 # Start with 1 mole of feed + liquid_fraction_remaining = 1.0 + GOR_stage = zeros(n_stages) + + p_sc = stages[end].p + T_sc = stages[end].T + + for i in 1:n_stages + V_i = sep_result.V_stages[i] + n_gas_i = liquid_fraction_remaining * V_i + n_liq_i = liquid_fraction_remaining * (1.0 - V_i) + + # GOR at this stage: gas volume at SC / oil volume at SC + if n_liq_i > 0 + V_gas_sc = n_gas_i * phase_molar_volume_from_Z(p_sc, T_sc, 1.0) + GOR_stage[i] = V_gas_sc + end + + liquid_fraction_remaining = n_liq_i + end + + # Stock tank oil properties + z_st = sep_result.oil_compositions[end] + props_st = flash_and_properties(eos, p_sc, T_sc, z_st) + V_mol_oil_sc = props_st.V_mol_l + ρ_oil_st = props_st.ρ_l + + # Stock tank oil volume per initial mole + V_oil_st = liquid_fraction_remaining * V_mol_oil_sc + + # Bo = oil volume at reservoir / stock tank oil volume + Bo = V_mol_oil_res / V_oil_st + + # Rs = total gas volume at SC / stock tank oil volume + V_gas_total_sc = sum(GOR_stage) + Rs = V_gas_total_sc / V_oil_st + + # Normalize GOR per stage by stock tank oil volume + GOR_stage ./= V_oil_st + + # Gas density at stock tank + z_gas_st = sep_result.gas_compositions[end] + mw_gas_st = mixture_molar_mass(z_gas_st, eos) + V_mol_gas_sc = phase_molar_volume_from_Z(p_sc, T_sc, 1.0) + ρ_gas_st = mw_gas_st / V_mol_gas_sc + + # API gravity + sg_oil = ρ_oil_st / WATER_DENSITY_SC # specific gravity relative to water + API = 141.5 / sg_oil - 131.5 + + return MSSResult(stages, z_feed, Bo, Rs, + sep_result.gas_compositions, + sep_result.oil_compositions, + GOR_stage, ρ_oil_st, ρ_gas_st, API) +end diff --git a/src/PVTExperiments/tables.jl b/src/PVTExperiments/tables.jl new file mode 100644 index 0000000..fa8710d --- /dev/null +++ b/src/PVTExperiments/tables.jl @@ -0,0 +1,358 @@ +""" + pvto_table(eos, z, T; p_range, n_rs, n_undersaturated, p_sc, T_sc, separator_stages) + +Generate a PVTO (live oil) table from DLE or MSS experiment data. + +# Arguments +- `eos`: Equation of state +- `z`: Overall composition (mole fractions) +- `T`: Reservoir temperature (K) +- `p_range`: Pressure range (p_max, p_min). Default: (50e6, 1e5) +- `n_rs`: Number of Rs (dissolved gas) levels. Default: 15 +- `n_undersaturated`: Number of undersaturated pressure steps per Rs level. Default: 5 +- `p_sc`: Standard condition pressure (Pa). Default: 101325.0 +- `T_sc`: Standard condition temperature (K). Default: 288.706 +- `separator_stages`: Optional separator stages. If provided, MSS is used for surface conditions. +""" +function pvto_table(eos, z, T; + p_range = (50e6, 1e5), + n_rs = 15, + n_undersaturated = 5, + p_sc = 101325.0, + T_sc = 288.706, + separator_stages = nothing + ) + z = collect(Float64, z) + z ./= sum(z) + + # Find saturation pressure + p_sat, _ = find_saturation_pressure(eos, z, T; p_min = p_range[2], p_max = p_range[1]) + + # Set up surface flash conditions + if isnothing(separator_stages) + surface_stages = [SeparatorStage(p_sc, T_sc)] + else + surface_stages = separator_stages + end + + # DLE-like liberation: step from p_sat to p_min + pressures_lib = collect(range(p_sat, p_range[2], length = n_rs + 1)) + + # Storage for the table + Rs_values = Float64[] + p_bub_values = Float64[] + p_table = Vector{Vector{Float64}}() + Bo_table = Vector{Vector{Float64}}() + mu_o_table = Vector{Vector{Float64}}() + + z_oil = copy(z) + n_oil = 1.0 + V_gas_sc_per_mol = phase_molar_volume_from_Z(p_sc, T_sc, 1.0) + + # Track liberation data + lib_data = [] + + for (i, p_lib) in enumerate(pressures_lib) + if i == 1 + # At bubble point + props = flash_and_properties(eos, p_lib, T, z_oil) + push!(lib_data, (p = p_lib, z_oil = copy(z_oil), n_oil = n_oil, props = props)) + else + # Flash and liberate gas + props = flash_and_properties(eos, p_lib, T, z_oil) + n_gas = n_oil * props.V + n_oil = n_oil * (1.0 - props.V) + z_oil = copy(props.x) + push!(lib_data, (p = p_lib, z_oil = copy(z_oil), n_oil = n_oil, props = props)) + end + end + + # Compute stock tank oil reference from final liquid + sep_res = separator_flash(eos, z_oil, surface_stages) + z_st_oil = sep_res.oil_compositions[end] + props_st = flash_and_properties(eos, p_sc, T_sc, z_st_oil) + V_mol_st_oil = props_st.V_mol_l + n_oil_at_st = n_oil + for V_i in sep_res.V_stages + n_oil_at_st *= (1.0 - V_i) + end + + # Build the table in reverse (from dead oil up to live oil) + for (idx, ld) in enumerate(reverse(lib_data)) + p_bubble = ld.p + z_oil_at_level = ld.z_oil + + # Flash this oil through separators to get surface products + sep = separator_flash(eos, z_oil_at_level, surface_stages) + z_st = sep.oil_compositions[end] + props_st_level = flash_and_properties(eos, p_sc, T_sc, z_st) + + # Liquid fraction through separator + liq_frac = 1.0 + gas_total_sc = 0.0 + for (j, V_j) in enumerate(sep.V_stages) + gas_at_j = liq_frac * V_j + gas_total_sc += gas_at_j * V_gas_sc_per_mol + liq_frac *= (1.0 - V_j) + end + V_oil_st = liq_frac * props_st_level.V_mol_l + + if V_oil_st ≤ 0.0 + continue + end + + Rs = max(0.0, gas_total_sc / V_oil_st) + + # Properties at bubble point + props_bub = flash_and_properties(eos, p_bubble, T, z_oil_at_level) + Bo_bub = props_bub.V_mol_l / (V_oil_st / liq_frac) + + # Saturated entry + p_entries = Float64[p_bubble] + Bo_entries = Float64[Bo_bub] + mu_entries = Float64[props_bub.μ_l] + + # Undersaturated entries: higher pressures + if p_bubble < p_range[1] + p_max_us = min(p_range[1], 2.0 * p_bubble) + dp = (p_max_us - p_bubble) / n_undersaturated + for k in 1:n_undersaturated + p_us = p_bubble + k * dp + props_us = flash_and_properties(eos, p_us, T, z_oil_at_level) + Bo_us = props_us.V_mol_l / (V_oil_st / liq_frac) + push!(p_entries, p_us) + push!(Bo_entries, Bo_us) + push!(mu_entries, props_us.μ_l) + end + end + + push!(Rs_values, Rs) + push!(p_bub_values, p_bubble) + push!(p_table, p_entries) + push!(Bo_table, Bo_entries) + push!(mu_o_table, mu_entries) + end + + # Sort by Rs + order = sortperm(Rs_values) + Rs_values = Rs_values[order] + p_bub_values = p_bub_values[order] + p_table = p_table[order] + Bo_table = Bo_table[order] + mu_o_table = mu_o_table[order] + + return PVTOTable(Rs_values, p_bub_values, p_table, Bo_table, mu_o_table) +end + +""" + pvdg_table(eos, z, T; p_range, n_points, p_sc, T_sc) + +Generate a PVDG (dry gas) table. + +# Arguments +- `eos`: Equation of state +- `z`: Gas composition (mole fractions) +- `T`: Temperature (K) +- `p_range`: Pressure range. Default: (50e6, 1e5) +- `n_points`: Number of pressure points. Default: 20 +- `p_sc`: Standard condition pressure (Pa). Default: 101325.0 +- `T_sc`: Standard condition temperature (K). Default: 288.706 +""" +function pvdg_table(eos, z, T; + p_range = (50e6, 1e5), + n_points = 20, + p_sc = 101325.0, + T_sc = 288.706 + ) + z = collect(Float64, z) + z ./= sum(z) + + pressures = collect(range(p_range[2], p_range[1], length = n_points)) + sort!(pressures) + + Bg_arr = zeros(n_points) + mu_g_arr = zeros(n_points) + + V_mol_sc = phase_molar_volume_from_Z(p_sc, T_sc, 1.0) + + for (i, p) in enumerate(pressures) + props = flash_and_properties(eos, p, T, z) + V_mol_gas = props.V_mol_v + Bg_arr[i] = V_mol_gas / V_mol_sc + mu_g_arr[i] = props.μ_v + end + + return PVDGTable(pressures, Bg_arr, mu_g_arr) +end + +""" + pvtg_table(eos, z, T; p_range, n_rv, n_undersaturated, p_sc, T_sc) + +Generate a PVTG (wet gas / gas condensate) table. + +# Arguments +- `eos`: Equation of state +- `z`: Overall composition (mole fractions) +- `T`: Temperature (K) +- `p_range`: Pressure range. Default: (50e6, 1e5) +- `n_rv`: Number of pressure (Rv) levels. Default: 15 +- `n_undersaturated`: Number of undersaturated entries per level. Default: 3 +- `p_sc`: Standard condition pressure (Pa). Default: 101325.0 +- `T_sc`: Standard condition temperature (K). Default: 288.706 +""" +function pvtg_table(eos, z, T; + p_range = (50e6, 1e5), + n_rv = 15, + n_undersaturated = 3, + p_sc = 101325.0, + T_sc = 288.706 + ) + z = collect(Float64, z) + z ./= sum(z) + + # Find dew point + p_sat, _ = find_saturation_pressure(eos, z, T; p_min = p_range[2], p_max = p_range[1]) + + # Pressure steps + pressures = collect(range(p_sat, p_range[2], length = n_rv + 1)) + + V_mol_sc = phase_molar_volume_from_Z(p_sc, T_sc, 1.0) + + p_values = Float64[] + Rv_sat_values = Float64[] + Rv_sub_table = Vector{Vector{Float64}}() + Bg_table = Vector{Vector{Float64}}() + mu_g_table = Vector{Vector{Float64}}() + + # CVD-like approach for gas condensate + z_gas = copy(z) + + for (i, p) in enumerate(pressures) + props = flash_and_properties(eos, p, T, z_gas) + + # Saturated gas properties + V_mol_gas = props.V_mol_v + Bg_sat = V_mol_gas / V_mol_sc + + # Rv: liquid content in gas at surface conditions + # Flash the vapor phase at surface conditions + if props.state == MultiComponentFlash.two_phase_lv + y = props.y + else + y = copy(z_gas) + end + + props_gas_sc = flash_and_properties(eos, p_sc, T_sc, y) + # Rv = volume of oil from gas at SC / volume of gas at SC + if props_gas_sc.state == MultiComponentFlash.two_phase_lv + V_oil_from_gas = (1.0 - props_gas_sc.V) * props_gas_sc.V_mol_l + V_gas_at_sc = props_gas_sc.V * props_gas_sc.V_mol_v + if V_gas_at_sc > 0 + Rv = V_oil_from_gas / V_gas_at_sc + else + Rv = 0.0 + end + else + Rv = 0.0 + end + + # Saturated entry + rv_entries = Float64[Rv] + bg_entries = Float64[Bg_sat] + mu_entries = Float64[props.μ_v] + + # Undersaturated entries: reduce Rv (leaner gas) + if Rv > 0 + for k in 1:n_undersaturated + frac = 1.0 - k / (n_undersaturated + 1) + # Create a leaner gas by mixing with pure gas component + z_lean = frac .* y .+ (1.0 - frac) .* z_gas + z_lean ./= sum(z_lean) + props_lean = flash_and_properties(eos, p, T, z_lean) + Bg_lean = props_lean.V_mol_v / V_mol_sc + + # Get Rv for lean gas + props_lean_sc = flash_and_properties(eos, p_sc, T_sc, z_lean) + if props_lean_sc.state == MultiComponentFlash.two_phase_lv + V_oil_l = (1.0 - props_lean_sc.V) * props_lean_sc.V_mol_l + V_gas_l = props_lean_sc.V * props_lean_sc.V_mol_v + Rv_lean = V_gas_l > 0 ? V_oil_l / V_gas_l : 0.0 + else + Rv_lean = 0.0 + end + + push!(rv_entries, Rv_lean) + push!(bg_entries, Bg_lean) + push!(mu_entries, props_lean.μ_v) + end + end + + # Add zero Rv entry + push!(rv_entries, 0.0) + # Dry gas Bg and viscosity + # Use the lightest component (index 1) as a proxy for dry gas. + # This assumes the mixture is ordered with the lightest component first. + z_first = zeros(length(z)) + z_first[1] = 1.0 + props_dry = flash_and_properties(eos, p, T, z_first) + push!(bg_entries, props_dry.V_mol_v / V_mol_sc) + push!(mu_entries, props_dry.μ_v) + + push!(p_values, p) + push!(Rv_sat_values, Rv) + push!(Rv_sub_table, rv_entries) + push!(Bg_table, bg_entries) + push!(mu_g_table, mu_entries) + + # Update gas composition (CVD-like) + if props.state == MultiComponentFlash.two_phase_lv + z_gas = copy(props.y) + end + end + + return PVTGTable(p_values, Rv_sat_values, Rv_sub_table, Bg_table, mu_g_table) +end + +""" + surface_densities(eos, z, T; p_sc, T_sc, separator_stages) + +Compute surface densities of oil and gas at standard conditions. + +# Arguments +- `eos`: Equation of state +- `z`: Overall composition (mole fractions) +- `T`: Reservoir temperature (K) +- `p_sc`: Standard condition pressure (Pa). Default: 101325.0 +- `T_sc`: Standard condition temperature (K). Default: 288.706 +- `separator_stages`: Optional separator stages for surface conditions. +""" +function surface_densities(eos, z, T; + p_sc = 101325.0, + T_sc = 288.706, + separator_stages = nothing + ) + z = collect(Float64, z) + z ./= sum(z) + + if isnothing(separator_stages) + surface_stages = [SeparatorStage(p_sc, T_sc)] + else + surface_stages = separator_stages + end + + # Flash through separator stages + sep = separator_flash(eos, z, surface_stages) + + # Oil density at stock tank + z_oil_st = sep.oil_compositions[end] + props_oil = flash_and_properties(eos, p_sc, T_sc, z_oil_st) + ρ_oil = props_oil.ρ_l + + # Gas density at stock tank + z_gas_st = sep.gas_compositions[end] + mw_gas = mixture_molar_mass(z_gas_st, eos) + V_mol_gas = phase_molar_volume_from_Z(p_sc, T_sc, 1.0) + ρ_gas = mw_gas / V_mol_gas + + return SurfaceDensities(ρ_oil, ρ_gas) +end diff --git a/src/PVTExperiments/types.jl b/src/PVTExperiments/types.jl new file mode 100644 index 0000000..f7ac414 --- /dev/null +++ b/src/PVTExperiments/types.jl @@ -0,0 +1,213 @@ +""" + SeparatorStage(p, T) + +A single separator stage defined by pressure `p` (Pa) and temperature `T` (K). +""" +struct SeparatorStage + p::Float64 + T::Float64 +end + +""" + CCEResult + +Results from a Constant Composition Expansion experiment. + +# Fields +- `T`: Temperature (K) +- `z`: Overall composition (mole fractions) +- `pressures`: Pressure steps (Pa) +- `Z_factors`: Two-phase compressibility factors +- `V_factors`: Vapor mole fraction at each step +- `relative_volume`: Relative volume V/V_sat +- `density_oil`: Oil (liquid) density at each step (kg/m³) +- `density_gas`: Gas (vapor) density at each step (kg/m³) +- `viscosity_oil`: Oil viscosity (Pa·s) +- `viscosity_gas`: Gas viscosity (Pa·s) +- `p_sat`: Bubble/dew point pressure (Pa) +- `is_bubblepoint`: true if bubble point, false if dew point +""" +struct CCEResult + T::Float64 + z::Vector{Float64} + pressures::Vector{Float64} + Z_factors::Vector{Float64} + V_factors::Vector{Float64} + relative_volume::Vector{Float64} + density_oil::Vector{Float64} + density_gas::Vector{Float64} + viscosity_oil::Vector{Float64} + viscosity_gas::Vector{Float64} + p_sat::Float64 + is_bubblepoint::Bool +end + +""" + DLEResult + +Results from a Differential Liberation Expansion experiment. + +# Fields +- `T`: Temperature (K) +- `z_init`: Initial overall composition (mole fractions) +- `pressures`: Pressure steps (Pa) (from p_sat down to final pressure) +- `Bo`: Oil formation volume factor (m³/m³ at standard conditions) +- `Rs`: Solution gas-oil ratio (m³/m³ at standard conditions) +- `density_oil`: Oil density at each step (kg/m³) +- `density_gas`: Gas density at each step (kg/m³) +- `viscosity_oil`: Oil viscosity (Pa·s) +- `viscosity_gas`: Gas viscosity (Pa·s) +- `Bg`: Gas formation volume factor (m³/m³ at standard conditions) +- `Z_gas`: Gas compressibility factor +- `gas_composition`: Gas composition at each liberation step +- `oil_composition`: Oil composition at each liberation step +- `p_sat`: Bubble point pressure (Pa) +""" +struct DLEResult + T::Float64 + z_init::Vector{Float64} + pressures::Vector{Float64} + Bo::Vector{Float64} + Rs::Vector{Float64} + density_oil::Vector{Float64} + density_gas::Vector{Float64} + viscosity_oil::Vector{Float64} + viscosity_gas::Vector{Float64} + Bg::Vector{Float64} + Z_gas::Vector{Float64} + gas_composition::Vector{Vector{Float64}} + oil_composition::Vector{Vector{Float64}} + p_sat::Float64 +end + +""" + CVDResult + +Results from a Constant Volume Depletion experiment. + +# Fields +- `T`: Temperature (K) +- `z_init`: Initial overall composition (mole fractions) +- `pressures`: Pressure steps (Pa) +- `Z_factors`: Two-phase Z factors +- `liquid_dropout`: Liquid volume fraction (fraction of cell volume) +- `gas_produced_cumulative`: Cumulative gas produced (mole fraction of initial) +- `density_gas`: Gas density at each step (kg/m³) +- `viscosity_gas`: Gas viscosity (Pa·s) +- `Z_gas`: Gas phase compressibility factor +- `gas_composition`: Gas composition at each step +- `p_sat`: Dew point pressure (Pa) +""" +struct CVDResult + T::Float64 + z_init::Vector{Float64} + pressures::Vector{Float64} + Z_factors::Vector{Float64} + liquid_dropout::Vector{Float64} + gas_produced_cumulative::Vector{Float64} + density_gas::Vector{Float64} + viscosity_gas::Vector{Float64} + Z_gas::Vector{Float64} + gas_composition::Vector{Vector{Float64}} + p_sat::Float64 +end + +""" + MSSResult + +Results from a Multistage Separator Test. + +# Fields +- `stages`: Separator stages (p, T pairs) +- `z_feed`: Feed composition (mole fractions) +- `Bo`: Oil formation volume factor from reservoir to stock tank (m³/m³) +- `Rs`: Solution GOR from reservoir to stock tank (m³/m³) +- `gas_composition`: Gas composition at each stage +- `oil_composition`: Oil composition exiting each stage +- `GOR_stage`: Gas-oil ratio at each separator stage (m³/m³) +- `density_oil_st`: Stock tank oil density (kg/m³) +- `density_gas_st`: Gas density at stock tank (kg/m³) +- `API_gravity`: API gravity of stock tank oil +""" +struct MSSResult + stages::Vector{SeparatorStage} + z_feed::Vector{Float64} + Bo::Float64 + Rs::Float64 + gas_composition::Vector{Vector{Float64}} + oil_composition::Vector{Vector{Float64}} + GOR_stage::Vector{Float64} + density_oil_st::Float64 + density_gas_st::Float64 + API_gravity::Float64 +end + +""" + PVTOTable + +Black oil PVTO (live oil) table for dissolved gas oil systems. + +# Fields +- `Rs`: Solution gas-oil ratio values (m³/m³ at surface conditions) +- `p_bub`: Bubble point pressures for each Rs (Pa) +- `p`: Pressure values for undersaturated oil at each Rs level (Pa) +- `Bo`: Oil formation volume factor (m³/m³) +- `mu_o`: Oil viscosity (Pa·s) +""" +struct PVTOTable + Rs::Vector{Float64} + p_bub::Vector{Float64} + p::Vector{Vector{Float64}} + Bo::Vector{Vector{Float64}} + mu_o::Vector{Vector{Float64}} +end + +""" + PVDGTable + +Black oil PVDG (dry gas) table. + +# Fields +- `p`: Pressure values (Pa) +- `Bg`: Gas formation volume factor (m³/m³) +- `mu_g`: Gas viscosity (Pa·s) +""" +struct PVDGTable + p::Vector{Float64} + Bg::Vector{Float64} + mu_g::Vector{Float64} +end + +""" + PVTGTable + +Black oil PVTG (wet gas / gas condensate) table. + +# Fields +- `p`: Pressure values (Pa) +- `Rv`: Vaporized oil-gas ratio at each pressure (m³/m³ at surface conditions) +- `Rv_sub`: Undersaturated Rv values at each pressure level +- `Bg`: Gas formation volume factor for each (p, Rv) entry (m³/m³) +- `mu_g`: Gas viscosity for each (p, Rv) entry (Pa·s) +""" +struct PVTGTable + p::Vector{Float64} + Rv::Vector{Float64} + Rv_sub::Vector{Vector{Float64}} + Bg::Vector{Vector{Float64}} + mu_g::Vector{Vector{Float64}} +end + +""" + SurfaceDensities + +Surface densities of oil and gas at standard conditions. + +# Fields +- `oil`: Oil density at surface conditions (kg/m³) +- `gas`: Gas density at surface conditions (kg/m³) +""" +struct SurfaceDensities + oil::Float64 + gas::Float64 +end diff --git a/src/PVTExperiments/utils.jl b/src/PVTExperiments/utils.jl new file mode 100644 index 0000000..c5bef40 --- /dev/null +++ b/src/PVTExperiments/utils.jl @@ -0,0 +1,330 @@ +""" + find_saturation_pressure(eos, z, T; p_min, p_max, n_points) + +Find the saturation pressure (bubble or dew point) for a given composition +and temperature using bisection. + +Returns `(p_sat, is_bubblepoint)`. +""" +function find_saturation_pressure(eos, z, T; + p_min = 1e5, + p_max = 100e6, + n_points = 50, + tol = 1e4 + ) + # Scan from high to low pressure to find where the phase transition occurs + pressures = range(p_max, p_min, length = n_points) + z_work = copy(z) + method = SSIFlash() + S = flash_storage(eos, (p = p_max, T = T, z = z_work), method = method) + K = initial_guess_K(eos, (p = p_max, T = T, z = z_work)) + + # Find the first pressure where two phases exist + p_upper = p_max + p_lower = p_min + found_two_phase = false + found_single = false + + for p in pressures + cond = (p = p, T = T, z = z_work) + K .= 1.0 + initial_guess_K!(K, eos, cond) + V, K_out, rep = flash_2ph!(S, K, eos, cond; extra_out = true, check = false) + if isnan(V) + if found_two_phase + p_lower = p + break + end + found_single = true + p_upper = p + else + if !found_two_phase && found_single + p_upper = pressures[max(1, findfirst(==(p), pressures) - 1)] + end + found_two_phase = true + p_lower = p + end + end + + if !found_two_phase + # Mixture is single-phase everywhere in the range + # Use label to determine phase type + label = single_phase_label(eos, (p = p_max, T = T, z = z_work)) + return (p_min, label < 0.5) + end + + # Bisection to refine + for _ in 1:60 # max bisection iterations + p_mid = 0.5 * (p_upper + p_lower) + cond = (p = p_mid, T = T, z = z_work) + K .= 1.0 + initial_guess_K!(K, eos, cond) + V, _, _ = flash_2ph!(S, K, eos, cond; extra_out = true, check = false) + if isnan(V) + p_upper = p_mid + else + p_lower = p_mid + end + if abs(p_upper - p_lower) < tol + break + end + end + + p_sat = 0.5 * (p_upper + p_lower) + + # Determine if bubble point or dew point + # Check composition: flash just below p_sat + cond = (p = p_sat * 0.99, T = T, z = z_work) + K .= 1.0 + initial_guess_K!(K, eos, cond) + V, K_out, _ = flash_2ph!(S, K, eos, cond; extra_out = true, check = false) + if isnan(V) + # Barely two-phase, use label + label = single_phase_label(eos, (p = p_sat, T = T, z = z_work)) + is_bubblepoint = label < 0.5 + else + is_bubblepoint = V < 0.5 + end + + return (p_sat, is_bubblepoint) +end + +""" + flash_and_properties(eos, p, T, z) + +Perform a flash calculation and extract all relevant properties. +Returns a NamedTuple with phase properties or `nothing` if flash fails. +""" +function flash_and_properties(eos, p, T, z) + z_work = copy(z) + result = flashed_mixture_2ph(eos, (p = p, T = T, z = z_work); method = SSIFlash()) + state = result.state + V = result.V + + if state == MultiComponentFlash.two_phase_lv + x = copy(result.liquid.mole_fractions) + y = copy(result.vapor.mole_fractions) + Z_L = result.liquid.Z + Z_V = result.vapor.Z + ρ_l = mass_density(eos, p, T, result.liquid) + ρ_v = mass_density(eos, p, T, result.vapor) + μ_l = lbc_viscosity(eos, p, T, result.liquid) + μ_v = lbc_viscosity(eos, p, T, result.vapor) + V_mol_l = molar_volume(eos, p, T, result.liquid) + V_mol_v = molar_volume(eos, p, T, result.vapor) + elseif state == MultiComponentFlash.single_phase_l + x = copy(z_work) + y = copy(z_work) + ph = result.liquid + Z_L = ph.Z + Z_V = ph.Z + ρ_l = mass_density(eos, p, T, ph) + ρ_v = ρ_l + μ_l = lbc_viscosity(eos, p, T, ph) + μ_v = μ_l + V_mol_l = molar_volume(eos, p, T, ph) + V_mol_v = V_mol_l + V = 0.0 + else + x = copy(z_work) + y = copy(z_work) + ph = result.vapor + Z_L = ph.Z + Z_V = ph.Z + ρ_l = mass_density(eos, p, T, ph) + ρ_v = ρ_l + μ_l = lbc_viscosity(eos, p, T, ph) + μ_v = μ_l + V_mol_l = molar_volume(eos, p, T, ph) + V_mol_v = V_mol_l + V = 1.0 + end + + return ( + state = state, + V = V, + x = x, + y = y, + Z_L = Z_L, + Z_V = Z_V, + ρ_l = ρ_l, + ρ_v = ρ_v, + μ_l = μ_l, + μ_v = μ_v, + V_mol_l = V_mol_l, + V_mol_v = V_mol_v + ) +end + +""" + mixture_molar_mass(z, eos) + +Compute the molar mass of a mixture given composition z. +""" +function mixture_molar_mass(z, eos) + mw = 0.0 + props = eos.mixture.properties + for (i, prop) in enumerate(props) + mw += z[i] * prop.mw + end + return mw +end + +""" + phase_molar_volume_from_Z(p, T, Z) + +Compute molar volume from compressibility factor: V = ZRT/p +""" +function phase_molar_volume_from_Z(p, T, Z) + return Z * R_GAS * T / p +end + +""" + separator_flash(eos, z_feed, stages::Vector{SeparatorStage}) + +Run a multistage separator train. Returns the compositions and properties at each stage. +""" +function separator_flash(eos, z_feed, stages::Vector{SeparatorStage}) + n_stages = length(stages) + nc = number_of_components(eos) + + gas_compositions = Vector{Vector{Float64}}(undef, n_stages) + oil_compositions = Vector{Vector{Float64}}(undef, n_stages) + V_stages = zeros(n_stages) + + z_current = copy(z_feed) + + for (i, stage) in enumerate(stages) + props = flash_and_properties(eos, stage.p, stage.T, z_current) + if props.state == MultiComponentFlash.two_phase_lv + gas_compositions[i] = copy(props.y) + oil_compositions[i] = copy(props.x) + V_stages[i] = props.V + z_current = copy(props.x) # liquid goes to next stage + else + # Single phase - no separation + gas_compositions[i] = copy(z_current) + oil_compositions[i] = copy(z_current) + V_stages[i] = props.V + # z_current stays the same + end + end + + return ( + gas_compositions = gas_compositions, + oil_compositions = oil_compositions, + V_stages = V_stages + ) +end + +# Print functions +function Base.show(io::IO, r::CCEResult) + println(io, "CCE Experiment Results") + println(io, "=" ^ 60) + @printf(io, "Temperature: %.2f K (%.2f °C)\n", r.T, r.T - 273.15) + @printf(io, "Saturation pressure: %.4e Pa (%.2f bar)\n", r.p_sat, r.p_sat / 1e5) + println(io, r.is_bubblepoint ? "Type: Bubble point" : "Type: Dew point") + println(io, "-" ^ 60) + @printf(io, "%14s %10s %10s %10s %12s\n", "P (Pa)", "V_rel", "V_frac", "ρ_oil", "ρ_gas") + @printf(io, "%14s %10s %10s %10s %12s\n", "", "", "", "(kg/m³)", "(kg/m³)") + println(io, "-" ^ 60) + for i in eachindex(r.pressures) + @printf(io, "%14.4e %10.4f %10.4f %10.2f %12.2f\n", + r.pressures[i], r.relative_volume[i], r.V_factors[i], + r.density_oil[i], r.density_gas[i]) + end +end + +function Base.show(io::IO, r::DLEResult) + println(io, "DLE Experiment Results") + println(io, "=" ^ 60) + @printf(io, "Temperature: %.2f K (%.2f °C)\n", r.T, r.T - 273.15) + @printf(io, "Bubble point pressure: %.4e Pa (%.2f bar)\n", r.p_sat, r.p_sat / 1e5) + println(io, "-" ^ 60) + @printf(io, "%14s %10s %10s %10s %10s\n", "P (Pa)", "Bo", "Rs", "ρ_oil", "Bg") + @printf(io, "%14s %10s %10s %10s %10s\n", "", "(m³/m³)", "(m³/m³)", "(kg/m³)", "(m³/m³)") + println(io, "-" ^ 60) + for i in eachindex(r.pressures) + @printf(io, "%14.4e %10.4f %10.4f %10.2f %10.6f\n", + r.pressures[i], r.Bo[i], r.Rs[i], + r.density_oil[i], r.Bg[i]) + end +end + +function Base.show(io::IO, r::CVDResult) + println(io, "CVD Experiment Results") + println(io, "=" ^ 60) + @printf(io, "Temperature: %.2f K (%.2f °C)\n", r.T, r.T - 273.15) + @printf(io, "Dew point pressure: %.4e Pa (%.2f bar)\n", r.p_sat, r.p_sat / 1e5) + println(io, "-" ^ 60) + @printf(io, "%14s %10s %12s %10s\n", "P (Pa)", "Z_gas", "Liq dropout", "Gas prod") + println(io, "-" ^ 60) + for i in eachindex(r.pressures) + @printf(io, "%14.4e %10.4f %12.4f %10.4f\n", + r.pressures[i], r.Z_gas[i], + r.liquid_dropout[i], r.gas_produced_cumulative[i]) + end +end + +function Base.show(io::IO, r::MSSResult) + println(io, "Multistage Separator Test Results") + println(io, "=" ^ 60) + @printf(io, "Number of stages: %d\n", length(r.stages)) + @printf(io, "Bo = %.4f m³/m³\n", r.Bo) + @printf(io, "Rs = %.4f m³/m³\n", r.Rs) + @printf(io, "Stock tank oil density: %.2f kg/m³\n", r.density_oil_st) + @printf(io, "API gravity: %.1f\n", r.API_gravity) + println(io, "-" ^ 60) + @printf(io, "%6s %14s %10s %12s\n", "Stage", "P (Pa)", "T (K)", "GOR (m³/m³)") + println(io, "-" ^ 60) + for i in eachindex(r.stages) + @printf(io, "%6d %14.4e %10.2f %12.4f\n", + i, r.stages[i].p, r.stages[i].T, r.GOR_stage[i]) + end +end + +function Base.show(io::IO, t::PVTOTable) + println(io, "PVTO Table") + println(io, "=" ^ 60) + @printf(io, "%10s %14s %12s %12s\n", "Rs", "P (Pa)", "Bo", "μ_o (Pa·s)") + println(io, "-" ^ 60) + for i in eachindex(t.Rs) + for j in eachindex(t.p[i]) + rs_str = j == 1 ? @sprintf("%10.4f", t.Rs[i]) : " " ^ 10 + @printf(io, "%s %14.4e %12.6f %12.4e\n", + rs_str, t.p[i][j], t.Bo[i][j], t.mu_o[i][j]) + end + println(io, "/") + end +end + +function Base.show(io::IO, t::PVDGTable) + println(io, "PVDG Table") + println(io, "=" ^ 60) + @printf(io, "%14s %12s %12s\n", "P (Pa)", "Bg", "μ_g (Pa·s)") + println(io, "-" ^ 60) + for i in eachindex(t.p) + @printf(io, "%14.4e %12.6f %12.4e\n", t.p[i], t.Bg[i], t.mu_g[i]) + end +end + +function Base.show(io::IO, t::PVTGTable) + println(io, "PVTG Table") + println(io, "=" ^ 60) + @printf(io, "%14s %10s %12s %12s\n", "P (Pa)", "Rv", "Bg", "μ_g (Pa·s)") + println(io, "-" ^ 60) + for i in eachindex(t.p) + for j in eachindex(t.Rv_sub[i]) + p_str = j == 1 ? @sprintf("%14.4e", t.p[i]) : " " ^ 14 + @printf(io, "%s %10.6f %12.6f %12.4e\n", + p_str, t.Rv_sub[i][j], t.Bg[i][j], t.mu_g[i][j]) + end + println(io, "/") + end +end + +function Base.show(io::IO, s::SurfaceDensities) + println(io, "Surface Densities") + @printf(io, " Oil: %.2f kg/m³\n", s.oil) + @printf(io, " Gas: %.4f kg/m³\n", s.gas) +end diff --git a/src/eos_types.jl b/src/eos_types.jl index 8985420..748ada2 100644 --- a/src/eos_types.jl +++ b/src/eos_types.jl @@ -88,9 +88,12 @@ struct RedlichKwong <: AbstractGeneralizedCubic end Specialized cubic equation of state for Søreide-Whitson EOS. -See "Peng-Robinson predictions for hydrocarbons, CO2, N2, and H₂S with pure -water and NaCl brine" by I. Søreide and C. Whitson, Fluid Phase Equilibria 77, -1992. +See "Peng-Robinson predictions for hydrocarbons, CO₂, N₂, and H₂S with pure +water and NaCl brine" by I. Søreide and C. Whitson, Fluid Phase Equilibria 77 (1992). + +Note that this EOS requires water to be present in the mixture by name (it is +after all the whole point of the Søreide-Whitson EOS). Special treatment is +used for H2S, CO2, and N2 when given by name. """ struct SoreideWhitson{T} <: AbstractPengRobinson A::NTuple{3, T} diff --git a/test/pvt_experiments.jl b/test/pvt_experiments.jl new file mode 100644 index 0000000..a18c411 --- /dev/null +++ b/test/pvt_experiments.jl @@ -0,0 +1,209 @@ +using Test +using MultiComponentFlash +import MultiComponentFlash.PVTExperiments as PVTExp + +@testset "PVTExperiments" begin + # Set up a simple oil mixture (methane + n-decane) + mixture = MultiComponentMixture(["Methane", "n-Decane"]) + eos = GenericCubicEOS(mixture, PengRobinson()) + z_oil = [0.5, 0.5] + T_res = 373.15 # 100°C + + # Gas condensate mixture (methane-heavy) + z_gas = [0.9, 0.1] + + @testset "Saturation Pressure" begin + p_sat, is_bp = PVTExp.find_saturation_pressure(eos, z_oil, T_res) + @test p_sat > 0.0 + @test is_bp # Should be bubble point for oil mixture + + p_sat_gas, is_bp_gas = PVTExp.find_saturation_pressure(eos, z_gas, T_res) + @test p_sat_gas > 0.0 + end + + @testset "Flash and Properties" begin + props = PVTExp.flash_and_properties(eos, 5e6, T_res, z_oil) + @test haskey(props, :V) + @test haskey(props, :ρ_l) + @test haskey(props, :ρ_v) + @test haskey(props, :μ_l) + @test haskey(props, :μ_v) + @test props.ρ_l > 0 + @test props.ρ_v > 0 + @test props.μ_l > 0 + @test props.μ_v > 0 + end + + @testset "CCE" begin + result = PVTExp.cce(eos, z_oil, T_res; n_points = 10) + @test result isa PVTExp.CCEResult + @test result.T == T_res + @test length(result.pressures) == 10 + @test all(result.density_oil .> 0) + @test all(result.density_gas .> 0) + @test all(result.viscosity_oil .> 0) + @test result.p_sat > 0 + # Relative volume should be 1.0 near saturation pressure + idx_sat = argmin(abs.(result.pressures .- result.p_sat)) + @test isapprox(result.relative_volume[idx_sat], 1.0, atol = 0.1) + # Test printing + io = IOBuffer() + show(io, result) + output = String(take!(io)) + @test contains(output, "CCE") + end + + @testset "DLE" begin + result = PVTExp.dle(eos, z_oil, T_res; n_points = 10) + @test result isa PVTExp.DLEResult + @test result.T == T_res + @test length(result.pressures) == 11 # n_points + 1 + @test result.p_sat > 0 + # Bo at bubble point should be > 1 + @test result.Bo[1] >= 1.0 + # Rs should decrease with pressure + @test result.Rs[1] >= result.Rs[end] + # All densities should be positive + @test all(result.density_oil .> 0) + # Test printing + io = IOBuffer() + show(io, result) + output = String(take!(io)) + @test contains(output, "DLE") + end + + @testset "CVD" begin + result = PVTExp.cvd(eos, z_gas, T_res; n_points = 10) + @test result isa PVTExp.CVDResult + @test result.T == T_res + @test length(result.pressures) == 11 # n_points + 1 + @test result.p_sat > 0 + # Liquid dropout starts at 0 (dew point) + @test result.liquid_dropout[1] ≈ 0.0 + # Gas Z-factors should be positive + @test all(result.Z_gas .> 0) + # Test printing + io = IOBuffer() + show(io, result) + output = String(take!(io)) + @test contains(output, "CVD") + end + + @testset "MSS" begin + stages = [ + PVTExp.SeparatorStage(3e6, 330.0), # High pressure separator + PVTExp.SeparatorStage(101325.0, 288.706) # Stock tank + ] + result = PVTExp.mss(eos, z_oil, T_res, stages) + @test result isa PVTExp.MSSResult + @test result.Bo > 0 + @test result.Rs >= 0 + @test result.density_oil_st > 0 + @test result.density_gas_st > 0 + @test length(result.gas_composition) == 2 + @test length(result.oil_composition) == 2 + # Test printing + io = IOBuffer() + show(io, result) + output = String(take!(io)) + @test contains(output, "Separator") + end + + @testset "PVTO Table" begin + table = PVTExp.pvto_table(eos, z_oil, T_res; n_rs = 5, n_undersaturated = 2) + @test table isa PVTExp.PVTOTable + @test length(table.Rs) == length(table.p_bub) + @test length(table.Rs) > 0 + # Rs should be non-negative and sorted + @test all(table.Rs .>= 0) + @test issorted(table.Rs) + # Bo should be positive + for bo_vec in table.Bo + @test all(bo_vec .> 0) + end + # Viscosity should be positive + for mu_vec in table.mu_o + @test all(mu_vec .> 0) + end + # Test printing + io = IOBuffer() + show(io, table) + output = String(take!(io)) + @test contains(output, "PVTO") + end + + @testset "PVDG Table" begin + table = PVTExp.pvdg_table(eos, z_gas, T_res; n_points = 10) + @test table isa PVTExp.PVDGTable + @test length(table.p) == 10 + @test all(table.Bg .> 0) + @test all(table.mu_g .> 0) + # Bg should decrease with pressure + @test table.Bg[1] > table.Bg[end] + # Test printing + io = IOBuffer() + show(io, table) + output = String(take!(io)) + @test contains(output, "PVDG") + end + + @testset "PVTG Table" begin + table = PVTExp.pvtg_table(eos, z_gas, T_res; n_rv = 5, n_undersaturated = 2) + @test table isa PVTExp.PVTGTable + @test length(table.p) > 0 + # Bg should be positive + for bg_vec in table.Bg + @test all(bg_vec .> 0) + end + # Test printing + io = IOBuffer() + show(io, table) + output = String(take!(io)) + @test contains(output, "PVTG") + end + + @testset "Surface Densities" begin + sd = PVTExp.surface_densities(eos, z_oil, T_res) + @test sd isa PVTExp.SurfaceDensities + @test sd.oil > 0 + @test sd.gas > 0 + # Oil should be denser than gas at surface + @test sd.oil > sd.gas + # Test printing + io = IOBuffer() + show(io, sd) + output = String(take!(io)) + @test contains(output, "Surface") + end + + @testset "High-Level Interface - Oil" begin + tables = generate_pvt_tables(eos, z_oil, T_res; + n_pvto = 5, n_pvdg = 10, n_undersaturated = 2) + @test tables.pvto !== nothing + @test tables.pvdg !== nothing + @test tables.surface_densities isa PVTExp.SurfaceDensities + @test tables.saturation_pressure > 0 + @test tables.is_bubblepoint == true + end + + @testset "High-Level Interface - Gas" begin + tables = generate_pvt_tables(eos, z_gas, T_res; + n_pvtg = 5, n_pvdg = 10) + @test tables.pvdg !== nothing + @test tables.surface_densities isa PVTExp.SurfaceDensities + @test tables.saturation_pressure > 0 + end + + @testset "High-Level Interface with Separator" begin + stages = [ + PVTExp.SeparatorStage(3e6, 330.0), + PVTExp.SeparatorStage(101325.0, 288.706) + ] + tables = generate_pvt_tables(eos, z_oil, T_res; + separator_stages = stages, + n_pvto = 5, n_pvdg = 10, n_undersaturated = 2) + @test tables.pvto !== nothing + @test tables.pvdg !== nothing + @test tables.surface_densities isa PVTExp.SurfaceDensities + end +end