MFlowCode · jasontruong2707 · Mar 26, 2026 · Mar 27, 2026 · Mar 27, 2026 · Mar 27, 2026
diff --git a/docs/documentation/case.md b/docs/documentation/case.md
@@ -386,6 +386,14 @@ Additional details on this specification can be found in [NACA airfoil](https://
 | `qvp`  ** | Real   | Stiffened-gas parameter $q'$ of fluid.         |
 | `sigma`   | Real   | Surface tension coefficient                    |
 | `G`       | Real   | Shear modulus of solid.                        |
+| `non_newtonian` | Logical | Enable non-Newtonian (Herschel-Bulkley) viscosity. See @ref sec-non-newtonian. |
+| `tau0`    | Real   | Yield stress \f$\tau_0\f$ (Herschel-Bulkley).  |
+| `K`       | Real   | Consistency index \f$K\f$ (Herschel-Bulkley).  |
+| `nn`      | Real   | Flow behavior index \f$n\f$ (Herschel-Bulkley).|
+| `hb_m`    | Real   | Papanastasiou regularization parameter \f$m\f$. |
+| `mu_min`  | Real   | Minimum viscosity clamp.                       |
+| `mu_max`  | Real   | Maximum viscosity clamp.                       |
+| `mu_bulk` | Real   | Bulk viscosity (non-Newtonian).                |
 
 Fluid material's parameters. All parameters except for sigma should be prepended with `fluid_pp(i)` where $i$ is the fluid index.
 
@@ -403,6 +411,38 @@ The parameters define material's property of compressible fluids that are used i
 When these parameters are undefined, fluids are treated as inviscid.
 Details of implementation of viscosity in MFC can be found in \cite Coralic15.
 
+#### Non-Newtonian Viscosity (Herschel-Bulkley) {#sec-non-newtonian}
+
+MFC supports non-Newtonian viscosity via the Herschel-Bulkley model with Papanastasiou regularization.
+Enable it per-fluid by setting ``fluid_pp(i)%%non_newtonian = 'T'`` (requires ``viscous = 'T'``).
+
+The effective viscosity is:
+
+\f[ \mu(\dot\gamma) = \frac{\tau_0}{\dot\gamma}\bigl(1 - e^{-m\,\dot\gamma}\bigr) + K\,\dot\gamma^{\,n-1} \f]
+
+where \f$\dot\gamma = \sqrt{2\,D_{ij}\,D_{ij}}\f$ is the shear rate computed from the strain-rate tensor.
+
+| Parameter | Type | Description |
+| ---: | :---: | :--- |
+| ``non_newtonian`` | Logical | Enable non-Newtonian viscosity for this fluid |
+| ``tau0``   | Real | Yield stress \f$\tau_0\f$. Set to 0 for a power-law fluid |
+| ``K``      | Real | Consistency index \f$K\f$ (must be > 0) |
+| ``nn``     | Real | Flow behavior index \f$n\f$: shear-thinning (\f$n<1\f$), Newtonian (\f$n=1\f$), shear-thickening (\f$n>1\f$) |
+| ``hb_m``   | Real | Papanastasiou regularization parameter \f$m\f$. Higher values approximate the true yield stress more closely (typical: 1000--10000) |
+| ``mu_min`` | Real | Minimum viscosity clamp (must be >= 0) |
+| ``mu_max`` | Real | Maximum viscosity clamp (must be > ``mu_min``) |
+| ``mu_bulk``| Real | Bulk viscosity for non-Newtonian fluids (optional, default 0 = inviscid bulk) |
+
+All parameters are prepended with ``fluid_pp(i)%%``.
+
+**Special cases:**
+- \f$\tau_0 = 0\f$: reduces to power-law model \f$\mu = K\,\dot\gamma^{\,n-1}\f$.
+- \f$n = 1,\, \tau_0 = 0\f$: reduces to Newtonian with \f$\mu = K\f$.
+- \f$n = 1,\, \tau_0 > 0\f$: reduces to Bingham plastic.
+
+> **Note:** When ``non_newtonian = 'T'``, the ``Re(1)``/``Re(2)`` parameters are ignored for that fluid.
+> The viscosity is computed from the Herschel-Bulkley model instead.
+
 - `fluid_pp(i)%%cv`, `fluid_pp(i)%%qv`, and `fluid_pp(i)%%qvp` define $c_v$, $q$, and $q'$ as parameters of $i$-th fluid that are used in stiffened gas equation of state.
 
 - `fluid_pp(i)%%G` is required for `hypoelasticity`.

diff --git a/docs/module_categories.json b/docs/module_categories.json
@@ -27,6 +27,8 @@
             "m_acoustic_src",
             "m_body_forces",
             "m_pressure_relaxation",
+            "m_hb_function",
+            "m_re_visc",
             "m_collisions"
         ]
     },

@@ -0,0 +1,117 @@
+#!/usr/bin/env python3
+"""
+2D Lid-Driven Cavity Flow with Herschel-Bulkley Non-Newtonian Fluid
+
+Re = 5000, n = 1.5 (shear-thickening) with mesh stretching near walls.
+
+HB model: mu = (tau0/gdot)*(1 - exp(-m*gdot)) + K * gdot^(n-1)
+Re_gen = rho * U^(2-n) * L^n / K = 1 * 1^0.5 * 1^1.5 / 0.0002 = 5000
+
+Mesh stretching: cosh-based clustering near all 4 walls (x_a, x_b, y_a, y_b).
-Re = 5000, n = 1.5 (shear-thickening) with mesh stretching near walls.
-
-HB model: mu = (tau0/gdot)*(1 - exp(-m*gdot)) + K * gdot^(n-1)
-Re_gen = rho * U^(2-n) * L^n / K = 1 * 1^0.5 * 1^1.5 / 0.0002 = 5000
-
-Mesh stretching: cosh-based clustering near all 4 walls (x_a, x_b, y_a, y_b).
+Re = 5000, n = 1.5 (shear-thickening) on a uniform grid.
+
+HB model: mu = (tau0/gdot)*(1 - exp(-m*gdot)) + K * gdot^(n-1)
+Re_gen = rho * U^(2-n) * L^n / K = 1 * 1^0.5 * 1^1.5 / 0.0002 = 5000
+
+(No mesh stretching parameters are enabled in this case file.)
-Re = 5000, n = 1.5 (shear-thickening) with mesh stretching near walls.
-
-HB model: mu = (tau0/gdot)*(1 - exp(-m*gdot)) + K * gdot^(n-1)
-Re_gen = rho * U^(2-n) * L^n / K = 1 * 1^0.5 * 1^1.5 / 0.0002 = 5000
-
-Mesh stretching: cosh-based clustering near all 4 walls (x_a, x_b, y_a, y_b).
+Re = 5000, n = 1.5 (shear-thickening) on a uniform grid.
+
+HB model: mu = (tau0/gdot)*(1 - exp(-m*gdot)) + K * gdot^(n-1)
+Re_gen = rho * U^(2-n) * L^n / K = 1 * 1^0.5 * 1^1.5 / 0.0002 = 5000
+
+(No mesh stretching parameters are enabled in this case file.)
+"""
+
+import json
+
+eps = 1e-6
+
+# HB model parameters
+tau0 = 0.0  # Yield stress (set to 0 for power-law fluid)
+K = 0.0002  # Consistency index (Re=5000: K = 1/5000)
+nn = 1.5  # Flow behavior index (shear-thickening)
+mu_min = 0.00002  # K * gdot_min^(n-1) = 0.0002 * (0.01)^0.5
+mu_max = 0.0632  # K * gdot_max^(n-1) = 0.0002 * (1e5)^0.5
+hb_m = 1000.0  # Papanastasiou regularization parameter
+mu_bulk = 0.0
+
+lid_velocity = 1.0  # Lid velocity (m/s)
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0.0,
+            "x_domain%end": 1.0,
+            "y_domain%beg": 0.0,
+            "y_domain%end": 1.0,
+            "m": 255,
+            "n": 255,
+            "p": 0,
+            "cfl_adap_dt": "T",
+            "cfl_target": 0.5,
+            "n_start": 0,
+            "t_stop": 50.0,
+            "t_save": 25.0,
+            # Simulation Algorithm Parameters
+            "num_patches": 1,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 2,
+            "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,
+            "bc_x%beg": -16,
+            "bc_x%end": -16,
+            "bc_y%beg": -16,
+            "bc_y%end": -16,
+            "bc_y%ve1": lid_velocity,
+            "viscous": "T",
+            # Formatted Database Files Structure Parameters
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "omega_wrt(3)": "T",
+            "fd_order": 4,
+            "parallel_io": "T",
+            # Patch 1: Base
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%x_centroid": 0.5,
+            "patch_icpp(1)%y_centroid": 0.5,
+            "patch_icpp(1)%length_x": 1.0,
+            "patch_icpp(1)%length_y": 1.0,
+            "patch_icpp(1)%vel(1)": 0,
+            "patch_icpp(1)%vel(2)": 0.0,
+            "patch_icpp(1)%pres": 1e3,
+            "patch_icpp(1)%alpha_rho(1)": 0.5,
+            "patch_icpp(1)%alpha(1)": 0.5,
+            "patch_icpp(1)%alpha_rho(2)": 0.5,
+            "patch_icpp(1)%alpha(2)": 0.5,
+            # Fluids Physical Parameters
+            # Fluid 1:
+            "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)%tau0": tau0,
+            "fluid_pp(1)%K": K,
+            "fluid_pp(1)%nn": nn,
+            "fluid_pp(1)%mu_max": mu_max,
+            "fluid_pp(1)%mu_min": mu_min,
+            "fluid_pp(1)%mu_bulk": mu_bulk,
+            "fluid_pp(1)%hb_m": hb_m,
+            # Fluid 2:
+            "fluid_pp(2)%gamma": 1.0 / (1.4 - 1.0),
+            "fluid_pp(2)%pi_inf": 0.0,
+            "fluid_pp(2)%Re(1)": 1.0 / K,
+            "fluid_pp(2)%non_newtonian": "T",
+            "fluid_pp(2)%tau0": tau0,
+            "fluid_pp(2)%K": K,
+            "fluid_pp(2)%nn": nn,
+            "fluid_pp(2)%mu_max": mu_max,
+            "fluid_pp(2)%mu_min": mu_min,
+            "fluid_pp(2)%mu_bulk": mu_bulk,
+            "fluid_pp(2)%hb_m": hb_m,
+        }
+    )
+)
@@ -0,0 +1,159 @@
+#!/usr/bin/env python3
+"""
+2D Poiseuille Flow with Herschel-Bulkley Non-Newtonian Fluid
+
+Pressure-driven channel flow between two no-slip walls, driven by a constant
+body force in the streamwise (x) direction. Periodic BCs in x, no-slip in y.
+
+HB model: mu = (tau0/gdot)*(1 - exp(-m*gdot)) + K * gdot^(n-1)
+  - tau0: yield stress
+  - K:    consistency index
+  - n:    flow behavior index (< 1 shear-thinning, > 1 shear-thickening)
+  - m:    Papanastasiou regularization parameter
+
+For Newtonian Poiseuille validation, set tau0=0, nn=1, K=mu.
+The analytical solution is: u(y) = (G/(2*mu)) * y * (H - y)
+where G = rho * g_x is the effective pressure gradient.
+"""
+
+import json
+import math
+
+# -- Channel geometry (square domain) --
+L = 1.0  # Channel length (streamwise, x)
+H = 1.0  # Channel height (wall-normal, y)
+
+# -- Grid resolution --
+Nx = 24  # Cells in x (streamwise, minimal — periodic)
+Ny = 81  # Cells in y (wall-normal)
+
+# -- Fluid properties --
+rho = 1.0  # Density
+p0 = 1e5  # Reference pressure (high for low Mach)
+gamma = 1.4  # Ratio of specific heats
+
+# Sound speed and CFL
+c = math.sqrt(gamma * p0 / rho)
+dx = L / (Nx + 1)
+cfl = 0.3
+dt = cfl * dx / c
+
+# -- Body force (pressure gradient substitute) --
+# G = rho * g_x acts as dp/dx driving force
+g_x = 0.5
+
+# -- HB non-Newtonian model parameters --
+tau0 = 0.0  # Yield stress (set 0 for power-law)
+K = 0.1  # Consistency index
+nn = 2.0  # Flow behavior index (< 1 = shear-thinning)
+hb_m = 1000.0  # Papanastasiou regularization parameter
+mu_min = 1e-4  # Minimum viscosity bound
+mu_max = 10.0  # Maximum viscosity bound
+mu_bulk = 0.0  # Bulk viscosity
+
+# Reference Re based on consistency index (used as baseline)
+Re_ref = 1.0 / K  # = 100
+
+# -- Time control --
+t_end = 10.0  # End time (allow flow to reach steady state)
+t_save = 5.0  # Save interval
+
+eps = 1e-6
+
+# Configuring case dictionary
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": 0.0,
+            "x_domain%end": L,
+            "y_domain%beg": 0.0,
+            "y_domain%end": H,
+            "m": Nx,
+            "n": Ny,
+            "p": 0,
+            "cfl_adap_dt": "T",
+            "cfl_target": cfl,
+            "n_start": 0,
+            "t_stop": t_end,
+            "t_save": t_save,
+            # Simulation Algorithm Parameters
+            "num_patches": 1,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 2,
+            "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,
+            # Boundary Conditions
+            # Periodic in x (streamwise), no-slip walls in y
+            "bc_x%beg": -1,
+            "bc_x%end": -1,
+            "bc_y%beg": -16,
+            "bc_y%end": -16,
+            # Viscous
+            "viscous": "T",
+            # Body Force (drives the flow like a pressure gradient)
+            "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",
+            "omega_wrt(3)": "T",
+            "fd_order": 4,
+            "parallel_io": "T",
+            # Patch 1: Entire channel domain (initially at rest)
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%x_centroid": L / 2.0,
+            "patch_icpp(1)%y_centroid": H / 2.0,
+            "patch_icpp(1)%length_x": L,
+            "patch_icpp(1)%length_y": H,
+            "patch_icpp(1)%vel(1)": 0.0,
+            "patch_icpp(1)%vel(2)": 0.0,
+            "patch_icpp(1)%pres": p0,
+            "patch_icpp(1)%alpha_rho(1)": rho * 0.5,
+            "patch_icpp(1)%alpha(1)": 0.5,
+            "patch_icpp(1)%alpha_rho(2)": rho * 0.5,
+            "patch_icpp(1)%alpha(2)": 0.5,
+            # Fluid 1: HB non-Newtonian fluid
+            "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0),
+            "fluid_pp(1)%pi_inf": 0.0,
+            "fluid_pp(1)%Re(1)": Re_ref,
+            "fluid_pp(1)%non_newtonian": "T",
+            "fluid_pp(1)%tau0": tau0,
+            "fluid_pp(1)%K": K,
+            "fluid_pp(1)%nn": nn,
+            "fluid_pp(1)%hb_m": hb_m,
+            "fluid_pp(1)%mu_min": mu_min,
+            "fluid_pp(1)%mu_max": mu_max,
+            "fluid_pp(1)%mu_bulk": mu_bulk,
+            # Fluid 2: same properties (single-phase effectively)
+            "fluid_pp(2)%gamma": 1.0 / (gamma - 1.0),
+            "fluid_pp(2)%pi_inf": 0.0,
+            "fluid_pp(2)%Re(1)": Re_ref,
+            "fluid_pp(2)%non_newtonian": "T",
+            "fluid_pp(2)%tau0": tau0,
+            "fluid_pp(2)%K": K,
+            "fluid_pp(2)%nn": nn,
+            "fluid_pp(2)%hb_m": hb_m,
+            "fluid_pp(2)%mu_min": mu_min,
+            "fluid_pp(2)%mu_max": mu_max,
+            "fluid_pp(2)%mu_bulk": mu_bulk,
+        }
+    )
+)