Innovative-Process-Applications/roller-compaction-ribbon-density

1"""2IPA Roller Compaction Synthetic Dataset Generator v1.03========================================================4Generates a physically plausible synthetic dataset based on:5  - Johanson rolling theory (Johanson, 1965) for nip angle & pressure6  - Heckel equation for powder densification under pressure7  - Empirical twin-feed-screw effect on ribbon uniformity8 9THIS IS SYNTHETIC EDUCATIONAL DATA. NOT REAL CUSTOMER OR LAB DATA.10Generated by Innovative Process Applications (IPA) for teaching DOE/RSM.11"""12 13import numpy as np14import pandas as pd15 16rng = np.random.default_rng(seed=42)  # reproducible17 18# ---- Physical constants & realistic ranges ----19# Grounded in published roller compaction literature for pharma excipients20# (microcrystalline cellulose / lactose blends on lab-to-pilot scale rolls)21 22N_RUNS = 600  # enough for DOE/RSM teaching, small enough to inspect by eye23 24# Process parameter ranges (realistic for IPA CL-series lab/pilot compactors)25ROLL_FORCE_KN_CM = (2.0, 14.0)     # kN/cm of roll width26ROLL_SPEED_RPM   = (1.0, 12.0)     # roll rotation speed27FEED_SCREW_RPM   = (10.0, 120.0)   # twin feed screw speed28ROLL_GAP_MM      = (1.5, 4.5)      # ribbon thickness target29 30# Material parameters (MCC-lactose blend, representative)31RHO_TRUE = 1.55           # true density, g/cc32RHO_BULK = 0.45           # tapped bulk density, g/cc33HECKEL_K = 0.018          # 1/MPa, Heckel compressibility constant34HECKEL_A = 0.62           # Heckel intercept (related to initial packing)35 36# Machine geometry (representative of CL2x6 class roller compactor)37ROLL_DIAMETER_MM = 150.038ROLL_WIDTH_MM    = 50.039 40def johanson_pressure_mpa(force_kn_per_cm, gap_mm):41    """42    Approximation of peak nip pressure from Johanson rolling theory.43    Pressure scales with force/(contact area), and contact area grows44    with sqrt(R * (gap_entry - gap)).45    """46    # Effective contact length proxy (mm): sqrt(R * delta_gap)47    # Assume entry gap ~ 3x exit gap for MCC-like materials48    contact_len_mm = np.sqrt(ROLL_DIAMETER_MM/2 * 2.0 * gap_mm)49    # Force per unit area: (kN/cm * 10 N/mm) / (contact_len_mm) -> MPa-ish50    pressure = (force_kn_per_cm * 100.0) / contact_len_mm51    return pressure  # MPa52 53def heckel_relative_density(pressure_mpa):54    """Heckel equation: ln(1/(1-D)) = K*P + A  =>  D = 1 - exp(-(K*P + A))"""55    return 1.0 - np.exp(-(HECKEL_K * pressure_mpa + HECKEL_A))56 57def feed_ratio_effect(feed_rpm, roll_rpm):58    """59    Twin feed screw effect: under-feeding starves the nip (low density,60    high variability); over-feeding causes slip and pre-compaction.61    Optimal ratio roughly 8-15 for twin screws on MCC.62    """63    ratio = feed_rpm / roll_rpm64    # Gaussian-like response centered at ratio=1165    optimality = np.exp(-((ratio - 11.0) ** 2) / (2 * 6.0**2))66    return optimality  # 0..167 68# ---- Generate runs ----69roll_force = rng.uniform(*ROLL_FORCE_KN_CM, N_RUNS)70roll_speed = rng.uniform(*ROLL_SPEED_RPM, N_RUNS)71feed_speed = rng.uniform(*FEED_SCREW_RPM, N_RUNS)72roll_gap   = rng.uniform(*ROLL_GAP_MM, N_RUNS)73 74# Physics75peak_pressure = johanson_pressure_mpa(roll_force, roll_gap)76base_rel_density = heckel_relative_density(peak_pressure)77feed_opt = feed_ratio_effect(feed_speed, roll_speed)78 79# Combine: good feed ratio lets you reach base density; poor ratio penalizes80relative_density = base_rel_density * (0.85 + 0.15 * feed_opt)81# Add realistic measurement noise (~1.5% relative)82relative_density += rng.normal(0, 0.015, N_RUNS)83relative_density = np.clip(relative_density, 0.40, 0.95)84 85ribbon_density = relative_density * RHO_TRUE               # g/cc86porosity = 1.0 - relative_density                          # fraction87 88# Ribbon uniformity: twin feed screws give lower CV across roll width89# when feed ratio is in the sweet spot; degrades otherwise90density_cv_pct = (2.0 + 4.5 * (1 - feed_opt)) + rng.normal(0, 0.4, N_RUNS)91density_cv_pct = np.clip(density_cv_pct, 1.0, 10.0)92 93throughput_kg_hr = (94    roll_speed * roll_gap * ROLL_WIDTH_MM * ribbon_density * 60 / 1000 * 0.8595) + rng.normal(0, 1.5, N_RUNS)96throughput_kg_hr = np.clip(throughput_kg_hr, 0, None)97 98df = pd.DataFrame({99    "run_id": np.arange(1, N_RUNS + 1),100    "roll_force_kN_per_cm": np.round(roll_force, 2),101    "roll_speed_rpm": np.round(roll_speed, 2),102    "feed_screw_rpm": np.round(feed_speed, 1),103    "roll_gap_mm": np.round(roll_gap, 2),104    "peak_pressure_MPa": np.round(peak_pressure, 1),105    "ribbon_rel_density": np.round(relative_density, 4),106    "ribbon_density_g_cc": np.round(ribbon_density, 4),107    "ribbon_porosity": np.round(porosity, 4),108    "density_CV_percent": np.round(density_cv_pct, 2),109    "throughput_kg_hr": np.round(throughput_kg_hr, 2),110})111 112df.to_csv("ribbon_density_v1.0.csv", index=False)113print(f"Wrote {len(df)} rows")114print(df.describe().round(3))115