6. Radical Semibatch Polymerization and Safety#
6.1. 馃摉 Introduction#
Semibatch processes are widely employed in the polymer industry due to their inherent flexibility, which enables precise control of reaction conditions and polymer properties. In this notebook, we will illustrate the basic principles needed to model these systems by examining a classical radical solution homopolymerization, based on the reaction scheme outlined in Notebook 3.
We assume a single-phase, well-stirred tank reactor where the mole balances for the species involved can be expressed as:
Here, \(n_i\), \(\dot{n}^{in}_{i}\) and \(r_i\) denote, respectively, the number of moles, the molar inflow, and the rate of production of species \(i\), while \(V\) represents the volume of the reaction mixture. The relevant species in our case include the initiator (\(I\)), monomer (\(M\)), optional solvent (\(S\)), radicals (\(R^{\cdot}_n\)), and dead polymer (\(P_n\)).
Note
One sometimes encounters semibatch models formulated with concentrations and volume as state variables. This can lead to unnecessary numerical issues (can you think of two examples?) and should be avoided. The state variables should preferably be extensive quantities like mass or moles, the latter having the advantage of working also for massless pseudo-species.
Let us further suppose that we are only interested in the average polymer chain length, rather than in the full molar mass distribution. In this scenario, the method of moments can be used to reduce the complexity of the model. Specifically, we focus on the first three moments of the radical and polymer distributions:
with \(j \in \{0,1,2\}\). As a result, the model will consist of molar balances for three species (\(I\), \(M\), \(S\)) and six pseudo-species (\(\mu_0\), \(\mu_1\), \(\mu_2\), \(\lambda_0\), \(\lambda_1\), and \(\lambda_2\)). The corresponding rates of formation are:
where \(k_t = k_{tc} + k_{td}\).
The volume of the reaction mixture can be evaluated from the calculated mole amounts of monomer, solvent and polymer repeating units (initiator neglected):
where \(MW_i\) and \(\rho_i\) represent, respectively, the molar mass and density of species \(i\).
Given that chain-growth reactions are strongly exothermal, the model would not be complete without an energy balance:
where \(m_t\) is the total mass of the reaction mixture, \(c_P\) is the average specific heat capacity of the reaction mixture and feed, \(\dot{m}_{in}\) is the total inlet mass flowrate, \(T_{in}\) is the mean feed temperature, \(\Delta H_{pol}\) is the polymerization enthalpy, \(U\) is the overall heat transfer coefficient, \(A\) is the heat transfer area, and \(T_j\) is the jacket temperature. We assume the jacket temperature is automatically adjusted by means of a simple proportional-only controller:
where \(T_{sp}(t)\) is the reactor temperature setpoint. Additionally, we also estimate the maximum temperature of the synthesis reaction (MTSR), which is an important safety quantity:
Enough theory! Let鈥檚 write some code to solve this semibatch homopolymerization model.
6.2. 馃М Model Implementation#
from enum import IntEnum
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
To improve code readability, we define an enum mapping the variable names to the corresponding row index in the state vector and related arrays.
class Ix(IntEnum):
I = 0
M = 1
S = 2
mu0 = 3
mu1 = 4
mu2 = 5
lbd0 = 6
lbd1 = 7
lbd2 = 8
T = 9
This model is more complex than the preceding ones, so it is essential to organize the code in a modular way. First, we implement a helper function arrhenius to describe the temperature dependence of the rate coefficients.
def arrhenius(T: float, k0: float, Ea: float, T0: float) -> float:
"""Arrhenius equation.
Parameters
----------
T : float
Temperature (掳C).
k0 : float
Value of rate coefficient at `T0` (#/s).
Ea : float
Energy of activation (J/mol).
T0 : float
Reference temperature (掳C).
Returns
-------
float
Value of rate coefficient at `T`.
"""
R = 8.314 # J/(mol路K)
return k0 * np.exp(-Ea / R * (1/(273.15 + T) - 1/(273.15 + T0)))
We can now implement a function to calculate the reaction rates for all species and moments.
def model_rates(C: np.ndarray, T: float, params: dict) -> np.ndarray:
"""Reaction rates.
Parameters
----------
C : np.ndarray
Concentration vector (mol/L).
T : float
Temperature (掳C).
params : dict
Model parameters.
Returns
-------
np.ndarray
Reaction rate vector (mol/(L路s)).
"""
# Rate coefficients
kinetics = params['kinetics']
kd = arrhenius(T, **kinetics['kd'])
kp = arrhenius(T, **kinetics['kp'])
kt = arrhenius(T, **kinetics['kt'])
f = kinetics['f']
alpha = kinetics['alpha']
ktd = kt*alpha
ktc = kt*(1 - alpha)
# Unpack concentrations
I = C[Ix.I]
M = C[Ix.M]
mu = C[Ix.mu0:Ix.mu2+1]
# Reaction rates
r = np.empty_like(C)
ri = 2*f*kd*I
r[Ix.I] = -kd*I
r[Ix.M] = -(ri + kp*M*mu[0])
r[Ix.S] = 0.0
r[Ix.mu0] = ri - 2*kt*mu[0]**2
r[Ix.mu1] = ri + kp*M*mu[0] - 2*kt*mu[0]*mu[1]
r[Ix.mu2] = ri + kp*M*(mu[0] + 2*mu[1]) - 2*kt*mu[0]*mu[2]
r[Ix.lbd0] = (ktc + 2*ktd)*mu[0]**2
r[Ix.lbd1] = 2*kt*mu[0]*mu[1]
r[Ix.lbd2] = 2*kt*mu[0]*mu[2] + 2*ktc*mu[1]**2
return r
Next, we implement a function to evaluate all auxiliary variables that we wish to track over time but that are not state variables. This function will be called both during and after the integration.
def model_aux(t: float, n: np.ndarray, T: float, params: dict) -> tuple:
"""Auxiliary model function.
Parameters
----------
t : float
Time (s).
n : np.ndarray
Components moles vector (mol).
T : float
Temperature (掳C).
params : dict
Model parameters.
Returns
-------
tuple
Auxiliary variables.
"""
# Unpack the moles vector (mol)
n_I = n[Ix.I]
n_M = n[Ix.M]
n_S = n[Ix.S]
n_mu1 = n[Ix.mu1]
n_ld1 = n[Ix.lbd1]
# Total mass of the reaction mixture (kg)
MW = params['MW']
m_t = n_I*MW['I'] + (n_M + n_mu1 + n_ld1)*MW['M'] + n_S*MW['S']
# Total volume of the reaction mixture (L)
rho = params['rho']
V = n_M*(MW['M']/rho['M']) + n_S*(MW['S']/rho['S']) \
+ (n_mu1 + n_ld1)*(MW['M']/rho['P'])
# Molar concentrations (mol/L)
C = n/(V + 1e-10)
# Reaction rates (L/(mol路s))
r = model_rates(C, T, params)
# Component molar feed rates (mol/s)
feed = params['feed']
ndot_in = np.zeros_like(n)
ndot_in[Ix.I] = feed['I'](t)
ndot_in[Ix.M] = feed['M'](t)
ndot_in[Ix.S] = feed['S'](t)
# Total mass feed rate (kg/s)
mdot_in = ndot_in[Ix.I]*MW['I'] + \
ndot_in[Ix.M]*MW['M'] + ndot_in[Ix.S]*MW['S']
# Temperature control
KP = params['control']['KP']
T_sp = params['control']['T_sp'](t)
Tj_min = params['control']['Tj_min']
Tj_max = params['control']['Tj_max']
Tj = T_sp - KP*(T - T_sp)
Tj = np.clip(Tj, Tj_min, Tj_max)
# MTSR
DHpol = params['heatbal']['DHpol']
cp = params['heatbal']['cp']
MTSR = T + n_M*(-DHpol)/(m_t*cp + 1e-10)
return m_t, V, C, r, ndot_in, mdot_in, Tj, MTSR
All parts are now in place to implement the derivative of the state vector. This function is relatively straightforward, as most of the heavy lifting is handled by the previously defined functions.
def model_xdot(t: float, x: np.ndarray, params: dict) -> np.ndarray:
"""Calculate the derivative of the state vector, dx/dt.
x = [n(I), n(M), n(S), n(mu_0)..n(mu_2), n(lambda_0)..n(lambda_2), T]
Parameters
----------
t : float
Time (s).
x : np.ndarray
State vector.
params : dict
Model parameters.
Returns
-------
np.ndarray
Time derivative of the state vector.
"""
# Unpack the state vector
n = x[:Ix.T]
T = x[Ix.T]
# Allocate the state vector derivative and respective views
xdot = np.zeros_like(x)
ndot = xdot[:Ix.T]
Tdot = xdot[Ix.T:]
# Evaluate the auxiliary variables
m_t, V, _, r, ndot_in, mdot_in, Tj, _ = model_aux(t, n, T, params)
# Material balances
ndot[:] = ndot_in + r*V
# Heat balance
# This could also be moved to own function or model_aux
heatbal = params['heatbal']
cp = heatbal['cp']
DHpol = heatbal['DHpol']
Tin = heatbal['Tin']
D = heatbal['D']
U = heatbal['U']
qdot_in = mdot_in*cp*(Tin - T)
qdot_rxn = DHpol*r[1]*V
A = 4*((V*1e-3)/D)
qdot_j = U*A*(T - Tj)
Tdot[:] = (qdot_in + qdot_rxn - qdot_j)/(m_t*cp + 1e-10)
return xdot
Finally, we perform the numerical integration using a suitable ODE solver. This system is very stiff, therefore we need to choose an implicit solver like LSODA.
def solve_model(n0: np.ndarray,
T0: float,
params: dict,
tend: float,
) -> tuple[np.ndarray, ...]:
"""Solve the dynamic model.
Parameters
----------
n0 : float
Initial component moles (mol).
T0 : float
Initial reactor temperature (掳C).
params : dict
Model parameters.
tend : float
Final time.
Returns
-------
tuple[np.ndarray, ...]
t, n, T, m_t, V, C, ndot_in, mdot_in, Tj, MTSR.
"""
# Initial conditions
x0 = np.hstack([n0, T0])
solution = solve_ivp(model_xdot,
t_span=(0.0, tend),
y0=x0,
args=(params,),
method='LSODA',
rtol=1e-4,
atol=1e-10,
max_step=600.0)
# Unpack the solution
t = solution.t
x = solution.y
n = x[:Ix.T, :]
T = x[Ix.T, :]
# Compute aux variables
m_t = np.empty_like(t)
V = np.empty_like(t)
C = np.empty_like(n)
ndot_in = np.empty_like(n)
mdot_in = np.empty_like(t)
Tj = np.empty_like(t)
MTSR = np.empty_like(t)
for i in range(t.size):
m_t[i], V[i], C[:, i], _, ndot_in[:, i], mdot_in[i], Tj[i], MTSR[i] = \
model_aux(t[i], n[:, i], T[i], params)
return t, n, T, m_t, V, C, ndot_in, mdot_in, Tj, MTSR
6.3. 鈻讹笍 Run Simulation#
We will use some educated guesses for an imaginary solution homopolymerization. Feel free to change them.
params = {
'MW': {
'I': 0.150, # kg/mol
'M': 0.100, # kg/mol
'S': 0.100, # kg/mol
},
'rho': {
'M': 0.90, # kg/L
'S': 0.80, # kg/L
'P': 1.10, # kg/L
},
'kinetics': {
'f': 0.5,
'alpha': 0.5, # ktd/(ktd + ktc)
'kd': {
'k0': 5e-4, # 1/s
'Ea': 50e3, # J/mol
'T0': 80. # 掳C
},
'kp': {
'k0': 2e3, # L/(mol路s)
'Ea': 20e3, # J/mol
'T0': 80. # 掳C
},
'kt': {
'k0': 5e7, # L/(mol路s)
'Ea': 10e3, # J/mol
'T0': 80. # 掳C
},
},
'control': {
'KP': 10.,
'Tj_min': 35., # 掳C
'Tj_max': 150., # 掳C
'T_sp': lambda t: np.interp(t, [0., 1800., 10*3600.], [20., 80., 80.])
},
'feed': {
'I': lambda t: 2e-2 if (t > 1800. and t < 8*3600.) else 0., # mol/s
'M': lambda t: 1.5e0 if (t > 1800. and t < 6*3600.) else 0., # mol/s
'S': lambda t: 2e0 if (t > 1800. and t < 6*3600.) else 0., # mol/s
},
'heatbal': {
'Tin': 25., # 掳C
'D': 2., # m
'U': 400., # W/(m虏路K)
'cp': 2100., # J/(kg路K)
'DHpol': -80e3 # J/mol
}
}
# Initial reactor content
n0 = np.zeros(len(Ix) - 1)
n0[Ix.I] = 0.0
n0[Ix.M] = 0.0
n0[Ix.S] = 1e4
# Initial temperature
T0 = params['control']['T_sp'](0.0)
# Simulation time
tend = 10*3600
Finally, we perform the integration; this should not take more than a fraction of a second, since there are only 10 states.
t, n, T, m_t, V, C, ndot_in, mdot_in, Tj, MTSR = solve_model(n0, T0, params, tend)
6.4. 馃搳 Plots#
Now, let鈥檚 visualize the model results.
6.4.1. Reactor Contents#
fig1, ax = plt.subplots(4, 1, sharex=True, figsize=(8, 8))
fig1.suptitle("Reactor Contents")
fig1.tight_layout()
fig1.align_ylabels()
t_hour = t/3600
# Feed I
ax[0].plot(t_hour, ndot_in[Ix.I, :], label="I")
ax[0].set_ylabel("Flow (mol/s)")
ax[0].grid(True)
ax[0].legend(loc='upper right')
# Feed others
ax[1].plot(t_hour, ndot_in[Ix.M, :], label="M")
ax[1].plot(t_hour, ndot_in[Ix.S, :], label="S")
ax[1].set_ylabel("Flow (mol/s)")
ax[1].grid(True)
ax[1].legend(loc='upper right')
# Volume
ax[2].plot(t_hour, V)
ax[2].set_ylabel("Volume (L)")
ax[2].set_ylim(0.0, None)
ax[2].grid(True)
# Mass
ax[3].plot(t_hour, m_t)
ax[3].set_ylabel("Mass (kg)")
ax[3].set_ylim(0.0, None)
ax[3].grid(True)
ax[-1].set_xlabel("Time (h)")
Text(0.5, 58.7222222222222, 'Time (h)')
6.4.2. Species Amounts & Concentrations#
fig2, ax = plt.subplots(3, 1, sharex=True, figsize=(8, 8))
fig2.suptitle("Species Amounts & Concentrations")
fig2.tight_layout()
fig2.align_ylabels()
# Number of moles
ax[0].plot(t_hour, n[Ix.I, :], label="I")
ax[0].plot(t_hour, n[Ix.M, :], label="M")
ax[0].plot(t_hour, n[Ix.S, :], label="S")
ax[0].plot(t_hour, n[Ix.mu1, :] + n[Ix.lbd1, :], label="P")
ax[0].set_ylabel("Moles (mol)")
ax[0].grid(True)
ax[0].legend(loc='upper right')
# Molar concentration I
ax[1].plot(t_hour, C[Ix.I, :], label="I")
ax[1].set_ylabel("Concentration (mol/L)")
ax[1].grid(True)
ax[1].legend(loc='upper right')
# Molar concentrations others
ax[2].plot(t_hour, C[Ix.M, :], label="M")
ax[2].plot(t_hour, C[Ix.S, :], label="S")
ax[2].plot(t_hour, C[Ix.mu1, :] + C[Ix.lbd1, :], label="P")
ax[2].set_ylabel("Concentration (mol/L)")
ax[2].grid(True)
ax[2].legend(loc='upper right')
ax[-1].set_xlabel("Time (h)")
Text(0.5, 58.7222222222222, 'Time (h)')
6.4.3. Polymer Chain Length and Dispersity#
fig3, ax = plt.subplots(2, 1, sharex=True, figsize=(8, 6))
fig3.suptitle("Polymer Chain Length and Dispersity")
fig3.tight_layout()
fig3.align_ylabels()
# DPn and DPw
DPn = (n[Ix.mu1, :] + n[Ix.lbd1, :])/(n[Ix.mu0, :] + n[Ix.lbd0, :] + 1e-10)
DPw = (n[Ix.mu2, :] + n[Ix.lbd2, :])/(n[Ix.mu1, :] + n[Ix.lbd1, :] + 1e-10)
ax[0].plot(t_hour, DPn, label="DPn")
ax[0].plot(t_hour, DPw, label="DPw")
ax[0].set_ylabel("Degree of polymerization")
ax[0].grid(True)
ax[0].legend(loc='best')
# PDI
PDI = DPw/(DPn + 1e-10)
ax[1].plot(t_hour, PDI)
ax[1].set_ylabel("PDI")
ax[1].set_ylim(1.0, None)
ax[1].grid(True)
ax[-1].set_xlabel("Time (h)")
Text(0.5, 36.72222222222221, 'Time (h)')
6.4.4. Temperature & Safety#
fig4, ax = plt.subplots(1, 1)
fig4.suptitle("Temperature & Safety")
fig4.tight_layout()
# Temperatures
ax.plot(t_hour, T, label="Reactor")
ax.plot(t_hour, Tj, label="Jacket")
ax.plot(t_hour, MTSR, label="MTSR")
ax.set_ylabel("Temperature (掳C)")
ax.grid(True)
ax.legend()
ax.set_xlabel("Time (h)")
Text(0.5, 23.52222222222222, 'Time (h)')
6.5. 馃攷 Questions#
Why is there a maximum in \([I]\) if the initiator feed rate is actually constant?
What happens to \(PDI\) when the parameter
alphais set to 0 or 1? Explain why this occurs.Why is the MTSR highest at the start of the feed period? What can be done to reduce it?
Why does the jacket temperature suddenly drop at 6 hours?
What changes occur in the jacket temperature between approximately 1 hour and 6 hours?
The current heat transfer model assumes a constant overall heat transfer coefficient, \(U\). What are your thoughts on this assumption?
Can you propose an empirical or semi-theoretical expression to describe the evolution of \(U\) during the run?
Implement the model proposed in Question 7.