Add Herschel-Bulkley non-Newtonian viscosity model

docs/documentation/case.md

@@ -1091,7 +1091,48 @@ This boundary condition can be used for fixed-temperature (isothermal) walls at
 
 
 
-### 19. GPU Performance (NVIDIA UVM)
+### 19. Non-Newtonian (Herschel-Bulkley) Viscosity {#sec-non-newtonian}
+
+| Parameter                         | Type    | Description                                                          |
+| ---:                              | :----:  | :---                                                                 |
+| `fluid_pp(i)%%non_newtonian`      | Logical | Enable Herschel-Bulkley non-Newtonian viscosity for fluid \f$i\f$.  |
+| `fluid_pp(i)%%K`                  | Real    | Consistency index \f$K\f$.                                          |
+| `fluid_pp(i)%%nn`                 | Real    | Flow index \f$n\f$ (\f$n<1\f$ shear-thinning, \f$n>1\f$ shear-thickening). |
+| `fluid_pp(i)%%tau0`               | Real    | Yield stress \f$\tau_0\f$; set to 0 for pure power-law.             |
+| `fluid_pp(i)%%hb_m`               | Real    | Papanastasiou regularization parameter \f$m\f$; required when `tau0 > 0`. |
+| `fluid_pp(i)%%mu_min`             | Real    | Lower viscosity clamp \f$\mu_{\min}\f$.                              |
+| `fluid_pp(i)%%mu_max`             | Real    | Upper viscosity clamp \f$\mu_{\max}\f$ (required).                  |
+| `fluid_pp(i)%%mu_bulk`            | Real    | Reserved; non-Newtonian bulk viscosity is not yet supported. The validator rejects it on a non-Newtonian fluid; on a Newtonian fluid it is accepted and ignored. |
+
+The effective dynamic viscosity is computed from the Papanastasiou-regularized Herschel-Bulkley model:
+
+\f[
+\mu_{\rm eff}(\dot\gamma) = \frac{\tau_0}{\dot\gamma}\!\left(1 - e^{-m\,\dot\gamma}\right) + K\,\dot\gamma^{n-1},
+\qquad
+\dot\gamma = \sqrt{2\,D_{ij}D_{ij}},
+\f]
+
+where \f$D_{ij} = \frac{1}{2}(\partial_i u_j + \partial_j u_i)\f$ is the strain-rate tensor and \f$\dot\gamma\f$ is the scalar shear rate.
+The result is clamped to \f$[\mu_{\min},\,\mu_{\max}]\f$.
+
+Special cases:
+
+- `tau0 = 0`: pure power-law fluid, \f$\mu_{\rm eff} = K\,\dot\gamma^{n-1}\f$.
+- `tau0 = 0`, `nn = 1`: Newtonian fluid with constant viscosity \f$\mu = K\f$.
+- `tau0 > 0`, `nn = 1`: Bingham plastic.
+
+Usage notes:
+
+- Requires `viscous = T`. `fluid_pp(i)%%Re(1)` must be set (use `1.0/K` to register the fluid as viscous; the HB model overrides \f$\mu_{\rm eff}\f$ cell-by-cell). `fluid_pp(i)%%Re(2)` (bulk viscosity) must not be set for a non-Newtonian fluid.
+- `mu_max` is required; `mu_min` is inactive if omitted (no lower clamp applied).
+- Positivity requirements: `K`, `nn`, and `mu_max` must be positive; `mu_min` and `hb_m` must be positive when set; `tau0` must be non-negative.
+- Requires `model_eqns = 2` or `3` and is incompatible with `igr`.
+- Supported only with `riemann_solver = 1` (HLL) or `riemann_solver = 2` (HLLC).
+- The HB parameters (`K`, `nn`, `tau0`, `hb_m`, `mu_min`, `mu_max`, `mu_bulk`) may only be set on a fluid with `non_newtonian = T`; the validator rejects them otherwise.
+- All HB parameters are non-dimensional (scaled by \f$\rho_{\rm ref} U_{\rm ref} L_{\rm ref}\f$), so \f$1/\mu_{\rm eff}\f$ equals the local effective Reynolds number.
+- For cylindrical geometry (`cyl_coord = T`) the shear rate uses the grid-direction strain components; curvature corrections to \f$\dot\gamma\f$ are not yet included.
+
+### 20. GPU Performance (NVIDIA UVM)
 
 | Parameter                  | Type    | Description                                              |
 | ---:                       | :---:   | :---                                                     |

docs/module_categories.json

@@ -15,6 +15,7 @@
         "category": "Physics Models",
         "modules": [
             "m_viscous",
+            "m_hb_function",
             "m_surface_tension",
             "m_bubbles",
             "m_bubbles_EE",

docs/references.bib

@@ -659,3 +659,14 @@ @article{Xu2019
   year    = {2019},
   doi     = {10.1063/1.5116374}
 }
+
+@article{Papanastasiou87,
+  author  = {Papanastasiou, Tasos C.},
+  title   = {Flows of Materials with Yield},
+  journal = {Journal of Rheology},
+  volume  = {31},
+  number  = {5},
+  pages   = {385--404},
+  year    = {1987},
+  doi     = {10.1122/1.549926}
+}

examples/2D_bingham_poiseuille_nn/README.md

@@ -0,0 +1,52 @@
+# 2D Bingham (Yield-Stress) Poiseuille Channel
+
+Validates the **yield-stress term** of MFC's Herschel-Bulkley non-Newtonian
+viscosity against a closed-form analytic Poiseuille profile. Demonstrates the
+Bingham regime (`nn = 1`, `tau0 > 0`): a rigid **plug** of uniform velocity forms
+near the centerline, where the shear stress falls below the yield stress.
+
+## Regime and parameters
+
+Single Papanastasiou-regularized Herschel-Bulkley fluid with unit flow index, so
+`K = mu` is a plain Newtonian consistency and the only non-Newtonian effect is the
+yield stress:
+
+| Parameter | Value | Role |
+|-----------|-------|------|
+| `fluid_pp(1)%K`   | `5.0e-2` | `n = 1` -> plain dynamic viscosity `mu` |
+| `fluid_pp(1)%nn`  | `1.0`    | flow index = 1 (Bingham, no power-law) |
+| `fluid_pp(1)%tau0`| `4.0e-3` | yield stress -> plug half-width `y0 = 0.4 H` |
+| `fluid_pp(1)%hb_m`| `1.0e4`  | sharp Papanastasiou yield regularization |
+| `fluid_pp(1)%mu_min`/`mu_max` | `1e-6` / `1.0` | viscosity clamp (rigid plug) |
+
+Driven by `g_x = 0.1`, `rho = 1`, `pres = 10`, giving `tau_w = rho*g*H = 1e-2`,
+`u_plug ~ 3.6e-3` and Mach ~1e-3. Channel `L_y = 0.2`, `H = 0.1`, no-slip walls,
+periodic in `x`.
+
+## Governing physics and analytic solution
+
+The shear stress is `tau = rho*g*(H - y)`; the fluid only flows where `|tau| > tau0`.
+With `n = 1`, `K = mu`, `tau_w = rho*g*H > tau0`:
+
+    plug half-width   : y0   = tau0/(rho*g)
+    sheared region    : u(y) = (1/(2*mu*rho*g)) *
+                               [ (tau_w - tau0)^2 - (rho*g*(H-y) - tau0)^2 ]   (|y-H| >= y0)
+    plug              : u_plug = (1/(2*mu*rho*g)) * (tau_w - tau0)^2            (|y-H| <  y0)
+
+The signature of a correct yield term is the flat plug of uniform velocity within
+`|y - H| < y0`. Requires `tau_w > tau0` for any flow.
+
+## How to run
+
+    ./mfc.sh run examples/2D_bingham_poiseuille_nn/case.py -n 2
+    python examples/2D_bingham_poiseuille_nn/compare_analytic.py
+
+## Validation result
+
+Relative L2 error vs. the analytic Bingham profile: **2.5%** (2-rank run, steady
+state confirmed). A plug forms at the centerline with half-width matching the
+analytic `y0 = tau0/(rho*g) = 0.4 H`.
+
+## References
+
+- Papanastasiou, T. C. (1987). Flows of materials with yield. *J. Rheol.* 31, 385.

examples/2D_bingham_poiseuille_nn/case.py

@@ -0,0 +1,132 @@
+#!/usr/bin/env python3
+"""
+2D Bingham (Herschel-Bulkley with flow index n = 1, yield stress tau0 > 0)
+Poiseuille channel.
+
+Validation case for the YIELD-STRESS term of the non-Newtonian (Herschel-Bulkley)
+viscosity. The companion examples/2D_poiseuille_nn validates the power-law term
+(tau0 = 0); this case isolates the yield term by fixing n = 1 (so K = mu is a plain
+Newtonian consistency) and turning on tau0 > 0. The signature of a correct yield
+term is a rigid PLUG of uniform velocity near the centerline, where the shear
+stress |tau| = rho*g*(H - y) drops below tau0 and the fluid stops yielding.
+
+A constant body acceleration g_x drives a fully-developed channel flow between two
+stationary no-slip walls (y = 0, L_y; half-height H = L_y/2, centerline y = H).
+The steady Bingham profile (n = 1, K = mu, tau_w = rho*g*H):
+
+  plug half-width from centerline   : y0   = tau0/(rho*g)
+  sheared region (0 <= y <= H - y0) : u(y) = (1/(2*mu*rho*g)) *
+                                              [ (tau_w - tau0)^2 - (rho*g*(H-y) - tau0)^2 ]
+  plug (H - y0 <= y <= H)           : u_plug = (1/(2*mu*rho*g)) * (tau_w - tau0)^2
+  upper half mirrors about y = H.
+
+Requires tau_w = rho*g*H > tau0 for any flow.
+
+Parameters (nondimensional MFC units):
+    rho  = 1.0, H = 0.1 (L_y = 0.2)
+    g_x  = 0.1          -> tau_w = rho*g*H = 1.0e-2
+    tau0 = 4.0e-3       -> y0 = tau0/(rho*g) = 4.0e-2 = 0.4 H  (clear plug + shear)
+    K    = mu = 5.0e-2  (n = 1 Bingham consistency; short diffusive time t_d = H^2 rho/mu = 0.2)
+    hb_m = 1.0e4        (sharp Papanastasiou yield; uncapped plug mu = tau0*hb_m = 40)
+    mu_min = 1e-6, mu_max = 1.0  (plug viscosity = 20*K, effectively rigid plug;
+          kept modest because the explicit viscous timestep scales as 1/mu_max)
+    pres = 10 -> sound speed ~3.74; u_plug ~ 3.6e-3 => Mach ~1e-3 (low Mach)
+    grid: m = 24 (x, periodic; >= 24 for a 2-rank y-split WENO5 decomposition),
+          n = 63 (y resolution of the plug; the viscous CFL dt ~ dy^2 rho/mu_max
+          makes a finer grid prohibitively slow for an explicit solver), L_x = L_y = 0.2
+    cfl_adap_dt, cfl_target = 0.3, t_stop = 1.5 (~7.5 diffusive times t_d = 0.2)
+
+Compare with compare_analytic.py.
+"""
+
+import json
+
+# Channel / fluid parameters
+L_x = 0.2
+L_y = 0.2
+rho = 1.0
+pres = 10.0
+K = 5.0e-2  # n = 1 -> K is the plain dynamic viscosity mu
+nn = 1.0
+tau0 = 4.0e-3
+g_x = 0.1
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0.0,
+            "x_domain%end": L_x,
+            "y_domain%beg": 0.0,
+            "y_domain%end": L_y,
+            "m": 24,
+            "n": 63,
+            "p": 0,
+            "cfl_adap_dt": "T",
+            "cfl_target": 0.3,
+            "n_start": 0,
+            "t_save": 0.15,
+            "t_stop": 1.5,
+            # Simulation Algorithm Parameters
+            "num_patches": 1,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 1,
+            "mpp_lim": "F",
+            "mixture_err": "T",
+            "time_stepper": 3,
+            "weno_order": 5,
+            "weno_eps": 1e-16,
+            "mapped_weno": "T",
+            "weno_Re_flux": "T",
+            "mp_weno": "T",
+            "weno_avg": "T",
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            # x: periodic channel; y: stationary no-slip walls
+            "bc_x%beg": -1,
+            "bc_x%end": -1,
+            "bc_y%beg": -16,
+            "bc_y%end": -16,
+            "viscous": "T",
+            # Constant body acceleration in +x: accel = g_x + k_x*sin(w_x*t - p_x)
+            "bf_x": "T",
+            "g_x": g_x,
+            "k_x": 0.0,
+            "w_x": 0.0,
+            "p_x": 0.0,
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "fd_order": 4,
+            "parallel_io": "T",
+            # Patch 1: full domain, quiescent
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%x_centroid": 0.5 * L_x,
+            "patch_icpp(1)%y_centroid": 0.5 * L_y,
+            "patch_icpp(1)%length_x": L_x,
+            "patch_icpp(1)%length_y": L_y,
+            "patch_icpp(1)%vel(1)": 0.0,
+            "patch_icpp(1)%vel(2)": 0.0,
+            "patch_icpp(1)%pres": pres,
+            "patch_icpp(1)%alpha_rho(1)": rho,
+            "patch_icpp(1)%alpha(1)": 1.0,
+            # Fluids Physical Parameters: single Bingham (HB with n = 1, tau0 > 0) fluid
+            "fluid_pp(1)%gamma": 1.0 / (1.4 - 1.0),
+            "fluid_pp(1)%pi_inf": 0.0,
+            "fluid_pp(1)%Re(1)": 1.0 / K,
+            "fluid_pp(1)%non_newtonian": "T",
+            "fluid_pp(1)%K": K,
+            "fluid_pp(1)%nn": nn,
+            "fluid_pp(1)%tau0": tau0,
+            "fluid_pp(1)%hb_m": 1.0e4,
+            "fluid_pp(1)%mu_min": 1e-6,
+            "fluid_pp(1)%mu_max": 1.0,
+        }
+    )
+)

examples/2D_bingham_poiseuille_nn/compare_analytic.py

@@ -0,0 +1,121 @@
+#!/usr/bin/env python3
+"""
+Compare the steady MFC velocity profile of the 2D Bingham Poiseuille channel
+(examples/2D_bingham_poiseuille_nn/case.py) against the closed-form analytic
+plug-flow solution, and measure the numerical plug half-width.
+
+For a Bingham fluid (Herschel-Bulkley with flow index n = 1, consistency K = mu,
+yield stress tau0) driven by body acceleration g in a channel of height L_y
+(half-height H = L_y/2, centerline y = H), with tau_w = rho*g*H > tau0:
+
+  plug half-width  : y0   = tau0/(rho*g)
+  sheared region   : u(y) = (1/(2*mu*rho*g)) *
+                            [ (tau_w - tau0)^2 - (rho*g*(H-y) - tau0)^2 ]   for |y-H| >= y0
+  plug             : u_plug = (1/(2*mu*rho*g)) * (tau_w - tau0)^2           for |y-H| <  y0
+
+The key yield-term signature is a flat PLUG of uniform velocity within |y-H| < y0.
+
+Reads the raw restart binaries (restart_data/lustre_*.dat) written by the run
+(conservative variables, Fortran order, x fastest), extracts the x-averaged u(y)
+from the last save, computes the analytic profile, and reports the relative L2
+error and the measured plug half-width (region near the centerline within 1% of
+u_max) versus y0. Steady state is confirmed by comparing the last two saves.
+
+Run:  ./build/venv/bin/python3 examples/2D_bingham_poiseuille_nn/compare_analytic.py
+"""
+
+import glob
+import os
+
+import numpy as np
+
+# Must match case.py
+RHO = 1.0
+G_X = 0.1
+K = 5.0e-2  # n = 1 -> mu = K
+TAU0 = 4.0e-3
+L_X = 0.2
+L_Y = 0.2
+H = 0.5 * L_Y
+
+TAU_W = RHO * G_X * H
+Y0 = TAU0 / (RHO * G_X)  # analytic plug half-width
+
+# Grid (m, n in case.py -> m+1, n+1 cells)
+NX = 25
+NY = 64
+# Conservative variables per cell: alpha_rho(1), mom_x, mom_y, E, alpha(1)
+NVAR = 5
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+RESTART_DIR = os.path.join(HERE, "restart_data")
+
+
+def u_analytic(y):
+    """Closed-form steady Bingham plug-flow profile."""
+    s = np.abs(y - H)  # distance from centerline
+    tau_local = RHO * G_X * s  # |shear stress| = rho*g*s: 0 at center, tau_w at walls
+    u_plug = (TAU_W - TAU0) ** 2 / (2.0 * K * RHO * G_X)
+    # Sheared region where s >= y0, i.e. tau_local >= tau0
+    u_shear = ((TAU_W - TAU0) ** 2 - np.maximum(tau_local - TAU0, 0.0) ** 2) / (2.0 * K * RHO * G_X)
+    return np.where(s < Y0, u_plug, u_shear)
+
+
+def read_profile(path):
+    """Return x-averaged u(y) from one restart binary (Fortran order, x fastest)."""
+    raw = np.fromfile(path, dtype=np.float64).reshape((NVAR, NY, NX))
+    u = raw[1] / raw[0]  # mom_x / alpha_rho -> u[y, x]
+    return u.mean(axis=1)  # average over x (flow is x-invariant)
+
+
+def main():
+    dats = sorted(
+        glob.glob(os.path.join(RESTART_DIR, "lustre_[0-9]*.dat")),
+        key=lambda p: int(os.path.basename(p).split("_")[1].split(".")[0]),
+    )
+    if not dats:
+        raise SystemExit(f"No restart files in {RESTART_DIR}. " "Run: ./mfc.sh run examples/2D_bingham_poiseuille_nn/case.py -n 2")
+
+    ycb = np.fromfile(os.path.join(RESTART_DIR, "lustre_y_cb.dat"), dtype=np.float64)
+    y = 0.5 * (ycb[:-1] + ycb[1:])
+
+    u_last = read_profile(dats[-1])
+    u_ana = u_analytic(y)
+
+    if len(dats) >= 2:
+        u_prev = read_profile(dats[-2])
+        drift = np.sqrt(np.sum((u_last - u_prev) ** 2) / np.sum(u_last**2))
+        print(f"Steady-state drift between last two saves (rel L2): {drift:.3e}")
+
+    err = np.sqrt(np.sum((u_last - u_ana) ** 2) / np.sum(u_ana**2))
+    umax = u_last.max()
+    i_peak = int(np.argmin(np.abs(y - H)))
+    print(f"tau_w = rho*g*H = {TAU_W:.4e}   tau0 = {TAU0:.4e}   (need tau_w > tau0: {TAU_W > TAU0})")
+    print(f"u_max numeric  = {umax:.6e}  (analytic plug {u_ana.max():.6e})")
+    print(f"u at walls     = {u_last[0]:.3e}, {u_last[-1]:.3e}  (analytic 0)")
+    print(f"u at centerline= {u_last[i_peak]:.6e}  (analytic {u_ana[i_peak]:.6e})")
+    print(f"Relative L2 error vs analytic Bingham profile: {err:.4e}")
+
+    # Measure numerical plug half-width: contiguous band around centerline where
+    # u >= 0.99 * u_max. Half-width = max distance from centerline in that band.
+    plug_mask = u_last >= 0.99 * umax
+    plug_y = y[plug_mask]
+    if plug_y.size:
+        meas_half = max(abs(plug_y.max() - H), abs(plug_y.min() - H))
+    else:
+        meas_half = 0.0
+    print(f"Plug half-width: measured (>=99% u_max) = {meas_half:.4e}   analytic y0 = {Y0:.4e}   (= {Y0 / H:.2f} H)")
+    print(f"  ratio measured/analytic = {meas_half / Y0:.3f}")
+
+    # Near-wall momentum balance for the sheared region: mu*|du/dy| + tau0 = rho*g*(H-y).
+    dudy = np.gradient(u_last, y)
+    print("Sheared-region balance mu|du/dy|+tau0 vs rho*g*(H-y):")
+    for frac in (0.05, 0.10, 0.20):
+        i = int(np.argmin(np.abs(y - frac * L_Y)))
+        tau_num = K * np.abs(dudy[i]) + TAU0
+        tau_th = RHO * G_X * (H - y[i])
+        print(f"  y={y[i]:.4f}  num={tau_num:.4e}  theory={tau_th:.4e}  ratio={tau_num / tau_th:.3f}")
+
+
+if __name__ == "__main__":
+    main()