Basic building blocks and CSTR

.gitignore

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

Project.toml

@@ -1,28 +1,19 @@
 name = "ProcessSimulator"
 uuid = "8886c03b-4dde-4be1-b6ee-87d056f985b8"
-authors = ["Avinash Subramanian"]
+authors = ["Chris Rackauckas", "Vinicius Viena Santana", "Sebastian Schmitt", "Andrés Riedemann", "Pierre Walker", "Avinash Subramanian"]
 version = "1.0.0-DEV"
 
 [deps]
-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]
-Clapeyron = "0.6"
-DifferentialEquations = "7"
-JSON = "0.21"
 ModelingToolkit = "9"
-NonlinearSolve = "3"
-OrdinaryDiffEq = "6"
 julia = "1.10"
 
 [extras]
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 
 [targets]
-test = ["SafeTestsets","Test"]
+test = ["Test", "SafeTestsets", "DifferentialEquations"]

src/ProcessSimulator.jl

@@ -1,26 +1,19 @@
 module ProcessSimulator
 
-
-using ModelingToolkit, JSON, OrdinaryDiffEq
+using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
-import ModelingToolkit: scalarize, equations, get_unknowns, defaults
-using Clapeyron
-
-include("utils.jl")
-
-include("Sources/MaterialSource.jl")
-
-include("Reactors/ReactionManager/KineticReaction.jl")
-
-include("Reactors/CSTR.jl")
-
-include("Valve/Valves.jl")
-
-include("HeatExchange/Jacket.jl")
+using ModelingToolkit: scalarize, equations, get_unknowns
 
-include("Separation/Flash.jl")
+# Base
+include("base/materials.jl")
+include("base/base_components.jl")
+include("base/utils.jl")
 
-include("Sources/Sourceutils.jl")
+# Fluid handling
+include("fluid_handling/compressors.jl")
+include("fluid_handling/heat_exchangers.jl")
 
+# Reactors
+include("reactors/CSTR.jl")
 
 end

src/base/base_components.jl

@@ -0,0 +1,143 @@
+
+@component function MaterialStream(ms::MaterialSource;phase="unknown",name)
+    @named c1 = MaterialConnector(ms)
+    @named c2 = MaterialConnector(ms)
+    mcons = [c1,c2]
+
+    vars = @variables begin
+        T(t),               [description="temperature"]         #, unit=u"K"]
+        ϱ(t),               [description="density"]             #, unit=u"mol m^-3"]
+        p(t),               [description="pressure"]            #, unit=u"Pa"]
+        h(t),               [description="molar enthalpy"]      #, unit=u"J mol^-1"]
+        n(t),               [description="molar flow "]         #, unit=u"mol s^-1"]
+        xᵢ(t)[1:ms.N_c],    [description="mole fractions"]      #, unit=u"mol mol^-1"]
+    end
+
+    eqs = [
+        # EOS 
+        ϱ ~ ms.molar_density(p,T,xᵢ;phase=phase),
+        h ~ ms.VT_enthalpy(ϱ,T,xᵢ),
+        1.0 ~ sum(collect(xᵢ)),
+        # Connector "1"
+        T ~ c1.T,
+        ϱ ~ c1.ϱ,
+        p ~ c1.p,
+        h ~ c1.h,
+        n ~ c1.n,
+        scalarize(xᵢ .~ c1.xᵢ)...,
+        # Connector "2"
+        T ~ c2.T,
+        ϱ ~ c2.ϱ,
+        p ~ c2.p,
+        h ~ c2.h,
+        n ~ -c2.n,
+        scalarize(xᵢ .~ c2.xᵢ)...,   
+    ]
+
+    return ODESystem(eqs, t, collect(Iterators.flatten(vars)), []; name, systems=mcons)
+end
+
+@connector function MaterialConnector(ms::MaterialSource;name)
+    vars = @variables begin
+        T(t),               [description="temperature", output=true]        #, unit=u"K"]
+        ϱ(t),               [description="density", output=true]            #, unit=u"mol m^-3"]
+        p(t),               [description="pressure"]                        #, unit=u"Pa"]
+        h(t),               [description="molar enthalpy", output=true]     #, unit=u"J mol^-1"]
+        xᵢ(t)[1:ms.N_c],    [description="mole fractions", output=true]     #, unit=u"mol mol^-1"]
+        n(t),               [description="total molar flow", connect=Flow]  #, unit=u"mol s^-1"]
+    end
+
+    return ODESystem(Equation[], t, collect(Iterators.flatten(vars)), []; name)
+end
+
+@component function HeatConnector(;name)
+    vars = @variables begin
+        Q(t),              [description="heat flux"]  #, unit=u"J s^-1"]
+    end
+
+    return ODESystem(Equation[], t, vars, []; name)
+end
+
+@component function WorkConnector(;name)
+    vars = @variables begin
+        W(t),              [description="power"]      #, unit=u"J s^1"]
+    end
+
+    return ODESystem(Equation[], t, vars, []; name)
+end
+
+@component function SimpleControlVolume(ms::MaterialSource;
+                                        N_mcons,
+                                        N_heats=0,
+                                        N_works=0,
+                                        name)
+    # Init 
+    mcons = [MaterialConnector(ms;name=Symbol("c$i")) for i in 1:N_mcons]
+    works = [WorkConnector(name=Symbol("w$i")) for i in 1:N_works]
+    heats = [HeatConnector(name=Symbol("q$i")) for i in 1:N_heats]
+
+    vars = @variables begin
+        ΔH(t), [description="Enthalpy difference inlets/outlets"]  #, unit=u"J/s"]
+        ΔE(t), [description="Added/removed heat or work"]          #, unit=u"J/s"]
+    end
+
+    eqs = Equation[
+        ΔH ~ sum([c.h*c.n for c in mcons])
+        ΔE ~ (isempty(heats) ? 0.0 : sum([q.Q for q in heats])) + (isempty(works) ? 0.0 : sum([w.W for w in works]))
+        # Energy balance
+        0.0 ~ ΔH + ΔE
+        # Mole balance
+        [0.0 ~ sum([c.xᵢ[j]*c.n for c in mcons]) for j in 1:ms.N_c][:]
+    ]
+
+    return ODESystem(eqs, t, vars, []; name, systems=[mcons...,works...,heats...])
+end
+
+@component function TPControlVolume(ms::MaterialSource;
+                                    N_mcons,
+                                    N_heats=0,
+                                    N_works=0,
+                                    phases=["unknown"],
+                                    reactive=false,
+                                    name)
+    # Init 
+    N_ph = length(phases)
+    mcons = [MaterialConnector(ms;name=Symbol("c$i")) for i in 1:N_mcons]
+    works = [WorkConnector(name=Symbol("w$i")) for i in 1:N_works]
+    heats = [HeatConnector(name=Symbol("q$i")) for i in 1:N_heats]
+
+    vars = @variables begin
+        T(t),                       [description="temperature", bounds=(0,Inf)]         #, unit=u"K"]
+        p(t),                       [description="pressure", bounds=(0,Inf)]            #, unit=u"Pa"]
+        ϱ(t)[1:N_ph],               [description="density", bounds=(0,Inf)]             #, unit=u"mol m^-3"]
+        (nᵢ(t))[1:N_ph,1:ms.N_c],   [description="molar holdup", bounds=(0,Inf)]        #, unit=u"mol"]
+        (xᵢ(t))[1:N_ph,1:ms.N_c],   [description="mole fractions", bounds=(0,1)]        #, unit=u"mol mol^-1"]
+        n(t),                       [description="total molar holdup", bounds=(0,Inf)]  #, unit=u"mol"]
+        U(t),                       [description="internal energy"]                     #, unit=u"J"]
+        ΔH(t),                      [description="enthalpy difference inlets/outlets"]  #, unit=u"J/s"]
+        ΔE(t),                      [description="added/removed heat or work"]          #, unit=u"J/s"]
+        V(t),                       [description="volume", bounds=(0,Inf)]             #, unit=u"m^3"]
+    end
+    reactive ? append!(vars, @variables begin
+        ΔnR(t)[1:ms.N_c],           [description="molar holdup change by reaction"]     #, unit=u"mol"])
+        ΔHᵣ(t),                     [description="enthalpy of reaction"]                #, unit=u"J/s"]
+    end) : nothing
+
+    eqs = [
+        ΔH ~ sum([c.h*c.n for c in mcons]),
+        ΔE ~ (isempty(heats) ? 0.0 : sum([q.Q for q in heats])) + (isempty(works) ? 0.0 : sum([w.W for w in works])),
+        # Energy balance
+        D(U) ~ ΔH + ΔE + (reactive ? ΔHᵣ : 0.0),
+        # Mole balance
+        [D(sum(collect(nᵢ[:,i]))) ~ sum([c.xᵢ[i]*c.n for c in mcons]) + 
+            (reactive ? ΔnR[i] : 0.0) for i in 1:ms.N_c]...,
+        n ~ sum(collect([nᵢ...])),
+        [xᵢ[j,i] ~ nᵢ[j,i]/sum(collect(nᵢ[j,:])) for i in 1:ms.N_c, j in 1:N_ph]...,
+        # Thermodynamic system properties
+        [ϱ[j] ~ ms.molar_density(p,T,collect(xᵢ[j,:]);phase=phases[j]) for j in 1:N_ph]...,
+        U ~ sum([ms.VT_internal_energy(ϱ[j],T,collect(xᵢ[j,:]))*sum(collect(nᵢ[j,:])) for j in 1:N_ph]),
+        V ~ n / sum(collect(ϱ)),
+    ]
+
+    return ODESystem(eqs, t, collect(Iterators.flatten(vars)), []; name, systems=[mcons...,works...,heats...])
+end

src/base/materials.jl

@@ -0,0 +1,51 @@
+abstract type AbstractMaterialSource end
+
+struct Reaction
+    ν::Vector{Float64}              # Stoichiometry
+    r::Function                     # Reaction rate (defined as f(p,T,xᵢ)) in mol/s
+    Δhᵣ::Function                   # Reaction enthalpy (defined as f(T)) in J/mol
+    Reaction(;ν, r, Δhᵣ) = new(ν, r, Δhᵣ)
+end
+
+struct MaterialSource <: AbstractMaterialSource
+    name::String                    # Name of the material source
+    components::Vector{String}      # Component names
+    N_c::Int                        # Number of components
+    Mw::Vector{Float64}             # Molar weight in kg/mol
+    pressure::Function              # Pressure function (defined as f(ϱ,T,xᵢ;kwargs...)) in Pa
+    molar_density::Function         # Molar density function (defined as f(p,T,xᵢ;kwargs...)) in mol/m³
+    VT_internal_energy::Function    # Internal energy function (defined as f(ϱ,T,xᵢ;kwargs...)) in J/mol
+    VT_enthalpy::Function           # Enthalpy function (defined as f(ϱ,T,xᵢ;kwargs...)) in J/mol
+    VT_entropy::Function            # Entropy function (defined as f(ϱ,T,xᵢ;kwargs...)) in J/(mol K)
+    tp_flash::Function              # Flash function (defined as f(p,T,xᵢ;kwargs...))
+    reaction::Vector{Reaction}      # Reaction struct 
+end
+
+function MaterialSource(components::Union{String,Vector{String}}; kwargs...)
+    components = components isa String ? [components] : components
+
+    # Check for mandatory keyword arguments
+    mandatory = [:Mw, :molar_density, :VT_enthalpy]
+    [haskey(kwargs, k) || throw(ArgumentError("Missing keyword argument $k")) for k in mandatory]
+
+    N_c = length(components)
+    length(kwargs[:Mw]) == N_c || throw(ArgumentError("Length of Mw must be equal to the number of components"))
+    name = haskey(kwargs, :name) ? kwargs[:name] : join(components, "_")
+
+    f_NA(field) = error("Function $field not defined in MaterialSource")
+
+    MaterialSource(
+        name,
+        components,
+        N_c,
+        kwargs[:Mw] isa Number ? [kwargs[:Mw]] : kwargs[:Mw],
+        get(kwargs, :pressure, (a,T,n;kws...) -> f_NA(:pressure)),
+        kwargs[:molar_density],
+        get(kwargs, :VT_internal_energy, (a,T,n;kws...) -> f_NA(:VT_internal_energy)),
+        kwargs[:VT_enthalpy],
+        get(kwargs, :VT_entropy, (a,T,n;kws...) -> f_NA(:VT_entropy)),
+        get(kwargs, :tp_flash, (a,T,n;kws...) -> f_NA(:tp_flash)),
+        get(kwargs, :reactions, Reaction[])
+    )
+end
+

src/base/utils.jl

@@ -0,0 +1,28 @@
+@component function Port(ms;phase="unknown", name)
+    @named c = MaterialConnector(ms)
+
+    vars = @variables begin
+        T(t),               [description="inlet temperature"]   #, unit=u"K"]
+        ϱ(t),               [description="inlet density"]       #, unit=u"mol m^-3"]
+        p(t),               [description="inlet pressure"]      #, unit=u"Pa"]
+        n(t),               [description="inlet molar flow"]    #, unit=u"mol s^-1"]
+        m(t),               [description="inlet mass flow"]     #, unit=u"kg s^-1"]
+        xᵢ(t)[1:ms.N_c],    [description="inlet mole fractions"]#, unit=u"mol mol^-1"] 
+    end
+
+    eqs = Equation[
+        # EOS 
+        ϱ ~ ms.molar_density(p,T,xᵢ;phase=phase),
+        1.0 ~ sum(collect(xᵢ)),
+        c.h ~ ms.VT_enthalpy(ϱ,T,xᵢ),
+        # Connector
+        T ~ c.T,
+        ϱ ~ c.ϱ,
+        p ~ c.p,
+        scalarize(xᵢ .~ c.xᵢ)...,
+        n ~ m / sum(ms.Mw[i] * xᵢ[i] for i in 1:ms.N_c),
+        0.0 ~ n + c.n,
+     ]
+
+    return ODESystem(eqs, t, collect(Iterators.flatten(vars)), []; name, systems=[c])
+end

src/fluid_handling/compressors.jl

@@ -0,0 +1,26 @@
+@component function SimpleAdiabaticCompressor(ms::MaterialSource; name)
+
+    # Subsystems
+    @named cv = SimpleControlVolume(ms;N_mcons=2,N_works=1)
+
+    # Variables
+    vars = @variables begin
+        W(t),               [description="power"]       #, unit=u"J s^1"]
+        T2_s(t),            [description="temperature"] #, unit=u"K"]
+        ϱ2_s(t),            [description="density"]     #, unit=u"mol m^-3"]
+    end
+
+    pars = @parameters begin
+        ηᴱ = 1.0,           [description="efficiency"]
+    end
+
+    # Equations
+    eqs = [
+        cv.w1.W ~ W,
+        ms.VT_entropy(cv.c1.ϱ,cv.c1.T,cv.c1.xᵢ) ~ ms.VT_entropy(cv.c2.ϱ,cv.c2.T,cv.c2.xᵢ),
+        ϱ2_s ~ ms.molar_density(cv.c2.p,T2_s,cv.c2.xᵢ),
+        ηᴱ ~ cv.w1.W / ((ms.VT_enthalpy(ϱ2_s,T2_s,cv.c2.xᵢ) - cv.c1.h)*cv.c1.n),
+    ]
+
+    return ODESystem(eqs, t, vars, pars; name, systems=[cv])
+end

src/fluid_handling/heat_exchangers.jl

@@ -0,0 +1,18 @@
+@component function SimpleIsobaricHeatExchanger(ms::MaterialSource; name)
+
+    # Subsystems
+    @named cv = SimpleControlVolume(ms;N_mcons=2,N_heats=1)
+
+    # Variables
+    vars = @variables begin
+        Q(t),               [description="heat flux"] #, unit=u"J s^-1"]
+    end 
+
+    # Equations
+    eqs = [
+        cv.q1.Q ~ Q,
+        cv.c1.p ~ cv.c2.p,
+    ]
+
+    return ODESystem(eqs, t, vars, []; systems=[cv], name)
+end

src/reactors/CSTR.jl

@@ -0,0 +1,34 @@
+@component function CSTR(ms::MaterialSource;name, i_reacts=1:length(ms.reaction), flowtype="")
+    # Subsystems
+    @named cv = TPControlVolume(ms;N_mcons=2,N_heats=1,N_works=1,phases=["liquid"],reactive=true)
+
+    # Variables
+    vars = @variables begin
+        Q(t),                   [description="heat flux"] #, unit=u"J s^-1"]
+    end
+
+    # Parameters
+    pars = @parameters begin
+        W=0.0,                  [description="work"]      #, unit=u"J s^1"]
+    end
+
+    # Equations
+    eqs = Equation[
+        cv.w1.W ~ W,
+        cv.q1.Q ~ Q,
+        [cv.c2.xᵢ[i] ~ cv.xᵢ[1,i] for i in 1:ms.N_c-1]...,
+        cv.c2.p ~ cv.p,
+        cv.c2.T ~ cv.T,
+        # Change of moles by reaction
+        [cv.ΔnR[i] ~ reac.r(cv.p,cv.T,collect(cv.xᵢ[1,:]))*reac.ν[i]*cv.n for i in 1:ms.N_c, reac in ms.reaction[i_reacts]]...,
+        # Enthalpy of reaction
+        [cv.ΔHᵣ ~ reac.r(cv.p,cv.T,collect(cv.xᵢ[1,:]))*reac.Δhᵣ(cv.T)*cv.n for reac in ms.reaction[i_reacts]]...,
+    ]
+    if flowtype == "const. mass"
+        push!(eqs,0.0 ~ cv.c1.n * sum(ms.Mw[i]*cv.c1.xᵢ[i] for i in 1:ms.N_c) + cv.c2.n * sum(ms.Mw[i]*cv.c2.xᵢ[i] for i in 1:ms.N_c))
+    elseif flowtype == "const. volume"
+        push!(eqs,D(cv.V) ~ 0.0)
+    end
+
+    return ODESystem(eqs, t, vars, pars; name, systems=[cv])
+end