Vini/refactorconnectors

.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,12 @@ authors = ["Avinash Subramanian"]
 version = "1.0.0-DEV"
 
 [deps]
-SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+Clapeyron = "7c7805af-46cc-48c9-995b-ed0ed2dc909a"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 
 [compat]
 julia = "1.10"

src/HeatExchange/Jacket.jl

@@ -0,0 +1,99 @@
+@component function Jacket(;substances_user,
+    Nc = length(substances_user),
+    phase,
+    thermal_fluid_model,
+    heat_transfer_coef,
+    name,
+   )
+
+
+
+properties = Dict(subs => load_component_properties(subs) for subs in substances_user)
+MWs = [properties[subs]["MW"] for subs in substances_user]
+
+pars = @parameters begin
+    U = heat_transfer_coef
+end
+
+systems = @named begin 
+    Out = matcon(; Nc = Nc)
+    In = matcon(; Nc = Nc)
+    EnergyCon = thermal_energy_connector()
+end 
+
+vars = @variables begin
+    LMTD(t), [description = "log mean temperature difference (-)"]
+    Tⱼ(t), [description = "Thermal fluid exiting temperature (K)"]
+    H(t), [description = "Thermal fluid exiting enthalpy (J/mol)"]
+    S(t), [description = "Thermal fluid exiting entropy (J/K/mol)"]
+    ρ(t), [description = "Thermal fluid exiting molar density (mol/m³)"]
+    ρʷ(t), [description = "Thermal fluid exiting mass density (kg/m³)"]
+    P_out(t), [description = "Thermal fluid exiting pressure (Pa)"]
+    F_out(t), [description = "Thermal fluid exiting molar flow rate (mol/s)"] 
+    Fʷ_out(t), [description = "Thermal fluid exiting mass flow rate (mol/s)"]
+    Q_out(t), [description = "Thermal fluid exiting volumetric flow rate (m³/s)"] 
+    Q̇(t), [description = "Heat transfer rate from/to the jacket"] 
+end
+
+
+enthalpy_eq = [H ~ enthalpy(thermal_fluid_model, P_out, Tⱼ, Out.z₁, phase = "unknown")]
+
+entropy_eq = [S ~ 0.0]
+
+densities_eq = [ρ ~ molar_density(thermal_fluid_model, P_out, Tⱼ, Out.z₁, phase = "unknown")
+             ρʷ ~ mass_density(thermal_fluid_model, P_out, Tⱼ, Out.z₁, phase = "unknown")]
+
+pressure_drop = [P_out ~ In.P]
+
+mass_balance = [In.Fʷ - Fʷ_out ~ 0.0
+                Q_out ~ Fʷ_out/ρʷ
+                Q_out ~ F_out/ρ]
+
+lmtd = [LMTD ~ log((EnergyCon.T - In.T)/(EnergyCon.T - Tⱼ))]
+
+heat_flux = [Q̇ ~ U*EnergyCon.A*(In.T - Tⱼ)/LMTD]
+
+energy_balance = [In.H*In.F - H*F_out - Q̇ ~ 0.0] # Quasi steady state assumption (Temperature change is much faster than reactor change)
+
+
+
+#Outlet connector equations:
+out_conn = [Out.P ~ P_out
+            Out.T ~ Tⱼ
+            Out.F ~ F_out
+            Out.Fʷ ~ Fʷ_out
+            Out.H ~ H
+            Out.S ~ S
+            Out.ρʷ ~ ρʷ
+            Out.ρ ~ ρ
+            scalarize(Out.z₁ .~ In.z₁)...
+            Out.MW[1] ~ MWs
+            EnergyCon.ϕᴱ ~ Q̇
+]
+
+
+if phase == :liquid
+    out_conn_phases = [
+                    scalarize(Out.z₂ .~ 0.0)...
+                    scalarize(Out.z₃ .~ In.z₁)...
+                    Out.MW[2] ~  0.0
+                    Out.MW[3] ~ MWs
+                    Out.α_g ~ 0.0]
+
+    elseif phase == :vapor
+        out_conn_phases = [
+        scalarize(Out.z₂ .~ In.z₁)...
+        scalarize(Out.z₃ .~ 0.0)...
+        Out.MW[2] ~ MWs
+        Out.MW[3] ~ 0.0
+        Out.α_g ~ 1.0]
+end
+
+eqs = [enthalpy_eq...; entropy_eq...; densities_eq...; pressure_drop...; mass_balance...; lmtd...; heat_flux...; energy_balance...; out_conn...; out_conn_phases...]
+
+ODESystem([eqs...;], t, collect(Iterators.flatten(vars)), collect(Iterators.flatten(pars)); name, systems = [Out, In, EnergyCon])
+
+end
+
+
+export Jacket

src/ProcessSimulator.jl

@@ -1,7 +1,22 @@
 module ProcessSimulator
 
-# Write your package code here.
-export Gibbs
-include("Reactors/Gibbs.jl")
+
+using ModelingToolkit, JSON, OrdinaryDiffEq
+using ModelingToolkit: t_nounits as t, D_nounits as D
+import ModelingToolkit: scalarize, equations, get_unknowns, defaults
+using Clapeyron
+
+include("utils")
+
+include("Sources/MaterialSource.jl")
+
+include("Reactors/ReactionManager/KineticReaction.jl")
+
+include("Reactors/CSTR.jl")
+
+include("HeatExchange/Jacket.jl")
+
+include("Sources/Sourceutils.jl")
+
 
 end

src/Reactors/CSTR.jl

@@ -0,0 +1,155 @@
+@component function CSTR(; substances_user, 
+    Nc = length(substances_user), 
+    phase, 
+    model,
+    Reaction,
+    ninports,  
+    Ac, 
+    height_out_port,
+    name,
+    guesses 
+   )
+
+
+
+#Constants
+gramsToKilograms = 10^(-3)
+Rᵍ = 8.314 # J/(mol K)
+Nri = Reaction.Nri
+
+#Properties of individual substances
+properties = Dict(subs => load_component_properties(subs) for subs in substances_user)
+MWs = [properties[subs]["MW"] for subs in substances_user]
+ΔH₀f = [properties[subs]["IGHF"]/10^3 for subs in substances_user] # (IG formation enthalpy) J/mol
+
+
+pars = @parameters begin 
+Af_r[1:Nri] = Reaction.Af_r, [description = "Arrhenius constant of each reaction at given temperature ()"]
+Coef_Cr[1:Nri, 1:Nc] = Reaction.Coef_Cr, [description = "Stoichiometric coefficients of each component in each reaction (-)"]
+Do_r[1:Nri, 1:Nc] = Reaction.Do_r, [description = "Forward order of the components (-) "]
+Ef_r[1:Nri] = Reaction.Ef_r, [description = "Activation energy of each reaction at given temperature ()"]
+A = Ac, [description = "Cross sectional area of the tank (m²)"]
+end
+
+#Ports creation
+
+InPorts = [matcon(; Nc = Nc, name = Symbol("InPorts$i")) for i in 1:ninports]
+
+OutPorts = @named begin
+    Out = matcon(; Nc = Nc) 
+    EnergyCon = thermal_energy_connector()
+end   
+
+vars = @variables begin
+    M(t), [description = "Mass holdup in the tank (kg)"]
+    N(t), [description = "Total molar holdup in the tank (kmol)"]
+    V(t), [description = "Volume holdup in the tank (m³)", guess = guesses[:V]]
+    (Nᵢ(t))[1:Nc], [description = "Molar holdup of each component in the tank (mol)"]
+    (Cᵢ(t))[1:Nc], [description = "Concentration of each component in the tank (mol/m³)", guess = guesses[:Cᵢ]]
+    ρ(t), [description = "Molar Density of the fluid in the tank (mol/m³)"]
+    ρʷ(t), [description = "Mass Density of the fluid in the tank (kg/m³)"]
+    MW(t), [description = "Molecular weight of fluid in the tank (kg/kmol)"]  
+    T(t), [description = "Temperature of vessel contents (K)"]  
+    P_out(t), [description = "Pressure at the outlet stream level (Pa)"]
+    H(t), [description = "Enthalpy holdup of the fluid in the tank (J)"] 
+    #S(t), [description = "Entropy holdup of the fluid in the tank (J/K)"]
+    F_out(t), [description = "Outlet molar flow rate (mol/s)", guess = guesses[:F_out]] 
+    Fʷ_out(t), [description = "Outlet molar flow rate (mol/s)", guess = guesses[:Fʷ_out]]
+    Q_out(t), [description = "Outlet volumetric flow rate (m³/s)"] # DoF
+    height(t), [description = "Liquid level in vessel measured from bottom of the tank (m)"]
+    P_atm(t), [description = "Tank pressure (Pa)"] # Equal to inlet pressures.
+    Q̇(t), [description = "Heat transfer rate (J/s)"] # Potential DoF
+    (r(t))[1:Nri], [description = "Rate of each reaction for each component (mol/s/m³)"]
+    (R(t))[1:Nc], [description = "Overall reaction rate (mol/s/m³)"]
+
+    (Cᵢ_in(t))[1:Nc, 1:ninports], [description = "Inlet concentration of each component (mol/m³)"] # DoF through inlet stream
+    (F_in(t))[1:ninports], [description = "Inlet molar flow rate (mol/s)"] # DoF through inlet stream
+    (Q_in(t))[1:ninports], [description = "Inlet volumetric flow rate(s) (m³/s)"]
+    (T_in(t))[1:ninports], [description = "Inlet temperature (K)"] # DoF through inlet stream
+    (h_in(t))[1:ninports], [description = "Inlet specific enthalpy (J/mol)"]
+    (ρ_in(t))[1:ninports], [description = "Inlet density (mol/m³)"]
+    (ρʷ_in(t))[1:ninports], [description = "Inlet density (mol/m³)"]
+
+end
+
+
+
+#Reaction equations
+reaction_rate = [r[i] ~ Af_r[i]*exp(-Ef_r[i]/(Rᵍ*T))*prod(scalarize((Cᵢ[:].^Do_r[i, :]))) for i in 1:Nri] # If there's an inert, the order is just zero, but it has to be written
+overall_reaction_rate = [R[i] ~ sum(scalarize(r[:].*Coef_Cr[:, i])) for i in 1:Nc]  # If there's an inert, the coefficient is just zero, but it has to be written
+
+
+#Inlet connector variables's equations
+atm_pressure = [P_atm ~ InPorts[1].P]
+mass_density_eqs = [ρʷ_in[j] ~ InPorts[j].ρʷ for j in 1:ninports]
+molar_density_eqs = [ρ_in[j] ~ InPorts[j].ρ for j in 1:ninports]
+inletenthalpy = [h_in[j] ~ InPorts[j].H for j in 1:ninports]
+inletconcentrations = [Cᵢ_in[i, j] ~ InPorts[j].z₁[i]*ρ_in[j] for j in 1:ninports for i in 1:Nc]
+inlettemperature_eqs = [T_in[j] ~ InPorts[j].T for j in 1:ninports]
+inletmolarflow_eqs = [F_in[j] ~ InPorts[j].F for j in 1:ninports]
+volumetricflow_eqs = [Q_in[j] ~ F_in[j] / ρ_in[j] for j in 1:ninports]
+
+
+#Outlet connector equations:
+out_conn = [Out.P ~ P_out
+            Out.T ~ T
+            Out.F ~ F_out
+            Out.Fʷ ~ Fʷ_out
+            Out.H ~ H/N
+            Out.S ~ 0.0
+            Out.ρʷ ~ ρʷ
+            Out.ρ ~ ρ
+            scalarize(Out.z₁ .~ Nᵢ/N)...
+            Out.MW[1] ~ MW
+]
+
+if phase == :liquid
+out_conn_phases = [
+                scalarize(Out.z₂ .~ 0.0)...
+                scalarize(Out.z₃ .~ Nᵢ/N)...
+                Out.MW[2] ~  0.0
+                Out.MW[3] ~ MW
+                Out.α_g ~ 0.0]
+
+elseif phase == :vapor
+    out_conn_phases = [
+    scalarize(Out.z₂ .~ Nᵢ/N)...
+    scalarize(Out.z₃ .~ 0.0)...
+    Out.MW[2] ~ MW
+    Out.MW[3] ~ 0.0
+    Out.α_g ~ 1.0]
+end
+
+
+#balances
+
+component_balance = [D(Nᵢ[i]) ~ sum(scalarize(Q_in[:].*Cᵢ_in[i, :])) - Q_out*Cᵢ[i] + R[i]*V for i in 1:Nc] #Neglectable loss to vapor phase head space
+energy_balance = [D(H) ~ sum(scalarize(Q_in.*ρ_in.*h_in)) - Q_out*ρ*H/N + Q̇]
+jacket_energy_balance = [Q̇ ~ -2.27*4184.0*(T - 288.7)*(1.0 - exp(-8440.26/(2.27*4184)))] 
+mass_volume_eq = [ρʷ*V ~ M, ρ*V ~ N] 
+mol_holdup = [N ~ sum(collect(Nᵢ))]
+mol_to_concentration = [scalarize(Nᵢ .~ Cᵢ*V)...]
+height_to_volume = [height*A ~ V]
+volumetricflow_to_molarflow = [Q_out ~ F_out/ρ, Q_out ~ sum(scalarize(Q_in))]
+volumetricflow_to_massflow = [Q_out ~ Fʷ_out/ρʷ]
+
+#Thermodynamic properties (outlet)
+pressure_out = [phase == :liquid ? P_out ~ P_atm : P_out ~ P_atm] #Estimation considering static pressure (May be off as tank is agitated and not static)
+density_eqs = [ρ ~ molar_density(model, P_out, T, Nᵢ, phase = "unknown")] 
+mass_density = [ρʷ ~ ρ*MW]
+globalEnthalpy_eq = [H ~ enthalpy(model, P_out, T, Nᵢ, phase = "unknown") + sum(scalarize(ΔH₀f.*Nᵢ))]
+molar_mass = [MW ~ sum(scalarize(MWs.*Nᵢ))/N*gramsToKilograms]
+
+
+eqs = [reaction_rate...; overall_reaction_rate...; atm_pressure...; mass_density_eqs...; molar_density_eqs...; inletenthalpy...; inletconcentrations...; inlettemperature_eqs...; inletmolarflow_eqs...; volumetricflow_eqs...; out_conn...;
+out_conn_phases...; component_balance...; energy_balance...; jacket_energy_balance...; mass_volume_eq...; mol_holdup...; mol_to_concentration...; height_to_volume...; volumetricflow_to_molarflow...; volumetricflow_to_massflow...;
+pressure_out...; density_eqs...; mass_density...; globalEnthalpy_eq...; molar_mass...]
+
+
+ODESystem([eqs...;], t, collect(Iterators.flatten(vars)),
+ collect(Iterators.flatten(pars)); name, systems = [InPorts...; OutPorts])
+
+end
+
+export CSTR
+

src/Reactors/ReactionManager/KineticReaction.jl

@@ -0,0 +1,24 @@
+# Definindo a estrutura
+struct KineticReactionNetwork{T}
+    Coef_Cr::Union{Vector{T}, Matrix{T}} # Stoichiometric coefficients of each component in each reaction (-)
+    Do_r::Union{Vector{T}, Matrix{T}} # Forward order of the components ()
+    substances_user::Vector{String} # Substances in the reaction network
+    Nc::Int # Number of components in the reaction network
+    Nri::Int # Number of reactions in the reaction network (r)
+    Af_r::T # Arrhenius constant of each reaction (s⁻¹)
+    Ef_r::T # Activation energy of each reaction (J.mol⁻¹)
+    name::String
+end
+
+# Construtor para aceitar palavras-chave
+function KineticReactionNetwork(; Coef_Cr::Union{Vector{T}, Matrix{T}}, Do_r::Union{Vector{T}, Matrix{T}},
+        substances_user::Vector{String}, Af_r::T, Ef_r::T, name::String) where {T}
+
+    Nc = length(substances_user)
+    Nri = Coef_Cr isa Matrix ? size(Coef_Cr, 1) : length(Coef_Cr)
+
+
+    return KineticReactionNetwork{T}(Coef_Cr, Do_r, substances_user, Nc, Nri, Af_r, Ef_r, name)
+end
+
+export KineticReactionNetwork