Professor, I have a vague image of non-Newtonian fluids as 'fluids with weird viscosity,' but could you explain them a bit more precisely?

In a Newtonian fluid, the shear stress $\tau$ and the shear rate $\dot{\gamma}$ have a linear relationship: $\tau = \mu \dot{\gamma}$. A fluid for which this linear relationship does not hold is a non-Newtonian fluid. There are many examples around us: blood, paint, ketchup, shampoo, polymer solutions, cement slurry, and so on.

Non-Newtonian Fluid

Q: Fluids with a yield stress seem difficult to handle numerically.

That's correct. Because the position of the yield surface is unknown, regularization techniques are often used. In the Papanastasiou model: $$ \eta(\dot{\gamma}) = \frac{\tau_y[1 - \exp(-m\dot{\gamma})]}{\dot{\gamma}} + K\dot{\gamma}^{n-1} $$ Increasing the parameter $m$ brings the behavior closer to an ideal Bingham plastic, but it also increases numerical stiffness. A typical range is $m = 100\text{--}1000\,\text{s}$.

Category: 流体解析（CFD） | Integrated 2026-04-06

CAE visualization for non newtonian theory - technical simulation diagram

非ニュートン流体

Theory and Physics

What is a Non-Newtonian Fluid?

🧑‍🎓

Professor, my image of a non-Newtonian fluid is just "a fluid with weird viscosity." Could you explain it a bit more properly?

🎓

For Newtonian fluids, shear stress $\tau$ and shear rate $\dot{\gamma}$ have a linear relationship $\tau = \mu \dot{\gamma}$. Fluids for which this linear relationship does not hold are non-Newtonian fluids. There are many around us. Blood, paint, ketchup, shampoo, polymer solutions, cement slurry, etc.

🎓

Non-Newtonian fluids are broadly classified as follows.

Type	Characteristics	Examples
Shear-thinning (Pseudoplastic)	Shear rate ↑ → Viscosity ↓	Paint, blood, polymer solutions
Shear-thickening (Dilatant)	Shear rate ↑ → Viscosity ↑	Cornstarch-water mixture
Bingham Plastic	Does not deform below yield stress	Toothpaste, chocolate
Viscoelastic Fluid	Possesses both viscosity and elasticity	Polymer melt, DNA solution
Thixotropy	Viscosity decreases over time	Paint, yogurt

🧑‍🎓

So ketchup being hard to get out is also non-Newtonian! Shaking it lowers the viscosity, making it easier to pour.

Constitutive Equation (Viscosity Model)

🎓

The core of non-Newtonian fluids lies in the formulation of the viscosity function $\eta(\dot{\gamma})$, which depends on shear rate. Let me introduce the main models.

🎓

1. Power-law Model

It is the simplest model.

$$ \eta(\dot{\gamma}) = K \dot{\gamma}^{n-1} $$

$K$ is the consistency index, $n$ is the power-law index. $n < 1$ for shear-thinning, $n > 1$ for shear-thickening. It reduces to a Newtonian fluid when $n = 1$.

🎓

2. Carreau Model

A model that improves upon the shortcomings of the Power-law (unrealistic behavior at low/high shear rates).

$$ \eta(\dot{\gamma}) = \eta_\infty + (\eta_0 - \eta_\infty)\left[1 + (\lambda\dot{\gamma})^2\right]^{(n-1)/2} $$

$\eta_0$ is the zero-shear viscosity, $\eta_\infty$ is the infinite-shear viscosity, $\lambda$ is the relaxation time.

🎓

3. Herschel-Bulkley Model

For fluids with yield stress (a generalization of Bingham plastic).

$$ \tau = \tau_y + K\dot{\gamma}^n \quad (\tau > \tau_y) $$

$\tau_y$ is the yield stress. For $\tau < \tau_y$, there is no flow (rigid behavior). It reduces to the Bingham model when $n = 1$.

🧑‍🎓

Fluids with yield stress seem numerically tricky to handle.

🎓

Exactly. Since the location of the yield surface is unknown, regularization techniques are often used. In the Papanastasiou model:

$$ \eta(\dot{\gamma}) = \frac{\tau_y[1 - \exp(-m\dot{\gamma})]}{\dot{\gamma}} + K\dot{\gamma}^{n-1} $$

Increasing the parameter $m$ brings it closer to ideal Bingham behavior, but numerically increases stiffness. $m = 100\text{--}1000\,\text{s}$ is a typical range.

Coffee Break Trivia

The Foundation of Non-Newtonian Fluid Theory—Bingham and Ostwald-de Waele (1906–1929)

Systematic research on non-Newtonian fluids began in the early 20th century. Eugene Bingham (1916) measured the "Yield Stress" of mud slurries and pastes and proposed the "Bingham Plastic" model, which behaves linearly after yielding. Meanwhile, Ostwald (1925) and de Waele (1923) independently formulated the "Power Law" fluid (τ=K·γ̇ⁿ), providing a unified description for shear-thinning (n<1: blood, polymer solutions) and shear-thickening (n>1: cornstarch-water mixture). These century-old models are still implemented as the foundation of CFD material model libraries and served as the starting point for more refined Cross, Casson, and Herschel-Bulkley models.

Physical Meaning of Each Term

Temporal Term $\partial(\rho\phi)/\partial t$: Imagine the moment you turn on a faucet. At first, water comes out in an unstable, spluttering manner, but after a while, it becomes a steady flow, right? This "period of change" is described by the temporal term. The pulsation of blood flow from a heartbeat, or the flow fluctuation each time an engine valve opens and closes—all are unsteady phenomena. So what is steady-state analysis? It looks only at "after sufficient time has passed and the flow has settled down"—meaning setting this term to zero. Since computational cost drops significantly, solving first in steady-state is a basic CFD strategy.
Convection Term $\nabla \cdot (\rho \mathbf{u} \phi)$: If you drop a leaf into a river, what happens? It gets carried downstream by the flow, right? This is "convection"—the effect where fluid motion transports things. Warm air from a heater reaching the far corner of a room is also because the "carrier," air, transports heat via convection. Here's the interesting part—this term contains "velocity × velocity," making it nonlinear. That is, as flow speed increases, this term rapidly strengthens, making control difficult. This is the root cause of turbulence. A common misconception: "Convection and conduction are similar" → They are completely different! Convection is carried by flow, conduction is transmitted by molecules. There's an order of magnitude difference in efficiency.
Diffusion Term $\nabla \cdot (\Gamma \nabla \phi)$: Have you ever put milk in coffee and left it? Even without stirring, after a while they naturally mix. That's molecular diffusion. Now, next question—honey and water, which flows more easily? Obviously water, right? Honey has high viscosity ($\mu$), so it flows poorly. When viscosity is large, the diffusion term becomes strong, and the fluid moves in a "thick" manner. In low Reynolds number flow (slow, viscous), diffusion is dominant. Conversely, in high Re number flow, convection overwhelms and diffusion plays a supporting role.
Pressure Term $-\nabla p$: When you push the plunger of a syringe, liquid shoots out forcefully from the needle tip, right? Why? Because the piston side is high pressure, the needle tip is low pressure—this pressure difference provides the force pushing the fluid. Dam discharge works on the same principle. On a weather map, where isobars are densely packed? That's right, strong winds blow. "Flow arises where there is a pressure difference"—this is the physical meaning of the pressure term in the Navier-Stokes equations. A point of confusion here: "Pressure" in CFD is often gauge pressure, not absolute pressure. If results go wrong immediately after switching to compressible analysis, mixing up absolute/gauge pressure might be the cause.
Source Term $S_\phi$: Warmed air rises—why? Because it becomes lighter (lower density) than its surroundings, so buoyancy pushes it up. This buoyancy is added to the equation as a source term. Other examples: chemical reaction heat from a gas stove flame, Lorentz force acting on molten metal in a factory's electromagnetic pump... These are all actions that "inject energy or force into the fluid from the outside," expressed by the source term. What happens if you forget the source term? In natural convection analysis, forgetting buoyancy means the fluid doesn't move at all—a physically impossible result where warm air doesn't rise in a room with the heater on in winter.

Assumptions and Applicability Limits

Continuum Assumption: Valid for Knudsen number Kn < 0.01 (molecular mean free path ≪ characteristic length)
Newtonian Fluid Assumption: Shear stress and strain rate have a linear relationship (viscosity model required for non-Newtonian fluids)
Incompressibility Assumption (for Ma < 0.3): Treat density as constant. For Mach number ≥ 0.3, compressibility effects must be considered.
Boussinesq Approximation (Natural Convection): Density variation is considered only in the buoyancy term; constant density is used in other terms.
Non-applicable Cases: Rarefied gas (Kn > 0.1), supersonic/hypersonic flow (shock capturing required), free surface flow (VOF/Level Set, etc. required)

Dimensional Analysis and Unit Systems

Variable	SI Unit	Notes / Conversion Memo
Velocity $u$	m/s	When converting from volumetric flow rate for inlet conditions, pay attention to cross-sectional area units.
Pressure $p$	Pa	Distinguish between gauge and absolute pressure. Use absolute pressure for compressible analysis.
Density $\rho$	kg/m³	Air: ~1.225 kg/m³@20°C, Water: ~998 kg/m³@20°C
Viscosity Coefficient $\mu$	Pa·s	Be careful not to confuse with kinematic viscosity $\nu = \mu/\rho$ [m²/s]
Reynolds Number $Re$	Dimensionless	$Re = \rho u L / \mu$. Indicator for laminar/turbulent transition.
CFL Number	Dimensionless	$CFL = u \Delta t / \Delta x$. Directly related to time-step stability.

Numerical Methods and Implementation

Numerical Methods for Non-Newtonian Fluids

🧑‍🎓

When solving non-Newtonian fluids with CFD, what's different from Newtonian fluids?

🎓

The fundamental difference is that viscosity depends on the velocity field. The viscous term in the Navier-Stokes equations becomes nonlinear.

$$ \nabla \cdot [\eta(\dot{\gamma})\dot{\boldsymbol{\gamma}}] \quad \text{where} \quad \dot{\gamma} = \sqrt{2\mathbf{D}:\mathbf{D}} $$

$\mathbf{D} = \frac{1}{2}(\nabla\mathbf{u} + \nabla\mathbf{u}^T)$ is the strain rate tensor.

🎓

As a solution method, Picard iteration (successive substitution method) becomes the basis.

1. Calculate $\dot{\gamma}$ from the previous step's velocity field

2. Update $\eta(\dot{\gamma})$

3. Solve the N-S equations with the updated viscosity to obtain a new velocity field

4. Repeat steps 1-3 until convergence

🧑‍🎓

So there's one more iteration loop compared to Newtonian fluids.

🎓

Correct. The difficulty of non-Newtonian fluid computation lies in this outer iteration loop being hard to converge. Special care is needed when shear-thinning is strong (small $n$) or when yield stress is present.

Generalized Reynolds Number

🎓

For non-Newtonian fluids, the definition of Reynolds number also changes. For pipe flow of Power-law fluids, the Metzner-Reed generalized Reynolds number is used.

$$ \text{Re}_{\text{MR}} = \frac{\rho U^{2-n} D^n}{K \cdot 8^{n-1} \left(\frac{3n+1}{4n}\right)^n} $$

🎓

Laminar-turbulent transition occurs around $\text{Re}_{\text{MR}} \approx 2100$ (almost the same critical value as Newtonian fluids). However, post-transition turbulent behavior differs significantly from Newtonian fluids.

🧑‍🎓

The generalized Re number formula is quite complex.

Discretization in the Finite Volume Method

🎓

Let me explain the discretization points to note when solving non-Newtonian fluids with CFD.

🎓

Face value interpolation of viscosity becomes important. When interpolating viscosity calculated at cell centers to faces:

Harmonic Mean: Recommended for interfaces where viscosity changes sharply. $\eta_f = \frac{2\eta_L \eta_R}{\eta_L + \eta_R}$
Arithmetic Mean: For cases where viscosity variation is gradual. $\eta_f = \frac{\eta_L + \eta_R}{2}$

🎓

For shear-thinning fluids, viscosity can differ by orders of magnitude between near-wall regions and the channel center. For example, in injection molding of polymer melts, $\eta$ can vary in the range of $10^1 \sim 10^4\,\text{Pa}\cdot\text{s}$. To handle this sharp change, the mesh near walls needs to be sufficiently refined.

Power-law $n$	Viscosity Ratio $\eta_{\text{center}}/\eta_{\text{wall}}$	Mesh Requirement
0.8	~3	Mild refinement sufficient
0.5	~30	10+ layers near wall recommended
0.2	~1000	Very fine mesh required

🧑‍🎓

The smaller $n$ is, the stricter the mesh requirement becomes. Convergence probably gets harder too.

🎓

Exactly. A practical technique is to lower the relaxation factor to 0.3-0.5 and also apply relaxation to the viscosity update.

Coffee Break Trivia

Numerical Stability of Non-Newtonian Fluid CFD—Singularity Handling for Yield-Stress Fluids (Bingham Fluids)

In CFD for Bingham fluids (yield stress τ_y>0), the discontinuity in material properties—"solid below yield stress, liquid above"—causes numerical instability. In the limit of shear rate γ̇→0, the viscosity of the Bingham model diverges to infinity, causing Newton's method to fail to converge. A practical solution is Papanastasiou's (1987) regularization approximation: τ = (τ_y/γ̇)(1-exp(-mγ̇)) + ηγ̇, where parameter m (typical values 10³–10⁵) balances convergence and accuracy. If m is too large, numerical stiffness causes divergence; if too small, the yield stress effect is underestimated. An appropriate m value is around 100 times τ_y/η (inverse of viscous time scale), and it is recommended to verify this along with mesh refinement.

Upwind Scheme

1st-order upwind: Large numerical diffusion but stable. 2nd-order upwind: Improved accuracy but risk of oscillations. Essential for high Reynolds number flows.

Central Differencing

2nd-order accurate, but numerical oscillations occur for Pe number > 2. Suitable for low Reynolds number, diffusion-dominated flows.

TVD Scheme (MUSCL, QUICK, etc.)

Maintains high accuracy while suppressing numerical oscillations via limiter functions. Effective for capturing shocks and steep gradients.

Finite Volume Method vs Finite Element Method

FVM: Naturally satisfies conservation laws. Mainstream in CFD. FEM: Advantageous for complex shapes and multiphysics. Mesh-free methods like SPH are also developing.

CFL Condition (Courant Number)

Explicit method: CFL ≤ 1 is the stability condition. Implicit method: Stable even for CFL > 1, but affects accuracy and iteration count. LES: CFL ≈ 1 recommended. Physical meaning: Information should not travel more than one cell per time step.

Residual Monitoring

Convergence is judged when residuals for continuity, momentum, and energy each drop by 3–4 orders of magnitude. The mass conservation residual is particularly important.

Relaxation Factor

Pressure: 0.2–0.3, Velocity: 0.5–0.7 are typical initial values. Lower the factor if diverging. Increase after convergence to accelerate.