Merge pull request #21 from ClapeyronThermo/vini/thermolangmuirmodel

Vini/thermolangmuirmodel

docs/src/tutorials/background.md

@@ -14,13 +14,13 @@ The first approach focuses on kinetics, where adsorption and desorption rates ar
 
 The heat of adsorption is a critical design parameter in adsorptive gas separation units. During adsorption, heat is released as adsorbate molecules transition to a lower energy state on the surface of the adsorbent compared to their higher energy state in the bulk gas phase. This exothermic process significantly impacts both the efficiency and operational conditions of adsorption systems. For a single component, the isosteric heat is given by:
 
-$Q_{st} = -T*(V_g - V_a)*\left( \frac{dP_i}{dT} \right)\rvert_{(N_i,A)}$ 
+$Q_{st} = T*(V_g - V_a)*\left( \frac{dP_i}{dT} \right)\rvert_{(N_i,A)}$ 
 
 where $Q_st$ is the isosteric heat of the component being adsorbed, $T$ is the temperature, $V_g$ is the molar volume of the component in gas phase, $V_a$ is the molar volume of the component in adsorbed phase, $N_i$ is the amount of component adsorbed of the component
 
 When the isotherm is of the form $N_i = f(T, P_i)$, one can write:
 
-$Q_{st, i} = -T \cdot (V_g - V_a) \cdot \left( \frac{\partial N_i / \partial T \mid P_i}{\partial N_i / \partial P \mid T} \right)$
+$Q_{st, i} = T \cdot (V_g - V_a) \cdot \left( \frac{\partial N_i / \partial T \mid P_i}{\partial N_i / \partial P \mid T} \right)$
 
 
 ## Multi component adsorption

docs/src/tutorials/getting_started.md

@@ -61,7 +61,7 @@ Below it is assumed that the ideal gas law is a good approximation to describe t
 import Langmuir: Rgas
 ΔH = map(P -> isosteric_heat(isotherm, P, 300.), P[2:end]) |> x -> round.(x, digits = 7)
 scatter(l_at_300[2:end], ΔH, size = (500, 250),  ylabel = "Isosteric heat (J/mol)", xlabel = "loading (mol/kg)", label = "Estimated isosteric heat with AD")
-plot!([first(l_at_300), last(l_at_300)], [-E, -E], label = "Expected value") 
+plot!([first(l_at_300), last(l_at_300)], [E, E], label = "Expected value") 
 ```
 
 ## Estimating properties in multicomponent adsorption.

docs/src/tutorials/tutorial.md

@@ -108,8 +108,8 @@ P_C₂₌_283K = sort(lp_ethylene[1][2][2:end])
 l_C₂₌_283K = sort(loading1_ethylene[2:end])
 ΔH_C₂₌_283K = map(p -> isosteric_heat(ethylene_isotherm, p, 283.0), P_C₂₌_283K)
 
-scatter(l_C₂_283K, -ΔH_C₂_283K, xlabel = "Loading (mol/kg)", ylabel = "Isosteric Heat (J/mol)", m = (4, :white, stroke(1, :lightslateblue)), markershape = :circle, label = "Ethane - 283.0 K", size = (600, 300))
-scatter!(l_C₂₌_283K, -ΔH_C₂₌_283K, xlabel = "Loading (mol/kg)", ylabel = "Isosteric Heat (J/mol)", markershape = :square, m = (3, :white, stroke(1, :springgreen2)), label = "Ethylene - 283.0 K")
+scatter(l_C₂_283K, ΔH_C₂_283K, xlabel = "Loading (mol/kg)", ylabel = "Isosteric Heat (J/mol)", m = (4, :white, stroke(1, :lightslateblue)), markershape = :circle, label = "Ethane - 283.0 K", size = (600, 300))
+scatter!(l_C₂₌_283K, ΔH_C₂₌_283K, xlabel = "Loading (mol/kg)", ylabel = "Isosteric Heat (J/mol)", markershape = :square, m = (3, :white, stroke(1, :springgreen2)), label = "Ethylene - 283.0 K")
 ```
 The plot shows the isosteric heat of adsorption for ethane and ethylene at 283.0 K as a function of loading. Ethane exhibits a steeper decline in adsorption heat, suggesting stronger initial interactions that weaken significantly as loading increases, whereas ethylene's decline is more gradual, indicating a slower reduction in adsorption strength. 
 

src/methods/isotherm_methods.jl

@@ -220,7 +220,7 @@ function isosteric_heat(model::IsothermModel, p, T; Vᵃ = zero(eltype(p)), Vᵍ
 
     ∂n_∂p, ∂n_∂T = _df
 
-    return -T*(Vᵍ - Vᵃ)*∂n_∂T/∂n_∂p
+    return T*(Vᵍ - Vᵃ)*∂n_∂T/∂n_∂p
 end
 
 #useful for creating pseudo langmuir models for multicomponent adsoption.

src/models/models.jl

@@ -292,3 +292,4 @@ include("unilan.jl")
 include("Toth.jl")
 include("multisite.jl")
 include("interpolation.jl")
+include("thermodynamic_langmuir.jl")

src/models/thermodynamic_langmuir.jl

@@ -0,0 +1,104 @@
+@with_metadata struct ThermodynamicLangmuir{T} <: IsothermModel{T}
+    (M::T, (0.0, Inf), "saturation loading")
+    (K₀::T, (0.0, Inf), "affinity parameter") 
+    (E::T, (-Inf, 0.0), "energy parameter")
+    (Bᵢᵩ::T, (-Inf, 0.0), "adsorbate-adsorbent interaction coefficient")
+end
+
+
+function gibbs_excess_free_energy(model::ThermodynamicLangmuir, T, θ)
+
+    _1 = one(eltype(T))
+    Bᵢᵩ = model.Bᵢᵩ
+    T⁻¹ = _1/T
+    θᵢ, θᵩ = θ 
+
+    τᵢᵩ = Bᵢᵩ*T⁻¹
+    Gᵢᵩ = exp(-0.3*τᵢᵩ)
+
+    gᴱ_RT = θᵢ*θᵩ*τᵢᵩ*(Gᵢᵩ - _1)/(θᵢ*Gᵢᵩ + θᵩ)
+
+    return gᴱ_RT
+end
+
+function activity_coefficient(model::ThermodynamicLangmuir, T, θ)
+
+         fun(θ) = gibbs_excess_free_energy(model, T, θ)
+
+         return exp.(gradient(fun, θ))
+end
+
+function loading_x0(model::ThermodynamicLangmuir, p, T)    
+
+    guess_model = from_vec(LangmuirS1, [model.M, model.K₀, model.E])
+
+    return loading(guess_model, p, T)
+end
+
+function isotherm_res(model::ThermodynamicLangmuir, q, p, T)
+    _1 = one(eltype(model))
+    M = model.M
+    K₀ = model.K₀
+    E = model.E
+    K = K₀*exp(-E/(Rgas(model)*T))
+    θᵢ = q/M
+    γᵢ, γᵩ = activity_coefficient(model, T, [θᵢ, _1 - θᵢ])
+    return q - M * K *p / (γᵢ/γᵩ + K*p) 
+end
+
+function henry_coefficient(model::ThermodynamicLangmuir, T)
+
+    _0 = zero(eltype(T))
+
+    q_0 = loading(model, _0, T)
+
+    f(∂q, ∂p) = isotherm_res(model, ∂q, ∂p, T)
+
+    _f,_df = fgradf2(f, q_0, _0)
+
+    ∂f0_∂q, ∂f0_∂P = _df
+
+    return -∂f0_∂P/∂f0_∂q
+end
+
+function isosteric_heat(model::ThermodynamicLangmuir, p, T; Vᵃ = zero(eltype(p)), Vᵍ = Rgas(model)*T/p)
+
+    q = loading(model, p, T)
+
+    f(∂p, ∂T) = isotherm_res(model, q, ∂p, ∂T)
+
+    _f,_df = fgradf2(f, p, T)
+
+    ∂n_∂p, ∂n_∂T = _df
+
+    return T*(Vᵍ - Vᵃ)*∂n_∂T/∂n_∂p
+
+end
+
+function loading_impl(model::ThermodynamicLangmuir, p, T)
+
+    _1 = one(eltype(p))
+    q0 = loading_x0(model, p, T)
+    M = model.M
+    K₀ = model.K₀
+    E = model.E
+    K = K₀*exp(-E/(Rgas(model)*T))
+
+    f0(q, P_T) = begin 
+    p, T = P_T
+    θᵢ = q/M
+    γᵢ, γᵩ = activity_coefficient(model, T, [θᵢ, _1 - θᵢ])
+    q - M * K *p / (γᵢ/γᵩ + K*p)    
+    end 
+    prob = Roots.ZeroProblem(f0, q0)
+    return Roots.solve(prob, p = (p, T))
+
+end
+
+
+function loading(model::ThermodynamicLangmuir, p, T)
+    return loading_impl(model, p, T)
+end
+
+export ThermodynamicLangmuir
+

src/utils.jl

@@ -126,4 +126,8 @@ function _eltype(x::T) where T
     else
         return eltyp
     end
+end
+
+@inline function gradient(f::F, x) where {F}
+    return ForwardDiff.gradient(f,x)
 end

test/runtests.jl

@@ -1,8 +1,12 @@
 using Langmuir
 import Langmuir: loading_ad, sp_res, to_vec, sp_res_numerical, isosteric_heat, Rgas, from_vec, fit, pressure, temperature, x0_guess_fit
+import Langmuir: gibbs_excess_free_energy, activity_coefficient
 import Langmuir: IsothermFittingProblem, DEIsothermFittingSolver
 using Test
 const LG = Langmuir
+
+
+
 #we test that definitions of loading and sp_res are consistent.
 
 function test_sp_res_loading(model, prange, T)
@@ -26,7 +30,7 @@ end
         test_sp_res_loading(cL_test, [1e5, 2e5, 3e5], 298.)
 
         # Isosteric heat
-        ΔH = isosteric_heat(cL_test, 101325.0, 298.) ≈ -cL_test.E
+        ΔH = isosteric_heat(cL_test, 101325.0, 298.) ≈ cL_test.E
     end
 end
 
@@ -73,7 +77,16 @@ end
     end
 end
 
-
+@testset "ThermodynamicLangmuir" begin
+    tlang = ThermodynamicLangmuir(2.468, 7.03e-10, -2540*Rgas(LangmuirS1{Float64}), -611.63)
+    analyical_gammas(T, θ) = (exp(θ[2]^2*tlang.Bᵢᵩ/T*(exp(-0.3*tlang.Bᵢᵩ/T) - 1)/(θ[1]*exp(-0.3*tlang.Bᵢᵩ/T) + θ[2])^2),
+    exp(θ[1]^2*-tlang.Bᵢᵩ/T*(exp(-0.3*-tlang.Bᵢᵩ/T) - 1)/(θ[1] + θ[2]*exp(-0.3*-tlang.Bᵢᵩ/T))^2)                                                                                         )
+    @test analyical_gammas(273.15, [0.4, 0.6])[1] ≈ activity_coefficient(tlang, 273.15, [0.4, 0.6])[1]
+    @test analyical_gammas(273.15, [0.4, 0.6])[2] ≈ activity_coefficient(tlang, 273.15, [0.4, 0.6])[2]
+
+    tlang2 = ThermodynamicLangmuir(2.468, 7.03e-10, -2540*Rgas(LangmuirS1{Float64}), 0.0)
+    @test loading(tlang2, 101325.0, 298.15) ≈ loading(LangmuirS1(tlang2.M, tlang2.K₀, tlang2.E), 101325.0, 298.15)
+end
 
 
 @testset "IAST" begin
@@ -116,3 +129,6 @@ end
     @test iast(models,p,T,y,IASTNestedLoop())[1] ≈ q_tot
     @test iast(models,p,T,y,IASTNestedLoop(),x0 = x0)[1] ≈ q_tot
 end
+
+
+