From 30a1333a8a11ffd4529115249e00ac8ae647d725 Mon Sep 17 00:00:00 2001 From: pw0908 Date: Thu, 1 Aug 2024 15:56:45 -0700 Subject: [PATCH 01/18] Add definition of activities --- src/methods/pT.jl | 25 ++++++++++++++++++++++++- src/models/Activity/equations.jl | 6 ++++++ 2 files changed, 30 insertions(+), 1 deletion(-) diff --git a/src/methods/pT.jl b/src/methods/pT.jl index 8c7248abe..0a5af8d46 100755 --- a/src/methods/pT.jl +++ b/src/methods/pT.jl @@ -455,6 +455,29 @@ function activity_coefficient(model::EoSModel,p,T,z=SA[1.]; phase=:unknown, thre return exp.((μ_mixt .- μ_pure) ./ R̄ ./ T) ./z end + +""" + activity(model::EoSModel,p,T,z=SA[1.0]; phase=:unknown, threaded=true, vol0=nothing) + +Calculates the activity, defined as: +```julia +log(a) = (μ_mixt - μ_pure) / R̄ / T +``` +where `μ_mixt` is the chemical potential of the mixture and `μ_pure` is the chemical potential of the pure components. + +Internally, it calls [`Clapeyron.volume`](@ref) to obtain `V` and +calculates the property via `VT_fugacity_coefficient(model,V,T,z)`. + +The keywords `phase`, `threaded` and `vol0` are passed to the [`Clapeyron.volume`](@ref) solver. +""" +function activity(model::EoSModel,p,T,z=SA[1.]; phase=:unknown, threaded=true, vol0=nothing) + pure = split_model(model) + μ_mixt = chemical_potential(model, p, T, z; phase, threaded, vol0) + μ_pure = gibbs_free_energy.(pure, p, T; phase, threaded, vol0) + R̄ = Rgas(model) + return exp.((μ_mixt .- μ_pure) ./ R̄ ./ T) +end + """ compressibility_factor(model::EoSModel, p, T, z=SA[1.]; phase=:unknown, threaded=true, vol0=nothing) @@ -618,6 +641,6 @@ export fundamental_derivative_of_gas_dynamics #volume properties export mass_density,molar_density, compressibility_factor #molar gradient properties -export chemical_potential, activity_coefficient, fugacity_coefficient +export chemical_potential, activity_coefficient, activity, fugacity_coefficient export chemical_potential_res export mixing, excess, gibbs_solvation \ No newline at end of file diff --git a/src/models/Activity/equations.jl b/src/models/Activity/equations.jl index f12da9dc8..d75d6055a 100755 --- a/src/models/Activity/equations.jl +++ b/src/models/Activity/equations.jl @@ -27,6 +27,12 @@ function activity_coefficient(model::ActivityModel,p,T,z) return exp.(Solvers.gradient(x->excess_gibbs_free_energy(model,p,T,x),z)/(R̄*T))::X end +function activity(model::ActivityModel,p,T,z=SA[1.]) + γ = activity_coefficient(model, p, T, z) + x = z ./ sum(z) + return γ .* x +end + function test_activity_coefficient(model::ActivityModel,p,T,z) X = gradient_type(p,T,z) return exp.(Solvers.gradient(x->excess_gibbs_free_energy(model,p,T,x),z)/(R̄*T))::X From 868eb5e01dcef230f9346d31f31689594933170d Mon Sep 17 00:00:00 2001 From: pw0908 Date: Thu, 1 Aug 2024 15:57:01 -0700 Subject: [PATCH 02/18] Begin work on reactive systems --- src/Clapeyron.jl | 6 +- .../property_solvers/reactive/reactive.jl | 65 ++++++++++++++++++ src/models/Reactive/reactive.jl | 67 +++++++++++++++++++ 3 files changed, 136 insertions(+), 2 deletions(-) create mode 100644 src/methods/property_solvers/reactive/reactive.jl create mode 100644 src/models/Reactive/reactive.jl diff --git a/src/Clapeyron.jl b/src/Clapeyron.jl index f878facfb..d7cd8dc2e 100755 --- a/src/Clapeyron.jl +++ b/src/Clapeyron.jl @@ -258,12 +258,14 @@ include("models/SAFT/SAFTVRMie/variants/SAFTVREMie.jl") include("models/SAFT/SAFTVRMie/variants/eSAFTVRMie.jl") include("models/SAFT/PCSAFT/variants/ePCSAFT.jl") # include("models/Electrolytes/ElectrolyteSAFT/eCPA.jl") - - include("methods/property_solvers/electrolytes/electrolytes.jl") include("methods/property_solvers/multicomponent/tp_flash/electrolyte_flash.jl") include("models/AnalyticalSLV/AnalyticalSLV.jl") +# Export reactive systems +include("models/Reactive/reactive.jl") +include("methods/property_solvers/reactive/reactive.jl") + #Unitful support, transition from dependency to ext if !isdefined(Base,:get_extension) include("../ext/ClapeyronUnitfulExt.jl") diff --git a/src/methods/property_solvers/reactive/reactive.jl b/src/methods/property_solvers/reactive/reactive.jl new file mode 100644 index 000000000..dbfce89fc --- /dev/null +++ b/src/methods/property_solvers/reactive/reactive.jl @@ -0,0 +1,65 @@ +# function pH(model::ReactiveModel,p,T,n0;z0=nothing) +# ν = stoichiometric_coefficient(model) # Exists +# ΔHf = model.params.ΔHf.values +# Sf = model.params.Sf.values +# ΔGf = ΔHf-T*Sf +# ΔrG = sum(ΔGf*ν,dims=1) +# Keq = exp.(-ΔrG/Rgas(model.model)/T) + +# if isnothing(z0) +# z0 = x0_equilibrium_conditions(model,p,T,n0) +# end + +# f!(F,z) = equilibrium_conditions(model,F,p,T,n0,z,ν,Keq) +# sol = Solvers.nlsolve(f!,z0,TrustRegion(Newton(), NWI()),NEqOptions(),ForwardDiff.Chunk{1}()) # Exists +# ξ = Solvers.x_sol(sol) +# n = n0+ν*ξ + +# a = activity(model.model,p,T,n) # Exists + +# hydronium_id = find_hydronium_index(model) + +# return -log10(a[hydronium_id]) +# end + +# function find_hydronium_index(model) +# idx = findfirst(model.components.=="hydronium") +# return idx +# end + +function equilibrate(model::ReactiveModel,p,T,n0;z0=nothing) + ν = stoichiometric_coefficient(model) # Exists + ΔHf = model.params.ΔHf.values + Sf = model.params.Sf.values + ΔGf = ΔHf-T*Sf + ΔrG = sum(ΔGf.*ν,dims=1) + Keq = exp.(-ΔrG/Rgas(model.eosmodel)/T) + + if isnothing(z0) + z0 = x0_equilibrium_conditions(model,p,T,n0) + end + + f!(F,z) = equilibrium_conditions(model,F,p,T,n0,z,ν,Keq) + sol = Solvers.nlsolve(f!,z0,TrustRegion(Newton(), NWI()),NEqOptions(),ForwardDiff.Chunk{1}()) # Exists + ξ = Solvers.x_sol(sol) + n = n0+ν*ξ + return n +end + +function equilibrium_conditions(model::ReactiveModel,F,p,T,n0,ξ,ν,Keq) + n = n0+ν*ξ + + a = activity(model.eosmodel,p,T,n) + + F[1:end] = sum(log.(a).*ν,dims=1) + log.(Keq) + return F +end + +function x0_equilibrium_conditions(model::ReactiveModel,p,T,n0) + + z0 = zeros(length(model.components)) + z0[1] = log10(p) + return z0 +end + +export equilibrate \ No newline at end of file diff --git a/src/models/Reactive/reactive.jl b/src/models/Reactive/reactive.jl new file mode 100644 index 000000000..3ca3a63c9 --- /dev/null +++ b/src/models/Reactive/reactive.jl @@ -0,0 +1,67 @@ +struct ReactiveParams{T} <: EoSParam where T <: Float64 + ΔHf::SingleParam{T} + Sf::SingleParam{T} +end + +struct ReactiveModel{T} + components::Vector{String} + reactions::Vector + params::ReactiveParams{Float64} + eosmodel::T +end + +function ReactiveModel(components::Vector{String}, reactions::Vector; + model::EoSModel, + userlocations=String[], + verbose=false) + params = getparams(components, ["properties/formation.csv"]; userlocations = userlocations, verbose = verbose) + ΔHf = params["Hf"] + Sf = params["Sf"] + + packagedparams = ReactiveParams(ΔHf,Sf) + + return ReactiveModel(components,reactions,packagedparams,model) +end + +function stoichiometric_coefficient(model::ReactiveModel) + ν = zeros(length(model.components),length(model.reactions)) + for (i,reaction) in enumerate(model.reactions) + species = [reaction[k][1] for k in 1:length(reaction)] + coeff = [reaction[k][2] for k in 1:length(reaction)] + for (j,component) in enumerate(model.components) + if component in species + ν[j,i] = coeff[findfirst(isequal(component),species)] + end + end + end + return ν +end + +function Base.show(io::IO, mime::MIME"text/plain", model::ReactiveModel) + print(io, typeof(model)) + length(model.eosmodel) == 1 && println(io, " with 1 component:") + length(model.eosmodel) > 1 && println(io, " with ", length(model.eosmodel), " components:") + show_pairs(io,model.components) + println(io, "") + + length(model.reactions) == 1 && println(io, "and 1 reaction:") + length(model.reactions) > 1 && println(io, "and ", length(model.reactions), " reactions:") + for reaction in model.reactions + species = [reaction[k][1] for k in 1:length(reaction)] + coeff = [reaction[k][2] for k in 1:length(reaction)] + reactants = "" + products = "" + for i in 1:length(species) + if coeff[i] < 0 + reactants *= string(abs(coeff[i]))*" "*species[i]*" + " + else + products *= string(abs(coeff[i]))*" "*species[i]*" + " + end + end + reactants = reactants[1:end-3] + products = products[1:end-3] + println(io, reactants, " <=> ", products) + end +end + +export ReactiveModel \ No newline at end of file From 48c926f69c1496cdb3e2176c24b8d3d8d2885014 Mon Sep 17 00:00:00 2001 From: pw0908 Date: Mon, 12 Aug 2024 19:28:14 -0700 Subject: [PATCH 03/18] Add database --- database/properties/formation.csv | 7 +++++++ database/properties/formation_aqueous.csv | 9 +++++++++ 2 files changed, 16 insertions(+) create mode 100755 database/properties/formation.csv create mode 100755 database/properties/formation_aqueous.csv diff --git a/database/properties/formation.csv b/database/properties/formation.csv new file mode 100755 index 000000000..ffd9b4775 --- /dev/null +++ b/database/properties/formation.csv @@ -0,0 +1,7 @@ +Clapeyron Database File,, +Enthalpies and Entropies of formation Single Params,, +species,Hf,Sf +ethanol,-234.7e3,283.0 +water,-241.83e3,188.84 +acetic acid,-438.1e3,282.84 +ethyl acetate,-445.0e3,363.0 \ No newline at end of file diff --git a/database/properties/formation_aqueous.csv b/database/properties/formation_aqueous.csv new file mode 100755 index 000000000..a839adf74 --- /dev/null +++ b/database/properties/formation_aqueous.csv @@ -0,0 +1,9 @@ +Clapeyron Database File,, +Enthalpies and Gibbs free energy of formation in aqueous system Single Params,, +species,Hf,Gf +hydroxide,-229.994e3,-157.244e3 +hydronium,-157.244e3,-237.129e3 +water,-157.244e3,-237.129e3 +acetate,-486.01e3,-369.31e3 +acetic acid,-485.76e3,-396.46e3 +sodium,-240.12e3,-261.905e3 \ No newline at end of file From 6b9e18b6ea4e1fbbe32b960eb31ac65d368f178c Mon Sep 17 00:00:00 2001 From: pw0908 Date: Mon, 12 Aug 2024 19:28:24 -0700 Subject: [PATCH 04/18] Minor typo --- src/methods/property_solvers/reactive/reactive.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/methods/property_solvers/reactive/reactive.jl b/src/methods/property_solvers/reactive/reactive.jl index dbfce89fc..fe82e4b88 100644 --- a/src/methods/property_solvers/reactive/reactive.jl +++ b/src/methods/property_solvers/reactive/reactive.jl @@ -51,7 +51,7 @@ function equilibrium_conditions(model::ReactiveModel,F,p,T,n0,ξ,ν,Keq) a = activity(model.eosmodel,p,T,n) - F[1:end] = sum(log.(a).*ν,dims=1) + log.(Keq) + F[1:end] = sum(log.(a).*ν,dims=1) - log.(Keq) return F end From e7b4161a7ae8afb1e0d87bf2ab533626f92951a1 Mon Sep 17 00:00:00 2001 From: pw0908 Date: Mon, 12 Aug 2024 19:28:43 -0700 Subject: [PATCH 05/18] Add aqueous reactive systems --- src/Clapeyron.jl | 2 + src/methods/pT.jl | 39 ++++++++++- .../property_solvers/reactive/reactive_aq.jl | 51 ++++++++++++++ src/models/Reactive/reactive_aq.jl | 68 +++++++++++++++++++ 4 files changed, 159 insertions(+), 1 deletion(-) create mode 100644 src/methods/property_solvers/reactive/reactive_aq.jl create mode 100644 src/models/Reactive/reactive_aq.jl diff --git a/src/Clapeyron.jl b/src/Clapeyron.jl index d7cd8dc2e..27c396eeb 100755 --- a/src/Clapeyron.jl +++ b/src/Clapeyron.jl @@ -264,7 +264,9 @@ include("models/AnalyticalSLV/AnalyticalSLV.jl") # Export reactive systems include("models/Reactive/reactive.jl") +include("models/Reactive/reactive_aq.jl") include("methods/property_solvers/reactive/reactive.jl") +include("methods/property_solvers/reactive/reactive_aq.jl") #Unitful support, transition from dependency to ext if !isdefined(Base,:get_extension) diff --git a/src/methods/pT.jl b/src/methods/pT.jl index 0a5af8d46..76918b78d 100755 --- a/src/methods/pT.jl +++ b/src/methods/pT.jl @@ -478,6 +478,43 @@ function activity(model::EoSModel,p,T,z=SA[1.]; phase=:unknown, threaded=true, v return exp.((μ_mixt .- μ_pure) ./ R̄ ./ T) end +function find_hydronium_index(model) + idx = findfirst(model.components.=="hydronium") + return idx +end + +function find_water_indx(model) + idx = findfirst(model.components.=="water") + return idx +end + +""" + aqueous_activity(model::EoSModel,p,T,z=SA[1.0]; phase=:unknown, threaded=true, vol0=nothing) + +Calculates the activity with the reference being infinite dilution in water, defined as: +```julia +log(a) = (μ_mixt - μ_inf) / R̄ / T +``` +where `μ_mixt` is the chemical potential of the mixture and `μ_inf` is the chemical potential of the components at infinite dilution in water. + +Internally, it calls [`Clapeyron.volume`](@ref) to obtain `V` and +calculates the property via `VT_fugacity_coefficient(model,V,T,z)`. + +The keywords `phase`, `threaded` and `vol0` are passed to the [`Clapeyron.volume`](@ref) solver. +""" +function aqueous_activity(model::EoSModel,p,T,z=SA[1.]; phase=:unknown, threaded=true, vol0=nothing) + μ_mixt = chemical_potential(model, p, T, z; phase, threaded, vol0) + + idx_w = find_water_indx(model) + zref = ones(length(model)) + zref[1:length(model) .!= idx_w] .*= 0.01801528 + zref ./= sum(zref) + + μ_ref = chemical_potential(model, p, T, zref; phase, threaded, vol0) + R̄ = Rgas(model) + return exp.((μ_mixt .- μ_ref) ./ R̄ ./ T) +end + """ compressibility_factor(model::EoSModel, p, T, z=SA[1.]; phase=:unknown, threaded=true, vol0=nothing) @@ -641,6 +678,6 @@ export fundamental_derivative_of_gas_dynamics #volume properties export mass_density,molar_density, compressibility_factor #molar gradient properties -export chemical_potential, activity_coefficient, activity, fugacity_coefficient +export chemical_potential, activity_coefficient, activity, aqueous_activity, fugacity_coefficient export chemical_potential_res export mixing, excess, gibbs_solvation \ No newline at end of file diff --git a/src/methods/property_solvers/reactive/reactive_aq.jl b/src/methods/property_solvers/reactive/reactive_aq.jl new file mode 100644 index 000000000..2dc063ad0 --- /dev/null +++ b/src/methods/property_solvers/reactive/reactive_aq.jl @@ -0,0 +1,51 @@ +function pH(model::ReactiveAqModel,p,T,n0;z0=nothing) + n = equilibrate(model, p, T, n0; z0 = z0) + + a = aqueous_activity(model.eosmodel,p,T,n) + + hydronium_id = find_hydronium_index(model) + + return -log10(a[hydronium_id]) +end + +function equilibrate(model::ReactiveAqModel,p,T,n0;z0=nothing) + T0 = 298.15 + ν = stoichiometric_coefficient(model) + ΔGf0 = model.params.ΔHf.values + ΔHf = model.params.ΔHf.values + ΔGf = ΔGf0/T0+ΔHf*(1/T-1/T0) + ΔrG = sum(ΔGf.*ν,dims=1) + Keq = exp.(-ΔrG/Rgas(model.eosmodel)) + + # println(Keq) + + if isnothing(z0) + z0 = x0_equilibrium_conditions(model,p,T,n0) + end + + f!(F,z) = equilibrium_conditions(model,F,p,T,n0,exp10.(z),ν,Keq) + sol = Solvers.nlsolve(f!,log10.(z0),TrustRegion(Newton(), NWI()),NEqOptions(),ForwardDiff.Chunk{1}()) # Exists + ξ = exp10.(Solvers.x_sol(sol)) + n = n0+ν*ξ + return n +end + +function equilibrium_conditions(model::ReactiveAqModel,F,p,T,n0,ξ,ν,Keq) + n = n0+ν*ξ + # println(n) + + a = aqueous_activity(model.eosmodel,p,T,n) + + F[1:end] = sum(log.(a).*ν,dims=1) - log.(Keq) + # println(F) + return F +end + +function x0_equilibrium_conditions(model::ReactiveAqModel,p,T,n0) + + z0 = zeros(length(model.components)) + z0[1] = log10(p) + return z0 +end + +export equilibrate \ No newline at end of file diff --git a/src/models/Reactive/reactive_aq.jl b/src/models/Reactive/reactive_aq.jl new file mode 100644 index 000000000..6716dc567 --- /dev/null +++ b/src/models/Reactive/reactive_aq.jl @@ -0,0 +1,68 @@ +struct ReactiveAqParams{T} <: EoSParam where T <: Float64 + ΔGf::SingleParam{T} + ΔHf::SingleParam{T} +end + +struct ReactiveAqModel{T} + components::Vector{String} + reactions::Vector + params::ReactiveParams{Float64} + eosmodel::T +end + +function ReactiveAqModel(components::Vector{String}, reactions::Vector; + model::EoSModel, + userlocations=String[], + verbose=false) + params = getparams(components, ["properties/formation_aqueous.csv"]; userlocations = userlocations, verbose = verbose) + ΔGf = params["Gf"] + ΔHf = params["Hf"] + + packagedparams = ReactiveParams(ΔGf,ΔHf) + + return ReactiveAqModel(components,reactions,packagedparams,model) +end + +function stoichiometric_coefficient(model::ReactiveAqModel) + ν = zeros(length(model.components),length(model.reactions)) + for (i,reaction) in enumerate(model.reactions) + species = [reaction[k][1] for k in 1:length(reaction)] + coeff = [reaction[k][2] for k in 1:length(reaction)] + for (j,component) in enumerate(model.components) + if component in species + ν[j,i] = coeff[findfirst(isequal(component),species)] + end + end + end + return ν +end + + +function Base.show(io::IO, mime::MIME"text/plain", model::ReactiveAqModel) + print(io, typeof(model)) + length(model.eosmodel) == 1 && println(io, " with 1 component:") + length(model.eosmodel) > 1 && println(io, " with ", length(model.eosmodel), " components:") + show_pairs(io,model.components) + println(io, "") + + length(model.reactions) == 1 && println(io, "and 1 reaction:") + length(model.reactions) > 1 && println(io, "and ", length(model.reactions), " reactions:") + for reaction in model.reactions + species = [reaction[k][1] for k in 1:length(reaction)] + coeff = [reaction[k][2] for k in 1:length(reaction)] + reactants = "" + products = "" + for i in 1:length(species) + if coeff[i] < 0 + reactants *= string(abs(coeff[i]))*" "*species[i]*" + " + else + products *= string(abs(coeff[i]))*" "*species[i]*" + " + end + end + reactants = reactants[1:end-3] + products = products[1:end-3] + println(io, reactants, " <=> ", products) + end +end + +export ReactiveAqModel \ No newline at end of file From 0b2824c0cc70ee62ee625900d6139a0912bcc2d6 Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Wed, 14 Aug 2024 01:24:25 -0400 Subject: [PATCH 06/18] reorder reactive, define reference_chemical_potential --- docs/src/properties/bulk.md | 13 +- src/methods/pT.jl | 146 ++++++++++++++---- .../property_solvers/reactive/reactive.jl | 58 ++----- .../property_solvers/reactive/reactive_aq.jl | 45 ------ src/models/Reactive/reactive.jl | 29 +++- src/models/Reactive/reactive_aq.jl | 62 ++------ 6 files changed, 174 insertions(+), 179 deletions(-) diff --git a/docs/src/properties/bulk.md b/docs/src/properties/bulk.md index c668f707e..36b59681d 100755 --- a/docs/src/properties/bulk.md +++ b/docs/src/properties/bulk.md @@ -76,6 +76,15 @@ Clapeyron.chemical_potential_res Clapeyron.fugacity_coefficient ``` +### Activity Coefficient +```@docs +Clapeyron.reference_chemical_potential +Clapeyron.reference_chemical_potential_type +Clapeyron.activity_coefficient +Clapeyron.activity +Clapeyron.aqueous_activity +``` + ### Mixing ```@docs Clapeyron.mixing @@ -83,6 +92,8 @@ Clapeyron.mixing ## Initial guess functions +These methods are considered internal, they don't support `Symbolics.jl` or `Unitful.jl` overloads. + ```@docs Clapeyron.lb_volume Clapeyron.T_scale @@ -97,4 +108,4 @@ Clapeyron.x0_psat Clapeyron.x0_saturation_temperature Clapeyron.antoine_coef Clapeyron.x0_crit_pure -``` \ No newline at end of file +``` diff --git a/src/methods/pT.jl b/src/methods/pT.jl index 76918b78d..7f01d7b1d 100755 --- a/src/methods/pT.jl +++ b/src/methods/pT.jl @@ -397,7 +397,7 @@ Returns the phase of a fluid at the conditions specified by `V`, `T` and `z`. Uses the phase identification parameter criteria from `Clapeyron.pip` returns `:liquid` if the phase is liquid (or liquid-like), `:vapour` if the phase is vapour (or vapour-like), and `:unknown` if the calculation of the phase identification parameter failed. - + Internally, it calls [`Clapeyron.volume`](@ref) to obtain `V` and calculates the property via `VT_enthalpy(model,V,T,z)`. The keywords `phase`, `threaded` and `vol0` are passed to the [`Clapeyron.volume`](@ref) solver. @@ -411,7 +411,7 @@ end fundamental_derivative_of_gas_dynamics(model::EoSModel, p, T, z=SA[1.]; phase=:gas, threaded=true, vol0=nothing)::Symbol Calculates the fundamental derivative of gas dynamics. - + Internally, it calls [`Clapeyron.volume`](@ref) to obtain `V` and calculates the property via `VT_enthalpy(model,V,T,z)`. The keywords `phase`, `threaded` and `vol0` are passed to the [`Clapeyron.volume`](@ref) solver. @@ -447,48 +447,91 @@ function fugacity_coefficient!(φ,model::EoSModel,p,T,z=SA[1.]; phase=:unknown, VT_fugacity_coefficient!(φ,model,V,T,z,p) end -function activity_coefficient(model::EoSModel,p,T,z=SA[1.]; phase=:unknown, threaded=true, vol0=nothing) - pure = split_model(model) - μ_mixt = chemical_potential(model, p, T, z; phase, threaded, vol0) - μ_pure = gibbs_free_energy.(pure, p, T; phase, threaded, vol0) - R̄ = Rgas(model) - return exp.((μ_mixt .- μ_pure) ./ R̄ ./ T) ./z +""" + activity_coefficient(model::EoSModel,p,T,z=SA[1.0];reference = :pure, phase=:unknown, threaded=true, vol0=nothing) + +Calculates the activity, defined as: +```julia +log(γ*z) = (μ_mixt - μ_ref) / R̄ / T +``` +where `μ_mixt` is the chemical potential of the mixture and `μ_ref` is the reference chemical potential for the model at `p`,`T` conditions, calculated via [`Clapeyron.reference_chemical_potential`](@ref). +Internally, it calls [`Clapeyron.volume`](@ref) to obtain `V` and +calculates the property via `VT_fugacity_coefficient(model,V,T,z)`. + +The keywords `phase`, `threaded` and `vol0` are passed to the [`Clapeyron.volume`](@ref) solver. + +If the `μ_ref` keyword argument is not provided, the `reference` keyword is used to specify the reference chemical potential.. +""" +function activity_coefficient(model::EoSModel,p,T,z=SA[1.]; + μ_ref = nothing, + reference = :pure, + phase=:unknown, + threaded=true, + vol0=nothing) + if model isa ActivityModel + return activity_coefficient(model,p,T,z) + end + if μ_ref == nothing + return activity_coefficient_impl(model,p,T,z,reference_chemical_potential(model,p,T,reference;phase,threaded),reference,phase,threaded,vol0) + else + return activity_coefficient_impl(model,p,T,z,μ_ref,reference,phase,threaded,vol0) + end end +function activity_coefficient_impl(model,p,T,z,μ_ref,reference,phase,threaded,vol0) + R̄ = Rgas(model) + μ_mixt = chemical_potential(model, p, T, z; phase, threaded, vol0) + return exp.((μ_mixt .- μ_ref) ./ R̄ ./ T) ./z +end -""" - activity(model::EoSModel,p,T,z=SA[1.0]; phase=:unknown, threaded=true, vol0=nothing) +""" + activity(model::EoSModel,p,T,z=SA[1.0];reference = :pure, phase=:unknown, threaded=true, vol0=nothing) Calculates the activity, defined as: ```julia -log(a) = (μ_mixt - μ_pure) / R̄ / T +log(a) = (μ_mixt - μ_ref) / R̄ / T ``` -where `μ_mixt` is the chemical potential of the mixture and `μ_pure` is the chemical potential of the pure components. - +where `μ_mixt` is the chemical potential of the mixture and `μ_ref` is the reference chemical potential for the model at `p`,`T` conditions, calculated via [`Clapeyron.reference_chemical_potential`](@ref). Internally, it calls [`Clapeyron.volume`](@ref) to obtain `V` and calculates the property via `VT_fugacity_coefficient(model,V,T,z)`. The keywords `phase`, `threaded` and `vol0` are passed to the [`Clapeyron.volume`](@ref) solver. + +If the `μ_ref` keyword argument is not provided, the `reference` keyword is used to specify the reference chemical potential.. """ -function activity(model::EoSModel,p,T,z=SA[1.]; phase=:unknown, threaded=true, vol0=nothing) - pure = split_model(model) - μ_mixt = chemical_potential(model, p, T, z; phase, threaded, vol0) - μ_pure = gibbs_free_energy.(pure, p, T; phase, threaded, vol0) +function activity(model::EoSModel,p,T,z=SA[1.]; + μ_ref = nothing, + reference = :pure, + phase=:unknown, + threaded=true, + vol0=nothing) + + if μ_ref == nothing + return activity_impl(model,p,T,z,reference_chemical_potential(model,p,T,reference,phase;threaded),reference,phase,threaded,vol0) + else + return activity_impl(model,p,T,z,μ_ref,reference,phase,threaded,vol0) + end +end + +function activity_impl(model,p,T,z,μ_ref,reference,phase,threaded,vol0) R̄ = Rgas(model) - return exp.((μ_mixt .- μ_pure) ./ R̄ ./ T) + μ_mixt = chemical_potential(model, p, T, z; phase, threaded, vol0) + return exp.((μ_mixt .- μ_ref) ./ R̄ ./ T) end function find_hydronium_index(model) - idx = findfirst(model.components.=="hydronium") + idx = findfirst(isequal("hydronium"),components) + idx == nothing && return 0 return idx end function find_water_indx(model) - idx = findfirst(model.components.=="water") + idx = findfirst(isequal("water"),components) + idx == nothing && return 0 return idx end -""" +""" aqueous_activity(model::EoSModel,p,T,z=SA[1.0]; phase=:unknown, threaded=true, vol0=nothing) Calculates the activity with the reference being infinite dilution in water, defined as: @@ -502,17 +545,58 @@ calculates the property via `VT_fugacity_coefficient(model,V,T,z)`. The keywords `phase`, `threaded` and `vol0` are passed to the [`Clapeyron.volume`](@ref) solver. """ -function aqueous_activity(model::EoSModel,p,T,z=SA[1.]; phase=:unknown, threaded=true, vol0=nothing) - μ_mixt = chemical_potential(model, p, T, z; phase, threaded, vol0) +function aqueous_activity(model::EoSModel,p,T,z=SA[1.]; + μ_ref = nothing, + reference = :aqueous, + phase=:unknown, + threaded=true, + vol0=nothing) + return activity(model,p,T,z;reference=reference,μ_ref = μ_ref,phase=phase,threaded=threaded,vol0=vol0) +end - idx_w = find_water_indx(model) - zref = ones(length(model)) - zref[1:length(model) .!= idx_w] .*= 0.01801528 - zref ./= sum(zref) - μ_ref = chemical_potential(model, p, T, zref; phase, threaded, vol0) - R̄ = Rgas(model) - return exp.((μ_mixt .- μ_ref) ./ R̄ ./ T) +""" + reference_chemical_potential_type(model)::Symbol + +Returns a symbol with the type of reference chemical potential used by the input `model`. + +""" +reference_chemical_potential_type(model) = :pure + +""" + reference_chemical_potential(model::EoSModel,p,T,reference; phase=:unknown, threaded=true, vol0=nothing) + +Returns a reference chemical potential. used in calculation of `activity` and actitivy_coefficient. there are two available references: +- `:pure`: the reference potential is a pure component at specified `T`, `p` and `phase` +- `:aqueous`: the chemical potential of the pure components at specified `T`, `p` and `phase` +- `:sat_pure_T`: the reference potential is the pure saturated liquid phase at specified `T`. +- `:zero`: the reference potential is equal to zero for all components (used for `ActivityModel`) +The keywords `phase`, `threaded` and `vol0` are passed to the [`Clapeyron.volume`](@ref) solver. +""" +function reference_chemical_potential(model::EoSModel,p,T,reference = reference_chemical_potential_type(model), phase=:unknown, threaded=true, vol0=nothing) + if reference == :pure + pure = split_model.(model) + return gibbs_free_energy.(pure, p, T; phase, threaded) + elseif reference == :aqueous + idx_w = find_water_indx(model) + if idx_w == 0 + throw(ArgumentError("There is no water in $(model).")) + end + zref = ones(length(model)) + zref[1:length(model) .!= idx_w] .*= 0.01801528 + zref ./= sum(zref) + return chemical_potential(model, p, T, zref; phase, threaded, vol0) + elseif reference == :sat_pure_T + pure = split_model.(model) + sat = saturation_pressure.(pure,T) + vl_pure = getindex.(sat,2) + return VT_gibbs_free_energy.(pure, vl_pure, T) + elseif reference == :zero + _0 = Base.promote_eltype(model,p,T) + return fill(_0,length(model)) + else + throw(ArgumentError("reference must be one of :pure, :aqueous, :sat_pure_T, :zero")) + end end """ @@ -678,6 +762,6 @@ export fundamental_derivative_of_gas_dynamics #volume properties export mass_density,molar_density, compressibility_factor #molar gradient properties -export chemical_potential, activity_coefficient, activity, aqueous_activity, fugacity_coefficient +export chemical_potential, activity_coefficient, activity, aqueous_activity, fugacity_coefficient,reference_chemical_potential,reference_chemical_potential_type export chemical_potential_res export mixing, excess, gibbs_solvation \ No newline at end of file diff --git a/src/methods/property_solvers/reactive/reactive.jl b/src/methods/property_solvers/reactive/reactive.jl index fe82e4b88..61a0bed02 100644 --- a/src/methods/property_solvers/reactive/reactive.jl +++ b/src/methods/property_solvers/reactive/reactive.jl @@ -1,65 +1,29 @@ -# function pH(model::ReactiveModel,p,T,n0;z0=nothing) -# ν = stoichiometric_coefficient(model) # Exists -# ΔHf = model.params.ΔHf.values -# Sf = model.params.Sf.values -# ΔGf = ΔHf-T*Sf -# ΔrG = sum(ΔGf*ν,dims=1) -# Keq = exp.(-ΔrG/Rgas(model.model)/T) - -# if isnothing(z0) -# z0 = x0_equilibrium_conditions(model,p,T,n0) -# end - -# f!(F,z) = equilibrium_conditions(model,F,p,T,n0,z,ν,Keq) -# sol = Solvers.nlsolve(f!,z0,TrustRegion(Newton(), NWI()),NEqOptions(),ForwardDiff.Chunk{1}()) # Exists -# ξ = Solvers.x_sol(sol) -# n = n0+ν*ξ - -# a = activity(model.model,p,T,n) # Exists - -# hydronium_id = find_hydronium_index(model) - -# return -log10(a[hydronium_id]) -# end - -# function find_hydronium_index(model) -# idx = findfirst(model.components.=="hydronium") -# return idx -# end - -function equilibrate(model::ReactiveModel,p,T,n0;z0=nothing) +function equilibrate(model::ReactiveEoSModel,p,T,n0;z0=nothing) ν = stoichiometric_coefficient(model) # Exists - ΔHf = model.params.ΔHf.values - Sf = model.params.Sf.values - ΔGf = ΔHf-T*Sf - ΔrG = sum(ΔGf.*ν,dims=1) - Keq = exp.(-ΔrG/Rgas(model.eosmodel)/T) - + Keq = ideal_Keq(model,T,z,ν) if isnothing(z0) z0 = x0_equilibrium_conditions(model,p,T,n0) end - - f!(F,z) = equilibrium_conditions(model,F,p,T,n0,z,ν,Keq) - sol = Solvers.nlsolve(f!,z0,TrustRegion(Newton(), NWI()),NEqOptions(),ForwardDiff.Chunk{1}()) # Exists - ξ = Solvers.x_sol(sol) + μ_ref_type = reference_chemical_potential_type(model) + μ_ref = reference_chemical_potential(model.eosmodel,p,T,μ_ref_type) + f!(F,z) = equilibrium_conditions(model,F,p,T,n0,exp10.(z),ν,Keq,μ_ref) + sol = Solvers.nlsolve(f!,log10.(z0),TrustRegion(Newton(), NWI()),NEqOptions(),ForwardDiff.Chunk{1}()) # Exists + ξ = exp10.(Solvers.x_sol(sol)) n = n0+ν*ξ return n end -function equilibrium_conditions(model::ReactiveModel,F,p,T,n0,ξ,ν,Keq) +function equilibrium_conditions(model::ReactiveModel,F,p,T,n0,ξ,ν,Keq,μ_ref) n = n0+ν*ξ - - a = activity(model.eosmodel,p,T,n) - + a = activity_impl(model,p,T,n,μ_ref,:unknown,:unknown,true,nothing) F[1:end] = sum(log.(a).*ν,dims=1) - log.(Keq) return F end -function x0_equilibrium_conditions(model::ReactiveModel,p,T,n0) - +function x0_equilibrium_conditions(model::ReactiveEoSModel,p,T,n0) z0 = zeros(length(model.components)) z0[1] = log10(p) return z0 end -export equilibrate \ No newline at end of file +export equilibrate diff --git a/src/methods/property_solvers/reactive/reactive_aq.jl b/src/methods/property_solvers/reactive/reactive_aq.jl index 2dc063ad0..693cc8889 100644 --- a/src/methods/property_solvers/reactive/reactive_aq.jl +++ b/src/methods/property_solvers/reactive/reactive_aq.jl @@ -1,51 +1,6 @@ function pH(model::ReactiveAqModel,p,T,n0;z0=nothing) n = equilibrate(model, p, T, n0; z0 = z0) - a = aqueous_activity(model.eosmodel,p,T,n) - hydronium_id = find_hydronium_index(model) - return -log10(a[hydronium_id]) end - -function equilibrate(model::ReactiveAqModel,p,T,n0;z0=nothing) - T0 = 298.15 - ν = stoichiometric_coefficient(model) - ΔGf0 = model.params.ΔHf.values - ΔHf = model.params.ΔHf.values - ΔGf = ΔGf0/T0+ΔHf*(1/T-1/T0) - ΔrG = sum(ΔGf.*ν,dims=1) - Keq = exp.(-ΔrG/Rgas(model.eosmodel)) - - # println(Keq) - - if isnothing(z0) - z0 = x0_equilibrium_conditions(model,p,T,n0) - end - - f!(F,z) = equilibrium_conditions(model,F,p,T,n0,exp10.(z),ν,Keq) - sol = Solvers.nlsolve(f!,log10.(z0),TrustRegion(Newton(), NWI()),NEqOptions(),ForwardDiff.Chunk{1}()) # Exists - ξ = exp10.(Solvers.x_sol(sol)) - n = n0+ν*ξ - return n -end - -function equilibrium_conditions(model::ReactiveAqModel,F,p,T,n0,ξ,ν,Keq) - n = n0+ν*ξ - # println(n) - - a = aqueous_activity(model.eosmodel,p,T,n) - - F[1:end] = sum(log.(a).*ν,dims=1) - log.(Keq) - # println(F) - return F -end - -function x0_equilibrium_conditions(model::ReactiveAqModel,p,T,n0) - - z0 = zeros(length(model.components)) - z0[1] = log10(p) - return z0 -end - -export equilibrate \ No newline at end of file diff --git a/src/models/Reactive/reactive.jl b/src/models/Reactive/reactive.jl index 3ca3a63c9..e6446a3b6 100644 --- a/src/models/Reactive/reactive.jl +++ b/src/models/Reactive/reactive.jl @@ -1,29 +1,42 @@ -struct ReactiveParams{T} <: EoSParam where T <: Float64 +abstract type ReactiveEoSModel <: EoSModel end + +struct ReactiveParams{T} <: ParametricEoSParam{T} ΔHf::SingleParam{T} Sf::SingleParam{T} end -struct ReactiveModel{T} +struct ReactiveModel{T} <: ReactiveEoSModel components::Vector{String} reactions::Vector params::ReactiveParams{Float64} eosmodel::T end -function ReactiveModel(components::Vector{String}, reactions::Vector; +function ReactiveModel(_components, reactions::Vector; model::EoSModel, userlocations=String[], verbose=false) + components = format_components(_components) params = getparams(components, ["properties/formation.csv"]; userlocations = userlocations, verbose = verbose) ΔHf = params["Hf"] Sf = params["Sf"] + packagedparams = ReactiveParams{Float64}(ΔHf,Sf) + return ReactiveModel(components,reactions,packagedparams,model) +end - packagedparams = ReactiveParams(ΔHf,Sf) +reference_chemical_potential_type(model::ReactiveEoSModel) = :pure - return ReactiveModel(components,reactions,packagedparams,model) +ideal_Keq(model::ReactiveEoSModel,T,z) = ideal_Keq(model,T,z,stoichiometric_coefficient(model)) + +function ideal_Keq(model::ReactiveModel,T,z,ν) + ΔHf = model.params.ΔHf.values + Sf = model.params.Sf.values + ΔGf = ΔHf-T*Sf + ΔrG = sum(ΔGf.*ν,dims=1) + Keq = exp.(-ΔrG/Rgas(model.eosmodel)/T) end -function stoichiometric_coefficient(model::ReactiveModel) +function stoichiometric_coefficient(model::ReactiveEoSModel) ν = zeros(length(model.components),length(model.reactions)) for (i,reaction) in enumerate(model.reactions) species = [reaction[k][1] for k in 1:length(reaction)] @@ -37,7 +50,7 @@ function stoichiometric_coefficient(model::ReactiveModel) return ν end -function Base.show(io::IO, mime::MIME"text/plain", model::ReactiveModel) +function Base.show(io::IO, mime::MIME"text/plain", model::ReactiveModel) print(io, typeof(model)) length(model.eosmodel) == 1 && println(io, " with 1 component:") length(model.eosmodel) > 1 && println(io, " with ", length(model.eosmodel), " components:") @@ -64,4 +77,4 @@ function Base.show(io::IO, mime::MIME"text/plain", model::ReactiveModel) end end -export ReactiveModel \ No newline at end of file +export ReactiveModel diff --git a/src/models/Reactive/reactive_aq.jl b/src/models/Reactive/reactive_aq.jl index 6716dc567..3cae5cc61 100644 --- a/src/models/Reactive/reactive_aq.jl +++ b/src/models/Reactive/reactive_aq.jl @@ -1,68 +1,36 @@ -struct ReactiveAqParams{T} <: EoSParam where T <: Float64 +struct ReactiveAqParams{T} <: ParametricEoSParam{T} ΔGf::SingleParam{T} ΔHf::SingleParam{T} end -struct ReactiveAqModel{T} +struct ReactiveAqModel{T} <: ReactiveEoSModel components::Vector{String} reactions::Vector params::ReactiveParams{Float64} eosmodel::T end -function ReactiveAqModel(components::Vector{String}, reactions::Vector; +reference_chemical_potential_type(model::ReactiveAqModel) = :aqueous + +function ReactiveAqModel(_components, reactions::Vector; model::EoSModel, userlocations=String[], verbose=false) + components = format_components(_components) params = getparams(components, ["properties/formation_aqueous.csv"]; userlocations = userlocations, verbose = verbose) ΔGf = params["Gf"] ΔHf = params["Hf"] - - packagedparams = ReactiveParams(ΔGf,ΔHf) - + packagedparams = ReactiveAqParams{Float64}(ΔGf,ΔHf) return ReactiveAqModel(components,reactions,packagedparams,model) end -function stoichiometric_coefficient(model::ReactiveAqModel) - ν = zeros(length(model.components),length(model.reactions)) - for (i,reaction) in enumerate(model.reactions) - species = [reaction[k][1] for k in 1:length(reaction)] - coeff = [reaction[k][2] for k in 1:length(reaction)] - for (j,component) in enumerate(model.components) - if component in species - ν[j,i] = coeff[findfirst(isequal(component),species)] - end - end - end - return ν -end - - -function Base.show(io::IO, mime::MIME"text/plain", model::ReactiveAqModel) - print(io, typeof(model)) - length(model.eosmodel) == 1 && println(io, " with 1 component:") - length(model.eosmodel) > 1 && println(io, " with ", length(model.eosmodel), " components:") - show_pairs(io,model.components) - println(io, "") - - length(model.reactions) == 1 && println(io, "and 1 reaction:") - length(model.reactions) > 1 && println(io, "and ", length(model.reactions), " reactions:") - for reaction in model.reactions - species = [reaction[k][1] for k in 1:length(reaction)] - coeff = [reaction[k][2] for k in 1:length(reaction)] - reactants = "" - products = "" - for i in 1:length(species) - if coeff[i] < 0 - reactants *= string(abs(coeff[i]))*" "*species[i]*" + " - else - products *= string(abs(coeff[i]))*" "*species[i]*" + " - end - end - reactants = reactants[1:end-3] - products = products[1:end-3] - println(io, reactants, " <=> ", products) - end +function ideal_Keq(model::ReactiveAqModel,T,z,ν) + T0 = 298.15 + ΔGf0 = model.params.ΔHf.values + ΔHf = model.params.ΔHf.values + ΔGf = ΔGf0/T0+ΔHf*(1/T-1/T0) + ΔrG = sum(ΔGf.*ν,dims=1) + Keq = exp.(-ΔrG/Rgas(model.eosmodel)) end -export ReactiveAqModel \ No newline at end of file +export ReactiveAqModel From 4cfb99ac532cce7ef2220aabf0f847bd96defa72 Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Wed, 14 Aug 2024 01:24:41 -0400 Subject: [PATCH 07/18] use reference_chemical_potential in more places --- .../multicomponent/solids/sle_solubility.jl | 11 ++++------ .../multicomponent/solids/slle_solubility.jl | 10 +++++----- src/models/Activity/equations.jl | 20 +++++++++++++++---- src/models/CompositeModel/CompositeModel.jl | 13 +++++++++--- src/models/CompositeModel/FluidCorrelation.jl | 12 ++++++++++- src/models/CompositeModel/GammaPhi.jl | 16 +++++++++++---- 6 files changed, 58 insertions(+), 24 deletions(-) diff --git a/src/methods/property_solvers/multicomponent/solids/sle_solubility.jl b/src/methods/property_solvers/multicomponent/solids/sle_solubility.jl index 14bdcada8..b7bd68f7c 100644 --- a/src/methods/property_solvers/multicomponent/solids/sle_solubility.jl +++ b/src/methods/property_solvers/multicomponent/solids/sle_solubility.jl @@ -65,8 +65,8 @@ function sle_solubility(model::CompositeModel,p,T,z;solute=nothing,x0=nothing) if isnothing(x0) x0 = x0_sle_solubility(model,p,T,z,idx_solv,idx_sol_l,ν_l,μsol) end - - data = (fluid_r,idx_liq_r,solid_r,idx_sol_r,idx_sol_l,idx_sol_s,idx_solv,μsol[1]) + μ_ref = reference_chemical_potential(model.fluid,p,T,reference_chemical_potential_type(model.fluid)) + data = (μ_ref,idx_sol_l,idx_solv,μsol[1]) # println(x0) f!(F,x) = obj_sle_solubility(F,model,p,T,z[idx_solv],exp10(x[1]),data,ν_l) results = Solvers.nlsolve(f!,x0,LineSearch(Newton()),NEqOptions(),ForwardDiff.Chunk{1}()) @@ -82,13 +82,13 @@ function sle_solubility(model::CompositeModel,p,T,z;solute=nothing,x0=nothing) end function obj_sle_solubility(F,model,p,T,zsolv,solu,data,ν_l) - fluid_r,idx_liq_r,solid_r,idx_sol_r,idx_sol_l,idx_sol_s,idx_solv,μsol = data + μ_ref,idx_sol_l,idx_solv,μsol = data z = zeros(typeof(solu),length(model.fluid)) z[.!(idx_solv)] .= solu z[idx_solv] .= zsolv R = Rgas(model.fluid) ∑z = sum(z) - γliq = activity_coefficient(model.fluid,p,T,z/∑z) + γliq = activity_coefficient(model.fluid,p,T,z/∑z,μ_ref = μ_ref) γl = @view(γliq[idx_sol_l]) zl = @view(z[idx_sol_l]) μliq = zero(eltype(γliq)) @@ -97,9 +97,6 @@ function obj_sle_solubility(F,model,p,T,zsolv,solu,data,ν_l) μliq_i = R*T*log(γl[i]*xli) μliq += μliq_i*ν_l[i] end - #μliq = R*T*log.(@view(γliq[idx_sol_l]).*@view(z[idx_sol_l]) ./ ∑z) - #solid_r,idx_sol_r = index_reduction(model.solid,idx_sol_s) - #μliq = dot(μliq,ν_l) F[1] = μliq - μsol return F end diff --git a/src/methods/property_solvers/multicomponent/solids/slle_solubility.jl b/src/methods/property_solvers/multicomponent/solids/slle_solubility.jl index ed2cbba2f..78a9469bc 100644 --- a/src/methods/property_solvers/multicomponent/solids/slle_solubility.jl +++ b/src/methods/property_solvers/multicomponent/solids/slle_solubility.jl @@ -17,10 +17,10 @@ function slle_solubility(model::CompositeModel,p,T) solid_r,idx_sol_r = index_reduction(model.solid,idx_sol) μsol = chemical_potential(solid_r,p,T,[1.]) - + μ_ref = reference_chemical_potential(model.fluid,p,T,:pure) x0 = x0_slle_solubility(model,p,T,μsol) - f!(F,x) = obj_slle_solubility(F,model.fluid,p,T,[exp10(x[1]),1-exp10(x[1])-exp10(x[2]),exp10(x[2])],[exp10(x[3]),1-exp10(x[3])-exp10(x[4]),exp10(x[4])],μsol) + f!(F,x) = obj_slle_solubility(F,model.fluid,p,T,[exp10(x[1]),1-exp10(x[1])-exp10(x[2]),exp10(x[2])],[exp10(x[3]),1-exp10(x[3])-exp10(x[4]),exp10(x[4])],μsol,μ_ref) results = Solvers.nlsolve(f!,x0,LineSearch(Newton())) sol = exp10.(Solvers.x_sol(results)) x1 = [sol[1],1-sol[1]-sol[2],sol[2]] @@ -28,11 +28,11 @@ function slle_solubility(model::CompositeModel,p,T) return (x1,x2) end -function obj_slle_solubility(F,liquid,p,T,x1,x2,μsol) +function obj_slle_solubility(F,liquid,p,T,x1,x2,μsol,μ_ref) nc = length(liquid) - γliq1 = activity_coefficient(liquid,p,T,x1) + γliq1 = activity_coefficient(liquid,p,T,x1,μ_ref = μ_ref) μliq1 = @. Rgas()*T*log(γliq1*x1) - γliq2 = activity_coefficient(liquid,p,T,x2) + γliq2 = activity_coefficient(liquid,p,T,x2,μ_ref = μ_ref) μliq2 = @. Rgas()*T*log(γliq2*x2) F[1] = (μliq1[end] - μsol[1])/Rgas()/T diff --git a/src/models/Activity/equations.jl b/src/models/Activity/equations.jl index d75d6055a..51ae0f52b 100755 --- a/src/models/Activity/equations.jl +++ b/src/models/Activity/equations.jl @@ -27,10 +27,22 @@ function activity_coefficient(model::ActivityModel,p,T,z) return exp.(Solvers.gradient(x->excess_gibbs_free_energy(model,p,T,x),z)/(R̄*T))::X end -function activity(model::ActivityModel,p,T,z=SA[1.]) - γ = activity_coefficient(model, p, T, z) - x = z ./ sum(z) - return γ .* x +function activity_coefficient_impl(model::ActivityModel,p,T,z,μ_ref,reference,phase,threaded,vol0) + #TODO: what to do if the reference is not pure? + return activity_coefficient(model,p,T,z) +end + +reference_chemical_potential_type(model::ActivityModel) = :zero + +function activity(model::ActivityModel,p,T,z) + γ = activity_coefficient(model,p,T,z) + ∑z = sum(z) + return γ .* z ./ ∑z +end + +function activity_impl(model::ActivityModel,p,T,z,μ_ref,reference,phase,threaded,vol0) + #TODO: what to do if the reference is not pure? + return activity(model,p,T,z) end function test_activity_coefficient(model::ActivityModel,p,T,z) diff --git a/src/models/CompositeModel/CompositeModel.jl b/src/models/CompositeModel/CompositeModel.jl index 8a052837f..f92c25440 100644 --- a/src/models/CompositeModel/CompositeModel.jl +++ b/src/models/CompositeModel/CompositeModel.jl @@ -62,7 +62,7 @@ julia> bubble_pressure(model,300.15,[0.9,0.1]) (2694.150594740186, 0.00016898441224336215, 0.9239727973658585, [0.7407077952279438, 0.2592922047720562]) #using a correlation-based fluid -julia> fluidmodel = CompositeModel(["octane","heptane"],liquid = RackettLiquid,saturation = DIPPR101Sat,gas = BasicIdeal); +julia> fluidmodel = CompositeModel(["octane","heptane"],liquid = RackettLiquid,saturation = DIPPR101Sat,gas = BasicIdeal); model2 = CompositeModel(["octane","heptane"],liquid = UNIFAC, fluid = fluidmodel) Composite Model (γ-ϕ) with 2 components: Activity Model: UNIFAC{PR{BasicIdeal, PRAlpha, NoTranslation, vdW1fRule}}("octane", "heptane") @@ -286,10 +286,17 @@ function volume_impl(model::CompositeModel,p,T,z,phase,threaded,vol0) end end -function activity_coefficient(model::CompositeModel,p,T,z=SA[1.]; phase = :unknown, threaded=true) - return activity_coefficient(model.fluid,p,T,z;phase,threaded) +function activity_coefficient(model::CompositeModel,p,T,z=SA[1.]; + μ_ref = nothing, + reference = :pure, + phase=:unknown, + threaded=true, + vol0=nothing) + return activity_coefficient(model.fluid,p,T,z;μ_ref,reference,phase,threaded,vol0) end +reference_chemical_potential_type(model::CompositeModel) = reference_chemical_potential_type(model.fluid) + function init_preferred_method(method::typeof(saturation_pressure),model::CompositeModel,kwargs) return init_preferred_method(saturation_pressure,model.fluid,kwargs) end diff --git a/src/models/CompositeModel/FluidCorrelation.jl b/src/models/CompositeModel/FluidCorrelation.jl index 512d0fcda..bee81018a 100644 --- a/src/models/CompositeModel/FluidCorrelation.jl +++ b/src/models/CompositeModel/FluidCorrelation.jl @@ -11,7 +11,17 @@ struct FluidCorrelation{V,L,Sat} <: RestrictedEquilibriaModel end __gas_model(model::FluidCorrelation) = model.gas -activity_coefficient(model::FluidCorrelation, p, T,z=SA[1.]; phase = :unknown, threaded=true) = FillArrays.Ones(length(model)) + +function activity_coefficient(model::FluidCorrelation,p,T,z=SA[1.]; + μ_ref = nothing, + reference = :pure, + phase=:unknown, + threaded=true, + vol0=nothing) + return FillArrays.Ones(length(model)) +end + +reference_chemical_potential_type(model::FluidCorrelation) = :zero function volume_impl(model::FluidCorrelation, p, T, z, phase, threaded, vol0) _0 = zero(p+T+first(z)) diff --git a/src/models/CompositeModel/GammaPhi.jl b/src/models/CompositeModel/GammaPhi.jl index cbd4e938b..197433ddb 100644 --- a/src/models/CompositeModel/GammaPhi.jl +++ b/src/models/CompositeModel/GammaPhi.jl @@ -19,10 +19,18 @@ function init_preferred_method(method::typeof(saturation_temperature),model::Gam return init_preferred_method(saturation_temperature,model.fluid,kwargs) end -function activity_coefficient(model::GammaPhi,p,T,z=SA[1.]; phase = :unknown, threaded=true) - return activity_coefficient(model.activity,p,T,z) +function activity_coefficient(model::GammaPhi,p,T,z=SA[1.]; + μ_ref = nothing, + reference = :pure, + phase=:unknown, + threaded=true, + vol0=nothing) + + return activity_coefficient(model.activity,p,T,z;μ_ref,reference,phase,threaded,vol0) end +reference_chemical_potential_type(model::GammaPhi) = reference_chemical_potential_type(model.activity) + function saturation_pressure(model::GammaPhi,T,method::SaturationMethod) return saturation_pressure(model.fluid,T,method) end @@ -198,7 +206,7 @@ function __eval_G_DETPFlash(wrapper::PTFlashWrapper{<:GammaPhi},p,T,x,equilibriu return gl,vl else throw(error("γ-ϕ Composite Model does not support VLE calculation with `DETPFlash`. if you want to calculate LLE equilibria, try using `DETPFlash(equilibrium = :lle)`")) - #= + #= vv = volume(model.fluid.model,p,T,x,phase = :v) gv = VT_gibbs_free_energy(model.fluid.model, vv, T, x) if gv > gl @@ -259,7 +267,7 @@ function tpd_obj(model::GammaPhi, p, T, di, isliquid, cache = tpd_neq_cache(mode lnγw = γ fx = @sum(w[i]*(lnγw[i] + log(w[i]) - di[i])) - sum(w) + 1 end - + obj = Solvers.ADScalarObjective(f,di,ForwardDiff.Chunk{2}()) optprob = OptimizationProblem(obj = obj,inplace = true) end From 359fa24ae69049e1c0472cd81719a82359acf38e Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Wed, 14 Aug 2024 01:26:35 -0400 Subject: [PATCH 08/18] add exception for activity model in the general activity definition --- src/methods/pT.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/methods/pT.jl b/src/methods/pT.jl index 7f01d7b1d..890b402b5 100755 --- a/src/methods/pT.jl +++ b/src/methods/pT.jl @@ -505,7 +505,9 @@ function activity(model::EoSModel,p,T,z=SA[1.]; phase=:unknown, threaded=true, vol0=nothing) - + if model isa ActivityModel + return activity(model,p,T,z) + end if μ_ref == nothing return activity_impl(model,p,T,z,reference_chemical_potential(model,p,T,reference,phase;threaded),reference,phase,threaded,vol0) else From 93a61b79829b32cfe38c476dce0ddf40b515db50 Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Wed, 14 Aug 2024 01:45:39 -0400 Subject: [PATCH 09/18] typo in reference_chemical_potential --- src/methods/pT.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/methods/pT.jl b/src/methods/pT.jl index 890b402b5..f091764ec 100755 --- a/src/methods/pT.jl +++ b/src/methods/pT.jl @@ -509,7 +509,7 @@ function activity(model::EoSModel,p,T,z=SA[1.]; return activity(model,p,T,z) end if μ_ref == nothing - return activity_impl(model,p,T,z,reference_chemical_potential(model,p,T,reference,phase;threaded),reference,phase,threaded,vol0) + return activity_impl(model,p,T,z,reference_chemical_potential(model,p,T,reference;phase,threaded),reference,phase,threaded,vol0) else return activity_impl(model,p,T,z,μ_ref,reference,phase,threaded,vol0) end @@ -575,7 +575,7 @@ Returns a reference chemical potential. used in calculation of `activity` and ac - `:zero`: the reference potential is equal to zero for all components (used for `ActivityModel`) The keywords `phase`, `threaded` and `vol0` are passed to the [`Clapeyron.volume`](@ref) solver. """ -function reference_chemical_potential(model::EoSModel,p,T,reference = reference_chemical_potential_type(model), phase=:unknown, threaded=true, vol0=nothing) +function reference_chemical_potential(model::EoSModel,p,T,reference = reference_chemical_potential_type(model); phase=:unknown, threaded=true, vol0=nothing) if reference == :pure pure = split_model.(model) return gibbs_free_energy.(pure, p, T; phase, threaded) From 7b5e4897f30f12b133d2656e5264c1bb9ee4c6a8 Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Thu, 15 Aug 2024 19:33:58 -0400 Subject: [PATCH 10/18] add missing parameters for joback --- database/ideal/JobackIdeal.csv | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/database/ideal/JobackIdeal.csv b/database/ideal/JobackIdeal.csv index 181745dce..43ea57e31 100755 --- a/database/ideal/JobackIdeal.csv +++ b/database/ideal/JobackIdeal.csv @@ -10,8 +10,8 @@ H2O,0.9725540967209461,0.05527786574464416,38.44807453383915,174.95,150.65,-242. =CH-~|~=CH~|~CH=,0.0129,-0.0006,46.0,24.96,8.73,37.97,48.53,-8.00,1.05E-1,-9.63E-5,3.56E-8,2.691,2.205,82.28,-0.242, =C<,0.0117,0.0011,38.0,24.14,11.14,83.99,92.36,-2.81E+1,2.08E-1,-3.06E-4,1.46E-7,3.063,2.138,0.0,0.0, =C=,0.0026,0.0028,36.0,26.15,17.78,142.14,136.70,2.74E+1,-5.57E-2,1.01E-4,-5.02E-8,4.720,2.661,0.0,0.0, -CH,0.0027,-0.0008,46.0,9.20,-11.18,79.30,77.71,2.45E+1,-2.71E-2,1.11E-4,-6.78E-8,2.322,1.155,0.0,0.0, -C,0.0020,0.0016,37.0,27.38,64.32,115.51,109.82,7.87,2.01E-2,-8.33E-6,1.39E-9,4.151,3.302,0.0,0.0, +CH~|~≡CH,0.0027,-0.0008,46.0,9.20,-11.18,79.30,77.71,2.45E+1,-2.71E-2,1.11E-4,-6.78E-8,2.322,1.155,0.0,0.0, +C~|~≡C−,0.0020,0.0016,37.0,27.38,64.32,115.51,109.82,7.87,2.01E-2,-8.33E-6,1.39E-9,4.151,3.302,0.0,0.0, ring-CH2-~|~CH2_hex~|~CH2_pent~|~cCH2_hex~|~cCH2_pent,0.0100,0.0025,48.0,27.15,7.75,-26.80,-3.68,-6.03,8.54E-2,-8.00E-6,-1.80E-8,0.490,2.398,307.53,-0.798, ring>CH-~|~CH_hex~|~CH_pent~|~cCH_hex~|~cCH_pent,0.0122,0.0004,38.0,21.78,19.88,8.67,40.99,-2.05E+1,1.62E-1,-1.60E-4,6.24E-8,3.243,1.942,-394.29,1.251, ring>C<,0.0042,0.0061,27.0,21.32,60.15,79.72,87.88,-9.09E+1,5.57E-1,-9.00E-4,4.69E-7,-1.373,0.644,0.0,0.0, @@ -43,6 +43,6 @@ O=CH- (aldehyde)~|~CH=O,0.0379,0.0030,82.0,72.24,36.90,-162.03,-143.48,3.09E+1,- -SH,0.0031,0.0084,63.0,63.56,20.09,-17.33,-22.99,3.53E+1,-7.58E-2,1.85E-4,-1.03E-7,2.360,6.884,0.0,0.0, -S- (non-ring),0.0119,0.0049,54.0,68.78,34.40,41.87,33.12,1.96E+1,-5.61E-3,4.02E-5,-2.76E-8,4.130,6.817,0.0,0.0, -S- (ring),0.0019,0.0051,38.0,52.10,79.93,39.10,27.76,1.67E+1,4.81E-3,2.77E-5,-2.11E-8,1.557,5.984,0.0,0.0, -C≡CH~|~CH#C,0.0,0.0,0.0,0.0,0.0,0.0,0.0,32.37,-0.007,0.00010267,-6.641e-08,0.0,0.0,0.0,0.0,0.0 -OCH3,0.0,0.0,0.0,0.0,0.0,0.0,0.0,45.0,-0.07128,0.000264,-1.515e-07,0.0,0.0,0.0,0.0,0.0 -OCH2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,24.591,0.0318,5.66e-05,-4.29e-08,0.0,0.0,0.0,0.0,0.0 +C≡CH~|~CH#C,0.0047,0.0008,83.0,36.58,53.14,194.81,187.53,32.37,-0.007,0.00010267,-6.641e-08,6.473,4.457,0.0,0.0, +OCH3,0.0309,0.0003,83.0,46.0,17.13,-208.67,-148.96,45.0,-0.07128,0.000264,-1.515e-07,2.096,4.783,670.38,-2.105, +OCH2,0.0357,0.0015,74.0,45.3, 33.5,-152.86,-96.58,24.591,0.0318,5.66e-05,-4.29e-08,3.778,4.636,216.25,-0.585, From 562fd3966b9b8ba6d8b09717823557794df71f95 Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Thu, 15 Aug 2024 21:05:35 -0400 Subject: [PATCH 11/18] add more joback functions --- src/models/ideal/JobackIdeal.jl | 235 ++++++++++++++++++++++++++------ test/test_models_others.jl | 12 ++ 2 files changed, 208 insertions(+), 39 deletions(-) diff --git a/src/models/ideal/JobackIdeal.jl b/src/models/ideal/JobackIdeal.jl index 6a9e22568..66dbe6db3 100755 --- a/src/models/ideal/JobackIdeal.jl +++ b/src/models/ideal/JobackIdeal.jl @@ -166,71 +166,228 @@ end ReidIdeal(model::JobackIdeal) = ReidIdeal(model.components,ReidIdealParam(model.params.coeffs,model.params.reference_state),model.references) +""" + JobackGC + +Module containing group contribution calculations using the joback method. the available functions are: + +- `JobackGC.T_c(model::JobackModel)`: critical temperature (in K) +- `JobackGC.P_c(model::JobackModel)`: critical pressure (in Pa) +- `JobackGC.V_c(model::JobackModel)`: critical volume (in m3/mol) +- `JobackGC.T_b(model::JobackModel)`: normal boiling point (in K) +- `JobackGC.H_form(model::JobackModel)`: enthalpy of formation at 298K, ideal gas (in J/mol) +- `JobackGC.G_form(model::JobackModel)`: gibbs energy of formation at 298K, ideal gas (in J/mol) +- `JobackGC.S_form(model::JobackModel)`: entropy of formation at 298K, ideal gas (in J/mol/K) +- `JobackGC.H_fusion(model::JobackModel)`: enthalpy of fusion (in J/mol, at 1 atm) +- `JobackGC.H_vap(model::JobackModel)`: molar enthalpy of vaporization (in J/mol, at normal boiling point) +- `JobackGC.C_p(model::JobackModel, T)`: ideal gas isobaric heat capacity (in J/mol/K) +""" +module JobackGC + using Clapeyron: JobackIdeal + using LinearAlgebra: dot +""" + JobackGC.T_b(model::JobackIdeal)::Vector{Float64} + +Given a `JobackIdeal` model, returns a vector containing normal boiling points (in K), for each component. +""" function T_b(model::JobackIdeal) + result = zeros(Float64,length(model)) n = model.groups.n_flattenedgroups T_b = model.params.T_b.values - res = 0.0 - @inbounds begin - ni = n[1] - res += ∑(T_b[j]*ni[j] for j in @groups(1)) + for i in 1:length(model) + result[i] = dot(n[i],T_b) + 198.2 end - res = res + 198.2 + return result end +""" + JobackGC.T_c(model::JobackIdeal)::Vector{Float64} + +Given a `JobackIdeal` model, returns a vector containing critical temperatures (in K), for each component. +""" function T_c(model::JobackIdeal,Tb=T_b(model)) n = model.groups.n_flattenedgroups T_c = model.params.T_c.values - res = 0.0 - @inbounds begin - ni = n[1] - ΣT_ci = ∑(T_c[j]*ni[j] for j in @groups(1)) - res += Tb * (0.584+0.965*ΣT_ci - ΣT_ci^2)^-1 + result = zeros(Float64,length(model)) + for i in 1:length(model) + ΣT_ci = dot(n[i],T_c) + result[i] = Tb[i] / (0.584+0.965*ΣT_ci - ΣT_ci^2) end - return res + return result end -#this style of writing is uglier but faster, -#ideally they should be of the same speed +""" + JobackGC.V_c(model::JobackIdeal)::Vector{Float64} + +Given a `JobackIdeal` model, returns a vector containing critical volumes (in m3/mol), for each component. +""" function V_c(model::JobackIdeal) n = model.groups.n_flattenedgroups V_c = model.params.V_c.values - res = 0.0 - groups = @groups(1) - ni = n[1] - @inbounds begin - ΣV_ci = zero(res) - for idx in 1:length(groups) - j = groups[idx] - ΣV_ci +=V_c[j]*ni[j] - end - res += 17.5 + ΣV_ci + result = zeros(Float64,length(model)) + for i in 1:length(model) + result[i] = dot(n[i],V_c) + 17.5 end - #result in mL/mol, converting to m3/mol - return res*1e-6 + result .*= 1e-6 + return result end +""" + JobackGC.P_c(model::JobackIdeal)::Vector{Float64} + +Given a `JobackIdeal` model, returns a vector containing critical pressures (in Pa), for each component. +""" function P_c(model::JobackIdeal) n = model.groups.n_flattenedgroups P_c = model.params.P_c.values N_a = model.params.N_a.values - groups = @groups(1) - ni = n[1] - @inbounds begin - ΣP_ci = 0.0 - ΣN_a = 0 - for idx in 1:length(groups) - j = groups[idx] - ΣP_ci += P_c[j]*ni[j] - ΣN_a += N_a[j]*ni[j] - end - res = (0.113 + 0.0032*ΣN_a - ΣP_ci)^-2 + result = zeros(Float64,length(model)) + for i in 1:length(model) + ΣP_ci = dot(n[i],P_c) + ΣN_a = dot(n[i],N_a) + result[i] = (0.113 + 0.0032*ΣN_a - ΣP_ci)^-2 end #result in bar, converting to Pa - return res*100000.0 + result .*= 1e5 + return result +end + +""" + G_form(model::JobackIdeal)::Vector{Float64} + +Given a `JobackIdeal` model, returns a vector containing gibbs energies of formation (J/mol, at ideal gas, 298K), for each component. +""" +function G_form(model::JobackIdeal) + n = model.groups.n_flattenedgroups + G_form = model.params.G_form.values + result = zeros(Float64,length(model)) + for i in 1:length(model) + result[i] = dot(n[i],G_form) + 53.88 + end + #result in kJ, converting to J + result .*= 1000 + return result end +""" + H_form(model::JobackIdeal)::Vector{Float64} + +Given a `JobackIdeal` model, returns a vector containing entalpies of formation (J/mol, at ideal gas, 298K), for each component. +""" +function H_form(model::JobackIdeal) + n = model.groups.n_flattenedgroups + H_form = model.params.H_form.values + result = zeros(Float64,length(model)) + for i in 1:length(model) + result[i] = dot(n[i],H_form) + 68.29 + end + #result in kJ, converting to J + result .*= 1000 + return result +end + +""" + S_form(model::JobackIdeal)::Vector{Float64} + +Given a `JobackIdeal` model, returns a vector containing entropies of formation (J/mol/K, at ideal gas, 298K), for each component. +""" +function S_form(model::JobackIdeal) + n = model.groups.n_flattenedgroups + H_form = model.params.H_form.values + G_form = model.params.G_form.values + result = zeros(Float64,length(model)) + for i in 1:length(model) + H_formᵢ = dot(n[i],H_form) + 68.29 + G_formᵢ = dot(n[i],G_form) + 53.88 + result[i] = (H_formᵢ - G_formᵢ)/298 + end + #result in kJ, converting to J + result .*= 1000 + return result +end + +""" + H_fusion(model::JobackIdeal)::Vector{Float64} + +Given a `JobackIdeal` model, returns a vector containing enthalpies of fusion (J/mol, at 1 atm), for each component. +""" +function H_fusion(model::JobackIdeal) + n = model.groups.n_flattenedgroups + H_fusion = model.params.H_fusion.values + result = zeros(Float64,length(model)) + for i in 1:length(model) + result[i] = dot(n[i],H_fusion) - 0.88 + end + #result in kJ, converting to J + result .*= 1000 + return result +end + +""" + H_vap(model::JobackIdeal)::Vector{Float64} + +Given a `JobackIdeal` model, returns a vector containing enthalpies of vaporization (J/mol, at normal boiling point), for each component. +""" +function H_vap(model::JobackIdeal) + n = model.groups.n_flattenedgroups + H_vap = model.params.H_vap.values + result = zeros(Float64,length(model)) + for i in 1:length(model) + result[i] = dot(n[i],H_vap) + 15.3 + end + #result in kJ, converting to J + result .*= 1000 + return result +end + +""" + C_p(model::JobackIdeal,T)::Vector + +Given a `JobackIdeal` model, returns a vector containing ideal gas isobaric heat capacities (J/mol/K), for each component. +""" +function C_p(model::JobackIdeal,T) + n = model.groups.n_flattenedgroups + a = model.params.a.values + b = model.params.b.values + c = model.params.c.values + d = model.params.d.values + result = zeros(Base.promote_eltype(T,Float64),length(model)) + for i in 1:length(model) + ai = dot(n[i],a) - 37.93 + bi = dot(n[i],b) + 0.210 + ci = dot(n[i],c) - 3.91e-4 + di = dot(n[i],d) + 2.06e-7 + result[i] = evalpoly(T, (ai,bi,ci,di)) + end + return result +end + +""" + Visc(model::JobackIdeal,T)::Vector + +Given a `JobackIdeal` model, returns a vector containing enthalpies of vaporization (J/mol, at normal boiling point), for each component. +""" +function Visc(model::JobackIdeal,T) + n = model.groups.n_flattenedgroups + ηa = model.params.eta_a.values + ηb = model.params.eta_b.values + Mw = model.params.Mw.values + result = zeros(Base.promote_eltype(T,Float64),length(model)) + for i in 1:length(model) + ηai = dot(n[i],ηa) + ηbi = dot(n[i],ηb) + Mwi = dot(n[i],Mw) + result[i] = Mwi*exp((ηai - 597.82)/T + ηbi - 11.202) + end + return result +end + +end #JobackGC module + function crit_pure(model::JobackIdeal) - return (T_c(model),P_c(model),V_c(model)) + return (JobackGC.T_c(model),JobackGC.P_c(model),JobackGC.V_c(model)) end -export JobackIdeal \ No newline at end of file + + +export JobackIdeal, JobackGC \ No newline at end of file diff --git a/test/test_models_others.jl b/test/test_models_others.jl index c92136054..9fff3cab0 100644 --- a/test/test_models_others.jl +++ b/test/test_models_others.jl @@ -81,6 +81,18 @@ end @test Clapeyron.crit_pure(system)[1] ≈ 500.2728274871347 rtol = 1e-6 @test Clapeyron.a_ideal(system,V,T,z) ≈ 9.210841420941021 rtol = 1e-6 @test Clapeyron.ideal_consistency(system,V,T,z) ≈ 0.0 atol = 1e-14 + + s0 = JobackIdeal("acetone") + @test Clapeyron.JobackGC.T_c(s0)[1] ≈ 500.5590 rtol = 1e-6 + @test Clapeyron.JobackGC.P_c(s0)[1] ≈ 48.025e5 rtol = 1e-6 + @test Clapeyron.JobackGC.V_c(s0)[1] ≈ 209.5e-6 rtol = 1e-6 + @test Clapeyron.JobackGC.T_b(s0)[1] ≈ 322.1100 rtol = 1e-6 + @test Clapeyron.JobackGC.H_form(s0)[1] ≈ −217.83e3 rtol = 1e-6 + @test Clapeyron.JobackGC.G_form(s0)[1] ≈ −154.54e3 rtol = 1e-6 + @test Clapeyron.JobackGC.C_p(s0,300)[1] ≈ 75.3264 rtol = 1e-6 + @test Clapeyron.JobackGC.H_fusion(s0)[1] ≈ 5.1250e3 rtol = 1e-6 + @test Clapeyron.JobackGC.H_vap(s0)[1] ≈ 29.0180e3 rtol = 1e-6 + @test Clapeyron.JobackGC.Visc(s0,300)[1] ≈ 0.0002942 rtol = 9e-4 end @testset "Reid" begin From 76427d0278b9574e3af3b573f823379f06635d6f Mon Sep 17 00:00:00 2001 From: longemen3000 Date: Thu, 15 Aug 2024 21:07:12 -0400 Subject: [PATCH 12/18] typo in joback doc --- src/models/ideal/JobackIdeal.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/models/ideal/JobackIdeal.jl b/src/models/ideal/JobackIdeal.jl index 66dbe6db3..dc966b680 100755 --- a/src/models/ideal/JobackIdeal.jl +++ b/src/models/ideal/JobackIdeal.jl @@ -181,6 +181,7 @@ Module containing group contribution calculations using the joback method. the a - `JobackGC.H_fusion(model::JobackModel)`: enthalpy of fusion (in J/mol, at 1 atm) - `JobackGC.H_vap(model::JobackModel)`: molar enthalpy of vaporization (in J/mol, at normal boiling point) - `JobackGC.C_p(model::JobackModel, T)`: ideal gas isobaric heat capacity (in J/mol/K) +- `JobackGC.Visc(model::JobackModel, T)`: liquid dynamic viscocity (in Pa*s) """ module JobackGC using Clapeyron: JobackIdeal @@ -365,7 +366,7 @@ end """ Visc(model::JobackIdeal,T)::Vector -Given a `JobackIdeal` model, returns a vector containing enthalpies of vaporization (J/mol, at normal boiling point), for each component. +Given a `JobackIdeal` model, returns a vector containing liquid dynamic viscocities (Pa*s), for each component. """ function Visc(model::JobackIdeal,T) n = model.groups.n_flattenedgroups From 5999feaec7ac2edd93e82d1821f65416b44cc229 Mon Sep 17 00:00:00 2001 From: pw0908 Date: Tue, 20 Aug 2024 19:01:46 -0700 Subject: [PATCH 13/18] Small fixes --- src/methods/pT.jl | 10 ++++++++-- src/methods/property_solvers/reactive/reactive.jl | 4 ++-- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/methods/pT.jl b/src/methods/pT.jl index f091764ec..8868ec61f 100755 --- a/src/methods/pT.jl +++ b/src/methods/pT.jl @@ -522,13 +522,19 @@ function activity_impl(model,p,T,z,μ_ref,reference,phase,threaded,vol0) end function find_hydronium_index(model) - idx = findfirst(isequal("hydronium"),components) + idx = findfirst(isequal("hydronium"),model.components) + idx == nothing && return 0 + return idx +end + +function find_hydroxide_index(model) + idx = findfirst(isequal("hydroxide"),model.components) idx == nothing && return 0 return idx end function find_water_indx(model) - idx = findfirst(isequal("water"),components) + idx = findfirst(isequal("water"),model.components) idx == nothing && return 0 return idx end diff --git a/src/methods/property_solvers/reactive/reactive.jl b/src/methods/property_solvers/reactive/reactive.jl index 61a0bed02..d7d15c556 100644 --- a/src/methods/property_solvers/reactive/reactive.jl +++ b/src/methods/property_solvers/reactive/reactive.jl @@ -1,6 +1,6 @@ function equilibrate(model::ReactiveEoSModel,p,T,n0;z0=nothing) ν = stoichiometric_coefficient(model) # Exists - Keq = ideal_Keq(model,T,z,ν) + Keq = ideal_Keq(model,T,n0,ν) if isnothing(z0) z0 = x0_equilibrium_conditions(model,p,T,n0) end @@ -15,7 +15,7 @@ end function equilibrium_conditions(model::ReactiveModel,F,p,T,n0,ξ,ν,Keq,μ_ref) n = n0+ν*ξ - a = activity_impl(model,p,T,n,μ_ref,:unknown,:unknown,true,nothing) + a = activity_impl(model.eosmodel,p,T,n,μ_ref,:unknown,:unknown,true,nothing) F[1:end] = sum(log.(a).*ν,dims=1) - log.(Keq) return F end From 77522739df77f09193fce05c1e9d548411f03907 Mon Sep 17 00:00:00 2001 From: pw0908 Date: Tue, 20 Aug 2024 19:02:12 -0700 Subject: [PATCH 14/18] Fix aqueous reference (not ideal but we'll see) --- .../property_solvers/reactive/reactive_aq.jl | 54 +++++++++++++++++++ src/models/Reactive/reactive_aq.jl | 8 +-- 2 files changed, 58 insertions(+), 4 deletions(-) diff --git a/src/methods/property_solvers/reactive/reactive_aq.jl b/src/methods/property_solvers/reactive/reactive_aq.jl index 693cc8889..feacd740f 100644 --- a/src/methods/property_solvers/reactive/reactive_aq.jl +++ b/src/methods/property_solvers/reactive/reactive_aq.jl @@ -1,6 +1,60 @@ +function equilibrate(model::ReactiveAqModel,p,T,n0;z0=nothing) + ν = stoichiometric_coefficient(model) # Exists + Keq = ideal_Keq(model,T,n0,ν) + println(Keq) + if isnothing(z0) + z0 = x0_equilibrium_conditions(model,p,T,n0) + end + f!(F,z) = equilibrium_conditions(model,F,p,T,n0,exp10.(z),ν,Keq) + sol = Solvers.nlsolve(f!,log10.(z0),TrustRegion(Newton(), NWI()),NEqOptions(),ForwardDiff.Chunk{1}()) # Exists + ξ = exp10.(Solvers.x_sol(sol)) + n = n0+ν*ξ + return n +end + function pH(model::ReactiveAqModel,p,T,n0;z0=nothing) n = equilibrate(model, p, T, n0; z0 = z0) a = aqueous_activity(model.eosmodel,p,T,n) hydronium_id = find_hydronium_index(model) + hydroxide_id = find_hydroxide_index(model) + + idx_w = find_water_indx(model) + + zref = ones(length(model.components)).*1e-30 + zref[idx_w] = 1. + zref[hydronium_id] = 0.01801528*1. + zref[hydroxide_id] = 0.01801528*1. + + zref ./= sum(zref) + + μref = chemical_potential(model.eosmodel,p,T,zref) + μmix = chemical_potential(model.eosmodel,p,T,n) + + a = exp.((μmix .- μref)./R̄./T) + return -log10(a[hydronium_id]) end + +function equilibrium_conditions(model::ReactiveAqModel,F,p,T,n0,ξ,ν,Keq) + n = n0+ν*ξ + + nreact = length(model.reactions) + + μmix = chemical_potential(model.eosmodel,p,T,n) + + for i in 1:nreact + idx_w = find_water_indx(model) + + zref = ones(length(model.components)).*1e-30 + zref[idx_w] = 1. + zref[ν[:,i] .!= 0 .&& model.eosmodel.charge .!= 0] .= 0.01801528*1. + + zref ./= sum(zref) + + μref = chemical_potential(model.eosmodel,p,T,zref) + + a = exp.((μmix .- μref)./R̄./T) + F[i] = sum(log.(a).*ν[:,i]) - log(Keq[i]) + end + return F +end \ No newline at end of file diff --git a/src/models/Reactive/reactive_aq.jl b/src/models/Reactive/reactive_aq.jl index 3cae5cc61..deee4431e 100644 --- a/src/models/Reactive/reactive_aq.jl +++ b/src/models/Reactive/reactive_aq.jl @@ -3,10 +3,10 @@ struct ReactiveAqParams{T} <: ParametricEoSParam{T} ΔHf::SingleParam{T} end -struct ReactiveAqModel{T} <: ReactiveEoSModel +struct ReactiveAqModel{T} <: ReactiveEoSModel where T <: EoSModel components::Vector{String} reactions::Vector - params::ReactiveParams{Float64} + params::ReactiveAqParams{Float64} eosmodel::T end @@ -26,9 +26,9 @@ end function ideal_Keq(model::ReactiveAqModel,T,z,ν) T0 = 298.15 - ΔGf0 = model.params.ΔHf.values + ΔGf0 = model.params.ΔGf.values ΔHf = model.params.ΔHf.values - ΔGf = ΔGf0/T0+ΔHf*(1/T-1/T0) + ΔGf = ΔGf0/T0-ΔHf*(1/T-1/T0) ΔrG = sum(ΔGf.*ν,dims=1) Keq = exp.(-ΔrG/Rgas(model.eosmodel)) end From 1caa61b85cbe4e157bee248c0c48a1ac6562170c Mon Sep 17 00:00:00 2001 From: pw0908 Date: Wed, 21 Aug 2024 09:01:39 -0700 Subject: [PATCH 15/18] Minor issue with GCMSABorn --- src/models/Electrolytes/Ion/GCMSABorn.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/models/Electrolytes/Ion/GCMSABorn.jl b/src/models/Electrolytes/Ion/GCMSABorn.jl index efd5608ba..4e9df2d42 100644 --- a/src/models/Electrolytes/Ion/GCMSABorn.jl +++ b/src/models/Electrolytes/Ion/GCMSABorn.jl @@ -52,6 +52,11 @@ function GCMSABorn(solvents,ions; RSPmodel=ConstW, userlocations=String[], verbo mix_segment!(groups,shapefactor.values,segment.values) sigma = params["sigma"] + + if typeof(sigma) <: PairParam + sigma = SingleParam("sigma",components,diagvalues(sigma.values)[:]) + end + sigma.values .*= 1E-10 gc_sigma = deepcopy(sigma) gc_sigma.values .^= 3 From 72131e95d88a6ee4fa6459de4d6b71990b3a44c4 Mon Sep 17 00:00:00 2001 From: pw0908 Date: Wed, 21 Aug 2024 16:47:06 -0700 Subject: [PATCH 16/18] Update database --- database/Electrolytes/properties/charges.csv | 1 + .../PCSAFT/gcPCPSAFT/homo/HomogcPCPSAFT_like.csv | 2 +- .../SAFT/SAFTgammaMie/SAFTgammaMie_groups.csv | 15 +++++++++++++++ database/properties/formation_aqueous.csv | 6 +++++- database/properties/molarmass_groups.csv | 3 +++ 5 files changed, 25 insertions(+), 2 deletions(-) diff --git a/database/Electrolytes/properties/charges.csv b/database/Electrolytes/properties/charges.csv index 95baa71e2..2d59876a6 100644 --- a/database/Electrolytes/properties/charges.csv +++ b/database/Electrolytes/properties/charges.csv @@ -5,6 +5,7 @@ water~|~water08,0, acetonitrile,0, methanol,0, carbon dioxide,0, +acetic acid,0, hydronium~|~H3O+,1, lithium~|~Li+,1, sodium~|~Na+,1, diff --git a/database/SAFT/PCSAFT/gcPCPSAFT/homo/HomogcPCPSAFT_like.csv b/database/SAFT/PCSAFT/gcPCPSAFT/homo/HomogcPCPSAFT_like.csv index 1cb2bc62c..387ade9cf 100644 --- a/database/SAFT/PCSAFT/gcPCPSAFT/homo/HomogcPCPSAFT_like.csv +++ b/database/SAFT/PCSAFT/gcPCPSAFT/homo/HomogcPCPSAFT_like.csv @@ -16,7 +16,7 @@ CH_pent~|~cCH_pent,13.01854,0.03314,7.719,1297.7,0.0,0,0,10.1007/s10765-023-0329 CH_arom~|~aCH,13.01854,0.42335,3.727,274.41,0.0,0,0,10.1007/s10765-023-03290-3 C_arom~|~aC,12.0107,0.15371,3.9622,527.2,0.0,0,0,10.1007/s10765-023-03290-3 CH=O,29.01754,1.5774,2.8035,242.99,2.4556,0,1,10.1007/s10765-023-03290-3 ->C=O,28.0097,1.223,2.8124,249.04,3.2432,0,1,10.1007/s10765-023-03290-3 +C=O,28.0097,1.223,2.8124,249.04,3.2432,0,1,10.1007/s10765-023-03290-3 OCH3,31.03322,1.6539,3.0697,196.05,1.3866,0,1,10.1007/s10765-023-03290-3 OCH2,30.02538,1.1349,3.2037,187.13,2.744,0,1,10.1007/s10765-023-03290-3 HCOO,45.01654,1.7525,2.9043,229.63,2.7916,0,2,10.1007/s10765-023-03290-3 diff --git a/database/SAFT/SAFTgammaMie/SAFTgammaMie_groups.csv b/database/SAFT/SAFTgammaMie/SAFTgammaMie_groups.csv index 378e13b6b..588857c08 100755 --- a/database/SAFT/SAFTgammaMie/SAFTgammaMie_groups.csv +++ b/database/SAFT/SAFTgammaMie/SAFTgammaMie_groups.csv @@ -34,6 +34,21 @@ acetone~|~propanone,"[""CH3COCH3""=>1]" propanol~|~1-propanol~|~propan-1-ol,"[""CH2OH""=>1,""CH2""=>1,""CH3""=>1]" methylamine~|~methanamine,"[""CH3"" => 1, ""NH2"" => 1]" ethylamine~|~ethanamine,"[""CH3"" => 1, ""CH2"" => 1, ""NH2"" => 1]" +hydronium,"[""H3O+""=>1]" sodium,"[""Na+""=>1]" +lithium,"[""Li+""=>1]" +potassium,"[""K+""=>1]" +rubidium,"[""Rb+""=>1]" +cesium,"[""Cs+""=>1]" +magnesium,"[""Mg2+""=>1]" +calcium,"[""Ca2+""=>1]" +strontium,"[""Sr2+""=>1]" +barium,"[""Ba2+""=>1]" +fluoride,"[""F-""=>1]" chloride,"[""Cl-""=>1]" +bromide,"[""Br-""=>1]" +iodide,"[""I-""=>1]" +hydroxide,"[""OH-""=>1]" acetate,"[""CH3""=>1,""COO-""=>1]" +acetic acid,"[""CH3""=>1,""COOH""=>1]" +bicarbonate,"[""HCO3-""=>1]" \ No newline at end of file diff --git a/database/properties/formation_aqueous.csv b/database/properties/formation_aqueous.csv index a839adf74..1d40e143b 100755 --- a/database/properties/formation_aqueous.csv +++ b/database/properties/formation_aqueous.csv @@ -6,4 +6,8 @@ hydronium,-157.244e3,-237.129e3 water,-157.244e3,-237.129e3 acetate,-486.01e3,-369.31e3 acetic acid,-485.76e3,-396.46e3 -sodium,-240.12e3,-261.905e3 \ No newline at end of file +sodium,-240.12e3,-261.905e3 +chloride,-167.159e3,-131.228e3 +carbon dioxide,-413.8e3,-385.98e3 +carbonate,-677.14e3,-527.81e3 +bicarbonate,-691.99e3,-586.77e3 \ No newline at end of file diff --git a/database/properties/molarmass_groups.csv b/database/properties/molarmass_groups.csv index ca3be54e1..b3961dbc0 100755 --- a/database/properties/molarmass_groups.csv +++ b/database/properties/molarmass_groups.csv @@ -85,6 +85,9 @@ N2,28.014 SF6,146.055 cO,16.00 OH (primary),17.008 +OH-,17.008 +H3O+,19.02 +HCO3-,61.01684 Na+,22.989 Cl-,35.453 K+,39.098 From b570a3c0b38be22e1402b4762b0627e3a508abdf Mon Sep 17 00:00:00 2001 From: pw0908 Date: Wed, 21 Aug 2024 16:47:22 -0700 Subject: [PATCH 17/18] Generalise ReactiveAq for multiple species --- .../property_solvers/reactive/reactive_aq.jl | 38 ++++++++++++++++--- 1 file changed, 32 insertions(+), 6 deletions(-) diff --git a/src/methods/property_solvers/reactive/reactive_aq.jl b/src/methods/property_solvers/reactive/reactive_aq.jl index feacd740f..f82c8688b 100644 --- a/src/methods/property_solvers/reactive/reactive_aq.jl +++ b/src/methods/property_solvers/reactive/reactive_aq.jl @@ -6,8 +6,10 @@ function equilibrate(model::ReactiveAqModel,p,T,n0;z0=nothing) z0 = x0_equilibrium_conditions(model,p,T,n0) end f!(F,z) = equilibrium_conditions(model,F,p,T,n0,exp10.(z),ν,Keq) - sol = Solvers.nlsolve(f!,log10.(z0),TrustRegion(Newton(), NWI()),NEqOptions(),ForwardDiff.Chunk{1}()) # Exists + sol = Solvers.nlsolve(f!,log10.(z0),TrustRegion(Newton(), NWI()),NEqOptions(),ForwardDiff.Chunk{length(model.reactions)}()) # Exists + println(sol) ξ = exp10.(Solvers.x_sol(sol)) + println(ξ) n = n0+ν*ξ return n end @@ -39,21 +41,45 @@ function equilibrium_conditions(model::ReactiveAqModel,F,p,T,n0,ξ,ν,Keq) n = n0+ν*ξ nreact = length(model.reactions) + nspecies = length(model.components) μmix = chemical_potential(model.eosmodel,p,T,n) - + idx_w = find_water_indx(model) for i in 1:nreact - idx_w = find_water_indx(model) + μref = zeros(nspecies) + + # 1. Charged Reactants + ireact_chrg = findall(x->ν[x,i]<0 && model.eosmodel.charge[x] != 0,1:nspecies) - zref = ones(length(model.components)).*1e-30 + zref = ones(nspecies).*1e-30 zref[idx_w] = 1. - zref[ν[:,i] .!= 0 .&& model.eosmodel.charge .!= 0] .= 0.01801528*1. + zref[ireact_chrg] .= 0.01801528*abs.(ν[ireact_chrg,i]) + zref ./= sum(zref) + + μref[ireact_chrg] = chemical_potential(model.eosmodel,p,T,zref)[ireact_chrg] + # 2. Charged Products + iprod_chrg = findall(x->ν[x,i]>0 && model.eosmodel.charge[x] != 0,1:nspecies) + zref = ones(nspecies).*1e-30 + zref[idx_w] = 1. + zref[iprod_chrg] .= 0.01801528*abs.(ν[iprod_chrg,i]) zref ./= sum(zref) - μref = chemical_potential(model.eosmodel,p,T,zref) + μref[iprod_chrg] = chemical_potential(model.eosmodel,p,T,zref)[iprod_chrg] + + # 3. Neutral species + ineutral = findall(x->model.eosmodel.charge[x] == 0,1:nspecies) + for j in ineutral + zref = ones(nspecies).*1e-30 + zref[j] = 0.01801528*abs(ν[j,i]) + zref[idx_w] = 1. + zref ./= sum(zref) + + μref[j] = chemical_potential(model.eosmodel,p,T,zref)[j] + end a = exp.((μmix .- μref)./R̄./T) + F[i] = sum(log.(a).*ν[:,i]) - log(Keq[i]) end return F From 1ecb7da570091fe0372604651c7d41feefaf2a6a Mon Sep 17 00:00:00 2001 From: pw0908 Date: Wed, 21 Aug 2024 16:50:33 -0700 Subject: [PATCH 18/18] Add reactive nb example --- examples/reactive.ipynb | 405 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 405 insertions(+) create mode 100644 examples/reactive.ipynb diff --git a/examples/reactive.ipynb b/examples/reactive.ipynb new file mode 100644 index 000000000..bfe7b86c6 --- /dev/null +++ b/examples/reactive.ipynb @@ -0,0 +1,405 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\u001b[32m\u001b[1m Activating\u001b[22m\u001b[39m project at `~/Library/CloudStorage/OneDrive-CaliforniaInstituteofTechnology/University/UROP/SAFT_codes/Clapeyron`\n" + ] + } + ], + "source": [ + "using Pkg, Revise\n", + "Pkg.activate(\".\")\n", + "using Clapeyron, GCIdentifier, ChemicalIdentifiers\n", + "using Plots" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "true" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "revise(Clapeyron)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Neutral Systems" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "species = [\"ethanol\",\"acetic acid\",\"ethyl acetate\",\"water\"]\n", + "components = get_groups_from_name.(species,UNIFAC)\n", + "reactions = [(\"ethanol\"=>-1,\"acetic acid\"=>-1,\"ethyl acetate\"=>1,\"water\"=>1)]\n", + "\n", + "eosmodel = UNIFAC(components)\n", + "model = ReactiveModel(species,reactions; model = eosmodel)\n", + "\n", + "T = LinRange(273.15,373.15,100)\n", + "ξ = zeros(100)\n", + "z01 = 4e-1\n", + "z02 = 4e-1\n", + "n0 = [0.5,0.5,0.0,0.0]\n", + "for i in 1:100\n", + " ξ[i] = equilibrate(model,1e5,T[i],n0; z0=[z01])[3]\n", + " z01 = ξ[i]\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt = plot(grid=:off,framestyle=:box,foreground_color_legend = nothing,legend_font=font(12), guidefontsize = 12, tickfontsize = 9)\n", + "plot!(plt,T,ξ/0.5,label=\"UNIFAC\",linewidth=2)\n", + "# plot!(plt,T,ξ[:,2]/0.5,label=\"PC-SAFT\",linewidth=2)\n", + "xlims!(plt,273.15,333.15)\n", + "xlabel!(plt,\"Temperature (K)\")\n", + "ylabel!(plt,\"Conversion\")\n", + "title!(plt,\"Ethanol + Acetic Acid ↔ Ethyl Acetate + Water\")\n", + "savefig(plt,\"ethanol_acetic_acid_ethyl_acetate_water.pdf\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Aqueous Systems" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "species = [\"water\",\"hydronium\",\"hydroxide\"]\n", + "reactions = [(\"water\"=>-2,\"hydronium\"=>1,\"hydroxide\"=>1)]\n", + "\n", + "eosmodel = SAFTVREMie([\"water\"],[\"hydronium\",\"hydroxide\"])\n", + "model = ReactiveAqModel(species,reactions; model = eosmodel)\n", + "\n", + "n0 = [1.0,0.0,0.0]\n", + "\n", + "T = LinRange(273.15,333.15,100)\n", + "pH = zeros(100)\n", + "z0 = 1e-7\n", + "for i in 1:100\n", + " pH[i] = Clapeyron.pH(model,1e5,T[i],n0; z0=[z0])\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt = plot(grid=:off,framestyle=:box,foreground_color_legend = nothing,legend_font=font(12), guidefontsize = 12, tickfontsize = 9)\n", + "plot!(plt,T,pH,label=\"SAFT-VRE Mie\",linewidth=2)\n", + "# plot!(plt,T,ξ[:,2]/0.5,label=\"PC-SAFT\",linewidth=2)\n", + "xlims!(plt,273.15,333.15)\n", + "xlabel!(plt,\"Temperature (K)\")\n", + "ylabel!(plt,\"pH\")\n", + "title!(plt,\"pH of pure water\")\n", + "# savefig(plt,\"pH.pdf\")\n", + "plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "species = [\"water\",\"hydronium\",\"hydroxide\",\"sodium\"]\n", + "reactions = [(\"water\"=>-2,\"hydronium\"=>1,\"hydroxide\"=>1)]\n", + "\n", + "eosmodel = SAFTVREMie([\"water\"],[\"hydronium\",\"hydroxide\",\"sodium\"])\n", + "model = ReactiveAqModel(species,reactions; model = eosmodel)\n", + "\n", + "salts = [(\"sodium hydroxide\",[\"sodium\"=>1,\"hydroxide\"=>1]),\n", + " (\"h3o.oh\",[\"hydronium\"=>1,\"hydroxide\"=>1])]\n", + "\n", + "m = 10 .^LinRange(-7,-1,100)\n", + "T = 298.15\n", + "\n", + "pH = zeros(100)\n", + "i = 1\n", + "\n", + "for i in 1:100\n", + " z0 = molality_to_composition(model.eosmodel, salts, [m[i],0.0])\n", + "\n", + " pH[i] = Clapeyron.pH(model,1e5,T,z0; z0=[1e-7])\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt = plot(grid=:off,framestyle=:box,foreground_color_legend = nothing,legend_font=font(12), guidefontsize = 12, tickfontsize = 9, xaxis=:log)\n", + "plot!(plt,m,pH,label=\"SAFT-VRE Mie\",linewidth=2)\n", + "xlabel!(plt,\"NaOH molality (mol/kg)\")\n", + "ylabel!(plt,\"pH\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "species = [\"water\",\"hydronium\",\"hydroxide\",\"chloride\"]\n", + "reactions = [(\"water\"=>-2,\"hydronium\"=>1,\"hydroxide\"=>1)]\n", + "\n", + "eosmodel = SAFTVREMie([\"water\"],[\"hydronium\",\"hydroxide\",\"chloride\"])\n", + "model = ReactiveAqModel(species,reactions; model = eosmodel)\n", + "\n", + "salts = [(\"hydrogen chloride\",[\"hydronium\"=>1,\"chloride\"=>1]),\n", + " (\"h3o.oh\",[\"hydronium\"=>1,\"hydroxide\"=>1])]\n", + "\n", + "m = 10 .^LinRange(-7,-1,100)\n", + "T = 298.15\n", + "\n", + "pH = zeros(100)\n", + "i = 1\n", + "\n", + "for i in 1:100\n", + " z0 = molality_to_composition(model.eosmodel, salts, [m[i],0.0])\n", + "\n", + " pH[i] = Clapeyron.pH(model,1e5,T,z0; z0=[1e-7])\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt = plot(grid=:off,framestyle=:box,foreground_color_legend = nothing,legend_font=font(12), guidefontsize = 12, tickfontsize = 9, xaxis=:log)\n", + "plot!(plt,m,pH,label=\"SAFT-VRE Mie\",linewidth=2)\n", + "xlabel!(plt,\"HCl molality (mol/kg)\")\n", + "ylabel!(plt,\"pH\")" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[1.0110054626220328e-14 1.7519664617156773e-5]\n", + "Results of solving non-linear equations\n", + "\n", + "* Algorithm:\n", + " Newton's method with default linsolve with Trust Region (Newton, eigen)\n", + "\n", + "* Candidate solution:\n", + " Final residual 2-norm: 2.61e-13\n", + " Final residual Inf-norm: 2.27e-13\n", + "\n", + " Initial residual 2-norm: 5.37e+01\n", + " Initial residual Inf-norm: 5.01e+01\n", + "\n", + "* Stopping criteria\n", + " |F(x')| = 2.27e-13 <= 1.00e-08 (true)\n", + "\n", + "* Work counters\n", + " Seconds run: 4.37e-02\n", + " Iterations: 3\n", + "\n" + ] + }, + { + "data": { + "text/plain": [ + "-1.1873918341936291" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "species = [\"water\",\"acetic acid\",\"hydronium\",\"hydroxide\",\"acetate\"]\n", + "reactions = [(\"water\"=>-2,\"hydronium\"=>1,\"hydroxide\"=>1),\n", + " (\"water\"=>-1,\"acetic acid\"=>-1,\"acetate\"=>1,\"hydronium\"=>1)]\n", + "\n", + "eosmodel = SAFTgammaEMie([\"water\",\"acetic acid\"],[\"hydronium\",\"hydroxide\",\"acetate\"])\n", + "model = ReactiveAqModel(species,reactions; model = eosmodel)\n", + "\n", + "salts = [(\"acetic acid\",[\"acetate\"=>1,\"hydronium\"=>1]),\n", + " (\"h3o.oh\",[\"hydronium\"=>1,\"hydroxide\"=>1])]\n", + "\n", + "z0 =[0.9,0.1,0.0,0.0,0.0]\n", + "\n", + "# z0 = molality_to_composition(model.eosmodel, salts, [0.0,1e-14,0.0], [0.98,0.02])\n", + "\n", + "Clapeyron.pH(model,1e5,298.15,z0; z0=[1e-7,1e-4])" + ] + }, + { + "cell_type": "code", + "execution_count": 49, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[1.0110054626220328e-14 4.302121165288985e-7]\n", + "Results of solving non-linear equations\n", + "\n", + "* Algorithm:\n", + " Newton's method with default linsolve with Trust Region (Newton, eigen)\n", + "\n", + "* Candidate solution:\n", + " Final residual 2-norm: 4.62e+00\n", + " Final residual Inf-norm: 4.14e+00\n", + "\n", + " Initial residual 2-norm: 1.11e+01\n", + " Initial residual Inf-norm: 9.98e+00\n", + "\n", + "* Stopping criteria\n", + " |F(x')| = 4.14e+00 <= 1.00e-08 (false)\n", + "\n", + "* Work counters\n", + " Seconds run: 1.66e+00\n", + " Iterations: 43\n", + "\n", + "[7.589023646643154e-43, 5.826776409944142e-7]\n" + ] + }, + { + "data": { + "text/plain": [ + "6-element Vector{Float64}:\n", + " 0.998998834444918\n", + " 0.0009994173221590057\n", + " 5.826776409944142e-7\n", + " 9.999999998e-11\n", + " 5.826776409944142e-7\n", + " 9.999999998e-11" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "species = [\"water\",\"carbon dioxide\",\"hydronium\",\"hydroxide\",\"bicarbonate\",\"sodium\"]\n", + "reactions = [(\"water\"=>-2,\"hydronium\"=>1,\"hydroxide\"=>1),\n", + " (\"water\"=>-2,\"carbon dioxide\"=>-1,\"bicarbonate\"=>1,\"hydronium\"=>1)]\n", + "\n", + "eosmodel = SAFTgammaEMie([\"water\",\"carbon dioxide\"],[\"hydronium\",\"hydroxide\",\"bicarbonate\",\"sodium\"])\n", + "model = ReactiveAqModel(species,reactions; model = eosmodel)\n", + "\n", + "salts = [(\"carbon dioxide\",[\"bicarbonate\"=>1,\"hydronium\"=>1]),\n", + " (\"sodium hydroxide\",[\"sodium\"=>1,\"hydroxide\"=>1]),\n", + " (\"h3o.oh\",[\"hydronium\"=>1,\"hydroxide\"=>1])]\n", + "\n", + "z0 =[0.999,0.001,0.0,1e-10,0.0,1e-10]\n", + "z0 /= sum(z0)\n", + "\n", + "# z0 = molality_to_composition(model.eosmodel, salts, [0.0,1e-5,0.0], [0.99,0.01])\n", + "\n", + "z = Clapeyron.equilibrate(model,1e5,298.15,z0; z0=[1e-7,1e-7])" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "6×2 Matrix{Float64}:\n", + " -2.0 -2.0\n", + " 0.0 -1.0\n", + " 1.0 1.0\n", + " 1.0 0.0\n", + " 0.0 1.0\n", + " 0.0 0.0" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "Clapeyron.stoichiometric_coefficient(model)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "species = [\"water\",\"hydronium\",\"hydroxide\"]\n", + "reactions = [(\"water\"=>-2,\"hydronium\"=>1,\"hydroxide\"=>1)]\n", + "\n", + "eosmodel = SAFTVREMie([\"water\"],[\"hydronium\",\"hydroxide\"])\n", + "model = ReactiveAqModel(species,reactions; model = eosmodel)\n", + "\n", + "Clapeyron.equilibrate(model,1e5,298.15,[1.0,0.0,0.0]; z0=[1e-7])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.9.4", + "language": "julia", + "name": "julia-1.9" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.9.4" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}