Bulk viscosity¶
The extra pressure a fluid feels while it is compressed or expanded.
Stokes's bulk viscous term: Pi = -zeta div(v), zeta the bulk viscosity and div(v) the velocity's divergence, the rate the fluid's volume grows per volume. A compressing fluid (div(v) < 0) pushes back harder than its equilibrium pressure, an expanding one less, so a shock or a blast's energy goes to heat; zeta = 0 is Stokes's hypothesis, a monatomic gas's. Source: Landau and Lifshitz, Fluid Mechanics, 2nd ed., section 15 (the stress tensor's bulk term); Batchelor, An Introduction to Fluid Dynamics, section 3.4. Checked against a reference in numpy, which has its own tests.
- Source
- Landau and Lifshitz, Fluid Mechanics, 2nd ed., section 15; Batchelor, An Introduction to Fluid Dynamics, section 3.4.
- Proof
- reference, invariant, limit.
bulk_viscosity.trait
namespace lusine.laws;
float<Pa> bulk_viscosity(float<Pa*s> viscosity, float<s^-1> divergence)
{
return -(viscosity * divergence);
}
The extra pressure a fluid feels while it is compressed or expanded, Pi = -zeta div(v).
import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
import numpy as np
from lusine import trait_run
F32 = np.float32
U = float(np.finfo(F32).eps) / 2 # binary32's unit roundoff, 2^-24
rng = np.random.default_rng(2026)
The kernel¶
The law as the library states it.
print(trait_run.source('laws/bulk_viscosity.trait'))
// lusine.laws.bulk_viscosity: the extra pressure a fluid feels while it is compressed or expanded.
// Stokes's bulk viscous term: Pi = -zeta div(v), zeta the bulk viscosity and div(v) the velocity's divergence, the
// rate the fluid's volume grows per volume. A compressing fluid (div(v) < 0) pushes back harder than its equilibrium
// pressure, an expanding one less, so a shock or a blast's energy goes to heat; zeta = 0 is Stokes's hypothesis, a
// monatomic gas's.
// Source: Landau and Lifshitz, Fluid Mechanics, 2nd ed., section 15 (the stress tensor's bulk term); Batchelor, An
// Introduction to Fluid Dynamics, section 3.4.
// Checked against a reference in numpy, which has its own tests.
namespace lusine.laws;
float<Pa> bulk_viscosity(float<Pa*s> viscosity, float<s^-1> divergence)
{
return -(viscosity * divergence);
}
The program¶
Each section runs the kernel over a lattice of items, one a case: the host sets every parameter, and main() binds the kernel over the items.
print(trait_run.law_program('laws/bulk_viscosity.trait', 'bulk_viscosity'))
from lusine.kernels.laws.bulk_viscosity import bulk_viscosity as reference
def law(**kw):
return trait_run.law('laws/bulk_viscosity.trait', 'bulk_viscosity', **kw)
namespace probe;
lusine.version version = { { number = 1 } };
struct item {
host float<Pa*s> viscosity;
host float<s^-1> divergence;
float<Pa> result_0;
};
lattice<item, 1> items(spacing = 1 m, boundary = wall);
void main() { items.result_0 = lusine.laws.bulk_viscosity(items.viscosity, items.divergence); }
Reference: compressions and expansions¶
Bulk viscosities to 10 Pa s and divergences either way: the kernel within a rounding of the reference.
z, dv = rng.uniform(0.0, 10.0, 400), rng.uniform(-50.0, 50.0, 400)
got, want = law(viscosity=z, divergence=dv), reference(z, dv)
err = np.max(np.abs(got - want) / np.maximum(np.abs(want), 1e-30))
assert err <= 4 * U, err / U
print(f'within {err / U:.1f} x 2^-24 over {len(z)} flows')
within 2.2 x 2^-24 over 400 flows
Invariant: always a heating¶
The work done against it, -Pi div(v) = zeta div(v)^2, is never negative, so it never cools the fluid.
pi = law(viscosity=z, divergence=dv)
assert np.all(-pi * dv >= 0.0)
print('never a cooling')
never a cooling
Limit: Stokes's hypothesis¶
With no bulk viscosity, a monatomic gas's, there is no extra pressure however fast the fluid is compressed.
pi = law(viscosity=np.zeros(3), divergence=np.array([-1e3, 0.0, 1e3]))
assert np.all(pi == 0.0)
print('zeta = 0: none')
zeta = 0: none