Copolymerization-Reaction-Modeling README and Notes

Source: gtancev/Copolymerization-Reaction-Modeling
Author: Georgi Tancev, PhD
Relevance to project: ⭐⭐⭐⭐⭐ — Core model for radical copolymerization simulations

What This Code Does

Simulates batch bulk radical copolymerization of two monomers (A and B) in a 15 L reactor. The system of ODEs tracks concentrations of initiator, both monomers, and both radical species over time (~2 h).

Inferred Monomer System (from hardcoded parameters)

Parameter	Monomer A	Monomer B
Molar mass	104 g/mol	100 g/mol
Density	0.90 kg/L	0.94 kg/L
kp₀	410 L/(mol·s)	930 L/(mol·s)
kt₀	2.4×10⁷ L/(mol·s)	9.2×10⁶ L/(mol·s)
Reactivity ratio	rA = 0.52	rB = 0.46

This matches styrene (A) / methyl methacrylate (B) nearly exactly:

Styrene MW = 104, rA ≈ 0.52; MMA MW = 100.12, rB ≈ 0.46 (literature values)
Initiator MI = 164 g/mol → AIBN (MW = 164.21), kd = 6.77×10⁻⁶ s⁻¹ (≈60°C)

Model Components

1. Kinetic Mechanism (`batch.m`)

Uses the terminal model (Mayo-Lewis). The ODE system:

dcI — initiator decomposition (1st order, kd)
dcA, dcB — monomer consumption via homo- and cross-propagation
dcA_R, dcB_R — radical balance for A-centered and B-centered radicals (quasi-steady-state not assumed — radicals are integrated explicitly)

Cross-propagation from reactivity ratio definitions:

kp_AB = kp_AA / rA
kp_BA = kp_BB / rB

Cross-termination by geometric mean:

kt_AB = sqrt(kt_AA * kt_BB)

2. Gel/Glass Effect

Rate constants are modified by polymer weight fraction (wp) via an empirical diffusion-limited formula:

kp = 1 / (1/kp0 + exp(C_eta * wp) / kp_D)
kt = 1 / (1/kt0 + exp(C_eta * wp) / kt_D) + C_RD * kp * (1 - wp)

Where C_eta = 25, C_RD = 180. This captures the Trommsdorff (gel) effect and glass effect empirically. More rigorous alternatives use the Vrentas-Duda free volume theory (see Ref. [1] below).

3. Post-processing (`cumulative.m`)

Computes using the method of moments:

Instantaneous copolymer composition (F_A) — Mayo-Lewis equation
Cumulative composition (F_A_c) via numerical integration (trapz)
Number- and weight-average chain lengths (n_N, n_W) and PDI
Termination split between combination (ktc) and disproportionation (ktd), with ratio C_t = 1000 (i.e., predominantly disproportionation)

Bugs and Code Quality Issues

Parameter duplication — rate constants, densities, and molar masses are hardcoded independently in batch.m, cumulative.m, and main.m. Changing a parameter in one place does not update the others. Any adaptation should centralize parameters.
No Arrhenius temperature dependence — all rate constants are fixed (implicitly ~60°C). Adding k = A * exp(-Ea/RT) would make the model temperature-aware.
f = 0.5 is fixed — initiator efficiency should ideally be conversion-dependent (decreases at high conversion due to cage effect).
No chain transfer — to monomer (Cm), to solvent, or to chain transfer agent (CTA). This limits applicability to systems without deliberate MWD control.
trapz integration in cumulative.m is first-order accurate and sensitive to the ODE time step. Consider switching to cumulative Simpson's rule or an analytic running integral.
No semibatch/continuous mode — the model is batch-only. Semibatch addition of monomer is a common industrial process for composition control.

Physical Limitations

Terminal model only — penultimate unit effect is not included. For rA×rB << 1 or rA×rB >> 1 (highly alternating or blocky), the penultimate model may be needed.
Bulk/solution only — no heterogeneous phase (see the companion emulsion repo for that).
No diffusion-controlled propagation at high conversion — at very high wp (glass effect), kp can also become diffusion-limited, but this requires the kp_D term and glass-effect parameters to be carefully tuned.
Quasi-bulk initiator decomposition — f is assumed constant; in real systems it can drop from ~0.5 to ~0.1 at high conversion.

Relevance and Adaptation Suggestions for This Project

High Relevance

This is the most directly applicable model for radical copolymerization simulation. The Mayo-Lewis formalism with gel effect is standard for batch bulk/solution copolymerization of styrene, acrylates, methacrylates.

Suggested Adaptations

Python port — convert to Python with scipy.integrate.solve_ivp (use method='Radau' or 'BDF' as ode15s equivalent). Example structure:

from scipy.integrate import solve_ivp
import numpy as np

def batch(t, c, params):
    cI, cA, cA_R, cB, cB_R = c
    # ... kinetics with params dict
    return [dcI, dcA, dcA_R, dcB, dcB_R]

sol = solve_ivp(batch, [0, 7200], c0, method='Radau', args=(params,), dense_output=True)

Parameterize for other monomer pairs — the structure is general. By changing rA, rB, kp_AA, kp_BB, kt_AA, kt_BB, MA, MB you can simulate any comonomer pair (e.g., styrene/BA, MMA/BA, VAc/ethylene).
Add chain transfer — add dcA -= ktr_A * cA * cR terms and a transfer constant Cm to control Mn independently of conversion.
Add temperature as a variable — couple to an energy balance ODE for non-isothermal simulation.
Connect to emulsion model — the batch.m kinetics can be embedded inside the particle phase of the emulsion model to simulate emulsion copolymerization.
Validate against mcPolymer — the mcPolymer Monte Carlo simulator in the mcPolymer/ directory can provide reference MWDs and composition distributions to validate this deterministic model.

Key Literature

[1] Modeling Insight into the Diffusion-Limited Cause of the Gel Effect in Free Radical Polymerization — G.A. O'Neil et al., Macromolecules 1999. Vrentas-Duda free volume theory for gel effect; more rigorous than the empirical exponential form used here.
[2] Free radical polymerizations associated with the Trommsdorff effect under semibatch reactor conditions — Ray et al., Polym. Eng. Sci. 1995. Gel and glass effects in batch and semibatch, MMA/PS systems.
[3] Pseudo-Homopolymerization Approach To Predict the MWD in Copolymerization — Zapata-González et al., I&EC Res. 2018. More efficient approach to MWD in copolymerization using the PHP+RQSSA method.
[4] Odian, G. Principles of Polymerization, 4th ed., Wiley 2004. Chapters 6–7 for copolymerization fundamentals and kinetics.
[5] Dotson, N.A. et al. Polymerization Process Modeling, VCH 1996. Standard reference for the moment method used in cumulative.m.

Notes written: 2026-06-26

7.1 KiB Raw Blame History Unescape Escape