Avi/streams cstr

.gitignore

@@ -3,4 +3,5 @@
 *.jl.mem
 /docs/Manifest.toml
 /docs/build/
-Manifest.toml
+Manifest.toml
+.vscode/**/*

.vscode/settings.json

@@ -0,0 +1 @@
+{}

Project.toml

@@ -4,7 +4,18 @@ authors = ["Avinash Subramanian"]
 version = "1.0.0-DEV"
 
 [deps]
-SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Clapeyron = "7c7805af-46cc-48c9-995b-ed0ed2dc909a"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+GCIdentifier = "b7ea765e-cbac-4e4a-9b0d-5427cc302506"
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
+UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
 julia = "1.10"

src/ProcessSimulator.jl

@@ -1,7 +1,18 @@
 module ProcessSimulator
 
-# Write your package code here.
-export Gibbs
-include("Reactors/Gibbs.jl")
 
+using ModelingToolkit, JSON , DifferentialEquations
+using ModelingToolkit: t_nounits as t, D_nounits as D
+import ModelingToolkit: scalarize, equations, get_unknowns, defaults
+using Clapeyron, Symbolics
+
+
+export load_component_properties, read_reidcp, my_model, enthalpy_simple, molar_density_simple, MaterialSource, SimpleCSTR, my_model, KineticReactionNetwork, Display
+include("utils")
+include("Reactors/ReactionManager/KineticReaction.jl")
+include("Reactors/SimplifiedCSTR.jl")
+
+export MaterialSource
+include("Sources/MaterialSource.jl")
+include("Sources/Sourceutils.jl")
 end

src/Reactors/Gibbs_Reactor/GibbsReactor_numeric.jl

@@ -0,0 +1,164 @@
+using Pkg
+Pkg.activate("GibbsReactor")
+
+using ModelingToolkit
+using Clapeyron
+using GCIdentifier
+using DelimitedFiles
+using LinearAlgebra
+
+#define components
+components_names = ["carbon dioxide", "carbon monoxide", "hydrogen", "water", "methane"]
+components_smiles = ["C(=O)=O", "[C-]#[O+]", "[HH]", "O=O", "C"]
+
+
+#Define the model
+my_model = ShomateIdeal(components_names;
+ userlocations = ["C:/Users/Vinic/OneDrive/Pos-Doc/Flowsheeting/GibbsReactor/GibbsReactor.jl", "myShomate.csv"])
+
+function readEnthalpyEntropyData(file_path, component_names)
+    # Read formation enthalpy from file
+    # file_path: File path of the file containing the formation enthalpy
+    # component_names: Names of the components
+    # Returns: Formation enthalpy in J/mol
+
+    # initialize the formation enthalpy and entropy vector
+    H₀_S₀ = zeros(length(component_names), 2)
+
+    # Read the TSV file
+    data = readdlm(file_path, '\t', header = true)
+
+    # Get the column of chemical names
+    chemical_names = data[1][:, 2]
+
+    # Find the index of the chemical name in the column
+    for i in eachindex(component_names)
+        index = findfirst(x -> x == component_names[i], chemical_names)
+
+        # If the chemical name was found, return the corresponding line
+        if index !== nothing
+            H₀_S₀[i, :] =  data[1][index, 3:end]
+        end
+    end
+
+    return H₀_S₀
+
+end
+
+H₀_S₀ = readEnthalpyEntropyData("Hf0_Sf0.tsv", components_names)
+
+
+# Define reaction coefficients:
+# Reaction 1: 5 CO + 3 H20 -> 4 CO2 + H2 + CH4
+# Reaction 2: CO + 3 H2 -> CH4 + H20 
+
+ν = [4  0; -5 -1; 1 -3; -3  1; 1  1] # Should be done automatically in the future 
+ν*[0.0; 0.0] 
+
+function GibbsFreeEnergy_TP(T, P, H₀_S₀, ξ, ν, model, f_feed)
+    # Calculate the Gibbs free energy of the system
+    # T: Temperature in K
+    # P: Pressure in Pa
+    # ξ: Extent of reaction
+    # ν: Stoichiometric coefficients
+    # my_model: Model of the system
+    # Returns: Gibbs free energy in J
+
+    f_prod = f_feed + ν*ξ #Mass balance 
+    ΔHrxn = dot(H₀_S₀[:, 1], ν*ξ) #Heat released or taken from Reaction
+    ΔSrxn = dot(H₀_S₀[:, 2], ν*ξ) #Standard entropy of reaction
+    H = enthalpy(model, P, T, f_prod) + ΔHrxn #Outlet stream enthalpy
+    S = entropy(model, P, T, f_prod) + ΔSrxn #Outlet stream entropy
+    G = H - T*S   
+    return G
+
+end
+
+
+
+function Entropy_TP(T, P, H₀_S₀, ξ, ν, model, f_feed)
+    # Calculate the entropy of the system (J)      
+
+    f_prod = f_feed + ν*ξ 
+
+    ΔSrxn = dot(H₀_S₀[:, 2], ν*ξ) #Standard entropy of reaction times how much was converted from reactants to products
+
+    S = entropy(model, P, T, f_prod) + ΔSrxn 
+
+    return S
+
+end
+
+
+
+function EnergyBalance(T, P, H₀_S₀, ξ, ν, model, p)
+    # Calculate the energy balance of the system
+    # res: Energy balance error
+    # ξ: Extent of reaction
+    # p: Parameters of the system
+    # Returns: Energy balance error (res) in J
+
+    f_feed, T_feed = p[1], p[2] 
+
+    # Calculate the energy balance
+
+    f_prod = f_feed + ν*ξ #mass balance
+
+    Hin = enthalpy(model, P, T_feed, f_feed)
+    Hout = enthalpy(model, P, T, f_prod)
+    ΔHrxnX = dot(H₀_S₀[:, 1], ν*ξ)
+
+    res = Hin + ΔHrxnX - Hout 
+
+    return res
+
+end
+
+#Testing the functions:
+GibbsFreeEnergy_TP(300, 10^5, H₀_S₀, [0.0, 0.0], ν, my_model, [5.0, 20.0, 5.0, 5.0, 5.0])
+Entropy_TP(300, 10^5, H₀_S₀, [0.0, 0.0], ν, my_model, [5.0, 20.0, 5.0, 5.0, 5.0])
+EnergyBalance(500, 10^5, H₀_S₀, [0.0, 0.0], ν, my_model, ([5.0, 20.0, 5.0, 5.0, 5.0], 500))
+
+#Minimizing S(T,ξ) with respect to ξ constrained to energy balance (adiabatic reactor):
+using OptimizationOptimJL, Optimization, OptimizationMOI, Ipopt
+
+f_feed = [5.0, 20.0, 5.0, 5.0, 5.0] #["carbon dioxide", "carbon monoxide", "hydrogen", "water", "methane"]
+T_feed = 500.0 #K
+P = 101325.0 #Pa
+
+#Target function
+Smin(ξ_T) = -Entropy_TP(ξ_T[end], P, H₀_S₀, ξ_T[1:2], ν, my_model, f_feed)
+
+
+#Energy balance constraint
+E(ξ_T) = EnergyBalance(ξ_T[end], P, H₀_S₀, ξ_T[1:2], ν, my_model, [f_feed, T_feed])
+
+#product flow rate positive
+function f_prod_cons(ξ_T, f_feed, ν)
+
+    res = f_feed + ν*ξ_T[1:2]
+
+    return res
+
+end
+
+
+#Initial guess
+ξ₀ = [0.000; 0.000]
+T₀ = 600.0
+
+#Testing the functions:
+Smin([ξ₀ ; T₀])
+E([ξ₀ ; T₀])
+f_prod_cons([ξ₀ ; T₀], f_feed, ν)
+
+
+#Lower and upper bounds
+lb = [0.00; zeros(size(components_names, 1))]
+ub = [0.00; Inf*ones(size(components_names, 1))]
+
+optf = OptimizationFunction((x, p) -> Smin(x), Optimization.AutoForwardDiff(), cons = (res, x, p) -> [E(x), f_prod_cons(x, f_feed, ν)])
+
+prob = OptimizationProblem(optf, [ξ₀; T₀], lcons = lb, ucons = ub)
+
+sol = solve(prob, Ipopt.Optimizer())

src/Reactors/Gibbs_Reactor/GibbsReactor_symbolic.jl

@@ -0,0 +1,83 @@
+using Symbolics, SymbolicUtils, ModelingToolkit
+using Reexport
+@reexport using ModelingToolkit
+
+using Clapeyron
+using GCIdentifier
+using DelimitedFiles
+using LinearAlgebra
+
+components_names = ["carbon dioxide", "carbon monoxide", "hydrogen", "water", "methane"]
+abspath(joinpath(@__DIR__, "myShomate.csv"))
+syngas_prop_model =  ShomateIdeal(components_names; userlocations = [abspath(joinpath(@__DIR__, "myShomate.csv"))])
+
+n_reactions = 2
+
+function readEnthalpyEntropyData(file_path, component_names)
+    # Read formation enthalpy from file
+    # file_path: File path of the file containing the formation enthalpy
+    # component_names: Names of the components
+    # Returns: Formation enthalpy in J/mol
+
+    # initialize the formation enthalpy and entropy vector
+    H₀_S₀ = zeros(length(component_names), 2)
+
+    # Read the TSV file
+    data = readdlm(file_path, '\t', header = true)
+
+    # Get the column of chemical names
+    chemical_names = data[1][:, 2]
+
+    # Find the index of the chemical name in the column
+    for i in eachindex(component_names)
+        index = findfirst(x -> x == component_names[i], chemical_names)
+
+        # If the chemical name was found, return the corresponding line
+        if index !== nothing
+            H₀_S₀[i, :] =  data[1][index, 3:end]
+        end
+    end
+
+    return H₀_S₀
+
+end
+
+H₀_S₀ = readEnthalpyEntropyData(abspath(joinpath(@__DIR__, "Hf0_Sf0.tsv")), components_names)
+
+@register_symbolic Entropy_TP1(T, P, H₀_S₀, ξ, ν, model, f_feed)
+
+function Entropy_TP(T, P, H₀_S₀, ξ_1, ξ_2, ν, model, f_feed)
+    # Calculate the entropy of the system (J)      
+
+    f_prod = f_feed + ν*ξ 
+
+    ΔSrxn = dot(H₀_S₀[:, 2], ν*ξ) #Standard entropy of reaction times how much was converted from reactants to products
+
+    S = entropy(model, P, T, f_prod) + ΔSrxn 
+
+    return S
+
+end
+
+function GibbsReactor(; name, n_reactions, v)
+end
+
+ν = [4  0; -5 -1; 1 -3; -3  1; 1  1]
+
+vars = @variables begin
+    (T = 600.0), [description = "Reactor temperature"]
+    (ξ_1 = 0.5), [description = "Extent of rxn 1"]
+    (ξ_2 = 0.5), [description = "Extent of rxn 2"]
+    (s = 0.0), [description = "Entropy of the system"]
+end
+
+pars = @parameters begin
+    (P = 101325.0), [description = "Reactor pressure"]
+    (T_feed = 500.0), [description = "Feed temperature"]
+    (f_feed[1:length(components_names)] = [5.0, 20.0, 5.0, 5.0, 5.0]), [description = "Feed molar fraction"]
+end
+
+# TODO: structural_simplify optimization system.
+eqns = [
+    s ~ Entropy_TP(T, P, H₀_S₀, ξ, ν, syngas_prop_model, f_feed)
+]