From 2a0e36a9f7cd9f6751b3d098560e0ef789a2d234 Mon Sep 17 00:00:00 2001 From: Castiel Date: Tue, 10 Dec 2024 11:22:23 +0100 Subject: [PATCH 1/2] mockup version that performs flash that calls V = flash_2ph(eos, conditions) and also compute the phase enthalpy --- src/eos.jl | 15 ++++++ src/isoenthalpic-flash/enthalpy_flash_util.jl | 47 +++++++++++++++++++ src/isoenthalpic-flash/isoenthalpic-script.jl | 44 +++++++++++++++++ src/stability.jl | 2 + 4 files changed, 108 insertions(+) create mode 100644 src/isoenthalpic-flash/enthalpy_flash_util.jl create mode 100644 src/isoenthalpic-flash/isoenthalpic-script.jl diff --git a/src/eos.jl b/src/eos.jl index 4d54857..9e0c5b0 100644 --- a/src/eos.jl +++ b/src/eos.jl @@ -326,6 +326,21 @@ function mixture_fugacities!(f, eos, cond, forces = force_coefficients(eos, cond return f end +## +function mixture_enthalpy!(H, eos, cond, forces = force_coefficients(eos, cond), scalars = force_scalars(eos, cond, forces)) + Z = mixture_compressibility_factor(eos, cond, forces, scalars) + # Cp -> Array with 4 Heat Capacity Coefficients for calculating + # Ex: Cp = [1 1 1 1] coeficients of the polynomial used to compute idealeEnthalpy of the component + # Hp phase enthalpy entha vector 1x n_phases + # for each phase + @inbounds for i in eachindex(H) + # Phase Enthalpy is the + H[i] = enthalpy_phase(eos, cond, i, Z, forces, scalars, T,T_r, Cp); + end + return H +end + + """ Specialization of solve_roots for cubics """ diff --git a/src/isoenthalpic-flash/enthalpy_flash_util.jl b/src/isoenthalpic-flash/enthalpy_flash_util.jl new file mode 100644 index 0000000..3d07d2d --- /dev/null +++ b/src/isoenthalpic-flash/enthalpy_flash_util.jl @@ -0,0 +1,47 @@ +function secant_update(T_curr, T_prev, R_curr, R_prev) + return T_curr - R_curr*(T_curr - T_prev)/(R_curr - R_prev) +end + +function secant_intialization(T_0, R_0, F_L, F_V, Cp_L, Cp_V) + return T_0 - (R_0/(F_L*Cp_L + F_V*Cp_V)) +end + +function regular_falsi_update(T_left, T_right, R_curr, R_left, R_right) + return T_left - R_curr*(T_left-T_right)/(R_left - R_right) +end + + +function enthalpy_ideal(T,T_r, Cp) +# Input: T -> Current Temperature, +# T_r -> Reference Temperature +# Cp -> Array with 4 Heat Capacity Coefficients +# # Ex: Cp = [1 1 1 1] coeficients of the polynomial used to compute idealeEnthalpy of the component +# Output: Ideal Component Enthalpy + p_exp = collect(1:4) + return sum((fill(T,4,).^p_exp - fill(T_r,4,).^p_exp ).*Cp') +end + +function enthalpy_departure(eos, cond, i, Z, forces, scalars) + # Input: T -> Current Temperature, + # T_r -> Reference Temperature + # Cp -> Array with 4 Heat Capacity Coefficients + # Output: Departure Enthalpy + # for each component + # Hc -> Array n_c components x 1 + R = 8.314 # J / mol·K + @inbounds for i in eachindex(c) + # d_component_fugacity_dt = d/dT(component_fugacity) with AD + Hc[i] = d_component_fugacity_dt(eos, cond, i, Z, forces, scalars) + end + return -sum(H_c)*R*(T^2) +end + +function enthalpy_phase(eos, cond, i, Z, forces, scalars, T,T_r, Cp) + # Input: eos -> Eqution of State + # T -> Current Temperature, + # T_r -> Reference Temperature + # Cp -> Array with 4 Heat Capacity Coefficients + # Output: Phase Enthalpy + return enthalpy_ideal(T,T_r, Cp) + enthalpy_departure(eos, cond, i, Z, forces, scalars) +end + diff --git a/src/isoenthalpic-flash/isoenthalpic-script.jl b/src/isoenthalpic-flash/isoenthalpic-script.jl new file mode 100644 index 0000000..5e73f17 --- /dev/null +++ b/src/isoenthalpic-flash/isoenthalpic-script.jl @@ -0,0 +1,44 @@ +include("enthalpy_flash_util.jl") +using MultiComponentFlash +# Define two species: One heavy and one light. +# The heavy component uses table lookup: +decane = MolecularProperty("n-Decane") +# The light component is given with explicit properties +mw = 0.0160428 # Molar mass (kg/mole) +P_c = 4.5992e6 # Critical pressure (Pa) +T_c = 190.564 # Critical temperature (°K) +V_c = 9.4118e-5 # Critical volume (m^3/mole) +ω = 0.22394 # Acentric factor +methane = MolecularProperty(mw, P_c, T_c, V_c, ω) +# or, equivialently, +# methane = MolecularProperty("Methane") +# Create a mixture +mixture = MultiComponentMixture((methane, decane)) +eos = GenericCubicEOS(mixture, PengRobinson()) +# Define conditions to flash at +p = 5e6 # 5 000 000 Pa, or 50 bar +T = 303.15 # 30 °C = 303.15 °K +z = [0.4, 0.6] # 1 mole methane per 9 moles of decane +conditions = (p = p, T = T, z = z) +# Perform a flash to get the vapor fraction +V = flash_2ph(eos, conditions) +cp = [1 1 1 1] + +# TT = cp + cp +# ref +# Tref = cp + +qq = enthalpy_phase(eos,2,1,cp) + + +# let index = 1 +# while index < 300 +# println("Counter is: $index") +# index += 1 +# end +# end + +# while index1 < 1000 +# println("Counter is: $index1") +# index1 += 1 +# end diff --git a/src/stability.jl b/src/stability.jl index fa922b3..5061cd0 100644 --- a/src/stability.jl +++ b/src/stability.jl @@ -28,6 +28,8 @@ function stability_2ph!(storage, K, eos, c; vapor = (p = p, T = T, z = y) # Update fugacities for current conditions used in both tests mixture_fugacities!(f_z, eos, c, forces) + # + mixture_enthalpy!(h_z, eos, c, forces) if check_vapor wilson_estimate!(K, eos, p, T) v = michelsen_test!(vapor, f_z, f_xy, vapor.z, z, K, eos, c, forces, Val(true); kwarg...) From c19d55f139d0958b7fc20d7c6f42b7e97deee8d9 Mon Sep 17 00:00:00 2001 From: Artur Castiel Date: Tue, 10 Dec 2024 14:22:26 +0100 Subject: [PATCH 2/2] updated the definition of Cp --- src/eos.jl | 6 ++-- src/isoenthalpic-flash/enthalpy_flash_util.jl | 30 +++++++++++++++---- 2 files changed, 29 insertions(+), 7 deletions(-) diff --git a/src/eos.jl b/src/eos.jl index 9e0c5b0..de77216 100644 --- a/src/eos.jl +++ b/src/eos.jl @@ -329,8 +329,10 @@ end ## function mixture_enthalpy!(H, eos, cond, forces = force_coefficients(eos, cond), scalars = force_scalars(eos, cond, forces)) Z = mixture_compressibility_factor(eos, cond, forces, scalars) - # Cp -> Array with 4 Heat Capacity Coefficients for calculating - # Ex: Cp = [1 1 1 1] coeficients of the polynomial used to compute idealeEnthalpy of the component +# Cp -> Matrix nc x 4. +# nc = number of components +# Each line has 4 coeficients for calculating the ideal enthalpy of a component +# # Hp phase enthalpy entha vector 1x n_phases # for each phase @inbounds for i in eachindex(H) diff --git a/src/isoenthalpic-flash/enthalpy_flash_util.jl b/src/isoenthalpic-flash/enthalpy_flash_util.jl index 3d07d2d..c426131 100644 --- a/src/isoenthalpic-flash/enthalpy_flash_util.jl +++ b/src/isoenthalpic-flash/enthalpy_flash_util.jl @@ -14,17 +14,37 @@ end function enthalpy_ideal(T,T_r, Cp) # Input: T -> Current Temperature, # T_r -> Reference Temperature -# Cp -> Array with 4 Heat Capacity Coefficients -# # Ex: Cp = [1 1 1 1] coeficients of the polynomial used to compute idealeEnthalpy of the component +# Cp -> Matrix nc x 4. The coeficients +# nc = number of components +# each component needs with 4 Heat Capacity Coefficients + +nc, _ = size(Cp) +h_c = 0.0; + +# Output: Ideal Mixture Enthalpy + p_exp = collect(1:4) + for i = 0:nc + # x_c -> fraction of c in a given phase + h_c = h_c + (x_c) * enthalpy_ideal_component(T,T_r, Cp[i,:]); + end + return hc +end + +function enthalpy_ideal_component(T,T_r, cp) +# Input: T -> Current Temperature, +# T_r -> Reference Temperature +# Cp -> Matrix 1 x 4. The coeficients +# nc = number of components +# each component needs with 4 Heat Capacity Coefficients # Output: Ideal Component Enthalpy p_exp = collect(1:4) - return sum((fill(T,4,).^p_exp - fill(T_r,4,).^p_exp ).*Cp') + return sum((fill(T,4,).^p_exp - fill(T_r,4,).^p_exp ).*cp') end + function enthalpy_departure(eos, cond, i, Z, forces, scalars) # Input: T -> Current Temperature, # T_r -> Reference Temperature - # Cp -> Array with 4 Heat Capacity Coefficients # Output: Departure Enthalpy # for each component # Hc -> Array n_c components x 1 @@ -40,7 +60,7 @@ function enthalpy_phase(eos, cond, i, Z, forces, scalars, T,T_r, Cp) # Input: eos -> Eqution of State # T -> Current Temperature, # T_r -> Reference Temperature - # Cp -> Array with 4 Heat Capacity Coefficients + # Cp -> Matrix with 4 Heat Capacity Coefficients # Output: Phase Enthalpy return enthalpy_ideal(T,T_r, Cp) + enthalpy_departure(eos, cond, i, Z, forces, scalars) end