Chapter 3: Deagglomeration and Dispersion Techniques

Learning Objectives

Compare and select appropriate mechanical dispersion methods
Understand cavitation and its role in ultrasonication
Design surfactant-based stabilization strategies
Optimize pH and ionic strength for electrostatic stabilization
Apply machine learning to dispersion process optimization
Select the right method combination for specific applications

3.1 Mechanical Dispersion Methods

Mechanical methods apply physical energy to break inter-particle bonds. The choice depends on agglomerate strength, particle fragility, and scale requirements.

graph TD A[Mechanical Dispersion] --> B[Ultrasonication] A --> C[Ball/Bead Milling] A --> D[High-Pressure Homogenization] A --> E[High-Shear Mixing] B --> B1[Probe sonicator] B --> B2[Bath sonicator] C --> C1[Planetary ball mill] C --> C2[Bead mill] C --> C3[Attritor] D --> D1[Microfluidizer] D --> D2[Valve homogenizer] style A fill:#f9f,stroke:#333 style B fill:#bbf style C fill:#bfb style D fill:#fbb

3.1.1 Ultrasonication

Ultrasonication uses high-frequency sound waves (20-100 kHz) to create cavitation bubbles in liquid. The collapse of these bubbles generates intense local forces that break agglomerates.

Cavitation Mechanism

Acoustic pressure creates microscopic bubbles that grow and violently collapse. At collapse, local conditions can reach:

Temperature: ~5000 K
Pressure: ~1000 atm
Jet velocity: ~100 m/s

These extreme conditions generate shear forces sufficient to break van der Waals bonds.

Example 1: Ultrasonication Parameter Optimization

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# ===================================
# Example 1: Ultrasonication Process Model
# ===================================

class UltrasonicDispersion:
    """Model for ultrasonic dispersion of nanoparticle agglomerates."""

    def __init__(self, d_agglomerate_nm=500, d_primary_nm=20,
                 viscosity_mPas=1.0, density_kg_m3=2200):
        self.d_agg = d_agglomerate_nm
        self.d_pri = d_primary_nm
        self.viscosity = viscosity_mPas * 1e-3
        self.density = density_kg_m3

    def deagglomeration_rate(self, power_W_mL, frequency_kHz, time_s):
        """
        Calculate particle size after ultrasonication.

        Parameters:
            power_W_mL: Power density (W/mL)
            frequency_kHz: Frequency (kHz)
            time_s: Treatment time (s)

        Returns:
            Final particle size (nm)
        """
        # Cavitation intensity (arbitrary units)
        # Higher power and lower frequency → stronger cavitation
        I_cav = power_W_mL * (40 / frequency_kHz)**0.5

        # Rate constant (empirical model)
        k = 0.1 * I_cav / (self.viscosity * 1000)

        # First-order decay toward primary particle size
        d_final = self.d_pri + (self.d_agg - self.d_pri) * np.exp(-k * time_s)

        return max(d_final, self.d_pri)

    def temperature_rise(self, power_W_mL, time_s, volume_mL=10,
                         heat_capacity_J_gK=4.2):
        """Calculate temperature rise during sonication."""
        # Approximate: 50-80% of power converted to heat
        efficiency = 0.6
        mass_g = volume_mL * self.density / 1000
        dT = efficiency * power_W_mL * volume_mL * time_s / (mass_g * heat_capacity_J_gK)
        return dT

    def energy_consumption(self, power_W_mL, time_s, volume_mL=10):
        """Total energy consumption (kJ)."""
        return power_W_mL * volume_mL * time_s / 1000

# Parameter study
sonic = UltrasonicDispersion(d_agglomerate_nm=500, d_primary_nm=20)

# Time evolution at different power levels
times = np.linspace(0, 600, 100)  # 0-10 min
powers = [0.5, 1.0, 2.0, 5.0]  # W/mL

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Plot 1: Size vs time at different powers
ax1 = axes[0, 0]
for P, color in zip(powers, plt.cm.viridis(np.linspace(0.2, 0.9, len(powers)))):
    sizes = [sonic.deagglomeration_rate(P, 20, t) for t in times]
    ax1.plot(times/60, sizes, linewidth=2, label=f'{P} W/mL', color=color)

ax1.axhline(y=20, color='red', linestyle='--', alpha=0.7)
ax1.text(8, 30, 'Primary particle size', fontsize=9, color='red')
ax1.set_xlabel('Time (minutes)', fontsize=11)
ax1.set_ylabel('Particle Size (nm)', fontsize=11)
ax1.set_title('Deagglomeration Kinetics (20 kHz)', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 10)

# Plot 2: Effect of frequency
ax2 = axes[0, 1]
frequencies = [20, 40, 100]  # kHz
P_fixed = 2.0

for f, color in zip(frequencies, ['blue', 'green', 'red']):
    sizes = [sonic.deagglomeration_rate(P_fixed, f, t) for t in times]
    ax2.plot(times/60, sizes, linewidth=2, label=f'{f} kHz', color=color)

ax2.axhline(y=20, color='gray', linestyle='--', alpha=0.7)
ax2.set_xlabel('Time (minutes)', fontsize=11)
ax2.set_ylabel('Particle Size (nm)', fontsize=11)
ax2.set_title(f'Frequency Effect at {P_fixed} W/mL', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0, 10)

# Plot 3: Temperature rise
ax3 = axes[1, 0]
for P, color in zip(powers, plt.cm.viridis(np.linspace(0.2, 0.9, len(powers)))):
    temps = [25 + sonic.temperature_rise(P, t) for t in times]
    ax3.plot(times/60, temps, linewidth=2, label=f'{P} W/mL', color=color)

ax3.axhline(y=60, color='orange', linestyle='--', alpha=0.7)
ax3.axhline(y=80, color='red', linestyle='--', alpha=0.7)
ax3.text(8, 62, 'Caution', fontsize=9, color='orange')
ax3.text(8, 82, 'Risk of degradation', fontsize=9, color='red')
ax3.set_xlabel('Time (minutes)', fontsize=11)
ax3.set_ylabel('Temperature (°C)', fontsize=11)
ax3.set_title('Temperature Rise During Sonication', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)
ax3.set_xlim(0, 10)
ax3.set_ylim(20, 100)

# Plot 4: Process optimization map
ax4 = axes[1, 1]
P_range = np.linspace(0.5, 5, 30)
t_range = np.linspace(60, 600, 30)
P_mesh, t_mesh = np.meshgrid(P_range, t_range)

# Calculate final size for each combination
final_size = np.zeros_like(P_mesh)
for i in range(len(t_range)):
    for j in range(len(P_range)):
        final_size[i, j] = sonic.deagglomeration_rate(P_mesh[i, j], 20, t_mesh[i, j])

contour = ax4.contourf(P_mesh, t_mesh/60, final_size, levels=20, cmap='RdYlGn_r')
ax4.contour(P_mesh, t_mesh/60, final_size, levels=[25, 30, 50], colors='white', linewidths=1.5)
plt.colorbar(contour, ax=ax4, label='Final Size (nm)')
ax4.set_xlabel('Power Density (W/mL)', fontsize=11)
ax4.set_ylabel('Time (minutes)', fontsize=11)
ax4.set_title('Process Window for Target Size', fontsize=12, fontweight='bold')

plt.tight_layout()
plt.savefig('ultrasonication.png', dpi=150, bbox_inches='tight')
plt.show()

# Print recommendations
print("\n=== Ultrasonication Process Recommendations ===")
for P in powers:
    t90 = None
    for t in range(10, 601, 10):
        if sonic.deagglomeration_rate(P, 20, t) < 30:
            t90 = t
            break
    dT = sonic.temperature_rise(P, t90 if t90 else 600)
    if t90:
        print(f"  {P} W/mL: Reach 30 nm in {t90/60:.1f} min (ΔT ≈ {dT:.0f}°C)")
    else:
        print(f"  {P} W/mL: Cannot reach 30 nm in 10 min")

Example Output:
=== Ultrasonication Process Recommendations ===
0.5 W/mL: Cannot reach 30 nm in 10 min
1.0 W/mL: Reach 30 nm in 6.7 min (ΔT ≈ 58°C)
2.0 W/mL: Reach 30 nm in 3.3 min (ΔT ≈ 58°C)
5.0 W/mL: Reach 30 nm in 1.3 min (ΔT ≈ 58°C)

Practical Considerations

Ice bath cooling is essential for temperature-sensitive samples
Pulsed mode (e.g., 5s on / 5s off) reduces heat buildup
Probe tip erosion can contaminate samples with metal particles
Over-sonication can cause re-agglomeration or particle damage

3.1.2 Ball Milling and Bead Milling

Method	Media Size	Final Size Range	Throughput	Contamination Risk
Planetary Ball Mill	3-10 mm	100 nm - 10 μm	Batch (g-kg)	High
Bead Mill	0.1-2 mm	20-500 nm	Continuous	Moderate
High-Energy Ball Mill	5-20 mm	10-100 nm	Batch (g)	Very High

Example 2: Bead Milling Process Simulation

import numpy as np
import matplotlib.pyplot as plt

# ===================================
# Example 2: Bead Milling Process Model
# ===================================

def bead_mill_size_reduction(d_initial_nm, bead_size_mm, speed_rpm,
                              solid_loading, time_min, media_fill=0.8):
    """
    Model particle size reduction in bead mill.

    Parameters:
        d_initial_nm: Initial particle/agglomerate size
        bead_size_mm: Grinding media diameter
        speed_rpm: Mill speed
        solid_loading: Solid content (weight fraction)
        time_min: Milling time
        media_fill: Volume fraction of beads

    Returns:
        Final particle size (nm)
    """
    # Stress intensity factor
    stress_intensity = (bead_size_mm * 1e-3)**3 * (speed_rpm / 1000)**2

    # Number of stress events depends on bead count
    n_beads = media_fill / (bead_size_mm * 1e-3)**3

    # Effective rate constant
    k = 0.01 * n_beads * stress_intensity / (solid_loading + 0.1)

    # Size reduction (limit at ~20 nm for typical oxide particles)
    d_min = 20  # nm, practical limit
    d_final = d_min + (d_initial_nm - d_min) * np.exp(-k * time_min)

    return max(d_final, d_min)

# Compare bead sizes
d_initial = 500  # nm
times = np.linspace(0, 120, 100)  # 0-2 hours

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Plot 1: Effect of bead size
ax1 = axes[0]
bead_sizes = [0.1, 0.3, 0.5, 1.0, 2.0]
colors = plt.cm.plasma(np.linspace(0.1, 0.9, len(bead_sizes)))

for bs, color in zip(bead_sizes, colors):
    sizes = [bead_mill_size_reduction(d_initial, bs, 2000, 0.2, t) for t in times]
    ax1.plot(times, sizes, linewidth=2, label=f'{bs} mm beads', color=color)

ax1.axhline(y=50, color='green', linestyle='--', alpha=0.7)
ax1.text(100, 55, 'Target: 50 nm', fontsize=9, color='green')
ax1.set_xlabel('Milling Time (min)', fontsize=11)
ax1.set_ylabel('Particle Size (nm)', fontsize=11)
ax1.set_title('Effect of Bead Size on Milling', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 120)
ax1.set_ylim(0, 550)

# Plot 2: Energy efficiency comparison
ax2 = axes[1]
# Specific energy consumption (kWh/kg) for different methods
methods = ['Ultrasonication', 'Bead Mill\n(0.3mm)', 'Bead Mill\n(1.0mm)',
           'High-Pressure\nHomogenizer', 'Planetary\nBall Mill']
energies = [5.0, 1.5, 3.0, 0.8, 8.0]
final_sizes = [25, 30, 60, 40, 80]

colors = plt.cm.RdYlGn_r(np.array(final_sizes) / 100)
bars = ax2.bar(methods, energies, color=colors, edgecolor='black')

# Add final size annotations
for bar, size in zip(bars, final_sizes):
    ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.2,
             f'{size} nm', ha='center', fontsize=9)

ax2.set_ylabel('Specific Energy (kWh/kg)', fontsize=11)
ax2.set_title('Energy Efficiency Comparison\n(reaching respective minimum size)', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig('bead_milling.png', dpi=150, bbox_inches='tight')
plt.show()

# Print selection guide
print("\n=== Bead Size Selection Guide ===")
print("Target Size    Recommended Bead Size    Notes")
print("-" * 55)
print("20-50 nm       0.1-0.3 mm               Yttria-stabilized ZrO2")
print("50-100 nm      0.3-0.5 mm               Standard ceramic beads")
print("100-500 nm     0.5-1.0 mm               Glass or steel beads")
print(">500 nm        1.0-2.0 mm               Cost-effective grinding")

3.2 Chemical Dispersion Methods

3.2.1 Surfactant Selection

Surfactants adsorb at particle surfaces, providing steric and/or electrostatic stabilization.

Surfactant Type	Examples	Stabilization	Best For
Anionic	SDS, sodium citrate	Electrostatic	Aqueous, low ionic strength
Cationic	CTAB, DTAB	Electrostatic	Negatively charged surfaces
Non-ionic	Tween, Triton, PEG	Steric	High ionic strength, biomedical
Polymeric	PVP, PVA, PAA	Steric + electrostatic	Long-term stability

Example 3: Surfactant Concentration Optimization

import numpy as np
import matplotlib.pyplot as plt

# ===================================
# Example 3: Surfactant Adsorption and Stability
# ===================================

def langmuir_adsorption(C_surf, K_ads, Gamma_max):
    """
    Langmuir adsorption isotherm.

    Parameters:
        C_surf: Surfactant concentration (mM)
        K_ads: Adsorption constant (1/mM)
        Gamma_max: Maximum surface coverage (molecules/nm²)

    Returns:
        Surface coverage (molecules/nm²)
    """
    return Gamma_max * K_ads * C_surf / (1 + K_ads * C_surf)

def surface_charge_density(Gamma, charge_per_molecule):
    """Calculate surface charge from adsorbed surfactant."""
    return Gamma * charge_per_molecule * 1.602e-19 * 1e18  # C/m²

def stability_from_coverage(Gamma, Gamma_max, type='electrostatic'):
    """
    Estimate stability based on surface coverage.

    Returns stability score 0-100.
    """
    coverage_fraction = Gamma / Gamma_max

    if type == 'electrostatic':
        # Need sufficient charge density
        if coverage_fraction < 0.3:
            return coverage_fraction / 0.3 * 50
        elif coverage_fraction < 0.7:
            return 50 + (coverage_fraction - 0.3) / 0.4 * 45
        else:
            return 95 + (coverage_fraction - 0.7) / 0.3 * 5
    else:  # steric
        # Need dense polymer brush
        if coverage_fraction < 0.5:
            return coverage_fraction / 0.5 * 60
        else:
            return 60 + (coverage_fraction - 0.5) / 0.5 * 40

# Surfactant parameters
surfactants = {
    'SDS (anionic)': {'K': 5.0, 'Gamma_max': 4.0, 'charge': -1, 'type': 'electrostatic'},
    'CTAB (cationic)': {'K': 8.0, 'Gamma_max': 3.5, 'charge': 1, 'type': 'electrostatic'},
    'Tween-80 (non-ionic)': {'K': 2.0, 'Gamma_max': 1.5, 'charge': 0, 'type': 'steric'},
    'PVP (polymer)': {'K': 0.5, 'Gamma_max': 0.8, 'charge': 0, 'type': 'steric'}
}

C_range = np.logspace(-2, 2, 100)  # 0.01 - 100 mM

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Plot 1: Adsorption isotherms
ax1 = axes[0, 0]
for name, params in surfactants.items():
    Gamma = [langmuir_adsorption(C, params['K'], params['Gamma_max']) for C in C_range]
    ax1.semilogx(C_range, Gamma, linewidth=2, label=name)

ax1.axhline(y=2.0, color='gray', linestyle='--', alpha=0.7)
ax1.text(50, 2.1, 'Monolayer ~2 mol/nm²', fontsize=9, color='gray')
ax1.set_xlabel('Surfactant Concentration (mM)', fontsize=11)
ax1.set_ylabel('Surface Coverage (molecules/nm²)', fontsize=11)
ax1.set_title('Adsorption Isotherms', fontsize=12, fontweight='bold')
ax1.legend(loc='lower right', fontsize=9)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0.01, 100)

# Plot 2: Stability vs concentration
ax2 = axes[0, 1]
for name, params in surfactants.items():
    stabilities = []
    for C in C_range:
        Gamma = langmuir_adsorption(C, params['K'], params['Gamma_max'])
        stab = stability_from_coverage(Gamma, params['Gamma_max'], params['type'])
        stabilities.append(stab)
    ax2.semilogx(C_range, stabilities, linewidth=2, label=name)

ax2.axhline(y=80, color='green', linestyle='--', alpha=0.7)
ax2.text(0.015, 82, 'Stable (>80)', fontsize=9, color='green')
ax2.axhline(y=50, color='orange', linestyle='--', alpha=0.7)
ax2.text(0.015, 52, 'Marginal (50-80)', fontsize=9, color='orange')
ax2.set_xlabel('Surfactant Concentration (mM)', fontsize=11)
ax2.set_ylabel('Stability Score', fontsize=11)
ax2.set_title('Dispersion Stability vs Concentration', fontsize=12, fontweight='bold')
ax2.legend(loc='lower right', fontsize=9)
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0.01, 100)
ax2.set_ylim(0, 100)

# Plot 3: Optimal concentration determination
ax3 = axes[1, 0]
# For SDS - detailed analysis
params = surfactants['SDS (anionic)']
C_detail = np.linspace(0.1, 20, 100)

Gamma_values = [langmuir_adsorption(C, params['K'], params['Gamma_max']) for C in C_detail]
stability_values = [stability_from_coverage(G, params['Gamma_max'], 'electrostatic')
                   for G in Gamma_values]
# Cost increases linearly
cost_values = C_detail / 20 * 100

# Efficiency = stability / cost
efficiency = np.array(stability_values) / (np.array(cost_values) + 10)
efficiency = efficiency / efficiency.max() * 100

ax3.plot(C_detail, stability_values, 'b-', linewidth=2, label='Stability')
ax3.plot(C_detail, cost_values, 'r--', linewidth=2, label='Relative Cost')
ax3.plot(C_detail, efficiency, 'g-', linewidth=2, label='Efficiency')

# Find optimal
opt_idx = np.argmax(efficiency)
ax3.axvline(x=C_detail[opt_idx], color='purple', linestyle=':', linewidth=2)
ax3.scatter([C_detail[opt_idx]], [efficiency[opt_idx]], s=100, c='purple', zorder=5)
ax3.text(C_detail[opt_idx] + 0.5, efficiency[opt_idx] + 5,
         f'Optimal: {C_detail[opt_idx]:.1f} mM', fontsize=10, color='purple')

ax3.set_xlabel('SDS Concentration (mM)', fontsize=11)
ax3.set_ylabel('Score / Cost', fontsize=11)
ax3.set_title('SDS: Finding Optimal Concentration', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)
ax3.set_xlim(0, 20)
ax3.set_ylim(0, 110)

# Plot 4: Surfactant selection guide
ax4 = axes[1, 1]
conditions = ['Low I\nAqueous', 'High I\nAqueous', 'Biomedical', 'Organic\nSolvent', 'Food\nGrade']
surfactant_scores = {
    'SDS': [95, 30, 20, 10, 40],
    'CTAB': [90, 25, 10, 15, 10],
    'Tween-80': [70, 85, 95, 60, 90],
    'PVP': [80, 90, 85, 70, 75]
}

x = np.arange(len(conditions))
width = 0.2
multiplier = 0

for name, scores in surfactant_scores.items():
    offset = width * multiplier
    ax4.bar(x + offset, scores, width, label=name, alpha=0.8)
    multiplier += 1

ax4.set_xticks(x + width * 1.5)
ax4.set_xticklabels(conditions)
ax4.set_ylabel('Suitability Score', fontsize=11)
ax4.set_title('Surfactant Selection by Application', fontsize=12, fontweight='bold')
ax4.legend(loc='upper right')
ax4.grid(True, alpha=0.3, axis='y')
ax4.set_ylim(0, 100)

plt.tight_layout()
plt.savefig('surfactant_selection.png', dpi=150, bbox_inches='tight')
plt.show()

3.2.2 Surface Functionalization

Covalent surface modification provides permanent stabilization that survives processing and storage.

Common Functionalization Strategies

Silane coupling: For oxides (SiO₂, TiO₂, Fe₃O₄) - APTES, MPTMS, OTS
Thiol SAMs: For noble metals (Au, Ag) - thiols, disulfides
Carboxylate coordination: For metal oxides - citrate, oleate
Polymer grafting: "Grafting-from" or "grafting-to" approaches

3.3 Physicochemical Methods

3.3.1 pH Control

pH determines surface charge through protonation/deprotonation of surface groups.

Example 4: pH Optimization for Electrostatic Stabilization

import numpy as np
import matplotlib.pyplot as plt

# ===================================
# Example 4: pH-Dependent Zeta Potential and Stability
# ===================================

def zeta_potential_vs_pH(pH, IEP, zeta_max=50):
    """
    Model zeta potential as function of pH.

    Parameters:
        pH: Solution pH
        IEP: Isoelectric point
        zeta_max: Maximum |ζ| (mV)

    Returns:
        Zeta potential (mV)
    """
    # Sigmoidal transition around IEP
    return zeta_max * np.tanh((IEP - pH) * 0.8)

def stability_criterion(zeta_mV):
    """
    Stability based on zeta potential.
    |ζ| > 30 mV generally indicates good stability.
    """
    abs_zeta = np.abs(zeta_mV)
    if abs_zeta > 40:
        return 'Excellent'
    elif abs_zeta > 30:
        return 'Good'
    elif abs_zeta > 20:
        return 'Moderate'
    else:
        return 'Poor'

# Common oxide materials and their IEPs
materials = {
    'SiO₂': {'IEP': 2.0, 'color': 'blue'},
    'TiO₂ (anatase)': {'IEP': 6.0, 'color': 'green'},
    'Al₂O₃': {'IEP': 9.0, 'color': 'red'},
    'Fe₃O₄': {'IEP': 6.5, 'color': 'orange'},
    'ZnO': {'IEP': 9.5, 'color': 'purple'}
}

pH_range = np.linspace(1, 13, 200)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Plot 1: Zeta potential vs pH
ax1 = axes[0]
for name, props in materials.items():
    zeta = [zeta_potential_vs_pH(pH, props['IEP']) for pH in pH_range]
    ax1.plot(pH_range, zeta, linewidth=2, label=name, color=props['color'])

ax1.axhline(y=30, color='gray', linestyle='--', alpha=0.7)
ax1.axhline(y=-30, color='gray', linestyle='--', alpha=0.7)
ax1.axhline(y=0, color='black', linestyle='-', linewidth=0.5)
ax1.fill_between(pH_range, -30, 30, alpha=0.1, color='red', label='Unstable region')
ax1.set_xlabel('pH', fontsize=11)
ax1.set_ylabel('Zeta Potential (mV)', fontsize=11)
ax1.set_title('Zeta Potential vs pH for Common Oxides', fontsize=12, fontweight='bold')
ax1.legend(loc='upper right', fontsize=9)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(1, 13)
ax1.set_ylim(-60, 60)

# Plot 2: pH stability windows
ax2 = axes[1]
y_pos = np.arange(len(materials))
bar_height = 0.6

for i, (name, props) in enumerate(materials.items()):
    IEP = props['IEP']

    # Stable at low pH (below IEP - 2)
    if IEP > 3:
        ax2.barh(i, IEP - 3, left=1, height=bar_height, color='green', alpha=0.7)

    # Unstable around IEP (IEP ± 2)
    ax2.barh(i, 4, left=max(1, IEP - 2), height=bar_height, color='red', alpha=0.3)

    # Stable at high pH (above IEP + 2)
    if IEP < 11:
        ax2.barh(i, 13 - (IEP + 2), left=IEP + 2, height=bar_height, color='green', alpha=0.7)

    # Mark IEP
    ax2.scatter([IEP], [i], marker='|', s=200, color='black', zorder=5)
    ax2.text(IEP, i + 0.35, f'IEP={IEP}', ha='center', fontsize=9)

ax2.set_yticks(y_pos)
ax2.set_yticklabels(materials.keys())
ax2.set_xlabel('pH', fontsize=11)
ax2.set_title('pH Stability Windows', fontsize=12, fontweight='bold')
ax2.set_xlim(1, 13)
ax2.grid(True, alpha=0.3, axis='x')

# Legend
from matplotlib.patches import Patch
legend_elements = [Patch(facecolor='green', alpha=0.7, label='Stable (|ζ| > 30 mV)'),
                   Patch(facecolor='red', alpha=0.3, label='Unstable (near IEP)')]
ax2.legend(handles=legend_elements, loc='lower right')

plt.tight_layout()
plt.savefig('pH_optimization.png', dpi=150, bbox_inches='tight')
plt.show()

# Print pH recommendations
print("\n=== pH Recommendations for Stable Dispersions ===")
print(f"{'Material':20s} {'IEP':>6s} {'Stable pH ranges':>25s}")
print("-" * 55)
for name, props in materials.items():
    IEP = props['IEP']
    ranges = []
    if IEP > 4:
        ranges.append(f"pH < {IEP-2:.0f}")
    if IEP < 10:
        ranges.append(f"pH > {IEP+2:.0f}")
    print(f"{name:20s} {IEP:>6.1f} {' or '.join(ranges):>25s}")

Example Output:
=== pH Recommendations for Stable Dispersions ===
Material IEP Stable pH ranges
-------------------------------------------------------
SiO₂ 2.0 pH > 4
TiO₂ (anatase) 6.0 pH < 4 or pH > 8
Al₂O₃ 9.0 pH < 7 or pH > 11
Fe₃O₄ 6.5 pH < 5 or pH > 9
ZnO 9.5 pH < 8

3.4 Process Optimization with Machine Learning

Example 5: Bayesian Optimization of Dispersion Process

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.optimize import minimize

# ===================================
# Example 5: Bayesian Optimization for Dispersion Process
# ===================================

class DispersionProcessOptimizer:
    """
    Optimize dispersion process using Bayesian optimization.

    Process parameters to optimize:
    - Ultrasonication power (W/mL)
    - Surfactant concentration (mM)
    - pH

    Objective: Minimize particle size while maintaining stability.
    """

    def __init__(self):
        self.bounds = {
            'power': (0.5, 5.0),      # W/mL
            'surfactant': (0.1, 10),   # mM
            'pH': (3, 11)
        }

        self.X_observed = []
        self.Y_observed = []

    def true_objective(self, power, surfactant, pH):
        """
        True objective function (unknown in practice).
        Lower is better (minimizing particle size × instability factor).
        """
        # Particle size depends mainly on power
        d_particle = 20 + 400 * np.exp(-0.5 * power)

        # Stability depends on surfactant and pH
        # Optimal surfactant around 2-5 mM
        surf_factor = 1 + 5 * np.exp(-((surfactant - 3) / 2)**2)

        # pH stability (optimal around 4 or 10, poor around 7)
        pH_factor = 1 + 3 * np.exp(-((pH - 7) / 2)**2)

        # Combined objective
        objective = d_particle * surf_factor * pH_factor / 100

        # Add noise
        return objective + np.random.normal(0, 0.5)

    def evaluate(self, params):
        """Evaluate at given parameters."""
        power, surfactant, pH = params
        return self.true_objective(power, surfactant, pH)

    def gp_predict(self, X_test):
        """Simple Gaussian process prediction."""
        if len(self.X_observed) < 2:
            return 50, 10  # Prior

        X_obs = np.array(self.X_observed)
        Y_obs = np.array(self.Y_observed)

        # Normalize
        X_mean, X_std = X_obs.mean(axis=0), X_obs.std(axis=0) + 1e-6
        X_norm = (X_obs - X_mean) / X_std
        X_test_norm = (X_test - X_mean) / X_std

        # Distance-based prediction
        dists = np.linalg.norm(X_norm - X_test_norm, axis=1)
        weights = np.exp(-dists**2 / 0.5)
        weights = weights / (weights.sum() + 1e-10)

        mean = weights @ Y_obs
        std = 5 * np.exp(-np.min(dists))

        return mean, std

    def expected_improvement(self, X_test):
        """Calculate expected improvement."""
        mean, std = self.gp_predict(X_test)
        best_y = min(self.Y_observed) if self.Y_observed else 100

        z = (best_y - mean) / (std + 1e-9)
        ei = (best_y - mean) * norm.cdf(z) + std * norm.pdf(z)
        return ei

    def optimize(self, n_iterations=20, n_initial=5):
        """Run Bayesian optimization."""
        # Initial random sampling
        np.random.seed(42)
        for _ in range(n_initial):
            params = [
                np.random.uniform(*self.bounds['power']),
                np.random.uniform(*self.bounds['surfactant']),
                np.random.uniform(*self.bounds['pH'])
            ]
            y = self.evaluate(params)
            self.X_observed.append(params)
            self.Y_observed.append(y)

        # Sequential optimization
        for i in range(n_iterations):
            # Find next point by maximizing EI
            best_ei = -1
            best_x = None

            for _ in range(50):  # Random search for EI maximum
                x_candidate = [
                    np.random.uniform(*self.bounds['power']),
                    np.random.uniform(*self.bounds['surfactant']),
                    np.random.uniform(*self.bounds['pH'])
                ]
                ei = self.expected_improvement(x_candidate)
                if ei > best_ei:
                    best_ei = ei
                    best_x = x_candidate

            # Evaluate and update
            y_new = self.evaluate(best_x)
            self.X_observed.append(best_x)
            self.Y_observed.append(y_new)

            if (i + 1) % 5 == 0:
                best_idx = np.argmin(self.Y_observed)
                print(f"Iteration {i+1}: Best objective = {self.Y_observed[best_idx]:.2f}")
                print(f"  Parameters: Power={self.X_observed[best_idx][0]:.2f}, "
                      f"Surf={self.X_observed[best_idx][1]:.2f}, "
                      f"pH={self.X_observed[best_idx][2]:.2f}")

        return self.X_observed, self.Y_observed

# Run optimization
print("=== Bayesian Optimization of Dispersion Process ===\n")
optimizer = DispersionProcessOptimizer()
X_history, Y_history = optimizer.optimize(n_iterations=20)

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Plot 1: Convergence
ax1 = axes[0, 0]
best_so_far = [min(Y_history[:i+1]) for i in range(len(Y_history))]
ax1.plot(range(1, len(Y_history)+1), Y_history, 'o-', alpha=0.5, label='Observed')
ax1.plot(range(1, len(Y_history)+1), best_so_far, 'r-', linewidth=2, label='Best so far')
ax1.axhline(y=min(Y_history), color='green', linestyle='--', alpha=0.7)
ax1.set_xlabel('Iteration', fontsize=11)
ax1.set_ylabel('Objective Value', fontsize=11)
ax1.set_title('Optimization Convergence', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Plot 2: Parameter evolution - Power
ax2 = axes[0, 1]
powers = [x[0] for x in X_history]
surfactants = [x[1] for x in X_history]
pHs = [x[2] for x in X_history]

ax2.scatter(range(1, len(powers)+1), powers, c=Y_history, cmap='RdYlGn_r', s=80)
ax2.set_xlabel('Iteration', fontsize=11)
ax2.set_ylabel('Power (W/mL)', fontsize=11)
ax2.set_title('Power Parameter Exploration', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)

# Plot 3: Parameter space (Power vs Surfactant)
ax3 = axes[1, 0]
scatter = ax3.scatter(powers, surfactants, c=Y_history, cmap='RdYlGn_r', s=100, edgecolor='black')
best_idx = np.argmin(Y_history)
ax3.scatter([powers[best_idx]], [surfactants[best_idx]], s=300, marker='*',
            c='gold', edgecolor='black', linewidth=2, label='Optimum', zorder=5)
plt.colorbar(scatter, ax=ax3, label='Objective')
ax3.set_xlabel('Power (W/mL)', fontsize=11)
ax3.set_ylabel('Surfactant (mM)', fontsize=11)
ax3.set_title('Parameter Space Exploration', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)

# Plot 4: Final results summary
ax4 = axes[1, 1]
best_params = X_history[best_idx]
param_names = ['Power\n(W/mL)', 'Surfactant\n(mM)', 'pH']
param_values = best_params
param_ranges = [(0.5, 5.0), (0.1, 10), (3, 11)]

# Normalize to 0-1 scale
norm_values = [(v - r[0]) / (r[1] - r[0]) for v, r in zip(param_values, param_ranges)]

bars = ax4.barh(param_names, norm_values, color=['blue', 'green', 'red'], alpha=0.7)

# Add actual values
for bar, val in zip(bars, param_values):
    ax4.text(bar.get_width() + 0.02, bar.get_y() + bar.get_height()/2,
             f'{val:.2f}', va='center', fontsize=11)

ax4.set_xlabel('Normalized Parameter Value', fontsize=11)
ax4.set_title(f'Optimal Parameters (Objective = {Y_history[best_idx]:.2f})', fontsize=12, fontweight='bold')
ax4.set_xlim(0, 1.2)
ax4.grid(True, alpha=0.3, axis='x')

plt.tight_layout()
plt.savefig('bayesian_optimization.png', dpi=150, bbox_inches='tight')
plt.show()

# Final summary
print("\n=== Optimization Results ===")
print(f"Best objective value: {Y_history[best_idx]:.2f}")
print(f"Optimal parameters:")
print(f"  - Ultrasonication power: {best_params[0]:.2f} W/mL")
print(f"  - Surfactant concentration: {best_params[1]:.2f} mM")
print(f"  - pH: {best_params[2]:.2f}")

3.5 Method Selection Guide

Criterion	Ultrasonication	Bead Milling	HPH	Chemical
Agglomerate strength	Weak-Medium	Strong	Medium	Any
Scale	mL-L	L-m³	L-m³	Any
Contamination	Low	Medium	Very low	None
Energy efficiency	Low	Medium	High	Very high
Particle damage	Medium	High	Medium	None

Recommended Approach

Start with chemical stabilization - surfactant or surface modification
Apply gentle mechanical dispersion - sonication or high-shear mixing
Optimize pH and ionic strength - maximize electrostatic repulsion
Use aggressive methods only if needed - bead milling for strong aggregates
Validate stability - measure zeta potential and particle size

Chapter Summary

Key Takeaways

Ultrasonication works via cavitation; effective for weak agglomerates but causes heating
Bead milling provides high energy; use smaller beads (0.1-0.3 mm) for finer particles
Surfactant selection depends on medium and application; optimize concentration
pH control is critical - stay away from IEP for electrostatic stability
Combined approaches (mechanical + chemical) usually work best
Machine learning can efficiently optimize multi-parameter processes

Exercises

Exercise 1: Ultrasonication protocol design

You have 500 nm TiO₂ agglomerates and need to reach <50 nm. Design a sonication protocol including power, time, and cooling strategy.

Exercise 2: Surfactant selection

Select an appropriate surfactant for dispersing 30 nm silica nanoparticles in (a) pure water, (b) 100 mM NaCl, (c) cell culture medium. Justify your choices.

Exercise 3: pH optimization

You need to disperse a mixture of TiO₂ (IEP = 6) and Al₂O₃ (IEP = 9) in water. What pH would you recommend? What challenges might arise?

← Chapter 2: Factors Affecting Agglomeration Chapter 4: Stability Evaluation →