Add Herschel-Bulkley non-Newtonian viscosity model

README.md

@@ -148,6 +148,7 @@ MFC is a SPEChpc benchmark candidate, part of the JSC JUPITER Early Access Progr
 	* Euler-Lagrange particle tracking
 	* Quadrature-based moment methods (QBMM)
 * Viscous effects (high-order accurate representations)
+	* Newtonian and non-Newtonian rheology (Herschel-Bulkley: power-law, Bingham, and yield-stress fluids)
 * Hypoelastic and hyperelastic material models
 * Ideal and stiffened gas equations of state
 * Body forces

docs/documentation/case.md

@@ -1101,7 +1101,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/documentation/equations.md

@@ -387,6 +387,14 @@ and similarly for all other components. Cylindrical coordinate formulations incl
 
 Input parameters: `Re_inv` (shear and volume Reynolds numbers per fluid).
 
+**Non-Newtonian viscosity (`fluid_pp(i)%%non_newtonian = .true.`):**
+
+Per-fluid Herschel-Bulkley rheology with Papanastasiou regularization (\cite Papanastasiou87) replaces the constant viscosity with a shear-rate-dependent effective viscosity, recomputed from the local strain-rate tensor at every time step:
+
+\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\,\mathbf{D}:\mathbf{D}}\f]
+
+This covers power-law shear-thinning/thickening (\f$\tau_0 = 0\f$), Bingham plastic (\f$n = 1,\ \tau_0 > 0\f$), and general Herschel-Bulkley yield-stress fluids. The mixture rule above applies with \f$1/\text{Re}_j = \mu_{{\rm eff},j}\f$ evaluated at the local shear rate; Newtonian and non-Newtonian fluids can be mixed. Supported with the HLL and HLLC Riemann solvers, including immersed boundaries. See @ref sec-non-newtonian "the case documentation" for parameters, constraints, and validated example cases.
+
 ---
 
 ## 5. Cylindrical Coordinates (`cyl_coord = .true.`) (\cite Wilfong26 Sec. 2.3)

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,
+        }
+    )
+)