Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

By the end of this section, you will be able to:

Explain what viscosity physically is — momentum diffusion between adjacent fluid layers — and why a single number $\mu$ captures it for many real fluids.
Build the 1D Newton law of viscosity $\tau = \mu \, \mathrm{d}u/\mathrm{d}y$ from the slope of a velocity profile.
Promote that scalar relation to the full $3 \times 3$ Cauchy stress tensor $\boldsymbol{\sigma}$ , and decompose it into pressure + deviatoric parts.
Compute the strain-rate tensor $\mathbf{D} = \tfrac{1}{2}(\nabla \mathbf{v} + (\nabla \mathbf{v})^{T})$ and use the Newtonian constitutive law $\boldsymbol{\sigma} = -p\,\mathbf{I} + 2\mu\,\mathbf{D}$ .
Connect the divergence of the stress tensor to the viscous term $\mu\,\nabla^2 \mathbf{v}$ in the Navier–Stokes momentum equation.
Implement the stress tensor in plain Python (finite differences) and again in PyTorch (exact autograd) for a worked Couette-flow case.

Intuition: Why a Tensor at All?

The key idea. Pressure pushes equally in all directions and needs only one number per point. Viscous forces depend on which face you push on and which direction you push. Two directions → a matrix → a tensor.

Imagine a tiny cube of fluid sitting in a flow. Each of its six faces feels a force from the neighbouring fluid. That force has both a normal component (squeeze or pull) and two shear components (sliding sideways). Three faces × three force components gives nine numbers — the entries of the stress tensor $\sigma_{ij}$ .

The first index of $\sigma_{ij}$ labels the face (the normal direction), the second labels the force component. So $\sigma_{xy}$ reads: "on the face whose normal points in $x$ , the $y$ -component of the force per unit area."

Why "tensor", not just "matrix"? Tensors are matrices that transform in a specific way when you rotate your axes. Pressure stays the same under any rotation. Shear stresses mix in a predictable pattern — the same pattern as

R\,\sigma\,R^{T}

. That well-defined behaviour under rotation is what makes

\boldsymbol{\sigma}

a physical object, not just an array of numbers.

Newton's Law of Viscosity (1D)

Start with the simplest possible setup — the one Newton himself analyzed in the Principia. Two large flat plates sandwich a thin layer of fluid. The bottom plate is fixed. The top plate moves with steady velocity $U$ . After a moment of transients, the fluid in between settles into a steady linear velocity profile:

u(y) = U\,\dfrac{y}{h}, \qquad 0 \le y \le h

Each fluid layer slides past the one below at a tiny relative speed. There must be a force per unit area between them, because a real fluid resists shearing — that's viscosity. Newton hypothesised that for a great many fluids (water, air, oil) that force is exactly proportional to the slope of the velocity profile:

\boxed{\;\tau \;=\; \mu \, \dfrac{\mathrm{d}u}{\mathrm{d}y}\;}

The constant of proportionality $\mu$ is the dynamic viscosity. Its units fall right out of dimensional analysis:

[\mu] = \dfrac{[\tau]}{[\mathrm{d}u/\mathrm{d}y]} = \dfrac{\mathrm{Pa}}{\mathrm{s}^{-1}} = \mathrm{Pa\cdot s}

Intuition: viscosity as momentum diffusion

Pick a horizontal surface inside the fluid. Molecules constantly hop across it. The molecules coming from above carry slightly more horizontal momentum than those coming from below, because the layer above moves slightly faster. Net transfer of horizontal momentum downward = a horizontal force pulling the lower fluid along. That microscopic momentum exchange is macroscopic viscous stress.

The viscosity $\mu$ measures how efficient that hopping is. Air is bad at it ( $\mu \approx 1.8 \times 10^{-5}$ Pa·s). Honey is great at it ( $\mu \approx 10$ Pa·s).

Fluid	μ at 20°C (Pa·s)	How "sticky"?
Air	1.8 × 10⁻⁵	Almost frictionless
Water	1.0 × 10⁻³	Easy to stir
Olive oil	≈ 0.08	Pours slowly
Glycerine	≈ 1.4	Like honey
Honey	≈ 10	Hangs from a spoon
Pitch	≈ 2 × 10⁸	Looks solid for years

Interactive: Couette Flow Explorer

Before doing any heavier math, get a feel for the relationship between the velocity profile, the shear rate, and the stress. Drag the sliders below and watch the slope, the arrows, and the numbers change together.

Loading interactive shear-flow visualizer…

Things to try

Double $U$ . The slope doubles, so $\tau$ doubles too — that is the linearity of Newton's law in action.
Halve the gap $h$ . The slope doubles even though $U$ is unchanged — squeezing the fluid into a thinner layer makes it work harder.
Switch to shear-thinning. Increase $U$ — the stress grows, but more slowly than linearly. That is why ketchup gets easier to pour the harder you shake it.

The Cauchy Stress Tensor

Real flows are three-dimensional, and forces inside the fluid have to be described on every possible face orientation. Cauchy's breakthrough (1822) was the realisation that you do not need a different rule for every plane. Knowing the force per unit area on three mutually perpendicular faces is enough — any other face is a weighted combination.

\boldsymbol{\sigma} = \begin{pmatrix} \sigma_{xx} & \sigma_{xy} & \sigma_{xz} \\ \sigma_{yx} & \sigma_{yy} & \sigma_{yz} \\ \sigma_{zx} & \sigma_{zy} & \sigma_{zz} \end{pmatrix}

Cauchy's stress theorem says that the traction vector $\mathbf{t}$ on a face with unit normal $\mathbf{n}$ is just $\mathbf{t} = \boldsymbol{\sigma}\,\mathbf{n}$ . Nine numbers describe stress at any orientation, at any point.

Symmetry: angular-momentum conservation

For a fluid that is not spinning with infinite angular acceleration (i.e. any real fluid), conservation of angular momentum on a vanishing element forces the stress tensor to be symmetric:

\sigma_{ij} = \sigma_{ji}

That cuts the nine entries down to six independent ones. Geometrically, the symmetric tensor diagonalises in some rotated frame — the principal directions of stress. In those axes there are no shear forces, only pure normal pulls or pushes.

Pressure + Deviatoric Decomposition

Every symmetric tensor splits uniquely into an isotropic part (a scalar times the identity) and a deviatoric traceless part:

\boldsymbol{\sigma} \;=\; \underbrace{-\,p\,\mathbf{I}}_{\text{isotropic (pressure)}} \;+\; \underbrace{\boldsymbol{\tau}}_{\text{deviatoric (viscous)}}

We define the pressure as minus one third of the trace (the negative sign keeps it positive for a fluid that is being squeezed):

p \;=\; -\tfrac{1}{3}\,\operatorname{tr}(\boldsymbol{\sigma}) \;=\; -\tfrac{1}{3}\,(\sigma_{xx} + \sigma_{yy} + \sigma_{zz})

Everything else — every off-diagonal entry and the diagonal's deviation from being all-equal — is the viscous stress tensor $\boldsymbol{\tau}$ . For an inviscid fluid $\boldsymbol{\tau} = \mathbf{0}$ and the Cauchy equation collapses to Euler's equation.

Why this split matters. Pressure does no work on the fluid's shape — it only changes its size (if it is compressible at all). The deviatoric part is what actually deforms fluid elements and dissipates kinetic energy into heat. The Newtonian assumption is a statement about that part alone.

The Strain-Rate Tensor

What does the fluid actually do at a point? The first thing we can ask is: how does the velocity change as we move infinitesimally? The Taylor expansion of $\mathbf{v}$ around a point introduces the velocity-gradient tensor $\mathbf{L} = \nabla \mathbf{v}$ with components $L_{ij} = \partial v_i / \partial x_j$ .

Split into deformation + rotation

Any matrix splits into a symmetric and an antisymmetric part:

\mathbf{L} \;=\; \underbrace{\tfrac{1}{2}\bigl(\mathbf{L} + \mathbf{L}^{T}\bigr)}_{\displaystyle \mathbf{D}\ \text{strain-rate}} \;+\; \underbrace{\tfrac{1}{2}\bigl(\mathbf{L} - \mathbf{L}^{T}\bigr)}_{\displaystyle \mathbf{W}\ \text{vorticity}}

The symmetric part $\mathbf{D}$ describes how a fluid element is stretching or shearing — the geometric distortion. The antisymmetric part $\mathbf{W}$ is a rigid-body rotation of the element. Rigid rotation costs no internal friction, so only $\mathbf{D}$ can produce viscous stress.

Strain rate D — fluid deforms

$D_{ij} = \tfrac{1}{2}\bigl(\partial_j v_i + \partial_i v_j\bigr)$

Diagonal entries: rate of stretching along that axis. Off-diagonal entries: rate of angular shear between two perpendicular axes. This is what costs energy.

Vorticity W — fluid spins

$W_{ij} = \tfrac{1}{2}\bigl(\partial_j v_i - \partial_i v_j\bigr)$

Encodes the local angular velocity $\boldsymbol{\omega} = \tfrac{1}{2}\nabla \times \mathbf{v}$ . No internal friction here.

Trace and incompressibility

The trace of $\mathbf{D}$ is $\operatorname{tr}(\mathbf{D}) = \partial_x u + \partial_y v + \partial_z w = \nabla \cdot \mathbf{v}$ — the divergence of the velocity. For an incompressible fluid this is zero, so $\mathbf{D}$ is traceless and the strain rate captures pure shearing, no volume change.

The Newtonian Constitutive Law

We have a kinematic object ( $\mathbf{D}$ ) and a dynamic object ( $\boldsymbol{\tau}$ ). A constitutive law connects them. Newton's assumption — generalised by Stokes to 3D — is the simplest possible: stress is a linear, isotropic function of strain rate.

The general (compressible) form

Isotropy plus linearity force the constitutive law to be (up to a sign convention):

\boldsymbol{\sigma} \;=\; -p\,\mathbf{I} \;+\; 2\mu\,\mathbf{D} \;+\; \lambda\,(\nabla \cdot \mathbf{v})\,\mathbf{I}

The first term is pressure. The second is the shearing response with viscosity $\mu$ . The third (volumetric viscosity $\lambda$ ) only matters if the fluid is being compressed or expanded — for almost all liquids we assume $\nabla \cdot \mathbf{v} = 0$ and the third term vanishes outright. Stokes' hypothesis $\lambda = -\tfrac{2}{3}\mu$ is the standard choice for compressible flows when one wants the mean normal stress to equal $-p$ .

The incompressible form

Setting $\nabla \cdot \mathbf{v} = 0$ leaves the clean, memorable form used throughout the rest of this chapter:

\boxed{\;\boldsymbol{\sigma} \;=\; -p\,\mathbf{I} \;+\; 2\mu\,\mathbf{D}\;}

Component-wise:

\sigma_{ij} \;=\; -p\,\delta_{ij} \;+\; \mu\,\Bigl(\dfrac{\partial v_i}{\partial x_j} \;+\; \dfrac{\partial v_j}{\partial x_i}\Bigr)

That single line contains all of viscous fluid mechanics for water-like fluids. Every more advanced model — turbulence closures, non-Newtonian rheology, viscoelasticity — replaces or generalises just this constitutive equation.

From Stress Tensor to Navier–Stokes

Cauchy's momentum equation says that the net force on a fluid element is the divergence of the stress tensor plus body forces:

\rho\,\dfrac{D\mathbf{v}}{Dt} \;=\; \nabla \cdot \boldsymbol{\sigma} \;+\; \rho\,\mathbf{g}

Plug in $\boldsymbol{\sigma} = -p\,\mathbf{I} + 2\mu\,\mathbf{D}$ and grind through the derivatives. Two facts make it collapse to something familiar:

$\nabla \cdot (-p\,\mathbf{I}) = -\nabla p$ — the divergence of a scalar-times-identity is just the gradient.
$\nabla \cdot (2\mu\,\mathbf{D}) = \mu\,\nabla^{2}\mathbf{v} + \mu\,\nabla(\nabla \cdot \mathbf{v}) = \mu\,\nabla^{2}\mathbf{v}$ when $\nabla \cdot \mathbf{v} = 0$ .

Substituting yields the incompressible Navier–Stokes momentum equation:

\rho\,\dfrac{D\mathbf{v}}{Dt} \;=\; -\nabla p \;+\; \mu\,\nabla^{2}\mathbf{v} \;+\; \rho\,\mathbf{g}

The big picture. The mysterious

\mu\,\nabla^{2}\mathbf{v}

term you met in Section 30.1 was never "just a Laplacian." It is literally the divergence of the deviatoric stress tensor — a force per unit volume that arises from layers of fluid rubbing past each other.

Worked Example: Couette Flow by Hand

Let's compute every tensor for the flow in the explorer. Numbers: $U = 2\ \mathrm{m/s}$ , $h = 1\ \mathrm{m}$ , $\mu = 1\ \mathrm{Pa\cdot s}$ . Velocity field: $\mathbf{v} = (Uy/h,\,0,\,0)$ .

📝 Click to expand: full hand calculation (try it yourself first)

Step 1 — velocity gradient L = ∇v

Each entry of $\mathbf{L}$ is $L_{ij} = \partial v_i/\partial x_j$ . Only $v_x = Uy/h$ is non-zero, and it depends only on $y$ . So a single partial survives: $\partial v_x/\partial y = U/h = 2/1 = 2\ \mathrm{s}^{-1}$ . Every other partial is zero. Therefore

\mathbf{L} = \begin{pmatrix} 0 & 2 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \ \mathrm{s}^{-1}.

Step 2 — strain rate D = ½(L + Lᵀ)

Take the symmetrization. The single non-zero entry $L_{xy} = 2$ splits across both off-diagonal positions:

\mathbf{D} = \tfrac{1}{2}\begin{pmatrix} 0 & 2 & 0 \\ 2 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} = \begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}\ \mathrm{s}^{-1}.

Step 3 — vorticity W = ½(L − Lᵀ)

\mathbf{W} = \tfrac{1}{2}\begin{pmatrix} 0 & 2 & 0 \\ -2 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} = \begin{pmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}\ \mathrm{s}^{-1}.

The angular velocity vector is $\boldsymbol{\omega} = \tfrac{1}{2}\nabla \times \mathbf{v} = (0, 0, -\tfrac{1}{2}\,U/h) = (0,0,-1)$ rad/s — every fluid element spins clockwise (looking down the z-axis) at 1 rad/s. Yet there is no viscous stress associated with that rotation, because $\mathbf{W}$ drops out of the constitutive law.

Step 4 — stress tensor σ = −pI + 2μD

Take the gauge pressure to be zero. Then

\boldsymbol{\sigma} = 2 \cdot 1 \cdot \begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} = \begin{pmatrix} 0 & 2 & 0 \\ 2 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}\ \mathrm{Pa}.

Step 5 — read off the shear stress on the plates

The top plate's outward normal is $\mathbf{n} = (0,1,0)$ . The traction it feels from the fluid is

\mathbf{t} = \boldsymbol{\sigma}\,\mathbf{n} = (2,\,0,\,0)\ \mathrm{Pa}.

So the fluid pulls the plate in the $+x$ direction with 2 N/m². If you keep the plate moving you need to push it forward with 2 N/m² × area. That is exactly Newton's 1D answer $\tau = \mu U/h = 2$ — but now it falls out of a full tensorial calculation. The interactive explorer above is producing the same number.

Step 6 — check the Navier–Stokes momentum equation

For this steady, fully developed flow the material derivative is zero. Body force is zero (or balanced by the unknown gauge pressure). The viscous term is $\mu \nabla^{2}\mathbf{v}$ . Compute it: $\nabla^{2} v_x = \partial^{2}(Uy/h)/\partial y^{2} = 0$ . Every term vanishes — momentum is conserved trivially. That is why Couette flow is the simplest exact solution of Navier–Stokes you can write down.

Python: Computing the Stress Tensor

Now lift the hand calculation into code. The Python version uses central differences — pedestrian but transparent. Every value matches what we computed analytically.

Stress tensor for Couette flow (NumPy + finite differences)

🐍python

Explanation(8)

Code(46)

1Import NumPy

We need NumPy for the small 3x3 matrices that hold the velocity gradient, strain rate, and stress tensors. Nothing else is required — the math is just elementary array operations.

5Pick a concrete flow: plane Couette

U = 2 m/s is the top-plate speed, h = 1 m is the gap, mu = 1 Pa·s is the dynamic viscosity (close to glycerine). Together they fix the shear rate gamma_dot = U/h = 2 s⁻¹ and the expected shear stress tau = mu * gamma_dot = 2 Pa. We will use these numbers as a hand-checkable target.

EXAMPLE

U = 2.0
h = 1.0
mu = 1.0
gamma_dot = U/h = 2.0 s^-1
tau (expected) = mu * gamma_dot = 2.0 Pa

10Velocity field u(y) = U y / h

For plane Couette flow the velocity points purely in the x-direction and depends only on y. The bottom plate at y = 0 is stationary, the top plate at y = h moves with speed U, and the linear interpolation between them is what a steady Newtonian fluid produces. The function returns (u, v, w) so downstream code can treat the field as 3D.

EXAMPLE

velocity([0.0, 0.5, 0.0]) = [U * 0.5 / 1.0, 0, 0] = [1.0, 0, 0] m/s

14Velocity-gradient tensor L[i,j] = ∂u_i/∂x_j

L is the 3x3 Jacobian of the velocity field. We approximate each column by a central difference: bump x_j by ±eps, recompute v, and divide by 2*eps. For our flow only column j = 1 (the y derivative of u) is nonzero, so we expect L to have a single non-zero entry at L[0,1] = U/h = 2.0.

EXAMPLE

Expected L =
[[ 0   U/h   0 ]      [[ 0   2   0 ]
 [ 0    0    0 ]  =   [ 0   0   0 ]
 [ 0    0    0 ]]      [ 0   0   0 ]]

22Strain-rate tensor D = ½(L + Lᵀ)

D is the symmetric part of L. It captures the actual *deformation* of fluid elements — stretching along principal axes — independently of any rigid rotation. For our shear flow, L is non-symmetric (only L[0,1] is non-zero), so its symmetrization spreads that single entry across both L[0,1] and L[1,0], each becoming U/(2h).

EXAMPLE

D = 0.5 * (L + L.T)
D = [[ 0      U/(2h)  0 ]
     [ U/(2h)  0     0 ]
     [ 0       0     0 ]]

At our numbers:  D[0,1] = D[1,0] = 1.0 s^-1

27Stress tensor σ = -p I + 2μ D (Newtonian, incompressible)

This is the constitutive law. Pressure p contributes only on the diagonal because it pushes equally in every direction. The viscous part 2μD lives in the off-diagonal entries whenever the fluid is shearing. With pressure = 0 (gauge) the diagonal of σ vanishes and we are left with a pure shear stress tensor whose off-diagonal entry σ[0,1] = 2μ * D[0,1] = μ * (U/h) = 2.0 Pa.

EXAMPLE

sigma = -0 * I + 2 * mu * D
      = 2 * 1.0 * [[0,   1,   0],
                   [1,   0,   0],
                   [0,   0,   0]]
      = [[0,   2,   0],
         [2,   0,   0],
         [0,   0,   0]]
sigma[0,1] = tau_xy = 2.0 Pa  ✓

32Pick the evaluation point

We sit at the middle of the gap, point = (0, 0.5, 0). Because Couette flow is *uniform* in x and z and the strain rate is constant in y, every interior point gives the same tensors — choosing the midpoint is just for readability.

EXAMPLE

p_point = [0, 0.5, 0]
u at p_point = U * 0.5 / 1.0 = 1.0 m/s

38Print the three tensors

All three tensors should match the analytic targets to ~7 decimals — finite differences with eps = 1e-6 are more than accurate enough for a linear velocity field. The final scalar σ[0,1] is what an engineer would actually measure with a couette-cell rheometer: the drag force per unit plate area.

EXAMPLE

Velocity gradient L =
[[0. 2. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]
Strain-rate D =
[[0. 1. 0.]
 [1. 0. 0.]
 [0. 0. 0.]]
Stress tensor sigma =
[[0. 2. 0.]
 [2. 0. 0.]
 [0. 0. 0.]]
Shear stress tau_xy = 2.0000 Pa

38 lines without explanation

1import numpy as np
2
3# Velocity field for steady plane Couette flow:
4#   u(y) = U * y / h,   v = 0,   w = 0
5U = 2.0   # top-plate velocity (m/s)
6h = 1.0   # gap between plates (m)
7mu = 1.0  # dynamic viscosity (Pa·s)
8
9def velocity(point):
10    """Return v = (u, v, w) at a 3D point."""
11    x, y, z = point
12    return np.array([U * y / h, 0.0, 0.0])
13
14def velocity_gradient(point, eps=1e-6):
15    """Numerical 3x3 Jacobian L[i,j] = du_i/dx_j via central differences."""
16    L = np.zeros((3, 3))
17    for j in range(3):
18        e = np.zeros(3); e[j] = eps
19        L[:, j] = (velocity(point + e) - velocity(point - e)) / (2 * eps)
20    return L
21
22def strain_rate(point):
23    """Symmetric part: D = 0.5 (L + L^T)."""
24    L = velocity_gradient(point)
25    return 0.5 * (L + L.T)
26
27def stress_tensor(point, pressure):
28    """Newtonian, incompressible:  sigma = -p I + 2 mu D."""
29    D = strain_rate(point)
30    return -pressure * np.eye(3) + 2 * mu * D
31
32# Evaluate everything at a single point inside the fluid:
33p_point = np.array([0.0, 0.5, 0.0])
34pressure = 0.0  # gauge pressure (no pressure gradient in this flow)
35
36L = velocity_gradient(p_point)
37D = strain_rate(p_point)
38sigma = stress_tensor(p_point, pressure)
39
40print("Velocity gradient L =")
41print(L)
42print("Strain-rate D =")
43print(D)
44print("Stress tensor sigma =")
45print(sigma)
46print(f"Shear stress tau_xy = {sigma[0, 1]:.4f} Pa")

Sanity check. The printed tau_xy = 2.0000 Pa agrees with

\;\mu\,(U/h) = 1\cdot 2 = 2

. The diagonal of

\boldsymbol{\sigma}

is all zeros because we chose gauge pressure

p = 0

; in a real engineering calculation you would replace it with the local pressure measured by a manometer.

PyTorch: Strain Rate via Autograd

The finite-difference version is fine for a flow we know analytically. For more complex flows — or as a building block inside a physics-informed neural network — we want exact derivatives. PyTorch's torch.autograd.functional.jacobian gives us $\mathbf{L} = \nabla \mathbf{v}$ in one call, to full machine precision.

Strain rate, vorticity & stress tensor with PyTorch autograd

🐍python

Explanation(7)

Code(33)

1Import PyTorch

PyTorch lets us replace the finite-difference approximation with *exact* derivatives (to machine precision) using reverse-mode autodiff. Same Couette-flow setup, same physics, but a cleaner pipeline: define the velocity as a differentiable function, then ask autograd for its Jacobian.

3Wrap the physical constants as tensors

U, h, mu live as 0-dim PyTorch tensors so that all subsequent algebra stays on the same tensor type. They have no requires_grad — they are parameters of the *flow*, not coordinates we differentiate with respect to.

EXAMPLE

U = tensor(2.0), h = tensor(1.0), mu = tensor(1.0)

7Velocity field as a torch function

velocity_field takes a 3-vector point and returns the stacked velocity vector (u, v, w). We use torch.stack rather than a Python list so the output is itself a differentiable tensor — autograd needs that to build a graph.

EXAMPLE

velocity_field(tensor([0.0, 0.5, 0.0])) = tensor([1.0, 0.0, 0.0])

15Exact Jacobian via torch.autograd.functional.jacobian

One line replaces the entire central-difference loop from the Python version. autograd.functional.jacobian sweeps every output component, differentiates it with respect to every input coordinate, and stacks the results into a (3, 3) tensor. For our flow it returns a sparse pattern with a single non-zero at L[0, 1] = 2.

EXAMPLE

L = velocity_gradient([0, 0.5, 0])
L = tensor([[0., 2., 0.],
            [0., 0., 0.],
            [0., 0., 0.]])

19Symmetric & antisymmetric decomposition

Any 3x3 matrix L splits uniquely into a symmetric part D = ½(L + Lᵀ) and an antisymmetric part W = ½(L - Lᵀ). D is the *strain-rate tensor* — the rate at which fluid elements stretch or shear. W is the *vorticity tensor* — the rate at which fluid elements rotate as rigid bodies. Only D contributes to viscous stress; W does not, because rigid rotation costs no internal friction. For pure shear, both halves are equally non-zero and equally important physically, but only D shows up in sigma.

EXAMPLE

L = [[0, 2, 0], [0, 0, 0], [0, 0, 0]]
D = 0.5 * (L + L.T) = [[0, 1, 0], [1, 0, 0], [0, 0, 0]]
W = 0.5 * (L - L.T) = [[0, 1, 0], [-1, 0, 0], [0, 0, 0]]
# notice W is antisymmetric: W[1,0] = -W[0,1]

22Stress tensor σ = -p I + 2μ D

Same constitutive law as before. The diagonal -p*I represents the isotropic pressure (compressive on the diagonal by convention). The deviatoric viscous part 2μD captures all of the shearing response. Because our pressure is 0 and our D has only off-diagonal entries, sigma is a pure shear-stress tensor.

EXAMPLE

sigma = -0 * I + 2 * mu * D
      = 2 * 1.0 * [[0, 1, 0], [1, 0, 0], [0, 0, 0]]
      = [[0, 2, 0], [2, 0, 0], [0, 0, 0]]

26Run on the midpoint and print

Same evaluation point (0, 0.5, 0) as the NumPy version. autograd gives us the four tensors to full floating-point precision — no finite-difference round-off. The final scalar sigma[0, 1] = 2.0 Pa is the shear stress that a real rheometer measures, and it now agrees *exactly* with the hand calculation tau = mu * U/h = 2.0.

EXAMPLE

L (velocity gradient) =
 tensor([[0., 2., 0.], [0., 0., 0.], [0., 0., 0.]])
D (strain rate) =
 tensor([[0., 1., 0.], [1., 0., 0.], [0., 0., 0.]])
W (vorticity tensor) =
 tensor([[ 0.,  1.,  0.], [-1.,  0.,  0.], [ 0.,  0.,  0.]])
sigma (stress) =
 tensor([[0., 2., 0.], [2., 0., 0.], [0., 0., 0.]])
tau_xy = 2.0 Pa

26 lines without explanation

1import torch
2
3U  = torch.tensor(2.0)
4h  = torch.tensor(1.0)
5mu = torch.tensor(1.0)
6
7def velocity_field(point):
8    """Plane Couette: u = U y/h, v = 0, w = 0."""
9    x, y, z = point[0], point[1], point[2]
10    u = U * y / h
11    v = torch.zeros_like(x)
12    w = torch.zeros_like(x)
13    return torch.stack([u, v, w])
14
15def velocity_gradient(point):
16    """Exact 3x3 Jacobian via reverse-mode autograd."""
17    return torch.autograd.functional.jacobian(velocity_field, point)
18
19def strain_and_stress(point, pressure):
20    L = velocity_gradient(point)
21    D = 0.5 * (L + L.T)                     # symmetric part
22    W = 0.5 * (L - L.T)                     # antisymmetric (vorticity) part
23    sigma = -pressure * torch.eye(3) + 2 * mu * D
24    return L, D, W, sigma
25
26point = torch.tensor([0.0, 0.5, 0.0], requires_grad=True)
27L, D, W, sigma = strain_and_stress(point, pressure=torch.tensor(0.0))
28
29print("L (velocity gradient) =\n", L)
30print("D (strain rate)        =\n", D)
31print("W (vorticity tensor)   =\n", W)
32print("sigma (stress)         =\n", sigma)
33print("tau_xy =", sigma[0, 1].item(), "Pa")

Why bother with autograd here? For Couette flow the finite-difference result is already exact (because the field is linear). But the moment you put a real velocity field on the right-hand side — say the output of a neural network or the result of a CFD interpolation — autograd's exactness keeps the strain-rate and stress tensors free of discretisation noise. That same machinery is what makes PINNs (physics-informed neural networks) able to enforce

\nabla \cdot \boldsymbol{\sigma} + \rho\mathbf{g} = \rho\,D\mathbf{v}/Dt

as a loss term.

Beyond Newton: Non-Newtonian Fluids

Many real fluids violate the linear stress–strain-rate relation. The interactive explorer above lets you toggle between three rheology models. Their $\tau$ -vs- $\dot{\gamma}$ curves look like this:

Class	τ vs γ̇	Examples
Newtonian	τ = μ γ̇ (straight line)	Water, air, gasoline
Shear-thinning (pseudoplastic)	n < 1, slope decreases	Blood, ketchup, paint, polymer melts
Shear-thickening (dilatant)	n > 1, slope increases	Corn-starch slurry, wet sand
Bingham plastic	yield stress τ₀ then linear	Toothpaste, mayonnaise, drilling mud
Viscoelastic	stress depends on rate AND history	Polymer solutions, silly putty

Important. The Navier–Stokes equations as written rely on the Newtonian constitutive law. For a non-Newtonian fluid you keep the conservation equations (mass, momentum) but swap

\boldsymbol{\sigma} = -p\,\mathbf{I} + 2\mu\,\mathbf{D}

for whatever model your fluid actually obeys. The rest of the derivation is unchanged.

Summary

Viscosity is microscopic momentum diffusion between fluid layers. One number $\mu$ captures it for Newtonian fluids.
The 1D Newton law $\tau = \mu\,\mathrm{d}u/\mathrm{d}y$ generalises to the Cauchy stress tensor $\boldsymbol{\sigma}$ — a symmetric $3 \times 3$ object that describes force per area on every face orientation.
The stress tensor splits into an isotropic pressure part $-p\,\mathbf{I}$ and a deviatoric viscous part $\boldsymbol{\tau}$ .
The velocity gradient $\mathbf{L} = \nabla \mathbf{v}$ splits into the symmetric strain-rate $\mathbf{D}$ (deformation) and the antisymmetric vorticity $\mathbf{W}$ (rotation). Only $\mathbf{D}$ produces viscous stress.
The Newtonian constitutive law $\boldsymbol{\sigma} = -p\,\mathbf{I} + 2\mu\,\mathbf{D}$ is a linear, isotropic relation between stress and strain rate.
Taking the divergence of $\boldsymbol{\sigma}$ in Cauchy's momentum equation produces the familiar $\mu\,\nabla^{2}\mathbf{v}$ viscous term in Navier–Stokes.
The interactive explorer and the worked Couette-flow example show every tensor collapsing to the same scalar answer $\tau = \mu U/h$ .
Real fluids include shear-thinning, shear-thickening, and yield-stress materials. Swapping the constitutive law is how those models are folded into the same momentum-conservation framework.

In one sentence:

"Viscosity turns the linear shear stress on a wall into a tensorial law for stress everywhere — and that tensor, once you take its divergence, is exactly the viscous term in Navier–Stokes."

Coming next: Section 30.5 looks at boundary conditions in fluid flow — no-slip, free-slip, inlet/outlet — and how they pin down a unique solution of the Navier–Stokes equations.