Skip to content

Commit e01e168

Browse files
sbryngelsonsbryngelson
committed
Add Herschel-Bulkley non-Newtonian viscosity model
1 parent 52e8af1 commit e01e168

39 files changed

+3318
-775
lines changed

‎docs/documentation/case.md‎

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -374,6 +374,14 @@ Additional details on this specification can be found in [NACA airfoil](https://
374374
| `qvp` ** | Real | Stiffened-gas parameter $q'$ of fluid. |
375375
| `sigma` | Real | Surface tension coefficient |
376376
| `G` | Real | Shear modulus of solid. |
377+
| `non_newtonian` | Logical | Enable non-Newtonian (Herschel-Bulkley) viscosity. See @ref sec-non-newtonian. |
378+
| `tau0` | Real | Yield stress \f$\tau_0\f$ (Herschel-Bulkley). |
379+
| `K` | Real | Consistency index \f$K\f$ (Herschel-Bulkley). |
380+
| `nn` | Real | Flow behavior index \f$n\f$ (Herschel-Bulkley).|
381+
| `hb_m` | Real | Papanastasiou regularization parameter \f$m\f$. |
382+
| `mu_min` | Real | Minimum viscosity clamp. |
383+
| `mu_max` | Real | Maximum viscosity clamp. |
384+
| `mu_bulk` | Real | Bulk viscosity (non-Newtonian). |
377385

378386
Fluid material's parameters. All parameters except for sigma should be prepended with `fluid_pp(i)` where $i$ is the fluid index.
379387

@@ -391,6 +399,38 @@ The parameters define material's property of compressible fluids that are used i
391399
When these parameters are undefined, fluids are treated as inviscid.
392400
Details of implementation of viscosity in MFC can be found in \cite Coralic15.
393401

402+
#### Non-Newtonian Viscosity (Herschel-Bulkley) {#sec-non-newtonian}
403+
404+
MFC supports non-Newtonian viscosity via the Herschel-Bulkley model with Papanastasiou regularization.
405+
Enable it per-fluid by setting ``fluid_pp(i)%%non_newtonian = 'T'`` (requires ``viscous = 'T'``).
406+
407+
The effective viscosity is:
408+
409+
\f[ \mu(\dot\gamma) = \frac{\tau_0}{\dot\gamma}\bigl(1 - e^{-m\,\dot\gamma}\bigr) + K\,\dot\gamma^{\,n-1} \f]
410+
411+
where \f$\dot\gamma = \sqrt{2\,D_{ij}\,D_{ij}}\f$ is the shear rate computed from the strain-rate tensor.
412+
413+
| Parameter | Type | Description |
414+
| ---: | :---: | :--- |
415+
| ``non_newtonian`` | Logical | Enable non-Newtonian viscosity for this fluid |
416+
| ``tau0`` | Real | Yield stress \f$\tau_0\f$. Set to 0 for a power-law fluid |
417+
| ``K`` | Real | Consistency index \f$K\f$ (must be > 0) |
418+
| ``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$) |
419+
| ``hb_m`` | Real | Papanastasiou regularization parameter \f$m\f$. Higher values approximate the true yield stress more closely (typical: 1000--10000) |
420+
| ``mu_min`` | Real | Minimum viscosity clamp (must be >= 0) |
421+
| ``mu_max`` | Real | Maximum viscosity clamp (must be > ``mu_min``) |
422+
| ``mu_bulk``| Real | Bulk viscosity for non-Newtonian fluids (optional, default 0 = inviscid bulk) |
423+
424+
All parameters are prepended with ``fluid_pp(i)%%``.
425+
426+
**Special cases:**
427+
- \f$\tau_0 = 0\f$: reduces to power-law model \f$\mu = K\,\dot\gamma^{\,n-1}\f$.
428+
- \f$n = 1,\, \tau_0 = 0\f$: reduces to Newtonian with \f$\mu = K\f$.
429+
- \f$n = 1,\, \tau_0 > 0\f$: reduces to Bingham plastic.
430+
431+
> **Note:** When ``non_newtonian = 'T'``, the ``Re(1)``/``Re(2)`` parameters are ignored for that fluid.
432+
> The viscosity is computed from the Herschel-Bulkley model instead.
433+
394434
- `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.
395435

396436
- `fluid_pp(i)%%G` is required for `hypoelasticity`.

‎docs/module_categories.json‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,9 @@
2626
"m_chemistry",
2727
"m_acoustic_src",
2828
"m_body_forces",
29-
"m_pressure_relaxation"
29+
"m_pressure_relaxation",
30+
"m_hb_function",
31+
"m_re_visc"
3032
]
3133
},
3234
{

‎examples/2D_lid_driven_cavity_nn/Lid_driven_cavity_Re_100_n_0.5.png‎

198 KB
Loading

‎examples/2D_lid_driven_cavity_nn/Lid_driven_cavity_Re_100_n_1.5.png‎

191 KB
Loading

‎examples/2D_lid_driven_cavity_nn/case.py‎

Lines changed: 117 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,117 @@
1+
#!/usr/bin/env python3
2+
"""
3+
2D Lid-Driven Cavity Flow with Herschel-Bulkley Non-Newtonian Fluid
4+
5+
Re = 5000, n = 1.5 (shear-thickening) with mesh stretching near walls.
6+
7+
HB model: mu = (tau0/gdot)*(1 - exp(-m*gdot)) + K * gdot^(n-1)
8+
Re_gen = rho * U^(2-n) * L^n / K = 1 * 1^0.5 * 1^1.5 / 0.0002 = 5000
9+
10+
Mesh stretching: cosh-based clustering near all 4 walls (x_a, x_b, y_a, y_b).
11+
"""
12+
13+
import json
14+
15+
eps = 1e-6
16+
17+
# HB model parameters
18+
tau0 = 0.0 # Yield stress (set to 0 for power-law fluid)
19+
K = 0.0002 # Consistency index (Re=5000: K = 1/5000)
20+
nn = 1.5 # Flow behavior index (shear-thickening)
21+
mu_min = 0.00002 # K * gdot_min^(n-1) = 0.0002 * (0.01)^0.5
22+
mu_max = 0.0632 # K * gdot_max^(n-1) = 0.0002 * (1e5)^0.5
23+
hb_m = 1000.0 # Papanastasiou regularization parameter
24+
mu_bulk = 0.0
25+
26+
lid_velocity = 1.0 # Lid velocity (m/s)
27+
28+
# Configuring case dictionary
29+
print(
30+
json.dumps(
31+
{
32+
# Logistics
33+
"run_time_info": "T",
34+
# Computational Domain Parameters
35+
"x_domain%beg": 0.0,
36+
"x_domain%end": 1.0,
37+
"y_domain%beg": 0.0,
38+
"y_domain%end": 1.0,
39+
"m": 255,
40+
"n": 255,
41+
"p": 0,
42+
"cfl_adap_dt": "T",
43+
"cfl_target": 0.5,
44+
"n_start": 0,
45+
"t_stop": 50.0,
46+
"t_save": 25.0,
47+
# Simulation Algorithm Parameters
48+
"num_patches": 1,
49+
"model_eqns": 2,
50+
"alt_soundspeed": "F",
51+
"num_fluids": 2,
52+
"mpp_lim": "F",
53+
"mixture_err": "T",
54+
"time_stepper": 3,
55+
"weno_order": 5,
56+
"weno_eps": 1e-16,
57+
"mapped_weno": "T",
58+
"weno_Re_flux": "T",
59+
"mp_weno": "T",
60+
"weno_avg": "T",
61+
"riemann_solver": 2,
62+
"wave_speeds": 1,
63+
"avg_state": 2,
64+
"bc_x%beg": -16,
65+
"bc_x%end": -16,
66+
"bc_y%beg": -16,
67+
"bc_y%end": -16,
68+
"bc_y%ve1": lid_velocity,
69+
"viscous": "T",
70+
# Formatted Database Files Structure Parameters
71+
"format": 1,
72+
"precision": 2,
73+
"prim_vars_wrt": "T",
74+
"omega_wrt(3)": "T",
75+
"fd_order": 4,
76+
"parallel_io": "T",
77+
# Patch 1: Base
78+
"patch_icpp(1)%geometry": 3,
79+
"patch_icpp(1)%x_centroid": 0.5,
80+
"patch_icpp(1)%y_centroid": 0.5,
81+
"patch_icpp(1)%length_x": 1.0,
82+
"patch_icpp(1)%length_y": 1.0,
83+
"patch_icpp(1)%vel(1)": 0,
84+
"patch_icpp(1)%vel(2)": 0.0,
85+
"patch_icpp(1)%pres": 1e3,
86+
"patch_icpp(1)%alpha_rho(1)": 0.5,
87+
"patch_icpp(1)%alpha(1)": 0.5,
88+
"patch_icpp(1)%alpha_rho(2)": 0.5,
89+
"patch_icpp(1)%alpha(2)": 0.5,
90+
# Fluids Physical Parameters
91+
# Fluid 1:
92+
"fluid_pp(1)%gamma": 1.0 / (1.4 - 1.0),
93+
"fluid_pp(1)%pi_inf": 0.0,
94+
"fluid_pp(1)%Re(1)": 1.0 / K,
95+
"fluid_pp(1)%non_newtonian": "T",
96+
"fluid_pp(1)%tau0": tau0,
97+
"fluid_pp(1)%K": K,
98+
"fluid_pp(1)%nn": nn,
99+
"fluid_pp(1)%mu_max": mu_max,
100+
"fluid_pp(1)%mu_min": mu_min,
101+
"fluid_pp(1)%mu_bulk": mu_bulk,
102+
"fluid_pp(1)%hb_m": hb_m,
103+
# Fluid 2:
104+
"fluid_pp(2)%gamma": 1.0 / (1.4 - 1.0),
105+
"fluid_pp(2)%pi_inf": 0.0,
106+
"fluid_pp(2)%Re(1)": 1.0 / K,
107+
"fluid_pp(2)%non_newtonian": "T",
108+
"fluid_pp(2)%tau0": tau0,
109+
"fluid_pp(2)%K": K,
110+
"fluid_pp(2)%nn": nn,
111+
"fluid_pp(2)%mu_max": mu_max,
112+
"fluid_pp(2)%mu_min": mu_min,
113+
"fluid_pp(2)%mu_bulk": mu_bulk,
114+
"fluid_pp(2)%hb_m": hb_m,
115+
}
116+
)
117+
)

‎examples/2D_poiseuille_nn/case.py‎

Lines changed: 159 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,159 @@
1+
#!/usr/bin/env python3
2+
"""
3+
2D Poiseuille Flow with Herschel-Bulkley Non-Newtonian Fluid
4+
5+
Pressure-driven channel flow between two no-slip walls, driven by a constant
6+
body force in the streamwise (x) direction. Periodic BCs in x, no-slip in y.
7+
8+
HB model: mu = (tau0/gdot)*(1 - exp(-m*gdot)) + K * gdot^(n-1)
9+
- tau0: yield stress
10+
- K: consistency index
11+
- n: flow behavior index (< 1 shear-thinning, > 1 shear-thickening)
12+
- m: Papanastasiou regularization parameter
13+
14+
For Newtonian Poiseuille validation, set tau0=0, nn=1, K=mu.
15+
The analytical solution is: u(y) = (G/(2*mu)) * y * (H - y)
16+
where G = rho * g_x is the effective pressure gradient.
17+
"""
18+
19+
import json
20+
import math
21+
22+
# -- Channel geometry (square domain) --
23+
L = 1.0 # Channel length (streamwise, x)
24+
H = 1.0 # Channel height (wall-normal, y)
25+
26+
# -- Grid resolution --
27+
Nx = 24 # Cells in x (streamwise, minimal — periodic)
28+
Ny = 81 # Cells in y (wall-normal)
29+
30+
# -- Fluid properties --
31+
rho = 1.0 # Density
32+
p0 = 1e5 # Reference pressure (high for low Mach)
33+
gamma = 1.4 # Ratio of specific heats
34+
35+
# Sound speed and CFL
36+
c = math.sqrt(gamma * p0 / rho)
37+
dx = L / (Nx + 1)
38+
cfl = 0.3
39+
dt = cfl * dx / c
40+
41+
# -- Body force (pressure gradient substitute) --
42+
# G = rho * g_x acts as dp/dx driving force
43+
g_x = 0.5
44+
45+
# -- HB non-Newtonian model parameters --
46+
tau0 = 0.0 # Yield stress (set 0 for power-law)
47+
K = 0.1 # Consistency index
48+
nn = 2.0 # Flow behavior index (< 1 = shear-thinning)
49+
hb_m = 1000.0 # Papanastasiou regularization parameter
50+
mu_min = 1e-4 # Minimum viscosity bound
51+
mu_max = 10.0 # Maximum viscosity bound
52+
mu_bulk = 0.0 # Bulk viscosity
53+
54+
# Reference Re based on consistency index (used as baseline)
55+
Re_ref = 1.0 / K # = 100
56+
57+
# -- Time control --
58+
t_end = 10.0 # End time (allow flow to reach steady state)
59+
t_save = 5.0 # Save interval
60+
61+
eps = 1e-6
62+
63+
# Configuring case dictionary
64+
print(
65+
json.dumps(
66+
{
67+
# Logistics
68+
"run_time_info": "T",
69+
# Computational Domain Parameters
70+
"x_domain%beg": 0.0,
71+
"x_domain%end": L,
72+
"y_domain%beg": 0.0,
73+
"y_domain%end": H,
74+
"m": Nx,
75+
"n": Ny,
76+
"p": 0,
77+
"cfl_adap_dt": "T",
78+
"cfl_target": cfl,
79+
"n_start": 0,
80+
"t_stop": t_end,
81+
"t_save": t_save,
82+
# Simulation Algorithm Parameters
83+
"num_patches": 1,
84+
"model_eqns": 2,
85+
"alt_soundspeed": "F",
86+
"num_fluids": 2,
87+
"mpp_lim": "F",
88+
"mixture_err": "T",
89+
"time_stepper": 3,
90+
"weno_order": 5,
91+
"weno_eps": 1e-16,
92+
"mapped_weno": "T",
93+
"weno_Re_flux": "T",
94+
"mp_weno": "T",
95+
"weno_avg": "T",
96+
"riemann_solver": 2,
97+
"wave_speeds": 1,
98+
"avg_state": 2,
99+
# Boundary Conditions
100+
# Periodic in x (streamwise), no-slip walls in y
101+
"bc_x%beg": -1,
102+
"bc_x%end": -1,
103+
"bc_y%beg": -16,
104+
"bc_y%end": -16,
105+
# Viscous
106+
"viscous": "T",
107+
# Body Force (drives the flow like a pressure gradient)
108+
"bf_x": "T",
109+
"g_x": g_x,
110+
"k_x": 0.0,
111+
"w_x": 0.0,
112+
"p_x": 0.0,
113+
# Formatted Database Files Structure Parameters
114+
"format": 1,
115+
"precision": 2,
116+
"prim_vars_wrt": "T",
117+
"omega_wrt(3)": "T",
118+
"fd_order": 4,
119+
"parallel_io": "T",
120+
# Patch 1: Entire channel domain (initially at rest)
121+
"patch_icpp(1)%geometry": 3,
122+
"patch_icpp(1)%x_centroid": L / 2.0,
123+
"patch_icpp(1)%y_centroid": H / 2.0,
124+
"patch_icpp(1)%length_x": L,
125+
"patch_icpp(1)%length_y": H,
126+
"patch_icpp(1)%vel(1)": 0.0,
127+
"patch_icpp(1)%vel(2)": 0.0,
128+
"patch_icpp(1)%pres": p0,
129+
"patch_icpp(1)%alpha_rho(1)": rho * 0.5,
130+
"patch_icpp(1)%alpha(1)": 0.5,
131+
"patch_icpp(1)%alpha_rho(2)": rho * 0.5,
132+
"patch_icpp(1)%alpha(2)": 0.5,
133+
# Fluid 1: HB non-Newtonian fluid
134+
"fluid_pp(1)%gamma": 1.0 / (gamma - 1.0),
135+
"fluid_pp(1)%pi_inf": 0.0,
136+
"fluid_pp(1)%Re(1)": Re_ref,
137+
"fluid_pp(1)%non_newtonian": "T",
138+
"fluid_pp(1)%tau0": tau0,
139+
"fluid_pp(1)%K": K,
140+
"fluid_pp(1)%nn": nn,
141+
"fluid_pp(1)%hb_m": hb_m,
142+
"fluid_pp(1)%mu_min": mu_min,
143+
"fluid_pp(1)%mu_max": mu_max,
144+
"fluid_pp(1)%mu_bulk": mu_bulk,
145+
# Fluid 2: same properties (single-phase effectively)
146+
"fluid_pp(2)%gamma": 1.0 / (gamma - 1.0),
147+
"fluid_pp(2)%pi_inf": 0.0,
148+
"fluid_pp(2)%Re(1)": Re_ref,
149+
"fluid_pp(2)%non_newtonian": "T",
150+
"fluid_pp(2)%tau0": tau0,
151+
"fluid_pp(2)%K": K,
152+
"fluid_pp(2)%nn": nn,
153+
"fluid_pp(2)%hb_m": hb_m,
154+
"fluid_pp(2)%mu_min": mu_min,
155+
"fluid_pp(2)%mu_max": mu_max,
156+
"fluid_pp(2)%mu_bulk": mu_bulk,
157+
}
158+
)
159+
)

‎examples/2D_poiseuille_nn/velocity_profile.png‎

393 KB
Loading

0 commit comments

Comments
 (0)