physicskit.fluids #

Returns:

dict – "k", time-averaged "E", the sampled "t" and "energy" series, and the final "omega".

spectrum()[source]#

Shell-summed energy spectrum \(E(k)\) of the current state, with \(\sum_k E(k)\) the energy.

Return type:

Returns:

k (ndarray of float, shape (n//2,)) – Integer shell wavenumbers.
E (ndarray of float, shape (n//2,))

vorticity()[source]#

Real-space vorticity field, shape (n, n).

Return type:: NDArray[double]

class physicskit.fluids.ForcedTurbulence3D(n=64, kf=2.5, epsilon=1.0, hyperviscosity_order=2, damping_at_kmax=20.0, seed=None)[source]#

Bases: object

Forced, statistically steady 3D turbulence on a \((2\pi)^3\)-periodic cube.

The incompressible Navier-Stokes equations are integrated in rotational form,

\[\partial_t\hat{\mathbf u} = \mathsf P\,\widehat{\mathbf u\times\boldsymbol\omega} - \nu_p k^{2p}\hat{\mathbf u} + \hat{\mathbf f},\]

where \(\mathsf P\) projects onto divergence-free fields, with 2/3-rule dealiasing, an exact integrating factor for the (hyper)viscosity and a second-order Runge-Kutta step. The forcing \(\hat{\mathbf f} = \varepsilon\,\hat{\mathbf u}/(2E_f)\) on the shell \(0 < k < k_f\), with \(E_f\) that shell’s energy, injects energy at exactly the rate \(\varepsilon\) (Lundgren 2003; Alvelius 1999 for random variants). In 3D the energy cascades to small scales, and in the inertial range Kolmogorov’s 1941 law \(E(k) = C_K\varepsilon^{2/3}k^{-5/3}\), \(C_K \approx 1.5\), holds. On a \(64^3\) grid that range spans only \(3 \lesssim k \lesssim 10\).

Parameters:

n (int)
kf (float) – Forcing cutoff.
epsilon (float)
hyperviscosity_order (int) – Power \(p\); 1 is ordinary viscosity.
damping_at_kmax (float) – \(\nu_p k_{\max}^{2p}\), the damping rate at the dealiasing cutoff.
seed (int | None) – Seed for the random initial field.

Examples

>>> flow = ForcedTurbulence3D(n=16, seed=0)
>>> out = flow.run(t_max=0.2, dt=0.01)
>>> bool(flow.max_divergence() < 1e-10)
True

dissipation_rate()[source]#

\(\sum_k 2\nu_p k^{2p}E(\mathbf k)\), the rate the (hyper)viscosity removes energy.

Return type:: float

kinetic_energy()[source]#

\(\tfrac12\langle|\mathbf u|^2\rangle\).

Return type:: float

max_divergence()[source]#

Largest \(|\nabla\cdot\mathbf u|\) on the grid (zero up to roundoff).

Return type:: float

run(t_max=10.0, dt=0.01, t_average=5.0, sample_every=10)[source]#

Advance and time-average the spectrum over the last t_average.

Return type:

Returns:

dict – "k", time-averaged "E", sampled "t", "energy" and "dissipation".

Parameters:

t_max (float)
dt (float)
t_average (float)
sample_every (int)

spectrum()[source]#

Shell-summed energy spectrum, with \(\sum_k E(k)\) the kinetic energy.

Return type:: tuple[NDArray[double], NDArray[double]]

exception physicskit.fluids.InvalidParameterError[source]#

Bases: FluidskitError, ValueError

A system or tool was constructed or called with a physically invalid parameter.

Examples include a non-positive Reynolds number or viscosity, a grid size incompatible with the doubly-periodic pseudo-spectral method’s FFT assumptions, a cylinder radius too large to fit its sampling grid, or an explicit time step that violates a CFL/Courant stability bound.

class physicskit.fluids.LatticeBoltzmannD2Q9(nx, ny, tau, body_force=(0.0, 0.0), solid=None, rho0=1.0, velocity=None)[source]#

Bases: object

A D2Q9 BGK lattice Boltzmann fluid on a periodic box with optional solid obstacles.

Both axes wrap periodically; channels, cylinders, and cavities are built by marking solid nodes, whose faces toward the fluid become no-slip walls by halfway bounce-back (the wall sits half a lattice spacing outside the last fluid node). A uniform body force (acceleration) body_force drives the flow, e.g. as the pressure gradient of a periodic channel.

Parameters:

nx (int) – Lattice size.
ny (int) – Lattice size.
tau (float) – BGK relaxation time, \(\tau > 1/2\); the viscosity is \(\nu = (\tau - 1/2)/3\) (lbm_viscosity()).
body_force (tuple of float, default=(0.0, 0.0)) – Acceleration \(\mathbf g\) acting on the fluid (force per unit mass), in lattice units.
solid (array_like of bool, shape (nx, ny), optional) – Solid (wall) nodes.
rho0 (float, default=1.0) – Initial uniform density.
velocity (tuple of array_like, optional) – Initial velocity field (ux, uy); zero by default. The populations start at the matching equilibrium.

Variables:

f (ndarray, shape (9, nx, ny)) – Particle distributions.
time (int) – Number of steps taken.

Examples

Flow between two plates driven by a body force (plane Poiseuille flow):

>>> solid = np.zeros((1, 22), dtype=bool)
>>> solid[:, [0, -1]] = True
>>> lb = LatticeBoltzmannD2Q9(1, 22, tau=0.8, body_force=(1e-6, 0.0), solid=solid)
>>> lb.step(6000)
>>> H, nu = 20.0, lb.viscosity
>>> round(float(lb.velocity[0][0, 11] / (1e-6 * H**2 / (8 * nu))), 2)
1.0

property density: NDArray[float64]#: Density \(\rho = \sum_i f_i\) (zero-padded meaning on solid nodes).

step(n_steps=1)[source]#

Advance n_steps collide-and-stream updates.

Parameters:: n_steps (int)
Return type:: None

property velocity: tuple[NDArray[float64], NDArray[float64]]#

Velocity \(\mathbf u = (\sum_i f_i\mathbf e_i + \tfrac12\rho\mathbf g)/\rho\), zero on solid nodes.

Returns:: ux, uy (ndarray, shape (nx, ny))

vorticity()[source]#

Vorticity \(\partial_x u_y - \partial_y u_x\) by periodic central differences.

Return type:: NDArray[double]
Returns:: ndarray, shape (nx, ny)

class physicskit.fluids.NavierStokes2D(n=64, length=6.283185307179586, nu=0.01)[source]#

Bases: object

A doubly periodic 2D incompressible Navier-Stokes solver, in vorticity-streamfunction form.

Builds the pseudo-spectral grid once at construction and exposes rhs(), step(), and simulate() against it, so that repeated calls (e.g. one per animation frame) do not re-derive the wavenumber grid on every call the way the module-level simulate_vorticity_streamfunction() convenience function must.

Parameters:

n (int) – Number of grid points along each axis.
length (float) – Physical domain size (the domain is [0, length) x [0, length)).
nu (float) – Kinematic viscosity. Must be positive; a truly inviscid flow is ill-posed on this grid (nothing controls the cascade of enstrophy to the grid scale), so nu=0 is rejected.

Variables:

n (int) – Grid resolution.
length (float) – Domain size.
nu (float) – Kinematic viscosity.
Y (X,) – Real-space coordinate grids.
K2 (KX, KY,) – Wavenumber grids.

Raises:

InvalidParameterError – If nu is not positive (see physicskit.fluids.core.grid.spectral_grid() for the conditions on n and length).

Examples

>>> import numpy as np
>>> solver = NavierStokes2D(n=48, length=2 * np.pi, nu=0.05)
>>> omega0 = np.sin(solver.X) * np.sin(solver.Y)
>>> result = solver.simulate(omega0, dt=0.01, steps=20)
>>> result["omega"].shape
(48, 48)
>>> bool(np.max(np.abs(result["omega"])) < np.max(np.abs(omega0)))
True

rhs(omega)[source]#

Evaluate the vorticity-transport right-hand side at omega.

Parameters:: omega (NDArray[double]) – Vorticity field.
Return type:: NDArray[double]
Returns:: ndarray of float, shape (n, n) – \(-(\mathbf{u}\cdot\nabla)\omega + \nu\nabla^2\omega\).

simulate(omega0, dt, steps)[source]#

Time-step an initial vorticity field forward by steps RK4 steps.

Parameters:

omega0 (NDArray[double]) – Initial vorticity field.
dt (float) – Time step.
steps (int) – Number of RK4 steps to advance.

Return type:

Returns:

dict – {"omega": final vorticity, "psi": final streamfunction, "u": ..., "v": ...}.

step(omega, dt)[source]#

Advance the vorticity field by one RK4 step of size dt.

Parameters:

omega (NDArray[double]) – Current vorticity field.
dt (float) – Step size.

Return type:

Returns:

ndarray of float, shape (n, n) – Vorticity field advanced by one step.

velocity(omega)[source]#

Recover the velocity field associated with a vorticity field.

Parameters:: omega (NDArray[double]) – Vorticity field.
Return type:: tuple[NDArray[double], NDArray[double]]
Returns:: u, v (ndarray of float, shape (n, n)) – Velocity components.

class physicskit.fluids.PlaneElasticityResult(displacement, strain, stress, centroids)[source]#

Bases: NamedTuple

Solution of plane_elasticity_solve().

Variables:

displacement (ndarray, shape (n_nodes, 2)) – Nodal displacements \((u_x, u_y)\).
strain (ndarray, shape (n_triangles, 3)) – Per-element \((\varepsilon_{xx}, \varepsilon_{yy}, \gamma_{xy})\).
stress (ndarray, shape (n_triangles, 3)) – Per-element \((\sigma_{xx}, \sigma_{yy}, \sigma_{xy})\).
centroids (ndarray, shape (n_triangles, 2)) – Element centroids, where the (constant) strain and stress apply best.

Parameters:

displacement (ndarray)
strain (ndarray)
stress (ndarray)
centroids (ndarray)

centroids: ndarray#: Alias for field number 3

displacement: ndarray#: Alias for field number 0

strain: ndarray#: Alias for field number 1

stress: ndarray#: Alias for field number 2

class physicskit.fluids.PointVortexSystem(positions, circulations, core=1e-06)[source]#

Bases: object

A system of interacting 2D point vortices, advected by their mutual Biot-Savart field.

Parameters:

positions (ArrayLike) – Initial vortex positions (x, y).
circulations (ArrayLike) – Vortex circulations \(\Gamma_i\); fixed for the system’s lifetime (point vortices carry their circulation with them unchanged).
core (float) – Regularization core radius (see point_vortex_velocities()).

Variables:

positions (ndarray of float, shape (n, 2)) – Current vortex positions.
circulations (ndarray of float, shape (n,)) – Vortex circulations.
core (float) – Regularization core radius.
n_vortices (int) – Number of vortices, n.

Raises:

InvalidParameterError – If positions and circulations have inconsistent lengths.

Examples

>>> import numpy as np
>>> system = PointVortexSystem(positions=[[1.0, 0.0], [-1.0, 0.0]], circulations=[1.0, 1.0])
>>> t, trajectory = system.trajectory(dt=0.01, n_steps=10)
>>> trajectory.shape
(11, 2, 2)

trajectory(dt=0.01, n_steps=1000)[source]#

Integrate the vortex trajectories forward with RK4.

Parameters:

dt (float) – Time step.
n_steps (int) – Number of steps to advance.

Return type:

Returns:

times (ndarray of float, shape (n_steps + 1,))
positions (ndarray of float, shape (n_steps + 1, n, 2)) – Vortex positions at each step, including the current positions.

velocities()[source]#

Velocity currently induced on each vortex.

Return type:: NDArray[double]
Returns:: ndarray of float, shape (n, 2)

class physicskit.fluids.PotentialFlow(U_inf=1.0, alpha=0.0)[source]#

Bases: object

A superposition of elementary potential-flow solutions.

Because Laplace’s equation is linear, the complex potential of any combination of uniform flow, sources/sinks, doublets, and point vortices is simply the sum of their individual complex potentials – exactly what this class accumulates. Call add_source(), add_doublet(), or add_vortex() to build up a flow on top of the base uniform stream (U_inf, alpha), then velocity(), streamfunction(), or pressure_coefficient() to evaluate it on a grid.

Parameters:

U_inf (float) – Free-stream speed used both as the base uniform flow and as the reference speed for pressure_coefficient().
alpha (float) – Free-stream angle of attack (radians).

Variables:

alpha (U_inf,) – As above.

Examples

>>> flow = PotentialFlow(U_inf=1.0)
>>> flow.add_doublet(strength=2 * np.pi * 1.0**2)  # uniform flow + doublet = cylinder
>>> X, Y = np.meshgrid(np.linspace(-3, 3, 5), np.linspace(-3, 3, 5))
>>> u, v = flow.velocity(X, Y)
>>> u.shape
(5, 5)

add_doublet(strength, x0=0.0, y0=0.0, alpha=0.0)[source]#

Add a doublet.

Parameters:

strength (float) – Doublet strength.
x0 (float) – Doublet location.
y0 (float) – Doublet location.
alpha (float) – Orientation angle (radians).

Return type:

add_source(strength, x0=0.0, y0=0.0)[source]#

Add a point source (or sink, for negative strength).

Parameters:

strength (float) – Source strength; negative for a sink.
x0 (float) – Source location.
y0 (float) – Source location.

Return type:

add_vortex(circulation, x0=0.0, y0=0.0)[source]#

Add a point vortex.

Parameters:

circulation (float) – Circulation \(\Gamma\); positive is counterclockwise.
x0 (float) – Vortex location.
y0 (float) – Vortex location.

Return type:

complex_potential(X, Y)[source]#

Evaluate the total complex potential \(W(z)=\phi+i\psi\).

Parameters:

X (ArrayLike) – Field points.
Y (ArrayLike) – Field points.

Return type:

Returns:

ndarray of complex – Complex potential at each point.

pressure_coefficient(X, Y)[source]#

Evaluate the Bernoulli pressure coefficient \(C_p = 1 - (u^2+v^2)/U_\infty^2\).

Parameters:

X (ArrayLike) – Field points.
Y (ArrayLike) – Field points.

Return type:

Returns:

ndarray of float – Pressure coefficient at each point.

streamfunction(X, Y)[source]#

Evaluate the streamfunction \(\psi = \mathrm{Im}(W)\).

Parameters:

X (ArrayLike) – Field points.
Y (ArrayLike) – Field points.

Return type:

Returns:

ndarray of float – Streamfunction values.

velocity(X, Y, delta=1e-06)[source]#

Evaluate the velocity field \(u - iv = dW/dz\).

Differentiates the complex potential with a central difference in X (equivalent, by the Cauchy-Riemann equations, to differentiating in either direction), which sidesteps deriving and summing a separate closed-form velocity expression for every element type.

Parameters:

X (ArrayLike) – Field points.
Y (ArrayLike) – Field points.
delta (float) – Finite-difference step used to differentiate the complex potential.

Return type:

Returns:

u, v (ndarray of float) – Velocity components.

class physicskit.fluids.RayleighBenard2D(Ra, Pr=1.0, aspect=np.float64(2.8284271247461903), nx=32, nz=32)[source]#

Bases: object

Nonlinear 2D Rayleigh-Bénard convection between stress-free, isothermal walls.

The layer \(0 < z < 1\) is extended to \(-1 < z < 1\) with \(\theta\), \(w\) and \(\omega\) odd in \(z\), which satisfies \(w = \theta = \partial_z u = 0\) at both walls and makes the problem periodic in \(z\) with period 2. It is then solved pseudo-spectrally like NavierStokes2D, with 2/3-rule dealiasing, an exact integrating factor for diffusion, RK4 for the remaining terms, and the odd symmetry re-imposed at every sample.

Parameters:

Ra (float) – Rayleigh number.
Pr (float) – Prandtl number.
aspect (float) – Box width over depth; the default fits one critical wavelength \(2\pi/a_c\) of the stress-free problem.
nx (int) – Grid points in \(x\) and over the doubled depth.
nz (int) – Grid points in \(x\) and over the doubled depth.

Variables:

omega (theta,) – Temperature perturbation and vorticity on the doubled domain.
z (x,) – Grid coordinates; \(z \in [-1, 1)\).

Examples

>>> rb = RayleighBenard2D(Ra=2000.0, nx=16, nz=16)
>>> rb.seed_mode(amplitude=1e-3)
>>> out = rb.run(t_max=0.2, dt=1e-3)
>>> bool(out["kinetic_energy"][-1] > out["kinetic_energy"][0])  # supercritical: grows
True

kinetic_energy()[source]#

Mean kinetic energy \(\tfrac12\langle u^2 + w^2\rangle\).

Return type:: float

nusselt()[source]#

Nusselt number \(Nu = 1 + \langle w\theta\rangle\) (conductive heat flux is 1).

Return type:: float

run(t_max, dt=0.001, sample_every=10)[source]#

Advance the convection with RK4.

Parameters:

t_max (float) – Time to integrate, in thermal diffusion times \(d^2/\kappa\).
dt (float)
sample_every (int) – Steps between samples.

Return type:

Returns:

dict – Sampled "t", "kinetic_energy", "nusselt" and the amplitude "theta_mode" of the seeded \(\sin(\pi z) \cos(2\pi x/\Gamma)\) mode.

seed_mode(amplitude=0.001, noise=0.0, seed=None)[source]#

Set \(\theta = A\sin(\pi z)\cos(2\pi x/\Gamma)\) plus optional odd noise, at rest.

Parameters:

amplitude (float)
noise (float) – Amplitude of additional random temperature perturbations.
seed (int | None)

Return type:

velocity(omega=None)[source]#

Velocity \((u, w)\) from the vorticity, via \(\nabla^2\psi = -\omega\), \(u = \partial_z\psi\), \(w = -\partial_x\psi\).

Return type:: tuple[NDArray[double], NDArray[double]]
Parameters:: omega (NDArray[float64] | None)

class physicskit.fluids.RayleighBenardWalls2D(Ra, Pr=1.0, aspect=None, boundaries='rigid', nx=32, nz=24)[source]#

Bases: object

2D Rayleigh-Bénard convection between rigid or stress-free walls, Fourier-Chebyshev.

The streamfunction \(\psi\) (\(u = \partial_z\psi\), \(w = -\partial_x\psi\), vorticity \(\omega = -\nabla^2\psi\)) and the temperature perturbation \(\theta\) are expanded in Fourier modes \(e^{ikx}\) across a periodic box and collocated on Chebyshev points in \(0 \le z \le 1\), where

\[\partial_t\nabla^2\psi = Pr\,\nabla^4\psi - Pr\,Ra\,\partial_x\theta + (\mathbf u\cdot\nabla)\omega, \qquad \partial_t\theta = \nabla^2\theta + w - (\mathbf u\cdot\nabla)\theta.\]

Diffusion is treated implicitly (Crank-Nicolson) and the rest explicitly (second-order Adams-Bashforth). The walls are isothermal, \(\theta = 0\), impermeable, \(\psi = 0\), and either no-slip, \(\partial_z\psi = 0\) ("rigid"), or stress-free, \(\partial_z^2\psi = 0\) ("free"), imposed by replacing the collocation equations at the two points next to each wall (Canuto et al., Spectral Methods, 2006, sec. 3.7).

Parameters:

Ra (float)
Pr (float)
aspect (float | None) – Box width; defaults to one critical wavelength \(2\pi/a_c\) for the chosen walls.
boundaries (str)
nx (int) – Fourier points across the box.
nz (int) – Chebyshev polynomial degree (nz + 1 points).

Examples

>>> rb = RayleighBenardWalls2D(Ra=2500.0, nx=16, nz=16)
>>> rb.seed_mode(1e-4)
>>> out = rb.run(t_max=0.3, dt=2e-3)
>>> bool(out["theta_mode"][-1] > out["theta_mode"][0])  # above 1708: grows
True

fields()[source]#

Real-space "u", "w", "theta" and "omega", shape (nx, nz + 1).

Return type:: dict[str, NDArray[double]]

nusselt()[source]#

Nusselt number \(1 + \langle w\theta\rangle\), averaged over the layer (Clenshaw-Curtis).

Return type:: float

run(t_max, dt=0.002, sample_every=10)[source]#

Advance with CN-AB2.

Return type:

Returns:

dict – Sampled "t", "nusselt" and "theta_mode", the amplitude of the first Fourier mode of \(\theta\) at mid-depth.

Parameters:

t_max (float)
dt (float)
sample_every (int)

seed_mode(amplitude=0.001)[source]#

Set \(\theta = A\sin(\pi z)\cos(2\pi x/\Gamma)\) with the fluid at rest.

Return type:: None
Parameters:: amplitude (float)

physicskit.fluids.acoustic_wave_1d(p0, length, c, t_max, dt=None, ends=('closed', 'closed'), rho0=1.0, u0=None, n_frames=100, probe=None)[source]#

Sound in a 1D pipe: staggered-grid leapfrog for the linearized Euler equations.

The pipe \(0 \le x \le L\) is split into n = len(p0) cells; pressure lives at cell centers \(x_i = (i + \tfrac12)\Delta x\) and velocity on the n + 1 faces. Each step updates

\[u_{i}^{n+1/2} = u_i^{n-1/2} - \frac{\Delta t}{\rho_0\Delta x}(p_i^n - p_{i-1}^n), \qquad p_i^{n+1} = p_i^n - \frac{\rho_0 c^2\Delta t}{\Delta x}(u_{i+1}^{n+1/2} - u_i^{n+1/2}).\]

Parameters:

p0 (ArrayLike) – Initial acoustic pressure at the cell centers.
length (float) – Pipe length.
c (float) – Sound speed.
t_max (float) – Final time, rounded to a whole number of steps.
dt (float | None) – Time step; defaults to 0.5 * dx / c. Must satisfy c * dt <= dx.
ends (tuple of {"closed", "open"}, default=("closed", "closed")) – Rigid (u = 0) or pressure-release (p = 0) at each end.
rho0 (float) – Ambient density.
u0 (ArrayLike | None) – Initial face velocities (at \(t = -\Delta t/2\)); zero by default.
n_frames (int) – Approximate number of recorded pressure frames.
probe (int | None) – Cell whose pressure is recorded at every step (for spectra); defaults to cell 0.

Returns:

times (ndarray, shape (k + 1,)) – Frame times.
frames (ndarray, shape (k + 1, n)) – Pressure at each frame.
trace (ndarray, shape (n_steps + 1,)) – Pressure at probe after every step, sampled every dt.
dt (float) – The time step used.

Examples

The fundamental of a closed-closed pipe returns to itself after one period \(2L/c\):

>>> x = (np.arange(200) + 0.5) / 200
>>> t, P, trace, dt = acoustic_wave_1d(np.cos(np.pi * x), 1.0, 1.0, t_max=2.0, dt=0.0025, n_frames=1)
>>> bool(np.max(np.abs(P[-1] - P[0])) < 1e-3)
True

physicskit.fluids.acoustic_wave_2d(p0, size, c, t_max, dt=None, rho0=1.0, n_frames=100, probe=(0, 0))[source]#

Sound in a rigid-walled rectangular room: 2D staggered-grid leapfrog.

The 2D version of acoustic_wave_1d(), with pressure at cell centers, \(u_x\) on the x-faces and \(u_y\) on the y-faces, and zero normal velocity on all four walls. Its standing modes are those of rectangular_room_mode_frequencies().

Parameters:

p0 (ArrayLike) – Initial pressure at cell centers.
size (tuple of float) – Room dimensions (Lx, Ly).
c (float) – Sound speed.
t_max (float) – Final time.
dt (float | None) – Time step; defaults to half the stability limit 1 / (c * sqrt(1/dx^2 + 1/dy^2)).
rho0 (float) – Ambient density.
n_frames (int) – Approximate number of recorded frames.
probe (tuple of int, default=(0, 0)) – Cell whose pressure is recorded at every step.

Returns:

times (ndarray)
frames (ndarray, shape (k + 1, nx, ny))
trace (ndarray, shape (n_steps + 1,))
dt (float)

Examples

>>> p0 = np.zeros((40, 30))
>>> p0[20, 15] = 1.0
>>> t, P, trace, dt = acoustic_wave_2d(p0, (4.0, 3.0), 340.0, t_max=0.01, n_frames=5)
>>> P.shape[1:], bool(abs(P[-1].sum() - 1.0) < 1e-9)  # rigid walls conserve the mean pressure
((40, 30), True)

physicskit.fluids.annulus_mesh(inner_radius, outer_radius, n_r, n_theta, theta_max=1.5707963267948966)[source]#

Structured triangulation of an annular sector \(a \le r \le b\), \(0 \le \theta \le \theta_{\max}\).

A quarter annulus (the default) with symmetry conditions on its two straight edges models a full pressurized cylinder at a quarter of the cost.

Parameters:

inner_radius (float) – Radii \(a < b\).
outer_radius (float) – Radii \(a < b\).
n_r (int) – Number of cells radially and around.
n_theta (int) – Number of cells radially and around.
theta_max (float) – Opening angle of the sector.

Returns:

nodes (ndarray, shape ((n_r + 1) * (n_theta + 1), 2))
triangles (ndarray of int, shape (2 * n_r * n_theta, 3))
boundaries (dict) – "inner" and "outer": arrays of shape (n_theta, 2) of boundary edges, each directed with the solid on its right (the convention of pressure_load()); "theta0" and "theta_max": node indices on the two straight edges.

Examples

>>> nodes, tris, bnd = annulus_mesh(1.0, 2.0, 4, 8)
>>> len(nodes), len(tris), bnd["inner"].shape
(45, 64, (8, 2))

physicskit.fluids.blasius_boundary_layer_thickness(x, U_inf, nu, eta_99=4.91)[source]#

Blasius laminar boundary-layer thickness \(\delta_{99}(x)\).

The height at which \(u/U_\infty\) first reaches 0.99 occurs, in similarity coordinates, at \(\eta \approx 4.91\) (a property of the Blasius profile from blasius_solve(), not re-derived here), so

\[\delta_{99}(x) = \eta_{99}\sqrt{\frac{\nu x}{U_\infty}} \propto \sqrt{x}.\]

The boundary layer’s parabolic-looking growth with downstream distance – much slower than the layer’s own streamwise extent – is exactly the “thin layer” assumption Ludwig Prandtl’s 1904 boundary-layer theory used to simplify Navier-Stokes into the boundary-layer equations Blasius then solved.

Parameters:

x (ArrayLike) – Downstream distance from the leading edge; must be positive.
U_inf (float) – Free-stream speed.
nu (float) – Kinematic viscosity.
eta_99 (float) – Similarity coordinate at which \(f'(\eta)=0.99\).

Return type:

Returns:

ndarray of float – Boundary-layer thickness \(\delta_{99}(x)\).

Raises:

InvalidParameterError – If U_inf or nu is not positive, or any entry of x is not positive.

Examples

>>> round(float(blasius_boundary_layer_thickness(x=1.0, U_inf=1.0, nu=1e-4)), 4)
0.0491

physicskit.fluids.blasius_skin_friction_coefficient(reynolds_x, fpp0=0.33206)[source]#

Local skin-friction coefficient of the Blasius boundary layer.

The wall shear stress \(\tau_w = \mu U_\infty f''(0)\sqrt{U_\infty/(\nu x)}\) gives a local skin-friction coefficient

\[c_f(x) = \frac{\tau_w}{\tfrac{1}{2}\rho U_\infty^2} = \frac{2 f''(0)}{\sqrt{Re_x}}, \qquad Re_x = \frac{U_\infty x}{\nu},\]

the classic \(c_f \propto Re_x^{-1/2}\) scaling of a laminar boundary layer.

Parameters:

reynolds_x (ArrayLike) – Local Reynolds number \(Re_x = U_\infty x/\nu\); must be positive.
fpp0 (float) – The Blasius constant \(f''(0)\), as returned by blasius_solve().

Return type:

Returns:

ndarray of float – Local skin-friction coefficient \(c_f(x)\).

Raises:

InvalidParameterError – If any entry of reynolds_x is not positive.

Examples

>>> round(float(blasius_skin_friction_coefficient(reynolds_x=1e4)), 5)
0.00664

physicskit.fluids.blasius_solve(eta_max=10.0, n_points=2001, tol=1e-10, max_iter=100)[source]#

Solve the Blasius laminar boundary-layer equation by shooting.

The Blasius similarity reduction of the steady, 2D, zero-pressure-gradient boundary-layer equations collapses them to the single third-order ODE

\[f''' + \tfrac{1}{2} f f'' = 0, \qquad f(0) = f'(0) = 0, \qquad f'(\infty) = 1,\]

for the dimensionless streamfunction \(f(\eta)\) of the similarity variable \(\eta = y\sqrt{U_\infty/(\nu x)}\); the streamwise velocity is \(u/U_\infty = f'(\eta)\). Since \(f''(0)\) (rather than the far-field condition) is the free parameter needed to start an initial-value integration, this is a two-point boundary value problem solved here by shooting: bisect on \(f''(0)\) until the resulting _blasius_integrate() solution satisfies \(f'(\eta_{max}) \approx 1\).

Parameters:

eta_max (float) – Similarity coordinate treated as “infinity”; large enough that \(f'\) has converged to 1 to within tol.
n_points (int) – Number of output points between eta=0 and eta_max.
tol (float) – Bisection tolerance on \(f'(\eta_{max}) - 1\).
max_iter (int) – Maximum number of bisection iterations.

Return type:

Returns:

dict – {"eta": ..., "f": ..., "fp": ..., "fpp": ...}, each an array of length n_points. fpp[0] is the classic Blasius constant, \(f''(0) \approx 0.332\).

Examples

>>> result = blasius_solve()
>>> round(float(result["fpp"][0]), 3)
0.332
>>> round(float(result["fp"][-1]), 4)
1.0

physicskit.fluids.couette_flow_velocity(y, U_wall, h)[source]#

Velocity profile of steady plane Couette flow.

The flow between two infinite parallel plates, the lower one at rest and the upper one (at y = h) moving at speed U_wall, with no imposed pressure gradient. The steady Navier-Stokes equations reduce to \(\mu\, d^2u/dy^2 = 0\), whose solution subject to no-slip at both walls is the linear profile below – shear stress \(\tau = \mu\,U_{wall}/h\) is exactly uniform across the gap.

Parameters:

y (ArrayLike) – Height above the lower (stationary) plate, 0 <= y <= h.
U_wall (float) – Speed of the upper, moving plate.
h (float) – Gap between the plates.

Return type:

Returns:

ndarray of float – Velocity \(u(y) = U_{wall}\,y/h\).

Raises:

InvalidParameterError – If h is not positive.

Examples

>>> couette_flow_velocity(y=[0.0, 0.5, 1.0], U_wall=2.0, h=1.0)
array([0., 1., 2.])

physicskit.fluids.dam_break_exact(x, t, h_left, h_right=0.0, g=9.81)[source]#

Exact solution of the dam break at \(x = 0\), released at \(t = 0\).

A rarefaction fan \(h = (2c_L - x/t)^2/(9g)\), \(u = \tfrac23(c_L + x/t)\) connects the reservoir to a uniform state \((h_m, u_m)\), which is joined to the tailwater by a bore of speed \(s = h_m u_m/(h_m - h_R)\). \(h_m\) solves

\[2(c_L - c_m) = (h_m - h_R)\sqrt{\frac{g(h_m + h_R)}{2h_mh_R}}\]

(Stoker, Water Waves, 1957). For a dry bed, \(h_R = 0\), the fan reaches the wet front \(x = 2c_Lt\) (Ritter 1892).

Parameters:

x (NDArray[double]) – Positions.
t (float) – Time after the release, positive.
h_left (float) – Reservoir depth.
h_right (float) – Tailwater depth, 0 <= h_right < h_left.
g (float)

Return type:

Returns:

h, u (ndarray of float)

Examples

>>> h, u = dam_break_exact(np.array([0.0]), t=1.0, h_left=1.0)
>>> round(float(h[0]), 6)  # Ritter: h = 4 h_L / 9 at the dam site
0.444444

physicskit.fluids.doublet_potential(X, Y, strength, x0=0.0, y0=0.0, alpha=0.0)[source]#

Complex potential of a doublet (a source-sink pair in the zero-separation limit).

Parameters:

X (ArrayLike) – Field points.
Y (ArrayLike) – Field points.
strength (float) – Doublet strength \(\kappa\) (the limit of source_strength * separation as the separation shrinks to zero).
x0 (float) – Doublet location.
y0 (float) – Doublet location.
alpha (float) – Orientation angle (radians) from sink-to-source.

Return type:

Returns:

ndarray of complex – Complex potential \(W(z) = \frac{\kappa e^{-i\alpha}}{2\pi (z-z_0)}\).

physicskit.fluids.elasticity_matrix(E, nu, plane='stress')[source]#

Constitutive matrix \(D\) with \((\sigma_{xx}, \sigma_{yy}, \sigma_{xy}) = D(\varepsilon_{xx}, \varepsilon_{yy}, \gamma_{xy})\).

Plane stress:

\[\begin{split}D = \frac{E}{1-\nu^2}\begin{pmatrix}1&\nu&0\\ \nu&1&0\\ 0&0&\tfrac{1-\nu}{2}\end{pmatrix}.\end{split}\]

Plane strain:

\[\begin{split}D = \frac{E}{(1+\nu)(1-2\nu)}\begin{pmatrix}1-\nu&\nu&0\\ \nu&1-\nu&0\\ 0&0&\tfrac{1-2\nu}{2}\end{pmatrix}.\end{split}\]

Engineering shear strain \(\gamma_{xy} = 2\varepsilon_{xy}\) (Zienkiewicz and Taylor, The Finite Element Method, 5th ed., vol. 1, §4.2).

Parameters:

E (float) – Young’s modulus.
nu (float) – Poisson’s ratio.
plane (str)

Return type:

Returns:

ndarray, shape (3, 3)

Examples

>>> elasticity_matrix(1.0, 0.0).tolist()
[[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 0.5]]

physicskit.fluids.energy_spectrum(u, v, length)[source]#

Isotropic (azimuthally averaged) kinetic energy spectrum of a 2D velocity field.

Computes \(\hat{u}(\mathbf{k})\) and \(\hat{v}(\mathbf{k})\) by FFT, forms the kinetic energy density \(\tfrac{1}{2}(|\hat{u}|^2+|\hat{v}|^2)\) at each wavevector, and sums it over the discrete annulus of wavevectors with \(k - \tfrac12 \le |\mathbf{k}| < k + \tfrac12\) (in grid-index units) to give \(E(k)\), normalized so that \(\sum_k E(k) \approx\) the domain-averaged kinetic energy per unit mass (Parseval’s theorem).

Parameters:

u (NDArray[double]) – Velocity components on a doubly periodic grid, e.g. from physicskit.fluids.systems.navier_stokes.NavierStokes2D.velocity().
v (NDArray[double]) – Velocity components on a doubly periodic grid, e.g. from physicskit.fluids.systems.navier_stokes.NavierStokes2D.velocity().
length (float) – Physical domain size.

Return type:

Returns:

k (ndarray of float, shape (n // 2,)) – Wavenumber bins (angular, \(2\pi/\lambda\) convention).
E (ndarray of float, shape (n // 2,)) – Kinetic energy spectral density at each wavenumber bin.

Raises:

InvalidParameterError – If u and v are not the same shape, or are not square.

See also

kolmogorov_reference_slope: The -5/3 reference line to compare E against.

Examples

>>> import numpy as np
>>> n, length = 64, 2 * np.pi
>>> x = np.linspace(0, length, n, endpoint=False)
>>> X, Y = np.meshgrid(x, x, indexing="ij")
>>> u, v = np.sin(Y), np.sin(X)  # a single-mode flow
>>> k, E = energy_spectrum(u, v, length)
>>> bool(np.argmax(E) == 1)  # all the energy sits at the single k=1 mode
True
>>> round(float(E.sum()), 6)  # = mean of (u**2 + v**2) / 2
0.5

physicskit.fluids.exact_riemann_solution(x, t, left=(1.0, 0.0, 1.0), right=(0.125, 0.0, 0.1), x0=0.0, gamma=1.4)[source]#

Exact solution of the Riemann problem for the 1D Euler equations of an ideal gas.

The star-region pressure \(p_*\) is the root of \(f_L(p_*) + f_R(p_*) + u_R - u_L = 0\), where \(f_K\) is the shock (Rankine-Hugoniot) or rarefaction (isentropic) relation for side \(K\), found by Newton iteration; the solution is then sampled in the similarity variable \((x - x_0)/t\) (Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd ed., 2009, ch. 4). Vacuum generation is not handled.

Parameters:

x (NDArray[double]) – Positions.
t (float) – Time, positive.
left (tuple[float, float, float]) – (rho, u, p) on either side of the initial discontinuity.
right (tuple[float, float, float]) – (rho, u, p) on either side of the initial discontinuity.
x0 (float) – Initial position of the discontinuity.
gamma (float)

Return type:

Returns:

dict – "rho", "u", "p", the specific internal energy "e" at x, and the star-region "p_star" and "u_star".

Examples

>>> out = exact_riemann_solution(np.array([0.0]), 0.2)
>>> round(float(out["p_star"]), 5), round(float(out["u_star"]), 5)  # Sod: Toro Table 4.3
(0.30313, 0.92745)

physicskit.fluids.flow_past_cylinder(U_inf, radius, circulation=0.0)[source]#

Build the classic potential flow past a circular cylinder, with optional lift.

Superposing a uniform stream with a doublet of strength \(\kappa = 2\pi U_\infty R^2\) places a circular streamline of radius radius exactly at the origin – since \(\psi=0\) on that circle can be shown to be exactly the streamfunction of this combination, the circle itself can be interpreted as a solid boundary (potential flow does not enforce a no-slip condition, only no penetration). Adding a point vortex at the same location keeps the circular streamline intact (a vortex centered on the circle induces purely tangential velocity on it) while breaking front-back symmetry, which is exactly what a real cylinder’s boundary layer does when it separates asymmetrically – e.g. from a spinning cylinder (the Magnus effect).

Parameters:

U_inf (float) – Free-stream speed.
radius (float) – Cylinder radius.
circulation (float) – Circulation \(\Gamma\) added at the cylinder center; 0.0 gives the symmetric, lift-free flow.

Return type:

PotentialFlow

Returns:

PotentialFlow – The assembled flow: uniform stream + doublet (+ vortex if circulation is nonzero).

Raises:

InvalidParameterError – If U_inf or radius is not positive.

See also

kutta_joukowski_lift: The lift force generated by circulation.

Examples

>>> flow = flow_past_cylinder(U_inf=1.0, radius=1.0, circulation=4 * np.pi)
>>> theta = np.linspace(0, 2 * np.pi, 9)
>>> X, Y = 1.0 * np.cos(theta), 1.0 * np.sin(theta)
>>> psi = flow.streamfunction(X, Y)
>>> bool(np.max(np.abs(psi - psi[0])) < 1e-8)  # the cylinder surface is one streamline
True

physicskit.fluids.froude_number(velocity, length, g=9.81)[source]#

The Froude number: inertial forces over gravitational forces.

\[Fr = \frac{U}{\sqrt{gL}}\]

Governs free-surface (gravity wave) flows the way the Reynolds number governs viscous ones: a ship model and its full-size counterpart make the same waves, relative to their length, only if Fr matches. William Froude’s 1868 ship-hull towing-tank experiments established exactly this scaling law.

Parameters:

velocity (float) – Characteristic flow speed U.
length (float) – Characteristic length scale L.
g (float) – Gravitational acceleration.

Return type:

Returns:

float – Froude number Fr.

Raises:

InvalidParameterError – If length or g is not positive.

Examples

>>> round(froude_number(velocity=3.0, length=9.81, g=9.81), 4)
0.3058

physicskit.fluids.gaussian_heat_solution(coords, t, alpha, sigma0, amplitude=1.0, center=None)[source]#

Closed-form spreading of a Gaussian under the heat equation in free space.

A Gaussian of width \(\sigma_0\) convolved with the heat kernel stays Gaussian, with \(\sigma^2(t) = \sigma_0^2 + 2\alpha t\) and its peak lowered so that the total heat is conserved:

\[u(\mathbf x, t) = A\left(\frac{\sigma_0^2}{\sigma^2(t)}\right)^{d/2} \exp\!\left(-\frac{|\mathbf x - \mathbf x_0|^2}{2\sigma^2(t)}\right)\]

in \(d\) dimensions (e.g. Crank, The Mathematics of Diffusion, 2nd ed., 1975, eq. 3.5a).

Parameters:

coords (ndarray or sequence of ndarray) – One coordinate array (1D) or one per axis (e.g. from numpy.meshgrid() with indexing="ij").
t (float) – Time.
alpha (float) – Diffusivity.
sigma0 (float) – Initial width.
amplitude (float) – Initial peak value \(A\).
center (sequence of float, optional) – Center \(\mathbf x_0\); the origin by default.

Return type:

Returns:

ndarray – \(u(\mathbf x, t)\).

Examples

>>> float(gaussian_heat_solution(np.array([0.0]), t=1.5, alpha=1.0, sigma0=1.0)[0])
0.5

physicskit.fluids.geostrophic_adjustment_steady(eta0, x, g=9.81, H=1.0, f=1.0)[source]#

Balanced end state of Rossby’s adjustment problem, from potential-vorticity conservation.

The linear potential vorticity \(\partial_x v - f\eta/H\) is conserved, and the final state is in geostrophic balance, \(fv = g\partial_x\eta\), \(u = 0\), so

\[\eta - L_d^2\,\partial_x^2\eta = \eta_0, \qquad L_d = \sqrt{gH}/f.\]

For a step of height \(2\eta_0\) this gives the front \(\eta = \eta_0\,\mathrm{sgn}(x)(1 - e^{-|x|/L_d})\).

Parameters:

eta0 (NDArray[double]) – Initial displacement on the periodic grid x.
x (NDArray[double])
g (float)
H (float)
f (float)

Return type:

Returns:

dict – "eta" and "v" of the adjusted state, and the Rossby radius "L_d".

Examples

>>> x = np.linspace(0, 100, 512, endpoint=False)
>>> out = geostrophic_adjustment_steady(np.cos(2 * np.pi * x / 100), x, g=1.0)
>>> round(float(out["eta"].max()), 4)  # 1 / (1 + (L_d k)^2)
0.9961

physicskit.fluids.heat_equation(u0, spacing, alpha, t_max, dt, method='crank_nicolson', bc='dirichlet', n_frames=1)[source]#

Evolve a field under the heat equation with an explicit or Crank-Nicolson scheme.

With \(L\) the finite-difference Laplacian of laplacian_matrix(),

explicit (FTCS): \(u^{n+1} = u^n + \alpha\Delta t\,L u^n\);
Crank-Nicolson: \((I - \tfrac12\alpha\Delta t L)u^{n+1} = (I + \tfrac12\alpha\Delta t L)u^n\).

Parameters:

u0 (ArrayLike) – Initial field on a 1D, 2D, or 3D grid. With bc="dirichlet" its boundary values are held fixed.
spacing (float or sequence of float) – Grid spacing, one value or one per axis.
alpha (float) – Diffusivity (thermal diffusivity \(k/\rho c_p\), or a diffusion coefficient).
t_max (float) – Final time, rounded to a whole number of steps per frame.
dt (float) – Time step.
method (str) – Time-stepping scheme.
bc (str) – Boundary condition on every face.
n_frames (int) – Number of recorded frames after the initial one.

Return type:

Returns:

times (ndarray, shape (n_frames + 1,))
frames (ndarray, shape (n_frames + 1, *u0.shape))

Raises:

InvalidParameterError – If alpha or dt is not positive, the method or boundary condition is unknown, or the explicit scheme’s stability bound \(\alpha\Delta t\sum_a \Delta x_a^{-2} \le 1/2\) is violated.

Examples

A sine mode between cold walls decays as \(e^{-\alpha k^2 t}\), up to the grid’s \(O(k^2\Delta x^2)\) error in the decay rate:

>>> x = np.linspace(0.0, 1.0, 101)
>>> t, u = heat_equation(np.sin(np.pi * x), 0.01, alpha=0.1, t_max=1.0, dt=0.01)
>>> round(float(u[-1, 50] / np.exp(-0.1 * np.pi**2)), 3)
1.0

physicskit.fluids.heat_equation_spectral(u0, spacing, alpha, times)[source]#

Exact solution of the periodic heat equation by Fourier modes.

\[\hat u(\mathbf k, t) = \hat u(\mathbf k, 0)\,e^{-\alpha |\mathbf k|^2 t}\]

Each Fourier mode decays independently at its own rate, with no time stepping and no time-step error: Fourier’s 1807/1822 solution method.

Parameters:

u0 (ArrayLike) – Initial field on a periodic grid of any dimension.
spacing (float or sequence of float) – Grid spacing, one value or one per axis (the period along each axis is n * spacing).
alpha (float) – Diffusivity.
times (ArrayLike) – Times at which to return the field.

Return type:

Returns:

ndarray, shape (len(times), *u0.shape)

Examples

The mean is conserved while a mode of wavenumber \(k\) decays:

>>> x = np.arange(64) * (2 * np.pi / 64)
>>> u = heat_equation_spectral(1.0 + np.cos(3 * x), 2 * np.pi / 64, 0.5, [0.0, 1.0])
>>> round(float(u[1].mean()), 12), round(float(u[1, 0] - 1.0), 6), round(float(np.exp(-4.5)), 6)
(1.0, 0.011109, 0.011109)

physicskit.fluids.ideal_gas_sound_speed(gamma, pressure, density)[source]#

Adiabatic sound speed of an ideal gas, \(c = \sqrt{\gamma p/\rho}\).

Laplace’s (1816) correction to Newton’s isothermal \(\sqrt{p/\rho}\), which had underestimated the measured speed by about 15%.

Parameters:

gamma (float) – Ratio of specific heats (1.4 for air).
pressure (float) – Ambient pressure in Pa.
density (float) – Ambient density in kg/m^3.

Return type:

Returns:

float – Sound speed in m/s.

Examples

>>> round(ideal_gas_sound_speed(1.4, 101325.0, 1.204), 1)  # air at 20 C
343.2

physicskit.fluids.kelvin_helmholtz_growth_rate(k, delta_u)[source]#

Inviscid linear growth rate of a Kelvin-Helmholtz-unstable vortex sheet.

For a vortex sheet (an infinitesimally thin shear layer with a velocity jump delta_u across it), linear stability analysis of the Euler equations gives a perturbation of wavenumber k growing as \(e^{\sigma t}\) with

\[\sigma(k) = \frac{k\,\Delta u}{2},\]

unstable at every wavenumber with no threshold velocity – unlike, say, the Rayleigh-Taylor instability below, which is stabilized at short wavelength by surface tension or diffusion. A shear layer of finite thickness shear_width (as built by kelvin_helmholtz_ic()) instead has a fastest-growing wavelength comparable to the layer thickness, with growth cut off entirely for \(k\,(\text{shear width})\) large (Michalke, 1964); this vortex-sheet limit is the thin-layer, small-k approximation to that fuller theory, and is what the shear-layer roll-up in kelvin_helmholtz_ic() grows at initially, before finite-thickness and nonlinear effects take over.

Parameters:

k (float) – Perturbation wavenumber (angular, \(2\pi/\lambda\)).
delta_u (float) – Velocity jump across the shear layer.

Return type:

Returns:

float – Growth rate \(\sigma\), in units of inverse time.

Examples

>>> round(kelvin_helmholtz_growth_rate(k=1.0, delta_u=2.0), 6)
1.0

physicskit.fluids.kelvin_helmholtz_ic(n, length, shear_width=0.1, perturbation_amplitude=0.05)[source]#

Vorticity initial condition for the Kelvin-Helmholtz shear-layer instability.

A thin vortex sheet at mid-domain, with a small sinusoidal ripple that seeds the instability: under physicskit.fluids.systems.navier_stokes.NavierStokes2D, the ripple’s amplitude grows – initially at the rate predicted by kelvin_helmholtz_growth_rate() – as the shear layer rolls up into a row of discrete “cat’s eye” vortices.

Parameters:

n (int) – Number of grid points along each axis.
length (float) – Physical domain size.
shear_width (float) – Thickness of the vortex sheet.
perturbation_amplitude (float) – Amplitude of the seeding sinusoidal ripple.

Return type:

Returns:

ndarray of float, shape (n, n) – Vorticity field.

Examples

>>> omega0 = kelvin_helmholtz_ic(64, 2 * 3.141592653589793)
>>> omega0.shape
(64, 64)
>>> bool(omega0.max() > 0)
True

physicskit.fluids.kolmogorov_kraichnan_spectrum(k, epsilon, C=6.0)[source]#

Inertial-range energy spectrum \(E(k) = C\varepsilon^{2/3}k^{-5/3}\).

Parameters:

k (NDArray[double]) – Wavenumbers.
epsilon (float) – Energy flux through the inertial range.
C (float) – Kolmogorov constant; about 1.5 in 3D, and about 6 for the 2D inverse cascade.

Return type:

Returns:

ndarray

Examples

>>> float(kolmogorov_kraichnan_spectrum(1.0, epsilon=1.0))
6.0

physicskit.fluids.kolmogorov_reference_slope(k, k0, E0)[source]#

A Kolmogorov \(k^{-5/3}\) reference line, anchored at one point.

Kolmogorov’s K41 theory predicts that, in the inertial range, the energy spectrum depends only on the wavenumber k and the (scale-independent) energy dissipation rate, and dimensional analysis alone then fixes the power law:

\[E(k) = E_0 \left(\frac{k}{k_0}\right)^{-5/3}.\]

Plotting this line through one point (k0, E0) of a measured spectrum (see energy_spectrum()) on log-log axes is the standard visual check for an inertial range: a genuine turbulent cascade runs parallel to this line over at least a decade of k.

Parameters:

k (NDArray[double]) – Wavenumbers at which to evaluate the reference line.
k0 (float) – Anchor wavenumber; must be positive.
E0 (float) – Spectral energy density at k0.

Return type:

Returns:

ndarray of float – Reference spectral density \(E_0(k/k_0)^{-5/3}\) at each k.

Raises:

InvalidParameterError – If k0 is not positive.

Examples

>>> import numpy as np
>>> k = np.array([1.0, 2.0, 4.0])
>>> ref = kolmogorov_reference_slope(k, k0=1.0, E0=1.0)
>>> round(float(ref[1]), 4)
0.315

physicskit.fluids.kutta_joukowski_lift(rho, U_inf, circulation)[source]#

The Kutta-Joukowski lift theorem: lift per unit span from bound circulation.

\[L' = -\rho\,U_\infty\,\Gamma\]

Any 2D body generating a net circulation \(\Gamma\) around itself in a stream of speed \(U_\infty\) experiences a lift force per unit span of exactly this magnitude, directed perpendicular to the free stream – true regardless of the body’s shape, a remarkable consequence of potential theory (Kutta 1902, Zhukovsky 1906) that underlies all of classical airfoil theory. Real airfoils select circulation via the Kutta condition (smooth flow off a sharp trailing edge); the cylinder with an added point vortex built by flow_past_cylinder() sets it directly. With this module’s convention (free stream toward \(+x\), \(\Gamma>0\) counterclockwise, as in point_vortex_potential()) the lift is signed along \(+y\): a clockwise circulation (\(\Gamma<0\)), which speeds up the flow over the top, lifts upward, as on a conventional airfoil.

Parameters:

rho (float) – Fluid density.
U_inf (float) – Free-stream speed.
circulation (float) – Circulation \(\Gamma\) (positive counterclockwise) bound to the body.

Return type:

Returns:

float – Lift force per unit span along \(+y\), \(L'\).

Raises:

InvalidParameterError – If rho or U_inf is not positive.

Examples

>>> round(kutta_joukowski_lift(rho=1.2, U_inf=10.0, circulation=-5.0), 1)
60.0

physicskit.fluids.lame_parameters(E, nu)[source]#

Lamé’s first parameter \(\lambda\) and shear modulus \(\mu\) from \(E\) and \(\nu\).

\[\lambda = \frac{E\nu}{(1+\nu)(1-2\nu)}, \qquad \mu = \frac{E}{2(1+\nu)}\]

Parameters:

E (float) – Young’s modulus.
nu (float) – Poisson’s ratio, in (-1, 1/2).

Return type:

tuple[float, float]

Returns:

lam, mu (float)

Examples

>>> lame_parameters(2.5, 0.25)
(1.0, 1.0)

physicskit.fluids.lame_thick_cylinder(r, inner_radius, outer_radius, p_in, p_out=0.0, E=None, nu=None, plane='strain')[source]#

Lamé’s solution for a thick-walled cylinder under internal and external pressure.

With \(a, b\) the inner and outer radii,

\[\sigma_{rr} = A - \frac{B}{r^2}, \quad \sigma_{\theta\theta} = A + \frac{B}{r^2}, \quad A = \frac{p_i a^2 - p_o b^2}{b^2 - a^2}, \quad B = \frac{(p_i - p_o)a^2 b^2}{b^2 - a^2}\]

(Lamé, Leçons sur la théorie mathématique de l’élasticité des corps solides, 1852; Timoshenko and Goodier, Theory of Elasticity, 3rd ed., §28), and the radial displacement

\[u_r = \frac{1+\nu}{E}\left[(1-2\nu)Ar + \frac{B}{r}\right] \;\text{(plane strain)}, \qquad u_r = \frac{1}{E}\left[(1-\nu)Ar + (1+\nu)\frac{B}{r}\right] \;\text{(plane stress)}.\]

Parameters:

r (ArrayLike) – Radius, \(a \le r \le b\).
inner_radius (float) – \(a\) and \(b\).
outer_radius (float) – \(a\) and \(b\).
p_in (float) – Internal pressure \(p_i\).
p_out (float) – External pressure \(p_o\).
E (float | None) – Elastic constants; needed only for the displacement.
nu (float | None) – Elastic constants; needed only for the displacement.
plane (str) – Idealization for the displacement (the in-plane stresses are the same for both).

Returns:

sigma_rr, sigma_tt (ndarray)
u_r (ndarray or None) – None when E or nu is not given.

Examples

The hoop stress at the bore of a cylinder with \(b = 2a\):

>>> srr, stt, _ = lame_thick_cylinder([1.0, 2.0], 1.0, 2.0, p_in=3.0)
>>> srr.tolist(), stt.tolist()
([-3.0, 0.0], [5.0, 2.0])

physicskit.fluids.laplacian_matrix(shape, spacing, bc='dirichlet')[source]#

Second-order finite-difference Laplacian on a uniform grid, as a sparse matrix.

Built as a Kronecker sum of 1D three-point stencils \((u_{i+1} - 2u_i + u_{i-1})/\Delta x^2\), acting on the grid flattened in C order.

Parameters:

shape (tuple of int) – Grid shape (1 to 3 axes).
spacing (float or sequence of float) – Grid spacing, one value or one per axis.
bc (str) – "dirichlet" zeroes the rows of every boundary node, so those values stay fixed under heat_equation(); "neumann" imposes zero normal gradient (an insulated wall) with a mirrored ghost node; "periodic" wraps each axis.

Return type:

csr_matrix

Returns:

scipy.sparse.csr_matrix, shape (N, N) – With N = prod(shape).

Examples

>>> L = laplacian_matrix((5,), 1.0, bc="periodic")
>>> L.toarray()[0].tolist()
[-2.0, 1.0, 0.0, 0.0, 1.0]

physicskit.fluids.lbm_equilibrium(rho, ux, uy)[source]#

D2Q9 equilibrium distributions \(f_i^{\text{eq}}(\rho, \mathbf u)\).

\[f_i^{\text{eq}} = w_i\rho\left[1 + 3\,\mathbf e_i\cdot\mathbf u + \tfrac92(\mathbf e_i\cdot\mathbf u)^2 - \tfrac32 u^2\right]\]

Its zeroth and first moments are exactly \(\rho\) and \(\rho\mathbf u\).

Parameters:

rho (ArrayLike) – Density and velocity fields (broadcastable to a common shape).
ux (ArrayLike) – Density and velocity fields (broadcastable to a common shape).
uy (ArrayLike) – Density and velocity fields (broadcastable to a common shape).

Return type:

Returns:

ndarray, shape (9, *shape)

Examples

>>> feq = lbm_equilibrium(1.2, 0.05, -0.02)
>>> round(float(feq.sum()), 12), round(float(feq @ D2Q9_VELOCITIES[:, 0]), 12)
(1.2, 0.06)

physicskit.fluids.lbm_relaxation_time(viscosity)[source]#

Relaxation time giving a target viscosity, \(\tau = \nu/c_s^2 + 1/2\).

Parameters:: viscosity (float) – Kinematic viscosity in lattice units, positive.
Return type:: float
Returns:: float

Examples

>>> lbm_relaxation_time(1 / 6)
1.0

physicskit.fluids.lbm_viscosity(tau)[source]#

Kinematic viscosity of the BGK lattice Boltzmann model, \(\nu = c_s^2(\tau - 1/2)\).

Parameters:: tau (float) – Relaxation time, \(\tau > 1/2\).
Return type:: float
Returns:: float – Viscosity in lattice units.

Examples

>>> lbm_viscosity(1.0)
0.16666666666666666

physicskit.fluids.mach_number(velocity, speed_of_sound)[source]#

The Mach number: flow speed over the local speed of sound.

\[M = \frac{U}{c}\]

The single number that decides whether compressibility can be ignored (\(M \ll 1\), the regime of every other module in this package except physicskit.fluids.systems.compressible_flow) or dominates the flow (\(M \gtrsim 1\), where normal_shock_relations() applies).

Parameters:

velocity (float) – Flow speed U.
speed_of_sound (float) – Local speed of sound c.

Return type:

Returns:

float – Mach number M.

Raises:

InvalidParameterError – If speed_of_sound is not positive.

Examples

>>> mach_number(velocity=340.0, speed_of_sound=340.0)
1.0

physicskit.fluids.normal_shock_relations(M1, gamma=1.4)[source]#

Ideal-gas normal shock relations as a function of upstream Mach number.

Solving the Rankine-Hugoniot conditions (rankine_hugoniot_jump_conditions()) for an ideal gas gives the downstream state entirely in terms of the upstream Mach number \(M_1 = u_1/c_1\):

\[\frac{p_2}{p_1} = 1 + \frac{2\gamma}{\gamma+1}(M_1^2-1), \qquad \frac{\rho_2}{\rho_1} = \frac{(\gamma+1)M_1^2}{(\gamma-1)M_1^2+2}, \qquad M_2^2 = \frac{1+\tfrac{\gamma-1}{2}M_1^2}{\gamma M_1^2-\tfrac{\gamma-1}{2}}.\]

Only \(M_1 \geq 1\) gives a physically admissible (entropy-increasing) shock; a supersonic upstream flow needs a mechanism, such as this jump, to return to subsonic downstream, and no analogous jump exists in reverse.

Parameters:

M1 (float) – Upstream Mach number; must be at least 1.
gamma (float) – Ratio of specific heats (1.4 for a diatomic ideal gas, e.g. air).

Return type:

dict[str, float]

Returns:

dict – {"p2_p1": ..., "rho2_rho1": ..., "T2_T1": ..., "M2": ...}.

Raises:

InvalidParameterError – If M1 is less than 1.

Examples

>>> jump = normal_shock_relations(M1=1.0)
>>> [round(v, 6) for v in (jump["p2_p1"], jump["rho2_rho1"], jump["M2"])]
[1.0, 1.0, 1.0]
>>> jump = normal_shock_relations(M1=2.0)
>>> round(jump["p2_p1"], 3), round(jump["rho2_rho1"], 3), round(jump["M2"], 3)
(4.5, 2.667, 0.577)

physicskit.fluids.pipe_mode_frequencies(length, c, n_modes=5, ends=('open', 'open'))[source]#

Standing-wave frequencies of a pipe with closed (rigid) or open (pressure-release) ends.

Parameters:

length (float) – Pipe length \(L\).
c (float) – Sound speed.
n_modes (int) – Number of modes, lowest first.
ends (tuple of {"closed", "open"}, default=("open", "open")) – Boundary condition at each end.

Return type:

Returns:

ndarray, shape (n_modes,) – \(f_n = nc/2L\) for like ends and \((2n-1)c/4L\) for unlike ends. (A pipe closed at both ends also has a trivial \(f = 0\) uniform-pressure mode, which is not listed.)

Examples

>>> pipe_mode_frequencies(1.0, 340.0, 3, ends=("closed", "open")).tolist()
[85.0, 255.0, 425.0]

physicskit.fluids.pipe_mode_shape(x, length, mode, ends=('open', 'open'))[source]#

Pressure profile of a pipe’s mode-th standing wave (unit amplitude).

Parameters:

x (ArrayLike) – Position along the pipe, 0 <= x <= length.
length (float) – Pipe length.
mode (int) – Mode number, starting at 1 (the order of pipe_mode_frequencies()).
ends (tuple of {"closed", "open"}, default=("open", "open")) – Boundary condition at each end.

Return type:

Returns:

ndarray – \(\sin(n\pi x/L)\) (open-open), \(\cos(n\pi x/L)\) (closed-closed), \(\cos((2n-1)\pi x/2L)\) (closed-open), or \(\sin((2n-1)\pi x/2L)\) (open-closed).

Examples

>>> pipe_mode_shape([0.0, 0.5, 1.0], 1.0, 1, ends=("closed", "closed")).round(12).tolist()
[1.0, 0.0, -1.0]

physicskit.fluids.plane_elasticity_solve(nodes, triangles, E, nu, *, plane='stress', fixed=None, fixed_values=None, forces=None, thickness=1.0)[source]#

Solve 2D linear elasticity with constant-strain triangles.

Each triangle contributes the stiffness \(K_e = t A_e B_e^{\mathsf T} D B_e\), with \(B_e\) the (constant) strain-displacement matrix of linear shape functions and \(D\) from elasticity_matrix() (Zienkiewicz and Taylor §4.2). The global system \(K\mathbf u = \mathbf f\) is solved for the free degrees of freedom after moving prescribed displacements to the right-hand side. The element passes the patch test: any uniform strain field is reproduced exactly.

Parameters:

nodes (ArrayLike) – Node coordinates.
triangles (ArrayLike) – Counter-clockwise node indices.
E (float) – Young’s modulus and Poisson’s ratio.
nu (float) – Young’s modulus and Poisson’s ratio.
plane (str) – Plane idealization.
fixed (ArrayLike | None) – Which displacement components are prescribed. At least enough must be fixed to remove rigid-body motion.
fixed_values (ArrayLike | None) – Prescribed displacement values (zero by default).
forces (ArrayLike | None) – External nodal forces, e.g. from pressure_load().
thickness (float) – Out-of-plane thickness (per unit length for plane strain).

Return type:

PlaneElasticityResult

Returns:

PlaneElasticityResult

Examples

Uniaxial tension of a plate: \(\sigma_{xx} = \sigma\) everywhere and \(u_x = \sigma x / E\). The loaded right edge is listed top to bottom, so the plate lies on its right, and a negative pressure pulls:

>>> nodes, tris = rectangle_mesh(2.0, 1.0, 4, 2)
>>> fixed = np.zeros((len(nodes), 2), dtype=bool)
>>> fixed[nodes[:, 0] == 0, 0] = True
>>> fixed[0, 1] = True
>>> right = np.flatnonzero(nodes[:, 0] == 2.0)
>>> F = pressure_load(nodes, np.column_stack([right[1:], right[:-1]]), -3.0)
>>> res = plane_elasticity_solve(nodes, tris, 100.0, 0.3, fixed=fixed, forces=F)
>>> bool(np.allclose(res.stress[:, 0], 3.0)), round(float(res.displacement[right, 0].mean()), 6)
(True, 0.06)

physicskit.fluids.plot_energy_spectrum(k, E, ax=None, show_kolmogorov=True, **plot_kwargs)[source]#

Log-log plot of a kinetic energy spectrum, with an optional Kolmogorov -5/3 reference line.

Parameters:

k (NDArray[double]) – Wavenumbers, e.g. from physicskit.fluids.utils.spectral_analysis.energy_spectrum().
E (NDArray[double]) – Spectral energy density at each k.
ax (Axes | None) – Axes to draw into; a new figure is created if omitted.
show_kolmogorov (bool) – If True, overlay a \(k^{-5/3}\) reference line anchored at the lowest plotted (nonzero-energy) wavenumber.
**plot_kwargs (Any) – Additional keyword arguments forwarded to ax.loglog for the measured spectrum.

Return type:

tuple[Figure, Axes]

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

physicskit.fluids.plot_pressure_coefficient(theta, Cp, ax=None, **plot_kwargs)[source]#

Plot the surface pressure coefficient \(C_p(\theta)\) around a body.

Parameters:

theta (NDArray[double]) – Angular position around the body (radians), typically [0, 2*pi].
Cp (NDArray[double]) – Pressure coefficient at each theta, e.g. from physicskit.fluids.systems.potential_flow.pressure_coefficient() evaluated on a circle.
ax (Axes | None) – Axes to draw into; a new figure is created if omitted.
**plot_kwargs (Any) – Additional keyword arguments forwarded to ax.plot.

Return type:

tuple[Figure, Axes]

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

physicskit.fluids.plot_shock_tube_profiles(x, rho, u, p)[source]#

Plot density, velocity, and pressure profiles from a shock-tube solution.

Parameters:

x (NDArray[double]) – Spatial grid.
rho (NDArray[double]) – Density, velocity, and pressure profiles, e.g. from physicskit.fluids.systems.compressible_flow.sod_shock_tube().
u (NDArray[double]) – Density, velocity, and pressure profiles, e.g. from physicskit.fluids.systems.compressible_flow.sod_shock_tube().
p (NDArray[double]) – Density, velocity, and pressure profiles, e.g. from physicskit.fluids.systems.compressible_flow.sod_shock_tube().

Return type:

tuple[Figure, NDArray]

Returns:

fig (matplotlib.figure.Figure)
axes (ndarray of matplotlib.axes.Axes, shape (3,))

physicskit.fluids.plot_streamlines(X, Y, u, v, ax=None, **streamplot_kwargs)[source]#

Plot velocity-field streamlines.

Parameters:

X (NDArray[double]) – Real-space coordinate grids.
Y (NDArray[double]) – Real-space coordinate grids.
u (NDArray[double]) – Velocity components.
v (NDArray[double]) – Velocity components.
ax (Axes | None) – Axes to draw into; a new figure is created if omitted.
**streamplot_kwargs (Any) – Additional keyword arguments forwarded to ax.streamplot.

Return type:

tuple[Figure, Axes]

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

physicskit.fluids.plot_vorticity_field(X, Y, omega, u=None, v=None, ax=None, **pcolormesh_kwargs)[source]#

Plot a vorticity field as a heatmap, optionally overlaid with velocity streamlines.

Parameters:

X (NDArray[double]) – Real-space coordinate grids.
Y (NDArray[double]) – Real-space coordinate grids.
omega (NDArray[double]) – Vorticity field.
u (NDArray[double] | None) – Velocity components; if given, streamlines are overlaid on the heatmap.
v (NDArray[double] | None) – Velocity components; if given, streamlines are overlaid on the heatmap.
ax (Axes | None) – Axes to draw into; a new figure is created if omitted.
**pcolormesh_kwargs (Any) – Additional keyword arguments forwarded to ax.pcolormesh.

Return type:

tuple[Figure, Axes]

Returns:

fig (matplotlib.figure.Figure)
ax (matplotlib.axes.Axes)

physicskit.fluids.point_vortex_potential(X, Y, circulation, x0=0.0, y0=0.0)[source]#

Complex potential of an isolated point vortex.

Parameters:

X (ArrayLike) – Field points.
Y (ArrayLike) – Field points.
circulation (float) – Circulation \(\Gamma\) (positive is counterclockwise).
x0 (float) – Vortex location.
y0 (float) – Vortex location.

Return type:

Returns:

ndarray of complex – Complex potential \(W(z) = \frac{-i\Gamma}{2\pi}\ln(z - z_0)\).

physicskit.fluids.point_vortex_velocities(positions, circulations, core=1e-06)[source]#

Velocity induced on each point vortex by every other vortex (Biot-Savart law).

The 2D Biot-Savart law: an isolated point vortex of circulation \(\Gamma_j\) at \(\mathbf{r}_j\) induces a purely azimuthal velocity \(\Gamma_j/(2\pi r)\) at distance r, and a system of n vortices simply superposes these fields (Laplace’s equation is linear, just as in physicskit.fluids.systems.potential_flow – indeed a point vortex here is point_vortex_potential(), evaluated at every other vortex’s location rather than on a field grid).

Parameters:

positions (ArrayLike) – Vortex positions (x, y).
circulations (ArrayLike) – Vortex circulations \(\Gamma_i\) (positive counterclockwise).
core (float) – Regularization (Rankine core) radius, keeping the induced velocity finite if two vortex cores coincide; negligible for any pair separated by much more than core.

Return type:

Returns:

ndarray of float, shape (n, 2) – Velocity at each vortex’s own position, induced by all the others.

Examples

A pair of counter-rotating vortices separated by distance d induces equal and opposite velocities of magnitude \(\Gamma/(2\pi d)\) on each other, translating the pair sideways rather than rotating it:

>>> import numpy as np
>>> Gamma, d = 1.0, 2.0
>>> positions = np.array([[0.0, d / 2], [0.0, -d / 2]])
>>> circulations = np.array([Gamma, -Gamma])
>>> vel = point_vortex_velocities(positions, circulations)
>>> bool(np.allclose(vel[0], vel[1]))  # both translate together
True
>>> round(float(vel[0, 0]), 6) == round(Gamma / (2 * np.pi * d), 6)
True

physicskit.fluids.poiseuille_flow_rate(dpdx, mu, h)[source]#

Volumetric flow rate per unit depth of plane Poiseuille flow.

The integral of poiseuille_flow_velocity() across the channel:

\[Q = -\frac{h^3}{12\mu}\frac{dp}{dx}\]

Parameters:

dpdx (float) – Imposed pressure gradient.
mu (float) – Dynamic viscosity.
h (float) – Gap between the plates.

Return type:

Returns:

float – Volumetric flow rate per unit depth, Q.

Raises:

InvalidParameterError – If mu or h is not positive.

Examples

>>> round(poiseuille_flow_rate(dpdx=-12.0, mu=1.0, h=1.0), 4)
1.0

physicskit.fluids.poiseuille_flow_velocity(y, dpdx, mu, h)[source]#

Velocity profile of steady plane Poiseuille flow.

The flow between two stationary infinite parallel plates driven by a constant imposed pressure gradient dpdx. The steady Navier-Stokes equations reduce to \(\mu\,d^2u/dy^2 = dp/dx\), whose solution subject to no-slip at both walls (y=0 and y=h) is the parabolic profile below, maximal at the channel centerline.

Parameters:

y (ArrayLike) – Height above the lower plate, 0 <= y <= h.
dpdx (float) – Imposed pressure gradient (negative drives flow in +x).
mu (float) – Dynamic viscosity.
h (float) – Gap between the plates.

Return type:

Returns:

ndarray of float – Velocity \(u(y) = -\frac{1}{2\mu}\frac{dp}{dx}\,y\,(h-y)\).

Raises:

InvalidParameterError – If mu or h is not positive.

See also

poiseuille_flow_rate: The volumetric flow rate driven by this profile.

Examples

>>> y = np.linspace(0, 1, 5)
>>> u = poiseuille_flow_velocity(y, dpdx=-8.0, mu=1.0, h=1.0)
>>> round(float(u[2]), 4)  # peak at centerline
1.0
>>> bool(u[0] == 0.0 and u[-1] == 0.0)  # no-slip at both walls
True

physicskit.fluids.polar_stress(stress, points)[source]#

Rotate Cartesian stresses to polar components about the origin.

\[\begin{split}\sigma_{rr} = \sigma_{xx}\cos^2\theta + \sigma_{yy}\sin^2\theta + 2\sigma_{xy}\sin\theta\cos\theta,\\ \sigma_{\theta\theta} = \sigma_{xx}\sin^2\theta + \sigma_{yy}\cos^2\theta - 2\sigma_{xy}\sin\theta\cos\theta,\\ \sigma_{r\theta} = (\sigma_{yy} - \sigma_{xx})\sin\theta\cos\theta + \sigma_{xy}(\cos^2\theta - \sin^2\theta).\end{split}\]

Parameters:

stress (ArrayLike) – \((\sigma_{xx}, \sigma_{yy}, \sigma_{xy})\).
points (ArrayLike) – Where each stress applies.

Return type:

tuple[NDArray[double], NDArray[double], NDArray[double]]

Returns:

sigma_rr, sigma_tt, sigma_rt (ndarray, shape (n,))

Examples

>>> [round(float(v[0]), 12) for v in polar_stress([[1.0, 2.0, 0.0]], [[0.0, 1.0]])]
[2.0, 1.0, 0.0]

physicskit.fluids.pressure_coefficient(u, v, U_inf)[source]#

Bernoulli’s pressure coefficient for steady, incompressible potential flow.

Along a streamline of a steady, incompressible, inviscid flow, Bernoulli’s equation \(p + \tfrac{1}{2}\rho|\mathbf{u}|^2 = \mathrm{const}\) gives the dimensionless pressure coefficient

\[C_p = \frac{p - p_\infty}{\tfrac{1}{2}\rho U_\infty^2} = 1 - \frac{u^2+v^2}{U_\infty^2}.\]

Parameters:

u (ArrayLike) – Velocity components.
v (ArrayLike) – Velocity components.
U_inf (float) – Free-stream speed.

Return type:

Returns:

ndarray of float – Pressure coefficient.

Raises:

InvalidParameterError – If U_inf is not positive.

Examples

>>> import numpy as np
>>> bool(np.isclose(pressure_coefficient(u=1.0, v=0.0, U_inf=1.0), 0.0))
True

physicskit.fluids.pressure_load(nodes, edges, pressure, thickness=1.0)[source]#

Consistent nodal forces from a uniform pressure on boundary edges.

Each directed edge \(a \to b\) carries a normal traction of magnitude pressure pushing to its right, i.e. into the solid when edges are listed with the solid on their right (as annulus_mesh() returns them). For a linear element the traction splits equally between the two end nodes: each receives \(\tfrac12 p\,t\,\ell\,\hat{\mathbf n}\).

Parameters:

nodes (ArrayLike)
edges (ArrayLike)
pressure (float)
thickness (float)

Return type:

Returns:

ndarray, shape (n_nodes, 2)

Examples

>>> F = pressure_load([[0, 0], [1, 0]], [[0, 1]], 2.0)
>>> F.tolist()  # pushes toward -y, the right of the edge direction +x
[[0.0, -1.0], [0.0, -1.0]]

physicskit.fluids.rankine_hugoniot_jump_conditions(rho1, u1, p1, rho2, u2, p2, gamma=1.4)[source]#

Residuals of the Rankine-Hugoniot jump conditions across a stationary discontinuity.

A steady discontinuity in a 1D inviscid compressible flow (states 1 upstream, 2 downstream, in the frame where the discontinuity itself is at rest) must conserve mass, momentum, and energy flux exactly:

\[\rho_1 u_1 = \rho_2 u_2, \qquad p_1 + \rho_1 u_1^2 = p_2 + \rho_2 u_2^2, \qquad h_1 + \tfrac{1}{2}u_1^2 = h_2 + \tfrac{1}{2}u_2^2,\]

where \(h=\gamma p/((\gamma-1)\rho)\) is specific enthalpy (an ideal gas is assumed only for the energy residual). Rankine (1870) first wrote the mass and momentum conditions; Hugoniot (1887, 1889) added the energy condition and the resulting relation between the shock’s pressure and density ratios now called the Hugoniot curve. This function returns each condition’s residual, primarily to verify a candidate jump (e.g. the output of normal_shock_relations()) rather than to solve for one directly.

Parameters:

rho1 (float) – Upstream density, velocity, and pressure.
u1 (float) – Upstream density, velocity, and pressure.
p1 (float) – Upstream density, velocity, and pressure.
rho2 (float) – Downstream density, velocity, and pressure.
u2 (float) – Downstream density, velocity, and pressure.
p2 (float) – Downstream density, velocity, and pressure.
gamma (float) – Ratio of specific heats, used only in the enthalpy of the energy residual; must match the one used to generate a candidate jump (e.g. normal_shock_relations()’s gamma).

Return type:

dict[str, float]

Returns:

dict – {"mass": ..., "momentum": ..., "energy": ...} residuals; all zero for an exactly satisfied jump.

Examples

>>> jump = normal_shock_relations(M1=2.0)
>>> residuals = rankine_hugoniot_jump_conditions(1.0, 2.0, 1.0 / 1.4, jump["rho2_rho1"], 2.0 / jump["rho2_rho1"], jump["p2_p1"] / 1.4)
>>> bool(max(abs(v) for v in residuals.values()) < 1e-10)
True

physicskit.fluids.rayleigh_benard_critical(boundaries='rigid')[source]#

Critical Rayleigh number and wavenumber for the onset of convection.

Minimizes rayleigh_benard_marginal_rayleigh() over \(a\).

Parameters:: boundaries (str) – No-slip or stress-free horizontal walls (both isothermal).
Return type:: tuple[float, float]
Returns:: Ra_c, a_c (float) – \((27\pi^4/4, \pi/\sqrt2)\) exactly for "free"; about (1707.76, 3.117) for "rigid".

Examples

>>> Ra_c, a_c = rayleigh_benard_critical("rigid")
>>> round(Ra_c, 2), round(a_c, 3)
(1707.76, 3.116)

physicskit.fluids.rayleigh_benard_growth_rate_free(Ra, a, Pr=1.0, m=1)[source]#

Exact linear growth rate of the stress-free mode \(\sin(m\pi z)e^{iax}\).

With \(q^2 = m^2\pi^2 + a^2\), the growth rate is the larger root of

\[\sigma^2 + (1 + Pr)q^2\sigma + Pr\left(q^4 - Ra\,\frac{a^2}{q^2}\right) = 0.\]

Parameters:

Ra (float)
a (float) – Horizontal wavenumber.
Pr (float)
m (int) – Vertical mode number.

Return type:

Returns:

float – Growth rate, in units of \(\kappa/d^2\); zero at the marginal \(Ra = q^6/a^2\).

Examples

>>> a = np.pi / np.sqrt(2)
>>> abs(rayleigh_benard_growth_rate_free(27 * np.pi**4 / 4, a)) < 1e-9
True

physicskit.fluids.rayleigh_benard_marginal_rayleigh(a, boundaries='rigid')[source]#

Rayleigh number at which the mode of horizontal wavenumber \(a\) is marginally stable.

For "free" (stress-free) walls this is \((\pi^2 + a^2)^3/a^2\). For "rigid" (no-slip) walls the lowest, even mode is \(W = \sum_{j=0}^2 A_j \cosh(q_j z)\) on \(-\tfrac12 < z < \tfrac12\), with \((q_j^2 - a^2)^3 = -Ra\,a^2\), and \(Ra\) is the smallest root of the \(3\times3\) determinant imposing \(W = W' = \theta = 0\) at the wall (Chandrasekhar 1961, sec. 15).

Parameters:

a (float) – Horizontal wavenumber (in units of \(1/d\)).
boundaries (str)

Return type:

Returns:

float

Examples

>>> round(rayleigh_benard_marginal_rayleigh(3.117), 1)
1707.8

physicskit.fluids.rayleigh_taylor_growth_rate(k, atwood_number, g=1.0)[source]#

Linear growth rate of the Rayleigh-Taylor instability.

For two inviscid, semi-infinite fluid layers (heavy density \(\rho_2\) on top of light density \(\rho_1\), gravity g pointing from the heavy into the light layer), a perturbation of wavenumber k on the interface grows as \(e^{\sigma t}\) with

\[\sigma(k) = \sqrt{A\,g\,k}, \qquad A = \frac{\rho_2-\rho_1}{\rho_2+\rho_1}\]

(Rayleigh, 1883; Taylor, 1950), unstable at every wavenumber in the idealized inviscid limit – real fluids are cut off at short wavelength by viscosity and surface tension, which simulate_rayleigh_taylor() supplies dynamically via its buoyancy diffusivity.

Parameters:

k (float) – Perturbation wavenumber (angular, \(2\pi/\lambda\)).
atwood_number (float) – The Atwood number \(A\), \(0 < A < 1\).
g (float) – Gravitational acceleration, pointing from the heavy layer into the light one.

Return type:

Returns:

float – Growth rate \(\sigma\), in units of inverse time.

Raises:

InvalidParameterError – If atwood_number is not in (0, 1) or k is negative.

Examples

>>> round(rayleigh_taylor_growth_rate(k=2.0, atwood_number=0.5, g=1.0), 6)
1.0

physicskit.fluids.rayleigh_taylor_ic(n, length, atwood_number, perturbation_amplitude=0.02, interface_position=None)[source]#

Vorticity and buoyancy initial condition for the Rayleigh-Taylor instability.

A denser fluid layer sits atop a lighter one (heavy-on-top is exactly the unstable arrangement), separated by an interface at mid-domain with a small sinusoidal ripple, evolved through simulate_rayleigh_taylor(). Vorticity starts at exactly zero: in the Boussinesq system, all vorticity generation comes from the baroclinic torque acting on the rippled density interface once gravity is switched on (see physicskit.fluids.core.timestepping.buoyant_vorticity_rhs()).

Parameters:

n (int) – Number of grid points along each axis.
length (float) – Physical domain size.
atwood_number (float) – The Atwood number \(A=(\rho_{heavy}-\rho_{light})/(\rho_{heavy}+\rho_{light})\), \(0 < A < 1\), setting the buoyancy jump across the interface.
perturbation_amplitude (float) – Amplitude of the seeding sinusoidal ripple on the interface.
interface_position (float | None) – Vertical position of the unperturbed interface; defaults to length / 2.

Return type:

Returns:

omega0 (ndarray of float, shape (n, n)) – Initial vorticity field (identically zero).
buoyancy0 (ndarray of float, shape (n, n)) – Initial buoyancy field: \(+A\) below the interface (light fluid), \(-A\) above it (heavy fluid), smoothed over one grid cell so the spectral derivatives it feeds resolve without ringing.

Raises:

InvalidParameterError – If atwood_number is not in (0, 1).

Examples

>>> omega0, buoyancy0 = rayleigh_taylor_ic(64, 2 * 3.141592653589793, atwood_number=0.3)
>>> bool(np.all(omega0 == 0.0))
True
>>> bool(buoyancy0.max() > 0 > buoyancy0.min())
True

physicskit.fluids.rectangle_mesh(Lx, Ly, nx, ny)[source]#

Structured triangulation of the rectangle \([0, L_x]\times[0, L_y]\).

Each of the nx * ny cells is split into two counter-clockwise triangles.

Parameters:

Lx (float) – Rectangle size.
Ly (float) – Rectangle size.
nx (int) – Number of cells along each axis.
ny (int) – Number of cells along each axis.

Return type:

tuple[NDArray[double], NDArray[int_]]

Returns:

nodes (ndarray, shape ((nx + 1) * (ny + 1), 2)) – Node coordinates; node i * (ny + 1) + j is at (i Lx/nx, j Ly/ny).
triangles (ndarray of int, shape (2 * nx * ny, 3))

Examples

>>> nodes, tris = rectangle_mesh(2.0, 1.0, 2, 1)
>>> nodes.shape, tris.shape
((6, 2), (4, 3))

physicskit.fluids.rectangular_room_mode_frequencies(Lx, Ly, c, n_max=3)[source]#

Mode frequencies of a rigid-walled rectangular room (2D).

\[f_{mn} = \frac{c}{2}\sqrt{\left(\frac{m}{L_x}\right)^2 + \left(\frac{n}{L_y}\right)^2}\]

with pressure \(\cos(m\pi x/L_x)\cos(n\pi y/L_y)\) (Kinsler et al., §9.2).

Parameters:

Lx (float) – Room dimensions.
Ly (float) – Room dimensions.
c (float) – Sound speed.
n_max (int) – Largest mode index along each axis.

Return type:

Returns:

ndarray, shape (k, 3) – Rows (m, n, f_mn) sorted by frequency, excluding (0, 0).

Examples

>>> rectangular_room_mode_frequencies(4.0, 3.0, 340.0, n_max=1)[:2].tolist()
[[1.0, 0.0, 42.5], [0.0, 1.0, 56.666666666666664]]

physicskit.fluids.reynolds_number(velocity, length, nu)[source]#

The Reynolds number: inertial forces over viscous forces.

\[Re = \frac{U L}{\nu}\]

Osborne Reynolds’ 1883 pipe-flow experiments identified this ratio as the single parameter controlling the transition from smooth (laminar) to chaotic (turbulent) flow, regardless of the pipe’s size or the fluid’s identity – the founding result of dimensional-similarity reasoning in fluid mechanics.

Parameters:

velocity (float) – Characteristic flow speed U.
length (float) – Characteristic length scale L.
nu (float) – Kinematic viscosity.

Return type:

Returns:

float – Reynolds number Re.

Raises:

InvalidParameterError – If length or nu is not positive.

Examples

>>> reynolds_number(velocity=2.0, length=0.5, nu=1e-6)
1000000.0

physicskit.fluids.rotating_shallow_water_1d(eta0, x, times, g=9.81, H=1.0, f=1.0)[source]#

Exact evolution of the linear rotating shallow-water equations on a periodic line.

Each Fourier mode of \((u, v, \eta)\), starting from rest, splits into a steady geostrophic part and inertia-gravity (Poincaré) waves of frequency \(\omega_k = \sqrt{f^2 + gHk^2}\), whose evolution is computed exactly:

\[\hat\eta_k(t) = \hat\eta_k(0)\,\frac{f^2 + gHk^2\cos\omega_k t}{\omega_k^2},\quad \hat u_k = -\frac{igk\hat\eta_k(0)}{\omega_k}\sin\omega_k t,\quad \hat v_k = \frac{igkf\hat\eta_k(0)}{\omega_k^2}(1 - \cos\omega_k t).\]

Parameters:

eta0 (NDArray[double]) – Initial surface displacement, fluid at rest.
x (NDArray[double]) – Uniform periodic grid.
times (NDArray[double]) – Output times.
g (float) – Gravity, mean depth and Coriolis parameter.
H (float) – Gravity, mean depth and Coriolis parameter.
f (float) – Gravity, mean depth and Coriolis parameter.

Return type:

Returns:

dict – "eta", "u", "v", each of shape (len(times), n).

Examples

>>> x = np.linspace(0, 10, 64, endpoint=False)
>>> out = rotating_shallow_water_1d(np.cos(2 * np.pi * x / 10), x, [0.0], g=1.0)
>>> bool(np.allclose(out["eta"][0], np.cos(2 * np.pi * x / 10)))
True

physicskit.fluids.shallow_water_1d(h0, u0, x, t_max, g=9.81, cfl=0.45)[source]#

Solve the 1D shallow-water equations with a finite-volume HLL scheme.

Second order in space (minmod-limited MUSCL reconstruction of \(h\) and \(hu\)) and time (Heun), with transmissive boundaries. Dry cells (\(h = 0\)) are allowed.

Parameters:

h0 (NDArray[double]) – Initial depth and velocity at the cell centres x.
u0 (NDArray[double]) – Initial depth and velocity at the cell centres x.
x (NDArray[double]) – Uniformly spaced cell centres.
t_max (float) – Final time.
g (float)
cfl (float) – Courant number.

Return type:

Returns:

dict – "h" and "u" at t_max, and the number of "steps".

Examples

>>> x = np.linspace(-1, 1, 200)
>>> out = shallow_water_1d(np.where(x < 0, 1.0, 0.5), np.zeros_like(x), x, t_max=0.1)
>>> bool(abs(out["h"].mean() - 0.75) < 1e-2)  # mass is conserved away from the ends
True

physicskit.fluids.simulate_rayleigh_taylor(omega0, buoyancy0, nu, kappa, g, dt, steps, length)[source]#

Time-step the Boussinesq Rayleigh-Taylor system with RK4.

Evolves the coupled vorticity-buoyancy equations of physicskit.fluids.core.timestepping.buoyant_vorticity_rhs() (scaled by the gravitational acceleration g) starting from rayleigh_taylor_ic(), through the roll-up of the rippled interface into the characteristic Rayleigh-Taylor mushroom plumes.

Parameters:

omega0 (NDArray[double]) – Initial vorticity field, as from rayleigh_taylor_ic().
buoyancy0 (NDArray[double]) – Initial buoyancy field, as from rayleigh_taylor_ic().
nu (float) – Kinematic viscosity; must be positive.
kappa (float) – Buoyancy (density) diffusivity; must be positive.
g (float) – Gravitational acceleration; must be positive.
dt (float) – Time step.
steps (int) – Number of RK4 steps to advance.
length (float) – Physical domain size.

Return type:

Returns:

dict – {"omega": ..., "buoyancy": ..., "psi": ..., "u": ..., "v": ...} at the final time.

Raises:

InvalidParameterError – If nu, kappa, or g is not positive.

See also

rayleigh_taylor_ic: Builds the initial condition consumed here.
rayleigh_taylor_growth_rate: The linear-theory prediction this simulation exceeds nonlinearly.

physicskit.fluids.simulate_vorticity_streamfunction(omega0, nu, dt, steps, length)[source]#

Time-step 2D incompressible vorticity-streamfunction Navier-Stokes with RK4.

A one-shot functional convenience wrapper around NavierStokes2D, for scripts that only need a single run rather than repeated stepping against the same grid.

Parameters:

omega0 (NDArray[double]) – Initial vorticity field on a doubly periodic [0, length)^2 domain.
nu (float) – Kinematic viscosity; must be positive.
dt (float) – Time step.
steps (int) – Number of RK4 steps to advance.
length (float) – Physical domain size.

Return type:

Returns:

dict – {"omega": final vorticity, "psi": final streamfunction, "u": ..., "v": ...}.

Raises:

InvalidParameterError – If nu is not positive.

See also

NavierStokes2D: Reusable object form, for repeated stepping on one grid.
physicskit.fluids.systems.instabilities.kelvin_helmholtz_ic: Shear-layer initial condition that grows via this solver.

Examples

>>> import numpy as np
>>> n = 48
>>> from physicskit.fluids.core.grid import spectral_grid
>>> X, Y, KX, KY, K2 = spectral_grid(n, 2 * np.pi)
>>> omega0 = np.sin(X) * np.sin(Y)
>>> result = simulate_vorticity_streamfunction(omega0, nu=0.05, dt=0.01, steps=20, length=2 * np.pi)
>>> result["omega"].shape
(48, 48)

physicskit.fluids.sod_shock_tube(nx=400, x0=0.5, t_final=0.2, gamma=1.4, cfl=0.5)[source]#

Solve the classic Sod shock-tube problem with a Lax-Friedrichs finite-volume scheme.

The Sod (1978) problem is the standard test Riemann problem for a compressible-flow solver: a diaphragm at x0 initially separates high pressure/density gas (left, at rest) from low pressure/density gas (right, at rest); removing the diaphragm at \(t=0\) produces, for \(t>0\), exactly three simple waves fanning out from x0 – a left-running rarefaction (expansion) fan, a right-running contact discontinuity (a density jump with continuous pressure and velocity), and a right-running shock satisfying rankine_hugoniot_jump_conditions() – against which any finite-volume scheme’s numerical diffusion and shock-capturing behavior can be judged. The domain is [0, 1]; the classic Sod initial condition is (rho, u, p) = (1, 0, 1) for x < x0 and (0.125, 0, 0.1) for x >= x0.

Parameters:

nx (int) – Number of finite-volume cells.
x0 (float) – Initial diaphragm location.
t_final (float) – Time to integrate to.
gamma (float) – Ratio of specific heats.
cfl (float) – Courant number; the time step is chosen adaptively each macro-chunk of the integration as cfl * dx / max(|u| + c), and must satisfy 0 < cfl <= 1 for the explicit scheme to be stable.

Return type:

Returns:

dict – {"x": ..., "rho": ..., "u": ..., "p": ...} at t_final, each an array of length nx.

Raises:

InvalidParameterError – If nx is smaller than 4, t_final is not positive, or cfl is not in (0, 1].

Examples

>>> result = sod_shock_tube(nx=200, t_final=0.15)
>>> result["rho"].shape
(200,)
>>> bool(result["rho"][0] > result["rho"][-1] > 0)  # left state denser than right
True

physicskit.fluids.source_potential(X, Y, strength, x0=0.0, y0=0.0)[source]#

Complex potential of a point source (or, for negative strength, a sink).

Parameters:

X (ArrayLike) – Field points.
Y (ArrayLike) – Field points.
strength (float) – Source strength (volume flow rate per unit depth); negative for a sink.
x0 (float) – Source location.
y0 (float) – Source location.

Return type:

Returns:

ndarray of complex – Complex potential \(W(z) = \frac{m}{2\pi}\ln(z - z_0)\).

physicskit.fluids.stokes_drag(mu, radius, velocity)[source]#

Stokes’ law: drag force on a sphere at low Reynolds number.

\[F_D = 6\pi\mu R v\]

Valid in the creeping-flow limit \(Re = 2R|v|/\nu \ll 1\), where the nonlinear advection term in Navier-Stokes is negligible next to viscous diffusion, leaving the linear Stokes equations \(\mu\nabla^2\mathbf{u}=\nabla p\) – George Stokes’ 1851 exact solution for uniform flow past a sphere.

Parameters:

mu (float) – Dynamic viscosity.
radius (float) – Sphere radius.
velocity (float) – Relative speed between the sphere and the fluid far away.

Return type:

Returns:

float – Drag force magnitude \(F_D\).

Raises:

InvalidParameterError – If mu or radius is not positive.

Examples

>>> round(stokes_drag(mu=1.0, radius=1.0, velocity=1.0), 4)
18.8496

physicskit.fluids.strouhal_number(frequency, length, velocity)[source]#

The Strouhal number: a dimensionless vortex-shedding frequency.

\[St = \frac{f L}{U}\]

For flow past a bluff body shedding a von_karman_vortex_street() at frequency f, St stays remarkably constant (near 0.2 for a circular cylinder) over a wide range of Reynolds number, letting a single scaled model predict the shedding frequency at full scale.

Parameters:

frequency (float) – Vortex-shedding frequency f.
length (float) – Characteristic length scale L (e.g. the body’s cross-stream width).
velocity (float) – Characteristic flow speed U.

Return type:

Returns:

float – Strouhal number St.

Raises:

InvalidParameterError – If velocity is not positive.

Examples

>>> round(strouhal_number(frequency=2.0, length=0.1, velocity=1.0), 4)
0.2

physicskit.fluids.uniform_flow_potential(X, Y, U_inf, alpha=0.0)[source]#

Complex potential of a uniform stream.

Parameters:

X (ArrayLike) – Field points.
Y (ArrayLike) – Field points.
U_inf (float) – Free-stream speed.
alpha (float) – Angle of attack (radians), measured from the +x axis.

Return type:

Returns:

ndarray of complex – Complex potential \(W(z) = U_\infty e^{-i\alpha} z\).

physicskit.fluids.von_karman_vortex_street(n_pairs, spacing_l, spacing_ratio=0.2805)[source]#

Build a finite, doubly staggered von Karman vortex street.

Two parallel rows of point vortices, offset from each other by half the along-row spacing spacing_l and by a row separation h = spacing_ratio * spacing_l. Every vortex in the upper row has circulation -1 (clockwise) and every vortex in the lower row +1 (counterclockwise) – the staggered pattern shed alternately from the top and bottom of a bluff body in a stream toward \(+x\). Such an infinite street translates rigidly, relative to the surrounding fluid, at \(U = \frac{\Gamma}{2l}\tanh(\pi h/l)\) toward \(-x\) (von Karman 1911-1912; Lamb, Hydrodynamics, Sec. 156). Von Karman showed that a doubly infinite row of this form is linearly stable to vortex-row perturbations at exactly one spacing ratio, VON_KARMAN_SPACING_RATIO; any other ratio (or any symmetric, unstaggered arrangement) grows unstable, which is why the observed spacing behind real bluff bodies clusters so consistently near this value.

Parameters:

n_pairs (int) – Number of vortex pairs (one from each row) to construct; must be at least 2 to form a discernible street.
spacing_l (float) – Spacing between same-row, same-sign vortices along the street.
spacing_ratio (float) – Row separation over along-row spacing, h / spacing_l.

Return type:

Returns:

positions (ndarray of float, shape (2*n_pairs, 2)) – Vortex positions.
circulations (ndarray of float, shape (2*n_pairs,)) – Vortex circulations, interleaved upper/lower: -1 for the upper row (even indices), +1 for the lower row (odd indices).

Raises:

InvalidParameterError – If n_pairs is smaller than 2 or spacing_l is not positive.

Examples

>>> positions, circulations = von_karman_vortex_street(n_pairs=6, spacing_l=1.0)
>>> positions.shape
(12, 2)
>>> circulations[:4]  # upper row clockwise, lower row counterclockwise
array([-1.,  1., -1.,  1.])

physicskit.fluids.weber_number(rho, velocity, length, surface_tension)[source]#

The Weber number: inertial forces over surface-tension forces.

\[We = \frac{\rho U^2 L}{\sigma}\]

Governs whether a free liquid surface (a droplet, a jet, a bubble) deforms and breaks up under its own inertia (\(We \gg 1\)) or is held together by surface tension (\(We \lesssim 1\)) – the parameter that decides, for instance, whether a jet of liquid breaks into droplets.

Parameters:

rho (float) – Fluid density.
velocity (float) – Characteristic flow speed U.
length (float) – Characteristic length scale L (e.g. a droplet diameter).
surface_tension (float) – Surface tension coefficient \(\sigma\).

Return type: