Mockup: Isenthalpic Flash

src/eos.jl

@@ -326,6 +326,23 @@ 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  -> 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)
+        # 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
 """

src/isoenthalpic-flash/enthalpy_flash_util.jl

@@ -0,0 +1,67 @@
+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  -> 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')
+end
+
+
+function enthalpy_departure(eos, cond, i, Z, forces, scalars) 
+    # Input: T   -> Current Temperature, 
+    #        T_r -> Reference Temperature 
+    # 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  -> 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
+

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

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...)