systems.builtin.deterministic.continuous.ContinuousCSTR

systems.builtin.deterministic.continuous.ContinuousCSTR(*args, **kwargs)

Continuous-time Continuous Stirred-Tank Reactor (CSTR) with cooling jacket.

Physical System:

A continuous flow reactor where reactant A converts to product B in an exothermic reaction. The CSTR is one of the most studied nonlinear systems in chemical engineering, exhibiting rich dynamics including: - Multiple steady states (multiplicity) - Sustained oscillations (limit cycles) - Bifurcations and hysteresis - Thermal runaway behavior

Unlike batch reactors, CSTRs operate continuously with: - Continuous feed stream entering at Cₐ,feed, T_feed - Continuous product stream leaving at reactor conditions - Perfect mixing assumption (uniform concentration and temperature) - External cooling/heating via jacket

The CSTR represents a fundamental model in: - Process control (benchmark nonlinear system) - Nonlinear dynamics (canonical example of multiplicity) - Chemical reaction engineering (industrial reactor design) - Bifurcation theory (illustrates saddle-node, Hopf bifurcations)

State Space:

State: x = [Cₐ, T] Concentration state: - Cₐ: Concentration of reactant A in reactor [mol/L] * Governed by material balance: in - out - reaction * 0 ≤ Cₐ ≤ Cₐ,feed (bounded by feed concentration) * Low Cₐ → high conversion (desired but challenging to control) * High Cₐ → low conversion (safe but inefficient)

Temperature state:
- T: Reactor temperature [K]
  * Governed by energy balance: in - out + generation - removal
  * Typically T > T_feed for exothermic reactions
  * Exhibits strong nonlinear coupling with concentration
  * Critical for safety (runaway prevention)
  * Small changes can cause large rate changes (Arrhenius)

Control: u = [T_jacket] - T_jacket: Cooling jacket temperature [K] * Primary manipulated variable for temperature control * Affects heat removal rate via UA·(T - T_jacket) * Typically T_jacket < T (cooling mode) * Can be T_jacket > T (heating mode for startup/cold days) * Physical constraints: chiller/heater capacity limits * Rate constraints: jacket dynamics (not modeled here)

Output: y = [Cₐ, T] - Full state measurement (common in modern plants) - In practice: * Cₐ: Online analyzer (GC, HPLC, NIR spectroscopy) * T: Thermocouple or RTD (fast, reliable) * Both have measurement noise and potential delays

Dynamics:

The continuous-time CSTR dynamics are:

dCₐ/dt = (F/V)·(Cₐ,feed - Cₐ) - r
dT/dt = (F/V)·(T_feed - T) + (-ΔH/ρCₚ)·r + UA/(VρCₚ)·(T_jacket - T)

Reaction Rate (Arrhenius kinetics): r = k₀·Cₐ·exp(-E/T) [mol/(L·s)]

where: - k₀: Pre-exponential factor [1/s] * Collision frequency in Arrhenius equation * Typical range: 10⁶-10¹² for liquid-phase reactions * Material and reaction-specific constant - E: Activation energy [K] (dimensionless Eₐ/R) * Energy barrier for reaction to occur * Typical range: 5,000-15,000 K for Eₐ/R * Higher E → more temperature-sensitive reaction * Physical activation energy Eₐ typically 40-120 kJ/mol - exp(-E/T): Arrhenius temperature dependence * Exponential sensitivity creates strong nonlinearity * 10°C change can double/triple reaction rate * Source of multiple steady states and instability

Physical Interpretation of Each Term:

Material Balance (dCₐ/dt): 1. (F/V)·(Cₐ,feed - Cₐ): Convective in/out - F/V = 1/τ: Inverse residence time [1/s] - τ = V/F: Average time molecule spends in reactor [s] - Positive when Cₐ < Cₐ,feed (dilution effect) - Acts as “restoring force” toward feed concentration - Time scale: τ (seconds to minutes)

-r: Consumption by reaction
- Always negative (reactant consumed)
- Exponentially increases with temperature
- Depends on current concentration (first-order)
- Time scale: 1/k (fast at high T, slow at low T)

At steady state: inflow rate = outflow rate + reaction rate

Energy Balance (dT/dt): 1. (F/V)·(T_feed - T): Convective heat in/out - Negative when T > T_feed (typical for exothermic) - Same time scale as material balance (1/τ) - Acts as “restoring force” toward feed temperature

(-ΔH/ρCₚ)·r: Heat generation from reaction
- Positive for exothermic reaction (ΔH < 0)
- Couples concentration to temperature
- Creates positive feedback: higher T → faster r → more heat
- This term causes thermal runaway if unchecked
- Magnitude: |ΔH|/(ρCₚ) is adiabatic temperature rise per mol/L
UA/(VρCₚ)·(T_jacket - T): Heat removal via jacket
- Negative when T > T_jacket (cooling)
- Only term controlled by manipulated variable
- UA: Overall heat transfer coefficient × area [J/(s·K)]
- Larger UA → better temperature control
- Time scale: VρCₚ/UA (thermal time constant)

At steady state: heat in + heat generated = heat out + heat removed

Nonlinear Coupling and Feedback:

The CSTR exhibits strong positive feedback that can lead to instability:

Thermal Feedback Loop (Runaway Mechanism): T ↑ → r ↑ (Arrhenius) → heat generation ↑ → T ↑ (positive feedback)

This loop is stabilized by:
- Convective cooling: higher T → more heat removal to feed
- Jacket cooling: higher T → more heat removal via jacket
- Reactant depletion: higher r → lower Cₐ → lower r (negative feedback)
Material-Thermal Coupling:
- High T → fast reaction → low Cₐ (depletion)
- Low Cₐ → slow reaction → less heat generation → lower T
- This coupling creates multiple possible steady states
Competition Between Time Scales:
- Residence time τ = V/F (convective transport)
- Reaction time 1/k (chemical kinetics)
- Thermal time VρCₚ/UA (heat transfer)
- Relative magnitudes determine stability and multiplicity

Parameters:

F : float, default=100.0 Volumetric flow rate [L/s] - Controls residence time τ = V/F - Higher F → shorter τ → lower conversion but more stable - Lower F → longer τ → higher conversion but less stable - Typical range: 10-1000 L/s depending on reactor size - Often kept constant (set by upstream/downstream constraints)

V : float, default=100.0 Reactor volume [L] - Combined with F to give residence time τ - Larger V → more material holdup → slower dynamics - Typical range: 100-10,000 L for industrial reactors - Design parameter (fixed once reactor is built)

C_A_feed : float, default=1.0 Feed concentration [mol/L] - Upper bound for reactor concentration - Higher feed → more product but more heat generation - Typical range: 0.1-10 mol/L - Often a disturbance variable (feed composition changes)

T_feed : float, default=350.0 Feed temperature [K] - Inlet stream temperature - Typical range: 280-360 K (ambient to pre-heated) - Can be manipulated for control but usually fixed - Feed pre-heating can improve conversion but reduces stability margin

k0 : float, default=7.2e10 Pre-exponential factor [1/s] - Arrhenius equation parameter - Determines reaction speed at given temperature - Typical range: 10⁶-10¹² for liquid phase - Reaction and catalyst specific - Obtained from kinetic experiments or literature

E : float, default=8750.0 Activation energy [K] (actually Eₐ/R, dimensionless) - Energy barrier for reaction to occur - Higher E → more temperature-sensitive - Typical range: 5,000-15,000 K for Eₐ/R - Physical Eₐ typically 40-120 kJ/mol - Key parameter determining multiplicity region - Strong influence on stability

delta_H : float, default=-5e4 Heat of reaction [J/mol] - Energy released (negative) or absorbed (positive) per mole reacted - Negative = exothermic (releases heat) - most common case - Positive = endothermic (absorbs heat) - rare, simpler control - Typical for exothermic: -20,000 to -200,000 J/mol - Magnitude determines thermal coupling strength - Larger |ΔH| → stronger feedback → more multiplicity/instability

rho : float, default=1000.0 Density [kg/L] - Fluid density (assumed constant) - Typical for aqueous solutions: 900-1,100 kg/L - Affects thermal mass (heat capacity of reactor contents)

Cp : float, default=0.239 Specific heat capacity [J/(kg·K)] - Heat required to raise 1 kg by 1 K - Typical for aqueous: 0.2-0.5 J/(kg·K) - Note: Using J/(kg·K) not kJ/(kg·K), hence small value - Combined with ρ gives volumetric heat capacity ρCₚ - Higher Cₚ → larger thermal inertia → slower temperature changes

UA : float, default=5e4 Overall heat transfer coefficient × area [J/(s·K)] - Lumped parameter combining: * Jacket-side film coefficient * Wall thermal conductivity * Reactor-side film coefficient * Heat transfer area - Typical range: 10³-10⁵ J/(s·K) - Higher UA → better temperature control, faster cooling - Limited by physical design (jacket size, flow rate, area) - Critical parameter for preventing runaway

Equilibria and Multiple Steady States:

The Hallmark of CSTR Dynamics

The CSTR is famous for exhibiting multiple steady states - a phenomenon called multiplicity. For a given set of operating conditions (feed conditions and jacket temperature), the reactor can have:

One steady state (unique solution)
Three steady states (two stable, one unstable)
Five steady states (rare, special parameter combinations)

Three Steady State Scenario (most common interesting case):

Low-Conversion State (Stable):
- Low temperature (T ≈ T_feed + 10-30 K)
- High reactant concentration (Cₐ ≈ 0.7-0.9·Cₐ,feed)
- Slow reaction rate (low k·exp(-E/T))
- Heat generation < Heat removal capacity
- Basin of attraction: States with low initial temperature
- Characteristics:
  - Easy to start up (cold startup naturally goes here)
  - Safe and stable
  - Economically poor (low conversion, wasted reactant)
  - Easy to control (large stability margins)
Intermediate State (Unstable):
- Moderate temperature
- Moderate concentration
- Saddle point in phase space
- Not physically realizable (unstable equilibrium)
- Acts as separatrix between basins of attraction
- Physical meaning: Transition point where thermal generation rate exactly balances heat removal rate, but balance is unstable
High-Conversion State (Stable):
- High temperature (T ≈ T_feed + 50-100 K)
- Low reactant concentration (Cₐ ≈ 0.05-0.3·Cₐ,feed)
- Fast reaction rate (high k·exp(-E/T))
- Large heat generation balanced by cooling
- Basin of attraction: States with high initial temperature
- Characteristics:
  - Desired operating point (high conversion = high profit)
  - Requires good startup procedure (must cross unstable intermediate)
  - Risk of runaway if cooling fails
  - Smaller stability margins (closer to instability boundary)
  - More challenging to control (strong nonlinearity)

Physical Intuition for Multiple Steady States:

Imagine heat generation curve vs heat removal curve: - Heat generation: S-shaped (Arrhenius kinetics) * Low T: generation small (slow reaction) * Medium T: generation increases rapidly (exponential activation) * High T: generation levels off (reactant depletion)

Heat removal: Linear in T (jacket cooling)
- Straight line: q_removal = UA·(T - T_jacket)/VρCₚ
- Plus convective: (F/V)·(T - T_feed)

Intersections of these curves = steady states: - If removal line is steep (large UA): unique high-conversion state - If removal line is shallow (small UA): can have 3 intersections - As T_jacket varies: intersections appear/disappear (bifurcations)

Stability of Steady States:

Linear stability analysis (eigenvalues of Jacobian): - Stable: Both eigenvalues have negative real parts (LHP) * Perturbations decay back to steady state * Can operate here sustainably - Unstable: At least one eigenvalue has positive real part (RHP) * Perturbations grow exponentially * Cannot operate here (physically unrealizable)

For CSTR: - Low-conversion state: Typically stable - Intermediate state: Always unstable (saddle point) - High-conversion state: Stable if cooling is adequate

Bifurcations (Qualitative changes in behavior):

Saddle-Node Bifurcation (Fold): As T_jacket decreases (more cooling):
- Initially: Only low-conversion state exists
- At bifurcation point: Two new states appear (intermediate + high)
- Further cooling: Three states coexist
- Hysteresis: Different paths for heating vs cooling
Hopf Bifurcation (Oscillations): At certain parameters, high-conversion state can lose stability via Hopf bifurcation → sustained oscillations (limit cycle)
- Temperature and concentration oscillate periodically
- Can occur with insufficient cooling or long residence time
- Indicates poor controllability

Finding Steady States:

Steady states satisfy: dCₐ/dt = 0, dT/dt = 0

This gives two coupled nonlinear algebraic equations: 1. (F/V)·(Cₐ,feed - Cₐ) = k₀·Cₐ·exp(-E/T) 2. (F/V)·(T_feed - T) + (-ΔH/ρCₚ)·k₀·Cₐ·exp(-E/T) = -UA/(VρCₚ)·(T_jacket - T)

Solution methods: - Numerical: Newton-Raphson, fsolve with multiple initial guesses - Graphical: Plot dCₐ/dt and dT/dt surfaces, find intersections - Continuation: Track solutions as parameters vary (AUTO, MATCONT)

See find_steady_states() method for implementation.

Control Objectives:

Setpoint Regulation (Most Common):
- Maintain reactor at high-conversion steady state
- Reject disturbances (feed composition, flow rate, ambient temperature)
- Controllers: PID, LQR, MPC
- Challenge: Strong nonlinearity, operate near instability
Startup Control:
- Transition from low-conversion to high-conversion state
- Must cross unstable intermediate state (separatrix crossing)
- Requires aggressive transient cooling
- Strategies:
  - Bang-bang cooling (maximum jacket cooling)
  - Optimal control (minimize time/energy)
  - Gain scheduling (change controller as state changes)
- Risk: Overshoot → runaway
Runaway Prevention (Safety Critical):
- Detect incipient runaway conditions
- Implement emergency cooling/shutdown
- Monitor dT/dt (temperature rate of change)
- Constraint: T < T_max (safety limit)
- Last resort: Emergency cooling, feed shutoff, depressurization
Economic Optimization:
- Maximize profit = revenue - costs
- Revenue: Product value = price × F × conversion
- Costs: Cooling energy, reactant waste, equipment wear
- Often operates close to instability boundary for maximum conversion
- Tradeoff: Higher conversion (profit) vs safety margin (risk)
Disturbance Rejection: Common disturbances:
- Feed concentration variations: Cₐ,feed(t)
- Feed temperature changes: T_feed(t)
- Flow rate fluctuations: F(t)
- Ambient temperature (affects jacket): T_ambient(t)
- Catalyst deactivation: k₀(t) decreases slowly

State Constraints:

Non-negativity: Cₐ ≥ 0
- Concentration cannot be negative (physical)
- Rarely active (implies complete conversion)
Concentration Bounds: 0 ≤ Cₐ ≤ Cₐ,feed
- Cannot exceed feed (dilution + reaction only decrease)
- Upper bound active only if no reaction occurs
Temperature Limits: T_min ≤ T ≤ T_max
- Lower limit: Prevent solidification/freezing (≈ 280 K)
- Upper limit: Safety, prevent runaway (≈ 450-500 K)
- Material limits: polymer degradation, wall integrity
- Most critical constraint for safety
Jacket Temperature Limits: T_jacket,min ≤ T_jacket ≤ T_jacket,max
- Chiller capacity: T_jacket,min ≈ 280 K
- Heater capacity: T_jacket,max ≈ 400 K
- Rate limit: |dT_jacket/dt| ≤ rate_max (jacket dynamics)

Time Scales and Dynamics:

Multiple Time Scales make CSTR dynamics rich and challenging:

Fast Scale - Reaction: t_reaction = 1/k
- At low T (350 K): t_reaction ≈ 10-100 s
- At high T (400 K): t_reaction ≈ 0.1-1 s
- Exponentially dependent on temperature
- Can be very fast at high conversion state
Medium Scale - Residence Time: t_residence = τ = V/F
- Typical: 10-1000 s (seconds to minutes)
- Time for complete turnover of reactor contents
- Natural time scale for concentration changes
- Design parameter
Slow Scale - Thermal: t_thermal = VρCₚ/UA
- Typical: 100-10,000 s (minutes to hours)
- Time constant for temperature response
- Limited by heat transfer through jacket
- Design parameter (limited by area, jacket design)

Stiffness: When time scales differ by orders of magnitude: - Fast reactions with slow heat transfer → stiff system - Requires implicit ODE solvers (Radau, BDF) - Small numerical errors in fast variables → large errors in slow variables

Integration Recommendations:

Solver Selection:

For most CSTR problems: - Moderate stiffness: RK45 (adaptive Runge-Kutta) works well - High stiffness: Use stiff solvers * scipy: Radau, BDF, LSODA (auto-switching) * Julia (DiffEqPy): Rosenbrock23, Rodas5

Tolerance Selection: - Standard simulation: rtol=1e-6, atol=1e-8 - High accuracy (optimization): rtol=1e-9, atol=1e-11 - Looser tolerances may miss important dynamics

Event Detection: For safety-critical applications, use event detection: - Detect T > T_max (runaway) - Detect dT/dt > threshold (incipient runaway) - Detect steady state (convergence)

Example Usage:

Create CSTR with default parameters

cstr = ContinuousCSTR(F=100.0, V=100.0)

Find all steady states for given jacket temperature

T_jacket_op = 350.0 steady_states = cstr.find_steady_states(T_jacket_op) print(f”Found {len(steady_states)} steady states:“) for i, (C_A, T) in enumerate(steady_states): … print(f” State {i+1}: Cₐ={C_A:.3f} mol/L, T={T:.1f} K”) … print(f” Conversion: {cstr.compute_conversion(C_A, 1.0)*100:.1f}%“)

Choose high-conversion operating point

if len(steady_states) >= 2: … # High-conversion state (highest T, lowest Cₐ) … steady_states_sorted = sorted(steady_states, key=lambda x: x[1], reverse=True) … C_A_op, T_op = steady_states_sorted[0] … x_op = np.array([C_A_op, T_op]) … u_op = np.array([T_jacket_op]) … else: … # Use provided or default operating point … x_op = np.array([0.1, 390.0]) … u_op = np.array([350.0])

Verify equilibrium

dx = cstr(x_op, u_op) print(f”check: ||dx/dt|| = {np.linalg.norm(dx):.2e}“)

Linearize at operating point

A, B = cstr.linearize(x_op, u_op) eigenvalues = np.linalg.eigvals(A) print(f”eigenvalues: {eigenvalues}“) print(f”Stable: {np.all(np.real(eigenvalues) < 0)}“)

Design LQR controller (emphasize temperature control)

Q = np.diag([1.0, 100.0]) # Penalize temperature error heavily R = np.array([[1.0]]) lqr_result = cstr.control.design_lqr(A, B, Q, R, system_type=‘continuous’) K = lqr_result[‘gain’]

Simulate with LQR control and disturbance

def lqr_controller(x, t): … # Add feed temperature disturbance at t=50s … if t > 50: … # Disturbance increases effective T_feed by changing heat balance … # Compensate by reducing jacket temperature … disturbance_compensation = -2.0 … else: … disturbance_compensation = 0.0 … … u_fb = -K @ (x - x_op) + u_op … return u_fb + disturbance_compensation

Perturb from equilibrium

x0 = x_op + np.array([0.05, -5.0]) # Small perturbation

result = cstr.simulate( … x0=x0, … controller=lqr_controller, … t_span=(0, 200), … dt=0.1, … method=‘Radau’ # Stiff solver … )

Plot results

import matplotlib.pyplot as plt fig, axes = plt.subplots(3, 1, figsize=(10, 8))

Concentration

axes[0].plot(result[‘time’], result[‘states’][:, 0]) axes[0].axhline(x_op[0], color=‘r’, linestyle=‘–’, label=‘Setpoint’) axes[0].set_ylabel(‘Cₐ [mol/L]’) axes[0].legend() axes[0].grid(True)

Temperature

axes[1].plot(result[‘time’], result[‘states’][:, 1]) axes[1].axhline(x_op[1], color=‘r’, linestyle=‘–’, label=‘Setpoint’) axes[1].set_ylabel(‘T [K]’) axes[1].legend() axes[1].grid(True)

Control action

axes[2].plot(result[‘time’], result[‘controls’][:, 0]) axes[2].axhline(u_op[0], color=‘r’, linestyle=‘–’, label=‘Nominal’) axes[2].set_ylabel(‘T_jacket [K]’) axes[2].set_xlabel(‘Time [s]’) axes[2].legend() axes[2].grid(True)

plt.tight_layout() plt.show()

Startup simulation: low → high conversion

Start at low-conversion steady state

if len(steady_states) >= 3: … C_A_low, T_low = steady_states_sorted[-1] # Lowest temperature state … x_low = np.array([C_A_low, T_low]) … else: … x_low = np.array([0.9, 360.0])

def startup_controller(x, t): … ’‘’Aggressive cooling for startup, then switch to regulator’’’ … if t < 100: … # Phase 1: Aggressive cooling to jump to high-conversion state … return np.array([330.0]) # Cold jacket … else: … # Phase 2: LQR regulation around high-conversion setpoint … return lqr_controller(x, t)

result_startup = cstr.simulate( … x0=x_low, … controller=startup_controller, … t_span=(0, 300), … dt=0.1, … method=‘Radau’ … )

Check if startup succeeded

final_state = result_startup[‘states’][-1, :] distance = np.linalg.norm(final_state - x_op) print(f”result:“) print(f” Final state: Cₐ={final_state[0]:.3f}, T={final_state[1]:.1f}“) print(f” Distance to target: {distance:.3f}“) print(f” Success: {distance < 5.0}“)

Phase portrait (requires multiple simulations)

Shows basins of attraction for multiple steady states

Bifurcation diagram (vary T_jacket)

T_jacket_range = np.linspace(320, 360, 20) bifurcation_data = {‘T_jacket’: [], ‘C_A’: [], ‘T’: [], ‘stable’: []}

for Tj in T_jacket_range: … states = cstr.find_steady_states(Tj) … for C_A, T in states: … # Check stability … A_local, _ = cstr.linearize(np.array([C_A, T]), np.array([Tj])) … eigs = np.linalg.eigvals(A_local) … is_stable = np.all(np.real(eigs) < 0) … … bifurcation_data[‘T_jacket’].append(Tj) … bifurcation_data[‘C_A’].append(C_A) … bifurcation_data[‘T’].append(T) … bifurcation_data[‘stable’].append(is_stable)

Plot bifurcation diagram

plt.figure(figsize=(10, 6)) stable = np.array(bifurcation_data[‘stable’]) plt.plot( … np.array(bifurcation_data[‘T_jacket’])[stable], … np.array(bifurcation_data[‘T’])[stable], … ‘b-’, linewidth=2, label=‘Stable’ … ) plt.plot( … np.array(bifurcation_data[‘T_jacket’])[~stable], … np.array(bifurcation_data[‘T’])[~stable], … ‘r–’, linewidth=2, label=‘Unstable’ … ) plt.xlabel(‘Jacket Temperature [K]’) plt.ylabel(‘Reactor Temperature [K]’) plt.title(‘CSTR Bifurcation Diagram’) plt.legend() plt.grid(True) plt.show()

Physical Insights:

Why Multiple Steady States Occur:

The interplay between heat generation (nonlinear, S-shaped) and heat removal (linear) creates the possibility of multiple intersections:

Low T region:
- Reaction slow (small exp(-E/T))
- Generation curve starts flat
- Removal line dominates → low-conversion stable
Medium T region:
- Reaction accelerates rapidly (exponential growth)
- Generation curve steepens dramatically
- Can intersect removal line 3 times
High T region:
- Reaction very fast but Cₐ depleted
- Generation curve plateaus (limited by Cₐ)
- Removal continues linearly → high-conversion stable

Thermal Runaway Mechanism:

Positive feedback loop if cooling insufficient: 1. Small T increase (disturbance or control action) 2. Reaction rate jumps (exponential Arrhenius) 3. More heat generated (exothermic) 4. Temperature rises further 5. Loop continues → runaway!

Stabilizing mechanisms: - Reactant depletion (limits generation at high T) - Jacket cooling (removes heat) - Feed cooling (convective heat removal)

Industrial Significance:

CSTRs are ubiquitous in chemical industry: - Polymerization reactors - Pharmaceutical synthesis - Wastewater treatment (biological reactors) - Fermentation processes

Multiple steady states have practical implications: - Startup: Complex procedure to reach desired state - Control: Must prevent unintended switching - Safety: Runaway risk at high-conversion state - Economics: Higher conversion → more profit but harder control

Control Challenges:

Nonlinearity:
- Linear controllers (PID) may perform poorly
- Gain scheduling or nonlinear control needed
- Operating point dependent behavior
Instability risk:
- High-conversion state often close to instability
- Small disturbances can cause large excursions
- Need fast, aggressive control action
Constraints:
- Temperature limits (safety)
- Jacket temperature limits (physical)
- Actuator saturation degrades performance
Multiple states:
- System can “jump” between states
- Hysteresis complicates control
- Need to prevent unintended transitions

Design Considerations:

For industrial CSTR design: - Safety first: Adequate cooling capacity (large UA) - Residence time: Balance conversion vs stability (choose V/F) - Operating point: High conversion but safe margin from instability - Instrumentation: Fast, reliable temperature measurement - Emergency systems: Backup cooling, feed shutoff, pressure relief

Comparison with Other Reactors:

CSTR vs Batch Reactor: - CSTR: Continuous operation, steady state, higher throughput - Batch: Transient operation, finite time, better for small volumes

CSTR vs Plug Flow Reactor (PFR): - CSTR: Back-mixed, uniform concentration/temperature - PFR: No back-mixing, concentration/temperature gradients - PFR: Generally more efficient but harder to control

CSTR vs Semi-Batch: - CSTR: Continuous in/out - Semi-batch: Batch with continuous feed, better heat management

Notes:

Extensions: More complex CSTR models can include: - Multiple reactions (A → B → C, parallel reactions) - Non-ideal mixing (RTD, compartment models) - Jacket dynamics (first-order lag in cooling) - Catalyst deactivation (slow time scale) - pH effects (additional state) - Gas-liquid reactions (mass transfer limitations)

Methods

Name	Description
compute_conversion	Compute fractional conversion of reactant A.
compute_damkohler_number	Compute Damköhler number Da = k·τ (reaction rate × residence time).
compute_residence_time	Compute residence time τ = V/F.
define_system	Define continuous-time CSTR dynamics.
find_steady_states	Find all steady states for a given jacket temperature.
setup_equilibria	Set up steady-state equilibrium if provided.

compute_conversion

systems.builtin.deterministic.continuous.ContinuousCSTR.compute_conversion(
    C_A,
    C_A_feed,
)

Compute fractional conversion of reactant A.

Parameters

Name	Type	Description	Default
C_A	float	Current reactor concentration [mol/L]	required
C_A_feed	float	Feed concentration [mol/L]	required

Returns

Name	Type	Description
	float	Conversion fraction X_A = (C_A_feed - C_A) / C_A_feed

Examples

>>> cstr = ContinuousCSTR()
>>> X = cstr.compute_conversion(C_A=0.1, C_A_feed=1.0)
>>> print(f"Conversion: {X*100:.1f}%")
Conversion: 90.0%

Notes

High conversion (X > 0.9) typically corresponds to high-temperature steady state with fast kinetics and strong exothermic heat generation.

compute_damkohler_number

systems.builtin.deterministic.continuous.ContinuousCSTR.compute_damkohler_number(
    T,
)

Compute Damköhler number Da = k·τ (reaction rate × residence time).

Parameters

Name	Type	Description	Default
T	float	Temperature [K]	required

Returns

Name	Type	Description
	float	Damköhler number [dimensionless]

Notes

Damköhler number measures reaction rate relative to flow rate: - Da << 1: Reaction slow, flow dominates, low conversion - Da >> 1: Reaction fast, kinetics dominate, high conversion - Da ≈ 1: Balanced, optimal efficiency

Examples

>>> cstr = ContinuousCSTR()
>>> Da_low = cstr.compute_damkohler_number(T=350.0)
>>> Da_high = cstr.compute_damkohler_number(T=400.0)
>>> print(f"Da(350K) = {Da_low:.2f}")
>>> print(f"Da(400K) = {Da_high:.2f}")

compute_residence_time

systems.builtin.deterministic.continuous.ContinuousCSTR.compute_residence_time()

Compute residence time τ = V/F.

Returns

Name	Type	Description
	float	Residence time [s]

Examples

>>> cstr = ContinuousCSTR(F=100.0, V=100.0)
>>> tau = cstr.compute_residence_time()
>>> print(f"Residence time: {tau} s")
Residence time: 1.0 s

Notes

Residence time is the average time a molecule spends in the reactor. - Longer τ (smaller F): More conversion, less stable - Shorter τ (larger F): Less conversion, more stable

define_system

systems.builtin.deterministic.continuous.ContinuousCSTR.define_system(
    F_val=100.0,
    V_val=100.0,
    C_A_feed_val=1.0,
    T_feed_val=350.0,
    k0_val=72000000000.0,
    E_val=8750.0,
    delta_H_val=-50000.0,
    rho_val=1000.0,
    Cp_val=0.239,
    UA_val=50000.0,
    x_ss=None,
    u_ss=None,
)

Define continuous-time CSTR dynamics.

Parameters

Name	Type	Description	Default
F_val	float	Volumetric flow rate [L/s]	`100.0`
V_val	float	Reactor volume [L]	`100.0`
C_A_feed_val	float	Feed concentration [mol/L]	`1.0`
T_feed_val	float	Feed temperature [K]	`350.0`
k0_val	float	Pre-exponential factor [1/s]	`72000000000.0`
E_val	float	Activation energy [K] (dimensionless Eₐ/R)	`8750.0`
delta_H_val	float	Heat of reaction [J/mol] (negative = exothermic)	`-50000.0`
rho_val	float	Density [kg/L]	`1000.0`
Cp_val	float	Specific heat capacity [J/(kg·K)]	`0.239`
UA_val	float	Overall heat transfer coefficient × area [J/(s·K)]	`50000.0`
x_ss	Optional[np.ndarray]	Steady-state [Cₐ, T] for equilibrium setup	`None`
u_ss	Optional[np.ndarray]	Steady-state [T_jacket] for equilibrium setup	`None`

find_steady_states

systems.builtin.deterministic.continuous.ContinuousCSTR.find_steady_states(
    T_jacket,
    T_range=(300.0, 500.0),
    n_points=100,
)

Find all steady states for a given jacket temperature.

Uses multiple initial guesses across temperature range to find all solutions to the steady-state equations.

Parameters

Name	Type	Description	Default
T_jacket	float	Jacket temperature [K]	required
T_range	tuple	Temperature range to search (T_min, T_max) [K]	`(300.0, 500.0)`
n_points	int	Number of initial guesses for root finding	`100`

Returns

Name	Type	Description
	list	List of (C_A, T) steady state tuples

Examples

>>> cstr = ContinuousCSTR()
>>> steady_states = cstr.find_steady_states(T_jacket=350.0)
>>> print(f"Found {len(steady_states)} steady states")
>>> for i, (C_A, T) in enumerate(steady_states):
...     X = cstr.compute_conversion(C_A, 1.0)
...     print(f"  State {i+1}: Cₐ={C_A:.3f}, T={T:.1f}, X={X*100:.1f}%")

Notes

This method finds steady states by solving: dCₐ/dt = 0 (material balance) dT/dt = 0 (energy balance)

For CSTR, there can be 1, 2, or 3 steady states depending on parameters. This method attempts to find all of them by using many different initial guesses across the temperature range.

For production code, consider: - scipy.optimize.fsolve with multiple guesses - Continuation methods (AUTO, MATCONT) - Homotopy methods for guaranteed finding of all solutions

setup_equilibria

systems.builtin.deterministic.continuous.ContinuousCSTR.setup_equilibria()

Set up steady-state equilibrium if provided.

Notes

CSTR can have multiple steady states! Only add user-provided equilibrium. Finding all equilibria requires solving nonlinear algebraic equations (see find_steady_states() method).

# systems.builtin.deterministic.continuous.ContinuousCSTR { #cdesym.systems.builtin.deterministic.continuous.ContinuousCSTR } ```python systems.builtin.deterministic.continuous.ContinuousCSTR(*args, **kwargs) ``` Continuous-time Continuous Stirred-Tank Reactor (CSTR) with cooling jacket. ## Physical System: {.doc-section .doc-section-physical-system} A continuous flow reactor where reactant A converts to product B in an exothermic reaction. The CSTR is one of the most studied nonlinear systems in chemical engineering, exhibiting rich dynamics including: - Multiple steady states (multiplicity) - Sustained oscillations (limit cycles) - Bifurcations and hysteresis - Thermal runaway behavior Unlike batch reactors, CSTRs operate continuously with: - Continuous feed stream entering at Cₐ,feed, T_feed - Continuous product stream leaving at reactor conditions - Perfect mixing assumption (uniform concentration and temperature) - External cooling/heating via jacket The CSTR represents a fundamental model in: - Process control (benchmark nonlinear system) - Nonlinear dynamics (canonical example of multiplicity) - Chemical reaction engineering (industrial reactor design) - Bifurcation theory (illustrates saddle-node, Hopf bifurcations) ## State Space: {.doc-section .doc-section-state-space} State: x = [Cₐ, T] Concentration state: - Cₐ: Concentration of reactant A in reactor [mol/L] * Governed by material balance: in - out - reaction * 0 ≤ Cₐ ≤ Cₐ,feed (bounded by feed concentration) * Low Cₐ → high conversion (desired but challenging to control) * High Cₐ → low conversion (safe but inefficient) Temperature state: - T: Reactor temperature [K] * Governed by energy balance: in - out + generation - removal * Typically T > T_feed for exothermic reactions * Exhibits strong nonlinear coupling with concentration * Critical for safety (runaway prevention) * Small changes can cause large rate changes (Arrhenius) Control: u = [T_jacket] - T_jacket: Cooling jacket temperature [K] * Primary manipulated variable for temperature control * Affects heat removal rate via UA·(T - T_jacket) * Typically T_jacket < T (cooling mode) * Can be T_jacket > T (heating mode for startup/cold days) * Physical constraints: chiller/heater capacity limits * Rate constraints: jacket dynamics (not modeled here) Output: y = [Cₐ, T] - Full state measurement (common in modern plants) - In practice: * Cₐ: Online analyzer (GC, HPLC, NIR spectroscopy) * T: Thermocouple or RTD (fast, reliable) * Both have measurement noise and potential delays ## Dynamics: {.doc-section .doc-section-dynamics} The continuous-time CSTR dynamics are: dCₐ/dt = (F/V)·(Cₐ,feed - Cₐ) - r dT/dt = (F/V)·(T_feed - T) + (-ΔH/ρCₚ)·r + UA/(VρCₚ)·(T_jacket - T) **Reaction Rate (Arrhenius kinetics)**: r = k₀·Cₐ·exp(-E/T) [mol/(L·s)] where: - k₀: Pre-exponential factor [1/s] * Collision frequency in Arrhenius equation * Typical range: 10⁶-10¹² for liquid-phase reactions * Material and reaction-specific constant - E: Activation energy [K] (dimensionless Eₐ/R) * Energy barrier for reaction to occur * Typical range: 5,000-15,000 K for Eₐ/R * Higher E → more temperature-sensitive reaction * Physical activation energy Eₐ typically 40-120 kJ/mol - exp(-E/T): Arrhenius temperature dependence * Exponential sensitivity creates strong nonlinearity * 10°C change can double/triple reaction rate * Source of multiple steady states and instability **Physical Interpretation of Each Term**: Material Balance (dCₐ/dt): 1. **(F/V)·(Cₐ,feed - Cₐ)**: Convective in/out - F/V = 1/τ: Inverse residence time [1/s] - τ = V/F: Average time molecule spends in reactor [s] - Positive when Cₐ < Cₐ,feed (dilution effect) - Acts as "restoring force" toward feed concentration - Time scale: τ (seconds to minutes) 2. **-r**: Consumption by reaction - Always negative (reactant consumed) - Exponentially increases with temperature - Depends on current concentration (first-order) - Time scale: 1/k (fast at high T, slow at low T) At steady state: inflow rate = outflow rate + reaction rate Energy Balance (dT/dt): 1. **(F/V)·(T_feed - T)**: Convective heat in/out - Negative when T > T_feed (typical for exothermic) - Same time scale as material balance (1/τ) - Acts as "restoring force" toward feed temperature 2. **(-ΔH/ρCₚ)·r**: Heat generation from reaction - Positive for exothermic reaction (ΔH < 0) - Couples concentration to temperature - Creates positive feedback: higher T → faster r → more heat - This term causes thermal runaway if unchecked - Magnitude: |ΔH|/(ρCₚ) is adiabatic temperature rise per mol/L 3. **UA/(VρCₚ)·(T_jacket - T)**: Heat removal via jacket - Negative when T > T_jacket (cooling) - Only term controlled by manipulated variable - UA: Overall heat transfer coefficient × area [J/(s·K)] - Larger UA → better temperature control - Time scale: VρCₚ/UA (thermal time constant) At steady state: heat in + heat generated = heat out + heat removed **Nonlinear Coupling and Feedback**: The CSTR exhibits strong positive feedback that can lead to instability: 1. **Thermal Feedback Loop** (Runaway Mechanism): T ↑ → r ↑ (Arrhenius) → heat generation ↑ → T ↑ (positive feedback) This loop is stabilized by: - Convective cooling: higher T → more heat removal to feed - Jacket cooling: higher T → more heat removal via jacket - Reactant depletion: higher r → lower Cₐ → lower r (negative feedback) 2. **Material-Thermal Coupling**: - High T → fast reaction → low Cₐ (depletion) - Low Cₐ → slow reaction → less heat generation → lower T - This coupling creates multiple possible steady states 3. **Competition Between Time Scales**: - Residence time τ = V/F (convective transport) - Reaction time 1/k (chemical kinetics) - Thermal time VρCₚ/UA (heat transfer) - Relative magnitudes determine stability and multiplicity ## Parameters: {.doc-section .doc-section-parameters} F : float, default=100.0 Volumetric flow rate [L/s] - Controls residence time τ = V/F - Higher F → shorter τ → lower conversion but more stable - Lower F → longer τ → higher conversion but less stable - Typical range: 10-1000 L/s depending on reactor size - Often kept constant (set by upstream/downstream constraints) V : float, default=100.0 Reactor volume [L] - Combined with F to give residence time τ - Larger V → more material holdup → slower dynamics - Typical range: 100-10,000 L for industrial reactors - Design parameter (fixed once reactor is built) C_A_feed : float, default=1.0 Feed concentration [mol/L] - Upper bound for reactor concentration - Higher feed → more product but more heat generation - Typical range: 0.1-10 mol/L - Often a disturbance variable (feed composition changes) T_feed : float, default=350.0 Feed temperature [K] - Inlet stream temperature - Typical range: 280-360 K (ambient to pre-heated) - Can be manipulated for control but usually fixed - Feed pre-heating can improve conversion but reduces stability margin k0 : float, default=7.2e10 Pre-exponential factor [1/s] - Arrhenius equation parameter - Determines reaction speed at given temperature - Typical range: 10⁶-10¹² for liquid phase - Reaction and catalyst specific - Obtained from kinetic experiments or literature E : float, default=8750.0 Activation energy [K] (actually Eₐ/R, dimensionless) - Energy barrier for reaction to occur - Higher E → more temperature-sensitive - Typical range: 5,000-15,000 K for Eₐ/R - Physical Eₐ typically 40-120 kJ/mol - Key parameter determining multiplicity region - Strong influence on stability delta_H : float, default=-5e4 Heat of reaction [J/mol] - Energy released (negative) or absorbed (positive) per mole reacted - Negative = exothermic (releases heat) - most common case - Positive = endothermic (absorbs heat) - rare, simpler control - Typical for exothermic: -20,000 to -200,000 J/mol - Magnitude determines thermal coupling strength - Larger |ΔH| → stronger feedback → more multiplicity/instability rho : float, default=1000.0 Density [kg/L] - Fluid density (assumed constant) - Typical for aqueous solutions: 900-1,100 kg/L - Affects thermal mass (heat capacity of reactor contents) Cp : float, default=0.239 Specific heat capacity [J/(kg·K)] - Heat required to raise 1 kg by 1 K - Typical for aqueous: 0.2-0.5 J/(kg·K) - Note: Using J/(kg·K) not kJ/(kg·K), hence small value - Combined with ρ gives volumetric heat capacity ρCₚ - Higher Cₚ → larger thermal inertia → slower temperature changes UA : float, default=5e4 Overall heat transfer coefficient × area [J/(s·K)] - Lumped parameter combining: * Jacket-side film coefficient * Wall thermal conductivity * Reactor-side film coefficient * Heat transfer area - Typical range: 10³-10⁵ J/(s·K) - Higher UA → better temperature control, faster cooling - Limited by physical design (jacket size, flow rate, area) - Critical parameter for preventing runaway ## Equilibria and Multiple Steady States: {.doc-section .doc-section-equilibria-and-multiple-steady-states} **The Hallmark of CSTR Dynamics** The CSTR is famous for exhibiting **multiple steady states** - a phenomenon called **multiplicity**. For a given set of operating conditions (feed conditions and jacket temperature), the reactor can have: 1. **One steady state** (unique solution) 2. **Three steady states** (two stable, one unstable) 3. **Five steady states** (rare, special parameter combinations) **Three Steady State Scenario** (most common interesting case): 1. **Low-Conversion State (Stable)**: - Low temperature (T ≈ T_feed + 10-30 K) - High reactant concentration (Cₐ ≈ 0.7-0.9·Cₐ,feed) - Slow reaction rate (low k·exp(-E/T)) - Heat generation < Heat removal capacity - **Basin of attraction**: States with low initial temperature - **Characteristics**: * Easy to start up (cold startup naturally goes here) * Safe and stable * Economically poor (low conversion, wasted reactant) * Easy to control (large stability margins) 2. **Intermediate State (Unstable)**: - Moderate temperature - Moderate concentration - **Saddle point** in phase space - Not physically realizable (unstable equilibrium) - Acts as separatrix between basins of attraction - **Physical meaning**: Transition point where thermal generation rate exactly balances heat removal rate, but balance is unstable 3. **High-Conversion State (Stable)**: - High temperature (T ≈ T_feed + 50-100 K) - Low reactant concentration (Cₐ ≈ 0.05-0.3·Cₐ,feed) - Fast reaction rate (high k·exp(-E/T)) - Large heat generation balanced by cooling - **Basin of attraction**: States with high initial temperature - **Characteristics**: * Desired operating point (high conversion = high profit) * Requires good startup procedure (must cross unstable intermediate) * Risk of runaway if cooling fails * Smaller stability margins (closer to instability boundary) * More challenging to control (strong nonlinearity) **Physical Intuition for Multiple Steady States**: Imagine heat generation curve vs heat removal curve: - Heat generation: S-shaped (Arrhenius kinetics) * Low T: generation small (slow reaction) * Medium T: generation increases rapidly (exponential activation) * High T: generation levels off (reactant depletion) - Heat removal: Linear in T (jacket cooling) * Straight line: q_removal = UA·(T - T_jacket)/VρCₚ * Plus convective: (F/V)·(T - T_feed) Intersections of these curves = steady states: - If removal line is steep (large UA): unique high-conversion state - If removal line is shallow (small UA): can have 3 intersections - As T_jacket varies: intersections appear/disappear (bifurcations) **Stability of Steady States**: Linear stability analysis (eigenvalues of Jacobian): - **Stable**: Both eigenvalues have negative real parts (LHP) * Perturbations decay back to steady state * Can operate here sustainably - **Unstable**: At least one eigenvalue has positive real part (RHP) * Perturbations grow exponentially * Cannot operate here (physically unrealizable) For CSTR: - Low-conversion state: Typically stable - Intermediate state: Always unstable (saddle point) - High-conversion state: Stable if cooling is adequate **Bifurcations** (Qualitative changes in behavior): 1. **Saddle-Node Bifurcation** (Fold): As T_jacket decreases (more cooling): - Initially: Only low-conversion state exists - At bifurcation point: Two new states appear (intermediate + high) - Further cooling: Three states coexist - **Hysteresis**: Different paths for heating vs cooling 2. **Hopf Bifurcation** (Oscillations): At certain parameters, high-conversion state can lose stability via Hopf bifurcation → sustained oscillations (limit cycle) - Temperature and concentration oscillate periodically - Can occur with insufficient cooling or long residence time - Indicates poor controllability **Finding Steady States**: Steady states satisfy: dCₐ/dt = 0, dT/dt = 0 This gives two coupled nonlinear algebraic equations: 1. (F/V)·(Cₐ,feed - Cₐ) = k₀·Cₐ·exp(-E/T) 2. (F/V)·(T_feed - T) + (-ΔH/ρCₚ)·k₀·Cₐ·exp(-E/T) = -UA/(VρCₚ)·(T_jacket - T) Solution methods: - Numerical: Newton-Raphson, fsolve with multiple initial guesses - Graphical: Plot dCₐ/dt and dT/dt surfaces, find intersections - Continuation: Track solutions as parameters vary (AUTO, MATCONT) See `find_steady_states()` method for implementation. ## Control Objectives: {.doc-section .doc-section-control-objectives} 1. **Setpoint Regulation** (Most Common): - Maintain reactor at high-conversion steady state - Reject disturbances (feed composition, flow rate, ambient temperature) - Controllers: PID, LQR, MPC - Challenge: Strong nonlinearity, operate near instability 2. **Startup Control**: - Transition from low-conversion to high-conversion state - Must cross unstable intermediate state (separatrix crossing) - Requires aggressive transient cooling - Strategies: * Bang-bang cooling (maximum jacket cooling) * Optimal control (minimize time/energy) * Gain scheduling (change controller as state changes) - Risk: Overshoot → runaway 3. **Runaway Prevention** (Safety Critical): - Detect incipient runaway conditions - Implement emergency cooling/shutdown - Monitor dT/dt (temperature rate of change) - Constraint: T < T_max (safety limit) - Last resort: Emergency cooling, feed shutoff, depressurization 4. **Economic Optimization**: - Maximize profit = revenue - costs - Revenue: Product value = price × F × conversion - Costs: Cooling energy, reactant waste, equipment wear - Often operates close to instability boundary for maximum conversion - Tradeoff: Higher conversion (profit) vs safety margin (risk) 5. **Disturbance Rejection**: Common disturbances: - Feed concentration variations: Cₐ,feed(t) - Feed temperature changes: T_feed(t) - Flow rate fluctuations: F(t) - Ambient temperature (affects jacket): T_ambient(t) - Catalyst deactivation: k₀(t) decreases slowly ## State Constraints: {.doc-section .doc-section-state-constraints} 1. **Non-negativity**: Cₐ ≥ 0 - Concentration cannot be negative (physical) - Rarely active (implies complete conversion) 2. **Concentration Bounds**: 0 ≤ Cₐ ≤ Cₐ,feed - Cannot exceed feed (dilution + reaction only decrease) - Upper bound active only if no reaction occurs 3. **Temperature Limits**: T_min ≤ T ≤ T_max - **Lower limit**: Prevent solidification/freezing (≈ 280 K) - **Upper limit**: Safety, prevent runaway (≈ 450-500 K) - Material limits: polymer degradation, wall integrity - **Most critical constraint for safety** 4. **Jacket Temperature Limits**: T_jacket,min ≤ T_jacket ≤ T_jacket,max - Chiller capacity: T_jacket,min ≈ 280 K - Heater capacity: T_jacket,max ≈ 400 K - Rate limit: |dT_jacket/dt| ≤ rate_max (jacket dynamics) ## Time Scales and Dynamics: {.doc-section .doc-section-time-scales-and-dynamics} **Multiple Time Scales** make CSTR dynamics rich and challenging: 1. **Fast Scale - Reaction**: t_reaction = 1/k - At low T (350 K): t_reaction ≈ 10-100 s - At high T (400 K): t_reaction ≈ 0.1-1 s - Exponentially dependent on temperature - Can be very fast at high conversion state 2. **Medium Scale - Residence Time**: t_residence = τ = V/F - Typical: 10-1000 s (seconds to minutes) - Time for complete turnover of reactor contents - Natural time scale for concentration changes - Design parameter 3. **Slow Scale - Thermal**: t_thermal = VρCₚ/UA - Typical: 100-10,000 s (minutes to hours) - Time constant for temperature response - Limited by heat transfer through jacket - Design parameter (limited by area, jacket design) **Stiffness**: When time scales differ by orders of magnitude: - Fast reactions with slow heat transfer → stiff system - Requires implicit ODE solvers (Radau, BDF) - Small numerical errors in fast variables → large errors in slow variables ## Integration Recommendations: {.doc-section .doc-section-integration-recommendations} **Solver Selection**: For most CSTR problems: - **Moderate stiffness**: RK45 (adaptive Runge-Kutta) works well - **High stiffness**: Use stiff solvers * scipy: Radau, BDF, LSODA (auto-switching) * Julia (DiffEqPy): Rosenbrock23, Rodas5 **Tolerance Selection**: - Standard simulation: rtol=1e-6, atol=1e-8 - High accuracy (optimization): rtol=1e-9, atol=1e-11 - Looser tolerances may miss important dynamics **Event Detection**: For safety-critical applications, use event detection: - Detect T > T_max (runaway) - Detect dT/dt > threshold (incipient runaway) - Detect steady state (convergence) ## Example Usage: {.doc-section .doc-section-example-usage} >>> # Create CSTR with default parameters >>> cstr = ContinuousCSTR(F=100.0, V=100.0) >>> >>> # Find all steady states for given jacket temperature >>> T_jacket_op = 350.0 >>> steady_states = cstr.find_steady_states(T_jacket_op) >>> print(f"Found {len(steady_states)} steady states:") >>> for i, (C_A, T) in enumerate(steady_states): ... print(f" State {i+1}: Cₐ={C_A:.3f} mol/L, T={T:.1f} K") ... print(f" Conversion: {cstr.compute_conversion(C_A, 1.0)*100:.1f}%") >>> >>> # Choose high-conversion operating point >>> if len(steady_states) >= 2: ... # High-conversion state (highest T, lowest Cₐ) ... steady_states_sorted = sorted(steady_states, key=lambda x: x[1], reverse=True) ... C_A_op, T_op = steady_states_sorted[0] ... x_op = np.array([C_A_op, T_op]) ... u_op = np.array([T_jacket_op]) ... else: ... # Use provided or default operating point ... x_op = np.array([0.1, 390.0]) ... u_op = np.array([350.0]) >>> >>> # Verify equilibrium >>> dx = cstr(x_op, u_op) >>> print(f"\nEquilibrium check: ||dx/dt|| = {np.linalg.norm(dx):.2e}") >>> >>> # Linearize at operating point >>> A, B = cstr.linearize(x_op, u_op) >>> eigenvalues = np.linalg.eigvals(A) >>> print(f"\nLinearized eigenvalues: {eigenvalues}") >>> print(f"Stable: {np.all(np.real(eigenvalues) < 0)}") >>> >>> # Design LQR controller (emphasize temperature control) >>> Q = np.diag([1.0, 100.0]) # Penalize temperature error heavily >>> R = np.array([[1.0]]) >>> lqr_result = cstr.control.design_lqr(A, B, Q, R, system_type='continuous') >>> K = lqr_result['gain'] >>> >>> # Simulate with LQR control and disturbance >>> def lqr_controller(x, t): ... # Add feed temperature disturbance at t=50s ... if t > 50: ... # Disturbance increases effective T_feed by changing heat balance ... # Compensate by reducing jacket temperature ... disturbance_compensation = -2.0 ... else: ... disturbance_compensation = 0.0 ... ... u_fb = -K @ (x - x_op) + u_op ... return u_fb + disturbance_compensation >>> >>> # Perturb from equilibrium >>> x0 = x_op + np.array([0.05, -5.0]) # Small perturbation >>> >>> result = cstr.simulate( ... x0=x0, ... controller=lqr_controller, ... t_span=(0, 200), ... dt=0.1, ... method='Radau' # Stiff solver ... ) >>> >>> # Plot results >>> import matplotlib.pyplot as plt >>> fig, axes = plt.subplots(3, 1, figsize=(10, 8)) >>> >>> # Concentration >>> axes[0].plot(result['time'], result['states'][:, 0]) >>> axes[0].axhline(x_op[0], color='r', linestyle='--', label='Setpoint') >>> axes[0].set_ylabel('Cₐ [mol/L]') >>> axes[0].legend() >>> axes[0].grid(True) >>> >>> # Temperature >>> axes[1].plot(result['time'], result['states'][:, 1]) >>> axes[1].axhline(x_op[1], color='r', linestyle='--', label='Setpoint') >>> axes[1].set_ylabel('T [K]') >>> axes[1].legend() >>> axes[1].grid(True) >>> >>> # Control action >>> axes[2].plot(result['time'], result['controls'][:, 0]) >>> axes[2].axhline(u_op[0], color='r', linestyle='--', label='Nominal') >>> axes[2].set_ylabel('T_jacket [K]') >>> axes[2].set_xlabel('Time [s]') >>> axes[2].legend() >>> axes[2].grid(True) >>> >>> plt.tight_layout() >>> plt.show() >>> >>> # Startup simulation: low → high conversion >>> # Start at low-conversion steady state >>> if len(steady_states) >= 3: ... C_A_low, T_low = steady_states_sorted[-1] # Lowest temperature state ... x_low = np.array([C_A_low, T_low]) ... else: ... x_low = np.array([0.9, 360.0]) >>> >>> def startup_controller(x, t): ... '''Aggressive cooling for startup, then switch to regulator''' ... if t < 100: ... # Phase 1: Aggressive cooling to jump to high-conversion state ... return np.array([330.0]) # Cold jacket ... else: ... # Phase 2: LQR regulation around high-conversion setpoint ... return lqr_controller(x, t) >>> >>> result_startup = cstr.simulate( ... x0=x_low, ... controller=startup_controller, ... t_span=(0, 300), ... dt=0.1, ... method='Radau' ... ) >>> >>> # Check if startup succeeded >>> final_state = result_startup['states'][-1, :] >>> distance = np.linalg.norm(final_state - x_op) >>> print(f"\nStartup result:") >>> print(f" Final state: Cₐ={final_state[0]:.3f}, T={final_state[1]:.1f}") >>> print(f" Distance to target: {distance:.3f}") >>> print(f" Success: {distance < 5.0}") >>> >>> # Phase portrait (requires multiple simulations) >>> # Shows basins of attraction for multiple steady states >>> >>> # Bifurcation diagram (vary T_jacket) >>> T_jacket_range = np.linspace(320, 360, 20) >>> bifurcation_data = {'T_jacket': [], 'C_A': [], 'T': [], 'stable': []} >>> >>> for Tj in T_jacket_range: ... states = cstr.find_steady_states(Tj) ... for C_A, T in states: ... # Check stability ... A_local, _ = cstr.linearize(np.array([C_A, T]), np.array([Tj])) ... eigs = np.linalg.eigvals(A_local) ... is_stable = np.all(np.real(eigs) < 0) ... ... bifurcation_data['T_jacket'].append(Tj) ... bifurcation_data['C_A'].append(C_A) ... bifurcation_data['T'].append(T) ... bifurcation_data['stable'].append(is_stable) >>> >>> # Plot bifurcation diagram >>> plt.figure(figsize=(10, 6)) >>> stable = np.array(bifurcation_data['stable']) >>> plt.plot( ... np.array(bifurcation_data['T_jacket'])[stable], ... np.array(bifurcation_data['T'])[stable], ... 'b-', linewidth=2, label='Stable' ... ) >>> plt.plot( ... np.array(bifurcation_data['T_jacket'])[~stable], ... np.array(bifurcation_data['T'])[~stable], ... 'r--', linewidth=2, label='Unstable' ... ) >>> plt.xlabel('Jacket Temperature [K]') >>> plt.ylabel('Reactor Temperature [K]') >>> plt.title('CSTR Bifurcation Diagram') >>> plt.legend() >>> plt.grid(True) >>> plt.show() ## Physical Insights: {.doc-section .doc-section-physical-insights} **Why Multiple Steady States Occur**: The interplay between heat generation (nonlinear, S-shaped) and heat removal (linear) creates the possibility of multiple intersections: 1. **Low T region**: - Reaction slow (small exp(-E/T)) - Generation curve starts flat - Removal line dominates → low-conversion stable 2. **Medium T region**: - Reaction accelerates rapidly (exponential growth) - Generation curve steepens dramatically - Can intersect removal line 3 times 3. **High T region**: - Reaction very fast but Cₐ depleted - Generation curve plateaus (limited by Cₐ) - Removal continues linearly → high-conversion stable **Thermal Runaway Mechanism**: Positive feedback loop if cooling insufficient: 1. Small T increase (disturbance or control action) 2. Reaction rate jumps (exponential Arrhenius) 3. More heat generated (exothermic) 4. Temperature rises further 5. Loop continues → runaway! Stabilizing mechanisms: - Reactant depletion (limits generation at high T) - Jacket cooling (removes heat) - Feed cooling (convective heat removal) **Industrial Significance**: CSTRs are ubiquitous in chemical industry: - Polymerization reactors - Pharmaceutical synthesis - Wastewater treatment (biological reactors) - Fermentation processes Multiple steady states have practical implications: - **Startup**: Complex procedure to reach desired state - **Control**: Must prevent unintended switching - **Safety**: Runaway risk at high-conversion state - **Economics**: Higher conversion → more profit but harder control **Control Challenges**: 1. **Nonlinearity**: - Linear controllers (PID) may perform poorly - Gain scheduling or nonlinear control needed - Operating point dependent behavior 2. **Instability risk**: - High-conversion state often close to instability - Small disturbances can cause large excursions - Need fast, aggressive control action 3. **Constraints**: - Temperature limits (safety) - Jacket temperature limits (physical) - Actuator saturation degrades performance 4. **Multiple states**: - System can "jump" between states - Hysteresis complicates control - Need to prevent unintended transitions **Design Considerations**: For industrial CSTR design: - **Safety first**: Adequate cooling capacity (large UA) - **Residence time**: Balance conversion vs stability (choose V/F) - **Operating point**: High conversion but safe margin from instability - **Instrumentation**: Fast, reliable temperature measurement - **Emergency systems**: Backup cooling, feed shutoff, pressure relief ## Comparison with Other Reactors: {.doc-section .doc-section-comparison-with-other-reactors} **CSTR vs Batch Reactor**: - CSTR: Continuous operation, steady state, higher throughput - Batch: Transient operation, finite time, better for small volumes **CSTR vs Plug Flow Reactor (PFR)**: - CSTR: Back-mixed, uniform concentration/temperature - PFR: No back-mixing, concentration/temperature gradients - PFR: Generally more efficient but harder to control **CSTR vs Semi-Batch**: - CSTR: Continuous in/out - Semi-batch: Batch with continuous feed, better heat management ## See Also: {.doc-section .doc-section-see-also} DiscreteCSTR : Discrete-time version for digital control ContinuousBatchReactor : Batch operation instead of continuous DiscreteBatchReactor : Discrete batch reactor ## Notes: {.doc-section .doc-section-notes} **Extensions**: More complex CSTR models can include: - Multiple reactions (A → B → C, parallel reactions) - Non-ideal mixing (RTD, compartment models) - Jacket dynamics (first-order lag in cooling) - Catalyst deactivation (slow time scale) - pH effects (additional state) - Gas-liquid reactions (mass transfer limitations) ## Methods | Name | Description | | --- | --- | | [compute_conversion](#cdesym.systems.builtin.deterministic.continuous.ContinuousCSTR.compute_conversion) | Compute fractional conversion of reactant A. | | [compute_damkohler_number](#cdesym.systems.builtin.deterministic.continuous.ContinuousCSTR.compute_damkohler_number) | Compute Damköhler number Da = k·τ (reaction rate × residence time). | | [compute_residence_time](#cdesym.systems.builtin.deterministic.continuous.ContinuousCSTR.compute_residence_time) | Compute residence time τ = V/F. | | [define_system](#cdesym.systems.builtin.deterministic.continuous.ContinuousCSTR.define_system) | Define continuous-time CSTR dynamics. | | [find_steady_states](#cdesym.systems.builtin.deterministic.continuous.ContinuousCSTR.find_steady_states) | Find all steady states for a given jacket temperature. | | [setup_equilibria](#cdesym.systems.builtin.deterministic.continuous.ContinuousCSTR.setup_equilibria) | Set up steady-state equilibrium if provided. | ### compute_conversion { #cdesym.systems.builtin.deterministic.continuous.ContinuousCSTR.compute_conversion } ```python systems.builtin.deterministic.continuous.ContinuousCSTR.compute_conversion( C_A, C_A_feed, ) ``` Compute fractional conversion of reactant A. #### Parameters {.doc-section .doc-section-parameters} | Name | Type | Description | Default | |----------|--------|---------------------------------------|------------| | C_A | float | Current reactor concentration [mol/L] | _required_ | | C_A_feed | float | Feed concentration [mol/L] | _required_ | #### Returns {.doc-section .doc-section-returns} | Name | Type | Description | |--------|--------|-------------------------------------------------------| | | float | Conversion fraction X_A = (C_A_feed - C_A) / C_A_feed | #### Examples {.doc-section .doc-section-examples} ```python >>> cstr = ContinuousCSTR() >>> X = cstr.compute_conversion(C_A=0.1, C_A_feed=1.0) >>> print(f"Conversion: {X*100:.1f}%") Conversion: 90.0% ``` #### Notes {.doc-section .doc-section-notes} High conversion (X > 0.9) typically corresponds to high-temperature steady state with fast kinetics and strong exothermic heat generation. ### compute_damkohler_number { #cdesym.systems.builtin.deterministic.continuous.ContinuousCSTR.compute_damkohler_number } ```python systems.builtin.deterministic.continuous.ContinuousCSTR.compute_damkohler_number( T, ) ``` Compute Damköhler number Da = k·τ (reaction rate × residence time). #### Parameters {.doc-section .doc-section-parameters} | Name | Type | Description | Default | |--------|--------|-----------------|------------| | T | float | Temperature [K] | _required_ | #### Returns {.doc-section .doc-section-returns} | Name | Type | Description | |--------|--------|----------------------------------| | | float | Damköhler number [dimensionless] | #### Notes {.doc-section .doc-section-notes} Damköhler number measures reaction rate relative to flow rate: - Da << 1: Reaction slow, flow dominates, low conversion - Da >> 1: Reaction fast, kinetics dominate, high conversion - Da ≈ 1: Balanced, optimal efficiency #### Examples {.doc-section .doc-section-examples} ```python >>> cstr = ContinuousCSTR() >>> Da_low = cstr.compute_damkohler_number(T=350.0) >>> Da_high = cstr.compute_damkohler_number(T=400.0) >>> print(f"Da(350K) = {Da_low:.2f}") >>> print(f"Da(400K) = {Da_high:.2f}") ``` ### compute_residence_time { #cdesym.systems.builtin.deterministic.continuous.ContinuousCSTR.compute_residence_time } ```python systems.builtin.deterministic.continuous.ContinuousCSTR.compute_residence_time() ``` Compute residence time τ = V/F. #### Returns {.doc-section .doc-section-returns} | Name | Type | Description | |--------|--------|--------------------| | | float | Residence time [s] | #### Examples {.doc-section .doc-section-examples} ```python >>> cstr = ContinuousCSTR(F=100.0, V=100.0) >>> tau = cstr.compute_residence_time() >>> print(f"Residence time: {tau} s") Residence time: 1.0 s ``` #### Notes {.doc-section .doc-section-notes} Residence time is the average time a molecule spends in the reactor. - Longer τ (smaller F): More conversion, less stable - Shorter τ (larger F): Less conversion, more stable ### define_system { #cdesym.systems.builtin.deterministic.continuous.ContinuousCSTR.define_system } ```python systems.builtin.deterministic.continuous.ContinuousCSTR.define_system( F_val=100.0, V_val=100.0, C_A_feed_val=1.0, T_feed_val=350.0, k0_val=72000000000.0, E_val=8750.0, delta_H_val=-50000.0, rho_val=1000.0, Cp_val=0.239, UA_val=50000.0, x_ss=None, u_ss=None, ) ``` Define continuous-time CSTR dynamics. #### Parameters {.doc-section .doc-section-parameters} | Name | Type | Description | Default | |--------------|------------------------|----------------------------------------------------|-----------------| | F_val | float | Volumetric flow rate [L/s] | `100.0` | | V_val | float | Reactor volume [L] | `100.0` | | C_A_feed_val | float | Feed concentration [mol/L] | `1.0` | | T_feed_val | float | Feed temperature [K] | `350.0` | | k0_val | float | Pre-exponential factor [1/s] | `72000000000.0` | | E_val | float | Activation energy [K] (dimensionless Eₐ/R) | `8750.0` | | delta_H_val | float | Heat of reaction [J/mol] (negative = exothermic) | `-50000.0` | | rho_val | float | Density [kg/L] | `1000.0` | | Cp_val | float | Specific heat capacity [J/(kg·K)] | `0.239` | | UA_val | float | Overall heat transfer coefficient × area [J/(s·K)] | `50000.0` | | x_ss | Optional\[np.ndarray\] | Steady-state [Cₐ, T] for equilibrium setup | `None` | | u_ss | Optional\[np.ndarray\] | Steady-state [T_jacket] for equilibrium setup | `None` | ### find_steady_states { #cdesym.systems.builtin.deterministic.continuous.ContinuousCSTR.find_steady_states } ```python systems.builtin.deterministic.continuous.ContinuousCSTR.find_steady_states( T_jacket, T_range=(300.0, 500.0), n_points=100, ) ``` Find all steady states for a given jacket temperature. Uses multiple initial guesses across temperature range to find all solutions to the steady-state equations. #### Parameters {.doc-section .doc-section-parameters} | Name | Type | Description | Default | |----------|--------|------------------------------------------------|------------------| | T_jacket | float | Jacket temperature [K] | _required_ | | T_range | tuple | Temperature range to search (T_min, T_max) [K] | `(300.0, 500.0)` | | n_points | int | Number of initial guesses for root finding | `100` | #### Returns {.doc-section .doc-section-returns} | Name | Type | Description | |--------|--------|--------------------------------------| | | list | List of (C_A, T) steady state tuples | #### Examples {.doc-section .doc-section-examples} ```python >>> cstr = ContinuousCSTR() >>> steady_states = cstr.find_steady_states(T_jacket=350.0) >>> print(f"Found {len(steady_states)} steady states") >>> for i, (C_A, T) in enumerate(steady_states): ... X = cstr.compute_conversion(C_A, 1.0) ... print(f" State {i+1}: Cₐ={C_A:.3f}, T={T:.1f}, X={X*100:.1f}%") ``` #### Notes {.doc-section .doc-section-notes} This method finds steady states by solving: dCₐ/dt = 0 (material balance) dT/dt = 0 (energy balance) For CSTR, there can be 1, 2, or 3 steady states depending on parameters. This method attempts to find all of them by using many different initial guesses across the temperature range. For production code, consider: - scipy.optimize.fsolve with multiple guesses - Continuation methods (AUTO, MATCONT) - Homotopy methods for guaranteed finding of all solutions ### setup_equilibria { #cdesym.systems.builtin.deterministic.continuous.ContinuousCSTR.setup_equilibria } ```python systems.builtin.deterministic.continuous.ContinuousCSTR.setup_equilibria() ``` Set up steady-state equilibrium if provided. #### Notes {.doc-section .doc-section-notes} CSTR can have multiple steady states! Only add user-provided equilibrium. Finding all equilibria requires solving nonlinear algebraic equations (see find_steady_states() method).