From 44daee861d46c394e3d62e53e07740581d8518ae Mon Sep 17 00:00:00 2001 From: se-schmitt Date: Wed, 9 Oct 2024 12:01:26 +0200 Subject: [PATCH 1/7] Implementation of base components and validation by simple thermodynamic cycle --- .gitignore | 3 +- Project.toml | 6 +- src/ProcessSimulator.jl | 15 ++- src/Reactors/Gibbs.jl | 3 - src/base/base_components.jl | 184 ++++++++++++++++++++++++++ src/base/materials.jl | 13 ++ src/base/utils.jl | 10 ++ src/fluid_handling/compressors.jl | 27 ++++ src/fluid_handling/heat_exchangers.jl | 17 +++ src/reactors/CSTR.jl | 8 ++ src/separation/distillation.jl | 19 +++ test/Reactor_tests/gibbs_tests.jl | 6 - test/base/simple_steady_state.jl | 93 +++++++++++++ test/reactors/simple_cstr.jl | 22 +++ test/runtests.jl | 4 +- 15 files changed, 413 insertions(+), 17 deletions(-) delete mode 100644 src/Reactors/Gibbs.jl create mode 100644 src/base/base_components.jl create mode 100644 src/base/materials.jl create mode 100644 src/base/utils.jl create mode 100644 src/fluid_handling/compressors.jl create mode 100644 src/fluid_handling/heat_exchangers.jl create mode 100644 src/reactors/CSTR.jl create mode 100644 src/separation/distillation.jl delete mode 100644 test/Reactor_tests/gibbs_tests.jl create mode 100644 test/base/simple_steady_state.jl create mode 100644 test/reactors/simple_cstr.jl diff --git a/.gitignore b/.gitignore index 8d27bf6..b415cf5 100644 --- a/.gitignore +++ b/.gitignore @@ -3,4 +3,5 @@ *.jl.mem /docs/Manifest.toml /docs/build/ -Manifest.toml \ No newline at end of file +Manifest.toml +.vscode/ \ No newline at end of file diff --git a/Project.toml b/Project.toml index f0d28c1..c1e8275 100644 --- a/Project.toml +++ b/Project.toml @@ -4,13 +4,15 @@ authors = ["Avinash Subramanian"] version = "1.0.0-DEV" [deps] -SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" [compat] julia = "1.10" [extras] +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" [targets] -test = ["Test"] +test = ["Test", "SafeTestsets", "DifferentialEquations"] diff --git a/src/ProcessSimulator.jl b/src/ProcessSimulator.jl index 29592a9..16658d8 100644 --- a/src/ProcessSimulator.jl +++ b/src/ProcessSimulator.jl @@ -1,7 +1,16 @@ module ProcessSimulator -# Write your package code here. -export Gibbs -include("Reactors/Gibbs.jl") +using ModelingToolkit +using ModelingToolkit: t_nounits as t, D_nounits as D +using ModelingToolkit: scalarize, equations, get_unknowns + +# Base +include("base/materials.jl") +include("base/base_components.jl") +include("base/utils.jl") + +# Fluid handling +include("fluid_handling/compressors.jl") +include("fluid_handling/heat_exchangers.jl") end diff --git a/src/Reactors/Gibbs.jl b/src/Reactors/Gibbs.jl deleted file mode 100644 index 3e90b2e..0000000 --- a/src/Reactors/Gibbs.jl +++ /dev/null @@ -1,3 +0,0 @@ -function Gibbs(x) - return 2x -end \ No newline at end of file diff --git a/src/base/base_components.jl b/src/base/base_components.jl new file mode 100644 index 0000000..66973b1 --- /dev/null +++ b/src/base/base_components.jl @@ -0,0 +1,184 @@ +# Base components + +@connector function State(ms::MaterialSource;N_c=1,n_bnds=(-Inf,Inf),phase="unknown",name) + vars = @variables begin + T(t), [description="temperature", unit="K", output=true] + p(t), [description="pressure", unit="Pa", output=true] + ϱ(t), [description="density", unit="mol/m³", output=true] + (nᵢ(t))[1:N_c], [description="molar flow", unit="mol/s", bounds=n_bnds, output=true] + n(t), [description="total molar flow", unit="mol/s", output=true] + end + + eqs = [ + # Density + ϱ ~ ms.molar_density(p,T,nᵢ;phase=phase) + n ~ sum(nᵢ) + ] + + pars = @parameters begin + end + + ODESystem(eqs, t, collect(Iterators.flatten(vars)), pars; name) +end + +""" + MaterialStream(ms::MaterialSource;N_c=1, phase="unknown", flowdir=:generic, name) + +Create a material stream connector (between points `A` and `B`). Constant pressure and temperature are assumed. +... +""" +@connector function MaterialStream(ms::MaterialSource;N_c=1, phase="unknown", flowdir=:AB, name) + if flowdir == :bidir + boundsA = boundsB = (-Inf,Inf) + elseif flowdir == :AB + boundsA = (-Inf,0) + boundsB = (0,Inf) + elseif flowdir == :BA + boundsA = (0,Inf) + boundsB = (-Inf,0) + else + error("Invalid type") + end + + @named A = State(ms;N_c=N_c,n_bnds=boundsA,phase=phase) + @named B = State(ms;N_c=N_c,n_bnds=boundsB,phase=phase) + + eqs = [ + scalarize(0.0 .~ A.nᵢ .+ B.nᵢ)... # Mass balance + A.T ~ B.T # Thermal equilibrium + A.p ~ B.p # Mechanical equilibrium + ] + + ODESystem(eqs, t; systems=[A,B], name) +end + +""" + EnergyStream(;name) + +Create an energy stream connector (between poinnts `A` and `B`). +... +""" +@connector function HeatStream(; name) + vars = @variables begin + QA(t), [description="heat A", unit="J/s", output=true] + QB(t), [description="heat B", unit="J/s", output=true] + end + + eqs = [ + 0.0 ~ QA + QB + ] + + ODESystem(eqs, t, collect(Iterators.flatten(vars)), []; name) +end + + +@connector function WorkStream(; name) + vars = @variables begin + WA(t), [description="work A", unit="J/s", output=true] + WB(t), [description="work B", unit="J/s", output=true] + end + + eqs = [ + 0.0 ~ WA + WB + ] + + ODESystem(eqs, t, collect(Iterators.flatten(vars)), []; name) +end + +""" + ControlVolume(ms::MaterialSource; ... name) + +Create a control volume with constant pressure and temperature and chemical equilibrium. +... +""", +@component function TPControlVolume(ms::MaterialSource; + N_states, + N_heats=0, + N_works=0, + reactions=Vector{ODESystem}[], + N_ph=1, phases=repeat(["unknown"],N_ph), + name) + # Init + N_c = length(ms.components) + states = [State(ms;N_c=N_c) for i in 1:N_states] + + vars = @variables begin + T(t), [description="temperature", unit="K", output=true, bounds=(0,Inf)] + p(t), [description="pressure", unit="Pa", output=true, bounds=(0,Inf)] + ϱ(t)[1:N_ph], [description="density", unit="mol/m³", output=true, bounds=(0,Inf)] + (nᵢ(t))[1:N_ph,1:N_c], [description="molar holdup", unit="mol", output=true, bounds=(0,Inf)] + n(t), [description="total molar holdup", unit="mol", output=true, bounds=(0,Inf)] + U(t), [description="internal energy", unit="J", output=true] + ΔH(t), [description="enthalpy difference inlets/outlets", unit="J/s", output=true] + ΔHᵣ(t), [description="enthalpy difference reactions", unit="J/s", output=true] + ΔE(t), [description="added/removed heat or work", unit="J/s", output=true] + end + if N_heats > 0 + @variables (Qs(t))[1:N_heats], [description="Heats added/removed", unit="J/s", output=true] + vars = vcat(vars, Qs) + else + Qs = 0.0 + end + if N_works > 0 + @variables (Ws(t))[1:N_works], [description="Work added/removed", unit="J/s", output=true] + vars = vcat(vars, Ws) + else + Ws = 0.0 + end + + eqs = [ + ΔH ~ sum([ms.VT_enthalpy(st.ϱ,st.T,st.nᵢ) for st in states]) + ΔHᵣ ~ 0.0 # TODO: Reactions -> sum([reaction.ΔH for reaction in reactions]) + ΔE ~ sum(Qs) + sum(Ws) + # Energy balance + D(U) ~ ΔH + ΔE + ΔHᵣ + # Mole balance + [D(sum(nᵢ[:,j])) ~ sum([st.nᵢ[j] for st in states]) for j in 1:N_c][:] + n ~ sum(nᵢ) + # TODO: Definition of chemical equilibrium and reactions outside of the control volume obejct? + # Thermodynamic system properties + [ϱ[i] ~ ms.molar_density(p,T,nᵢ[i,:];phase=phases[i]) for i in 1:N_ph] + U ~ sum([ms.VT_internal_energy(ϱ[i],T,nᵢ[i,:]) for i in 1:N_ph]) + ] + + ODESystem([eqs...], t, collect(Iterators.flatten(vars)), []; name, systems=sys) +end + +@component function SimpleControlVolume(ms::MaterialSource; + N_states, + N_heats=0, + N_works=0, + N_ph=1, phases=repeat(["unknown"],N_ph), + name) + # Init + N_c = length(ms.components) + states = [State(ms;N_c=N_c,name=Symbol("s$i")) for i in 1:N_states] + + vars = @variables begin + ΔH(t), [description="Enthalpy difference inlets/outlets", unit="J/s", output=true] + ΔE(t), [description="Added/removed heat or work", unit="J/s", output=true] + end + if N_heats > 0 + @variables (Qs(t))[1:N_heats], [description="Heats added/removed", unit="J/s", output=true] + vars = vcat(vars, Qs) + else + Qs = zeros(1) + end + if N_works > 0 + @variables (Ws(t))[1:N_works], [description="Work added/removed", unit="J/s", output=true] + vars = vcat(vars, Ws) + else + Ws = zeros(1) + end + + eqs = [ + ΔH ~ sum([ms.VT_enthalpy(st.ϱ,st.T,st.nᵢ)*sign(st.nᵢ[1]) for st in states]) + ΔE ~ sum([Q for Q in Qs]) + sum([W for W in Ws]) + # Energy balance + 0.0 ~ ΔH + ΔE + # Mole balance + [0.0 ~ sum([st.nᵢ[j] for st in states]) for j in 1:N_c][:] + ] + + ODESystem([eqs...], t, collect(Iterators.flatten(vars)), []; name, systems=states) +end \ No newline at end of file diff --git a/src/base/materials.jl b/src/base/materials.jl new file mode 100644 index 0000000..3a807df --- /dev/null +++ b/src/base/materials.jl @@ -0,0 +1,13 @@ +abstract type AbstractMaterialSource end + +struct MaterialSource <: AbstractMaterialSource + name::String # Name of the material source + components::Vector{String} # Component names + Mw::Vector{Float64} # Molar weight in kg/mol + pressure::Function # Pressure function (defined as f(ϱ,T,nᵢ;kwargs...)) + molar_density::Function # Molar density function (defined as f(p,T,nᵢ;kwargs...)) + VT_internal_energy::Function # Internal energy function (defined as f(ϱ,T,nᵢ;kwargs...)) + VT_enthalpy::Function # Enthalpy function (defined as f(ϱ,T,nᵢ;kwargs...)) + VT_entropy::Function # Entropy function (defined as f(ϱ,T,nᵢ;kwargs...)) + tp_flash::Function # Flash function (defined as f(p,T,nᵢ;kwargs...)) +end \ No newline at end of file diff --git a/src/base/utils.jl b/src/base/utils.jl new file mode 100644 index 0000000..4e3a819 --- /dev/null +++ b/src/base/utils.jl @@ -0,0 +1,10 @@ +function connect_states(state_1,state_2; is_stream=false) + ndir = is_stream ? -1 : 1 + con_eqs = [ + state_1.T ~ state_2.T + state_1.p ~ state_2.p + # state_1.ϱ ~ state_2.ϱ + scalarize(state_1.nᵢ .~ ndir * state_2.nᵢ)... + ] + return con_eqs +end \ No newline at end of file diff --git a/src/fluid_handling/compressors.jl b/src/fluid_handling/compressors.jl new file mode 100644 index 0000000..f6bfb58 --- /dev/null +++ b/src/fluid_handling/compressors.jl @@ -0,0 +1,27 @@ +@component function SimpleAdiabaticCompressor(ms::MaterialSource; name) + + # Subsystems + @named cv = SimpleControlVolume(ms;N_states=2,N_works=1) + @named cv_s = SimpleControlVolume(ms;N_states=2,N_works=1) + sys = [cv,cv_s] + + # Variables + vars = @variables begin + ηᴱ(t), [description="efficiency", unit="-", output=true] + end + + pars = [] + + # Equations + eqs = [ + cv.s1.p ~ cv_s.s1.p, + cv.s1.T ~ cv_s.s1.T, + cv.s2.p ~ cv_s.s2.p, + scalarize(cv.s1.nᵢ .~ cv_s.s1.nᵢ)..., + ms.VT_entropy(cv.s1.ϱ,cv.s1.T,cv.s1.nᵢ) ~ ms.VT_entropy(cv_s.s2.ϱ,cv_s.s2.T,cv_s.s2.nᵢ), + ηᴱ ~ cv.Ws[1] / cv_s.Ws[1] + ] + eqs = Equation[eqs...] + + return ODESystem([eqs...], t, vars, pars; name, systems=sys) +end \ No newline at end of file diff --git a/src/fluid_handling/heat_exchangers.jl b/src/fluid_handling/heat_exchangers.jl new file mode 100644 index 0000000..17f69cd --- /dev/null +++ b/src/fluid_handling/heat_exchangers.jl @@ -0,0 +1,17 @@ +@component function SimpleIsobaricHeatExchanger(ms::MaterialSource; name) + + # Subsystems + @named cv = SimpleControlVolume(ms;N_states=2,N_heats=1) + sys = [cv] + + # Variables + vars = @variables begin + end + + # Equations + eqs = [ + cv.s1.p ~ cv.s2.p, + ] + + return ODESystem(eqs, t, vars, []; systems=sys, name) +end \ No newline at end of file diff --git a/src/reactors/CSTR.jl b/src/reactors/CSTR.jl new file mode 100644 index 0000000..595763a --- /dev/null +++ b/src/reactors/CSTR.jl @@ -0,0 +1,8 @@ +@component function CSTR(ms:MaterialSource,reac::Reaction;name) + # Subsystems + @named cv = SimpleControlVolume(ms;) + + # Stoichiometry + + return ODESystem() +end \ No newline at end of file diff --git a/src/separation/distillation.jl b/src/separation/distillation.jl new file mode 100644 index 0000000..7a7c41d --- /dev/null +++ b/src/separation/distillation.jl @@ -0,0 +1,19 @@ +@component function DistillationColumn(ms::MaterialSource; N_stages, i_feeds, name) + # Subsystems + stages = [SimpleStage(ms; name="stage$i", add_flows=sum(i.==[1,N_stages])+sum(i.==i_feeds)) for i in 1:N_stages] + + # Connect stages + # ... + + return ODESystem(; name) +end + +@component function SimpleStage(ms::MaterialSource; name, add_flows) + # Subsystem + @named cv = TPControlVolume(ms; N_states=4) + + # VLE + # ... + + return ODESystem(; name) +end \ No newline at end of file diff --git a/test/Reactor_tests/gibbs_tests.jl b/test/Reactor_tests/gibbs_tests.jl deleted file mode 100644 index d8c9d28..0000000 --- a/test/Reactor_tests/gibbs_tests.jl +++ /dev/null @@ -1,6 +0,0 @@ -using ProcessSimulator -using Test - -y = Gibbs(2.0) - -@test y == 4.0 \ No newline at end of file diff --git a/test/base/simple_steady_state.jl b/test/base/simple_steady_state.jl new file mode 100644 index 0000000..e959d43 --- /dev/null +++ b/test/base/simple_steady_state.jl @@ -0,0 +1,93 @@ +using ProcessSimulator +using ModelingToolkit, DifferentialEquations +using ModelingToolkit: t_nounits as t + +const PS = ProcessSimulator + +# Create material sources +Mw = 0.004 # kg/mol +R = 2.1e3*Mw # J/(mol K) +cₚ = 5.2e3*Mw # J/(mol K) +cᵥ = cₚ - R + +matsource = PS.MaterialSource( + "helium", + ["helium"],[0.004], + (ϱ,T,n) -> ϱ*R*T, + (p,T,n;kwargs...) -> p/(R*Mw*T), + (ϱ,T,n) -> cᵥ*T, + (ϱ,T,n) -> cₚ*T, + (ϱ,T,n) -> cᵥ*log(T) + R*log(1/ϱ), + (p,T,n) -> NaN +) + +# Create flowsheet +# Components +comp_12 = PS.SimpleAdiabaticCompressor(matsource; name=:comp_12) +heat_22⁺ = PS.SimpleIsobaricHeatExchanger(matsource; name=:heat_22⁺) +heat_2⁺3 = PS.SimpleIsobaricHeatExchanger(matsource; name=:heat_2⁺2) +turb_34 = PS.SimpleAdiabaticCompressor(matsource; name=:turb_34) +heat_44⁺ = PS.SimpleIsobaricHeatExchanger(matsource; name=:heat_44⁺) +heat_4⁺1 = PS.SimpleIsobaricHeatExchanger(matsource; name=:heat_4⁺1) +sys = [comp_12,heat_22⁺,heat_2⁺3,turb_34,heat_44⁺,heat_4⁺1] + +# Connections +con_eqs = Equation[] +append!(con_eqs,PS.connect_states(comp_12.cv.s2,heat_22⁺.cv.s1; is_stream=true)) +append!(con_eqs,PS.connect_states(heat_22⁺.cv.s2,heat_2⁺3.cv.s1; is_stream=true)) +append!(con_eqs,PS.connect_states(heat_2⁺3.cv.s2,turb_34.cv.s1; is_stream=true)) +append!(con_eqs,PS.connect_states(turb_34.cv.s2,heat_44⁺.cv.s1; is_stream=true)) +append!(con_eqs,PS.connect_states(heat_44⁺.cv.s2,heat_4⁺1.cv.s1; is_stream=true)) + +# Additional connections (internal heat exchanger) +append!(con_eqs,[heat_22⁺.cv.Qs[1] ~ -heat_4⁺1.cv.Qs[1]]) + +@named flowsheet_ = ODESystem(con_eqs, t, [], []; systems=sys) + +giv = [ + comp_12.cv.s1.nᵢ[1] => 1.0, # T1 + comp_12.cv.s1.T => 298.15, # T1 + comp_12.cv.s1.p => 1e5, # p1 + comp_12.cv.s2.p => 1e6, # p2 + heat_22⁺.cv.s2.T => 298.15, # T2⁺ + turb_34.cv.s1.T => 253.15, # T3 + turb_34.cv.s2.p => 1e5, # p4 + heat_44⁺.cv.s2.T => 253.15, # T4⁺ + comp_12.ηᴱ => 1.0, # ηᴱ + turb_34.ηᴱ => 1.0 # ηᴱ +] +unk = [ + heat_44⁺.cv.Qs[1], + comp_12.cv.Ws[1], + turb_34.cv.Ws[1] +] + +flowsheet,idx = structural_simplify(flowsheet_, (first.(giv), first.(unk))) + +u0 = [ + comp_12.cv.s2.T => 1.0, + comp_12.cv_s.s2.T => 1.0, + (comp_12.cv.Ws)[1] => -1.0, + turb_34.cv.s2.T => 1.0, + turb_34.cv_s.s2.T => 1.0, + (turb_34.cv.Ws)[1] => 1.0, + heat_4⁺1.cv.s2.T => 1.0, +] + +prob = SteadyStateProblem(flowsheet,giv,u0) +sol = solve(prob) + +ε_KM = abs(sol[heat_44⁺.cv.Qs[1]])/abs(sol[turb_34.cv.Ws[1]] + sol[comp_12.cv.Ws[1]]) + +(T₁,T₂,T₂₊,T₃,T₄,T₄₊) = ( + sol.ps[comp_12.cv.s1.T], + sol[comp_12.cv.s2.T], + sol.ps[heat_22⁺.cv.s2.T], + sol.ps[turb_34.cv.s1.T], + sol[turb_34.cv.s2.T], + sol.ps[heat_44⁺.cv.s2.T] +) + +@test round(T₂,digits=2) ≈ 755.58 +@test round(T₄,digits=2) ≈ 99.89 +@test round(ε_KM,digits=3) ≈ 0.504 \ No newline at end of file diff --git a/test/reactors/simple_cstr.jl b/test/reactors/simple_cstr.jl new file mode 100644 index 0000000..11d5f2c --- /dev/null +++ b/test/reactors/simple_cstr.jl @@ -0,0 +1,22 @@ +using ProcessSimulator +using ModelingToolkit, DifferentialEquations +using ModelingToolkit: t_nounits as t + +const PS = ProcessSimulator + +# Create material sources +Mw = 0.004 # kg/mol +R = 2.1e3*Mw # J/(mol K) +cₚ = 5.2e3*Mw # J/(mol K) +cᵥ = cₚ - R + +matsource = PS.MaterialSource( + "helium", + ["helium"],[0.004], + (ϱ,T,n) -> ϱ*R*T, + (p,T,n;kwargs...) -> p/(R*Mw*T), + (ϱ,T,n) -> cᵥ*T, + (ϱ,T,n) -> cₚ*T, + (ϱ,T,n) -> cᵥ*log(T) + R*log(1/ϱ), + (p,T,n) -> NaN +) \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index c7b573c..cc92390 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,6 +2,6 @@ using ProcessSimulator using Test using SafeTestsets -@safetestset "Gibbs reactor" begin - include("Reactor_tests/gibbs_tests.jl") +@safetestset "Simple Steady State" begin + include("base/simple_steady_state.jl") end From d6631cf58fc122f352a2091196bc5d0403efdf05 Mon Sep 17 00:00:00 2001 From: se-schmitt Date: Mon, 14 Oct 2024 16:04:19 +0200 Subject: [PATCH 2/7] Rewrite components as functions --- src/base/base_components.jl | 140 +++++++++--------------------- src/base/utils.jl | 10 --- src/fluid_handling/compressors.jl | 22 +++-- 3 files changed, 50 insertions(+), 122 deletions(-) diff --git a/src/base/base_components.jl b/src/base/base_components.jl index 66973b1..12588ac 100644 --- a/src/base/base_components.jl +++ b/src/base/base_components.jl @@ -1,96 +1,56 @@ -# Base components +@component function MaterialFlow(ms::MaterialSource;N_c=1,phase="unknown",name) -@connector function State(ms::MaterialSource;N_c=1,n_bnds=(-Inf,Inf),phase="unknown",name) - vars = @variables begin - T(t), [description="temperature", unit="K", output=true] - p(t), [description="pressure", unit="Pa", output=true] - ϱ(t), [description="density", unit="mol/m³", output=true] - (nᵢ(t))[1:N_c], [description="molar flow", unit="mol/s", bounds=n_bnds, output=true] - n(t), [description="total molar flow", unit="mol/s", output=true] - end + @named s = TϱState(N_c=N_c) + @named c = PHxConnector(ms;N_c=N_c,phase=phase) + sub = [s,c] eqs = [ - # Density - ϱ ~ ms.molar_density(p,T,nᵢ;phase=phase) - n ~ sum(nᵢ) + c.p ~ ms.pressure(s.ϱ,s.T,c.xᵢ) + c.h ~ ms.VT_enthalpy(s.ϱ,s.T,c.xᵢ) + 1.0 ~ sum([xᵢ for xᵢ in c.xᵢ]) + scalarize(c.xᵢ .~ s.nᵢ ./ c.n) ] - pars = @parameters begin + ODESystem(eqs, t, [], []; name, systems=sub) +end + +@component function TϱState(;N_c=1, name) + vars = @variables begin + T(t), [description="temperature", output=true] #, unit=u"K"] + ϱ(t), [description="density", output=true] #, unit=u"mol m^-3"] + nᵢ(t)[1:N_c], [description="mole fractions", output=true] #, unit=u"mol s^-1"] end - ODESystem(eqs, t, collect(Iterators.flatten(vars)), pars; name) + return ODESystem(Equation[], t, collect(Iterators.flatten(vars)), []; name) end -""" - MaterialStream(ms::MaterialSource;N_c=1, phase="unknown", flowdir=:generic, name) - -Create a material stream connector (between points `A` and `B`). Constant pressure and temperature are assumed. -... -""" -@connector function MaterialStream(ms::MaterialSource;N_c=1, phase="unknown", flowdir=:AB, name) - if flowdir == :bidir - boundsA = boundsB = (-Inf,Inf) - elseif flowdir == :AB - boundsA = (-Inf,0) - boundsB = (0,Inf) - elseif flowdir == :BA - boundsA = (0,Inf) - boundsB = (-Inf,0) - else - error("Invalid type") +@connector function PHxConnector(ms::MaterialSource;N_c=1,phase="unknown",name) + vars = @variables begin + p(t), [description="pressure"] #, unit=u"Pa"] + h(t), [description="molar enthalpy", connect=Stream] #, unit=u"J mol^-1"] + xᵢ(t)[1:N_c]=1/N_c, [description="mole fractions", connect=Stream] #, unit=u"mol mol^-1"] + n(t), [description="total molar flow", connect=Flow] #, unit=u"mol s^-1"] end - @named A = State(ms;N_c=N_c,n_bnds=boundsA,phase=phase) - @named B = State(ms;N_c=N_c,n_bnds=boundsB,phase=phase) - - eqs = [ - scalarize(0.0 .~ A.nᵢ .+ B.nᵢ)... # Mass balance - A.T ~ B.T # Thermal equilibrium - A.p ~ B.p # Mechanical equilibrium - ] - - ODESystem(eqs, t; systems=[A,B], name) + return ODESystem(Equation[], t, collect(Iterators.flatten(vars)), []; name) end -""" - EnergyStream(;name) - -Create an energy stream connector (between poinnts `A` and `B`). -... -""" -@connector function HeatStream(; name) +@connector function HeatConnector(;name) vars = @variables begin - QA(t), [description="heat A", unit="J/s", output=true] - QB(t), [description="heat B", unit="J/s", output=true] + Q(t)#, [description="heat flux", connect=Flow] #, unit=u"J s^-1"] end - eqs = [ - 0.0 ~ QA + QB - ] - - ODESystem(eqs, t, collect(Iterators.flatten(vars)), []; name) + return ODESystem(Equation[], t, collect(Iterators.flatten(vars)), []; name) end - -@connector function WorkStream(; name) +@connector function WorkConnector(;name) vars = @variables begin - WA(t), [description="work A", unit="J/s", output=true] - WB(t), [description="work B", unit="J/s", output=true] + W(t)#, [description="power", connect=Flow] #, unit=u"J s^1"] end - eqs = [ - 0.0 ~ WA + WB - ] - - ODESystem(eqs, t, collect(Iterators.flatten(vars)), []; name) + return ODESystem(Equation[], t, collect(Iterators.flatten(vars)), []; name) end -""" - ControlVolume(ms::MaterialSource; ... name) - -Create a control volume with constant pressure and temperature and chemical equilibrium. -... -""", @component function TPControlVolume(ms::MaterialSource; N_states, N_heats=0, @@ -100,7 +60,9 @@ Create a control volume with constant pressure and temperature and chemical equi name) # Init N_c = length(ms.components) - states = [State(ms;N_c=N_c) for i in 1:N_states] + flows = [MaterialFlow(ms;N_c=N_c,name=Symbol("s$i")) for i in 1:N_states] + 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", unit="K", output=true, bounds=(0,Inf)] @@ -113,23 +75,11 @@ Create a control volume with constant pressure and temperature and chemical equi ΔHᵣ(t), [description="enthalpy difference reactions", unit="J/s", output=true] ΔE(t), [description="added/removed heat or work", unit="J/s", output=true] end - if N_heats > 0 - @variables (Qs(t))[1:N_heats], [description="Heats added/removed", unit="J/s", output=true] - vars = vcat(vars, Qs) - else - Qs = 0.0 - end - if N_works > 0 - @variables (Ws(t))[1:N_works], [description="Work added/removed", unit="J/s", output=true] - vars = vcat(vars, Ws) - else - Ws = 0.0 - end eqs = [ ΔH ~ sum([ms.VT_enthalpy(st.ϱ,st.T,st.nᵢ) for st in states]) ΔHᵣ ~ 0.0 # TODO: Reactions -> sum([reaction.ΔH for reaction in reactions]) - ΔE ~ sum(Qs) + sum(Ws) + Δ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 + ΔHᵣ # Mole balance @@ -152,32 +102,22 @@ end name) # Init N_c = length(ms.components) - states = [State(ms;N_c=N_c,name=Symbol("s$i")) for i in 1:N_states] + flows = [MaterialFlow(ms;N_c=N_c,name=Symbol("s$i")) for i in 1:N_states] + 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="J/s", output=true] ΔE(t), [description="Added/removed heat or work", unit="J/s", output=true] end - if N_heats > 0 - @variables (Qs(t))[1:N_heats], [description="Heats added/removed", unit="J/s", output=true] - vars = vcat(vars, Qs) - else - Qs = zeros(1) - end - if N_works > 0 - @variables (Ws(t))[1:N_works], [description="Work added/removed", unit="J/s", output=true] - vars = vcat(vars, Ws) - else - Ws = zeros(1) - end eqs = [ - ΔH ~ sum([ms.VT_enthalpy(st.ϱ,st.T,st.nᵢ)*sign(st.nᵢ[1]) for st in states]) - ΔE ~ sum([Q for Q in Qs]) + sum([W for W in Ws]) + ΔH ~ sum([f.c.h*f.c.n for f in flows]) + Δ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([st.nᵢ[j] for st in states]) for j in 1:N_c][:] + [0.0 ~ sum([f.s.nᵢ[j] for f in flows]) for j in 1:N_c][:] ] ODESystem([eqs...], t, collect(Iterators.flatten(vars)), []; name, systems=states) diff --git a/src/base/utils.jl b/src/base/utils.jl index 4e3a819..e69de29 100644 --- a/src/base/utils.jl +++ b/src/base/utils.jl @@ -1,10 +0,0 @@ -function connect_states(state_1,state_2; is_stream=false) - ndir = is_stream ? -1 : 1 - con_eqs = [ - state_1.T ~ state_2.T - state_1.p ~ state_2.p - # state_1.ϱ ~ state_2.ϱ - scalarize(state_1.nᵢ .~ ndir * state_2.nᵢ)... - ] - return con_eqs -end \ No newline at end of file diff --git a/src/fluid_handling/compressors.jl b/src/fluid_handling/compressors.jl index f6bfb58..2e844c4 100644 --- a/src/fluid_handling/compressors.jl +++ b/src/fluid_handling/compressors.jl @@ -2,26 +2,24 @@ # Subsystems @named cv = SimpleControlVolume(ms;N_states=2,N_works=1) - @named cv_s = SimpleControlVolume(ms;N_states=2,N_works=1) - sys = [cv,cv_s] + sys = [cv] # Variables vars = @variables begin - ηᴱ(t), [description="efficiency", unit="-", output=true] + T2_s(t), [description="temperature", output=true] + ϱ2_s(t), [description="density", output=true] end - pars = [] + pars = @parameters begin + ηᴱ, [description="efficiency", output=true] + end # Equations eqs = [ - cv.s1.p ~ cv_s.s1.p, - cv.s1.T ~ cv_s.s1.T, - cv.s2.p ~ cv_s.s2.p, - scalarize(cv.s1.nᵢ .~ cv_s.s1.nᵢ)..., - ms.VT_entropy(cv.s1.ϱ,cv.s1.T,cv.s1.nᵢ) ~ ms.VT_entropy(cv_s.s2.ϱ,cv_s.s2.T,cv_s.s2.nᵢ), - ηᴱ ~ cv.Ws[1] / cv_s.Ws[1] + ms.VT_entropy(cv.f1.s.ϱ,cv.f1.s.T,cv.f1.c.xᵢ) ~ ms.VT_entropy(cv.f2.s.ϱ,cv.f2.s.T,cv.f2.c.xᵢ) + ϱ2_s ~ ms.molar_density(cv.f2.c.p,T2_s,cv.f2.c.xᵢ) + ηᴱ ~ cv.Ws[1] / (ms.VT_enthalpy(ϱ2_s,T2_s,cv.f2.c.xᵢ) - ms.VT_enthalpy(cv.f1.s.ϱ,cv.f1.s.T,cv.f1.c.xᵢ)) ] - eqs = Equation[eqs...] - return ODESystem([eqs...], t, vars, pars; name, systems=sys) + return ODESystem(eqs, t, vars, pars; name, systems=sys) end \ No newline at end of file From b0abd88af181208a8eae9201fc8b1df722444da5 Mon Sep 17 00:00:00 2001 From: se-schmitt Date: Sat, 19 Oct 2024 12:38:30 +0200 Subject: [PATCH 3/7] Base components with test --- src/base/base_components.jl | 157 +++++++++++++------------ src/base/materials.jl | 40 ++++++- src/base/utils.jl | 41 +++++++ src/fluid_handling/compressors.jl | 19 ++-- src/fluid_handling/heat_exchangers.jl | 11 +- test/base/simple_steady_state.jl | 158 +++++++++++++++----------- test/reactors/simple_cstr.jl | 19 ---- 7 files changed, 267 insertions(+), 178 deletions(-) diff --git a/src/base/base_components.jl b/src/base/base_components.jl index 12588ac..3fc77b5 100644 --- a/src/base/base_components.jl +++ b/src/base/base_components.jl @@ -1,124 +1,141 @@ -@component function MaterialFlow(ms::MaterialSource;N_c=1,phase="unknown",name) - @named s = TϱState(N_c=N_c) - @named c = PHxConnector(ms;N_c=N_c,phase=phase) - sub = [s,c] +@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"] + s(t), [description="molar entropy"] #, unit=u"J mol^-1 K^-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 = [ - c.p ~ ms.pressure(s.ϱ,s.T,c.xᵢ) - c.h ~ ms.VT_enthalpy(s.ϱ,s.T,c.xᵢ) - 1.0 ~ sum([xᵢ for xᵢ in c.xᵢ]) - scalarize(c.xᵢ .~ s.nᵢ ./ c.n) + # EOS + ϱ ~ ms.molar_density(p,T,xᵢ;phase=phase), + h ~ ms.VT_enthalpy(ϱ,T,xᵢ), + s ~ ms.VT_entropy(ϱ,T,xᵢ), + 1.0 ~ sum([xᵢ for xᵢ in xᵢ]), + # Connector "1" + T ~ c1.T, + ϱ ~ c1.ϱ, + p ~ c1.p, + h ~ c1.h, + s ~ c1.s, + n ~ c1.n, + scalarize(xᵢ .~ c1.xᵢ)..., + # Connector "2" + T ~ c2.T, + ϱ ~ c2.ϱ, + p ~ c2.p, + h ~ c2.h, + s ~ c2.s, + n ~ -c2.n, + scalarize(xᵢ .~ c2.xᵢ)..., ] - ODESystem(eqs, t, [], []; name, systems=sub) + return ODESystem(eqs, t, collect(Iterators.flatten(vars)), []; name, systems=mcons) end -@component function TϱState(;N_c=1, name) - vars = @variables begin - T(t), [description="temperature", output=true] #, unit=u"K"] - ϱ(t), [description="density", output=true] #, unit=u"mol m^-3"] - nᵢ(t)[1:N_c], [description="mole fractions", output=true] #, unit=u"mol s^-1"] - end - - return ODESystem(Equation[], t, collect(Iterators.flatten(vars)), []; name) -end - -@connector function PHxConnector(ms::MaterialSource;N_c=1,phase="unknown",name) +@component function MaterialConnector(ms::MaterialSource;name) 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", connect=Stream] #, unit=u"J mol^-1"] - xᵢ(t)[1:N_c]=1/N_c, [description="mole fractions", connect=Stream] #, unit=u"mol mol^-1"] + s(t), [description="molar entropy", connect=Stream] #, unit=u"J mol^-1 K^-1"] + xᵢ(t)[1:ms.N_c], [description="mole fractions", connect=Stream] #, 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 -@connector function HeatConnector(;name) +@component function HeatConnector(;name) vars = @variables begin - Q(t)#, [description="heat flux", connect=Flow] #, unit=u"J s^-1"] + Q(t), [description="heat flux"] #, unit=u"J s^-1"] end - return ODESystem(Equation[], t, collect(Iterators.flatten(vars)), []; name) + return ODESystem(Equation[], t, vars, []; name) end -@connector function WorkConnector(;name) +@component function WorkConnector(;name) vars = @variables begin - W(t)#, [description="power", connect=Flow] #, unit=u"J s^1"] + W(t), [description="power"] #, unit=u"J s^1"] end - return ODESystem(Equation[], t, collect(Iterators.flatten(vars)), []; name) + return ODESystem(Equation[], t, vars, []; name) end -@component function TPControlVolume(ms::MaterialSource; - N_states, - N_heats=0, - N_works=0, - reactions=Vector{ODESystem}[], - N_ph=1, phases=repeat(["unknown"],N_ph), - name) +@component function SimpleControlVolume(ms::MaterialSource; + N_mcons, + N_heats=0, + N_works=0, + name) # Init - N_c = length(ms.components) - flows = [MaterialFlow(ms;N_c=N_c,name=Symbol("s$i")) for i in 1:N_states] + 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", unit="K", output=true, bounds=(0,Inf)] - p(t), [description="pressure", unit="Pa", output=true, bounds=(0,Inf)] - ϱ(t)[1:N_ph], [description="density", unit="mol/m³", output=true, bounds=(0,Inf)] - (nᵢ(t))[1:N_ph,1:N_c], [description="molar holdup", unit="mol", output=true, bounds=(0,Inf)] - n(t), [description="total molar holdup", unit="mol", output=true, bounds=(0,Inf)] - U(t), [description="internal energy", unit="J", output=true] - ΔH(t), [description="enthalpy difference inlets/outlets", unit="J/s", output=true] - ΔHᵣ(t), [description="enthalpy difference reactions", unit="J/s", output=true] - ΔE(t), [description="added/removed heat or work", unit="J/s", output=true] + Δ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 = [ - ΔH ~ sum([ms.VT_enthalpy(st.ϱ,st.T,st.nᵢ) for st in states]) - ΔHᵣ ~ 0.0 # TODO: Reactions -> sum([reaction.ΔH for reaction in reactions]) + 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 - D(U) ~ ΔH + ΔE + ΔHᵣ + 0.0 ~ ΔH + ΔE # Mole balance - [D(sum(nᵢ[:,j])) ~ sum([st.nᵢ[j] for st in states]) for j in 1:N_c][:] - n ~ sum(nᵢ) - # TODO: Definition of chemical equilibrium and reactions outside of the control volume obejct? - # Thermodynamic system properties - [ϱ[i] ~ ms.molar_density(p,T,nᵢ[i,:];phase=phases[i]) for i in 1:N_ph] - U ~ sum([ms.VT_internal_energy(ϱ[i],T,nᵢ[i,:]) for i in 1:N_ph]) + [0.0 ~ sum([c.xᵢ[j]*c.n for c in mcons]) for j in 1:ms.N_c][:] ] - ODESystem([eqs...], t, collect(Iterators.flatten(vars)), []; name, systems=sys) + return ODESystem(eqs, t, vars, []; name, systems=[mcons...,works...,heats...]) end -@component function SimpleControlVolume(ms::MaterialSource; - N_states, - N_heats=0, - N_works=0, - N_ph=1, phases=repeat(["unknown"],N_ph), - name) +@component function TPControlVolume(ms::MaterialSource; + N_mcons, + N_heats=0, + N_works=0, + N_ph=1, phases=repeat(["unknown"],N_ph), + name) # Init - N_c = length(ms.components) - flows = [MaterialFlow(ms;N_c=N_c,name=Symbol("s$i")) for i in 1:N_states] + 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="J/s", output=true] - ΔE(t), [description="Added/removed heat or work", unit="J/s", output=true] + 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"] + 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"] + end + + pars = @parameters begin + V, [description="volume", bounds=(0,Inf)] #, unit=u"m^3"] end eqs = [ - ΔH ~ sum([f.c.h*f.c.n for f in flows]) + Δ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 + D(U) ~ ΔH + ΔE + ΔHᵣ # Mole balance - [0.0 ~ sum([f.s.nᵢ[j] for f in flows]) for j in 1:N_c][:] + [D(sum(nᵢ[:,j])) ~ sum([c.xᵢ[j]*c.n for c in mcons]) for j in 1:ms.N_c][:] + n ~ sum(nᵢ) + # Thermodynamic system properties + [ϱ[i] ~ ms.molar_density(p,T,nᵢ[i,:];phase=phases[i]) for i in 1:N_ph] + U ~ sum([ms.VT_internal_energy(ϱ[i],T,nᵢ[i,:]) for i in 1:N_ph]) ] - ODESystem([eqs...], t, collect(Iterators.flatten(vars)), []; name, systems=states) + return ODESystem(eqs, t, collect(Iterators.flatten(vars)), pars; name, systems=[mcons...,works...,heats...]) end \ No newline at end of file diff --git a/src/base/materials.jl b/src/base/materials.jl index 3a807df..c707401 100644 --- a/src/base/materials.jl +++ b/src/base/materials.jl @@ -3,11 +3,39 @@ abstract type AbstractMaterialSource 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,nᵢ;kwargs...)) - molar_density::Function # Molar density function (defined as f(p,T,nᵢ;kwargs...)) - VT_internal_energy::Function # Internal energy function (defined as f(ϱ,T,nᵢ;kwargs...)) - VT_enthalpy::Function # Enthalpy function (defined as f(ϱ,T,nᵢ;kwargs...)) - VT_entropy::Function # Entropy function (defined as f(ϱ,T,nᵢ;kwargs...)) - tp_flash::Function # Flash function (defined as f(p,T,nᵢ;kwargs...)) + pressure::Function # Pressure function (defined as f(ϱ,T,xᵢ;kwargs...)) + molar_density::Function # Molar density function (defined as f(p,T,xᵢ;kwargs...)) + VT_internal_energy::Function # Internal energy function (defined as f(ϱ,T,xᵢ;kwargs...)) + VT_enthalpy::Function # Enthalpy function (defined as f(ϱ,T,xᵢ;kwargs...)) + VT_entropy::Function # Entropy function (defined as f(ϱ,T,xᵢ;kwargs...)) + tp_flash::Function # Flash function (defined as f(p,T,xᵢ;kwargs...)) +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)) + ) end \ No newline at end of file diff --git a/src/base/utils.jl b/src/base/utils.jl index e69de29..70fda72 100644 --- a/src/base/utils.jl +++ b/src/base/utils.jl @@ -0,0 +1,41 @@ +@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([xᵢ for xᵢ in xᵢ]), + c.h ~ ms.VT_enthalpy(ϱ,T,xᵢ), + c.s ~ ms.VT_entropy(ϱ,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 + +function mconnect(c1,c2) + Equation[ + c1.T ~ c2.T, + c1.ϱ ~ c2.ϱ, + c1.p ~ c2.p, + c1.h ~ c2.h, + c1.s ~ c2.s, + scalarize(c1.xᵢ .~ c2.xᵢ)..., + 0.0 ~ c1.n + c2.n, + ] +end \ No newline at end of file diff --git a/src/fluid_handling/compressors.jl b/src/fluid_handling/compressors.jl index 2e844c4..4f80798 100644 --- a/src/fluid_handling/compressors.jl +++ b/src/fluid_handling/compressors.jl @@ -1,25 +1,26 @@ @component function SimpleAdiabaticCompressor(ms::MaterialSource; name) # Subsystems - @named cv = SimpleControlVolume(ms;N_states=2,N_works=1) - sys = [cv] + @named cv = SimpleControlVolume(ms;N_mcons=2,N_works=1) # Variables vars = @variables begin - T2_s(t), [description="temperature", output=true] - ϱ2_s(t), [description="density", output=true] + 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 - ηᴱ, [description="efficiency", output=true] + ηᴱ = 1.0, [description="efficiency"] end # Equations eqs = [ - ms.VT_entropy(cv.f1.s.ϱ,cv.f1.s.T,cv.f1.c.xᵢ) ~ ms.VT_entropy(cv.f2.s.ϱ,cv.f2.s.T,cv.f2.c.xᵢ) - ϱ2_s ~ ms.molar_density(cv.f2.c.p,T2_s,cv.f2.c.xᵢ) - ηᴱ ~ cv.Ws[1] / (ms.VT_enthalpy(ϱ2_s,T2_s,cv.f2.c.xᵢ) - ms.VT_enthalpy(cv.f1.s.ϱ,cv.f1.s.T,cv.f1.c.xᵢ)) + 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=sys) + return ODESystem(eqs, t, vars, pars; name, systems=[cv]) end \ No newline at end of file diff --git a/src/fluid_handling/heat_exchangers.jl b/src/fluid_handling/heat_exchangers.jl index 17f69cd..f2affdd 100644 --- a/src/fluid_handling/heat_exchangers.jl +++ b/src/fluid_handling/heat_exchangers.jl @@ -1,17 +1,18 @@ @component function SimpleIsobaricHeatExchanger(ms::MaterialSource; name) # Subsystems - @named cv = SimpleControlVolume(ms;N_states=2,N_heats=1) - sys = [cv] + @named cv = SimpleControlVolume(ms;N_mcons=2,N_heats=1) # Variables vars = @variables begin - end + Q(t), [description="heat flux"] #, unit=u"J s^-1"] + end # Equations eqs = [ - cv.s1.p ~ cv.s2.p, + cv.q1.Q ~ Q, + cv.c1.p ~ cv.c2.p, ] - return ODESystem(eqs, t, vars, []; systems=sys, name) + return ODESystem(eqs, t, vars, []; systems=[cv], name) end \ No newline at end of file diff --git a/test/base/simple_steady_state.jl b/test/base/simple_steady_state.jl index e959d43..6074479 100644 --- a/test/base/simple_steady_state.jl +++ b/test/base/simple_steady_state.jl @@ -10,84 +10,104 @@ R = 2.1e3*Mw # J/(mol K) cₚ = 5.2e3*Mw # J/(mol K) cᵥ = cₚ - R -matsource = PS.MaterialSource( - "helium", - ["helium"],[0.004], - (ϱ,T,n) -> ϱ*R*T, - (p,T,n;kwargs...) -> p/(R*Mw*T), - (ϱ,T,n) -> cᵥ*T, - (ϱ,T,n) -> cₚ*T, - (ϱ,T,n) -> cᵥ*log(T) + R*log(1/ϱ), - (p,T,n) -> NaN +# Perfect gas equations +matsource = PS.MaterialSource("helium"; + Mw = 0.004, + molar_density = (p, T, x; kwargs...) -> p/(R*Mw*T), + VT_internal_energy = (ϱ, T, x) -> cᵥ*T , + VT_enthalpy = (ϱ, T, x) -> cₚ*T, + VT_entropy = (ϱ, T, x) -> cᵥ*log(T) + R*log(1/ϱ) ) -# Create flowsheet -# Components -comp_12 = PS.SimpleAdiabaticCompressor(matsource; name=:comp_12) -heat_22⁺ = PS.SimpleIsobaricHeatExchanger(matsource; name=:heat_22⁺) -heat_2⁺3 = PS.SimpleIsobaricHeatExchanger(matsource; name=:heat_2⁺2) -turb_34 = PS.SimpleAdiabaticCompressor(matsource; name=:turb_34) -heat_44⁺ = PS.SimpleIsobaricHeatExchanger(matsource; name=:heat_44⁺) -heat_4⁺1 = PS.SimpleIsobaricHeatExchanger(matsource; name=:heat_4⁺1) -sys = [comp_12,heat_22⁺,heat_2⁺3,turb_34,heat_44⁺,heat_4⁺1] - -# Connections -con_eqs = Equation[] -append!(con_eqs,PS.connect_states(comp_12.cv.s2,heat_22⁺.cv.s1; is_stream=true)) -append!(con_eqs,PS.connect_states(heat_22⁺.cv.s2,heat_2⁺3.cv.s1; is_stream=true)) -append!(con_eqs,PS.connect_states(heat_2⁺3.cv.s2,turb_34.cv.s1; is_stream=true)) -append!(con_eqs,PS.connect_states(turb_34.cv.s2,heat_44⁺.cv.s1; is_stream=true)) -append!(con_eqs,PS.connect_states(heat_44⁺.cv.s2,heat_4⁺1.cv.s1; is_stream=true)) - -# Additional connections (internal heat exchanger) -append!(con_eqs,[heat_22⁺.cv.Qs[1] ~ -heat_4⁺1.cv.Qs[1]]) - -@named flowsheet_ = ODESystem(con_eqs, t, [], []; systems=sys) - -giv = [ - comp_12.cv.s1.nᵢ[1] => 1.0, # T1 - comp_12.cv.s1.T => 298.15, # T1 - comp_12.cv.s1.p => 1e5, # p1 - comp_12.cv.s2.p => 1e6, # p2 - heat_22⁺.cv.s2.T => 298.15, # T2⁺ - turb_34.cv.s1.T => 253.15, # T3 - turb_34.cv.s2.p => 1e5, # p4 - heat_44⁺.cv.s2.T => 253.15, # T4⁺ - comp_12.ηᴱ => 1.0, # ηᴱ - turb_34.ηᴱ => 1.0 # ηᴱ -] -unk = [ - heat_44⁺.cv.Qs[1], - comp_12.cv.Ws[1], - turb_34.cv.Ws[1] +# Compressor +@named inlet = PS.Port(matsource) +@named comp = PS.SimpleAdiabaticCompressor(matsource) +@named outlet = PS.Port(matsource) + +con_eqs_comp = [ + PS.mconnect(inlet.c, comp.cv.c1)..., + PS.mconnect(comp.cv.c2, outlet.c)... ] +@named compressor_ = ODESystem(con_eqs_comp, t, [], []; systems=[inlet,comp,outlet]) -flowsheet,idx = structural_simplify(flowsheet_, (first.(giv), first.(unk))) +inp_comp = [ + inlet.m => 1.0, + inlet.T => 300., + inlet.p => 1e5, + outlet.p => 1e6 +] +out_comp = [comp.W] +compressor,idx = structural_simplify(compressor_, (first.(inp_comp), out_comp)) -u0 = [ - comp_12.cv.s2.T => 1.0, - comp_12.cv_s.s2.T => 1.0, - (comp_12.cv.Ws)[1] => -1.0, - turb_34.cv.s2.T => 1.0, - turb_34.cv_s.s2.T => 1.0, - (turb_34.cv.Ws)[1] => 1.0, - heat_4⁺1.cv.s2.T => 1.0, +u0_comp = [ + (u => 1. for u in unknowns(compressor))..., + comp.ηᴱ => 0.8 ] -prob = SteadyStateProblem(flowsheet,giv,u0) -sol = solve(prob) +prob_comp = SteadyStateProblem(compressor,inp_comp,u0_comp) +sol_comp = solve(prob_comp) -ε_KM = abs(sol[heat_44⁺.cv.Qs[1]])/abs(sol[turb_34.cv.Ws[1]] + sol[comp_12.cv.Ws[1]]) +@test sol_comp[comp.W] ≈ 2.3933999e6 rtol=1e-5 -(T₁,T₂,T₂₊,T₃,T₄,T₄₊) = ( - sol.ps[comp_12.cv.s1.T], - sol[comp_12.cv.s2.T], - sol.ps[heat_22⁺.cv.s2.T], - sol.ps[turb_34.cv.s1.T], - sol[turb_34.cv.s2.T], - sol.ps[heat_44⁺.cv.s2.T] +# Create flowsheet +@named comp_12 = PS.SimpleAdiabaticCompressor(matsource) +@named s2 = PS.MaterialStream(matsource) +@named heat_22⁺ = PS.SimpleIsobaricHeatExchanger(matsource) +@named s2⁺ = PS.MaterialStream(matsource) +@named heat_2⁺3 = PS.SimpleIsobaricHeatExchanger(matsource) +@named s3 = PS.MaterialStream(matsource) +@named turb_34 = PS.SimpleAdiabaticCompressor(matsource) +@named s4 = PS.MaterialStream(matsource) +@named heat_44⁺ = PS.SimpleIsobaricHeatExchanger(matsource) +@named s4⁺ = PS.MaterialStream(matsource) +@named heat_4⁺1 = PS.SimpleIsobaricHeatExchanger(matsource) + +systems = [comp_12,heat_22⁺,heat_2⁺3,turb_34,heat_44⁺,heat_4⁺1] +streams = [inlet,s2,s2⁺,s3,s4,s4⁺,outlet] + +con_eqs = vcat( + PS.mconnect(inlet.c, comp_12.cv.c1)..., + PS.mconnect(comp_12.cv.c2, s2.c1)..., + PS.mconnect(s2.c2, heat_22⁺.cv.c1)..., + PS.mconnect(heat_22⁺.cv.c2, s2⁺.c1)..., + PS.mconnect(s2⁺.c2, heat_2⁺3.cv.c1)..., + PS.mconnect(heat_2⁺3.cv.c2, s3.c1)..., + PS.mconnect(s3.c2, turb_34.cv.c1)..., + PS.mconnect(turb_34.cv.c2, s4.c1)..., + PS.mconnect(s4.c2, heat_44⁺.cv.c1)..., + PS.mconnect(heat_44⁺.cv.c2, s4⁺.c1)..., + PS.mconnect(s4⁺.c2, heat_4⁺1.cv.c1)..., + PS.mconnect(heat_4⁺1.cv.c2, outlet.c)... ) +append!(con_eqs,[0.0 ~ heat_2⁺3.Q + heat_4⁺1.Q]) + +@named flowsheet_ = ODESystem(con_eqs, t, [], []; systems=[systems...,streams...]) + +inp = [ + inlet.n => 1.0, + inlet.T => 298.15, + inlet.p => 1e5, + s2.p => 1e6, + s2⁺.T => 298.15, + s3.T => 253.15, + outlet.T => 298.15, + outlet.p => 1e5, +] +out = [ + heat_44⁺.Q, + comp_12.W, + turb_34.W, +] + +flowsheet,idx = structural_simplify(flowsheet_, (first.(inp), out)) + +u0 = [u => 1. for u in unknowns(flowsheet)] + +prob = SteadyStateProblem(flowsheet,inp,u0) +sol = solve(prob) + +ε_KM = abs(sol[heat_44⁺.Q])/abs(sol[turb_34.W] + sol[comp_12.W]) -@test round(T₂,digits=2) ≈ 755.58 -@test round(T₄,digits=2) ≈ 99.89 +@test round(sol[s2.T],digits=2) ≈ 755.58 +@test round(sol[s4.T],digits=2) ≈ 99.89 @test round(ε_KM,digits=3) ≈ 0.504 \ No newline at end of file diff --git a/test/reactors/simple_cstr.jl b/test/reactors/simple_cstr.jl index 11d5f2c..62635af 100644 --- a/test/reactors/simple_cstr.jl +++ b/test/reactors/simple_cstr.jl @@ -1,22 +1,3 @@ using ProcessSimulator using ModelingToolkit, DifferentialEquations using ModelingToolkit: t_nounits as t - -const PS = ProcessSimulator - -# Create material sources -Mw = 0.004 # kg/mol -R = 2.1e3*Mw # J/(mol K) -cₚ = 5.2e3*Mw # J/(mol K) -cᵥ = cₚ - R - -matsource = PS.MaterialSource( - "helium", - ["helium"],[0.004], - (ϱ,T,n) -> ϱ*R*T, - (p,T,n;kwargs...) -> p/(R*Mw*T), - (ϱ,T,n) -> cᵥ*T, - (ϱ,T,n) -> cₚ*T, - (ϱ,T,n) -> cᵥ*log(T) + R*log(1/ϱ), - (p,T,n) -> NaN -) \ No newline at end of file From 0aa4d7b9c521db3873751d1db40647b6bf703b7a Mon Sep 17 00:00:00 2001 From: se-schmitt Date: Sat, 19 Oct 2024 13:01:41 +0200 Subject: [PATCH 4/7] Use MTK connectors --- src/base/base_components.jl | 12 ++++++------ src/base/utils.jl | 12 ------------ test/base/simple_steady_state.jl | 28 ++++++++++++++-------------- 3 files changed, 20 insertions(+), 32 deletions(-) diff --git a/src/base/base_components.jl b/src/base/base_components.jl index 3fc77b5..e8d4c6a 100644 --- a/src/base/base_components.jl +++ b/src/base/base_components.jl @@ -41,14 +41,14 @@ return ODESystem(eqs, t, collect(Iterators.flatten(vars)), []; name, systems=mcons) end -@component function MaterialConnector(ms::MaterialSource;name) +@connector function MaterialConnector(ms::MaterialSource;name) vars = @variables begin - T(t), [description="temperature"] #, unit=u"K"] - ϱ(t), [description="density"] #, unit=u"mol m^-3"] + 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", connect=Stream] #, unit=u"J mol^-1"] - s(t), [description="molar entropy", connect=Stream] #, unit=u"J mol^-1 K^-1"] - xᵢ(t)[1:ms.N_c], [description="mole fractions", connect=Stream] #, unit=u"mol mol^-1"] + h(t), [description="molar enthalpy", output=true] #, unit=u"J mol^-1"] + s(t), [description="molar entropy", output=true] #, unit=u"J mol^-1 K^-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 diff --git a/src/base/utils.jl b/src/base/utils.jl index 70fda72..e480f01 100644 --- a/src/base/utils.jl +++ b/src/base/utils.jl @@ -27,15 +27,3 @@ return ODESystem(eqs, t, collect(Iterators.flatten(vars)), []; name, systems=[c]) end - -function mconnect(c1,c2) - Equation[ - c1.T ~ c2.T, - c1.ϱ ~ c2.ϱ, - c1.p ~ c2.p, - c1.h ~ c2.h, - c1.s ~ c2.s, - scalarize(c1.xᵢ .~ c2.xᵢ)..., - 0.0 ~ c1.n + c2.n, - ] -end \ No newline at end of file diff --git a/test/base/simple_steady_state.jl b/test/base/simple_steady_state.jl index 6074479..7ee645e 100644 --- a/test/base/simple_steady_state.jl +++ b/test/base/simple_steady_state.jl @@ -25,8 +25,8 @@ matsource = PS.MaterialSource("helium"; @named outlet = PS.Port(matsource) con_eqs_comp = [ - PS.mconnect(inlet.c, comp.cv.c1)..., - PS.mconnect(comp.cv.c2, outlet.c)... + connect(inlet.c, comp.cv.c1), + connect(comp.cv.c2, outlet.c) ] @named compressor_ = ODESystem(con_eqs_comp, t, [], []; systems=[inlet,comp,outlet]) @@ -66,18 +66,18 @@ systems = [comp_12,heat_22⁺,heat_2⁺3,turb_34,heat_44⁺,heat_4⁺1] streams = [inlet,s2,s2⁺,s3,s4,s4⁺,outlet] con_eqs = vcat( - PS.mconnect(inlet.c, comp_12.cv.c1)..., - PS.mconnect(comp_12.cv.c2, s2.c1)..., - PS.mconnect(s2.c2, heat_22⁺.cv.c1)..., - PS.mconnect(heat_22⁺.cv.c2, s2⁺.c1)..., - PS.mconnect(s2⁺.c2, heat_2⁺3.cv.c1)..., - PS.mconnect(heat_2⁺3.cv.c2, s3.c1)..., - PS.mconnect(s3.c2, turb_34.cv.c1)..., - PS.mconnect(turb_34.cv.c2, s4.c1)..., - PS.mconnect(s4.c2, heat_44⁺.cv.c1)..., - PS.mconnect(heat_44⁺.cv.c2, s4⁺.c1)..., - PS.mconnect(s4⁺.c2, heat_4⁺1.cv.c1)..., - PS.mconnect(heat_4⁺1.cv.c2, outlet.c)... + connect(inlet.c, comp_12.cv.c1), + connect(comp_12.cv.c2, s2.c1), + connect(s2.c2, heat_22⁺.cv.c1), + connect(heat_22⁺.cv.c2, s2⁺.c1), + connect(s2⁺.c2, heat_2⁺3.cv.c1), + connect(heat_2⁺3.cv.c2, s3.c1), + connect(s3.c2, turb_34.cv.c1), + connect(turb_34.cv.c2, s4.c1), + connect(s4.c2, heat_44⁺.cv.c1), + connect(heat_44⁺.cv.c2, s4⁺.c1), + connect(s4⁺.c2, heat_4⁺1.cv.c1), + connect(heat_4⁺1.cv.c2, outlet.c) ) append!(con_eqs,[0.0 ~ heat_2⁺3.Q + heat_4⁺1.Q]) From b88132a6137e9af31a01b8a9b3bbad9a85f14ebf Mon Sep 17 00:00:00 2001 From: se-schmitt Date: Tue, 22 Oct 2024 10:34:09 +0200 Subject: [PATCH 5/7] Implementation of CSTR (not validated yet) --- src/ProcessSimulator.jl | 3 ++ src/base/base_components.jl | 49 ++++++++++++---------- src/base/materials.jl | 24 +++++++---- src/base/utils.jl | 3 +- src/reactors/CSTR.jl | 32 +++++++++++++-- test/reactors/simple_cstr.jl | 79 ++++++++++++++++++++++++++++++++++++ 6 files changed, 155 insertions(+), 35 deletions(-) diff --git a/src/ProcessSimulator.jl b/src/ProcessSimulator.jl index 16658d8..2ae7646 100644 --- a/src/ProcessSimulator.jl +++ b/src/ProcessSimulator.jl @@ -13,4 +13,7 @@ include("base/utils.jl") include("fluid_handling/compressors.jl") include("fluid_handling/heat_exchangers.jl") +# Reactors +include("reactors/CSTR.jl") + end diff --git a/src/base/base_components.jl b/src/base/base_components.jl index e8d4c6a..49937e8 100644 --- a/src/base/base_components.jl +++ b/src/base/base_components.jl @@ -9,7 +9,6 @@ ϱ(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"] - s(t), [description="molar entropy"] #, unit=u"J mol^-1 K^-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 @@ -18,14 +17,12 @@ # EOS ϱ ~ ms.molar_density(p,T,xᵢ;phase=phase), h ~ ms.VT_enthalpy(ϱ,T,xᵢ), - s ~ ms.VT_entropy(ϱ,T,xᵢ), - 1.0 ~ sum([xᵢ for xᵢ in xᵢ]), + 1.0 ~ sum(collect(xᵢ)), # Connector "1" T ~ c1.T, ϱ ~ c1.ϱ, p ~ c1.p, h ~ c1.h, - s ~ c1.s, n ~ c1.n, scalarize(xᵢ .~ c1.xᵢ)..., # Connector "2" @@ -33,7 +30,6 @@ ϱ ~ c2.ϱ, p ~ c2.p, h ~ c2.h, - s ~ c2.s, n ~ -c2.n, scalarize(xᵢ .~ c2.xᵢ)..., ] @@ -47,7 +43,6 @@ end ϱ(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"] - s(t), [description="molar entropy", output=true] #, unit=u"J mol^-1 K^-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 @@ -102,39 +97,49 @@ end N_mcons, N_heats=0, N_works=0, - N_ph=1, phases=repeat(["unknown"],N_ph), + phases=["unknown"], 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"] - 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"] + 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"] end + !isempty(ms.reaction) ? 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"] + τ(t), [description="residence time (m/Δm_dot)"] #, unit=u"s"] + end) : nothing pars = @parameters begin V, [description="volume", bounds=(0,Inf)] #, unit=u"m^3"] end 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])) + Δ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 + ΔHᵣ + D(U) ~ ΔH + ΔE + (isempty(ms.reaction) ? 0.0 : ΔHᵣ), # Mole balance - [D(sum(nᵢ[:,j])) ~ sum([c.xᵢ[j]*c.n for c in mcons]) for j in 1:ms.N_c][:] - n ~ sum(nᵢ) + [D(sum(collect(nᵢ[:,i]))) ~ + sum([c.xᵢ[i]*c.n for c in mcons])*(isempty(ms.reaction) ? 1.0 : τ) + + (isempty(ms.reaction) ? 0.0 : ΔnR[i]) 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 - [ϱ[i] ~ ms.molar_density(p,T,nᵢ[i,:];phase=phases[i]) for i in 1:N_ph] - U ~ sum([ms.VT_internal_energy(ϱ[i],T,nᵢ[i,:]) for i in 1:N_ph]) + [ϱ[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(nᵢ[j,:])) for j in 1:N_ph]) ] return ODESystem(eqs, t, collect(Iterators.flatten(vars)), pars; name, systems=[mcons...,works...,heats...]) diff --git a/src/base/materials.jl b/src/base/materials.jl index c707401..9e0b603 100644 --- a/src/base/materials.jl +++ b/src/base/materials.jl @@ -1,16 +1,24 @@ 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...)) - molar_density::Function # Molar density function (defined as f(p,T,xᵢ;kwargs...)) - VT_internal_energy::Function # Internal energy function (defined as f(ϱ,T,xᵢ;kwargs...)) - VT_enthalpy::Function # Enthalpy function (defined as f(ϱ,T,xᵢ;kwargs...)) - VT_entropy::Function # Entropy function (defined as f(ϱ,T,xᵢ;kwargs...)) + 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...) @@ -36,6 +44,8 @@ function MaterialSource(components::Union{String,Vector{String}}; kwargs...) 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, :tp_flash, (a,T,n;kws...) -> f_NA(:tp_flash)), + get(kwargs, :reactions, Reaction[]) ) -end \ No newline at end of file +end + diff --git a/src/base/utils.jl b/src/base/utils.jl index e480f01..676335d 100644 --- a/src/base/utils.jl +++ b/src/base/utils.jl @@ -13,9 +13,8 @@ eqs = Equation[ # EOS ϱ ~ ms.molar_density(p,T,xᵢ;phase=phase), - 1.0 ~ sum([xᵢ for xᵢ in xᵢ]), + 1.0 ~ sum(collect(xᵢ)), c.h ~ ms.VT_enthalpy(ϱ,T,xᵢ), - c.s ~ ms.VT_entropy(ϱ,T,xᵢ), # Connector T ~ c.T, ϱ ~ c.ϱ, diff --git a/src/reactors/CSTR.jl b/src/reactors/CSTR.jl index 595763a..44f29f7 100644 --- a/src/reactors/CSTR.jl +++ b/src/reactors/CSTR.jl @@ -1,8 +1,32 @@ -@component function CSTR(ms:MaterialSource,reac::Reaction;name) +@component function CSTR(ms::MaterialSource;name, i_reacts=1:length(ms.reaction)) # Subsystems - @named cv = SimpleControlVolume(ms;) + @named cv = TPControlVolume(ms;N_mcons=2,N_heats=1,N_works=1,phases=["liquid"]) - # Stoichiometry + # Variables + vars = @variables begin + Q(t), [description="heat flux"] #, unit=u"J s^-1"] + end - return ODESystem() + # 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]]..., + # Change of enthalpy by reaction + [cv.ΔHᵣ ~ reac.r(cv.p,cv.T,collect(cv.xᵢ[1,:]))*reac.Δhᵣ(cv.T) for reac in ms.reaction[i_reacts]]..., + # Residence time (assuming ̇vᵢₙ = ̇vₒᵤₜ, V = const.) + cv.τ ~ cv.V / sum(collect(cv.c1.n .* cv.c1.xᵢ ./ cv.c1.ϱ)), + cv.c1.n ./ cv.c1.ϱ ~ cv.c2.n ./ cv.c2.ϱ + ] + + return ODESystem(eqs, t, vars, pars; name, systems=[cv]) end \ No newline at end of file diff --git a/test/reactors/simple_cstr.jl b/test/reactors/simple_cstr.jl index 62635af..1694eeb 100644 --- a/test/reactors/simple_cstr.jl +++ b/test/reactors/simple_cstr.jl @@ -1,3 +1,82 @@ using ProcessSimulator using ModelingToolkit, DifferentialEquations using ModelingToolkit: t_nounits as t + +const PS = ProcessSimulator + +# Material Parameters +ϱ = [14.8,55.3,13.7,24.7].*1e3 # mol/m^3 +cₚ = [146.440,75.312,192.464,81.588] # J/mol/K + +# Reaction +ν = [-1.0, -1.0, 1.0, 0.0] +Δhᵣ = -83.734392e6 # J/mol + +matsource = PS.MaterialSource(["propylene oxide (A)","water (B)","propylene glycol (C)","methanol (M)"]; + Mw = [0.05808, 0.018015, 0.07609, 0.03204], + molar_density = (p, T, x; kwargs...) -> sum([ϱ[i]*x[i] for i in 1:4]), + VT_enthalpy = (ϱ, T, x) -> sum([cₚ[i]*x[i] for i in 1:4])*T, + VT_internal_energy = (ϱ, T, x) -> sum([cₚ[i]*x[i] for i in 1:4])*T, + reactions = [PS.Reaction( + ν = ν, + r = (p,T,x) -> -2.73e-4*exp(9059*(1/297-1/T))*x[1], + Δhᵣ = x -> Δhᵣ, + )], +) + +F = [36.3e3, 453.6e3, 0, 45.4e3]./3600 # mol/s +TF = 297.04 # K + +UA = 8440.06 # J/K/s +T_coolant = 288.71 # K +cₚ_coolant = 75.312 # J/mol/K + + +# Create flowsheet +@named inlet = PS.Port(matsource) +@named cstr = PS.CSTR(matsource) +@named outlet = PS.Port(matsource) + +# Connect the flowsheet +eqs = [ + connect(inlet.c,cstr.cv.c1), + connect(cstr.cv.c2,outlet.c), +] + +@named flowsheet_ = ODESystem(eqs, t, [], [], systems=[inlet,cstr,outlet]) + +pars = [ + cstr.cv.V => 1.89, +] + +inp = [ + inlet.T => 297.04, + inlet.p => 1e5, + inlet.n => sum(F), + inlet.xᵢ[1] => F[1]/sum(F), + inlet.xᵢ[2] => F[2]/sum(F), + inlet.xᵢ[3] => F[3]/sum(F), + cstr.cv.q1.Q => 0.0, + outlet.p => 1e5, +] +out = [ + cstr.cv.T, + cstr.cv.xᵢ[1,1], + cstr.cv.xᵢ[1,2], + cstr.cv.xᵢ[1,3], +] + +u0 = [ + cstr.cv.nᵢ[1,1] => 0.0, + cstr.cv.nᵢ[1,2] => 55.3e3*Dict(pars...)[cstr.cv.V], + cstr.cv.nᵢ[1,3] => 0.0, + cstr.cv.nᵢ[1,4] => 0.0, + cstr.cv.T => 297., + cstr.cv.ΔnR[1] => 0.0, + cstr.cv.ΔnR[2] => 0.0, + cstr.cv.ΔnR[3] => 0.0, +] + +flowsheet,idx = structural_simplify(flowsheet_,(first.(inp),out)) + +prob = ODEProblem(flowsheet, u0, (0.0, 1.4e4), collect(Iterators.flatten([pars,inp]))) \ No newline at end of file From fe5e3e6d704e504cdea5954c74feb6c2da07aa0c Mon Sep 17 00:00:00 2001 From: se-schmitt Date: Thu, 24 Oct 2024 13:41:12 +0200 Subject: [PATCH 6/7] CSTR implementation working --- src/base/base_components.jl | 13 ++-- src/reactors/CSTR.jl | 13 ++-- test/reactors/simple_cstr.jl | 127 ++++++++++++++++++++++++++++------- 3 files changed, 115 insertions(+), 38 deletions(-) diff --git a/src/base/base_components.jl b/src/base/base_components.jl index 49937e8..376d561 100644 --- a/src/base/base_components.jl +++ b/src/base/base_components.jl @@ -98,6 +98,7 @@ end N_heats=0, N_works=0, phases=["unknown"], + reactive=false, name) # Init N_ph = length(phases) @@ -116,10 +117,9 @@ end ΔH(t), [description="enthalpy difference inlets/outlets"] #, unit=u"J/s"] ΔE(t), [description="added/removed heat or work"] #, unit=u"J/s"] end - !isempty(ms.reaction) ? append!(vars, @variables begin + 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"] - τ(t), [description="residence time (m/Δm_dot)"] #, unit=u"s"] end) : nothing pars = @parameters begin @@ -130,16 +130,15 @@ end Δ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 + (isempty(ms.reaction) ? 0.0 : ΔHᵣ), + 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])*(isempty(ms.reaction) ? 1.0 : τ) + - (isempty(ms.reaction) ? 0.0 : ΔnR[i]) for i in 1:ms.N_c]..., + [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(nᵢ[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]) ] return ODESystem(eqs, t, collect(Iterators.flatten(vars)), pars; name, systems=[mcons...,works...,heats...]) diff --git a/src/reactors/CSTR.jl b/src/reactors/CSTR.jl index 44f29f7..e17f29f 100644 --- a/src/reactors/CSTR.jl +++ b/src/reactors/CSTR.jl @@ -1,6 +1,6 @@ @component function CSTR(ms::MaterialSource;name, i_reacts=1:length(ms.reaction)) # Subsystems - @named cv = TPControlVolume(ms;N_mcons=2,N_heats=1,N_works=1,phases=["liquid"]) + @named cv = TPControlVolume(ms;N_mcons=2,N_heats=1,N_works=1,phases=["liquid"],reactive=true) # Variables vars = @variables begin @@ -21,11 +21,12 @@ 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]]..., - # Change of enthalpy by reaction - [cv.ΔHᵣ ~ reac.r(cv.p,cv.T,collect(cv.xᵢ[1,:]))*reac.Δhᵣ(cv.T) for reac in ms.reaction[i_reacts]]..., - # Residence time (assuming ̇vᵢₙ = ̇vₒᵤₜ, V = const.) - cv.τ ~ cv.V / sum(collect(cv.c1.n .* cv.c1.xᵢ ./ cv.c1.ϱ)), - cv.c1.n ./ cv.c1.ϱ ~ cv.c2.n ./ cv.c2.ϱ + # 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]]..., + # Constant volume/mass + # 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), # const. masss + cv.V ~ cv.n / cv.ϱ[1] + cv.c1.n / cv.c1.ϱ + cv.c2.n / cv.c2.ϱ, # const. volume + # cv.c2.n ~ cv.c1.n/cv.c1.ϱ / cv.V * cv.n ] return ODESystem(eqs, t, vars, pars; name, systems=[cv]) diff --git a/test/reactors/simple_cstr.jl b/test/reactors/simple_cstr.jl index 1694eeb..2211db4 100644 --- a/test/reactors/simple_cstr.jl +++ b/test/reactors/simple_cstr.jl @@ -5,32 +5,39 @@ using ModelingToolkit: t_nounits as t const PS = ProcessSimulator # Material Parameters -ϱ = [14.8,55.3,13.7,24.7].*1e3 # mol/m^3 -cₚ = [146.440,75.312,192.464,81.588] # J/mol/K +ϱ = [14.8,55.3,13.7,24.7].*1e3 # mol/m^3 +cₚ = [35,18,46,19.5]*4.184 # J/mol/K # Reaction ν = [-1.0, -1.0, 1.0, 0.0] -Δhᵣ = -83.734392e6 # J/mol +Δhᵣ = 20013*4.184 # J/mol matsource = PS.MaterialSource(["propylene oxide (A)","water (B)","propylene glycol (C)","methanol (M)"]; - Mw = [0.05808, 0.018015, 0.07609, 0.03204], + Mw = [0.05808, 0.01801, 0.07609, 0.03204], molar_density = (p, T, x; kwargs...) -> sum([ϱ[i]*x[i] for i in 1:4]), VT_enthalpy = (ϱ, T, x) -> sum([cₚ[i]*x[i] for i in 1:4])*T, VT_internal_energy = (ϱ, T, x) -> sum([cₚ[i]*x[i] for i in 1:4])*T, reactions = [PS.Reaction( ν = ν, - r = (p,T,x) -> -2.73e-4*exp(9059*(1/297-1/T))*x[1], - Δhᵣ = x -> Δhᵣ, + r = (p,T,x) -> 2.73e-4*exp(9059*(1/297-1/T))*x[1], + Δhᵣ = (T) -> Δhᵣ, )], ) +# Parameters and boundary conditions +V_cstr = 1.89 # m^3 F = [36.3e3, 453.6e3, 0, 45.4e3]./3600 # mol/s -TF = 297.04 # K -UA = 8440.06 # J/K/s -T_coolant = 288.71 # K -cₚ_coolant = 75.312 # J/mol/K +# Initial conditions +T0 = 297. +n0 = 55.3e3*V_cstr +m0_inlet = F' * matsource.Mw +# Heat exchanger +UA = 7262*4184/3600 # J/K/s +T1_coolant = 288.71 # K +cₚ_coolant = 18*4.184 # J/mol/K +n_coolant = 453.6e3/3600 # mol/s # Create flowsheet @named inlet = PS.Port(matsource) @@ -41,42 +48,112 @@ cₚ_coolant = 75.312 # J/mol/K eqs = [ connect(inlet.c,cstr.cv.c1), connect(cstr.cv.c2,outlet.c), + cstr.Q ~ -n_coolant*cₚ_coolant*(cstr.cv.T-T1_coolant)*(1-exp(-UA/(n_coolant*cₚ_coolant))), ] @named flowsheet_ = ODESystem(eqs, t, [], [], systems=[inlet,cstr,outlet]) pars = [ - cstr.cv.V => 1.89, + cstr.cv.V => V_cstr, ] inp = [ - inlet.T => 297.04, + inlet.T => 297.0, inlet.p => 1e5, inlet.n => sum(F), inlet.xᵢ[1] => F[1]/sum(F), inlet.xᵢ[2] => F[2]/sum(F), inlet.xᵢ[3] => F[3]/sum(F), - cstr.cv.q1.Q => 0.0, outlet.p => 1e5, ] -out = [ - cstr.cv.T, - cstr.cv.xᵢ[1,1], - cstr.cv.xᵢ[1,2], - cstr.cv.xᵢ[1,3], -] u0 = [ + cstr.cv.U => matsource.VT_internal_energy(NaN, T0, [0,1,0,0])*n0, cstr.cv.nᵢ[1,1] => 0.0, - cstr.cv.nᵢ[1,2] => 55.3e3*Dict(pars...)[cstr.cv.V], + cstr.cv.nᵢ[1,2] => n0, cstr.cv.nᵢ[1,3] => 0.0, cstr.cv.nᵢ[1,4] => 0.0, + inlet.m => m0_inlet, cstr.cv.T => 297., - cstr.cv.ΔnR[1] => 0.0, - cstr.cv.ΔnR[2] => 0.0, - cstr.cv.ΔnR[3] => 0.0, + cstr.cv.c2.n => 0.0, + outlet.m => 0., ] -flowsheet,idx = structural_simplify(flowsheet_,(first.(inp),out)) +# u0 = [ +# cstr.cv.U => matsource.VT_internal_energy(NaN, T0, [0,1,0,0])*n0, +# cstr.cv.nᵢ[1,1] => 0.0, +# cstr.cv.nᵢ[1,2] => n0, +# cstr.cv.nᵢ[1,3] => 0.0, +# cstr.cv.nᵢ[1,4] => 0.0, +# cstr.cv.T => 297., +# inlet.m => m0_inlet, +# outlet.m => 0., +# ] + +flowsheet,idx = structural_simplify(flowsheet_,(first.(inp),[])) + +# prob_steady = SteadyStateProblem(flowsheet, vcat(pars,inp), u0) +# sol_steady = solve(prob_steady) + +prob = ODEProblem(flowsheet, u0, (0, 4).*3600., vcat(pars,inp))#; guesses = [outlet.m => -m0_inlet]) +@time sol = solve(prob, QNDF(), abstol = 1e-6, reltol = 1e-6) + + +# -------------------------- Plotting and evaluation ------------------------- # +using Plots, PrettyTables +pythonplot() +# gr() +# plotlyjs() +closeall() + +t_h = sol.t./3600 + +ps = [plot(;framestyle=:box,xlabel="t / h",xlims=(0.0,maximum(sol.t)*1.02/3600)) for i in 1:6] +plot!(ps[1],t_h,sol[cstr.cv.T],label="T",ylabel="T / K") +[plot!(ps[2],t_h, sol[cstr.cv.nᵢ[1,i]]./V_cstr./1e3;label="x[$i]",ylabel="Cᵢ / kmol/m³") for i in 1:4] +plot!(ps[3],t_h,sol[cstr.cv.n]./1e3,label="n",xlabel="t / h",ylabel="n / kmol") +plot!(ps[4],t_h,matsource.reaction[1].r.(NaN,sol[cstr.cv.T],sol[cstr.cv.xᵢ]),label="r",xlabel="t / h",ylabel="r / kmol") +plot!(ps[5],t_h,sol[cstr.cv.ΔnR[3]]./1e3,label="Q",xlabel="t / h",ylabel="ΔnR / kmol") +plot!(ps[6],t_h,sol[cstr.cv.c1.h],label="h_in",xlabel="t / h",ylabel="ΔE / Js") +plot!(ps[6],t_h,sol[cstr.cv.c2.h],label="h_out") +plot!(ps[6],t_h,sol[cstr.cv.ΔHᵣ],label="ΔHᵣ") +plot!(ps[6],t_h,abs.(sol[cstr.cv.q1.Q]),label="|Q|") +fig = plot(ps...;layout=(2,3),size=(1200,800)) +display(fig) + +# Plot Fogler (Figures E13-3.1+3.2) +plot() +plts_Fogler = [plot(;framestyle=:box,xlabel="t / h",xlims=(0,4)) for i in 1:2] +plot!(plts_Fogler[1],t_h, sol[cstr.cv.nᵢ[1,1]]./V_cstr./1e3; ylabel="Cₐ / kmol/m³") +plot!(plts_Fogler[2],t_h, sol[cstr.cv.T]; ylabel="T / K") +fig_Fogler = plot(plts_Fogler...;layout=(1,2)) +display(fig_Fogler) + +# Plots masss and volume +V_sim = sol[cstr.cv.n] ./ sol[cstr.cv.ϱ[1]] +m_sim = [(ni * matsource.Mw)[1] for ni in sol[cstr.cv.nᵢ]] + +fig_Vm = plot( + plot(t_h,V_sim./V_sim[1];xlabel="t / h",ylabel="V / m³",framestyle=:box), + plot(t_h,m_sim;xlabel="t / h",ylabel="m / kg",framestyle=:box); + layout=(2,1) +) +display(fig_Vm) + +# Table +nᵢ_sol = [sol[cstr.cv.nᵢ[1,i]] for i in 1:4] +header =[ + "Variable", "Initial Fogler", "Initial Sim.", "Final Fogler", "Final Sim." +] +data = [ + "Ca" 0.0 nᵢ_sol[1][1]/V_cstr/1e3 0.658258 nᵢ_sol[1][end]/V_cstr/1e3; + "Cb" 55.3 nᵢ_sol[2][1]/V_cstr/1e3 34.06019 nᵢ_sol[2][end]/V_cstr/1e3; + "Cc" 0.0 nᵢ_sol[3][1]/V_cstr/1e3 2.247301 nᵢ_sol[3][end]/V_cstr/1e3; + "Nb" 104.517 nᵢ_sol[2][1]/1e3 64.37375 nᵢ_sol[2][end]/1e3; + "Nc" 0.0 nᵢ_sol[3][1]/1e3 4.2474 nᵢ_sol[3][end]/1e3; + "Nm" 0.0 nᵢ_sol[4][1]/1e3 6.868166 nᵢ_sol[4][end]/1e3; + "Qr2" 39920*4184/3600 sol[cstr.cv.q1.Q][1] 205900 sol[cstr.cv.q1.Q][end]; + "T" 297.0 sol[cstr.cv.T][1] 331.4976 sol[cstr.cv.T][end]; +] -prob = ODEProblem(flowsheet, u0, (0.0, 1.4e4), collect(Iterators.flatten([pars,inp]))) \ No newline at end of file +pretty_table(data; header=header) \ No newline at end of file From b3093e624466caec63eb19c7a276ea5e5e1da5b1 Mon Sep 17 00:00:00 2001 From: se-schmitt Date: Thu, 24 Oct 2024 15:44:39 +0200 Subject: [PATCH 7/7] Changed volume to variable, added `flowtype` kw and test for CSTR --- src/base/base_components.jl | 10 ++-- src/reactors/CSTR.jl | 11 +++-- test/reactors/simple_cstr.jl | 95 +++++++----------------------------- test/runtests.jl | 6 ++- 4 files changed, 32 insertions(+), 90 deletions(-) diff --git a/src/base/base_components.jl b/src/base/base_components.jl index 376d561..095509c 100644 --- a/src/base/base_components.jl +++ b/src/base/base_components.jl @@ -116,16 +116,13 @@ end 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 - pars = @parameters begin - V, [description="volume", bounds=(0,Inf)] #, unit=u"m^3"] - end - 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])), @@ -138,8 +135,9 @@ end [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]) + 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)), pars; name, systems=[mcons...,works...,heats...]) + return ODESystem(eqs, t, collect(Iterators.flatten(vars)), []; name, systems=[mcons...,works...,heats...]) end \ No newline at end of file diff --git a/src/reactors/CSTR.jl b/src/reactors/CSTR.jl index e17f29f..dba50af 100644 --- a/src/reactors/CSTR.jl +++ b/src/reactors/CSTR.jl @@ -1,4 +1,4 @@ -@component function CSTR(ms::MaterialSource;name, i_reacts=1:length(ms.reaction)) +@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) @@ -23,11 +23,12 @@ [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]]..., - # Constant volume/mass - # 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), # const. masss - cv.V ~ cv.n / cv.ϱ[1] + cv.c1.n / cv.c1.ϱ + cv.c2.n / cv.c2.ϱ, # const. volume - # cv.c2.n ~ cv.c1.n/cv.c1.ϱ / cv.V * cv.n ] + 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 \ No newline at end of file diff --git a/test/reactors/simple_cstr.jl b/test/reactors/simple_cstr.jl index 2211db4..4f5cb5e 100644 --- a/test/reactors/simple_cstr.jl +++ b/test/reactors/simple_cstr.jl @@ -41,7 +41,7 @@ n_coolant = 453.6e3/3600 # mol/s # Create flowsheet @named inlet = PS.Port(matsource) -@named cstr = PS.CSTR(matsource) +@named cstr = PS.CSTR(matsource;flowtype="const. mass") @named outlet = PS.Port(matsource) # Connect the flowsheet @@ -53,9 +53,7 @@ eqs = [ @named flowsheet_ = ODESystem(eqs, t, [], [], systems=[inlet,cstr,outlet]) -pars = [ - cstr.cv.V => V_cstr, -] +pars = [] inp = [ inlet.T => 297.0, @@ -79,81 +77,22 @@ u0 = [ outlet.m => 0., ] -# u0 = [ -# cstr.cv.U => matsource.VT_internal_energy(NaN, T0, [0,1,0,0])*n0, -# cstr.cv.nᵢ[1,1] => 0.0, -# cstr.cv.nᵢ[1,2] => n0, -# cstr.cv.nᵢ[1,3] => 0.0, -# cstr.cv.nᵢ[1,4] => 0.0, -# cstr.cv.T => 297., -# inlet.m => m0_inlet, -# outlet.m => 0., -# ] - flowsheet,idx = structural_simplify(flowsheet_,(first.(inp),[])) -# prob_steady = SteadyStateProblem(flowsheet, vcat(pars,inp), u0) -# sol_steady = solve(prob_steady) - -prob = ODEProblem(flowsheet, u0, (0, 4).*3600., vcat(pars,inp))#; guesses = [outlet.m => -m0_inlet]) -@time sol = solve(prob, QNDF(), abstol = 1e-6, reltol = 1e-6) - - -# -------------------------- Plotting and evaluation ------------------------- # -using Plots, PrettyTables -pythonplot() -# gr() -# plotlyjs() -closeall() - -t_h = sol.t./3600 - -ps = [plot(;framestyle=:box,xlabel="t / h",xlims=(0.0,maximum(sol.t)*1.02/3600)) for i in 1:6] -plot!(ps[1],t_h,sol[cstr.cv.T],label="T",ylabel="T / K") -[plot!(ps[2],t_h, sol[cstr.cv.nᵢ[1,i]]./V_cstr./1e3;label="x[$i]",ylabel="Cᵢ / kmol/m³") for i in 1:4] -plot!(ps[3],t_h,sol[cstr.cv.n]./1e3,label="n",xlabel="t / h",ylabel="n / kmol") -plot!(ps[4],t_h,matsource.reaction[1].r.(NaN,sol[cstr.cv.T],sol[cstr.cv.xᵢ]),label="r",xlabel="t / h",ylabel="r / kmol") -plot!(ps[5],t_h,sol[cstr.cv.ΔnR[3]]./1e3,label="Q",xlabel="t / h",ylabel="ΔnR / kmol") -plot!(ps[6],t_h,sol[cstr.cv.c1.h],label="h_in",xlabel="t / h",ylabel="ΔE / Js") -plot!(ps[6],t_h,sol[cstr.cv.c2.h],label="h_out") -plot!(ps[6],t_h,sol[cstr.cv.ΔHᵣ],label="ΔHᵣ") -plot!(ps[6],t_h,abs.(sol[cstr.cv.q1.Q]),label="|Q|") -fig = plot(ps...;layout=(2,3),size=(1200,800)) -display(fig) - -# Plot Fogler (Figures E13-3.1+3.2) -plot() -plts_Fogler = [plot(;framestyle=:box,xlabel="t / h",xlims=(0,4)) for i in 1:2] -plot!(plts_Fogler[1],t_h, sol[cstr.cv.nᵢ[1,1]]./V_cstr./1e3; ylabel="Cₐ / kmol/m³") -plot!(plts_Fogler[2],t_h, sol[cstr.cv.T]; ylabel="T / K") -fig_Fogler = plot(plts_Fogler...;layout=(1,2)) -display(fig_Fogler) - -# Plots masss and volume -V_sim = sol[cstr.cv.n] ./ sol[cstr.cv.ϱ[1]] -m_sim = [(ni * matsource.Mw)[1] for ni in sol[cstr.cv.nᵢ]] - -fig_Vm = plot( - plot(t_h,V_sim./V_sim[1];xlabel="t / h",ylabel="V / m³",framestyle=:box), - plot(t_h,m_sim;xlabel="t / h",ylabel="m / kg",framestyle=:box); - layout=(2,1) -) -display(fig_Vm) +prob = ODEProblem(flowsheet, u0, (0, 2).*3600., vcat(inp))#; guesses = [outlet.m => -m0_inlet]) +sol = solve(prob, QNDF(), abstol = 1e-6, reltol = 1e-6) -# Table -nᵢ_sol = [sol[cstr.cv.nᵢ[1,i]] for i in 1:4] -header =[ - "Variable", "Initial Fogler", "Initial Sim.", "Final Fogler", "Final Sim." -] -data = [ - "Ca" 0.0 nᵢ_sol[1][1]/V_cstr/1e3 0.658258 nᵢ_sol[1][end]/V_cstr/1e3; - "Cb" 55.3 nᵢ_sol[2][1]/V_cstr/1e3 34.06019 nᵢ_sol[2][end]/V_cstr/1e3; - "Cc" 0.0 nᵢ_sol[3][1]/V_cstr/1e3 2.247301 nᵢ_sol[3][end]/V_cstr/1e3; - "Nb" 104.517 nᵢ_sol[2][1]/1e3 64.37375 nᵢ_sol[2][end]/1e3; - "Nc" 0.0 nᵢ_sol[3][1]/1e3 4.2474 nᵢ_sol[3][end]/1e3; - "Nm" 0.0 nᵢ_sol[4][1]/1e3 6.868166 nᵢ_sol[4][end]/1e3; - "Qr2" 39920*4184/3600 sol[cstr.cv.q1.Q][1] 205900 sol[cstr.cv.q1.Q][end]; - "T" 297.0 sol[cstr.cv.T][1] 331.4976 sol[cstr.cv.T][end]; -] +(Tmax, iTmax) = findmax(sol[cstr.cv.T]) + +@test Tmax ≈ 356.149 atol=1e-3 +@test sol.t[iTmax] ≈ 2823.24 atol=1e-2 + +if isinteractive() + # Plots + using Plots -pretty_table(data; header=header) \ No newline at end of file + ps = [plot(;framestyle=:box,xlabel="t / h",xlims=(0.0,2.0)) for i in 1:2] + plot!(ps[1],sol.t/3600,sol[cstr.cv.T],label="T",ylabel="T / K") + [plot!(ps[2],sol.t/3600, sol[cstr.cv.xᵢ[1,i]];label=matsource.components[i],ylabel="xᵢ / mol/mol") for i in 1:4] + plot(ps...;layout=(1,2),size=(800,400)) +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index cc92390..cb00bab 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,6 +2,10 @@ using ProcessSimulator using Test using SafeTestsets -@safetestset "Simple Steady State" begin +@safetestset "Base components" begin include("base/simple_steady_state.jl") end + +@safetestset "Reactors" begin + include("reactors/simple_cstr.jl") +end