Add table generator + improve docs (#27)

* Add table generator + improve docs

* Update Project.toml

Author: moyner

Project.toml

@@ -1,25 +1,27 @@
 name = "MultiComponentFlash"
 uuid = "35e5bd01-9722-4017-9deb-64a5d32478ff"
+version = "1.1.19"
 authors = ["Olav Møyner <olav.moyner@gmail.com>"]
-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"]

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
+```

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

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

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

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