physicskit.fluids #

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.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:

None

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:

None

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:

None

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.

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.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.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.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.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_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.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.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.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.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_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.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.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: