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}^{\mathrm{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}_i\)), and dead polymer (\(P_i\)).
State variable choice
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 \(m \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_{\mathrm{t}} = k_{\mathrm{tc}} + k_{\mathrm{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_{\mathrm{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}_{\mathrm{in}}\) is the total inlet mass flowrate, \(T_{\mathrm{in}}\) is the mean feed temperature, \(\Delta H_{\mathrm{pol}}\) is the polymerization enthalpy, \(U\) is the overall heat transfer coefficient, \(A\) is the heat transfer area, and \(T_{\mathrm{j}}\) is the jacket temperature. We assume the jacket temperature is automatically adjusted by means of a simple proportional-only controller:
where \(T_{\mathrm{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’s write some code to solve this semibatch homopolymerization model.
6.2. 🧮 Model Implementation#
6.2.1. Architecture#
This model is more complex than the preceding ones, so it is essential to organize the code in a modular way. The diagram below provides a high-level overview of the model architecture: solve_model advances the system by calling SciPy’s solve_ivp, which repeatedly evaluates model_xdot to compute the state derivatives.
Within each evaluation, model_xdot calls model_aux, which delegates the calculation of all kinetic contributions to model_rates, with rate constants obtained from arrhenius. After the integration is complete, solve_model calls model_aux once more to reconstruct the time evolution of the auxiliary variables using the computed state-variable solution as input.
graph TD
solve_model[solve_model]
solve_ivp[solve_ivp]
model_xdot[model_xdot]
model_aux[model_aux]
model_rates[model_rates]
solve_model --> solve_ivp
solve_model --> model_aux
solve_ivp --> model_xdot
model_xdot --> model_aux
model_aux --> model_rates
model_rates --> arrhenius
6.2.2. Code 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
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 formation rates for all species and moments.
def model_rates(C: np.ndarray, T: float, params: dict) -> np.ndarray:
"""Component formation rates.
Parameters
----------
C : np.ndarray
Concentration vector (mol/L).
T : float
Temperature (°C).
params : dict
Model parameters.
Returns
-------
np.ndarray
Component formation 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)
# Concentrations
I = C[Ix.I]
M = C[Ix.M]
mu = C[Ix.mu0:Ix.mu2+1]
# Formation 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, recipe: dict, params: dict
) -> tuple:
"""Auxiliary model function.
Parameters
----------
t : float
Time (s).
n : np.ndarray
Components moles vector (mol).
T : float
Temperature (°C).
recipe : dict
Polymerization recipe.
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)
# Component formation rates (mol/(L·s))
r = model_rates(C, T, params)
# Component molar feed rates (mol/s)
feed = recipe['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']
Tj_min = params['control']['Tj_min']
Tj_max = params['control']['Tj_max']
T_sp = recipe['T_sp'](t)
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, recipe: dict, 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.
recipe : dict
Polymerization recipe.
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, recipe, params)
# Material balances
ndot[:] = ndot_in + r*V
# Heat balance
# This could also be moved to own function or model_aux
cp = params['heatbal']['cp']
DHpol = params['heatbal']['DHpol']
D = params['heatbal']['D']
U = params['heatbal']['U']
Tin = recipe['feed']['Tin']
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(recipe: dict, params: dict) -> tuple[np.ndarray, ...]:
"""Solve the dynamic model.
Parameters
----------
recipe : dict
Polymerization recipe.
params : dict
Model parameters.
Returns
-------
tuple[np.ndarray, ...]
t, n, T, m_t, V, C, ndot_in, mdot_in, Tj, MTSR.
"""
# Initial conditions
n0 = recipe['n0']
T0 = recipe['T0']
x0 = np.hstack([n0, T0])
solution = solve_ivp(model_xdot,
t_span=(0.0, recipe['tend']),
y0=x0,
args=(recipe, 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], recipe, params)
return t, n, T, m_t, V, C, ndot_in, mdot_in, Tj, MTSR
6.3. 📊 Results and Visualization#
6.3.1. Runing the Simulation#
We will use some educated guesses for an imaginary solution homopolymerization. Feel free to change them.
# Model parameters
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': 100., # °C
},
'heatbal': {
'D': 2., # m
'U': 400., # W/(m²·K)
'cp': 2100., # J/(kg·K)
'DHpol': -80e3 # J/mol
}
}
# Polymerization recipe
recipe = {
'tend': 10*3600., # s
'n0': np.zeros(len(Ix) - 1),
'T0': 20.0, # °C
'T_sp': lambda t: np.interp(t, [0., 1800., 10*3600.], [20., 80., 80.]),
'feed': {
'Tin': 25.0, # °C
'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
}
}
recipe['n0'][Ix.S] = 1e4
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(recipe, params)
Now, let’s visualize the model results.
6.3.2. Reactor Contents#
fig1, ax = plt.subplots(4, 1, sharex=True, figsize=(7, 4*2))
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)");
6.3.3. Species Amounts & Concentrations#
fig2, ax = plt.subplots(3, 1, sharex=True, figsize=(7, 3*2))
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)");
6.3.4. Polymer Chain Length and Dispersity#
fig3, ax = plt.subplots(2, 1, sharex=True, figsize=(7, 2*2))
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)");
6.3.5. Temperature & Safety#
fig4, ax = plt.subplots(figsize=(7, 5))
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)");
6.4. 🔎 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.