examples for reactor and for steady state analysis

examples/example_reactor.py

@@ -0,0 +1,94 @@
+#########################################################################################
+##
+##                             PathSim Chemical Reactor Example
+##
+#########################################################################################
+
+# IMPORTS ===============================================================================
+import numpy as np
+import matplotlib.pyplot as plt
+
+from pathsim import Simulation, Connection
+from pathsim.blocks import ODE, Source, Scope
+from pathsim.solvers import ESDIRK32, ESDIRK43, GEAR52A
+
+# CSTR WITH CONSECUTIVE REACTIONS INITIAL VALUE PROBLEM =================================
+
+# Initial conditions
+Ca_0 = 1.0    # Initial concentration of A
+Cb_0 = 0.0    # Initial concentration of B
+T_0 = 300.0   # Initial temperature
+
+
+# System parameters
+Tc = 280.0    # Coolant temperature
+tau = 1.0     # Residence time
+k1_0 = 1e4    # Rate constant 1
+k2_0 = 1e3    # Rate constant 2
+E1 = 5e4      # Activation energy 1
+E2 = 5.5e4    # Activation energy 2
+dH1 = -5e4    # Reaction enthalpy 1
+dH2 = -5.2e4  # Reaction enthalpy 2
+rho = 1000.0  # Density
+Cp = 4.184    # Heat capacity
+U = 1000.0    # Heat transfer coefficient
+V = 0.1       # Reactor volume
+A = 0.1       # Heat transfer area
+R = 8.314     # Gas constant
+
+
+# Define system blocks
+Sco    = Scope(labels=['Ca', 'Cb', 'T'])
+Src_Ca = Source(lambda t: 2.0 + np.sin(0.5*t))
+Src_T  = Source(lambda t: 280.0 * (1 - 0.8 * np.exp(-0.6*t)))
+
+def reaction_rates(x, u, t):
+
+    #unpack states
+    Ca, Cb, T = x
+
+    #unpack inputs
+    Ca_in, T_in = u
+
+    # Concentration dynamics
+    dCa_dt = (Ca_in - Ca)/tau - k1_0*np.exp(-E1/(R*T))*Ca
+    dCb_dt = -Cb/tau + k1_0*np.exp(-E1/(R*T))*Ca - k2_0*np.exp(-E2/(R*T))*Cb
+
+    # Temperature dynamics
+    dT_dt = (T_in - T)/tau + \
+            (-dH1/(rho*Cp))*k1_0*np.exp(-E1/(R*T))*Ca + \
+            (-dH2/(rho*Cp))*k2_0*np.exp(-E2/(R*T))*Cb - \
+            U*A*(T-Tc)/(V*rho*Cp)
+
+    return np.array([dCa_dt, dCb_dt, dT_dt])
+
+CSTR = ODE(reaction_rates, np.array([Ca_0, Cb_0, T_0]))
+
+# Main system blocks and connections
+blocks = [CSTR, Src_Ca, Src_T, Sco]
+connections = [
+    Connection(CSTR, Sco),         # Ca output
+    Connection(CSTR[1], Sco[1]),   # Cb output
+    Connection(CSTR[2], Sco[2]),   # T output
+    Connection(Src_Ca, CSTR[0]),
+    Connection(Src_T, CSTR[1])  
+]
+
+# Initialize simulation
+Sim = Simulation(
+    blocks,
+    connections,
+    dt=0.001,
+    log=True,
+    Solver=GEAR52A,
+    tolerance_lte_abs=1e-6,
+    tolerance_lte_rel=1e-4
+)
+
+# Run simulation for 10 seconds
+Sim.run(20, False)  
+
+# Plot results
+Sco.plot(".-", lw=1.5)
+
+plt.show()

examples/example_steadystate.py

@@ -0,0 +1,64 @@
+#########################################################################################
+##
+##           example of a simple feedback system with steady state analysis
+##
+#########################################################################################
+
+# IMPORTS ===============================================================================
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+from pathsim import Simulation, Connection
+
+from pathsim.blocks import (
+    Source, 
+    Integrator, 
+    Amplifier, 
+    Adder, 
+    Scope
+    )
+
+from pathsim.solvers import RKBS32
+
+
+# 1st ORDER SYSTEM ======================================================================
+
+#simulation timestep
+dt = 0.01
+
+#step delay
+tau = 1.5
+
+#blocks that define the system
+Src = Source(lambda t: int(t>tau))
+Int = Integrator()
+Amp = Amplifier(-0.5)
+Add = Adder()
+Sco = Scope(labels=["step", "response"])
+
+blocks = [Src, Int, Amp, Add, Sco]
+
+#the connections between the blocks
+connections = [
+    Connection(Src, Add[0], Sco[0]),
+    Connection(Amp, Add[1]),
+    Connection(Add, Int),
+    Connection(Int, Amp, Sco[1])
+    ]
+
+#initialize simulation with the blocks, connections, timestep and logging enabled
+Sim = Simulation(blocks, connections, dt=dt, log=True)
+
+#run the simulation for some time
+Sim.run(2*tau)
+
+#then force to steady state
+Sim.steadystate(reset=False)
+
+#then run some more 
+Sim.run(tau, reset=False)
+
+Sco.plot()
+
+plt.show()