Initial commit

NOTES.md

@@ -0,0 +1,187 @@
+# Emulsion-Polymerization-Modeling — Analysis Notes
+
+**Source:** [gtancev/Emulsion-Polymerization-Modeling](https://github.com/gtancev/Emulsion-Polymerization-Modeling)  
+**Author:** Georgi Tancev, PhD  
+**Relevance to project:** ⭐⭐⭐⭐ — Heterogeneous-phase framework; pairs with the copolymerization model
+
+---
+
+## What This Code Does
+
+Simulates **seeded emulsion polymerization** (Stage II — monomer droplets still present) in a 10 L batch reactor. Two model variants are run:
+
+- **Model 1 (`emulsion1.m`)** — Fixed n̄ = 0.5 (Smith-Ewart Case II)
+- **Model 2 (`emulsion2.m`)** — Dynamic n̄ from radical entry/exit balance
+
+The state variables are:
+- `n_M` — total moles of monomer (both phases)
+- `n_M_p` — moles of monomer inside particles
+
+---
+
+## Inferred Monomer/System (from parameters)
+
+| Parameter | Value | Interpretation |
+|-----------|-------|----------------|
+| `k_p = 715` L/(mol·s) | — | Consistent with **styrene** at ~60°C (lit. ~480–715) |
+| `M_M = 0.1` kg/mol | — | Styrene MW = 104, close approximation |
+| `M_I = 0.164` kg/mol | — | **AIBN** (MW = 164.21) |
+| `k_d = 5.55e-6` s⁻¹ | — | AIBN at ~60°C |
+| `d_p_0 = 50 nm` | — | Typical seed particle diameter |
+| `phi_M = 0.5` | — | Monomer volume fraction in swollen particle |
+
+The system is a **styrene-like seeded emulsion polymerization** in water with AIBN initiator and a nonionic-like emulsifier (Langmuir adsorption with K_E = 100 L/mol).
+
+---
+
+## Model Components
+
+### Phase Structure
+
+Three phases are tracked implicitly:
+1. **Water phase** — initiator, emulsifier, free radicals (c_R computed at quasi-steady-state: `c_R = sqrt(2f·kd·cI/kt)`)
+2. **Particle phase** — polymer + swollen monomer (volume tracked via v_pol and v_p)
+3. **Droplet phase** — reservoir of monomer (present while `n_M_d > 0`)
+
+### Monomer Partitioning
+
+While monomer droplets exist (Stage II), monomer in particles is maintained at the thermodynamic swelling equilibrium:
+```
+dn_M_P = n_s * c_M * 1/(1-phi) * dv_pol
+```
+Once droplets are exhausted, `dn_M_P = dn_M` (Stage III, all remaining monomer is in particles).
+
+### Model 1 — Fixed n̄ (`emulsion1.m`)
+
+Smith-Ewart **Case II**: n̄ = 0.5 (at any instant, each particle contains exactly 0 or 1 radical with equal probability). Valid for small particles where termination upon entry of a second radical is instantaneous.
+
+```matlab
+n_bar = 0.5;
+dn_M = -k_p * c_M * n_bar / N_A * n_s;
+```
+
+### Model 2 — Dynamic n̄ (`emulsion2.m`)
+
+Steady-state radical balance inside particles:
+```
+n̄ = ρ / (2ρ + k_des)
+```
+where:
+- `ρ = A_p * k_e * c_R * N_A` — radical capture rate (1/s), proportional to particle surface area
+- `k_des` — radical desorption rate constant (1/s)
+
+This is a simplified **Nomura/Ugelstad** expression that includes desorption as an exit mechanism. It predicts n̄ < 0.5 when desorption is significant.
+
+### Emulsifier Partitioning (in `main.m` only)
+
+The code computes min/max bounds on emulsifier amount using a **Langmuir adsorption isotherm**:
+```
+F_E = F_E_inf * K_E * E_w / (1 + K_E * E_w)
+```
+where F_E_inf is the maximum surface coverage (from `a_E`, the area per emulsifier molecule). This is used to check whether the emulsifier loading falls between minimum stabilization (20% surface coverage) and CMC saturation. **This is only calculated, not used dynamically in the ODE.**
+
+---
+
+## Critical Bug: Typo in Function Names
+
+**The function files contain a typo that will cause MATLAB to fail:**
+
+- File: `emulsion1.m` → Function declared as: `function [ dm ] = emuslion1(...)` (missing 'l', 's' transposed)
+- File: `emulsion2.m` → Function declared as: `function [ dm ] = emuslion2(...)` (same typo)
+
+In MATLAB, the function name inside the file **must match the filename**. `main.m` calls `@emulsion1` and `@emulsion2` (correct names), but the functions are defined as `emuslion1` and `emuslion2`. This will throw an error on some MATLAB versions. **Fix by correcting the function declarations:**
+
+```matlab
+% emulsion1.m line 1: change
+function [ dm ] = emuslion1( t,n,n_s,n_M_0,V_p_0 )
+% to:
+function [ dm ] = emulsion1( t,n,n_s,n_M_0,V_p_0 )
+```
+
+Same fix needed in `emulsion2.m`.
+
+---
+
+## Physical Limitations
+
+1. **Seeded emulsion only (no nucleation)** — particle number `n_s` is fixed. No homogeneous, micellar, or coagulative nucleation is modeled. This restricts simulation to Stage II and III only.
+
+2. **Monodisperse particles** — all particles are assumed identical (same size, same n̄). Real emulsions have a particle size distribution (PSD).
+
+3. **Fixed radical concentration in water** — `c_R` is computed once from the initial initiator concentration and held constant. In reality, `cI` decays and `c_R` drops with it; this is a significant approximation for long reaction times.
+
+4. **No temperature model** — all rate constants are isothermal. The polymerization of styrene is exothermic (~70 kJ/mol); temperature rise would accelerate kd and kp, changing kinetics.
+
+5. **No gel/glass effect** — unlike the copolymerization model, emulsion1/2 use constant `k_p`. At high conversion inside particles (phi → 0), diffusion limitation on kp should be included.
+
+6. **No coagulation or Ostwald ripening** — particle stability and size evolution by mass transfer are not modeled.
+
+7. **n̄ expression in Model 2 is simplified** — the full Smith-Ewart recursion for n̄ accounts for termination, entry, and exit across all radical numbers (n = 0, 1, 2, ...). The simplified expression `ρ/(2ρ + k_des)` holds in the limit of small n̄ and pseudo-bulk behavior but is not exact for intermediate cases.
+
+---
+
+## Outputs Computed in `main.m`
+
+For both models, the following are plotted vs. conversion or time:
+1. Conversion vs. time
+2. Total particle volume and polymer volume vs. conversion
+3. Monomer volume fraction in particles (φ) vs. conversion
+4. Particle diameter vs. conversion
+5. n̄ vs. time/conversion
+6. Number-average (n_N) and weight-average (n_W) chain lengths
+
+The chain lengths are computed from the ratio of characteristic time scales:
+- `tau_p = 1/(kp * c_M)` — monomer consumption time
+- `tau_in = 1/ρ` — radical entry interval
+
+---
+
+## Adaptation Suggestions for This Project
+
+### 1. Fix the function name typo (immediate)
+See above — must match file name and function declaration.
+
+### 2. Fix the decaying initiator (important)
+Replace the fixed `c_R` with an ODE for initiator:
+```matlab
+% Add to state vector:
+dcI = -kd * cI;
+c_R = sqrt(2*f*kd*cI/kt);  % updated at each step
+```
+
+### 3. Couple with copolymerization kinetics
+The most powerful adaptation: replace the single-monomer propagation (`k_p * c_M * n_bar`) with the Mayo-Lewis terminal model (from `Copolymerization-Reaction-Modeling/batch.m`) operating inside the particle phase. This would give a full **emulsion copolymerization** model.
+
+### 4. Python port
+Convert to Python using `scipy.integrate.solve_ivp`. The particle volume calculations and phase-switching logic (Stage II → III) translate straightforwardly.
+
+```python
+def emulsion_odes(t, n, params):
+    n_M, n_M_p = n
+    # compute phase volumes, phi, c_M, n_bar
+    # return [dn_M, dn_M_p]
+```
+
+### 5. Add particle size distribution (PSD)
+Replace the single particle with a population balance — discretize particle sizes and track radical occupancy per size class.
+
+### 6. Compare n̄ models
+The two models (fixed vs. dynamic n̄) differ in their treatment of exit kinetics. Comparing their predictions vs. conversion reveals the importance of desorption for the specific monomer system. Small, water-soluble radicals (e.g., from vinyl acetate) have large k_des; hydrophobic monomers (styrene) have small k_des, so Model 1 and Model 2 should agree well for styrene.
+
+---
+
+## Key Literature
+
+- [1] [Free Radical Entry in Emulsion Polymerization: A Comprehensive Theory](https://consensus.app/papers/details/8827fafe862952529adba620738c1138/?utm_source=claude_code) — Herrera-Ordonez, *Macromol. React. Eng.* 2025. Critical review of radical entry models; shows diffusion/adsorption is rate-limiting (not aqueous-phase propagation as in older theories).
+
+- [2] [Calculations of the Average Number of Radicals per Particle](https://consensus.app/papers/details/cf3fcb89b1115796a8c5b3c05c827b09/?utm_source=claude_code) — Tiwari, 2021. Smith-Ewart recursion and pseudo-bulk vs. 0-1 kinetics — directly relevant to the n̄ approximation used here.
+
+- [3] [Particle by Particle Kinetic Monte Carlo Tracking in Miniemulsion](https://consensus.app/papers/details/b23fc44c9c905b37830dc7d44e9c1fc9/?utm_source=claude_code) — Marien et al., *Macromolecules* 2019. kMC approach showing the limits of the average-particle (Smith-Ewart) assumption; important context for Model 1 vs. Model 2.
+
+- [4] Gilbert, R.G. *Emulsion Polymerization: A Mechanistic Approach*, Academic Press, 1995. The canonical reference for emulsion polymerization kinetics.
+
+- [5] Nomura, M.; Harada, M. *J. Appl. Polym. Sci.* 1981, 26, 17. Original simplified n̄ expression similar to what is used in Model 2.
+
+---
+
+*Notes written: 2026-06-26*

emulsion1.m

@@ -0,0 +1,45 @@
+function [ dm ] = emuslion1( t,n,n_s,n_M_0,V_p_0 )
+
+k_p = 715; % L/(mol*s)
+n_bar = 0.5;
+N_A = 6.022*10^23; % 1/mol
+M_M = 0.1; % kg/mol
+p_M = 0.94; % kg/L
+p_P = 1.1; % kg/L
+
+%%
+
+n_M = n(1); % amount of total monomer available
+n_M_p = n(2); % amount of monomer in particles
+n_M_d = n_M-n_M_p; % amount of monomer in dispersed phase
+
+%%
+
+v_pol = V_p_0+(n_M_0-n_M)*M_M/(n_s*p_P);
+v_p = v_pol+(n_M_p/n_s)*(M_M/p_M);
+phi = (v_p-v_pol)/v_p;
+c_M = n_M_p/(n_s*v_p);
+
+%%
+
+dn_M = - k_p*c_M*n_bar/N_A*n_s;
+dv_pol = - dn_M*M_M/(p_P*n_s);
+
+%%
+    
+    if n_M_d > 0
+        
+        dn_M_P = n_s.*c_M.*1./(1-phi).*dv_pol;
+    
+    else
+
+        dn_M_P = dn_M;
+        
+    end
+    
+%%
+
+dm = [dn_M; dn_M_P];
+
+end
+

emulsion2.m

@@ -0,0 +1,60 @@
+function [ dm ] = emuslion2( t,n,n_s,n_M_0,V_p_0,V )
+
+k_p = 715; % L/(mol*s)
+k_des = 1e-1; % 1/s
+k_e = 1e-6; % dm/s
+k_t = 9.8e6; % L/(mol*s)
+k_d = 5.55e-6; % 1/s
+
+N_A = 6.022*10^23; % 1/mol
+M_M = 0.1; % kg/mol
+p_M = 0.94; % kg/L
+p_P = 1.10; % kg/L
+w_I = 0.01;
+m_M_0 = n_M_0*M_M; % kg
+m_I = w_I*m_M_0; % kg
+M_I = 0.164; % kg/mol
+f = 0.5;
+c_I = m_I/(M_I*V); % mol/L
+c_R = sqrt(2*f*k_d*c_I/k_t); % mol/L
+
+%%
+
+n_M = n(1); % amount of total monomer available
+n_M_p = n(2); % amount of monomer in particles
+n_M_d = n_M-n_M_p; % amount of monomer in dispersed phase
+
+%%
+
+v_pol = V_p_0+(n_M_0-n_M)*M_M/(n_s*p_P); 
+v_p = v_pol+(n_M_p/n_s)*(M_M/p_M);
+phi = (v_p-v_pol)/v_p;
+c_M = n_M_p/(n_s*v_p);
+
+diam = (6*v_p/pi)^(1/3);
+A_p = 4*pi*(diam/2)^2;
+rho = A_p*k_e*c_R*N_A;
+n_bar = (rho/(2*rho+k_des));
+
+%%
+
+dn_M = - k_p*c_M*n_bar/N_A*n_s;
+dv_pol = - dn_M*M_M/(p_P*n_s);
+
+%%
+    
+    if n_M_d > 0
+        
+        dn_M_P = n_s*c_M*1/(1-phi)*dv_pol;
+    
+    else
+
+        dn_M_P = dn_M;
+        
+    end
+    
+%%
+
+dm = [dn_M; dn_M_P];
+
+end

main.m

@@ -0,0 +1,213 @@
+%% Polymerization Reaction & Colloid Engineering
+% George Tancev
+
+clear all; close all; clc;
+options = optimset('Display','none','MaxFunEvals',1000);
+
+%% data
+
+V = 10; % L
+d_p_0 = 50e-8; % dm
+w_s = 0.05;
+c_M_0 = 2; % mol/L
+n_M_0 = c_M_0*V;
+w_I_0 = 0.01;
+a_E = 50*(10^(-9))^2; % dm^2
+N_A = 6.022*10^23; % 1/mol
+
+M_M = 0.1; % kg/mol
+M_I = 0.164; % kg/mol
+p_M = 0.94; % kg/L
+p_P = 1.1; % kg/L
+p_W = 1.0; % kg/L
+c_M = p_M/M_M;
+phi_M = 0.5;
+w_I = 0.01;
+f = 0.5;
+
+k_p = 715; % L/(mol*s)
+k_des = 1e-1; % 1/s
+k_e = 1e-6; % dm/s
+k_t = 9.8*10^6; % L/(mol*s)
+k_d = 5.55e-6; % 1/s
+
+cmc = 40e-3; % mol/L
+K_E = 100; % L/mol
+
+m_M_0 = n_M_0*M_M; % initial mass of monomer (kg)
+V_M_0 = m_M_0/p_M; % initial volume of monomer (L)
+
+%% a)
+
+mass_water = @(m_W_0)(V-p_M*m_M_0-p_W*m_W_0-p_P*(w_s*(m_M_0+m_W_0)/(1-w_s)));
+
+m_W_0 = fsolve(mass_water,0.5,options); % mass of water (kg)
+V_W_0 = m_W_0/p_W; % initial volume of water (L)
+m_tot = m_M_0+m_W_0; % total mass of solution (kg) without seeds
+m_s = w_s*(m_tot)/(1-w_s); % mass of seeds (kg)
+
+V_p_0 = (4/3)*pi*(d_p_0/2)^3; % volume of a seed
+m_1_s = V_p_0*p_P; % mass of one seed without monomer
+
+n_s = m_s/m_1_s; % number of seeds
+V_s_tot = n_s*2*V_p_0; % total volume of swollen seed
+d_p_m = (12*V_p_0/pi)^(1/3); % diameter of swollen seed
+A_s_tot = n_s*4*pi*(d_p_m/2)^2; % total surface of all particles at t = 0
+A_1_s = 4*pi*(d_p_m/2)^2; % surface of one particle at t = 0
+V_1_s = 2*V_p_0; % volume of one particle at t = 0
+
+F_E_inf = (A_1_s/V_1_s)/(a_E*6.02*10^23); % maximum amount of emulsifier adsorbed in mol/vol
+
+F_E_min = 0.2*F_E_inf; % min. adsorbed amount per particle
+E_w_min = (F_E_min/F_E_inf)/(K_E*(1-F_E_min/F_E_inf));
+
+E_w_max = cmc;
+F_E_max = F_E_inf*(K_E*E_w_max)/(1+K_E*E_w_max); % max. adsorbed amount per particle
+
+n_0_min = E_w_min*V_W_0+F_E_min*n_s*V_1_s; % minimum amount in mol
+n_0_max = E_w_max*V_W_0+F_E_max*n_s*V_1_s; % maximum amount in mol
+
+m_I = w_I*m_M_0;
+c_I = m_I/(M_I*V);
+c_R = sqrt(2*f*k_d*c_I/k_t);
+
+%% b)
+
+m0 = [n_M_0 n_s*V_p_0*c_M];
+tspan = 0:1:3600;
+options2 = odeset('RelTol',1e-6,'AbsTol',1e-10);
+[t,m] = ode15s(@(t,n)emulsion1( t,n,n_s,n_M_0,V_p_0 ),tspan,m0,options2);
+
+v_pol = V_p_0+(n_M_0-m(:,1))./n_s*(M_M/p_P);
+v_p = v_pol+(m(:,2)./n_s)*(M_M/p_M);
+
+phi = (v_p-v_pol)./v_p; % phi
+diam = (6.*v_p./pi).^(1/3); % diameter over time (nm)
+V_pol = n_s*v_pol; % total volume of polymer
+V_p = n_s*v_p;
+conv1 = (n_M_0-m(:,1))/n_M_0;
+A_p = 4.*pi.*(diam./2).^2;
+rho = A_p.*k_e.*c_R.*N_A; % 1/s
+n_bar = (rho./(2.*rho));
+
+c_M_1 = m(:,2)./(n_s.*v_p);
+c_R_1 = n_bar./(N_A.*v_p);
+tau_p1 = 1./(k_p.*c_M_1);
+tau_in1 = 1./rho;
+n_N = (tau_in1./tau_p1);
+n_W = 2.*n_N;
+
+figure(1);
+subplot(6,1,1);
+plot(t./3600,conv1);
+title('conversion vs. time');
+xlabel('time / [h]');
+ylabel('conversion');
+
+subplot(6,1,2);
+plot(conv1,V_pol,conv1,V_p);
+axis([0 1 0 2.5]);
+title('total volume of polymer and particles vs. conversion');
+xlabel('x_M');
+ylabel('V / [L]');
+legend('V_{pol}','V_p','Location','best');
+
+subplot(6,1,3);
+plot(conv1,phi);
+axis([0 1 0 0.55]);
+title('\phi vs. conversion');
+xlabel('x_M');
+ylabel('\phi');
+
+subplot(6,1,4);
+plot(conv1,diam.*10^8);
+axis([0 1 50 100]);
+title('diameter of a particle vs. conversion');
+xlabel('x_M');
+ylabel('d / [nm]');
+
+subplot(6,1,5);
+plot(t/3600,n_bar);
+axis([0 1 0 1]);
+title('n vs. conversion');
+xlabel('x_M');
+ylabel('n');
+
+subplot(6,1,6);
+plot(conv1,n_N,conv1,n_W);
+title('number and weight average vs. conversion');
+xlabel('x_M');
+ylabel('n_N, n_W');
+legend('n_N','n_W','Location','best');
+
+%% c)
+
+m0 = [n_M_0 n_s*V_p_0*c_M];
+tspan2 = 0:1:3600;
+[t2,m2] = ode15s(@(t,n)emulsion2( t,n,n_s,n_M_0,V_p_0,V ),tspan2,m0,options2);
+
+v_pol2 = V_p_0+(n_M_0-m2(:,1))./n_s*(M_M/p_P);
+v_p2 = v_pol2+(m2(:,2)./n_s)*(M_M/p_M);
+
+phi = (v_p2-v_pol2)./v_p2; % phi
+diam2 = (6*v_p2./pi).^(1/3); % diameter over time (dm)
+V_pol = n_s*v_pol2; % total volume of polymer
+V_p = n_s*v_p2;
+conv2 = (n_M_0-m2(:,1))/n_M_0;
+A_p2 = 4.*pi.*(diam2./2).^2;
+rho = A_p2.*k_e.*c_R.*N_A; % 1/s
+n_bar = (rho./(2.*rho+k_des));
+
+c_M_2 = m2(:,2)./(n_s.*v_p2);
+c_R_2 = n_bar./(N_A.*v_p2);
+tau_p = 1./(k_p.*c_M_2);
+tau_in = 1./rho;
+tau_out = 1./(k_des.*n_bar);
+
+n_N = 1./(tau_p./tau_in+tau_p./tau_out);
+n_W = 2.*n_N;
+
+figure(2);
+subplot(6,1,1);
+plot(t2./3600,conv2);
+title('conversion vs. time');
+xlabel('time / [h]');
+ylabel('conversion');
+
+subplot(6,1,2);
+plot(conv2,V_pol,conv2,V_p);
+axis([0 1 0 2.5]);
+title('total volume of polymer and particles vs. conversion');
+xlabel('x_M');
+ylabel('V / [L]');
+legend('V_{pol}','V_p','Location','best');
+
+subplot(6,1,3);
+plot(conv2,phi);
+axis([0 1 0 0.55]);
+title('\phi vs. conversion');
+xlabel('x_M');
+ylabel('\phi');
+
+subplot(6,1,4);
+plot(conv2,diam2.*10^8);
+axis([0 1 50 100]);
+title('diameter of a particle vs. conversion');
+xlabel('x_M');
+ylabel('d / [nm]');
+
+subplot(6,1,5);
+plot(t2/3600,n_bar);
+axis([0 1 0 0.5]);
+title('n vs. conversion');
+xlabel('time / [h]');
+ylabel('n');
+
+subplot(6,1,6);
+plot(conv2,n_N,conv2,n_W);
+title('number and weight average vs. conversion');
+xlabel('x_M');
+ylabel('n_N, n_W');
+legend('n_N','n_W','Location','best');
+
+%%