Generalized-Newtonian rheology

Many process fluids are non-Newtonian: their viscosity is not a constant but a function of the local deformation rate. Shear-thinning broths (filamentous fermentation, xanthan, high-cell-density cultures), yield-stress pastes, and polymer solutions all fall outside the Newtonian assumption that underlies the base solver. Nassu models this class through the generalized-Newtonian framework, in which the fluid keeps the Newtonian tensorial constitutive form but the viscosity becomes an apparent viscosity \(\eta(\dot\gamma)\) evaluated from the local shear rate.

The construction reuses machinery the solver already has: the shear-rate scalar is exactly the strain-rate magnitude \(|S|\) that the Smagorinsky model already computes, and the apparent viscosity enters only through the relaxation frequency \(\omega\), exactly where the LES eddy viscosity enters. No new streaming, collision, or communication path is introduced.

Constitutive framework

A generalized-Newtonian fluid keeps the Newtonian tensorial form of the constitutive law but lets the proportionality depend on the local deformation rate:

(1)\[ \tau_{\alpha\beta} = 2\,\eta(\dot\gamma)\, S_{\alpha\beta}, \qquad S_{\alpha\beta} = \tfrac{1}{2}\!\left(\partial_\alpha u_\beta + \partial_\beta u_\alpha\right) \]

where \(\eta(\dot\gamma)\) is the apparent dynamic viscosity and \(\dot\gamma\) is the shear-rate magnitude, the scalar invariant of the rate-of-strain tensor:

(2)\[ \dot\gamma \;=\; \sqrt{2\,S_{\alpha\beta}S_{\alpha\beta}} \;\equiv\; |S| \]

This is identical to the magnitude already fed to the Smagorinsky model (see the LES chapter, Eq. smagorinsky), written in the polymer-rheology convention Bird et al.[1] where the shear-rate tensor is \(\dot\gamma_{\alpha\beta}=2S_{\alpha\beta}\). In simple shear \(u_x=\dot\gamma\,y\) only \(S_{xy}=S_{yx}=\dot\gamma/2\) are nonzero, so \(|S|=\dot\gamma\): the scalar coincides with the engineering shear rate, which is why viscometric data measured in simple shear insert directly into \(\eta(\dot\gamma)\).

Because \(S_{\alpha\beta}\) is recovered algebraically from the local non-equilibrium stress \(\Pi^{\mathrm{neq}}_{\alpha\beta}\) (LBM chapter, Eq. strain), \(\dot\gamma\) is available at every node with no velocity gradients and no neighbour communication - the same locality that makes LBM-LES cheap is reused to evaluate a shear-dependent molecular viscosity.

Viscosity models

All formulas give the dynamic apparent viscosity \(\eta\); the mapping to the lattice relaxation follows in the next section. Physical \(\eta\) is converted to lattice units with the standard viscosity factor before that step.

Power-law (Ostwald-de Waele)

(3)\[ \eta(\dot\gamma) = K\,\dot\gamma^{\,n-1} \]

with consistency index \(K\) and flow-behaviour index \(n\): \(n<1\) shear-thinning, \(n>1\) shear-thickening, \(n=1\) Newtonian (\(\eta=K\)). Both tails are unphysical (\(\eta\to\infty\) or \(\eta\to0\)), so the power-law is a mid-shear fit that must be clipped (see stability below) Chhabra and Richardson[2].

Carreau-Yasuda

(4)\[ \eta(\dot\gamma) = \eta_\infty + (\eta_0 - \eta_\infty)\, \Big[\,1 + (\lambda\dot\gamma)^{a}\,\Big]^{\frac{n-1}{a}} \]

with zero-shear plateau \(\eta_0\), infinite-shear plateau \(\eta_\infty\), relaxation time \(\lambda\), thinning index \(n\), and Yasuda transition-width \(a\) (the original Carreau model is \(a=2\)) Bird et al.[1]. The model is inherently bounded, \(\eta_\infty\le\eta\le\eta_0\), so it needs little or no clipping - the numerically friendliest choice.

Herschel-Bulkley (yield stress)

(5)\[ \eta(\dot\gamma) = \frac{\tau_0}{\dot\gamma} + K\,\dot\gamma^{\,n-1} \qquad(\text{ideal},\ \dot\gamma>0) \]

with yield stress \(\tau_0\); \(\tau_0=0\) recovers the power-law, and \(\tau_0>0,\,n=1\) is the Bingham plastic \(\eta=\tau_0/\dot\gamma+\mu_p\). The term \(\tau_0/\dot\gamma\to\infty\) as \(\dot\gamma\to0\), so the ideal model cannot be evaluated in unyielded (plug) regions and requires a regularization. Nassu uses the Papanastasiou exponential form Papanastasiou[3], which is smooth and finite everywhere:

(6)\[ \eta(\dot\gamma) = K\,\dot\gamma^{\,n-1} \;+\; \frac{\tau_0\,\big(1 - e^{-m\,\dot\gamma}\big)}{\dot\gamma} \]

The stress-growth exponent \(m\) (units of time) controls the sharpness of the yielded/unyielded transition; \(m\to\infty\) recovers the ideal model. As \(\dot\gamma\to0\) the regularized term tends to the finite value \(m\tau_0\). A bi-viscosity regularization (a large Newtonian plateau below a critical \(\dot\gamma_c\)) is a cheaper \(C^0\) alternative.

Mapping to the LBM relaxation

In lattice units (\(\Delta x=\Delta t=1\), \(c_s^2=1/3\)) the relaxation follows the standard viscosity relation:

(7)\[ \tau = \frac{1}{2} + \frac{\eta(\dot\gamma)}{\rho\,c_s^2\,\Delta t} = \frac{1}{2} + 3\,\nu(\dot\gamma), \qquad \omega = \frac{1}{\tau}, \qquad \nu=\eta/\rho \]

Stability window. The RR-BGK operator requires \(\omega\in(0,2)\), i.e. \(\tau>1/2\), i.e. \(\nu>0\): as \(\nu\to0\), \(\omega\to2\) (marginal); as \(\nu\to\infty\), \(\omega\to0\) (over-damped). Because the power-law and Herschel-Bulkley models drive \(\eta\) out of this window at their tails, the apparent viscosity is hard-clipped to \([\eta_{\min},\eta_{\max}]\) (equivalently \(\omega\) to a sub-interval of \((0,2)\) with margin) after evaluating the model and before forming \(\omega\). Carreau-Yasuda is self-clipping; the clip guards both the divergence end and the vanishing-viscosity (\(\omega\to2\)) end.

Locality, and the honest limit on “no iteration”. The strain rate is read locally from \(\Pi^{\mathrm{neq}}\), directly analogous to the Dong (2008) explicit-\(\omega\) Smagorinsky route (LES chapter). However, since \(S_{\alpha\beta}=-\tfrac{\omega}{2\rho c_s^2\Delta t}Q_{\alpha\beta}\), the recovered \(\dot\gamma=|S|=B/\omega\) still depends on \(\omega\) (with \(B\) a known, \(\omega\)-independent quantity from \(Q_{\alpha\beta}\)). The self-consistency condition

(8)\[ c_s^2\!\left(\frac{1}{\omega}-\frac12\right) = \frac{1}{\rho}\,\eta\!\left(\frac{B}{\omega}\right) \]

is closed-form only for the Newtonian and Smagorinsky (\(\nu_{\mathrm{SGS}}\propto|S|\)) cases; for a general \(\eta(\dot\gamma)\) it is transcendental. Nassu evaluates \(\eta\) at a lagged strain rate (the stored previous-step \(|S|\)), which makes the per-node update fully explicit and local. This is standard generalized-Newtonian LBM practice Gabbanelli et al.[4], Boyd et al.[5]. The locality is the real advantage; an unconditional closed-form update is not claimed.

Coupling with LES

With an SGS model active the effective diffusivity is the sum of the molecular apparent viscosity and the eddy viscosity, both functions of the same resolved \(|S|\):

(9)\[ \nu_{\mathrm{total}}(\dot\gamma) = \nu_{\mathrm{gn}}(\dot\gamma) + \nu_{\mathrm{SGS}}(\dot\gamma), \qquad \nu_{\mathrm{SGS}} = C_S^2\,(\Delta x_n)^2\,|S| \]

so a single strain-rate evaluation feeds both, and \(\omega\) is formed from \(\nu_{\mathrm{total}}\) via Eq. rheo_tau. Additivity of a molecular viscosity and a modelled eddy viscosity is the standard eddy-viscosity closure. Two modeling caveats hold: the Smagorinsky constant \(C_S\) is calibrated on a Newtonian cascade and should be validated, not assumed, on a shear-thinning background; and because \(\eta(\dot\gamma)\) is nonlinear, evaluating the model at the resolved strain rate (\(\eta(\overline{\dot\gamma})\) rather than \(\overline{\eta(\dot\gamma)}\)) is a modeling assumption, not an identity.

Multiblock behaviour

The apparent viscosity is folded into \(\omega_{\mathrm{LES}}\) at the relaxation hook, so every interface-communication site that converts the stored strain rate to a stress through the effective viscosity (\(\mathrm{stress}=2\rho\nu S\) with \(\nu\) from \(\omega_{\mathrm{LES}}\)) remains correct with no change. The one site that assumes a constant molecular viscosity is the coarse-to-fine analytic stress inversion (the LES shear-stress matching, build_S_from_stress), whose closed form solves a quadratic valid only for shear-rate-independent \(\nu_0\). Under a generalized-Newtonian viscosity this inversion must instead evaluate the apparent viscosity at the (lagged) interface strain rate. This is the single genuine variation point introduced by the feature, and continuity of the strain rate and viscosity across a 2:1 refinement interface is therefore a mandatory validation target.

Neutral (Newtonian) limit

Every model collapses exactly to the constant-viscosity solver, so enabling the module with a degenerate parameter set is bit-identical to a Newtonian run:

• Power-law \(n=1\): \(\eta=K\) constant.

• Carreau-Yasuda \(n=1\) (or \(\eta_0=\eta_\infty\)): the exponent \((n-1)/a=0\) gives \(\eta=\eta_0\) constant.

• Herschel-Bulkley \(\tau_0=0,\,n=1\): \(\eta=K\) constant (and \(\tau_0=0\) zeroes the Papanastasiou term for any \(m\)).

In a constant-\(\eta\) state the local update falls through to exactly the existing Newtonian \(\omega\) (no strain-rate read, no clip), and a regression asserts bit-identical output against the Newtonian baseline (subject to the known IBM atomicAdd tolerance where IBM is present).

Validation targets

The fully-developed closed-form profiles below are the quantitative references for the validation cases. Plane channel of half-height \(H\) (walls at \(y=\pm H\)), constant favourable pressure gradient \(G=-\mathrm{d}p/\mathrm{d}x>0\), no-slip walls; pipe of radius \(R\).

T1 - Power-law plane channel Bird et al.[6], Gabbanelli et al.[4]:

(10)\[ u(y) = \frac{n}{n+1}\left(\frac{G}{K}\right)^{1/n} \Big[\,H^{\frac{n+1}{n}} - |y|^{\frac{n+1}{n}}\Big] \]

T2 - Power-law pipe (generalized Hagen-Poiseuille) Bird et al.[6]:

(11)\[ u(r) = \frac{n}{n+1}\left(\frac{G}{2K}\right)^{1/n} \Big[\,R^{\frac{n+1}{n}} - r^{\frac{n+1}{n}}\Big] \]

T3 - Bingham plane channel (yield-stress plug, plug half-width \(y_0=\tau_0/G\), plastic viscosity \(\mu_p\)) Chhabra and Richardson[2]:

(12)\[\begin{split} u(y) = \begin{cases} \dfrac{G}{2\mu_p}(H-y_0)^2, & |y|\le y_0 \\[8pt] \dfrac{1}{\mu_p}\!\left[\dfrac{G}{2}(H^2-y^2) - \tau_0(H-|y|)\right], & y_0\le|y|\le H \end{cases} \end{split}\]

T4 - Herschel-Bulkley plane channel (general \(n\); reduces to T3 at \(n=1,\,K=\mu_p\)) Chhabra and Richardson[2]:

(13)\[ u(y) = \frac{n}{(n+1)\,G\,K^{1/n}} \Big[\,(GH-\tau_0)^{\frac{n+1}{n}} - (G|y|-\tau_0)^{\frac{n+1}{n}}\Big], \quad y_0\le|y|\le H \]

Each reduces to standard Newtonian Poiseuille at \(n=1\), \(\tau_0=0\). Plane Couette flow (uniform shear stress) yields a linear profile independent of \(n\) and \(\tau_0\) and is therefore only a null/consistency check, not a discriminating test.

---

[1] (1,2)

R. Byron Bird, Robert C. Armstrong, and Ole Hassager. Dynamics of Polymeric Liquids, Volume 1: Fluid Mechanics. John Wiley & Sons, 2 edition, 1987. ISBN 978-0-471-80245-7.

[2] (1,2,3)

Rajendra P. Chhabra and John F. Richardson. Non-Newtonian Flow and Applied Rheology: Engineering Applications. Butterworth-Heinemann, 2 edition, 2008. ISBN 978-0-7506-8532-0.

[3]

Tasos C. Papanastasiou. Flows of materials with yield. Journal of Rheology, 31(5):385–404, 1987. doi:10.1122/1.549926.

[4] (1,2)

Susana Gabbanelli, German Drazer, and Joel Koplik. Lattice Boltzmann method for non-Newtonian (power-law) fluids. Physical Review E, 72(4):046312, 2005. doi:10.1103/PhysRevE.72.046312.

[5]

Joshua Boyd, James M. Buick, and Simon Green. Analysis of the Casson and Carreau-Yasuda non-Newtonian blood models in steady and oscillatory flows using the lattice Boltzmann method. Physics of Fluids, 19(9):093103, 2007. doi:10.1063/1.2772250.

[6] (1,2)

R. Byron Bird, Warren E. Stewart, and Edwin N. Lightfoot. Transport Phenomena. John Wiley & Sons, 2 edition, 2002. ISBN 978-0-471-41077-5.