""" Rheology intro animation: Newtonian fluid vs. shear-thinning polymer melt (PBT-like). Companion to entries/mathematics/rheology-pbt-extrusion.html. Run: python3 create_rheology_gif.py (writes GIF / PNG / optional MP4 alongside this file) ------------------------------------------------------------------------------ The physics this animation teaches ------------------------------------------------------------------------------ Rheology studies how materials deform and flow under applied stress. In simple shear a material is sheared between two plates; we track three quantities: tau shear stress [Pa] force per unit area gdot shear rate [1/s] du/dy, the velocity gradient eta viscosity [Pa s] resistance to flow, eta = tau / gdot Newtonian fluid (water, simple solvents at fixed T): tau = eta * gdot, eta = constant. So tau grows LINEARLY with gdot and the ratio tau/gdot never changes. Shear-thinning polymer melt (Ostwald-de Waele power law): tau = K * gdot**n, eta_app(gdot) = tau / gdot = K * gdot**(n-1). For 0 < n < 1 the exponent n-1 < 0, so apparent viscosity DECREASES as the shear rate increases. Molecularly: long entangled chains orient with the flow at high shear, offering less resistance. This is a POWER-LAW drop, not an exponential one. Temperature: viscosity also drops with temperature. We use a normalized Arrhenius-type factor about a reference temperature T_ref, eta_app(gdot, T) = K_ref * exp[ (Ea/R) * (1/T - 1/T_ref) ] * gdot**(n-1). This is a PEDAGOGICAL model. The constants below are illustrative, not a calibrated fit to any real PBT grade. """ import numpy as np import matplotlib matplotlib.use("Agg") # headless: render frames without a display import matplotlib.pyplot as plt from matplotlib.patches import FancyArrowPatch, Rectangle, Circle, FancyBboxPatch from matplotlib.animation import FuncAnimation, PillowWriter try: from matplotlib.animation import FFMpegWriter except Exception: # pragma: no cover - extremely rare FFMpegWriter = None from pathlib import Path # -------------------------------------------------------------------------- # Pedagogical material parameters (NOT a calibrated PBT fit) # -------------------------------------------------------------------------- K_ref = 1000.0 # Pa s^n, consistency index at T_ref n = 0.45 # shear-thinning index, 0 < n < 1 eta_newtonian = 300.0 # Pa s, constant Newtonian viscosity Ea = 50_000.0 # J/mol, apparent activation energy for flow R = 8.314 # J/(mol K), gas constant T_ref = 533.15 # K (~260 C), reference temperature T_cold = 513.15 # K (~240 C) T_hot = 553.15 # K (~280 C) # Shear-rate axis: log-spaced over [1e-1, 1e3] s^-1 GDOT = np.logspace(-1, 3, 400) # -------------------------------------------------------------------------- # Rheological model functions # -------------------------------------------------------------------------- def eta_polymer(gdot, T): """Apparent viscosity of the polymer melt, eta_app(gdot, T) [Pa s]. Power-law shear thinning times a normalized Arrhenius temperature factor. """ arrhenius = np.exp(Ea / R * (1.0 / T - 1.0 / T_ref)) return K_ref * arrhenius * gdot ** (n - 1.0) def tau_polymer(gdot, T): """Shear stress of the polymer melt, tau = eta_app * gdot [Pa].""" return eta_polymer(gdot, T) * gdot def tau_newtonian(gdot): """Newtonian shear stress, tau = eta * gdot [Pa].""" return eta_newtonian * gdot # -------------------------------------------------------------------------- # Color palette (<= 4 main colors) and global style # -------------------------------------------------------------------------- C_NEWT = "#1f6fb2" # Newtonian: blue C_POLY = "#e8820c" # polymer melt (reference T): orange C_COLD = "#5e3c99" # cold curve: purple C_HOT = "#d7301f" # hot curve: red C_INK = "#222222" # text / structure C_FAINT = "#888888" plt.rcParams.update({ "font.size": 15, "font.family": "DejaVu Sans", "axes.titlesize": 17, "axes.labelsize": 16, "axes.edgecolor": "#444444", "xtick.labelsize": 12, "ytick.labelsize": 12, "mathtext.fontset": "cm", "figure.facecolor": "white", "axes.facecolor": "white", "savefig.facecolor": "white", }) # -------------------------------------------------------------------------- # Frame / scene bookkeeping # -------------------------------------------------------------------------- fps = 12 duration_sec = 16 n_frames = fps * duration_sec # 192 frames scene_edges = [0, 36, 72, 120, 156, 192] # scene boundaries in frames def scene_progress(frame, start, end): """Fractional progress in [0,1] of `frame` within scene [start, end).""" return np.clip((frame - start) / (end - start), 0.0, 1.0) def smoothstep(x): """Smooth 0->1 ease (zero slope at both ends).""" x = np.clip(x, 0.0, 1.0) return x * x * (3.0 - 2.0 * x) # -------------------------------------------------------------------------- # Figure scaffold: a single 1280x720 canvas reused for every frame. # We use one full-canvas axes for drawings/text and add plot axes per scene. # -------------------------------------------------------------------------- DPI = 100 fig = plt.figure(figsize=(12.8, 7.2), dpi=DPI) # Background axes covering the whole figure for titles, schematics, labels. bg = fig.add_axes([0, 0, 1, 1]) bg.set_xlim(0, 1) bg.set_ylim(0, 1) bg.axis("off") def _text(x, y, s, size=16, color=C_INK, ha="left", va="center", weight="normal", style="normal", alpha=1.0): """Convenience wrapper for background-axes text in figure coordinates.""" return bg.text(x, y, s, size=size, color=color, ha=ha, va=va, weight=weight, style=style, alpha=alpha, zorder=10) # -------------------------------------------------------------------------- # Helper: simple-shear cell (two plates + sliding fluid layers) # -------------------------------------------------------------------------- def draw_shear_cell(ax, progress, mode="newtonian"): """Draw a sheared fluid between a fixed bottom plate and a moving top plate. progress in [0,1] animates the top plate sliding right and the per-layer velocity arrows growing (longer near the top: a linear velocity gradient). `mode` only sets the accent color of the fluid block. """ ax.clear() ax.set_xlim(0, 10) ax.set_ylim(0, 10) ax.axis("off") accent = C_NEWT if mode == "newtonian" else C_POLY x0, x1 = 1.2, 8.0 y_bot, y_top = 2.0, 8.0 shift = 1.6 * progress # how far the top plate has slid # Fluid block (light fill) -- the sheared material. ax.add_patch(Rectangle((x0, y_bot), x1 - x0, y_top - y_bot, facecolor=accent, alpha=0.10, edgecolor="none", zorder=1)) # Fixed bottom plate (hatched bar) and moving top plate. ax.add_patch(Rectangle((x0 - 0.3, y_bot - 0.45), x1 - x0 + 0.6, 0.45, facecolor="#cfcfcf", edgecolor=C_INK, zorder=3)) ax.add_patch(Rectangle((x0 - 0.3 + shift, y_top), x1 - x0 + 0.6, 0.45, facecolor="#cfcfcf", edgecolor=C_INK, zorder=3)) # Fluid layers: horizontal lines that skew with height (top moves most). n_layers = 7 for i in range(n_layers): frac = i / (n_layers - 1) # 0 at bottom, 1 at top y = y_bot + frac * (y_top - y_bot) dx = shift * frac # linear velocity profile ax.plot([x0 + dx, x1 + dx], [y, y], color=accent, lw=1.4, alpha=0.55, zorder=2) # Velocity arrow: length grows with height -> shear rate = du/dy. arrow_len = 0.25 + 2.4 * frac * (0.4 + 0.6 * progress) if arrow_len > 0.05: ax.add_patch(FancyArrowPatch( (x1 + dx + 0.15, y), (x1 + dx + 0.15 + arrow_len, y), arrowstyle="-|>", mutation_scale=12, color=accent, lw=2.0, zorder=4)) # Annotate the moving / fixed plates. ax.text((x0 + x1) / 2 + shift, y_top + 0.7, "top plate moves →", ha="center", va="bottom", fontsize=12.5, color=C_INK) ax.text((x0 + x1) / 2, y_bot - 0.95, "fixed bottom plate", ha="center", va="top", fontsize=12.5, color=C_INK) # Velocity-gradient bracket label. ax.text(x1 + 1.9, (y_bot + y_top) / 2, r"$u(y)$" "\n" r"$\dot{\gamma}=du/dy$", ha="left", va="center", fontsize=13, color=accent) # -------------------------------------------------------------------------- # Helper: cartoon polymer chains, oriented by `alignment` in [0,1] # -------------------------------------------------------------------------- def draw_polymer_chains(ax, alignment, seed=0): """Draw squiggly polymer chains; alignment 0 = tangled, 1 = flow-aligned. Each chain is a sine squiggle with a random phase. At low alignment the chains take random orientations; at high alignment they rotate toward the horizontal (flow) direction and stretch out, the molecular picture behind shear thinning. """ ax.clear() ax.set_xlim(0, 10) ax.set_ylim(0, 10) ax.axis("off") a = smoothstep(alignment) rng = np.random.default_rng(seed) # fixed seed -> stable chains across frames n_chains = 9 t = np.linspace(-1.0, 1.0, 60) for _ in range(n_chains): cx = rng.uniform(1.8, 8.2) cy = rng.uniform(2.0, 8.0) length = (1.1 + 0.4 * rng.random()) * (1.0 + 0.9 * a) # stretch out amp = (0.45 + 0.35 * rng.random()) * (1.0 - 0.75 * a) # uncrinkle phase = rng.uniform(0, 2 * np.pi) waves = rng.uniform(2.0, 3.5) # Local chain shape: a stretched axis (xs) + transverse squiggle (ys). xs = length * t ys = amp * np.sin(waves * np.pi * t + phase) # Orientation: random when tangled, -> horizontal as it aligns. theta0 = rng.uniform(0, 2 * np.pi) theta = (1.0 - a) * theta0 # collapse toward 0 rad (flow direction) ct, st = np.cos(theta), np.sin(theta) X = cx + ct * xs - st * ys Y = cy + st * xs + ct * ys ax.plot(X, Y, color=C_POLY, lw=2.2, alpha=0.85, solid_capstyle="round") # Flow-direction velocity arrows, growing with alignment (= shear rate). for yk in (3.0, 5.0, 7.0): L = 0.6 + 2.4 * a ax.add_patch(FancyArrowPatch((0.6, yk), (0.6 + L, yk), arrowstyle="-|>", mutation_scale=12, color=C_FAINT, lw=1.8)) # -------------------------------------------------------------------------- # Helper: simplified twin-screw extruder schematic # -------------------------------------------------------------------------- def draw_extruder(ax, progress): """Two intermeshing rotating screws in a barrel; melt moves left to right.""" ax.clear() ax.set_xlim(0, 10) ax.set_ylim(0, 10) ax.axis("off") bx0, bx1, by0, by1 = 0.6, 9.4, 3.2, 6.8 # Barrel outline. ax.add_patch(FancyBboxPatch((bx0, by0), bx1 - bx0, by1 - by0, boxstyle="round,pad=0.02,rounding_size=0.25", facecolor="#f1f1ec", edgecolor=C_INK, lw=1.8, zorder=1)) # Two screw shafts drawn as a series of flights (slanted ticks) that # appear to rotate/convey as `progress` advances. phase = 2 * np.pi * progress for cy in (by0 + 1.0, by1 - 1.0): ax.plot([bx0 + 0.4, bx1 - 0.4], [cy, cy], color=C_FAINT, lw=1.2, zorder=2) nf = 16 for k in range(nf): xf = bx0 + 0.6 + (bx1 - bx0 - 1.2) * k / (nf - 1) s = 0.5 * np.sin(phase + k * 0.9) # flight tilt animates ax.plot([xf - 0.18, xf + 0.18], [cy - 0.45 + s * 0.25, cy + 0.45 + s * 0.25], color=C_POLY, lw=2.4, solid_capstyle="round", zorder=3) # Melt moving left -> right: a few markers advancing with progress. for j in range(5): mx = bx0 + 0.8 + ((j / 5.0 + progress) % 1.0) * (bx1 - bx0 - 1.6) ax.add_patch(Circle((mx, (by0 + by1) / 2), 0.16, facecolor=C_POLY, edgecolor="none", alpha=0.7, zorder=4)) ax.text((bx0 + bx1) / 2, by1 + 0.5, "twin-screw extruder: two intermeshing screws convey the melt →", ha="center", va="bottom", fontsize=12.5, color=C_INK) # -------------------------------------------------------------------------- # Helper: clean up a Matplotlib plot axis (spines, grid) # -------------------------------------------------------------------------- def setup_clean_axis(ax): ax.spines["top"].set_visible(False) ax.spines["right"].set_visible(False) ax.grid(True, which="both", color="#dddddd", lw=0.6, alpha=0.7) ax.set_axisbelow(True) def save_last_frame(fig, filename): """Save the current figure state as a PNG poster of the final frame.""" fig.savefig(filename, dpi=DPI) # -------------------------------------------------------------------------- # Per-scene renderers. Each receives local progress p in [0,1]. # We create/destroy plot axes per frame so the single canvas can morph. # -------------------------------------------------------------------------- # Track dynamically-added axes so we can clear them each frame. _dyn_axes = [] def _clear_dynamic(): for ax in _dyn_axes: ax.remove() _dyn_axes.clear() bg.clear() bg.set_xlim(0, 1) bg.set_ylim(0, 1) bg.axis("off") def _header(title, subtitle=None): _text(0.5, 0.95, title, size=22, ha="center", weight="bold") if subtitle: _text(0.5, 0.895, subtitle, size=14, ha="center", color=C_FAINT) def scene1(p): """What is shear? Animated plates + growing velocity arrows.""" _header("1. What is shear?") ax = fig.add_axes([0.05, 0.08, 0.62, 0.74]); _dyn_axes.append(ax) draw_shear_cell(ax, smoothstep(p), mode="newtonian") _text(0.70, 0.70, "Shear = neighboring layers\nsliding past one another", size=15) _text(0.70, 0.50, r"Shear rate $\dot{\gamma}$: speed gradient", size=15, color=C_NEWT) _text(0.70, 0.38, r"Shear stress $\tau$: force per area", size=15, color=C_NEWT) _text(0.70, 0.22, r"$\eta = \tau/\dot{\gamma}$ (viscosity)", size=16) def scene2(p): """Newtonian fluid: tau = eta * gdot, linear, marker sweeps the line.""" _header("2. Newtonian fluid", r"$\tau = \eta\,\dot{\gamma}$, $\eta$ constant") # Left: shear cell. axL = fig.add_axes([0.03, 0.08, 0.40, 0.70]); _dyn_axes.append(axL) draw_shear_cell(axL, 1.0, mode="newtonian") # Right: tau vs gdot on LINEAR axes. axR = fig.add_axes([0.55, 0.16, 0.40, 0.64]); _dyn_axes.append(axR) setup_clean_axis(axR) g = np.linspace(0, 1000, 200) axR.plot(g, tau_newtonian(g), color=C_NEWT, lw=3) axR.set_xlabel(r"$\dot{\gamma}\ \mathrm{(s^{-1})}$") axR.set_ylabel(r"$\tau\ \mathrm{(Pa)}$") axR.set_xlim(0, 1000) axR.set_ylim(0, tau_newtonian(1000) * 1.05) # Marker sweeping along the line as gdot increases. gm = 1000 * smoothstep(p) axR.plot([gm], [tau_newtonian(gm)], "o", color=C_NEWT, ms=11, markeredgecolor="white") axR.text(0.05, 0.90, r"$\eta=\tau/\dot{\gamma}=300\ \mathrm{Pa\,s}$ (fixed)", transform=axR.transAxes, fontsize=13, color=C_NEWT) _text(0.5, 0.045, "Double the shear rate, double the stress. The ratio stays fixed.", size=15, ha="center") def scene3(p): """Polymer melt: chains align (left) + log-log eta drop (right).""" _header("3. Polymer melt: chains align, viscosity drops") a = smoothstep(p) # alignment / marker position driver # Left: polymer chains becoming aligned. axL = fig.add_axes([0.03, 0.10, 0.40, 0.66]); _dyn_axes.append(axL) draw_polymer_chains(axL, a, seed=3) if a < 0.5: _text(0.23, 0.085, "low shear: chains tangled;\nhigh apparent viscosity", size=13, ha="center", color=C_POLY) else: _text(0.23, 0.085, "higher shear: chains align;\napparent viscosity drops", size=13, ha="center", color=C_POLY) # Right: eta_app vs gdot on LOG-LOG axes. axR = fig.add_axes([0.55, 0.17, 0.40, 0.60]); _dyn_axes.append(axR) setup_clean_axis(axR) axR.set_xscale("log"); axR.set_yscale("log") axR.plot(GDOT, eta_polymer(GDOT, T_ref), color=C_POLY, lw=3, label="polymer melt") axR.plot(GDOT, np.full_like(GDOT, eta_newtonian), color=C_NEWT, lw=2.5, ls="--", label="Newtonian") axR.set_xlabel(r"$\dot{\gamma}\ \mathrm{(s^{-1})}$") axR.set_ylabel(r"$\eta_{\mathrm{app}}\ \mathrm{(Pa\cdot s)}$") axR.set_ylim(10, 1e4) # Direct labels (no tiny legend). axR.text(2e2, eta_newtonian * 1.25, "Newtonian $\\eta=300$", color=C_NEWT, fontsize=12, ha="center") axR.text(2e0, eta_polymer(2e0, T_ref) * 1.5, r"$\eta_{\mathrm{app}}=K\dot{\gamma}^{\,n-1}$", color=C_POLY, fontsize=13) # Marker sweeping from low to high gdot along the polymer curve. gm = 10 ** (-1 + 4 * a) axR.plot([gm], [eta_polymer(gm, T_ref)], "o", color=C_POLY, ms=11, markeredgecolor="white", zorder=5) _text(0.5, 0.05, r"$0