Chapter 4: Stability Evaluation Methods

4.1 Zeta Potential Measurement

Zeta potential (ζ) is the electrical potential at the slipping plane of a particle in solution. It is the primary indicator of electrostatic stability.

Example 1: Zeta Potential Analysis and pH Titration

import numpy as np
import matplotlib.pyplot as plt

# ===================================
# Example 1: Zeta Potential Analysis
# ===================================

class ZetaPotentialAnalyzer:
    """Analyze zeta potential data and predict stability."""

    def __init__(self, sample_name, iep=None):
        self.sample_name = sample_name
        self.iep = iep
        self.measurements = []

    def add_measurement(self, pH, zeta_mV, conductivity_mS=None):
        """Add a zeta potential measurement."""
        self.measurements.append({
            'pH': pH,
            'zeta': zeta_mV,
            'conductivity': conductivity_mS
        })

    def fit_sigmoidal(self):
        """Fit sigmoidal curve to pH titration data."""
        if len(self.measurements) < 3:
            return None

        pHs = np.array([m['pH'] for m in self.measurements])
        zetas = np.array([m['zeta'] for m in self.measurements])

        # Estimate IEP from sign change
        for i in range(len(zetas) - 1):
            if zetas[i] * zetas[i+1] < 0:
                # Linear interpolation
                self.iep = pHs[i] - zetas[i] * (pHs[i+1] - pHs[i]) / (zetas[i+1] - zetas[i])
                break

        return self.iep

    def stability_at_pH(self, pH):
        """Predict stability at given pH."""
        if not self.measurements:
            return None

        # Find closest measurements
        pHs = np.array([m['pH'] for m in self.measurements])
        zetas = np.array([m['zeta'] for m in self.measurements])

        # Interpolate
        zeta_interp = np.interp(pH, pHs, zetas)
        abs_zeta = abs(zeta_interp)

        if abs_zeta > 40:
            return {'rating': 'Excellent', 'score': 100, 'zeta': zeta_interp}
        elif abs_zeta > 30:
            return {'rating': 'Good', 'score': 80, 'zeta': zeta_interp}
        elif abs_zeta > 20:
            return {'rating': 'Moderate', 'score': 50, 'zeta': zeta_interp}
        else:
            return {'rating': 'Poor', 'score': 20, 'zeta': zeta_interp}

    def plot_titration(self, ax=None):
        """Plot pH titration curve."""
        if ax is None:
            fig, ax = plt.subplots(figsize=(10, 6))

        pHs = [m['pH'] for m in self.measurements]
        zetas = [m['zeta'] for m in self.measurements]

        ax.plot(pHs, zetas, 'o-', markersize=8, linewidth=2, label=self.sample_name)
        ax.axhline(y=0, color='black', linewidth=0.5)
        ax.axhline(y=30, color='green', linestyle='--', alpha=0.7)
        ax.axhline(y=-30, color='green', linestyle='--', alpha=0.7)

        if self.iep:
            ax.axvline(x=self.iep, color='red', linestyle=':', alpha=0.7)
            ax.text(self.iep + 0.1, max(zetas) * 0.9, f'IEP = {self.iep:.1f}', color='red')

        ax.fill_between(pHs, -30, 30, alpha=0.1, color='red')
        ax.set_xlabel('pH', fontsize=11)
        ax.set_ylabel('Zeta Potential (mV)', fontsize=11)
        ax.set_title(f'Zeta Potential Titration: {self.sample_name}', fontsize=12, fontweight='bold')
        ax.legend()
        ax.grid(True, alpha=0.3)

        return ax

# Simulate titration data for different materials
np.random.seed(42)

silica = ZetaPotentialAnalyzer('SiO₂ (20 nm)')
for pH in np.linspace(2, 10, 15):
    # SiO₂ IEP ≈ 2, strongly negative at higher pH
    zeta = -50 * (1 - np.exp(-(pH - 2) / 2)) + np.random.normal(0, 2)
    silica.add_measurement(pH, zeta)
silica.fit_sigmoidal()

titania = ZetaPotentialAnalyzer('TiO₂ (30 nm)')
for pH in np.linspace(2, 10, 15):
    # TiO₂ IEP ≈ 6
    zeta = 40 * np.tanh((6 - pH) * 0.8) + np.random.normal(0, 2)
    titania.add_measurement(pH, zeta)
titania.fit_sigmoidal()

alumina = ZetaPotentialAnalyzer('Al₂O₃ (50 nm)')
for pH in np.linspace(2, 12, 15):
    # Al₂O₃ IEP ≈ 9
    zeta = 50 * np.tanh((9 - pH) * 0.6) + np.random.normal(0, 2)
    alumina.add_measurement(pH, zeta)
alumina.fit_sigmoidal()

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

# Plot 1: All titration curves
ax1 = axes[0]
for sample in [silica, titania, alumina]:
    sample.plot_titration(ax1)
ax1.set_xlim(2, 12)
ax1.set_ylim(-70, 70)

# Plot 2: Stability windows
ax2 = axes[1]
pHs = np.linspace(2, 12, 100)

for sample, color in [(silica, 'blue'), (titania, 'green'), (alumina, 'red')]:
    scores = [sample.stability_at_pH(pH)['score'] for pH in pHs]
    ax2.plot(pHs, scores, linewidth=2, label=sample.sample_name, color=color)

ax2.axhline(y=80, color='gray', linestyle='--', alpha=0.7)
ax2.text(11, 82, 'Good stability', fontsize=9, color='gray')
ax2.set_xlabel('pH', fontsize=11)
ax2.set_ylabel('Stability Score', fontsize=11)
ax2.set_title('Predicted Stability vs pH', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)
ax2.set_xlim(2, 12)
ax2.set_ylim(0, 110)

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

# Print IEP values and recommended pH ranges
print("\n=== Zeta Potential Analysis Summary ===")
for sample in [silica, titania, alumina]:
    print(f"\n{sample.sample_name}:")
    print(f"  IEP = {sample.iep:.1f}")

    # Find stable pH ranges
    stable_ranges = []
    in_stable = False
    start_pH = None

    for pH in np.linspace(2, 12, 100):
        stab = sample.stability_at_pH(pH)
        if stab['score'] >= 80 and not in_stable:
            in_stable = True
            start_pH = pH
        elif stab['score'] < 80 and in_stable:
            in_stable = False
            stable_ranges.append(f"pH {start_pH:.1f}-{pH:.1f}")

    if in_stable:
        stable_ranges.append(f"pH {start_pH:.1f}-12.0")

    print(f"  Stable pH ranges: {', '.join(stable_ranges)}")

4.2 DLVO Theory

DLVO theory (Derjaguin-Landau-Verwey-Overbeek) describes colloidal stability as the balance between van der Waals attraction and electrostatic repulsion.

Example 2: DLVO Interaction Energy Calculation

import numpy as np
import matplotlib.pyplot as plt

# ===================================
# Example 2: DLVO Theory Implementation
# ===================================

class DLVOCalculator:
    """Calculate DLVO interaction energy profiles."""

    # Physical constants
    eps_0 = 8.854e-12      # Vacuum permittivity (F/m)
    kB = 1.38e-23          # Boltzmann constant (J/K)
    e = 1.602e-19          # Elementary charge (C)
    NA = 6.022e23          # Avogadro's number

    def __init__(self, R_nm, A_H, psi_mV, eps_r=78.5, T=298):
        """
        Initialize DLVO calculator.

        Parameters:
            R_nm: Particle radius (nm)
            A_H: Hamaker constant (J)
            psi_mV: Surface potential (mV)
            eps_r: Relative permittivity
            T: Temperature (K)
        """
        self.R = R_nm * 1e-9
        self.A_H = A_H
        self.psi = psi_mV * 1e-3
        self.eps_r = eps_r
        self.T = T

    def debye_length(self, ionic_strength_M):
        """Calculate Debye length (m)."""
        I = ionic_strength_M * 1000  # mol/m³
        kappa_inv = np.sqrt(
            (self.eps_r * self.eps_0 * self.kB * self.T) /
            (2 * self.NA * self.e**2 * I)
        )
        return kappa_inv

    def v_vdw(self, h):
        """Van der Waals attraction (J)."""
        h = np.maximum(h, 1e-12)  # Avoid division by zero
        return -self.A_H * self.R / (12 * h)

    def v_elec(self, h, ionic_strength_M):
        """Electrostatic repulsion (J)."""
        h = np.maximum(h, 1e-12)
        kappa = 1 / self.debye_length(ionic_strength_M)

        V = (2 * np.pi * self.eps_r * self.eps_0 * self.R *
             self.psi**2 * np.log(1 + np.exp(-kappa * h)))
        return V

    def v_total(self, h, ionic_strength_M):
        """Total DLVO interaction (J)."""
        return self.v_vdw(h) + self.v_elec(h, ionic_strength_M)

    def energy_barrier(self, ionic_strength_M, h_range=None):
        """Find maximum energy barrier."""
        if h_range is None:
            h_range = np.linspace(0.3e-9, 50e-9, 1000)

        V = self.v_total(h_range, ionic_strength_M)
        V_kT = V / (self.kB * self.T)

        max_idx = np.argmax(V_kT)
        return {
            'barrier_kT': V_kT[max_idx],
            'distance_nm': h_range[max_idx] * 1e9,
            'stable': V_kT[max_idx] > 15
        }

    def plot_profile(self, ionic_strengths_M, ax=None):
        """Plot DLVO energy profiles."""
        if ax is None:
            fig, ax = plt.subplots(figsize=(10, 6))

        h_range = np.linspace(0.3e-9, 30e-9, 500)

        for I in ionic_strengths_M:
            V_total = self.v_total(h_range, I)
            V_kT = V_total / (self.kB * self.T)
            lambda_D = self.debye_length(I) * 1e9
            ax.plot(h_range * 1e9, V_kT, linewidth=2,
                    label=f'I = {I*1000:.0f} mM (λD = {lambda_D:.1f} nm)')

        ax.axhline(y=0, color='black', linewidth=0.5)
        ax.axhline(y=15, color='green', linestyle='--', alpha=0.7)
        ax.text(25, 17, 'Stability threshold (15 kT)', fontsize=9, color='green')

        ax.set_xlabel('Separation Distance (nm)', fontsize=11)
        ax.set_ylabel('Interaction Energy (kT)', fontsize=11)
        ax.set_title(f'DLVO Profile: R = {self.R*1e9:.0f} nm, ψ₀ = {self.psi*1000:.0f} mV',
                     fontsize=12, fontweight='bold')
        ax.legend()
        ax.grid(True, alpha=0.3)
        ax.set_xlim(0, 30)

        return ax

# Create DLVO calculator for silica nanoparticles
silica_dlvo = DLVOCalculator(
    R_nm=15,           # 30 nm diameter
    A_H=0.83e-20,      # SiO₂-SiO₂ in water
    psi_mV=-40         # Typical zeta potential at pH 7
)

# Plot effect of ionic strength
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

ax1 = axes[0]
ionic_strengths = [0.001, 0.01, 0.05, 0.1, 0.5]
silica_dlvo.plot_profile(ionic_strengths, ax1)
ax1.set_ylim(-50, 50)

# Plot 2: Stability diagram
ax2 = axes[1]
psi_range = np.linspace(10, 60, 50)
I_range = np.logspace(-3, 0, 50)
psi_mesh, I_mesh = np.meshgrid(psi_range, I_range)

barriers = np.zeros_like(psi_mesh)
for i in range(len(I_range)):
    for j in range(len(psi_range)):
        calc = DLVOCalculator(R_nm=15, A_H=0.83e-20, psi_mV=psi_mesh[i, j])
        result = calc.energy_barrier(I_mesh[i, j])
        barriers[i, j] = result['barrier_kT']

contour = ax2.contourf(psi_mesh, I_mesh * 1000, barriers, levels=20, cmap='RdYlGn')
ax2.contour(psi_mesh, I_mesh * 1000, barriers, levels=[15], colors='black', linewidths=2)
plt.colorbar(contour, ax=ax2, label='Energy Barrier (kT)')
ax2.set_yscale('log')
ax2.set_xlabel('Surface Potential |ψ₀| (mV)', fontsize=11)
ax2.set_ylabel('Ionic Strength (mM)', fontsize=11)
ax2.set_title('DLVO Stability Map (Barrier = 15 kT contour)', fontsize=12, fontweight='bold')

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

# Print stability analysis
print("\n=== DLVO Stability Analysis (30 nm SiO₂, ψ = -40 mV) ===")
print(f"{'I (mM)':>10s} {'λD (nm)':>10s} {'Barrier (kT)':>15s} {'Stable?':>10s}")
print("-" * 50)
for I in ionic_strengths:
    result = silica_dlvo.energy_barrier(I)
    lambda_D = silica_dlvo.debye_length(I) * 1e9
    stable = "Yes" if result['stable'] else "No"
    print(f"{I*1000:>10.0f} {lambda_D:>10.1f} {result['barrier_kT']:>15.1f} {stable:>10s}")

4.3 Dynamic Light Scattering (DLS)

DLS measures particle size by analyzing fluctuations in scattered light intensity caused by Brownian motion.

Example 3: DLS Data Analysis

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

# ===================================
# Example 3: DLS Data Analysis
# ===================================

def stokes_einstein(D, T=298, eta=0.001):
    """
    Calculate hydrodynamic diameter from diffusion coefficient.

    Parameters:
        D: Diffusion coefficient (m²/s)
        T: Temperature (K)
        eta: Viscosity (Pa·s)

    Returns:
        Hydrodynamic diameter (nm)
    """
    kB = 1.38e-23
    d = kB * T / (3 * np.pi * eta * D)
    return d * 1e9  # Convert to nm

def lognormal_distribution(d, d_mean, sigma):
    """Log-normal particle size distribution."""
    return (1 / (d * sigma * np.sqrt(2 * np.pi))) * \
           np.exp(-(np.log(d) - np.log(d_mean))**2 / (2 * sigma**2))

def polydispersity_index(d_mean, d_std):
    """Calculate PDI from mean and standard deviation."""
    return (d_std / d_mean)**2

class DLSAnalyzer:
    """Analyze DLS measurement results."""

    def __init__(self, d_mean_nm, pdi, measurement_type='intensity'):
        """
        Initialize with measured values.

        Parameters:
            d_mean_nm: Z-average diameter (nm)
            pdi: Polydispersity index
            measurement_type: 'intensity', 'volume', or 'number'
        """
        self.d_z = d_mean_nm
        self.pdi = pdi
        self.sigma = np.sqrt(pdi)  # Approximate relationship

    def generate_distribution(self, d_range=None):
        """Generate size distribution."""
        if d_range is None:
            d_range = np.linspace(1, self.d_z * 3, 500)

        distribution = lognormal_distribution(d_range, self.d_z, max(0.1, self.sigma))
        distribution = distribution / np.trapz(distribution, d_range)

        return d_range, distribution

    def intensity_to_number(self, d_range, intensity_dist):
        """Convert intensity distribution to number distribution."""
        # Intensity ∝ d⁶ (Rayleigh scattering)
        number_dist = intensity_dist / d_range**6
        number_dist = number_dist / np.trapz(number_dist, d_range)
        return number_dist

    def quality_assessment(self):
        """Assess dispersion quality from DLS parameters."""
        if self.pdi < 0.1:
            quality = "Excellent (monodisperse)"
            score = 100
        elif self.pdi < 0.2:
            quality = "Good (narrow distribution)"
            score = 80
        elif self.pdi < 0.3:
            quality = "Moderate (polydisperse)"
            score = 60
        elif self.pdi < 0.5:
            quality = "Poor (broad distribution)"
            score = 40
        else:
            quality = "Very poor (aggregated?)"
            score = 20

        return {'quality': quality, 'score': score, 'pdi': self.pdi, 'd_z': self.d_z}

# Simulate DLS results for different samples
samples = {
    'Well-dispersed': {'d_z': 30, 'pdi': 0.08},
    'Moderately dispersed': {'d_z': 45, 'pdi': 0.25},
    'Partially aggregated': {'d_z': 150, 'pdi': 0.45},
    'Heavily aggregated': {'d_z': 500, 'pdi': 0.6}
}

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

# Plot 1: Size distributions (intensity-weighted)
ax1 = axes[0, 0]
for name, params in samples.items():
    analyzer = DLSAnalyzer(params['d_z'], params['pdi'])
    d, dist = analyzer.generate_distribution()
    ax1.plot(d, dist, linewidth=2, label=f"{name} (Z-avg={params['d_z']} nm)")

ax1.set_xlabel('Diameter (nm)', fontsize=11)
ax1.set_ylabel('Intensity (%)', fontsize=11)
ax1.set_title('Intensity-Weighted Size Distributions', fontsize=12, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 300)

# Plot 2: Number vs Intensity distribution comparison
ax2 = axes[0, 1]
analyzer = DLSAnalyzer(50, 0.2)  # Example with moderate PDI
d, intensity_dist = analyzer.generate_distribution()
number_dist = analyzer.intensity_to_number(d, intensity_dist)

ax2.plot(d, intensity_dist / intensity_dist.max(), 'b-', linewidth=2, label='Intensity')
ax2.plot(d, number_dist / number_dist.max(), 'r--', linewidth=2, label='Number')
ax2.axvline(x=50, color='blue', linestyle=':', alpha=0.7)
ax2.axvline(x=d[np.argmax(number_dist)], color='red', linestyle=':', alpha=0.7)
ax2.set_xlabel('Diameter (nm)', fontsize=11)
ax2.set_ylabel('Normalized Intensity', fontsize=11)
ax2.set_title('Intensity vs Number Distribution', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0, 150)

# Plot 3: Quality assessment
ax3 = axes[1, 0]
names = list(samples.keys())
pdi_values = [samples[n]['pdi'] for n in names]
colors = plt.cm.RdYlGn_r(np.array(pdi_values))

bars = ax3.bar(names, pdi_values, color=colors, edgecolor='black')
ax3.axhline(y=0.1, color='green', linestyle='--', alpha=0.7)
ax3.axhline(y=0.3, color='orange', linestyle='--', alpha=0.7)
ax3.axhline(y=0.5, color='red', linestyle='--', alpha=0.7)
ax3.text(3.5, 0.12, 'Monodisperse', fontsize=9, color='green')
ax3.text(3.5, 0.32, 'Acceptable', fontsize=9, color='orange')
ax3.text(3.5, 0.52, 'Aggregated', fontsize=9, color='red')
ax3.set_ylabel('Polydispersity Index (PDI)', fontsize=11)
ax3.set_title('Dispersion Quality Assessment', fontsize=12, fontweight='bold')
ax3.set_xticklabels(names, rotation=15, ha='right', fontsize=9)
ax3.grid(True, alpha=0.3, axis='y')

# Plot 4: Stability monitoring over time
ax4 = axes[1, 1]
times = [0, 1, 7, 14, 30]  # days
stable_sizes = [30, 31, 32, 33, 35]
unstable_sizes = [30, 45, 80, 150, 400]

ax4.semilogy(times, stable_sizes, 'go-', linewidth=2, markersize=8, label='Stable dispersion')
ax4.semilogy(times, unstable_sizes, 'ro-', linewidth=2, markersize=8, label='Unstable dispersion')
ax4.axhline(y=50, color='gray', linestyle='--', alpha=0.7)
ax4.text(25, 55, 'Aggregation threshold', fontsize=9, color='gray')
ax4.set_xlabel('Time (days)', fontsize=11)
ax4.set_ylabel('Z-Average Diameter (nm)', fontsize=11)
ax4.set_title('Stability Monitoring by DLS', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(True, alpha=0.3, which='both')
ax4.set_xlim(0, 32)

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

# Print quality assessment
print("\n=== DLS Quality Assessment ===")
for name, params in samples.items():
    analyzer = DLSAnalyzer(params['d_z'], params['pdi'])
    assessment = analyzer.quality_assessment()
    print(f"{name:25s}: Z-avg = {assessment['d_z']:>5.0f} nm, "
          f"PDI = {assessment['pdi']:.2f}, {assessment['quality']}")

4.4 Sedimentation Analysis

Sedimentation tests directly measure dispersion stability by monitoring particle settling over time.

Example 4: Stokes Sedimentation Model

import numpy as np
import matplotlib.pyplot as plt

# ===================================
# Example 4: Sedimentation Analysis
# ===================================

def stokes_settling_velocity(d_nm, rho_p, rho_f=1000, eta=0.001):
    """
    Calculate Stokes settling velocity.

    Parameters:
        d_nm: Particle diameter (nm)
        rho_p: Particle density (kg/m³)
        rho_f: Fluid density (kg/m³)
        eta: Fluid viscosity (Pa·s)

    Returns:
        Settling velocity (m/s)
    """
    d = d_nm * 1e-9
    g = 9.81
    delta_rho = rho_p - rho_f
    v = (delta_rho * g * d**2) / (18 * eta)
    return v

def settling_time(height_m, d_nm, rho_p, rho_f=1000, eta=0.001):
    """Calculate time for particle to settle given height."""
    v = stokes_settling_velocity(d_nm, rho_p, rho_f, eta)
    if v <= 0:  # Particle floats
        return np.inf
    return height_m / v

def brownian_diffusion_length(d_nm, t_s, T=298, eta=0.001):
    """
    Calculate characteristic Brownian diffusion length.

    Parameters:
        d_nm: Particle diameter (nm)
        t_s: Time (s)
        T: Temperature (K)
        eta: Viscosity (Pa·s)

    Returns:
        RMS displacement (m)
    """
    kB = 1.38e-23
    d = d_nm * 1e-9
    D = kB * T / (3 * np.pi * eta * d)
    return np.sqrt(2 * D * t_s)

# Compare settling for different particle sizes
d_range = np.logspace(0, 4, 100)  # 1 nm to 10 μm
rho_p = 2200  # kg/m³ (silica)
h = 0.05  # 5 cm settling height

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

# Plot 1: Settling velocity vs size
ax1 = axes[0]
v_settle = [stokes_settling_velocity(d, rho_p) for d in d_range]
ax1.loglog(d_range, np.array(v_settle) * 1e9, 'b-', linewidth=2.5, label='Settling velocity')

# Add Brownian velocity (characteristic velocity = sqrt(2D/t) where t = 1 hour)
t_ref = 3600  # 1 hour
v_brownian = [brownian_diffusion_length(d, t_ref) / t_ref for d in d_range]
ax1.loglog(d_range, np.array(v_brownian) * 1e9, 'r--', linewidth=2, label='Brownian velocity (1h)')

# Find crossover
crossover_idx = np.argmin(np.abs(np.array(v_settle) - np.array(v_brownian)))
ax1.axvline(x=d_range[crossover_idx], color='purple', linestyle=':', linewidth=2)
ax1.text(d_range[crossover_idx] * 1.2, 1e-2,
         f'Crossover ~{d_range[crossover_idx]:.0f} nm', fontsize=10, color='purple')

ax1.set_xlabel('Particle Diameter (nm)', fontsize=11)
ax1.set_ylabel('Velocity (nm/s)', fontsize=11)
ax1.set_title('Settling vs Brownian Motion', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3, which='both')
ax1.set_xlim(1, 10000)

# Plot 2: Settling time vs size
ax2 = axes[1]
t_settle = [settling_time(h, d, rho_p) for d in d_range]

ax2.loglog(d_range, np.array(t_settle) / 3600, 'b-', linewidth=2.5)  # Convert to hours

# Add time markers
ax2.axhline(y=1, color='green', linestyle='--', alpha=0.7)
ax2.axhline(y=24, color='orange', linestyle='--', alpha=0.7)
ax2.axhline(y=168, color='red', linestyle='--', alpha=0.7)
ax2.text(5000, 1.5, '1 hour', fontsize=9, color='green')
ax2.text(5000, 30, '1 day', fontsize=9, color='orange')
ax2.text(5000, 200, '1 week', fontsize=9, color='red')

ax2.set_xlabel('Particle Diameter (nm)', fontsize=11)
ax2.set_ylabel('Settling Time (hours)', fontsize=11)
ax2.set_title(f'Time to Settle 5 cm (SiO₂ in water)', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3, which='both')
ax2.set_xlim(1, 10000)
ax2.set_ylim(0.1, 1e6)

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

# Print settling times for key sizes
print("\n=== Sedimentation Analysis (5 cm height, SiO₂ in water) ===")
print(f"{'Diameter (nm)':>15s} {'Settling velocity':>20s} {'Time to settle':>20s}")
print("-" * 60)
for d in [10, 50, 100, 500, 1000, 5000]:
    v = stokes_settling_velocity(d, rho_p)
    t = settling_time(h, d, rho_p)
    if t < 3600:
        time_str = f"{t:.0f} s"
    elif t < 86400:
        time_str = f"{t/3600:.1f} hours"
    elif t < 86400 * 365:
        time_str = f"{t/86400:.1f} days"
    else:
        time_str = f"{t/(86400*365):.1f} years"
    print(f"{d:>15.0f} {v*1e9:>20.3f} nm/s {time_str:>20s}")

4.5 Comprehensive Stability Assessment

Example 5: Multi-Parameter Stability Score

import numpy as np
import matplotlib.pyplot as plt

# ===================================
# Example 5: Comprehensive Stability Assessment
# ===================================

class StabilityAssessment:
    """Comprehensive stability evaluation combining multiple metrics."""

    def __init__(self, sample_name):
        self.name = sample_name
        self.metrics = {}

    def add_zeta_potential(self, zeta_mV):
        """Add zeta potential measurement."""
        abs_zeta = abs(zeta_mV)
        if abs_zeta > 40:
            score = 100
        elif abs_zeta > 30:
            score = 80
        elif abs_zeta > 20:
            score = 50
        else:
            score = 20

        self.metrics['zeta'] = {
            'value': zeta_mV,
            'score': score,
            'weight': 0.3
        }

    def add_dls(self, d_z_nm, pdi, d_target_nm=None):
        """Add DLS measurement."""
        # Size score (closer to target is better)
        if d_target_nm:
            size_ratio = d_z_nm / d_target_nm
            if size_ratio < 1.2:
                size_score = 100
            elif size_ratio < 1.5:
                size_score = 70
            elif size_ratio < 2:
                size_score = 40
            else:
                size_score = 10
        else:
            size_score = 100 if d_z_nm < 100 else (200 - d_z_nm) / 2

        # PDI score
        if pdi < 0.1:
            pdi_score = 100
        elif pdi < 0.2:
            pdi_score = 80
        elif pdi < 0.3:
            pdi_score = 60
        else:
            pdi_score = max(0, 100 - pdi * 100)

        combined_score = 0.5 * size_score + 0.5 * pdi_score

        self.metrics['dls'] = {
            'value': {'d_z': d_z_nm, 'pdi': pdi},
            'score': combined_score,
            'weight': 0.35
        }

    def add_dlvo_barrier(self, barrier_kT):
        """Add DLVO energy barrier."""
        if barrier_kT > 25:
            score = 100
        elif barrier_kT > 15:
            score = 80
        elif barrier_kT > 5:
            score = 40
        else:
            score = 10

        self.metrics['dlvo'] = {
            'value': barrier_kT,
            'score': score,
            'weight': 0.2
        }

    def add_sedimentation(self, settling_time_days):
        """Add sedimentation stability."""
        if settling_time_days > 365:
            score = 100
        elif settling_time_days > 90:
            score = 80
        elif settling_time_days > 30:
            score = 60
        elif settling_time_days > 7:
            score = 40
        else:
            score = 10

        self.metrics['sedimentation'] = {
            'value': settling_time_days,
            'score': score,
            'weight': 0.15
        }

    def overall_score(self):
        """Calculate weighted overall stability score."""
        total_weight = sum(m['weight'] for m in self.metrics.values())
        weighted_sum = sum(m['score'] * m['weight'] for m in self.metrics.values())
        return weighted_sum / total_weight

    def stability_rating(self):
        """Get stability rating."""
        score = self.overall_score()
        if score >= 85:
            return 'Excellent'
        elif score >= 70:
            return 'Good'
        elif score >= 50:
            return 'Moderate'
        else:
            return 'Poor'

    def report(self):
        """Print detailed stability report."""
        print(f"\n=== Stability Assessment: {self.name} ===")
        print("-" * 50)

        for metric_name, data in self.metrics.items():
            print(f"\n{metric_name.upper()}:")
            print(f"  Value: {data['value']}")
            print(f"  Score: {data['score']:.0f}/100")
            print(f"  Weight: {data['weight']*100:.0f}%")

        print("\n" + "=" * 50)
        print(f"OVERALL SCORE: {self.overall_score():.0f}/100")
        print(f"RATING: {self.stability_rating()}")
        print("=" * 50)

# Compare different samples
samples = []

# Sample 1: Well-optimized dispersion
s1 = StabilityAssessment("Optimized SiO₂")
s1.add_zeta_potential(-42)
s1.add_dls(35, 0.08, d_target_nm=30)
s1.add_dlvo_barrier(35)
s1.add_sedimentation(1000)
samples.append(s1)

# Sample 2: Moderate quality
s2 = StabilityAssessment("Moderate TiO₂")
s2.add_zeta_potential(-28)
s2.add_dls(80, 0.22, d_target_nm=50)
s2.add_dlvo_barrier(12)
s2.add_sedimentation(60)
samples.append(s2)

# Sample 3: Poor stability
s3 = StabilityAssessment("Unstable Fe₃O₄")
s3.add_zeta_potential(-15)
s3.add_dls(250, 0.45, d_target_nm=30)
s3.add_dlvo_barrier(3)
s3.add_sedimentation(5)
samples.append(s3)

# Print reports
for sample in samples:
    sample.report()

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

# Plot 1: Radar chart (simplified as bar chart)
ax1 = axes[0]
metric_names = ['Zeta\nPotential', 'DLS\nQuality', 'DLVO\nBarrier', 'Sedimentation']
x = np.arange(len(metric_names))
width = 0.25

for i, sample in enumerate(samples):
    scores = [
        sample.metrics['zeta']['score'],
        sample.metrics['dls']['score'],
        sample.metrics['dlvo']['score'],
        sample.metrics['sedimentation']['score']
    ]
    ax1.bar(x + i * width, scores, width, label=sample.name, alpha=0.8)

ax1.set_xticks(x + width)
ax1.set_xticklabels(metric_names)
ax1.set_ylabel('Score', fontsize=11)
ax1.set_title('Component Stability Scores', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3, axis='y')
ax1.set_ylim(0, 110)

# Plot 2: Overall comparison
ax2 = axes[1]
names = [s.name for s in samples]
overall_scores = [s.overall_score() for s in samples]
ratings = [s.stability_rating() for s in samples]
colors = ['green' if r == 'Excellent' else 'yellow' if r == 'Good'
          else 'orange' if r == 'Moderate' else 'red' for r in ratings]

bars = ax2.barh(names, overall_scores, color=colors, edgecolor='black', alpha=0.8)
ax2.axvline(x=85, color='green', linestyle='--', alpha=0.7)
ax2.axvline(x=70, color='orange', linestyle='--', alpha=0.7)
ax2.axvline(x=50, color='red', linestyle='--', alpha=0.7)

for bar, score, rating in zip(bars, overall_scores, ratings):
    ax2.text(bar.get_width() + 1, bar.get_y() + bar.get_height()/2,
             f'{score:.0f} ({rating})', va='center', fontsize=10)

ax2.set_xlabel('Overall Stability Score', fontsize=11)
ax2.set_title('Comprehensive Stability Comparison', fontsize=12, fontweight='bold')
ax2.set_xlim(0, 110)
ax2.grid(True, alpha=0.3, axis='x')

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

Chapter Summary