diff --git a/docs/Project.toml b/docs/Project.toml index a99833c..134947b 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,5 +1,7 @@ [deps] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +Clapeyron = "7c7805af-46cc-48c9-995b-ed0ed2dc909a" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" MultiComponentFlash = "35e5bd01-9722-4017-9deb-64a5d32478ff" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/docs/make.jl b/docs/make.jl index d3551f5..7227a2e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,12 +1,15 @@ using MultiComponentFlash using Documenter -using BenchmarkTools, Plots +using BenchmarkTools, CairoMakie + +# CairoMakie.activate!(type = "svg") DocMeta.setdocmeta!(MultiComponentFlash, :DocTestSetup, :(using MultiComponentFlash); recursive=true) makedocs(; modules=[MultiComponentFlash], - warnonly = true, + warnonly = false, + checkdocs=:exports, authors="Olav Møyner ", repo="https://github.com/moyner/MultiComponentFlash.jl/blob/{commit}{path}#{line}", sitename="MultiComponentFlash.jl", @@ -19,7 +22,8 @@ makedocs(; "Home" => "index.md", "Examples" => Any[ "Basic usage" => "examples/basics.md", - "Advanced usage" => "examples/advanced.md" + "Advanced usage" => "examples/advanced.md", + "EOS reference validation" => "examples/eos_validation.md", ], "API" => Any[ "Mixtures" => "api/mixtures.md", diff --git a/docs/src/examples/advanced.md b/docs/src/examples/advanced.md index a72d4d5..3297c0b 100644 --- a/docs/src/examples/advanced.md +++ b/docs/src/examples/advanced.md @@ -86,9 +86,11 @@ between kernel work items. We create a three-component mixture and flash for a range of pressure and temperature conditions: -```julia -using MultiComponentFlash, Plots -ns = 1000 +```@example phase-diagram +using MultiComponentFlash, CairoMakie +CairoMakie.activate!(type = "svg") + +ns = 100 ubar = 1e5 # Pressure range p0 = 1*ubar @@ -110,20 +112,29 @@ cond = (p = p0, T = T0, z = z) m = SSIFlash() S = flash_storage(eos, cond, method = m) K = initial_guess_K(eos, cond) -data = zeros(ns, ns) -for ip = 1:ns - for iT = 1:ns - c = (p = p[ip], T = T[iT], z = z) - data[ip, iT] = flash_2ph!(S, K, eos, c, NaN, method = m) +data = zeros(length(T), length(p)) +for (iT, temperature) in pairs(T) + for (ip, pressure) in pairs(p) + c = (p = pressure, T = temperature, z = z) + V = flash_2ph!(S, K, eos, c, NaN, method = m) + data[iT, ip] = V end end -contour(p./ubar, T .- 273.15, data, levels = 10, fill=(true,cgrad(:hot))) -ylabel!("Pressure [Bar]") -xlabel!("T [°Celsius]") +fig = Figure(size = (760, 480)) +ax = Axis(fig[1, 1]; + xlabel = "Temperature [°C]", + ylabel = "Pressure [bar]", + title = "Vapor fraction") +contours = contourf!(ax, T .- 273.15, p./ubar, data; + levels = range(0.0, 1.0, length = 11), colormap = :hot) +Colorbar(fig[1, 2], contours; label = "Vapor mole fraction") +fig ``` -![Phase diagram](../assets/phase_diagram_simple.png) +The colored region is the two-phase envelope. Unfilled states are stable +single-phase conditions, for which the two-phase solver reports no intermediate +vapor fraction. ### PVT table generation @@ -132,3 +143,11 @@ There is experimental support for generating simulator input blackoil tables (e. ```@docs generate_pvt_tables ``` + +### Coupling to simulators and other utilities + +```@docs +cubic_benchmark +FlashedMixture2Phase +FlashedPhase +``` diff --git a/docs/src/examples/eos_validation.md b/docs/src/examples/eos_validation.md new file mode 100644 index 0000000..db086ec --- /dev/null +++ b/docs/src/examples/eos_validation.md @@ -0,0 +1,175 @@ +# EOS reference validation + +This example compares Peng-Robinson 1976 (PR76), Soave-Redlich-Kwong (SRK), and +Redlich-Kwong (RK) with +[Clapeyron.jl](https://github.com/ClapeyronThermo/Clapeyron.jl). It checks +compressibility and fugacity coefficients for a set of cases. + +## Comparison matrix + +```@example eos-reference +using Clapeyron, MultiComponentFlash, Test + +names = ["carbon dioxide", "methane", "decane"] +properties = ( + MolecularProperty(0.0440, 7.38e6, 304.1, 9.412e-5, 0.224), + MolecularProperty(0.0160, 4.60e6, 190.6, 9.863e-5, 0.011), + MolecularProperty(0.1420, 2.10e6, 617.7, 6.098e-4, 0.488), +) +conditions = ( + (p = 1.0e6, T = 300.0, z = [0.5, 0.3, 0.2]), + (p = 5.0e6, T = 350.0, z = [0.2, 0.7, 0.1]), + (p = 2.0e7, T = 500.0, z = [0.1, 0.2, 0.7]), +) +interactions = ( + zeros(3, 3), + [0.0 0.08 0.03; 0.08 0.0 0.015; 0.03 0.015 0.0], +) + +function compare_cubic_eos(names, properties, conditions, interactions) + relative_error(value, reference) = + abs(value - reference)/max(abs(reference), eps(Float64)) + Tc = [property.T_c for property in properties] + Pc = [property.p_c for property in properties] + Mw = [1000property.mw for property in properties] + acentricfactor = [property.ω for property in properties] + comparisons = NamedTuple[] + + for k in interactions + mixture = MultiComponentMixture(properties; names = names, A_ij = k) + parameters = (; Tc, Pc, Mw, acentricfactor, k, l = zeros(size(k))) + models = ( + (label = "PR76", eos = PengRobinson(), + reference = Clapeyron.PR(names; + userlocations = parameters, verbose = false)), + (label = "SRK", eos = SoaveRedlichKwong(), + reference = Clapeyron.SRK(names; + userlocations = parameters, verbose = false)), + (label = "RK", eos = RedlichKwong(), + reference = Clapeyron.RK(names; + userlocations = parameters, verbose = false)), + ) + + for model in models, state in conditions + eos = GenericCubicEOS(mixture, model.eos) + for (phase, reference_phase) in ((:liquid, :l), (:vapor, :v)) + state_with_phase = (; + p = state.p, T = state.T, z = state.z, phase) + forces = force_coefficients(eos, state_with_phase) + scalars = force_scalars(eos, state_with_phase, forces) + Z = mixture_compressibility_factor( + eos, state_with_phase, forces, scalars) + phi = [exp(MultiComponentFlash.component_fugacity_coefficient( + eos, state_with_phase, i, Z, forces, scalars)) + for i in eachindex(state.z)] + + reference_volume = Clapeyron.volume( + model.reference, state.p, state.T, state.z; + phase = reference_phase) + reference_Z = state.p*reference_volume/ + (Clapeyron.Rgas(model.reference)*state.T*sum(state.z)) + reference_phi = Clapeyron.fugacity_coefficient( + model.reference, state.p, state.T, state.z; + phase = reference_phase) + + push!(comparisons, (; + eos = model.label, + Z, + reference_Z, + Z_relative_error = relative_error(Z, reference_Z), + phi, + reference_phi, + phi_relative_error = relative_error.(phi, reference_phi))) + end + end + end + return comparisons +end + +comparisons = compare_cubic_eos( + names, properties, conditions, interactions) + +@test length(comparisons) == 36 +@test all(row -> isapprox(row.Z, row.reference_Z; + rtol = 1.0e-7, atol = 1.0e-10), comparisons) +@test all(row -> isapprox(row.phi, row.reference_phi; + rtol = 1.0e-7, atol = 1.0e-10), comparisons) + +length(comparisons) +``` + +## Parity plots + +The dashed line indicates exact agreement. Fugacity coefficients use logarithmic +axes because they span a wider range than ``Z``. + +```@example eos-reference +using CairoMakie +# CairoMakie.activate!(type = "svg") + +eos_names = ["PR76", "SRK", "RK"] +colors = [:dodgerblue3, :darkorange2, :seagreen4] + +fig = Figure(size = (920, 410)) +z_axis = Axis(fig[1, 1]; + xlabel = "Clapeyron Z", + ylabel = "MultiComponentFlash Z", + title = "Compressibility parity", + aspect = DataAspect()) +phi_axis = Axis(fig[1, 2]; + xlabel = "Clapeyron fugacity coefficient", + ylabel = "MultiComponentFlash fugacity coefficient", + title = "Fugacity-coefficient parity", + xscale = log10, + yscale = log10, + aspect = DataAspect()) + +for (eos_name, color) in zip(eos_names, colors) + rows = filter(row -> row.eos == eos_name, comparisons) + reference_z = [row.reference_Z for row in rows] + calculated_z = [row.Z for row in rows] + reference_phi = [value for row in rows for value in row.reference_phi] + calculated_phi = [value for row in rows for value in row.phi] + scatter!(z_axis, reference_z, calculated_z; + color = color, markersize = 9, label = eos_name) + scatter!(phi_axis, reference_phi, calculated_phi; + color = color, markersize = 9, label = eos_name) +end + +z_limits = extrema(vcat( + [row.reference_Z for row in comparisons], + [row.Z for row in comparisons])) +phi_limits = extrema(vcat( + [value for row in comparisons for value in row.reference_phi], + [value for row in comparisons for value in row.phi])) +z_line = collect(z_limits) +phi_line = collect(phi_limits) +lines!(z_axis, z_line, z_line; color = :black, linestyle = :dash) +lines!(phi_axis, phi_line, phi_line; color = :black, linestyle = :dash) +axislegend(z_axis; position = :lt) +fig +``` + +The maximum relative error makes the remaining differences visible. + +```@example eos-reference +maximum_z_error = [maximum(row.Z_relative_error + for row in comparisons if row.eos == eos_name) for eos_name in eos_names] +maximum_phi_error = [maximum(maximum(row.phi_relative_error) + for row in comparisons if row.eos == eos_name) for eos_name in eos_names] + +error_figure = Figure(size = (720, 420)) +error_axis = Axis(error_figure[1, 1]; + xlabel = "Equation of state", + ylabel = "Maximum relative error", + title = "Worst case across all states, roots, and BIC matrices", + yscale = log10, + xticks = (1:length(eos_names), eos_names)) +scatterlines!(error_axis, 1:length(eos_names), maximum_z_error; + markersize = 12, label = "Z") +scatterlines!(error_axis, 1:length(eos_names), maximum_phi_error; + markersize = 12, + label = "Fugacity coefficients") +axislegend(error_axis; position = :lt) +error_figure +``` diff --git a/docs/src/index.md b/docs/src/index.md index 1d8c42f..6407bf9 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -11,9 +11,12 @@ This package implements several equations of state for multicomponent vapor-liqu The following equations of state (EOS) are implemented as a class of generic cubics: * [Peng-Robinson](https://doi.org/10.1021/i160057a011) +* Corrected Peng-Robinson for large acentric factors * [Redlich-Kwong](https://doi.org/10.1021/cr60137a013) * [Soave-Redlich-Kwong](https://doi.org/10.1016/0009-2509(72)80096-4) * [Zudkevitch-Joffe](https://doi.org/10.1002/aic.690160122) +* [Soreide-Whitson](https://doi.org/10.1016/0378-3812(92)85105-H) for + water and brine mixtures The code is fully type stable, easy to use and fairly performant, with additional options to avoid allocations if you need to perform many flashes. The main implementation goal is to have a compact, performant and easy to use code suitable for integration in simulators of multiphase flow. @@ -38,15 +41,17 @@ If you want to use the former features, please know that they might be subject t * The module currently only supports cubic equations of state. These are limited in accuracy for longer chains of molecules without extensive tuning for a specific mixture. * Flash is limited to two-phase pressure-temperature (pT) flash. * The flash algorithms are limited to the basics - there are many strategies that could be implemented -* Peng-Robinson is the only equation of state that has been thoroughly validated. +* Cubic equations remain approximations whose practical accuracy depends on + component data and fitted binary-interaction coefficients. See + [EOS reference validation](@ref) for implementation-level comparisons. # Other packages Julia packages: -* [Clapeyron.jl](https://github.com/ypaul21/Clapeyron.jl) supports more advanced equations of state in addition to the cubics (SAFT-type and empirical EOS) and has a much larger API for thermodynamical properties. At the time of writing, this package does not support a full flash. +* [Clapeyron.jl](https://github.com/ypaul21/Clapeyron.jl) supports more advanced equations of state in addition to the cubics (SAFT-type and empirical EOS) and has a much larger API for thermodynamical properties. It is also used by this package's independent [EOS reference validation](@ref) suite. Julia interfaces to other useful packages: * [PyThermo.jl](https://github.com/stillyslalom/PyThermo.jl) is an interface to the [Thermo](https://pypi.org/project/thermo/) package. * [CoolProp.jl](https://github.com/CoolProp/CoolProp.jl) is an interface to the [CoolProp](http://www.coolprop.org/) library. # Contact -You can use the [Github webpage](https://github.com/moyner/MultiComponentFlash.jl) or drop me a line at [Olav Møyner](mailto:olav.moyner@gmail.no). +You can use the [Github webpage](https://github.com/moyner/MultiComponentFlash.jl) or drop me a line at [Olav Møyner](mailto:olav.moyner@sintef.no). diff --git a/src/PVTExperiments/PVTExperiments.jl b/src/PVTExperiments/PVTExperiments.jl index 83d0ef1..83761bc 100644 --- a/src/PVTExperiments/PVTExperiments.jl +++ b/src/PVTExperiments/PVTExperiments.jl @@ -21,5 +21,5 @@ module PVTExperiments include("tables.jl") include("interface.jl") - export generate_pvt_tables, PVTTableSet + export generate_pvt_tables end diff --git a/src/eos.jl b/src/eos.jl index d8264ef..78b9e71 100644 --- a/src/eos.jl +++ b/src/eos.jl @@ -159,18 +159,17 @@ function force_coefficients(eos::AbstractCubicEOS, cond; static_size = false) return coeff end -function get_force_coefficients(forces, eos, cond) - if forces_per_phase(eos) - phase = get_phase(cond) - if phase == :liquid - return forces.liquid - elseif phase == :vapor - return forces.vapor - else - error("Forces per phase are only supported for liquid and vapor phases, not $phase.") - end +@inline get_force_coefficients(forces, eos, cond) = forces + +@inline function get_force_coefficients( + forces::NamedTuple{(:liquid, :vapor)}, eos, cond) + phase = get_phase(cond) + if phase == :liquid + return forces.liquid + elseif phase == :vapor + return forces.vapor else - return forces + error("Forces per phase are only supported for liquid and vapor phases, not $phase.") end end @@ -363,7 +362,7 @@ function solve_cubic_positive_roots(a, b, c) M = R^2 - Q^3 single_root = M > 0 if single_root - # Single real roots + # Single real root S = -sign(R)*(abs(R) + sqrt(M))^(1/3) if S == 0 T = 0 @@ -371,9 +370,16 @@ function solve_cubic_positive_roots(a, b, c) T = Q/S end return S + T - a/3 + elseif iszero(Q) + # Q = R = 0: the polynomial has one triple root. The trigonometric + # expression below is otherwise undefined (acos(0/0)). + return -a/3 else # Three real roots - theta = acos(R/sqrt(Q^3)) + # Clamp guards against a round-off excursion outside [-1, 1] at a + # repeated root, where acos is still well defined analytically. + acos_arg = clamp(R/sqrt(Q^3), -one(R), one(R)) + theta = acos(acos_arg) r1 = -(2*sqrt(Q)*cos(theta/3)) - a/3 r2 = -(2*sqrt(Q)*cos((theta + 2*pi)/3)) - a/3 r3 = -(2*sqrt(Q)*cos((theta - 2*pi)/3)) - a/3 diff --git a/src/eos/soreide_whitson.jl b/src/eos/soreide_whitson.jl index ebcea0f..1113830 100644 --- a/src/eos/soreide_whitson.jl +++ b/src/eos/soreide_whitson.jl @@ -87,5 +87,5 @@ function soreide_whitson_hc_aqueous_bic(sw::SoreideWhitson, acf_i, T_ri) A_0 = a_0 + a_mw_0*sign(acf_i)*abs(acf_i)^(-0.1) A_1 = a_1 + a_mw_1*acf_i A_2 = a_2 + a_mw_2*acf_i - return A_0*(1 + α_0*c_sw) + A_1*T_ri*(1 + α_1*c_sw) + A_2*T_ri*(1 + α_2*c_sw) + return A_0*(1 + α_0*c_sw) + A_1*T_ri*(1 + α_1*c_sw) + A_2*T_ri^2*(1 + α_2*c_sw) end diff --git a/src/eos/zudkevitch_joffe.jl b/src/eos/zudkevitch_joffe.jl index ce94180..a94aef2 100644 --- a/src/eos/zudkevitch_joffe.jl +++ b/src/eos/zudkevitch_joffe.jl @@ -1,5 +1,5 @@ # ZudkevitchJoffe -function weight_ai(eos::GenericCubicEOS{ZudkevitchJoffe}, cond, i) +function weight_ai(eos::GenericCubicEOS{E}, cond, i) where E<:ZudkevitchJoffe zj = eos.type mix = eos.mixture T = cond.T @@ -7,7 +7,7 @@ function weight_ai(eos::GenericCubicEOS{ZudkevitchJoffe}, cond, i) return eos.ω_a*zj.F_a(T, i)*T_r^(-0.5) end -function weight_bi(eos::GenericCubicEOS{ZudkevitchJoffe}, cond, i) +function weight_bi(eos::GenericCubicEOS{E}, cond, i) where E<:ZudkevitchJoffe zj = eos.type T = cond.T return eos.ω_b*zj.F_b(T, i) diff --git a/src/static.jl b/src/static.jl index 7be38de..89a7409 100644 --- a/src/static.jl +++ b/src/static.jl @@ -93,8 +93,6 @@ end return max_root end -@inline get_force_coefficients(forces, eos::GenericCubicEOS, cond) = forces - """Immutable force coefficients for accelerator kernels.""" @inline function static_force_coefficients(eos::GenericCubicEOS{E, R, N}, cond, ::Type{T}) where {E, R, N, T} diff --git a/src/utils.jl b/src/utils.jl index c387497..b7edea4 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -24,7 +24,7 @@ function single_phase_label(mixture::MultiComponentMixture, cond) return Float64(cond.T > T_c) end -single_phase_label(eos::AbstractEOS,cond) = single_phase_label(eos.mixture,cond) +single_phase_label(eos::AbstractEOS, cond) = single_phase_label(eos.mixture,cond) """ lbc_viscosity(eos, p, T, ph; ) diff --git a/test/eos_equations.jl b/test/eos_equations.jl new file mode 100644 index 0000000..e2e2f47 --- /dev/null +++ b/test/eos_equations.jl @@ -0,0 +1,133 @@ +const MCF = MultiComponentFlash + +@testset "Cubic EOS defining equations" begin + property = MolecularProperty(0.050, 5.0e6, 400.0, 1.0e-4, 0.80) + mixture = MultiComponentMixture((property,); names = ["test component"]) + cond = (p = 2.0e6, T = 320.0, z = [1.0]) + Tr = cond.T/property.T_c + + @testset "Static and temperature coefficients" begin + pr = GenericCubicEOS(mixture, PengRobinson()) + prc = GenericCubicEOS(mixture, PengRobinsonCorrected()) + srk = GenericCubicEOS(mixture, SoaveRedlichKwong()) + rk = GenericCubicEOS(mixture, RedlichKwong()) + + @test collect(MCF.static_coefficients(PengRobinson())) ≈ + [0.457235529, 0.077796074, 1 + sqrt(2), 1 - sqrt(2)] + @test collect(MCF.static_coefficients(SoaveRedlichKwong())) ≈ + [0.4274802327, 0.08664035, 0.0, 1.0] + + kappa_pr = 0.37464 + 1.54226property.ω - 0.26992property.ω^2 + alpha_pr = (1 + kappa_pr*(1 - sqrt(Tr)))^2 + @test MCF.weight_ai(pr, cond, 1) ≈ 0.457235529*alpha_pr + + kappa_pr78 = 0.379642 + 1.48503property.ω - + 0.164423property.ω^2 + 0.016666property.ω^3 + alpha_pr78 = (1 + kappa_pr78*(1 - sqrt(Tr)))^2 + @test MCF.weight_ai(prc, cond, 1) ≈ 0.457235529*alpha_pr78 + + kappa_srk = 0.48 + 1.574property.ω - 0.176property.ω^2 + alpha_srk = (1 + kappa_srk*(1 - sqrt(Tr)))^2 + @test MCF.weight_ai(srk, cond, 1) ≈ 0.4274802327*alpha_srk + @test MCF.weight_ai(rk, cond, 1) ≈ 0.4274802327/sqrt(Tr) + end + + @testset "Zudkevitch-Joffe modifiers" begin + Fa(T, i) = 1 + 1.0e-3T + 0.1i + Fb(T, i) = 0.8 + 5.0e-4T + 0.05i + zj = GenericCubicEOS(mixture, ZudkevitchJoffe(; F_a = Fa, F_b = Fb)) + @test MCF.weight_ai(zj, cond, 1) ≈ + 0.4274802327*Fa(cond.T, 1)/sqrt(Tr) + @test MCF.weight_bi(zj, cond, 1) ≈ 0.08664035*Fb(cond.T, 1) + + rk = GenericCubicEOS(mixture, RedlichKwong()) + zj_default = GenericCubicEOS(mixture, ZudkevitchJoffe()) + @test MCF.weight_ai(zj_default, cond, 1) == MCF.weight_ai(rk, cond, 1) + @test MCF.weight_bi(zj_default, cond, 1) == MCF.weight_bi(rk, cond, 1) + end + + @testset "Generalized cubic polynomial" begin + for eos_type in (PengRobinson(), SoaveRedlichKwong(), RedlichKwong()) + eos = GenericCubicEOS(mixture, eos_type) + A = 0.31 + B = 0.047 + polynomial = MCF.cubic_polynomial(eos, A, B) + for Z in (0.09, 0.8, 1.3) + cubic_residual = Z^3 + polynomial[1]*Z^2 + + polynomial[2]*Z + polynomial[3] + eos_residual = 1 - 1/(Z - B) + + A/((Z + eos.m_1*B)*(Z + eos.m_2*B)) + denominator = (Z - B)*(Z + eos.m_1*B)*(Z + eos.m_2*B) + @test cubic_residual ≈ eos_residual*denominator + end + end + end +end + +@testset "Soreide-Whitson correlations" begin + molality = 2.5 + sw = SoreideWhitson(["Water", "generic hydrocarbon"]; molality = molality) + acentric = 0.35 + Tr = 1.20 + + A0 = sw.A[1] + sw.A_mw[1]*sign(acentric)*abs(acentric)^(-0.1) + A1 = sw.A[2] + sw.A_mw[2]*acentric + A2 = sw.A[3] + sw.A_mw[3]*acentric + expected_hc_bic = A0*(1 + sw.alphas[1]*molality) + + A1*Tr*(1 + sw.alphas[2]*molality) + + A2*Tr^2*(1 + sw.alphas[3]*molality) + @test MCF.soreide_whitson_hc_aqueous_bic(sw, acentric, Tr) ≈ + expected_hc_bic + + expected_n2_bic = -1.70235*(1 + 0.025587molality^0.75) + + 0.44338*(1 + 0.08126molality^0.75)*Tr + expected_co2_bic = -0.31092*(1 + 0.15587molality^0.7505) + + 0.23580*(1 + 0.17837molality^0.979)*Tr - + 21.2566exp(-6.7222Tr - molality) + @test MCF.soreide_whitson_n2_aqueous_bic(sw, acentric, Tr) ≈ expected_n2_bic + @test MCF.soreide_whitson_h2s_aqueous_bic(sw, acentric, Tr) ≈ + -0.20441 + 0.234267Tr + @test MCF.soreide_whitson_co2_aqueous_bic(sw, acentric, Tr) ≈ + expected_co2_bic + + water = MolecularProperty("Water") + hydrocarbon = MolecularProperty(0.050, 5.0e6, 400.0, 1.0e-4, acentric) + interaction = [0.0 0.12; 0.12 0.0] + mixture = MultiComponentMixture((water, hydrocarbon); + names = ["Water", "generic hydrocarbon"], A_ij = interaction) + sw = SoreideWhitson(mixture; molality = molality) + eos = GenericCubicEOS(mixture, sw) + cond = (p = 1.0e6, T = 350.0, z = [0.4, 0.6]) + water_Tr = cond.T/water.T_c + alpha_half = 1 + 0.4530*(1 - (1 - 0.0103molality^1.1)*water_Tr) + + 0.0034*(water_Tr^(-3) - 1) + @test MCF.weight_ai(eos, cond, 1) ≈ 0.457235529*alpha_half^2 + + liquid = MCF.set_phase(cond, :liquid) + vapor = MCF.set_phase(cond, :vapor) + expected = MCF.soreide_whitson_hc_aqueous_bic( + sw, acentric, cond.T/hydrocarbon.T_c) + @test MCF.binary_interaction(eos, 1, 2, liquid) ≈ expected + @test MCF.binary_interaction(eos, 2, 1, liquid) ≈ expected + @test MCF.binary_interaction(eos, 1, 2, vapor) == 0.12 + + storage = flash_storage(eos, cond) + @test MCF.get_force_coefficients(storage.forces, eos, liquid) === + storage.forces.liquid + @test MCF.get_force_coefficients(storage.forces, eos, vapor) === + storage.forces.vapor + vapor_fraction, K, report = flash_2ph(eos, cond; extra_out = true) + @test report.converged + @test isfinite(vapor_fraction) + @test all(isfinite, K) +end + +@testset "Cubic root degeneracies" begin + @test MCF.solve_cubic_positive_roots(-3.0, 3.0, -1.0) ≈ 1.0 + @test MCF.solve_cubic_positive_roots(0.0, 0.0, 0.0) ≈ 0.0 + @test sort(collect(MCF.solve_cubic_positive_roots(-4.0, 5.0, -2.0))) ≈ + [1.0, 1.0, 2.0] + + root = MCF.solve_cubic_positive_roots(0.0, 1.0, 1.0) + @test root^3 + root + 1 ≈ 0.0 atol = 1.0e-14 +end diff --git a/test/runtests.jl b/test/runtests.jl index 1e17c26..32047cf 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,5 @@ include("test_setup.jl") +include("eos_equations.jl") @testset "Rachford-Rice" begin test_rachford_rice()