Physics | Hydrodynamics » Michal Bejger - Introduction, equations of hydrodynamics

Introduction: equations of hydrodynamics Michał Bejger N. Copernicus Center, Warsaw 1 / 32 Outline ? Derivation of the equations: ? from the Boltzmann equation, ? by other methods, ? Regimes of applicability, ? some remarks about the advection equation. 2 / 32 What is a fluid? 3 / 32 The difference between fluid, and other phases, gas and solid: ? fluid flows, ? not resists deformation (not elastic), though may be viscous, ? forms free surface (gases don’t) 3 / 32 Fluid Source:description https://doksi.net ? system composed of many bodies, that can be described as a continuum, mean free path λ  fluid element size  system size In general λ = (nσ)−1 ? Forces between particles are short-range (saturation, screening). Some examples: Air water space λ 10−5 cm 10−9 cm 1015 cm 3 / 32 3 / 32 Euler Lagrange viewpoints Source: and https://doksi.net Eulerian approach: An observer looks at fluid motion from a specific location in space through which the

fluid flows as time passes. Lagrangian approach: An observer looks at fluid motion following an individual fluid parcel as it moves through space and time. 4 / 32 Vector operators: Some definitions 5 / 32 Vector operators: gradient Gradient represents a rate of change of a scalar field (derivative) in space. It is a ”vector” (technically a differential form): ∇f = ∂f ∂f e1 + · · · + en ∂x1 ∂xn 5 / 32 Vector operators: divergence Divergence is the volume density of the outward flux from an infinitesimal volume around a given point. The result is a scalar: For F = F i ei , ∇ · F = ∇i F i = ∂Fy ∂Fz ∂Fx ∂F i + + = ∂x i ∂x ∂y ∂z 5 / 32 Vector operators: curl Curl describes the infinitesimal rotation of a 3-dimensional vector field. The result is a vector: ∇ × F = e i ijk ∂j Fk       ∂Fy ∂Fz ∂Fx ∂Fx ∂Fy ∂Fz − i+ − j+ − k ∂y ∂z ∂z ∂x ∂x ∂y 5 / 32 Material/substantial derivative Fluid

properties are generally functions of position and time. We will denote ? ∂/∂t the rate of change w.rt time at some fixed position ? D/Dt the rate of change w.rt time while traveling with a fluid element. D Dt is called a material/substantial derivative. For some quantity f (vector field, scalar field etc.), it equals Df ∂f = + u · ∇f , Dt ∂t Why like that? Imagine we follow a change in f over a short time δt. Fluid element moved from r at t to r + uδt at t + δt Df f (r + uδt, t + δt) − f (r, t) = lim δt0 Dt δt 6 / 32 Continuity equation Mass conservation: the rate of change of fluid mass inside a volume equals to the net rate of fluid flow into the volume. Often expressed as ∂ρ + ∇ · (ρu) = 0 ∂t 7 / 32 Continuity equation Consider volume V enclosed by a surface S, with n normal vector. Total mass of fluid in V is Z Z Z d ρdV = − (ρu)·ndS ρdV and dt V S V (mass flux across S). Z d ρdV dt V Z = V fixed V ∂ρ dV ∂t Z = − Stokes

∇ · (ρu)dV V True for all V the continuity equation: ∂ρ + ∇ · (ρu) = 0 ∂t also written as Dρ + ρ∇ · u = 0 Dt 7 / 32 Continuity equation Example: advection equation from continuity equation. Advection operator u · ∇ = ui ∂ ∂x i We obtain the advection equation ∂ρ + u · ∇(ρ) = 0 ∂t from the continuity equation ∂ρ + ∇ · (ρu) = 0 ∂t by assuming ∇ · u = 0 (u solenoidal, incompressible, divergence-free vector; u ≡ ∇ × A) 7 / 32 Momentum equation (equation of motion) Momentum conservation: the rate of change of total fluid momentum in some volume equals to the sum of forces acting on the volume. For an inviscid fluid, it’s called the Euler equation: ∂u + u∇ · u = −∇P + ρf ∂t 8 / 32 Momentum equation (equation of motion) Consider volume V moving with the fluid (no flow through the boundary S). Fluid momentum is Z Z d ρudV ρudV . Its rate of change, dt V V equals to the force acting on V . Two types of forces

- body and surface: Z Z ρfdV and − PndS S | V {z } | {z } body force, e.g , gravity d dt Z surface force Z ρudV V (inviscid fluid, force normal to S) Z Z Du dV = ρfdV − PndS ρdV cons. V Dt V S Z Z = ρfdV − ∇PdV = Stokes ρ V V 8 / 32 Momentum equation (equation of motion) Consider volume V moving with the fluid (no flow through the boundary S). Fluid momentum is Z Z d ρudV ρudV . Its rate of change, dt V V equals to the force acting on V . Z Z Z Z Du d ρudV = ρ dV = ρfdV − PndS ρdV cons. V dt V Dt V S Z Z = ρfdV − ∇PdV Stokes V V We recover the Euler equation:   Du ∂u ρ =ρ + u∇ · u = −∇P + ρf Dt ∂t 8 / 32 Momentum equation (equation of motion) In case of viscosity, momentum equation is called the Navier-Stokes equation. Consider the i−th component of the surface force; Z Z i − Pn dS becomes σij nj dS S S with the stress tensor σij . For simple liquids and gases, with µ called the dynamical viscosity, and the second

viscosity −2/3µ:     1  ∂u  ∂uj 1   i σij = −Pδij + 2µ  + − (∇ · u) δij  , ∂xi 3  2 ∂xj  | {z } Strain rate tensor Z Z ∂ σij dV (Stokes again) S V ∂xj Assuming µ = const., we finally arrive at   Du 1 ρ = −∇P + ρf + µ ∇2 u + ∇(∇ · u) . Dt 3 σij nj dS = 8 / 32 Conserved quantities: energy Energy conservation: the rate of change of the total energy equals the rate at which work and other sources of energy act at the fluid element:  D ρ Dt    1 2 DU u + = −∇·(Pu)+ρu·f+ρ−∇·F. 2 Dt 9 / 32 Conserved quantities: energy Let’s multiply the momentum equation by u:   Du D 1 2 ρu · =ρ u = −u · ∇P + ρu · f Dt Dt 2 Rate of change of the kinetic energy equals the rate at which work is done by forces acting at the fluid element. Similarly, equation for the total energy, kinetic u2 /2 + internal U is  Z  Z Z d 1 2 u + U dV = ρu · fdV − u · (Pn)dS dt V 2 V S Z Z + ρdV −

F · ndS. | V {z } | S {z } generated at rate  heat flux across S 9 / 32 Conserved quantities: energy Let’s multiply the momentum equation by u:   Du D 1 2 ρu · =ρ u = −u · ∇P + ρu · f Dt Dt 2 Rate of change of the kinetic energy equals the rate at which work is done by forces acting at the fluid element. Similarly, equation for the total energy, kinetic u2 /2 + internal U is     DU D 1 2 ρ u + = −∇ · (Pu) + ρu · f + ρ − ∇ · F. Dt 2 Dt Equation for internal energy U only (subtracting the kinetic part): p Dρ 1 DU = 2 +  − ∇ · F. Dt ρ Dt ρ Since V = ρ−1 , we recover the first law of thermodynamics: dU = −PdV + δQ. |{z} added heat 9 / 32 Virial Source: theorem https://doksi.net Velocity is the rate of change of position following the fluid: u = Euler equation for gravitational field: ρ D 2r = −∇P+ρf Dt 2 = f=−∇ψ Dr . Dt −∇P−ρ∇ψ Let’s multiply by r and integrate over a volume V : Z Z Z D 2r ρdV = − r

· ∇PdV − r · ∇ψρdV . r· Dt 2 V V V Z Z Z  2 D 2r d Dr Dr LHS: r· ρdV = r · ρdV − ρdV 2 Dt dt V Dt Dt V V Z 1 d2 = |r|2 ρdV − 2T . 2 dt 2 V Z 1 Total kinetic energy T = ρu2 dV 2 V 10 / 32 Virial Source: theorem https://doksi.net Pressure term: ∇·r=3 Z − Z r·∇PdV = − V V z }| { Z Z Z ∇·(rP)dV + P∇ · rdV = − Pr · ndS +3 PdV V V } | S {z =0 Z Also, define total gravitational energy: − r·∇ψρdV = V Z and moment of inertia: I = 1 2 Z ψρ(r)dV = Ψ, V ρr2 dV , V to arrive at the scalar form of the virial theorem: Z 1 d 2I = 2T + 3 PdV + Ψ. 2 dt 2 V 10 / 32 Distribution function Particle’s distribution function is a function of positions, velocities and time, f (x, y , z, vx , vy , vz , t), which gives the number of particles per unit volume in single-particle phase space. Z Number density n(x, y , z, t) = f dvx dvy dvz Z Total number of particles N(t) = n dx dy dz. (x−µ)2 1 √ e − 2σ2 σ 2π Maxwell-Boltzmann

velocity probability distribution: !  m 3/2 m(vx2 + vy2 + vz2 ) f = exp − 2πkT 2kT Normal distribution: f (x, µ, σ) = describes particle speeds in idealized gases where the particles move freely inside a stationary container without interacting with one another, except for very brief collisions in which they exchange energy and momentum with each other. 11 / 32 Boltzmann’s Equation - Derivation For an ensemble of non-interacting point-like particles in momentum phase-space; Liouville 0 s theorem states that the distribution function f is conserved along the trajectories in phase space of positions q and momenta p: ∂f ∂f ∂qi ∂f ∂pi ∂f ∂f ∂f df = + + = + q̇i + ṗi = 0 dt ∂t ∂qi ∂t ∂pi ∂t ∂t ∂qi ∂pi ṗ i = F i is the “external” force, e.g, gravity, depending only on macroscopic quantities. Collision integral. In case of collisions, phase space evolves:   df ∂f ∂f ∂f ∂f = + q̇i + ṗi = 6= 0 dt ∂t ∂qi ∂pi ∂t col. 12 /

32 Boltzmann’s equation: collision integral The form of the collision integral can be approximated by assuming: ? ? ? ? Only short-range two-particle interactions should be important, Collisions are elastic, Collision time should be negligible. No correlations between the incoming and outgoing particles. This is Stosszahlansatz, (molecular chaos) assumption.  ∂f ∂t  ZZ = col. g dσ (g , Ω)[f (p0 A , t)f (p0 B , t)−f (pA , t)f (pB , t)] d Ω d 3 pA , dΩ where g = |pB − pA | = |p0 B − p0 A | and d σ/d Ω is the differential cross-section. In case of conserved quantities (mass, momentum, energy) the collision integral doesn’t provide a contribution to averaged quantities. 13 / 32 Momenta of the distribution function Z mfd 3 p Z Momentum ρu = mvfd 3 p Z m Internal energy ρ = (v − u)2 fd 3 p 2 Density ρ = (v = p/m, u is mean flow) 14 / 32 Continuity equation We start with df ∂f ∂f ∂f = + q̇i + ṗi = dt ∂t ∂qi ∂pi  ∂f

∂t  col. Let’s multiply by a particle mass m and integrate in momentum space: 15 / 32 Continuity equation df ∂f ∂f ∂f We start with = + q̇i + ṗi = dt ∂t ∂qi ∂pi  ∂f ∂t  col. Let’s multiply by a particle mass m and integrate in momentum space:   Z Z Z Z ∂f 3 ∂f ∂f 3 ∂f 3 m d p + mvi d p+ m ṗi d p = m d 3p ∂t ∂qi ∂pi ∂t col. | {z } =0, local mass cons. with ṗi = Fi , ∂ ∂ρ + ∂t ∂qi Z mfv i d 3 p − m | Z S n · FfdS = 0 {z } f 0 for V ∞ We recover the continuity equation for a mean fluid velocity, u i = hv i i of a fluid element: ∂ρ ∂(ρhv i i) + =0 ∂t ∂qi 15 / 32 Momentum equation Let’s multiply the Boltzmann equation by mv:   Z Z Z Z ∂f ∂f 3 ∂f 3 ∂f mv j d 3 p+ mv j v i d p+ mv j F i d p = mv j d 3p ∂t ∂qi ∂pi ∂t col. | {z } =0, local mom. cons Z ∂f 3 ∂(ρhv j i) ∂(ρhv i v j i) pj + + Fi d p=0 ∂t ∂qi ∂pi V Z Z Z Z ∂p j f 3 ∂p j 3 ∂f 3 d p= d p− fd p = −

fp i ni dS −ρδij pj ∂pi V ∂pi V ∂pi V | S {z } =0 j So we obtain, i j ∂(ρhv i) ∂(ρhv v i) + − ρF i δij = 0 ∂t ∂qi 16 / 32 Momentum equation ∂(ρhv j i) ∂(ρhv i v j i) + − ρF i δij = 0 ∂t ∂qi By subtracting the continuity equation multiplied by hv j i,   i ∂(ρhv j i) ∂(ρhv i v j i) j j ∂ρ j ∂(ρhv i) + − ρF − hv i + hv i =0 ∂t ∂qi ∂t ∂qi we obtain Navier-Stokes equation: ∂(hv j i) 1 ∂σ ij ∂(hv j i) + hv i i = Fj + , ∂t ∂qi ρ ∂qi with the stress tensor σ ij = −ρhv i v j i + ρhv i ihv j i. In case of collisions, σ ij = −Pδij + viscous stress tensor 16 / 32 Closure relations Since hydrodynamic equations are under-determined, one needs to close the set of equations with the stress tensor, a function of ρ, P and : ? Euler equations correspond to Maxwell distribution, σij = −Pδij (neglect viscous stresses and heat conduction), ? Chapman–Enskog theory: deviation from the

Maxwell–Boltzmann distribution in the equilibrium is small, and the first order corrections are     ∂uj 1 1 ∂ui − (∇ · u) δij , + σij = 2µ 2 ∂xj ∂xi 3 and the heat conduction term in the energy equation, F = κ∇T with shear viscosity µ ∝ √ mkT and thermal conductivity κ = 52 cV µ. 17 / 32 Dissipation in astrophysics An estimate for the effects of viscosity, heat conduction and other effects is provided by some dimensionless numbers: ? Knudsen number (are we in the fluid regime?): Kn = mean free path λ = size of the system l ? The Reynolds number ( dimensionless Navier-Stokes): Re = ρul ul inertial forces = = viscous forces µ ν ? Péclet number (how the heat is transported): Pe = advective transport rate ul = diffusive transport rate α ? Prandtl number (momentum-to-thermal diffusivity): Pr = viscous diffusion rate ν = thermal diffusion rate α (kinematic viscosity ν = µ/rho, thermal diffusivity α = k/(ρcP ) 18 / 32 Viscous

transport in disks Viscosity in action - what makes the accretion so efficient? α-viscosity prescription: ? ν ∝ αcs H, ? stress tensor in the disk, torque Tr φ ∝ ρνr ∂Ω ∂r ∝ −αP (at least for Shakura-Sunyaev disks). Perhaps magnetic field plays a role? Magnetorotational instability (Balbus-Hawley) 19 / 32 Simplifying approximations We have several simple flows at our disposals: ? Incompressible: ∇ · u = 0, equivalent to ρ = const. ? Anelastic: ∇ · (ρu) = 0, equivalent to ∂ρ/∂t = 0 (see continuity equation) ? Barotropic: P = P(ρ) instead of P = P(ρ, T . ), which is called baroclinic. It can result from the equation of state (degeneracy, eg for white dwarfs and neutron stars) or may be incidental, e.g when convection establishes s = const. throughout a star ? Adiabatic: heat transfer is neglected. Simplified energy relation d ρs/dt = 0 20 / 32 Inviscid and adiabatic flow When neglecting viscosity and heat conduction, the momentum and

energy equations are simplified. The Navier-Stokes equation reduces to the Euler equation:   Du 1 ρ = −∇P + ρf+µ ∇2 u + ∇(∇ · u) Dt 3  ρ D Dt  1 2 u 2  DU + Dt  = −∇ · (Pu) + ρu · f+ρ − ∇ · F The energy equation can be rewritten as a conservation of entropy: Dρs + ∇ · (ρsu) = 0. Dt 21 / 32 Laminar and Turbulent Flows ? Flow in which the kinetic energy decays due to the action of viscosity is called laminar flow. ? Large Re > 100-1000 turbulent flows (often seen in astrophysical settings). ? Problem: Resolution of turbulent flows down to the length scale where viscosity becomes important is not feasible (MRI etc.) 22 / 32 Plasma hydrodynamics ? Astrophysical fluids are quite often charged (ionized), so the neutrality condition is sometimes not satisfied, ? There is a limit for the mean free path λ: Coulomb scattering, Debye length (screening), s kB /e 2 P λD = ne /Te + i Zi2 ni /Ti for l ≥ λD fluid may be considered

neutral, ? If macroscopic process slower than the plasma frequency timescale (Langmuir waves, electron density fluctuations), s 4πne e 2 ωpe = me ? Effect of macroscopic magnetic fields (magnetospheres of pulsars) changes the mean free path (preferred directions, magnetic field gyration). 23 / 32 Magnetohydrodynamics In case of magnetic fields, even when the fluid is electrically neutral: ∂ρ + ∇ · (ρu) = 0, ∂t ∂u 1 2 1 + u∇ · u = −∇(P + B ) + ρf + j×B ∂t 8π 4π | {z } Lorentz force Maxwell equations: ∇ · B = 0, ∇ · E = ρe , 1 ∂E ∂B , ∇ × B = 4πj + 2 ∂t c ∂t + Equation for fluid energy evolution. ∇×E=− Ideal MHD approximation (primes denote the fluid co-moving quantities): ? j0 = σE0 σ(E + u × B) 1 ? ∂E ∂t 0 j = 4π ∇ × B ? Magnetic induction equation: magnetic diffusivity η. ∂B ∂t = ∇ × (u × B) + η∇2 B, with 24 / 32 Partial differential equations 25 / 32 Types of equations For linear

equations in two dimensions (x and t, say), one may classify the general one a ∂2F ∂2F ∂F ∂2F ∂F + c +e + fF + g = 0 + 2b +d 2 2 ∂x ∂x∂t ∂t ∂x ∂t For example: b2 < ac Elliptic b2 > ac Hyperbolic b2 = ac Parabolic Laplace equation Wave equation Diffusion ∂2F ∂x 2 ∂2F ∂x 2 ∂2F ∂x 2 − − ∂2F ∂t 2 =0 1 ∂2F c 2 ∂t 2 =0 + ∂F ∂t =0 25 / 32 Naming conventions Differences between elliptic and hyperbolic/parabolic equations: ? elliptic equations have boundary conditions which are specified around a closed boundary, ? derivatives are with respect to spatial variables (e.g, Laplace or Poisson), ? for hyperbolic - boundary conditions for time variable are initial conditions. 26 / 32 Hyperbolic equation - wave equation 2 ∂2y 2∂ y Consider a following equation, − c = 0. ∂t 2 ∂x 2    ∂ ∂ ∂ ∂ It can be rewritten as +c −c y = 0, ∂t ∂x ∂t ∂x or re-casted as a system of two first-order

equations, ∂z ∂z +c =0 ∂t ∂x ∂y ∂y −c =z ∂t ∂x This are examples of advection equations, ∂u ∂u +c =0 ∂t ∂x with boundary conditions u(x, t) = u0 (x) for t = t0 . 27 / 32 Advection equation: simple attempt Let’s try a centred difference for the space derivative (subscript j) and Euler’s method for the time derivative (superscript n): ∂u ∂u = −c ∂t ∂x ujn+1 − ujn δt = −c n n uj+1 − uj−1 2δx As expected, time part is 1-st order accurate: ujn+1 − ujn ≈ δt ∂u ∂t n j 1 ∂2u + δt 2 2 ∂t 2 n + . j Space part is 2-nd order: n n uj+1 − uj−1 ≈ 2δx ∂u ∂t n j 1 ∂3u + δx 3 3 ∂t 3 n + . j 28 / 32 Advection equation: simple attempt Substituting the Taylor expansions in the approximation: δt ∂u ∂t n j ∂2u 1 + δt 2 2 ∂t 2 n + ··· ≈ j cδt 2δx 2δx ∂u ∂t n j ∂3u 1 + δx 3 3 ∂t 3 n + . j - the truncation errors are 2nd order in time and 3rd order in space. How about

the stability in time? Test solution, plane wave vjn = v n exp(ikxj ) (von Neumann stability condition): v n+1 exp(ikxj ) = v n exp(ikxj ) − cδt n v (exp(ikxj+1 ) − exp(ikxj−1 )) , 2δx For small deviation, v n v n + δv n ,   cδt n+1 δv = 1−i sin(kδx) δv n δx | {z } growth factor S S = 1 − iα, so square of norm |S|2 = 1 + α2 ≤ 1 for all wave vectors k Problem. 28 / 32 Advection equation: Lax method Minor improvement: n n n n ujn+1 − (uj+1 − uj−1 )/2 uj+1 − uj−1 ∂u ∂u = −c = −c ∂t ∂x δt 2δx Let’s check the stability, again using the von Neumann condition:   cδt sin(kδx) δv n , δv n+1 = cos(kδx) − i δx | {z } growth factor S so the growth factor 2  2 !  cδt cδt sin2 (kδx) = 1 − sin2 (kδx) 1 − . |S|2 = cos2 (kδx) + δx δx This means stability for all k as long as δx ≥c δt Courant–Friedrichs–Lewy condition: the information propagation speed on the grid must be greater than all the physical speeds of

the problem. 29 / 32 Advection equation: two more examples Other schemes that result in the CFL condition: ? Upwind differencing (with backwards Euler in time): ujn+1 − ujn δt = −c n ujn − uj−1 δx ? Staggered leapfrog (centered differences in both space and time): ujn+1 − ujn−1 δt = −c n n uj+1 − uj−1 δx second order in time and space, but prone to the mesh drift instability; grid points of odd j + n and even j + n decoupled 30 / 32 Advection equation: two more examples Mesh drift instability 30 / 32 Advection equation with Lax-Wendroff Consider again an advection equation, with a following two-step method: n+1/2 First step is Lax method for uj+1/2 : n+1/2 n − ujn )/2 uj+1/2 − (uj+1 δt/2 = −c n uj+1 − ujn δx Second step, the quantities at tn+1 are calculated using the centered expression: n+1/2 n+1/2 uj+1/2 − uj−1/2 ujn+1 − ujn = −c δt δx n+1/2 n+1/2 (values uj+1/2 , uj−1/2 auxiliary; no mesh drift instability) 31

/ 32 Further reading. ? „An introduction to astrophysical fluid dynamics”, Michael J. Thompson ? „An Introduction to Astrophysical Hydrodynamics”, Steven N. Shore ? Amusing weblog: fuckyeahfluiddynamics.tumblrcom 32 / 32