Fluid-Dynamics

The FluidExtensions class bridges VectorField (∇·, ∇×) and ScalarField (∇, ∇²) to classical fluid dynamics — Navier–Stokes equations, Bernoulli's principle, continuity, vorticity, drag/lift, dimensionless numbers, viscous pipe flow, and hydrostatics.

Namespace: CSharpNumerics.Physics.FluidDynamics

🌀 Navier–Stokes Equations (Incompressible)

$\rho\!\left(\frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v}\cdot\nabla)\mathbf{v}\right) = -\nabla p + \mu\nabla^2\mathbf{v} + \mathbf{f}, \qquad \nabla\cdot\mathbf{v} = 0$

Individual terms and the full residual:

var velocity = new VectorField(
    r => -r.y, r => r.x, r => 0);  // rigid-body rotation
var pressure = new ScalarField(
    r => 0.5 * 1000 * (r.x * r.x + r.y * r.y));
var point = (1.0, 0.0, 0.0);

// Convective acceleration: (v·∇)v
Vector conv = velocity.ConvectiveAcceleration(point);

// Viscous term: μ∇²v
Vector visc = velocity.ViscousTerm(dynamicViscosity: 1e-3, point);

// Pressure gradient force: −∇p
Vector gradP = pressure.PressureGradientForce(point);

// Full residual: ∂v/∂t = (−∇p + μ∇²v + f)/ρ − (v·∇)v
// For steady-state flow this should be ≈ 0
Vector residual = velocity.NavierStokesResidual(
    pressure, density: 1000, dynamicViscosity: 1e-3, point);

// With body force (e.g. gravity)
Vector residualG = velocity.NavierStokesResidual(
    pressure, density: 1000, dynamicViscosity: 1e-3, point,
    bodyForce: new Vector(0, 0, -9810));

// Incompressibility check: ∇·v = 0
double divV = velocity.IncompressibilityResidual(point);  // should be ≈ 0

🏛️ Euler Equations (Inviscid Flow)

Navier–Stokes with $\mu = 0$:

Vector euler = velocity.EulerEquationResidual(
    pressure, density: 1000, point);

// With body force
Vector eulerG = velocity.EulerEquationResidual(
    pressure, density: 1000, point,
    bodyForce: new Vector(0, 0, -9810));

⚖️ Bernoulli's Principle

$p + \tfrac{1}{2}\rho v^2 + \rho g h = \text{const}$

Along a streamline for steady, inviscid, incompressible flow:

// Bernoulli constant
double B = 101325.0.BernoulliConstant(density: 1000, velocity: 5.0, height: 10);

// Pressure at a second point
double p2 = 101325.0.BernoulliPressure(
    density: 1000, v1: 2.0, v2: 8.0, h1: 0, h2: 3.0);

// Dynamic pressure: q = ½ρv²
double q = 1.225.DynamicPressure(speed: 50);

// Stagnation (total) pressure: p₀ = p + ½ρv²
double p0 = 101325.0.StagnationPressure(density: 1.225, speed: 50);

🔄 Continuity Equation

$\frac{\partial\rho}{\partial t} + \nabla\cdot(\rho\mathbf{v}) = 0$

// Mass flux as a vector field: ρv (constant density)
VectorField massFlux = velocity.MassFlux(density: 1000);

// Mass flux with spatially varying density
var rhoField = new ScalarField(r => 1000 + 10 * r.x);
VectorField massFlux2 = velocity.MassFlux(rhoField);

// Continuity residual: ∇·(ρv) — should be 0 for incompressible flow
double cont = velocity.ContinuityResidual(density: 1000, point);

// Volume flow rate: Q = A·v
double Q = 0.01.VolumeFlowRate(speed: 3.0);   // 0.03 m³/s

// Speed from continuity: A₁v₁ = A₂v₂
double v2 = 0.01.ContinuitySpeed(speed1: 3.0, area2: 0.005);  // 6.0 m/s

🌪️ Vorticity & Circulation

$\boldsymbol{\omega} = \nabla\times\mathbf{v}$

// Vorticity: ω = ∇×v
Vector omega = velocity.Vorticity(point);

// Enstrophy density: ½|ω|²
double enstrophy = velocity.Enstrophy(point);

// Helicity density: v·ω
double helicity = velocity.HelicityDensity(point);

🌊 Stream Function (2-D)

For incompressible 2-D flow, $v_x = \partial\psi/\partial y$, $v_y = -\partial\psi/\partial x$:

// Create velocity field from a stream function
var psi = new ScalarField(r => r.x * r.y);  // ψ = xy → saddle flow
VectorField v = psi.VelocityFromStreamFunction();
// vx = ∂ψ/∂y = x,  vy = −∂ψ/∂x = −y → strain flow

// The resulting field is automatically divergence-free
double div = v.Divergence((1, 1, 0));  // ≈ 0

🪂 Drag & Lift

// Drag force: F_D = ½ρv²C_D A
double drag = 0.47.DragForce(density: 1.225, speed: 30, referenceArea: 0.01);

// Lift force: F_L = ½ρv²C_L A
double lift = 1.2.LiftForce(density: 1.225, speed: 60, referenceArea: 20);

// Terminal velocity: v_t = √(2mg / (ρC_D A))
double vTerm = 0.5.TerminalVelocity(
    dragCoefficient: 0.47, density: 1.225, referenceArea: 0.01);

🔢 Dimensionless Numbers

// Reynolds number: Re = ρvL/μ
double Re = 1000.0.ReynoldsNumber(speed: 2.0, characteristicLength: 0.1,
    dynamicViscosity: 1e-3);  // Re = 200 000

// Mach number: Ma = v/c (defaults to speed of sound in air)
double Ma = 300.0.MachNumber();           // Ma ≈ 0.875
double Ma2 = 1500.0.MachNumber(1480);    // in water

// Froude number: Fr = v/√(gL)
double Fr = 5.0.FroudeNumber(characteristicLength: 2.0);

// Strouhal number: St = fL/v
double St = 50.0.StrouhalNumber(characteristicLength: 0.1, speed: 25);

// Weber number: We = ρv²L/σ
double We = 1000.0.WeberNumber(speed: 1.0, characteristicLength: 0.001,
    surfaceTension: 0.072);

🚰 Viscous / Pipe Flow

Hagen–Poiseuille flow in circular pipes:

// Volume flow rate: Q = πR⁴ΔP / (8μL)
double Q = 0.005.PoiseuilleFlowRate(
    pressureDrop: 1000, dynamicViscosity: 1e-3, length: 1.0);

// Velocity profile: v(r) = (ΔP/4μL)(R² − r²)
double vCenter = 0.005.PoiseuilleVelocity(
    pressureDrop: 1000, dynamicViscosity: 1e-3, length: 1.0,
    radialDistance: 0);    // maximum at centre

double vWall = 0.005.PoiseuilleVelocity(
    pressureDrop: 1000, dynamicViscosity: 1e-3, length: 1.0,
    radialDistance: 0.005);  // zero at wall

Stokes drag on a sphere in creeping flow:

// F = 6πμRv
double F = 1e-3.StokesDrag(radius: 0.01, speed: 0.1);

🏊 Hydrostatics

// Hydrostatic pressure at depth: p = p₀ + ρgh
double p = 10.0.HydrostaticPressure(density: 1025);  // 10 m depth in seawater

// Buoyant force (Archimedes): F_b = ρ_fluid · V · g
double Fb = 1025.0.BuoyantForce(displacedVolume: 0.5);  // ≈ 5025 N

⚡ Fluid Energy & Momentum

var v = new Vector(3, 4, 0);

// Kinetic energy density: e_k = ½ρ|v|²
double ek = v.KineticEnergyDensity(density: 1000);  // 12 500 J/m³

// Momentum density: ρv
Vector mom = v.MomentumDensity(density: 1000);  // (3000, 4000, 0) kg/(m²·s)

🌀 Kármán Vortex Street

The KarmanVortexStreetExtensions class provides formulas for analysing periodic vortex shedding behind bluff bodies (cylinders). It covers the Roshko empirical Strouhal–Reynolds correlation, vortex street geometry (von Kármán stability ratio), regime classification, and the idealised wake-drag formula.

Formula	Method	Description
$f = \mathrm{St} \cdot U / D$	VortexSheddingFrequency	Shedding frequency from Strouhal number
$\mathrm{St} \approx 0.198(1 - 19.7/\mathrm{Re})$	RoshkoStrouhalNumber	Roshko (1954) empirical correlation
$f = \mathrm{St}(\mathrm{Re}) \cdot U / D$	VortexSheddingFrequencyFromReynolds	Shedding frequency from Reynolds number
$\mathrm{Ro} = f D^2 / \nu$	RoshkoNumber	Roshko number (= St · Re)
$a = U / f$	VortexStreamwiseSpacing	Longitudinal spacing between vortices
$h/a = \tfrac{1}{\pi}\mathrm{arccosh}(\sqrt{2})$	StableSpacingRatio	Von Kármán stability criterion ≈ 0.281
$h = (h/a) \cdot a$	VortexLateralSpacing	Cross-stream spacing between rows
$U_v \approx 0.875\,U$	VortexConvectionVelocity	Downstream vortex convection speed
—	IsPeriodicSheddingRegime	True if 47 < Re < 2×10⁵
—	CylinderWakeRegime	Wake regime string from Re
$C_d$ (wake momentum)	VortexStreetDragCoefficient	Idealised wake-drag from vortex spacing

using CSharpNumerics.Physics.FluidDynamics;

// Roshko's Strouhal–Reynolds correlation for a circular cylinder
double Re = 1000;
double St = Re.RoshkoStrouhalNumber();           // ≈ 0.194

// Vortex shedding frequency
double U = 5.0, D = 0.05;
double f = St.VortexSheddingFrequency(U, D);     // ≈ 19.4 Hz

// Or directly from Reynolds number
double f2 = Re.VortexSheddingFrequencyFromReynolds(U, D);

// Roshko number Ro = f·D²/ν
double nu = 1e-5;
double Ro = f.RoshkoNumber(D, nu);

// Vortex street geometry
double a = U.VortexStreamwiseSpacing(f);         // streamwise spacing
double h = a.VortexLateralSpacing();             // lateral spacing ≈ 0.281·a

// Vortex convection velocity
double Uv = U.VortexConvectionVelocity();        // ≈ 0.875·U

// Regime classification
bool shedding = Re.IsPeriodicSheddingRegime();   // true
string regime = Re.CylinderWakeRegime();         // "Turbulent transition"

// Wake-drag coefficient from vortex street
double Cd = h.VortexStreetDragCoefficient(a, U, Uv, D);

🔬 Viscous Flow Model (IViscousFlowModel)

The IViscousFlowModel interface provides an abstraction for pipe flow physics used by simulation engines. The default implementation ViscousFlowModel encapsulates Hagen–Poiseuille physics.

Method	Description
DrivingForce(dP/dx, ρ)	Body force per unit mass $f = -(dP/dx) / \rho$
CylindricalDiffusion(ν, d²v/dr², dv/dr, r)	$\nu [d^2v/dr^2 + (1/r)(dv/dr)]$
SymmetryAxisDiffusion(ν, d²v/dr²)	$2\nu \, d^2v/dr^2$ (L'Hôpital at $r=0$)

using CSharpNumerics.Physics.FluidDynamics;
using CSharpNumerics.Physics.FluidDynamics.Interfaces;

IViscousFlowModel flow = new ViscousFlowModel();

double driving = flow.DrivingForce(pressureGradient: -100, density: 1000);
// 0.1 m/s² driving force

double diff = flow.CylindricalDiffusion(nu: 1e-6, d2vdr2: 50, dvdr: 10, r: 0.01);