diff --git a/docs/src/tutorials/background.md b/docs/src/tutorials/background.md index da32c1a..f3bf3e6 100644 --- a/docs/src/tutorials/background.md +++ b/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 diff --git a/docs/src/tutorials/getting_started.md b/docs/src/tutorials/getting_started.md index 0ea46e6..64ca098 100644 --- a/docs/src/tutorials/getting_started.md +++ b/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. diff --git a/docs/src/tutorials/tutorial.md b/docs/src/tutorials/tutorial.md index 4974364..a207d49 100644 --- a/docs/src/tutorials/tutorial.md +++ b/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. diff --git a/src/methods/isotherm_methods.jl b/src/methods/isotherm_methods.jl index a2d0a5c..3c8f4bf 100644 --- a/src/methods/isotherm_methods.jl +++ b/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. diff --git a/src/models/models.jl b/src/models/models.jl index 874dbca..1de0d68 100644 --- a/src/models/models.jl +++ b/src/models/models.jl @@ -292,3 +292,4 @@ include("unilan.jl") include("Toth.jl") include("multisite.jl") include("interpolation.jl") +include("thermodynamic_langmuir.jl") diff --git a/src/models/thermodynamic_langmuir.jl b/src/models/thermodynamic_langmuir.jl new file mode 100644 index 0000000..b71c659 --- /dev/null +++ b/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 + diff --git a/src/utils.jl b/src/utils.jl index 7f93c4f..8632d5f 100644 --- a/src/utils.jl +++ b/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 \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 936b01a..f90c4ad 100644 --- a/test/runtests.jl +++ b/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 + + +