Hydrodynamics Rami Vainio Spring 2010 Hydrodynamics – Spring 2010 53376 Hydrodynamics 10 credit points Lectures: Mon 10–12 and Fri 12–14 in Physicum D117 Lecturer: Rami Vainio (phone: 191 50676; email: rami.vainio@helsinki.fi) Exercises: Thu 14–16 in Exactum C220. ∼4 problems per week. Submit to the box in 2nd floor A-wing lobby by Tue, 6 pm. Assistant: Heli Hietala (heli.hietala@helsinki.fi) Student seminar: Each student has to give a seminar presentation about a topic related to the course (20+5 min) write an essay (of 5–8 pages including figures) on it. The list of topics with links to additional information will become available at the course website on the third week of the course; the students should reserve their topic from the lecturer by Fri, 19.02.2010. Topics are given out on “first come first served” basis. If you want to propose a topic of your own, contact the lecturer. Grading: final exam 50% + exercises 25% + essay 15% + seminar presentation 10%. Hydrodynamics – Spring 2010 1 Course outline The course will cover the basic theory of hydrodynamics and present a number of examples especially related to astrophysics. However, the scope of the course is general and aimed at atmospheric physicists and geophysicicts alike. We will follow the textbook (Ch. 1–9 in full, 14–16 selectively) A.R. Chourdhuri: The Physics of Fluids and Plasmas – An Introduction for Astrophysicists. Cambridge University Press, 1998 Additional material on non-astrophysical applications will be added, where appropriate. Preliminary Outline 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. Introduction Boltzmann equation Derivation of the hydrodynamic equations Ideal fluids Viscous flows Gas dynamics Linear waves and instabilities Turbulence Rotation in hydrodynamics Magnetohydrodynamics Hydrodynamics – Spring 2010 2 Introduction Hydrodynamics – Spring 2010 3 Scope of hydrodynamics Hydrodynamics = theory of fluid dynamics Fluids (fluidit): • Liquids (nesteet) • Gases (kaasut) Usually plasmas (ionized gases) are treated separately from neutral fluids. However, • hydrodynamics is a valid macroscopic theory for weakly magnetized plasmas as well. • magnetohydrodynamics (MHD) covers the macroscopic evolution of magnetized plasmas. • microscopic theories of neutral gases and plasmas are very different, as charged particles interact via long-range forces. On this course we will not address the microscopic treatment of plasmas but refer to the course Plasma physics (5 cp, fall term) as the appropriate introduction to plasmas. Hydrodynamics is a continuum theory: material treated at a macroscopic level without regard to the motion of individual particles. Intuitively clear : liquids (like water) and gases (like air) under normal conditions can be treated as continua. But what about dilute gases like the solar wind and the interstellar medium (a few particles per cc)? What are the necessary assumptions for the continuum hypothesis to be valid? To shed light on that, we will derive the hydrodynamic equations from a more general theory. Hydrodynamics – Spring 2010 4 Dynamical theories Our aim: develop dynamical theory for fluids. Dynamical theory: physical theory allowing the study of time evolution of a system. Familiar examples: Classical mechanics, classical electrodynamics, quantum mechanics. For a dynamical theory, we need a way to describe the state of the system. E.g., • classical mechanics: generalized coordinates qi and momenta pi ∀i = 1, ..., N • classical electrodynamics: electric and magnetic field E(x) and B(x) ∀x • quantum mechanics: wave function ψ(x) ∀x Thus, the state is prescribed by giving the numerical values of a set of variables. We also need a set of equations which tells us how the variables change with time. Once we know the state at some instant, we shall be able to calculate the state at any other instant. E.g., • classical mechanics: Hamilton’s (or Lagrange’s or Newton’s) equations • classical electrodynamics: Maxwell’s equations • quantum mechanics: time-dependent Schrödinger equation Thus: what are the appropriate state variables in fluid mechanics and how do they evolve in time? The answer depends on the level of theory we need to apply. Hydrodynamics – Spring 2010 5 Different levels of dynamical theories Looking at a fluid consisting of N particles (with log N ≫ 1), we may study the system at different levels of theory. Level 0: N quantum particles. State described by ψ(x1, x2, ..., xN ), time-evolution by ∂ψ = Ĥψ; i~ ∂t N ∂2 ~2 X Ĥ = − + V ({xi}) 2m i=1 ∂x2i Level 1: N classical particles. State described by {qs, ps}, time-evolution by              q̇s = ∂H ∂ps             ∂H ṗs = − ∂qs H({qs, ps}, t) = X r q̇r pr − L; L({qs, q̇s}, t) = T − V ; ps = ∂L ∂ q̇s Level 2: Kinetic description. State described by f (x, u), time-evolution by ∂f ∂f + ẋ · ∇f + u̇ · = C(f ); ∂t ∂u ẋ = u; u̇ = F (x, t) m Level 3: Continuum description. State described by, e.g., ρ(x), T (x) and v(x), time-evolution by the hydrodynamic equations. Hydrodynamics – Spring 2010 6 Level 0 → Level 1 Is it always possible to replace the quantum description by the classical description? Only if particle density is low enough so that quantum interference is unimportant. Size of the wave packet ∼ de Broglie wavelength = λ = h/p r √ Typical particle momentum p ∼ m κBT /m = mκBT Interparticle distance ∼ n−1/3 Thus, λ ≪ n−1/3 if hn1/3 n[cm−3]1/3 −7 1≫√ = 4.36 · 10 r mκBT m[amu] T [K] Clearly, the condition is well satisfied for the dilute gases in space. How about denser fluids? For air in standard temperature (273 K) and pressure (1013 hPa), the righthand side is about 0.015. For water at 20◦C it is about 0.19, so the condition is met reasonably well. However, with liquids you need to consider the molecular size as well. If the condition is met, then the wave packets of the particles evolve according to their respective Schrödinger equations in an isolated manner, i~∂ψi(xi, t)/∂t = Ĥiψi, Ĥi = −(~2/2m)∇2i + V̄ (xi) and the expectation values hxii(t) and hpii(t) evolve according to the classical equations of motion (Ehrenfest’s theorem). Hydrodynamics – Spring 2010 7 Level 1 → Level 2 Level 1: prescribe the system of N classical particles with their positions,(x1, ..., xN ), and velocities, (u1, ..., uN ). Evolution governed by Newton’s laws. Phase space Γ of the system composed of 6N dimensions, (x1, ..., xN , u1, ..., uN ) giving the coordinates. Each point in Γ corresponds to a different state of the system, evolution of the system traced by a curve in Γ. For large N (∼ NA), solving 6N first-order (or 3N second order) equations of motion not possible in practice. → Statistical description required. Level 2: introduce the distribution function f (x, u, t) = density of particles in a six-dimensional “µ-space” (x, u). Single point in Γ is mapped to N points in µ (each particle corresponding to one point given by its position and velocity coordinates). Mapping illustrated below by an analogue of three particles in a 6-D Γ-space (x1, x2, x3, u1, u2, u3) and 2-D µ-space (x, u) (i.e., particles in 1-D motion). u2 x2 Instead of following the motion of all the N points, we follow the evolution of the distribution function. It is given by the Boltzmann equation. Γ µ u u1 u3 u1 u2 u3 x1 x3 x1 x2 x x3 Hydrodynamics – Spring 2010 8 Level 2 → Level 3 A single-component gas in thermodynamic equilibrium fully described by two thermodynamic variables. A fluid in motion is not in TE as a whole. However, a small element of fluid in its own rest frame can be in approximate TE . → State of the fluid element prescribed by two thermodynamic variables (e.g., density ρ and temperature T ) and velocity v . Prescribe all fluid elements → continuum: ρ(x), T (x) and v(x) prescribe the whole system. [Other choices often used as well, e.g., ρ(x), pressure p(x) and v(x).] Hydrodynamic equations govern the time-evolution of these variables as a function of position and, hence, consitute a dynamical theory of a fluid as continuum. They can be derived from the Boltzmann equation. Note: Each number needed to prescribe the system corresponds to one dimension of phase space of the system. As the five variables ρ, T and vi prescribing the system need to be given in all points of configuration space, the phase space of the hydrodynamic system has actually an infinite number of dimensions. Hydrodynamics – Spring 2010 9 Hamiltonian description of Level 1 Let us study the transition “Level 1 → Level 2” in more detail. Before that, let us recall the Hamiltonian formulation of classical mechanics. Consider a dynamical system prescribed by (qs, ps; s = 1, ..., n) and governed by the Hamiltonian H(qs, ps, t). Thus, ∂H ṗs = − ∂qs q̇s = ∂H . ∂ps For example, for a conservative system of N classical particles H = T + V = (2m)−1 X s p2s + V (qs), where qs are the (n = 3N ) Cartesian spatial coordinates of the particles, so q̇s = ∂H = ps/m ⇒ ps = mq̇s = component of mechanical momentum ∂ps ṗs = − ∂H ∂V =− ∂qs ∂qs Newtonian equation of motion Note: generalised coordinates qs do not have to have dimensions of length. If they don’t, dimension of ps are not those of mechanical momentum. (Always, however, [qs][ps] = [h]). Hydrodynamics – Spring 2010 10 Ensembles in phase space Ensemble = a set of many replicas of the same system, which are identical in all other respects apart from being in different states at an instant of time. ⇒ Each member of an ensemble can be represented by a point in phase space Γ at an instant of time and their evolution corresponds to different trajectories in Γ. Ass. Ensemble points are distributed densely and smoothly in Γ ⇒ A density of ensemble points, ρens(qs, ps, t), can be defined (at any location in Γ). Example. Consider a set of N non-interacting identical 1-D harmonic oscillators. The system has a Hamiltonian N 1 2 mω02 2 X Hs, where Hs(qs, ps) := H= p + q = Es. 2m s 2 s s=1 As the oscillators do not interact, each one conserves its total energy Es. The microstate of the system, described by the phases of the oscillators, is unknown. A statistically representative ensemble of this system, therefore, has the density of ensemble points given by N ω0 !N Y δ(Hs(qs, ps) − Es), ρens(qs, ps, t) = 2π s=1 n n normalized so that s d qs d ps ρens = 1, i.e., interpreted as a probability density for finding the system in a given microstate (qs, ps). ´ Q Hydrodynamics – Spring 2010 11 Liouville’s theorem Liouville’s theorem: Dρens = 0, Dt where D/Dt denotes a time-derivative along the path in phase space of any member of the ensemble. Proof : Dρens ρens(qs + δqs, ps + δps, t + δt) − ρens(qs, ps, t) = lim , δt→0 Dt δt where (qs, ps) and (qs + δqs, ps + δps) represent points occupied by one member of the ensemble on the phase-space trajectory of the system at t and t + δt, respectively. Taylor-expand: ρens(qs + δqs, ps + δps, t + δt) = ρens(qs, ps, t) + X s δqs ∂ρens ∂ρens X ∂ρens + δps + δt s ∂qs ∂ps ∂t to get Dρens ∂ρens X ∂ρens X ∂ρens = + q˙s + p˙s s s Dt ∂t ∂qs ∂ps   ∂ρens X ∂ ∂ ṗs  X ∂ X ∂ q̇s (q˙sρens) + (p˙sρens) − ρens  + + = s ∂ps s ∂qs s ∂qs ∂t ∂ps   ∂ρens X ∂ ∂ ∂H ∂ ∂H  X ∂ X = (q˙sρens) + (p˙sρens) − ρens  − + s ∂ps s ∂qs s ∂qs ∂ps ∂t ∂ps ∂qs } | {z =0 Hydrodynamics – Spring 2010 12 Since the number of members in the ensemble is conserved we must have, for an arbitrary phase space volume VΓ, ˆ ˆ ∂ ρensv Γ · ds, ρensdV = − ∂t VΓ ∂VΓ where v Γ = (q̇s, ṗs) is the 2n-dimensional velocity vector in Γ space and ds is a surface element (with outward normal) of the 2n − 1 dimensional hypersurface ∂VΓ bounding VΓ. p vΓ ds q The surface integral can be converted to volume integral, so ˆ ˆ   ∂ρens ∂ ∂ X  dV = − (q̇sρens) + (ṗsρens) dV ∂ps VΓ ∂t VΓ s ∂qs and as the volume is arbitrary, we have   ∂ρens X  ∂ Dρens ∂ + . 0= (q̇sρens) + (ṗsρens) = s ∂qs ∂t ∂ps Dt Hydrodynamics – Spring 2010  13 Conservation of phase space volume Ass. Ensemble points lie initially the phase space volume dnqs dnps. p Ass. As the system evolves in time, the ensemble points move to another volume dnqs′ dnp′s. The number of points is conserved, i.e., ρensdnqs dnps = ρ′ensdnqs′ dnp′s and since according to Liouville’s theorem ρens = ρ′ens, we get dnqs dnps = dnqs′ dnp′s, q i.e., the phase space volume occupied by the points is conserved. Hydrodynamics – Spring 2010 14 Distribution function in µ-space Let us now return to the statistical description (level 2) of the fluid consisting of N similar classical particles. We noted that • the state of the system is represented by a point in the 6N -D phase space Γ and its timeevolution there is traced by a curve. • the state of the system can be mapped to a 6-D µ-space (x, u), where it is represented by N points with coordinates (xi, ui) and its evolution traced by N curves. The distribution function is defined as f (x, u, t) = lim + δV →0 f δN δV δN , δV where δN is the number of particles in a small µspace volume δV around the point (x, u). The limit δV → 0+ means that δV is very small compared to the overall volume of the fluid but still large enough so that a large number of particles is inside the volume. Only if this limit exists, can we pass from Level 1 to Level 2. Hydrodynamics – Spring 2010 . . u u+du . u . . . . .. . . .. . . ... . . . . . . . . . . .. . . . . . . . . .. . .. . . . . .. . . . ... . . . .. . . .. . .. . . . . . x+dx x . . . . . . . . . . . . . ... .. . . . . . . . . . . . . δV,. δN . . . .. . . . . . . . . . . .. . . . . . . . . . . x . 15 Collisionless Boltzmann equation Our N particles admit Hamiltonian description. Can we obtain an equation similar to Liouville’s equation for the distribution function f in µ-space? If the N particles are each governed by a Hamiltonian of form H = H(x, p, t) with p = mu, then mu̇i = − ∂H ; ∂xi mẋi = ∂H ; ∂ui i = 1, 2, 3 and Df ∂f ∂f ∂f + u̇i = 0. (summation over i = 1, 2, 3implied!) = + ẋi Dt ∂t ∂xi ∂ui This case corresponds to N non-interacting particles moving in an external potential and the equation is termed collisionless Boltzmann equation. If particles of the fluid get in close distances so that they interact, this treatment does not apply. Then, the potential energy of each particle is not only a function of the particle position but also of the distance to the neighboring particles. In that case, the Hamiltonian formulation in µ-space becomes impossible. (No problems arise, of course, in the Hamiltonian formulation in Γ-space.) In neutral gases, however, the treatment of interactions is relatively easy, as the particles interact only when at close distances, i.e., only when they collide. The effect of the collisions can be added on the RHS of the equation. → Boltzmann equation. Hydrodynamics – Spring 2010 16 Boltzmann equation Hydrodynamics – Spring 2010 17 Collisions in a dilute gas Collisionless Boltzmann equation satisfied by f (x, u, t) if interactions between particles are neglected. Consider a dilute fluid , in which the total volume of the fluid particles is much smaller than the volume occupied by the fluid, na3 ≪ 1; n = number density, a = radius of a particle. Neutral particles → no long-range forces → interactions between the particles assumed to occur only when they collide, i.e., when the separation between the particles is not much larger than d = 2a. A particle moves freely (in a straight line) between two collisions. The average distance between two collisions is the mean free path, 1 . λ=√ 2 nπd2 Clearly, in a dilute fluid, λ ≫ d. Thus, binary collisions are much more common than collisions with three or more particle interacting at the same time. Also, particle trajectories become zig-zag paths. Collisions produce changes in f (x, u, t): 1. some particle with initial velocity u may have other velocities after collisions ⇒ δf (u) < 0 2. some particles with other initial velocities may have velocity u after collisions ⇒ δf (u) > 0 Thus, Hydrodynamics – Spring 2010 Df 3 3 d x d u = −Cout + Cin Dt 18 Binary collisions Consider two particles with the same mass colliding. Ass. Initial velocities of the particles are u and u1 and final velocities are u′ and u′1. Conservation of momentum and energy: u u1 u + u1 = u′ + u′1 2 1 2 |u| + 12 |u1|2 = 21 |u′|2 + 12 |u′1|2 In addition, assuming a central force field between the particles, the relative velocity after the collision, u′rel = u′ −u′1, lies in the plane of the initial relative velocity, urel = u − u1, and initial radius vector between the particles, r = x − x1. u’1 u’ Five conditions, six unknowns (primed vectors). The sixth condition has to come from the nature of interaction between the particles. A statistical description of collisions sufficient. Thus, the sixth condition comes from the prescription of the differential scattering cross-section. Hydrodynamics – Spring 2010 19 Differential scattering cross-section Consider two beams of similar particles colliding; densities n1 and n and velocities u1 and u, respectively. Particles in the second beam experience a flux of I = |u − u1|n1 of particles from the first beam. What is the number of collisions, δnc, per unit time and unit volume, which deflect the particles from the second beam to the solid angle dΩ? u1 u’ u dΩ δnc = σ(u, u1|u′, u′1) n |u − u1|n1 dΩ, u’1 where σ(u, u1|u′, u′1) is the differential scattering cross section. Conservation of momentum and energy and the requirement that the scattering is in the solid angle dΩ determine the (primed) final velocities. If molecular processes are assumed reversible, we can write for the reverse cross-section σ(u′, u′1|u, u1) = σ(u, u1|u′, u′1). Hydrodynamics – Spring 2010 20 Collision integral Consider a beam with velocities in d3u around u. Its number density is n = f (x, u, t) d3u. Similarly, a beam with velocities in d3u1 around u1 has a number density of n1 = f (x, u1, t) d3u1. Thus, δnc = σ(u, u1|u′, u′1)|u − u1| f (x, u, t) f (x, u1, t) dΩ d3u d3u1 gives the number of particles scattered to the solid angle dΩ (around u′) per unit volume and unit time. Since Cout is the total number of collisions per unit time in µ-space volume d3x d3u, it corresponds to ˆ ˆ Cout = d3x d3u d3u1 dΩ σ(u, u1|u′, u′1)|u − u1| f (x, u, t) f (x, u1, t). For Cin, we need to consider reverse collisions between particles with velocities in d3u′ and d3u′1, which scatter into d3u and d3u1, respectively. Analogously, δn′c = σ(u′, u′1|u, u1)|u′ − u′1| f (x, u′, t) f (x, u′1, t) dΩ d3u′ d3u′1. Conservation laws ⇒ |u′ − u′1| = |u − u1|. Volume occupied by the particles conserved in scattering; thus, conservation of 2-particle phase space volume implies d3u′ d3u′1 = d3u d3u1, provided that the scattering can be described by a potential allowing a Hamiltonian treatment. Hydrodynamics – Spring 2010 21 Thus, assuming microscopic reversibility, δn′c = σ(u, u1|u′, u′1)|u − u1| f (x, u′, t) f (x, u′1, t) dΩ d3u d3u1 and Cin = d3u d3x Recall that ˆ d3u1 ˆ dΩ σ(u, u1|u′, u′1)|u − u1| f (x, u′, t) f (x, u′1, t). Df 3 3 d x d u = −Cout + Cin Dt Hence, we finally have ∂f F ∂f + u · ∇f + · = ∂t m ∂u ˆ d3u1 ˆ dΩ σ(Ω)|u − u1| (f ′f1′ − f f1) ≡ C[f ], where f = f (x, u, t); f1 = f (x, u1, t); f ′ = f (x, u′, t); f1′ = f (x, u′1, t), ẋ = u and F = mu̇ have been substituted. Here, F includes the external forces (like gravity) and not the interparticle forces modeled by collisions. Here, it has also been assumed that σ = σ(Ω), which holds if the differential cross section only depends on the angle between the relative particle velocities u − u1 and u′ − u′1. Note that when performing the integrations over Ω (specifying the direction of u′) and u1, one needs to use the conservation laws to give u′ and u′1 as functions of u (fixed argument) and u1 (integration variable). Boltzmann equation is a non-linear integro-differential equation for the distribution function. Hydrodynamics – Spring 2010 22 Evaluation of the collision integral Collision integral contains an integral over u1 measured in the laboratory frame. The integral is easier to evaluate in a coordinate system centered at u. Let us denote urel ≡ u1 − u. Thus, ˆ ˆ C[f ](x, u, t) = d3urel dΩ σ(θ)|urel|[f (u + u′)f (u + u′1) − f (u)f (u + urel)], where dΩ = sin θ dθ dϕ is measured in the same frame, so that the differential cross-section becomes a function of scattering angle θ, only. In this frame, the conservation laws read urel = u′ + u′1 2 1 2 |urel | = 12 |u′|2 + 12 |u′1|2 and, therefore, 21 |urel|2 = 12 |u′ +u′1|2 = 21 |u′|2 +u′ ·u′1 + 21 |u′1|2 ⇒ u′ ·u′1 = 0, so unless u′1 = 0 (head-on collision), the angle between the particle trajectories is 90◦. Thus, also |u′rel| = |u′1 −u′| = |urel| and u′ = |urel| cos θ(e3 cos θ + e⊥ sin θ); u′1 = |urel| sin θ(e3 sin θ − e⊥ cos θ), where e3 k urel, e⊥ = e1 cos ϕ + e2 sin ϕ and we take, e.g., e1 k urel × u, e2 = e3 × e1. u’ θ urel u’1 Hydrodynamics – Spring 2010 e1 u’ u’1 e ϕ u’1 u’ e 2 e3 23 Note, finally, that π − 2θ = (urel, u′rel) = θCM so cos θCM = − cos 2θ. Since σ(θ) dΩ = σCM(θCM) dΩCM, σ can also be evaluated in the center-of-mass (CM) frame, where it is easiest to calculate. As urel k u measured in the CM frame, ϕ remains the same in the two coordinate systems. Thus, σ(θ) sin θ dθ = σCM(θCM) sin θCMdθCM = σCM(θCM) d cos θCM = σCM(θCM) d cos 2θ = σCM(θCM) 2 sin 2θ dθ = 4 cos θ σCM(θCM) sin θ dθ giving σ(θ) = 4 cos θ σCM(π − 2θ). Clearly, the integral over Ω extends over the forward half-space, only, and this gets mapped to the full 4π sr in the CM frame. In summary: C[f ] = ˆ 3 d urel ˆ π/2 0 dθ ˆ 0 2π dϕ 2 sin 2θ σCM(π − 2θ)|urel|[f (u + u′)f (u + u′1) − f (u)f (u + urel)], u′ = |urel| cos θ(e3 cos θ + e⊥ sin θ); e3 k urel; u′1 = |urel| sin θ(e3 sin θ − e⊥ cos θ), e⊥ = e1 cos ϕ + e2 sin ϕ; e1 k urel × u, e2 = e3 × e1 Hydrodynamics – Spring 2010 24 Maxwellian distribution Maxwellian distribution was first derived by Maxwell using the isotropy argument. ´ +∞ Consider velocities in the x-direction, distributed according to F (ux) normalized as −∞ F (ux) dux = 1. However, as there’s nothing special about the x-direction, velocities in the y - and z -direction are also distributed according to the same function. Thus, the probablity of finding a particle in duxduy duz is 1 f (u) duxduy duz = F (ux) dux · F (uy ) duy · F (uz ) duz n However, if f is an isotropic distribution, it should only depend on the magnitude u (or, e.g., its square). Thus, 1 f (u) := f¯(u2x + u2y + u2z ) = F (ux)F (uy )F (uz ). n Thus, a sum of u2i ’s in the argument of f¯ has to correspond to a product of F (ui)’s and this is possible only if F (ux) is an exponential function of u2x, i.e., 2 2 2 2 2 F (ux) = A1/3e−Bux ⇒ f (u) = nAe−B(ux+uy +uz ) = nAe−Bu . Constant A (as a function of B > 0) can be found by normalizing the F -distribution to unity. To find out the value of B , we need to apply the average kinetic energy of the particles, h 21 mu2i = 23 κBT , which then gives    3/2 2 mu m   −  exp  f (u) = n  2πκBT 2κBT Hydrodynamics – Spring 2010 25 Maxwellian as a solution of Boltzmann equation Maxwellian is a steady-state solution of BE , as will be demonstrated below. Consider a uniform gas without external forces. Thus, BE in steady-state reads ˆ ˆ 0 = d3u1 dΩ σ(Ω)|u − u1| (f ′f1′ − f f1). In order for f to be a solution of this equation, we require f f1 = f ′f1′ ⇔ log f (u) + log f (u1) = log f (u′) + log f (u′1). ⋆ Let χ(u) be a quantity conserved in collisions. Then χ(u) + χ(u1) = χ(u′) + χ(u′1). Thus, the most general form of f satisfying (⋆) is log f (u) = C0 + X r χr (u), where χr (u) are all the independent quantities conserved in collisions. Conservation of energy and all three components of momentum implies log f (u) = C0 + C1u2 + C2xux + C2y uy + C2z uz = −B(u − u0)2 + log A ∴ f (u) = A exp{−B(u − u0)2} "streaming" Maxwellian Hydrodynamics – Spring 2010 26 Boltzmann’s H-theorem Microscopic processes at molecular level presumably reversible, but macroscopic processes not. When a velocity distribution relaxes to a Maxwellian as a result of collisions, that is an irreversible process. (Reverse process extremely improbable.) When deriving BE, we assumed reversibility of microscopic processes. Yet, an equation leading to an irreversible evolution of the distribution equation resulted. How come? Let us construct a quantity H= ˆ d3u f log f. (H) H is a dimensionless quantity per unit volume (log is dimensionless.) H -theorem (Boltzmann). If the distribution function f appearing in (H) evolves according to BE, then H can never increase with time, i.e., dH ≤ 0. dt Proof. Differentiating (H) with respect to time ˆ dH ∂f = d3u (1 + log f ) dt ∂t Hydrodynamics – Spring 2010 27 For a uniform gas with F = 0, we have ˆ ˆ ˆ dH = d3u d3u1 dΩ σ(Ω) |u − u1| (f ′f1′ − f f1) (1 + log f ). dt Integration wrt. both u and u1 is performed, these dummy variables can be interchanged : ˆ ˆ ˆ dH = d3u d3u1 dΩ σ(Ω) |u − u1| (f ′f1′ − f f1) (1 + log f1). dt 1 1 2 (1)+ 2 (2) ˆ ˆ ˆ dH 1 = ⇒ d3u d3u1 dΩ σ(Ω) |u − u1| (f ′f1′ − f f1) [2 + log(f f1)]. dt 2 Reversibility of molecular interactions assumed, so for collisions{u′, u′1} → {u, u1} ˆ ˆ ˆ dH 1 = d3u′ d3u′1 dΩ σ(Ω) |u′ − u′1| (f f1 − f ′f1′ ) [2 + log(f ′f1′ )], dt 2 and since |u′ − u′1| = |u − u1| and d3u′ d3u′1 = d3u d3u1 we have ˆ ˆ ˆ 1 dH =− d3u d3u1 dΩ σ(Ω) |u − u1| (f ′f1′ − f f1) [2 + log(f ′f1′ )]. dt 2 (1) (2) (3) (4) 1 1 2 (3)+ 2 (4) ˆ ˆ ˆ dH 1 f f1 = ⇒ d3u d3u1 dΩ σ(Ω) |u − u1| (f ′f1′ − f f1) log ′ ′ . dt 4 f f1 For any positive α and β , we have (α − β) log(α/β) ≤ 0, which proves the theorem. Hydrodynamics – Spring 2010  28 Clearly, since a Maxwellian has C[f ] = 0, it also has dH/dt = 0. If we begin with a non-equilibrium distribution, H keeps decreasing until the gas relaxes to the equilibrium distribution (Maxwellian) when H attains a minimum value. We started by assuming reversibility, but ended up with the quantity H , which has a non-symmetric evolution in time. The arrow of time has been picked up. H -theorem is not the same thing as the increase of entropy S . Here, H is defined only for a uniform gas (entropy can be defined for more complicated systems); the usual definition of entropy is valid only for systems in TE, whereas H is defined for non-equilibrium systems. For dilute gases in TE we can, however, relate the two by S = −κBH + constant. Finally, we note a paradox associated with the H -theorem: suppose, at one instant of time, we reverse all the velocities of the particles in a gas. Then, the particles will follow their original paths backwards. If we had dH/dt < 0 before, now we should have dH/dt > 0! (Loschmidt’s paradox or reversibility paradox.) The resolution of the paradox lies in the assumption that the molecules before the collisions are always uncorrelated. (This is the so-called assumption of molecular chaos, and it is hidden in the statistical formulation of the interactions through the scattering cross-section.) We cannot, however, assume that the molecules would be uncorrelated after the collisions. Thus, in the reverse situation, the molecules are not uncorrelated before the collisions, and the assumptions behind the BE do not apply. So, in fact, the arrow of time in BE is picked up because of the assumption of uncorrelated velocities before the collisions. Hydrodynamics – Spring 2010 29 Conservation equation Consider a quantity χ(x, u), carried by the particles, conserved in binary collisions, i.e., χ + χ1 = χ′ + χ′1 (X) with the usual abbreviations. The total amount of χ per unit volume in the gas is ˆ ˆ f (x, u, t) = n(x, t)hχi(x, t). nχ(x, t) = d3u χ(x, u) f (x, u, t) = n(x, t) d3u χ(x, u) n(x, t) Aim: find an evolution equation for nχ, when f evolves according to BE. Multiply BE by χ and integrate over velocities: ˆ ˆ Df 3 χ d u = C[f ] · χ d3u. Dt Consider, first, the RHS. Perform the same tricks as for 1 + log f in place of χ when deriving the H -theorem to get ˆ ˆ ˆ ˆ 1 d3u d3u1 dΩ σ(Ω) |u − u1| (f ′f1′ − f f1)(χ + χ1 − χ′ − χ′1) C[f ] · χ d3u = 4 ´ so, because of (X), C[f ] · χ d3u = 0. Thus, for any quantity conserved in binary collisions, 0= Hydrodynamics – Spring 2010 ˆ Df 3 χd u = Dt ˆ   ∂f Fi ∂f  3 ∂f  + ui χd u + ∂t ∂xi m ∂ui 30 yielding ∂ 0= ∂t ˆ ∂ f χ d3u + ∂xi 1 − m ˆ 1 ∂Fi f χ d3u − ∂ui m ˆ f χui d3u − ˆ ˆ ∂χ 3 1 ui f d u+ ∂xi m ˆ ∂ (f χFi) d3u ∂ui ∂χ f Fi d3u. ∂ui The last term on the first line = 0, because it can be converted to a surface integral at |u| → ∞, ˆ ˆ ∂ (f χFi) d3u = f χFiui u dΩ = 0, ∂ui |u|→∞ where it has been assumed that f tends to zero more rapidly than χFiui u increases, as u → ∞. ´ Thus, recalling that d3u Qf = nhQi for any quantity Q, we get * ∂ ∂ ∂χ + n * ∂χ + n * ∂Fi + (nhχi) + Fi (nhuiχi) − n ui − − χ =0 ∂t ∂xi ∂xi m ∂ui m ∂ui This tells us how the density nχ = nhχi of any quantity χ, conserved in binary collisions, evolves in time. It will be used repeteadly when deriving the hydrodynamic equations. Hydrodynamics – Spring 2010 31 Derivation of HD equations Hydrodynamics – Spring 2010 32 The moment equations The equation * ∂ ∂ ∂χ + n * ∂χ + n * ∂Fi + Fi (nhuiχi) − n ui − − χ =0 (nhχi) + ∂t ∂xi ∂xi m ∂ui m ∂ui prescribes the recipe of how to move from microscopic (molecular quantity χ) theory to macroscopic (amount of χ per unit volume, nhχi) theory. Below, we will limit ourselves to forces independent of velocity. Conservation of mass. First, take χ = m. Thus, ∂ ∂ (nm) + (nmhuii) = 0, i.e., ∂t ∂xi ∂ρ + ∇ · (ρv) = 0, ∂t where ρ = nm is the mass density and v = hui is the average velocity of the particles. Conservation of momentum. Next, consider χ = muj . Thus, * ∂ ∂ ∂ ∂ ∂uj + = (ρvj ) + (nmhuiuj i) − n Fi (ρhuiuj i) − nFj 0 = (nmhuj i) + ∂t ∂xi ∂u ∂t ∂x i | {z i} =δij Hydrodynamics – Spring 2010 33 Noting that u − v is the particle velocity with respect to the average flow, we define ˆ Pij ≡ m d3u (ui − vi)(uj − vj )f = ρh(ui − vi)(uj − vj )i = ρ(huiuj i − huivj i − hviuj i + hvivj i) = ρ(huiuj i − vivj − vivj + vivj ) = ρ(huiuj i − vivj ). Thus, ρhuiuj i = Pij + ρvivj and ∂ ∂Pij ρ ∂ (ρvj ) + (ρvivj ) = − + Fj , i.e., ∂t ∂xi ∂xi m ρ ∂(ρv) + ∇ · (ρvv) = −∇ · P + F , ∂t m where P = Pij eiej is a second-order tensor. (This is clear, since it is defined through a dyadic product of vectors.) Using the conservation of mass, we can simplify the LHS       ∂ρ ∂v ∂v ∂(ρv) + ∇ · (ρvv) = ρ  + v · ∇v  + v  + ∇ · (ρv) = ρ  + v · ∇v  ∂t ∂t ∂t ∂t Note: on this course, we mark the dyadic product of vectors without any symbol, just placing the vectors side by side. Thus, ˆ P=m Hydrodynamics – Spring 2010 d3u (u − v)(u − v)f. 34 Conservation of energy. In case the gas is mono-atomic, translational kinetic energy is conserved in collisions, so we take χ = 21 m|u − v|2, i.e., the kinetic energy in the rest frame of the gas. Thus, ∂ ∂qi ∂ (ρǫ) + (ρǫvi) + + Pij Λij = 0, i.e., ∂t ∂xi ∂xi ∂(ρǫ) + ∇ · (ρǫv) + ∇ · q + P : Λ = 0 ∂t where ǫ = 21 h|u − v|2i internal energy per unit mass q = 12 ρh(u − v)|u − v|2i energy flux   1 ∂vi ∂vj    rate of strain tensor + Λ = Λij eiej with Λij =  2 ∂xj ∂xi Again, we can use   ∂ ∂ ∂ǫ (ρǫ) + (ρǫvi) = ρ  + v · ∇ǫ . ∂t ∂xi ∂t Note: we denote the double inner product of two tensors with a colon: A : B = Aij Bji . In the case of P and Λ, both tensors are symmetric, so Pij = Pji and Λij = Λji. Thus, P : Λ = Pij Λij . Hydrodynamics – Spring 2010 35 Summary of moment equations We have derived the moment equations (termed so because they are obtained as velocity moments of the Boltzmann equation) ∂ρ ∂(ρvi) =0 + ∂t ∂xi ∂vj 1 ∂Pij Fj ∂vj + vi =− + ; ∂t ∂xi ρ ∂xi m   ∂ǫ ∂qi ∂ǫ  ρ  + vi =− − Pij Λij ; ∂t ∂xi ∂xi j = 1, 2, 3   1 ∂vi ∂vj    Λij =  + 2 ∂xj ∂xi so there are altogether five equations. The unknowns are velocity moments of f : ˆ ˆ ˆ 1 d3u uj f ; Pij = m d3u (ui − vi)(uj − vj )f ; ρ = m d3u f ; vj = n ˆ ˆ 1 1 ǫ= d3u (ui − vi)(ui − vi)f ; qj = d3u (uj − vj )(ui − vi)(ui − vi)f 2n 2 so there are altogether 1+3+6+1+3=14 unknowns. Thus, the moment equations do not yet constitute a dynamical theory of fluids. We could, of course, introduce more moment equations by taking higher-order moments of BE. However, they would always introduce higher order moments of the distribution function because of the term u·∇f in BE. Thus, a way of truncating the chain has to be developed in order to obtain a dynamical theory. Hydrodynamics – Spring 2010 36 Excursion to stellar dynamics Boltzmann equation describes the evolution of point-like masses interacting via binary collisions. Could BE be applied to the evolution of a large group of stars, like a cluster or a galaxy? In short, no. Stars interact via gravity, which is a long-range force. However, in many systems, collisions between stars are extremely rare. For galaxies, for example, the collisional relaxation time can be estimated to be longer than the age of the universe. If collisions are unimportant, the collisionless Boltzmann equation can be applied. In this case, contributions from the collision integral vanish identically, and moment equations derived for the BE apply (Jeans 1922). In this case, of course, the force F is not external, but should include the contribution from self-gravity of the system. The moment equations can be applied to astronomical data. As an example, consider the matter in the neighborhood of the Sun. Assume: 1. that the system is in a steady state and 2. that the variations of quantities is much faster in the direction perpendicular to the galactic plane than within the directions in the plane. Thus, the uz -moment equation gives 0= ∂ 1 d d ∂ (ρ⋆vz ) + (n⋆hu2z i) (ρ⋆huiuz i) − n⋆Fz = (ρ⋆hu2z i) − gz ρ⋆ ⇒ gz = ∂t ∂xi dz n⋆ dz Hydrodynamics – Spring 2010 37 Oort limit Oort (1932) estimated the gravitational field at different heights from the galactic plane using statistics on stars and equation 1 d gz (z) = (n⋆hu2z i) n⋆ dz along with the assumption hu2z i = constant. Now, g fulfills the Poisson equation, −4πGρ = ∇ · g = dgz , dz so the amount of matter producing the gravitational field in the solar neighborhood can be estimated. It turns out to be ρ ≈ (7 − 10) × 10−21 kg m−3. Calculating the the density by estimating the amount of material in the visible stars gives ρ⋆ ≈ (4 − 6) × 10−21 kg m−3. The difference, ρ − ρ⋆ ≈ (1 − 6) × 10−21 kg m−3, gives, therefore, an upper limit to the density of the interstellar gas, known as the Oort limit. Hydrodynamics – Spring 2010 38 Towards the dynamical theory of fluids Recall the moment equations: ∂ρ ∂(ρvi) + =0 ∂t ∂xi ∂vj ∂vj 1 ∂Pij Fj + vi =− + ; ∂t ∂xi ρ ∂xi m  j = 1, 2, 3  ∂ǫ ∂qi ∂ǫ  ρ  + vi =− − Pij Λij , ∂t ∂xi ∂xi   1 ∂vi ∂vj   and  + where Λij =  2 ∂xj ∂xi ˆ ˆ ˆ 1 ρ = m d3u f ; vj = d3u uj f ; Pij = m d3u (ui − vi)(uj − vj )f ; n ˆ ˆ 1 1 d3u (ui − vi)(ui − vi)f ; qj = d3u (uj − vj )(ui − vi)(ui − vi)f ǫ= 2n 2 Our task is to derive enough relations between the fourteen unknowns, ρ, vj , Pij , ǫ and qj so that the system of equations becomes closed and a dynamical theory for fluids is established. Hydrodynamics – Spring 2010 39 Zero-order approximation Collisions are the only way of transmitting momentum and energy in a fluid consisting of neutral particles and, thus, essential for fluid-like behaviour. We noted that collisions drive the distribution towards a Maxwellian equilibrium distribution, with a possible non-zero mean velocity. Thus, start by assuming that the distribution function is of form 3/2  m    f (x, u, t) = n(x, t)  2πκBT (x, t)   2  m{u − v(x, t)}  , exp − 2κBT (x, t) i.e., a Maxwellian with local values of the mean velocity, density and temperature. Start by calculating Pij (x, t) Pij = m ˆ d3u (ui − vi)(uj − vj )f U =u−v =  mn  3/2 m  2πκBT ˆ  − d3U UiUj exp  2  mU   = p δij , 2κBT where p = nκBT is the pressure of the fluid. Likewise, ˆ ˆ 1 1 3 κ T B ǫ= ; qj = d3U U 2 f = d3U Uj U 2f = 0, 2n 2 m 2 where ǫ is the internal energy (per unit mass) for a monoatomic gas. Finally, Pij Λij = pΛii = p Hydrodynamics – Spring 2010 ∂vi = p∇ · v. ∂xi 40 Thus, the zero-order moment equations can be written in the simplified form as: ∂ρ + ∇ · (ρv) = 0 ∂t ∂v 1 F + v · ∇v = − ∇p + ∂t ρ m   ∂ǫ ρ  + v · ∇ǫ = −p∇ · v, ∂t i.e., as five equations for the six quantities ρ, v , p and ǫ. However, the three thermodynamic quantities can be expressed as a function of number density and temperature, i.e., ρ = mn, p = nκBT and ǫ = 23 κBT /m. Thus, we conclude that the number of independent variables equals the number of equations. We have ended up with a dynamical theory of a fluid. This theory, however, lacks some important features of a theory for real fluids 1. Since q = 0, there is no transport of internal energy by other means than convection. I.e., heat conduction is missing from the theory. 2. Since P is diagonal, there is no viscosity in our fluid (proof will follow later). Thus, momentum transport inside the fluid is incomplete and relative motions of fluid layers are not damped out. In short, our description of the fluid lacks a proper description of transport phenomena. What went wrong in the derivation? The starting point, a local Maxwellian distribution, is too restrictive: if there is a temperature gradient in the fluid, particles arriving to a position from the direction of the gradient have slightly higher energies than those arriving from the opposite direction. Transport phenomena are clearly associated with (small) departures from Maxwellian distribution. Hydrodynamics – Spring 2010 41 Approximating C[f ] Consider where f (0) f (x, u, t) = f (0)(x, u, t) + g(x, u, t), is a Maxwellian distribution constructed with the moments ˆ ˆ ˆ m 1 d3U U 2 f d3u uf ; T = n = d3u f ; v = n 3nκB (U = u − v) and g is a small correction. Start from the collision integral ˆ ˆ C[f ] = d3u1 dΩ |u − u1|σ(Ω)(f ′f1′ − f f1). Assume that those distributions over which we will integrate, i.e., f ′, f ′1 and f1, are all well represented ′ (0) ′ ′ (0) (0) (0) ′ by Maxwellians f (0) , f1 and f1 , and use the property of a Maxwellian that f (0) f1 = f (0)f1 . Thus, ˆ d3u1 ˆ (0) (0) ˆ ˆ (0) dΩ |u − u1|σ(Ω)(f (0)f1 − f f1 ) = (f (0) − f ) d3u1 dΩ |u − u1|σ(Ω)f1 ˆ ˆ = −g(x, u, t) dΩ σ(Ω) d3u1|urel| f (0)(x, u1, t) = −σtot nh|urel|i g(x, u, t), C[f ] ≈ where h|urel|i(|u|, T ) is the average velocity between the colliding particles and h|urel|inσtot, thus, represents the average scattering rate of particles of speed |u|. Hydrodynamics – Spring 2010 42 BGK Equation We got C[f ] ≈ −σtot nhureli (f − f (0)) This motivates replacing collision integral by −(f − f (0))/τ , giving (Bhatnagar, Gross & Krook 1954) (0) F ∂ f − f ∂  f = − +u·∇+ · . ∂t m ∂u τ   The BGK equation describes a distribution relaxing towards a Maxwellian with a relaxation time τ taken to be constant. The relaxation time can, for example, be replaced by the mean scattering time for hard-sphere molecules (of diameter d),  1/2 πκ T 1 B  = 4nd2  , τ m but actually it is a parameter of the theory. We shall see that transport phenomena can be qualitatively described the BGK equation, but the values of the transport coefficients are not exact. The reason lies in the fact that Boltzmann collision integral actually demands 1/τ ∝ h|urel|i, which is a function of u. However, even if values of the coefficients are not accurate, the model provides us with the framework under which the transport processes can be identified. Hydrodynamics – Spring 2010 43 Departure from Maxwellian First, estimate the order of magnitude of the departure from Maxwellian. If the system is in an approximate steady state and there are no external forces in the system, then |u|τ (0) λ (0) |u|f (0) |g| ≈ ⇒ |g| ≈ f ≈ f , L τ L L where L is the gradient scale height of the system. Thus, departures from a Maxwellian will be small, if mean free path ≪ gradient scale height of the system. By using α = λ/L as a small expansion parameter, we can write f = f (0) + αf (1) + α2f (2) + ... where f (i)’s are of the same order. Substituting this expansion to BGK equation allows the estimation of the successive terms. Here, we consider only the first-order expansion: g = αf (1)   ∂ F ∂  (0) = −τ  + u · ∇ + · f , ∂t m ∂u where, e.g., ∂f (0) ∂n ∂f (0) ∂T ∂f (0) ∂vi ∂f (0) . = + + ∂t ∂t ∂n ∂t ∂T ∂t ∂vi Here, many terms can be eliminated by using the zero-order moment equations (recalling that p = nκBT ). This gives (exercise)      1 ∂T  m 1 5 m g = −τ  U2 −  + Λij UiUj − δij U 2 f (0) Ui T ∂xi 2κBT 2 κBT 3 Hydrodynamics – Spring 2010 44 Heat flux (0) 3 Consider now the moments calculated using f = f + g . Recall, that by construction n = d u f, ´ 3 ´ 3 nv = d u uf and 3nκTT = m d U U 2f , but other moments, in particular q and the off-diagonal elements of Pij are of interest. Let us start from the heat flux m qi = 2 ˆ m = −τ 2 ´ d3U UiU 2g      m 5 1 1 ∂T  m Uj d3U UiU 2  U2 −  + Λjk Uj Uk − δjk U 2 f (0) T ∂xj 2κBT 2 κBT 3 ˆ   5 m m 1 ∂T U 2 −  f (0) (no summation over i) d3U Ui2 U 2  = −τ 2 T ∂xi 2κBT 2 ˆ where odd terms of the integrand have been omitted. Since the integral has the same value for all i = 1, 2, 3, we can replace it with third of the sum over i. Thus, q = −K∇T ˆ   5 κ2BT m 5 τm (0) 3 4 2  d UU U − f = τn , K= 6T 2κBT 2 2 m where K is the coefficient of thermal conductivity. Thus, Fourier’s law of heat conduction is recovered. Hydrodynamics – Spring 2010 45 P tensor Consider next the P tensor, defined by Pij = m πij = m ˆ ´ d3u UiUj f = nκBT δij + πij , where d3U UiUj g      1 1 ∂T 5 m m U2 −  + Λkl Uk Ul − δkl U 2 f (0) = −τ m d3U UiUj  Uk  T ∂xk 2κBT 2 κBT 3 ˆ   τ m2 1 Λkl d3U UiUj Uk Ul − δkl U 2 f (0), =− κBT 3 ˆ and odd terms have again been eliminated. Clearly, πii = 0 (integrand odd for k 6= l and terms equal for k = l), so πij is a traceless tensor whose elements depend linearly on Λkl . Thus, πij = −2µ(Λij − 13 δij tr Λ) = −2µ(Λij − 13 δij ∇ · v), where the proportionality constant −2µ can be obtained, e.g., from τ m2 Λkl −2µΛ12 = π12 = − κBT   1 d3U U1U2 Uk Ul − δkl U 2 f (0) 3 ˆ τ m2 =− (Λ12 + Λ21) d3U U12U22f (0) = −2Λ12τ nκBT κBT ˆ giving µ = τ nκBT . We shall see later, that µ is the coefficient of (shear) viscosity. Hydrodynamics – Spring 2010 46 Comparison of µ with experiment Let us first consider the value of viscosity, µ = τ nκBT , obtained by using the mean hard-sphere scattering rate  1/2 1 πκ T B  = 4nd2  τ m as the inverse relaxation time in BGK model. Thus, 1  mκBT 1/2 . µ= 2 4d π   (Note the error d ← a in Chourdhuri.) A rigorous derivation utilizing the full collision integral yields √   5 πmκBT 5  mκBT 1/2 µ= = , 16 σtot 16d2 π so the BGK model underestimates the rigorous result by one part in 16, only! Notes: • Theory: µ independent of n. Experiment: confirmed (for gases). Why? More molecules → more agents of momentum transport but shorter mean free path. √ • Theory: µ ∝ T . Experiment: µ increases faster in gases. Why? Higher temperature → effective radius of molecule smaller → higher µ. Liquids: µ usually decreases with temperature. Hydrodynamics – Spring 2010 47 Comparison of K with experiment Heat conduction in the BGK model is 5 κB K= µ . 2 m 3 The internal energy per unit mass is ǫ = 2 κBT /m giving the specific heat per unit mass as cV = Thus, 3 κB . 2m 5 K = . µcV 3 A more rigorous model (full Boltzmann treatment) gives this number as 5/2 for hard-sphere molecules. Experimental values for monoatomic inert gases are Hydrodynamics – Spring 2010 Gas He Ne Ar Kr Xe K µcV 2.45 2.52 2.48 2.54 2.58 48 First-order moment equations After now expressing Pij and qi as functions of n, T and vi, we may finally write down the moment equations incorporating the transport phenomena. Thus, treating µ as a constant,   1 ∂p ∂  ∂Pij Λij − δij ∇ · v  = − 2µ ∂xi ∂xj ∂xi 3 =   1 ∂ ∂  ∂vj ∂vi  2µ ∂ ∂vk ∂p ∂p 2  + ∇ vj +  −µ + = − µ ∇ · v ∂xj ∂xi ∂xi ∂xj 3 ∂xj ∂xk ∂xj 3 ∂xj so     ∂v ρF 2 1   ρ + v · ∇v = −∇p + µ ∇ v + 3 ∇(∇ · v) + ∂t m Likewise, giving   (P)   1 Pij Λij = pδij Λij + πij Λij = p∇ · v − 2µ Λ : Λ − (∇ · v)2 3     ∂ǫ 1 ρ  + v · ∇ǫ = −p∇ · v + ∇ · (K∇T ) + 2µ Λ : Λ − (∇ · v)2 . ∂t 3 Along with the conservation of mass ∂ρ + ∇ · (ρv) = 0 ∂t equations (P) and (E) form a closed system, i.e., a dynamical theory of fluids. Hydrodynamics – Spring 2010 (E) (M) 49 Hydrodynamic equations Typically, we make some simplifications to the first order moment equations, when treating fluids • • • • Neglect the spatial variation of µ, as done already when deriving the equation of motion (P) Neglect the effects of compressibility (∇ · v ) in viscous forces in (P) Neglect the effect of viscous production of heat in the energy equation (E) Simplify the notation by writing F ← F /m. Thus, we get the following hydrodynamic equations: ∂ρ + ∇ · (ρv) = 0 ∂t ∂v 1 µ + v · ∇v = − ∇p + F + ∇2v ∂t ρ ρ  (Navier–Stokes equation)  ∂ǫ ρ  + v · ∇ǫ = −p∇ · v + ∇ · (K∇T ) ∂t We finally recall the starting point of our derivation, that the distribution remains close to a local Maxwellian. This was shown to be the case when the mean free path was much shorter than the gradient scales of the system. Thus, we expect the hydrodynamic equations to break down if this condition is not met. Hydrodynamics – Spring 2010 50 Ideal fluids Hydrodynamics – Spring 2010 51 Macroscopic derivation of the HD equations Above, we derived the hydrodynamic equations (for dilute gases) starting from molecular dynamics. Next, we will consider the macroscopic derivation of hydrodynamic equations, treating fluids as continua, which establishes that HD equations are valid for dense fluids as well. We shall, again, be aiming at finding a dynamical theory for fluids. Thus, taking time derivatives is in a key role. When treating continua, one introduces two different kinds of time derivatives: Eulerian: This is the (partial) time derivative ∂/∂t taken at a fixed point in space, i.e., regarding x constant. Lagrangian: This is the (total) time derivative d/dt taken at a position following a fluid element moving at velocity v . Consider a quantity Q(x, t) related to the fluid (e.g., temperature, density, component of velocity). In an infinitesimal time δt, fluid element that is at x at time t will move to x + v δt. Thus, Q(x + v δt, t + δt) − Q(x, t) dQ ≡ lim δt→0 dt δt ∂Q (x, t) + v δt · ∇Q(x, t) Taylor: Q(x + v δt, t + δt) ≈ Q(x, t) + δt ∂t ∴ Hydrodynamics – Spring 2010 dQ ∂Q = + v · ∇Q. dt ∂t 52 Equation of continuity Consider the amount of mass enclosed in a fixed volume V : ˆ MV (t) = d3x ρ(x, t). V Thus, ˆ ˆ d dMV ∂ρ = d3x ρ(x, t) V fixed = d3x . dt dt V ∂t V On the other hand, if mass is neither created nor destroyed, the temporal change of MV can only be due to the flow of mass through the surface ∂V bounding the volume. Thus, ˛ ˆ dMV =− ρv · ds = − d3x ∇ · (ρv) dt ∂V V so ˆ V and since V is arbitrary   ∂ρ d3x  + ∇ · (ρv) = 0 ∂t ∂ρ + ∇ · (ρv) = 0. ∂t (M) Thus, dρ ∂ρ + v · ∇ρ + ρ∇ · v ⇒ = −ρ∇ · v. ∂t dt Clearly, the fluid is incompressible (i.e., dρ/dt = 0) iff ∇ · v = 0. 0= Hydrodynamics – Spring 2010 53 The equation of motion Consider a fluid element of volume δV and, thus, mass ρ δV . The acceleration of the fluid element is dv/dt so Newton’s second law applied to the fluid element states that dv = δF = δF body + δF surf . ρ δV dt v −Pij ds j e i ds n The force is broken in two components: Body force δF body = ρ δV F contains forces that affect all the δ Fbody molecules of the fluid element; usually an external force in a neutral fluid. Surface force dF surf = −P · ds contains forces that are excerted on the fluid element through its surface element ds. As both dF surf and ds are vectors, the proportionality relation requires P to be a second-rank tensor. The total surface force acting on the element is ˛ ˆ ∂Pij ∂Pij 3 d x=− δV. (δFsurf )i = − Pij dsj = − ∂x ∂x j j ∂δV δV Thus, ρ dvi ∂Pij =− + ρFi dt ∂xj (P) which is identical to the result from the microscopic derivation (but without explicit relation of Pij to the fluid properties). The symmetry of Pij follows from the conservation of angular momentum. Hydrodynamics – Spring 2010 54 Surface forces for static fluids. Euler equation If the fluid is in static equilibrium, then it is experimentally established that dF surf k n ⇒ dF surf = −p ds, where p is a scalar pressure. Thus, Pij = p δij . (SE) Although this simple law holds for fluids at rest, it is no longer valid for fluids in motion (unless v is the same constant everywhere). This becomes clear if one considers a shear flow, v = y/τ ex, and takes the fluid element to have a surface in the xz plane. If (SE) holds, then the surface force is in ey direction and it is unable to damp out the relative motion between the fluid above and below the xz plane, contrary to experimental reality. We will first, however, restrict our attention to fluids where (SE) holds regardless of the motional state of the fluid. Thus, the equation of motion becomes ρ i.e., ∂p dvi =− + ρFi, dt ∂xi 1 ∂v + v · ∇v = − ∇p + F . ∂t ρ This equation is called the Euler equation and the fluids obeying it ideal fluids. Hydrodynamics – Spring 2010 55 First law of thermodynamics The first law of thermodynamics states dQ = dU + p dV, where dQ is the amount of heat added to the system dU is the change of the internal energy of the system p dV is the work done by the system (V is the system volume) Thermodynamic variables can be classified to extensive and intensive. If the variable V obeys V1+2 = V1 + V2 when combining two distinct thermodynamic systems into one, it is called extensive. If not, it’s called intensive. Clearly, Q, U and V are extensive, but p is intensive. Although the fluid is not in general in TE, the first law of TD can be applied to a small fluid element of mass δm. Thus, define the extensive variables as quantities per unit mass multiplied by δm: dQ = δm dq dU = δm dǫ   dρ 1 dV = δm d   = −δm 2 ρ ρ p ∴ dq = dǫ − 2 dρ ρ Hydrodynamics – Spring 2010 56 Energy equation Divide dq = dǫ − (p/ρ2) dρ by dt, and use dρ/dt = −ρ∇ · v to get p dρ dq dǫ = 2 + , dt ρ dt dt dǫ = −p∇ · v − L, dt where L = −ρ dq/dt is the heat loss rate per unit volume. ρ Heat flows from hotter to colder parts of the system at a rate proportional to the temperature difference. Thus, the heat flux can be given as q = −K∇T, where K is the coefficient of thermal conductivity. Rate of heat loss by conduction of the fluid element through its surface is, therefore, ˛ ˆ q · ds = (∇ · q) d3x. ∂ δV δV The corresponding heat loss rate per unit volume is L = ∇ · q = −∇ · (K∇T ) giving Hydrodynamics – Spring 2010   ∂ǫ ρ  + v · ∇ǫ = −p∇ · v + ∇ · (K∇T ) ∂t (E) 57 Macroscopic derivation of HD – a summary The continuum hypothesis states that matter is distributed continuously inside a body. The body, which in HD consists of a fluid, can be broken into small fictious (fluid) elements, each having a velocity v and carrying physical quantities like density ρ and temperature T . Continuum hypothesis yields immediately ∂ρ + ∇ · (ρv) = 0. ∂t Newton II applied on fluid elements dv = −∇ · P + ρF , dt where P = Pij eiej is a second-rank tensor representing the surface forces (per unit area) excerted on the fluid element by the fluid surrounding it. In a (local) thermal equilibrium, Pij = p δij is known empirically to apply. ρ By assuming LTE and considering the first law of thermodynamics, the energy equation for the fluid could be derived dǫ ρ = −p∇ · v − ∇ · q, dt where ǫ is the internal energy per unit mass and q gives the heat flux. We need to be able to express all the fourteen variables ρ, vi, Pij , ǫ and qi as functions of five variables, in order to obtain a dynamical theory for fluids. In the macroscopic theory, Pij , ǫ and qi have to be linked to vi and two thermodynamic variables via constitutive equations, which are essentially empirical relations. The relations Pij = p δij and qi = −K∂T /∂xi are examples of such relations, but depending on a material and the physical setup, more complicated relations may have to be applied. Hydrodynamics – Spring 2010 58 Vorticity Consider the Euler equation for a system, where the body force is conservative, i.e., F = −∇ϕ: 1 ∂v + v · ∇v = − ∇p − ∇ϕ. ∂t ρ Since for an arbitrary vector field A, A · ∇A = 21 ∇(A · A) − A × (∇ × A), we have 1 ∂v = v × (∇ × v) − ∇p − ∇(ϕ + 12 v 2). ∂t ρ Taking now a curl of both sides we get ∇ρ × ∇p ∂ω = ∇ × (v × ω) + , ∂t ρ2 where ω ≡ ∇ × v is the vorticity of the fluid. Hydrodynamics – Spring 2010 59 Vorticity equation in incompressible fluids For fluids with constant density (incompressible fluids), the equation for vorticity simplifies to ∂ω = ∇ × (v × ω). ∂t (V) The assumption incompressibility is clearly a good one for liquids like water. However, it applies very well also for gases, provided that they move much slower than compressions spread around. The spreading of compressions occurs at sound speed, so for fluid velocities much less than the speed of sound, incomressibility is a very good approximation. Incompressible fluids obey ∇·v =0 Any vector field can be solved, if its divergence and curl are specified. Thus, the vorticity ω =∇×v (I) (W) prescribes the velocity field in an incompressible fluid and equations (V), (I) and (W) constitute a dynamical theory of incompressible ideal fluids. Hydrodynamics – Spring 2010 60 Redundancy of energy equation. Barotropic fluids For incompressible ideal (real) fluids, Euler (Navier-Stokes) equation along with ∇ · v = 0 and ρ = constant provides a closed set of equations. (Four variables, i.e., v and p, and four equations.) Thus, for incompressible fluids, energy equation is redundant. In many astrophysical situations, however, incompressibility is too restrictive an assumption: many such systems • have large density variations because of the large system size and strong gravitational fields, and/or • exhibit flow speeds comparable to or exceeding the speed of sound. Thus, the energy equation has to be solved together with the momentum and continuity equations. The study of such fluids is referred to as gas dynamics. In certain problems, however, it is possible to assume a functional relation between p and ρ: p = p(ρ). This a called a barotropic relation and fluids obeying it barotropic fluids. Energy equation is redundant for barotropic fluids—the full dynamical theory now consists of the continuity equation and Euler (Navier–Stokes) equation. As ∇ρ k ∇p, the vorticity equation for ideal barotropic fluids becomes ∂ω = ∇ × (v × ω). ∂t Hydrodynamics – Spring 2010 61 Hydrodynamic equations in conservative form Equation of continuity has the form ∂(mass density) + ∇ · (mass flux) = 0. ∂t It is in a standard conservative form, expressing the fact that mass in a fixed region of space can only change because of its flux through the boundary of the region. As established before (using the equation of continuity) ∂(ρv) dv + ∇ · (ρvv) = ρ = −∇ · P + ρF . ∂t dt In the absense of (external) body forces, conservation of momentum can be stated as ∂(ρv) + ∇ · (ρvv + P) = 0, ∂t where ρv is the momentum density of the fluid and T = ρvv + P can, thus, be regarded as the momentum flux. Using this momentum equation in conservative form, it is easy to show that in the absense of external forces, the total amount of momentum in a fluid is conserved (excerise). Hydrodynamics – Spring 2010 62 The energy equation can also be cast into a conservative form. To do that, we first note that energy density is a sum of kinetic and internal energy densities, i.e., ρ( 12 v 2 + ǫ). In the absense of external forces (and neglecting work done by viscous forces) we get (exercise)      ∂  2 1 1 2 ρǫ + 2 ρv + ∇ · ρv 2 v + w − K∇T = 0, ∂t establishing that energy flux is  1 2 2v  ρv + w − K∇T, where w = ǫ + p/ρ is enthalpy per unit mass. Clearly, by integrating the conservative energy equation over the whole system, one can establish that the total energy of the fluid is conserved (under the assumptions used to derive the energy equation). The conservative form of HD equations is especially useful in numerical work. Hydrodynamics – Spring 2010 63 Hydrostatics In hydrostatics, we are interested in finding the static equilibrium solutions of the HD equations. Thus, we set v = 0 and ∂/∂t = 0 to get 0 = −∇p + ρF 0 = ∇ · (K∇T ) from the momentum and energy equations, respectively. Note that we have three unknowns (ρ, p, T ) and four equations here, so the system is actually overdetermined. This means that the momentum equation does not have a solution for an arbitrary force field F (x). For an incompressible (ρ constant) or a barotropic fluid with ρ = ρ(p), we can use ρ−1∇p = ∇P , where ˆ dp . P(p) ≡ ρ Thus, the force balance states ∇P = F . Clearly, the solution exists only if the force is conservative, i.e., ∇P = F = −∇ϕ ⇒ P = P0 − ϕ where P0 is a constant. For incompressible or barotropic fluids the energy equation is redundant. Hydrodynamics – Spring 2010 64 Examples Incompressible fluid : P0 − ϕ(x) = P = ˆ dp p = ⇒ p(x) = p0 − ρ ϕ(x) ρ ρ Consider a homogeneous gravitational field F = −gez , i.e., ϕ = gz . Thus, p(z) = p0 − gρz. As we know from experience in everyday life, hydrostatic pressure increases linearly with depth (−z ) in an incompressible fluid. Isothermal ideal gas: consider p = RρT , where T is constant and R = κB/m. Thus, P0 − ϕ(x) = RT ˆ   dp ϕ(x)  −  = RT ln p ⇒ p(x) = p0 exp  p RT so for ϕ = gz z! RT κBT p(z) = p0 exp − , where H = = H g mg Here, H is the pressure scale height of an isothermal atmosphere. Note that the density obeys the same law of stratification, ρ = ρ0e−z/H . Hydrodynamics – Spring 2010 65 Principle of Archimedes Consider an incompressible fluid in static equilibrium in a homogeneous gravitational field g = −gez . Thus, p(z) = p0 − gρ z. Assume that a solid body of volume V and mass M is immersed in the fluid, and that the body is at rest with respect to the fluid. What is the total force acting on the body? The force acting on the surface of the body resulting from the pressure of the fluid is1 ˛ ˛ ˆ ˆ ˆ ∂p 3 dx F surf = − p ds = − (pI) · ds = − ∇ · (pI) d3x = − ∇p d3x = −ez ∂z ∂V ∂V V V V ˆ gρ d3x = gρV ez = −Mfluid g = ez V This states the principle of Archimedes: A body immersed in a fluid experiences a buoyant force equal to the weight of the fluid it displaces. In addition, the body is affected by the gravitational force M g so the resultant force acting on the body is R = (M − Mfluid) g. Note that once the body body starts to move as a result of this force, the force becomes modified. 1 When transforming the surface integral into volume integral, the pressure field of the fluid is continued to the region inside the solid body. Hydrodynamics – Spring 2010 66 Polytropic atmosphere The assumption of constant T does not generally hold throughout the entire Earth’s atmosphere. Thus, replace this assumption by a polytropic equation of state, p = p0(ρ/ρ0)γ , where γ is the polytropic index of the gas. Isothermal state corresponds to γ = 1. Adiabatic state corresponds to γ = cp/cV , i.e., to γ = 53 for a monoatomic gas and γ = gas. However, other values, such as γ = 1.2, are often used in practical calculations. We have ρ = ρ0(p/p0)1/γ so P= Using ϕ = gz gives, thus, ˆ 7 5 for a diatomic dp p0γ  p 1−1/γ = ρ ρ0(γ − 1) p0   p0γ  p 1−1/γ + gz = const. ρ0(γ − 1) p0 At z = 0, take p = p0 (sea-level pressure). Thus, const. = γp0/[ρ0(γ − 1)] and   γ 1 " " z # γ−1 z # γ−1 p z# p = p0 1 − , ρ = ρ0 1 − and T = = T0 1 − h h Rρ h " where now h = Hydrodynamics – Spring 2010 γ p0 γ kB T0 = γ − 1 gρ0 γ − 1 mg gives the atmospheric thickness (for γ > 1). 67 Courtesy of Windows to the Universe, http://www.windows.ucar.edu Hydrodynamics – Spring 2010 68 Solar corona In case of the atmosphere of the Earth, the assumption of homogeneous gravitational field is not a bad one. However, for extended atmospheres of stars like the Sun, we have to discard this assumption. Eclipse photos ⇒ Sun has a thin outer atmosphere called the corona. Corona extends to large distances from the Sun ⇒ global hydrostatic balance considered. Corona has very little mass ⇒ F =− GM⊙ er . r2 Despite the relatively non-spherical appearance, the assumption of spherical symmetry used as a starting point. Thus, GM⊙ GM⊙ p dp = −ρ 2 = − 2 dr r r RT   1 d  2 dT  Kr 0 = ∇ · (K∇T ) = 2 r dr dr Hydrodynamics – Spring 2010 69 Consider, first, the energy equation. Observations tell us that corona is very hot (1–2 MK) so it is in a plasma state. Thus, it can be shown that K ∝ T 5/2. Therefore,   d  2 5/2 dT  r T 0= dr dr r2 T 5/2 ⇒ dT = −C1 dr ⇒ 2 7 T 7/2 = C0 + C1 r so, assuming that T → 0 as r → ∞, we get r0 !2/7 . r T (r) = T0 This derivation assumes that nothing is heating the gas above the coronal base at r = r0, so the temperature is obtained simply by applying heat conduction from the coronal base kept at T = T0. Now, the momentum equation becomes GM⊙  r 2/7 p dp =− 2 dr r r0 RT0   yielding, using p = p0 at r = r0,      7GM⊙  r0 !5/7 − 1 . p = p0 exp  5RT0r0 r  − 75 (GM⊙/RT0r0)  This has an asymptotic value p∞ = p0 exp , which is much larger than the interstellar pressure. Thus, Parker (1958) concluded that corona cannot be in a static eqilibrium but has to expand as a solar wind . Hydrodynamics – Spring 2010 70 Steady flows. Streamlines After considering hydrostatics we turn to the next simplest HD problem, the study of steady flows.2 For that purpose, introduce the concept of a streamline. Streamline. A streamline is defined as a line everywhere tangential to the velocity field v(x, t). Thus, dx = v(x, t) dτ ⇒ dτ = dx dy dz = = , vx vy vz where τ is a parameter, gives the equations of the streamline, when integrated regarding t as a constant. Note that the path of a fluid element, obtained by integrating dx = v(x, t) dt, is given by a streamline in a steady flow, i.e., if ∂v/∂t = 0. 2 By steady flows we mean flows that are independent of time. Hydrodynamics – Spring 2010 71 Example. Calculate the streamlines of the velocity field v = ω(xey − yex) + V ez . −     dy dx = ⇒ ωy ωx dy dφ dz = = ⇒ ωx ω V x = R cos φ x dx + y dy = 0 ⇒ x2 + y 2 = R2 ⇒   y = R sin φ     0 x = R cos ω z−z V φ z − z0 = ω V ⇒  0  y = R sin ω z−z V The streamlines are, thus, helixes around the z -axis. This is, of course, motion consisting of rigid rotation √ around and translation along the z -axis (Λ = 0). The speed of each fluid element is constant |v| = ω 2R2 + V 2. Example. What about v = az(xey − yex + Lez ) ? − dx dy = ⇒ azy azx     x = R cos φ x dx + y dy = 0 ⇒   y = R sin φ dy dφ dz = = ⇒ azx az azL     0 x = R cos z−z z − z0 L φ= ; ⇒  y = R sin z−z0 L L Streamlines still the same helixes, but now the fluid element has an accelerating speed |v| = √ |az| R2 + L2. This is not rigid-body motion, as although Λxx = Λyy = Λxy = 0, Λzz ∂vz = = aL; ∂z Hydrodynamics – Spring 2010 Λxz   1 ∂vz ∂vx  = − 21 ay; =  + 2 ∂x ∂z Λyz   1 ∂vz ∂vy  1 = 2 ax. =  + 2 ∂y ∂z 72 Bernoulli’s principle for steady flows Consider conservative external forces, F = −∇Φ. Thus, 1 − ∇p − ∇Φ = v · ∇v = ∇( 12 v 2) − v × (∇ × v). ρ Take an integral along a streamline ˆ     1   = 0 + dl · ∇( 21 v 2) − |v × (∇ × v) ∇p + ∇Φ  {z } ρ ⊥dl ˆ dp 1 2 v + P + Φ = constant P = ⇒ 2 ρ which states Bernoulli’s principle for steady flows. For example, considering an incompressible fluid in homogeneous gravitational field, we get 1 2 2v + p + gh = constant, ρ where h is the height from some horizontal level. Let us consider some examples of the application of Bernoulli’s principle. Hydrodynamics – Spring 2010 73 Flow in a horizontal pipe with variable cross section Consider incompressible flow in a horizontal pipe, whence h = constant. Thus, 1 2 2v + p/ρ = constant. Furthermore, conservation of mass implies ∇ · v = 0. Integrate the equation over an arbitrary part of the pipe, allowing the ends of the volume to have a different cross-section. Thus, ˛ v · ds = 0. S The only non-zero contributions to the integral come from the ends (as v · n = 0 elsewhere), whence ˆ ˆ (v · n) ds. (v · n) ds + 0= S1 S2 The flow is, to a good approximation, perpendicular to the cross-sectional area of the tube, so this yields v1A1 = v2A2, where Ai is the area of the surface S1. Thus, p2 − p1 = 12 ρ(v12 − v22) = 12 ρ(A22/A21 − 1)v22 A1 A2 and the pressure is, therefore, the higher the larger the local cross-sectional area of the tube. Hydrodynamics – Spring 2010 74 Flow through a small valve at the bottom of a tank Consider a tank with a small valve at the bottom. Let the level A represent the free surface in the tank and level B be the bottom level. Let the height of the free surface with respect to the bottom be H . Bernoulli: 2 2 + ρgH = pB + 12 ρvB . pA + 21 ρvA Consider first the case, when the vent is closed and there’s no flow. Thus, pA + ρgH = pB, the standard result for hydrostatic pressure. Now, open the valve. Both pA and pB are now equal to the atmospheric pressure. For a small valve, the velocity is extremely small at the free surface (as vAAA = vBAB and AA ≫ AB), so we can take vA = 0 to a good approximation. Thus, A ρgH = 12 ρvB2 or vB2 = 2gH, which is known as Torricelli’s theorem. B Note that the velocity is independent of the direction of the outlet, so if that points upwards, it is easy to see that according to Torricelli’s theorem the water will reach the height H again. In reality, because of viscosity, the water will not quite reach the level of the surface. Hydrodynamics – Spring 2010 75 Example. Use Torricelli’s theorem to estimate, how long it takes for a tank of 1 m height to get empty through a valve in the bottom, if the diameter of the valve is 1% of the diameter of the tank. Let H(t) be the height of the surface from the bottom of the tank. The vertical component of velocity at the surface is dH/dt < 0. Conservation of mass implies q dH Torricelli = −ABvB ≈ −AB 2gH AA dt √ √ dH AB √ AB √ √ ⇒ =− 2g dt ⇒ 2( H0 − H) = 2g t(H) AA AA H Thus, for the tank to get empty (H = 0) it takes t(0) = v u Au u t v u u 2u t A 2H0 = 100 AB g Hydrodynamics – Spring 2010 A B 2m ≈ 1 h 15 min. 9.81 m s−2 76 Kelvin’s theorem on vorticity Consider vorticity ω = ∇ × v in an ideal fluid. Thus, ∂ω = ∇ × (v × ω). ∂t Circulation is defined as Γ(t) = ˛ C v · dl, where C is a closed material curve, i.e., a closed curve traversing of adjacent fluid elements and following their motion. Using Stokes’s law, we can write ˆ ˆ ω · ds, Γ(S(t), t) = (∇ × v) · ds = S S where S is (any) surface enclosed by C . Thus, circulation is equal to the flux of vorticity. Kelvin’s theorem on vorticity states that in an ideal fluid, dΓ Γ(S1, t1) − Γ(S, t) ≡ lim = 0, t1 →t dt t1 − t where S1 is the surface consisting of the positions, at time t1, of those fluid elements that were on surface S at time t. Hydrodynamics – Spring 2010 77 Proof. We will consider a slightly more ´general theorem, stating that any solenoidal (∇ · Q = 0) vector field, Q, fullfills dΓ/dt = 0, where Γ = S Q · ds, if Q evolves in time as ∂Q = ∇ × (v × Q). ∂t Calculate Γ(S1, t1) − Γ(S, t) dΓ = lim , dt δt→0 δt where t1 = t + δt. For that, consider ˆ  ˆ  ∂Q Q(x, t) + δt Q(x, t + δt) · ds ≈ Γ(S1, t1) = (x, t) · ds ∂t S1 S1 ˆ ˆ ∂Q (x, t) · ds, ≈ Q(x, t) · ds + δt ∂t S1 S where terms of higher than first order in δt have been neglected. Thus,   ˆ 1  ∂Q Γ(S1, t + δt) − Γ(S, t)  ≈  (x, t) · ds Q(x, t) · ds + Q(x, t) · (−ds) + δt δt S1 ∂t S S ˆ ˆ ˆ ˛ ˆ []= Q(x, t) · ds − Q(x, t) · ds = (∇ · Q) d3x − Q(x, t) · ds = − Q(x, t) · ds, ˆ ˆ ∂V CC1 V CC1 CC1 where CC1 is the surface drawn by the material curve as it moves from the position of C to the position of C1, and V is the volume enclosed by ∂V = S1 ∪ S ∪ CC1. Hydrodynamics – Spring 2010 78 Calculate, next, Thus, ˆ CC1 ´ Q · ds. Here, CC1 Q · ds = ˆ CC 1 Q · (−δt v(x, t) × dl) = δt ˛ (v × Q) · dl = δt ˆ [∇ × (v × Q)] · ds C S so ds = −δt v(x, t) × dl. S1 C1 S dΓ = dt   ∂Q  − ∇ × (v × Q) · ds. S ∂t ˆ This holds for any solenoidal vector field. Since we assumed the temporal evolution of Q to satisfy ∂Q = ∇ × (v × Q), ∂t C v(x,t) δt dl ds the theorem is proven. Hydrodynamics – Spring 2010 79 Ramifications of Kelvin’s theorem – Helmholtz’s theorems Like for any solenoidal field (cf. magnetic induction B ), we can introduce the concept of flux tubes of vorticity. These vortex tubes are tubes consisting of bundles of field lines of ω . The flux of vorticity across the crosssection of the tube is, therefore, constant and equal to the circulation calculated along any closed curve surrounding the vortex tube. C3 ω C2 Helmholtz’s theorems state 1. The strength of a vortex tube (i.e., the circulation associated with it) does not vary in time 2. Fluid elements lying on a vortex line at some instant continue to lie on that vortex line. More simply, vortex lines move with the fluid. 3. Fluid elements initially free of vorticity remain free of vorticity. All these theorems are corollaries of Kelvin’s theorem. C1 Also, because of the solenoidal property of vorticity, ∇ · ω = 0, vortex lines and tubes must appear as closed loops, extend to infinity or start/end at solid boundaries. Hydrodynamics – Spring 2010 80 Potential flow Consider an incompressional ideal fluid with ω = 0 everywhere at some instant of time. Because of Helmholtz theorems, we then know that ω = 0 everywhere at all instants of time! The fluid is said to be irrotational . Thus, as ∇ × v = 0 everywhere at all times, we know that v has to be expressible as a gradient of a potential, v = −∇φ. Furthermore, as the fluid was assumed incompressible, we have 0 = −∇ · v = ∇ · (∇φ) = ∇2φ, i.e., velocity potential φ fulfills the Laplace equation. Since the normal component of velocity, vn = −en · ∇φ, has to vanish at a fixed boundary of the fluid, the appropriate boundary condition is en · ∇φ = ∂φ = 0. ∂n It is well known that Laplace equation with such Neumann-type boundary conditions has a unique solution. Notice that as potential theory (i.e., Laplace’s equation) is applied also in electro- and magnetostatic problems, many familiar results of electrodynamics can be immediately borrowed to fluid mechanics of irrotational flows. Below we will consider some of the classical solutions of the irrotational, inviscid and incompressible flows solved by potential theory. Hydrodynamics – Spring 2010 81 Sources and sinks In analogy with the electrostatic potential of a point charge, φES = (q/4πε0)r−1, we express the velocity potential as φ = (m/4π)r−1, where r denotes the distance from the point ◦, and m is constant. Clearly, ∇2φ = 0, when r 6= 0. Thus, v = −∇φ = m er , 4πr2 i.e., the flow is radially outwards from ◦, and the mass flow [kg s−1] through any spherical surface centered around ◦ is 4πr2ρ0vr = mρ0. Thus, the singular point ◦ is a source and m represents its strength. A sink is a negative source (denoted by •). If a source is a point of creation of the fluid, then sink is a point of annihilation. The flow is directed radially inwards and the potential (for a sink of strength m) is given by φ = −(m/4π)r−1. Hydrodynamics – Spring 2010 82 Flow near an infinite plate The method of images, familiar from electrostatics, can be used to determine the velocity potential near an infinite plate. Consider a point source of strength m at x = −aex and the plate at x = 0. The effect of the plate can be obtained by placing an image source of the same strength at x = aex. Thus,   1 m 1 + φ= 4π  r1 r2  r2 r1 m source x O m image with r1 = [(x + a)2 + y 2 + z 2]1/2 and r2 = [(x − a)2 + y 2 + z 2]1/2 being the distance from the source and the image, respectively. Therefore, vx = −   ∂φ m dr1 m dr2 m x + a x − a  = + = + 3  ∂x 4πr12 dx 4πr22 dx 4π r13 r2 Thus, at x = 0, vx = 0 and there’s no flow along the x direction. The surface of the plate is a stream surface. Hydrodynamics – Spring 2010 83 Flow past a cylinder Consider the flow past a fixed infinitely long cylinder (of radius a) with its axis perpendicular to the flow direction. At large distances from the cylinder, assume the flow to be uniform, i.e., v = −U ex. Determine the velocity field around the cylinder. Clearly, we have a translational symmetry (along the axis of the cylinder) in the problem. Thus, it is two-dimensional. Furthermore, the boundary condition at the surface of the cylinder is most easily formulated in polar coordinates (r, θ): ∂φ = 0 at r = a. ∂r It is, therefore, natural to make use of polar coordinates, where the Laplacian can be given as 1 ∂  ∂φ  1 ∂ 2φ r + 2 2. ∇ φ(r, θ) = r ∂r ∂r r ∂θ  2 Thus, write  ∂  ∂φ  ∂ 2φ r r =− 2 ∂r ∂r ∂θ and separate the variables, i.e., write φ = R(r)Θ(θ). Thus,   1 d2Θ r d  dR  . r = −K = − R dr dr Θ dθ2  Hydrodynamics – Spring 2010  84 Solve, first, the latter equation √ −θ K ′′ Θ − KΘ = 0 ⇒ Θ = a−e √ θ K + a+ e and note that periodicity Θ(θ) = Θ(θ + 2π) has to be required. Thus, needs to be adopted so K = −n2. Hence, √ K = i n, where n is an integer, Θn(θ) = Cn cos nθ + Dn sin nθ is the general solution for the angular part. Thus, the general solution for the Laplace equation is φ(r, θ) = where ∞ X n=0 Rn(r)[Cn cos nθ + Dn sin nθ],   d dRn  r r − n2Rn = 0. dr dr Proceeding to the radial part, we note first that n = 0 gives Θ0 = C0 and r dR0 = B0 = constant. dr which satisfies the boundary condition R′(a) = 0 only for B0 = 0 ⇒ R0(r) = A0 = constant. As constants can always be added to potentials without changing their gradients, we can set A0 = 0. Hydrodynamics – Spring 2010 85 For n > 0, writing y = ln r we get d/dr = (dy/dr) d/dy = (1/r) d/dy ⇒ r d/dr = d/dy so d2Rn − n 2 Rn = 0 2 dy This gives and ⇔ Rn = Aneny + Bne−ny = Anrn +  Bn . rn  Bn  r n + φ(r, θ) = (Cn cos nθ + Dn sin nθ) rn n=1 ∞ X   ∞ ∂φ Bn X n rn−1 − n+1  (Cn cos nθ + Dn sin nθ). = ∂r n=1 r The inner boundary condition then gives Bn = a2n ∀n = 1, 2, 3, ... so φ(r, θ) = ∞ X n=1    rn + 2n  a   (Cn cos nθ + Dn sin nθ) rn At r → ∞, v = −U ex so φ = U x = U r cos θ there. Hence, Dn = 0 ∀n, C1 = U and Cn = 0 for n = 2, 3, ... and finally   2 a  a2  φ(r, θ) = U cos θ r +  = U x + U cos θ, r r with ∇ = er ∂r + eθ r−1∂θ , a2 v = −U ex + U 2 (cos θ er + sin θ eθ ). r Hydrodynamics – Spring 2010 86 Streamlines of the flow past a cylinder Hydrodynamics – Spring 2010 87 Velocity field around a moving cylinder We showed that for a liquid streaming uniformly around a fixed cylinder, the velocity is a2 v = −U ex + U 2 (cos θ er + sin θ eθ ) r Make a Galilean transformation into a coordinate system moving with constant velocity −U ex. In that frame of reference, the fluid is at rest at r → ∞ and the cylinder moves at velocity U ex. Thus, a2 v = v + U ex = U 2 (cos θ er + sin θ eθ ) r ′ is the fluid velocity in that frame. Here the radial distance r and the angle θ are measured from and around the axis of the moving cylinder. Thus, |v ′| = U (a/r)2 is independent on the angle measured around the cylinder. This equation actually applies even if the cylinder is in a non-uniform 2-D motion (i.e., if U is a function of time but stays perpendicular to the axis of the cylinder). This is because we rassumed the fluid to be incompressible. An incompressible fluid has a an infinite sound speed (given by ∂p/∂ρ), implying that the fluid reacts to the changes in the motional state of the cylinder instantaneously. Hydrodynamics – Spring 2010 88 Equation of motion of the cylinder. d’Alembert’s paradox The kinetic energy of the fluid per unit length of the cylinder is ˆ ∞ 2 v ′ 2πr dr = 21 πρa2U 2 = 21 M ′U 2, Kfluid = 21 ρ a where M ′ = πρa2 is the mass of fluid displaced by a unit length of the cylider. If M is the mass of the cylinder per unit length, then the total kinetic energy associated with the unit length of the system is Ktot = 12 (M + M ′)U 2. It is now possible to to deduce the equation of motion of the cylinder. If F is the force acting on the unit length of the cylinder, then the rate of work done by the cylinder is F · U = (M + M ′) dU dU · U ⇒ F = (M + M ′) dt dt Looks just like Newton II for the cylinder in vacuum, but with M replaced by M + M ′, which is referred to as the effective mass of the cylinder. One also notes that if dU /dt = 0 there is no resultant force acting on the cylinder. This is contrary to our everyday experience and known as d’Alemert’s paradox. The result arises, of course, since we neglected viscosity from our calculation, and the fluid may flow freely on the surface of the cylinder in our model. Hydrodynamics – Spring 2010 89 Cylinder immersed in an ideal fluid in a gravitational field Consider an infite cylinder immersed in a fluid with its axis perpendicular to the vertical direction. Assume that a gravitational field of g = −gez is acting in the system. Since the fluid is ideal, we can use the hydrostatic result for the force F = −(M − M ′) g ez . Assume that the cylinder is intitally at rest. Thus, U = −U ez and the equation of motion for the cylinder is dU = (M − M ′) g (M + M ′) dt giving a result that the cylinder experiences a constant acceleration, M − M′ dU g. = dt M + M′ This, again, is against our everyday experience that states that a body falling in a fluid will approach asymptotically a terminal velocity. The discrepancy is again, of course, due to the neglect of viscosity. Hydrodynamics – Spring 2010 90 Stream function The equation of continuity for an incompressible fluid in two dimensions gives ∂vx ∂vy + = 0. ∂x ∂y In general, a dx + b dy is an exact differential iff ∂a/∂y = ∂b/∂x. Thus, vy dx − vx dy is an exact differential, say dψ , and ˆ ˆ ψ = dψ = (vy dx − vx dy) is called the stream function with ∂ψ = vy ; ∂x ∂ψ = −vx. ∂y Recall that the streamlines are given by dx dy = ⇒ vx vy 0 = vy dx − vx dy = dψ, i.e., ψ is constant along streamlines. Hydrodynamics – Spring 2010 91 Relation between stream function and velocity potential If the flow is irrotational,   ∂vy ∂vx  0 = ∇ × v = ez  , − ∂x ∂y then using vx = −∂ψ/∂y and vy = ∂ψ/∂x gives ∂ 2ψ ∂ 2ψ + 2 = ∇2ψ, 0= 2 ∂x ∂y i.e., ψ satisfies Laplace’s equation3 just as φ in v = −∇φ. We get ∂ψ ∂φ = vy = − ∂x ∂y and ∂ψ ∂φ = −vx = ∂y ∂x so the functions ψ and φ fulfill the Cauchy–Riemann conditions. As ∇ψ · ∇φ = −vy vx + vxvy = 0 the curves ψ = const. (stream lines) and φ = const. (equipotential lines) always intersect at right angles. Because of the Cauchy–Riemann conditions, the functions ψ and φ are conjugate functions, i.e., the real and imaginary parts of a single analytic complex function, called the complex potential . 3 Functions satisying Laplace’s equation are called harmonic. Hydrodynamics – Spring 2010 92 Complex potential Because they are conjugate function, the stream function and the velocity potential can be combined to a complex potential by treating the xy plane as the complex plane: f (z) = φ(x, y) + iψ(x, y), z = x + iy The derivative of this analytic function, ∂φ ∂ψ df = +i = −(vx − ivy ) ≡ −w dz ∂x ∂x where w is called the complex velocity. Thus, given the complex potential one can easily determine the the complex velocity and the flow velocity. The fact that a 2-D irrotational flow can be obtained as a derivative of an analytic complex potential is important from the point of view of applications, because it allows us to produce solutions to complicated flow problems from simple solutions using conformal mapping. Hydrodynamics – Spring 2010 93 Some simple flows Consider, first, f (z) = Az n = A(reiϕ)n = Arn(cos nϕ + i sin nϕ), where A ∈ R is a constant. Thus, φ = Arn cos nϕ; ψ = Arn sin nϕ; Examples: n = 1; A = U ; df = −nAz n−1. w = − dz φ = U r cos ϕ = U x; ψ = U r sin ϕ = U y; w = −U rectilinear flow in the − x direction n = 2 : φ = Ar2 cos 2ϕ = A(x2 − y 2); n=1 2 ψ = Ar sin 2ϕ = 2xy; n=2 vx − ivy = w = −2Az = −2A(x + iy) ⇒ vx = −2Ax; vy = 2Ay hyperbolic flow around a rectangle n = −1; A = U a2 : φ+1 + φ−1 = φcyl a2 a2 a2 φ = U cos ϕ; ψ = −U sin ϕ; w = U 2 r r z   a2   ⇒ f (z) = U z +  = complex potential of a flow around a cylinder z   a2   ψ = U sin ϕ r −  = stream function of a flow around a cylinder r Hydrodynamics – Spring 2010 94 Conformal mapping Assume that an analytic complex function g(z) = u(x, y) + iv(x, y) is given, and that the derivative of g is non-zero in a region of the xy plane. This function serves as a mapping between two complex regions, i.e., (x, y) ↔ (u, v). Such a mapping is called a conformal mapping.4 Assume that the potential φ̄(u, v) is a known solution of a simple flow problem, i.e., that it is a harmonic function of u and v and satisfies certain boundary conditions in a region of the uv plane. Thus, φ(x, y) = φ̄{u(x, y), v(x, y)} is a harmonic function of x and y and represents a solution to flow problem in the region of the xy plane that maps to the region of the flow in uv plane. This may be applied, of course, also to the stream function. Thus, if the streamlines in the uv -plane are given by ψ̄(u, v) = const. the relation ψ(x, y) = ψ̄{u(x, y), v(x, y)} = const. gives them in the xy -plane. 4 Note that u and v are just a transformed coordinates in the complex plane and have nothing to do with velocity. Hydrodynamics – Spring 2010 95 Example of conformal mapping Consider the rectilinear flow in the region u ∈ R, v ∈ [−π, π], described by φ̄(u, v) = −V u; ψ̄(u, v) = −V v. Consider a conformal mapping given by the analytic function g(w) as z = g(w) = w + ew = u + iv + eu(cos v + i sin v) = u + eu cos v + i(v + eu sin v)     x = u + eu cos v ∴  y = v + eu sin v The stream lines in the uv plane are given by v = const. Thus, in xy plane the streamlines are obtained in parametric form by fixing the value of v and regarding u as a parameter in the equations of the conformal mapping. E.g., take v = 0 ⇒ y = 0 and x = u + eu, i.e., the x axis. v = ±π gives y = ±π and x = u − eu, i.e., a horizontal flow line but with x → −∞ with both u → ±∞. As ∂x/∂u = 1 − eu = 0 ⇔ u = 0 ⇔ x = −1, the streamlines thus make a 180◦ turn at x = −1. At other values of v , we get x u + eu cos v = → cot v, y v + eu sin v as u → ∞, i.e., the streamlines become straight lines at large u. This solution describes the flow from a channel. Hydrodynamics – Spring 2010 96 The equipotential lines are obtained by keeping u fixed and letting v run from −π to π in     x = u + eu cos v y    = v + eu sin v For −u ≫ 1 they are just vertical lines at x = u between y = −π and π . At u ≫ 1, they become circles. π v y π u −π Hydrodynamics – Spring 2010 x −π 97 Viscous flows Hydrodynamics – Spring 2010 98 Newtonian fluids In our microscopic derivation of the pressure tensor, we found Pij = p δij + πij , where πij was found to be a traceless tensor linearly dependent on Λij = 21 (∂ivj + ∂j vi). In ideal fluids, we assumed πij = 0. This assumption was pointed out on several occasions to be erroneous and lead to unphysical results. Derive, next, the physical result for πij from macroscopic considerations. Internal friction of a fluid, called viscosity, opposes relative motion of adjacent layers of fluid . For a velocity field with shear, like the one in the Figure to the right, we know that the slower-moving fluid below the xz plane will act on the faster moving fluid above the xz plane (and vice versa) to slow down (speed up) its motion. The force has to act across the xz plane (normal in the y direction), but the force has to be in the x direction to level out the speed differences. As the surface force through dS = dS ey is (dFsurf )i = −p dSi − πij dSj , y x we note that the x component of the force has to be due to πxy . The shear force is expected to be larger for a larger velocity gradient. Hydrodynamics – Spring 2010 99 Newton postulated that the shear force is proportional to the velocity gradient, i.e., πxy = −µ ∂vx , ∂y (N) where µ is the coefficient of viscosity. Fluids obeying the proportionality relation between the shear stress and the velocity gradient are known as Newtonian fluids. Shear stress in the general situation has to be more complicated than the simple relation above. This can be seen, e.g., by considering a fluid rotating rigidly. Consider, for example, a body rotating with angular velocity Ω = Ωez : v = Ω × x = Ω(xey − yex) and −µ∂y vx = µ∂xvy = µΩ, so clearly the expression (N) is non-zero although there is no shear involved with rigid rotation. However, µ(∂y vx + ∂xvy ) = 0. In general, we can write     ∂vi 1 ∂vi ∂vj  1  ∂vi ∂vj  1 +  = Λij + ǫijk ωk   =  + − 2 ∂xj 2 ∂xj ∂xi 2 ∂xj ∂xi where ωk = ǫkrs∂r vs is the vorticity. Here, Λij is related to the rate of deformation and 12 ǫijk ωk to rotation. Thus, the shear stress of a Newtonian fluid should depend on the symmetric part of the velocity gradient, Λij , rather than the velocity gradient itself. The most general second rank tensor depending linearly on symmetric combinations of velocity gradients is   ∂vk ∂vi ∂vj  +b  + δij = 2aΛij + bΛkk δij . a ∂xj ∂xi ∂xk Hydrodynamics – Spring 2010 100 Therefore, a Newtonian fluid is a fluid obeying the constitutive equation  Pij = p δij − 2µ Λij − 31 δij Λkk      ∂vi ∂vj 2 ∂vk  ∂vk    δij − µ  + − δij − ζΛkk δij = p − ζ ∂xk ∂xj ∂xi 3 ∂xk i.e., P = p̄I + π where p̄ = p − ζtr Λ and π = −2µ(Λ − 13 I tr Λ) Here π is a traceless tensor corresponding to viscous stress due to velocity shear. The coefficient µ is the dynamic viscosity (or shear viscocity or viscosity) already obtained through our microscopic derivation. The first term, i.e., p̄I, is the spherical part of the P tensor. Here p̄ = 31 Pkk = p − ζ∇ · v = p + ζ 1 dρ . ρ dt The coefficient ζ is called the bulk viscosity (or volume viscosity or second viscosity). It is clearly related to viscous forces involving compressive motions. Bulk viscosity is difficult to determine experimentally, and important only in situations when compression is important (e.g., for shock waves and damping of sound waves). According to our microscopic derivation, bulk viscocity vanishes for dilute gases, but large values are experimentally obtaind for many liquids (Bhatia & Singh). Liquid Glycerol Water Methyl alcohol Toluene Benzene ζ/µ 1.1 2.5 3.2 13 100 Hydrodynamics – Spring 2010 101 Navier–Stokes equation The surface force density for P = p̄I + π, where p̄ = p − ζtr Λ and π = −2µ(Λ − 13 I tr Λ), is   ∂p ∂ ∂  ∂vi ∂vj 2 ∂vk  ∂Pij   =− +ζ (∇ · v) + µ + − δij − ∂xj ∂xi ∂xi ∂xj ∂xj ∂xi 3 ∂xk =− ∂p ∂ + µ∇2vi + (ζ + 13 µ) (∇ · v) ∂xi ∂xi i.e., −∇ · P = −∇p + (ζ + 31 µ)∇(∇ · v) + µ∇2v . Here, the viscocity coefficients are assumed to be constants. Thus, the equation of motion for a Newtonian fluid is ρ dv = ρF − ∇p + (ζ + 31 µ)∇(∇ · v) + µ∇2v. dt This is, of course, the Navier–Stokes equation. The only difference with respect to the microscopic result derived for dilute gases is the inclusion of the bulk viscocity term. As noted already in connection with the microscopic derivation, for many problems it is a reasonable approximation to neglect compressibility effects. Hence, a simpler version of the Navier–Stokes equation is obtained as 1 ∂v + v · ∇v = F − ∇p + ν∇2v, ∂t ρ where ν = µ/ρ is called the kinematic viscosity. Hydrodynamics – Spring 2010 102 Navier–Stokes vs. Euler equation The N–S equation in form ∂v 1 + v · ∇v = F − ∇p + ν∇2v ∂t ρ differs from Euler equation only by virtue of one extra term, ν∇2v , but this changes everything! The additional term contains second derivatives in spatial variables, whereas the Euler equation only contains first-order terms. Thus, more boundary conditions are required to solve the N–S equation than the Euler equation. For Euler equation, one takes the normal component of the velocity to vanish at fixed solid boundaries and obtains a unique solution. It is impossible, in general, to impose more boundary conditions, which led to the several unphysical results highlighted in the previous section. The N–S equation with the second derivatives now allows us to impose the viscous-fluid boundary condition, v = 0 at the solid surfaces, which allows one, for example, to consider the effect of friction on bodies immersed in fluids. But why do we actually impose v = 0 at solid surfaces? Experience tells us that this is the way to go. For example, blades of a fan collect dust. It is also very difficult to sweep surfaces clean from fine dust by blowing from one’s mouth. Hydrodynamics – Spring 2010 103 Vorticity equation for Newtonian fluids Taking a curl of the N–S equation gives ∂ω = ∇ × (v × ω) + ν∇2ω ∂t for incompressible or barotropic fluids with conservative external forces. Kelvin’s theorem on vorticity does not apply anymore and instead we get ˆ ˛ ˆ   ∂ω dΓ  = − ∇ × (v × ω) · ds = ν (∇2ω) · ds = ν (∇2v) · dl dt S C S ∂t where C = ∂S and the last form holds only for incompressible fluids. Thus, a finite viscosity allows the production and decay of vorticity. This allows one to overcome some problems that are associated with ideal fluids. For example, if a stick is moved through a fluid, like water, we can clearly see the generation of vortices behind the moving stick. This can only be understood if the fluid has non-zero viscosity. For incompressible viscous fluids, energy equation is redundant. A dynamical theory of the fluid is provided by the vorticity equation, definition of vorticity, ω = ∇ × v , and the continuity equation, ∇ · v = 0. Hydrodynamics – Spring 2010 104 Navier–Stokes equation in curvilinear coordinates In many flow problems, the symmetries of the boundaries or solid bodies immersed in the flow suggest the use of some curvilinear cooridnates system, e.g., cylindrical or spherical coordinates. Thus, it is necessary to express the terms of the N–S equation in these coordinate systems. Differential operators acting on vectors and tensors are quite complicated in curviliear cooridnates. The N–S equation contains the terms v · ∇v and ν∇2v , which need special care. It is usually easiest to use the identities v · ∇v = 21 ∇v 2 − v × (∇ × v) ∇2v = ∇(∇ · v) − ∇ × (∇ × v) before writing down the N–S equation in component form. In cylindrical coordinates, for example, one gets   ∂vr ∂ 2 vr 1 ∂ 2 vr ∂ 2 vr 1 ∂vr 1 ∂p 2 ∂vθ vr  ∂vr vθ ∂vr ∂vr vθ2 + + vr + + vz − = Fr − +ν 2 + 2 2 + − − 2 ∂t ∂r r ∂θ ∂z r ρ ∂r ∂r r ∂θ ∂z 2 r ∂r r2 ∂θ r   ∂v ∂v v ∂v v vr 1 ∂p 1 ∂ 2v ∂ 2 vθ 1 ∂vθ 2 ∂vr vθ  ∂ 2v ∂vθ + vr θ + θ θ + vz θ − θ = Fθ − + ν  2θ + 2 2θ + + + − 2 ∂t ∂r r ∂θ ∂z r ρr ∂θ ∂r r ∂θ ∂z 2 r ∂r r2 ∂θ r   ∂vz ∂ 2 vz 1 ∂ 2 vz ∂ 2 vz 1 ∂vz  1 ∂p ∂vz vθ ∂vz ∂vz . + + vr + + vz = Fz − +ν 2 + 2 2 + ∂t ∂r r ∂θ ∂z ρ ∂z ∂r r ∂θ ∂z 2 r ∂r Hydrodynamics – Spring 2010 105 Flow through a circular pipe Consider the steady-state flow in a straight pipe with circular cross-section (radius a). Take the flow to be incompressible and external forces to vanish. Take the flow velocity to be along the axis of the tube, i.e., v = vez . For symmetry reasons v cannot depend on the angle measured around the cylider and because of the incompressibility (uniform cross section), it cannot depend on z . Thus, v = v(r)ez and Therefore, ∂ dv ∂v = + v · ∇v = v (v ez ) = 0 dt ∂t ∂z   µ ∂  ∂v  ∇p = µ∇2v = r ez . r ∂r ∂r Since ∇p k ez , pressure is a function z , only. Since v depends on r, only, we have   dp µ d  dv  = r , dz r dr dr which is of course possible only if both sides are equal to a constant. From the left-hand side, this constant is seen to be (p2 − p1)/l = −∆p/l, where p1 and p2 are the pressure at z = 0 and z = l, respectively. Hydrodynamics – Spring 2010 106 We get   ∆p dv ∆p 2 ∆p 2 d  dv  r = − r ⇒ r = c1 − r ⇒ v = c0 + c1 ln r − r, dr dr µl dr 2µl 4µl where ci are constants of integration. As the velocity must remain finite at r = 0, c1 = 0, and because v = 0 at r = a, we get ∆p 2 a 4µl c0 = ⇒ v= ∆p 2 (a − r2) 4µl The maximum speed is, therefore, vmax = ∆p a2/4µl. The average speed is ´a ´a 2 2 v 2πr dr v vmax z max 0 r(a − r )dr 0 ´a vav = ´ a = 2 = a 2 0 2πr dr 0 r dr and the rate of mass flow through the cross-sectional area of the tube, or the discharge, is Q= ˆ 0 a ρvmax πa2 πa4∆p = . ρvz 2πr dr = 2 8νl Thus, the discharge is proportional to the fourth power of the radius of the tube. This result is known as the Hagen–Poiseuille formula. The Hagen–Poiseuille formula provides means to determine the viscosity of a fluid by measuring Q vs. ∆p. Note that the result is experimentally valid only up to a certain velocity vmax after which the flow becomes turbulent. When this occurs depends on the Reynolds number of the flow. Hydrodynamics – Spring 2010 107 Couette flow Consider incompressible fluid between two parallel plates at x = 0 and x = h. Assume that the upper plate is moving at speed V along the z axis while the lower plate is at rest. Now, the flow between the plates can be assumed to be along the z axis and independent on coordinates y and z . Thus, v = v(x) ez and dv ∂v = + v · ∇v = 0 dt ∂t so ∂ 2v 2 ∇p = µ∇ v = µ 2 ez , ∂x where, again, the pressure can depend only on z . Thus, d2v dp =µ 2 dz dx where both sides must be equal to a constant. This integrates to give x2 dp v = c0 + c1 x + , 2µ dz where dp/dz is a constant. Hydrodynamics – Spring 2010 108 Since v = 0 at x = 0, c0 = 0, and since v = V at x = h, we get the result h2 dp hc1 = V − 2µ dz ⇒ V x x2 − hx dp + v= h 2µ dz known as the generalized Couette flow. (Simple Couette flow is obtained by setting p = const.) Exercise: sketch the velocity profile for dp/dz < 0 and dp/dz > 0. What is the condition for zero discharge? The drag force excerted by the fluid on the plates (per unit area of the plate) is given by   V µ 2x − h dp  dv = ∓ + , ∓µ dx x=(h±h)/2 h 2 dz x=(h±h)/2 where the upper [lower] signs correspond to the upper [lower] plate. The force per unit area on the lower plate is, then, V µ h dp − h 2 dz and on the upper plate V µ h dp − . − h 2 dz In simple Couette flow (dp/dz = 0), the drag forces are the same but oppositely directed. Hydrodynamics – Spring 2010 109 Scaling and Reynolds number Suppose a new submarine is designed and we want to find out the flow pattern around it when it is moving through water. Is it possible to do laboratory experiments with a miniatyre model to find out, e.g., the drag force excerted by the water on the submarine moving at different speeds? In short, yes. Consider the fluid flow around geometrically similar objects. Let L be a typical size of the object and V be a typical velocity of the fluid flow. Thus, L/V can be taken as the typical time scale of the problem. Let us scale all variables by their typical values, x = x′L; v = v′V ; t = t′L/V ; ω = ω ′V /L We can see that the primed quantities are all dimensionless. The vorticity equation can be written in a dimensionless form as 1 ′2 ′ ∂ω ′ ′ ′ ′ = ∇ × (v × ω ) + ∇ ω ∂t′ Re where ∇′ = ∂/∂x′ = L∇ and Re = LV ν is a dimensionless number known as the Reynolds number of the flow. For two geometrically similar flows, the dimensionless flow patterns are identical, if their Reynolds numbers are the same. Hydrodynamics – Spring 2010 110 Examples of scaling • Uniform flow past a cylinder or a sphere: take L = 2a and V = U (relative speed of the fluid and the object) 2aU . ⇒ Re = ν Thus, for miniature modeling of submarines we need to consider faster-than-real flows or fluids less viscous than water (at ∼ 10◦C) in order to achieve a Reynolds number representative of the true situation. • Hagen–Poiseuille flow: take L = 2a and V = vave = ∆p a2/8µl ∆p a3 . ⇒ Re = ρl 4ν 2 Experimentally, laminar flow is obtained as long Re . 2400. If geometrical similarity is maintained (i.e., a/l kept constant) and the Reynold’s number is kept constant, then ∆p a2 ρ 4ν 2 is also constant. This dimensionless quantity equals Eu/(4Re)2, where Eu ≡ ∆p/ρV 2 is the Euler number of the flow. Thus, Hagen–Poiseuille flow is characterised by two dimensionless numbers picked from Re, Eu and a/l. Hydrodynamics – Spring 2010 111 Slow viscous flow past solid bodies For incompressional (and barotropic) fluids, we derived the dimensionless equation 1 ′2 ′ ∂ω ′ ′ ′ ′ ∇ ω. = ∇ × (v × ω ) + ∂t′ Re Clearly, at the limit of Re ≪ 1, in steady state the equation reduces to ∇2 ω = 0 which is a much simpler equation than N–S because it is linear. Assuming incompressibility, ∇ · v = 0, we can write v = ∇ × ψ , where ψ is the vector potential of the flow. In a two-dimensional flow, we can take ψ = −ψ(x, y)ez , i.e., v=− ∂ψ ∂ψ ex + ey ∂y ∂x showing that ψ is just the streamfunction. Use, next, ω = ∇ × (∇ × ψ) = ∇(∇ · ψ) − ∇2ψ , where the first term vanishes in the 2-D case. Thus, the stream function fullfils the bi-harmonic equation ∇4 ψ = 0 for a viscous fluid. The correct boundary condition at the surfaces is ∇ψ = 0. Hydrodynamics – Spring 2010 112 Stokes flow past a sphere For axisymmetric meridional flows, like the flow past a sphere, we can use the Stokes streamfunction Ψ defined, in spherical coordinates, by ψ = (Ψ/r sin θ) eϕ: vr = 1 ∂Ψ 1 ∂(sin θ ψϕ) = 2 ; r sin θ ∂θ r sin θ ∂θ vθ = − Thus, (exercise)  1 ∂(rψϕ) 1 ∂Ψ =− r ∂r r sin θ ∂r 2 2 sin θ ∂  1 ∂   ∂ 2  Ψ = 0. Ê Ψ ≡  2 + 2 ∂r r ∂θ sin θ ∂θ  Far from the sphere, the velocity components for a uniform flow along the z axis are ∂Ψ 1 ∂Ψ ⇒ = rU sin2 θ ⇒ Ψ = c(θ) + 12 U r2 sin2 θ r sin θ ∂r ∂r c′(θ) 1 ∂Ψ c′(θ) + U r2 sin θ cos θ = = 2 + U cos θ ⇒ c′(θ) = 0 U cos θ = vr = 2 2 r sin θ ∂θ r sin θ r sin θ −U sin θ = vθ = − ⇒ Ψ = 12 U r2 sin2 θ at r ∼ ∞ The boundary conditions on the sphere are vr = vθ = 0 at r = a, and the boundary conditions at infitinity are Ψ ∼ 12 U r2 sin2 θ. Let us try a solution of form Ψ = f (r) sin2 θ. Thus,  2  2 2 2 d d 2 2 2f    2 2 2 2 ′′  f (r) = 0. ÊΨ = f − 2  sin θ ⇒ Ê Ψ =  2 − 2  f (r) sin θ ⇒ D̂ f ≡  2 − 2  r dr r dr r   Hydrodynamics – Spring 2010 113 Seek for solutions of form f = rn. Thus, D̂f = [n(n − 1) − 2]rn−2 and D̂2f = [n(n − 1) − 2][(n − 2)(n − 3) − 2]rn−4 so n must fulfill [n(n − 1) − 2][(n − 2)(n − 3) − 2] = 0 ⇔ n = 1± √ √ 1± 1+8 1+8 = −1, 2 ∨ n = 2 + = 1, 4 2 2 Thus, n = −1, 1, 2, 4 and f (r) = Ar−1 + Br + Cr2 + Dr4. Now, at infinity, f ∼ 21 U r2 so D = 0 and C = 21 U . At the surface, f ′(a) = 0 = f (a) so     −Aa −2 + B + Ua =0 Aa−1 + Ba + 12 U a2 = 0        A = 41 U a3 ⇒  B = − 43 U a Hence, finally, Ψ = 12 U (r2 − 23 ar + 21 a3r−1) sin2 θ and 1 ∂Ψ = U (1 − 32 ar−1 + 12 a3r−3) cos θ 2 r sin θ ∂θ 1 ∂Ψ = −U (1 − 43 ar−1 − 41 a3r−3) sin θ vθ = − r sin θ ∂r vr = Hydrodynamics – Spring 2010 114 Streamlines of Stokes flow In each meridional plane, ϕ = constant, 1 ∂Ψ dΨ 1 ∂Ψ dr r dθ dr + 2 r dθ = 0 ⇒ =0 = ⇒ vr vθ r sin θ ∂r r sin θ ∂θ r sin θ so the streamlines are contours of the Stokes streamfunction, i.e., (r2 − 23 ar + 21 a3r−1) sin2 θ = b2, where b is the distance of the streamline from the z axis at r → ∞. Thus, the flow is symmetric with respect to equator, θ = 12 π . Hydrodynamics – Spring 2010 115 Stresses. Stokes’s law Stresses (incompressible flow): −P = −pI + 2µΛ, where Λ = 12 [∇v + (∇v)T], ∇v is the velocity gradient tensor, and where p can be solved from the linearized equation of motion (consistent with Stokes flow), 0 = −∇p + µ∇2v = −∇p − µ∇ × ω. Vorticity:   3 sin θ 1 ∂(rvθ ) 1 ∂vr   eϕ = − U a  − eϕ ω= r ∂r r ∂θ 2 r2 Pressure distribution: Write the equation of motion in component form as 1 ∂ 3U aµ cos θ ∂p 3U aµ cos θ (sin θωϕ) = = −µ ⇒ p = p (θ) − 0 ∂r r sin θ ∂θ r3 2r2 3U aµ sin θ 1∂ 1 ∂p (rωϕ) = =µ ⇒ p′0(θ) = 0 3 r ∂θ r ∂r 2r 3U aµ cos θ ⇒ p = p0 − 2r2 Hydrodynamics – Spring 2010 116 Viscous stresses: The components of the viscous stress −π = 2µΛ = µ[∇v + (∇v)T] are most easily given in spherical coordinates. For axisymmetric meridional flows, v = vr (r, θ)er + vθ (r, θ)eθ , so the gradient operator can be given as ∇ = er ∂r + r−1eθ ∂θ . Thus, ∇v = (er ∂r + r−1eθ ∂θ )[vr (r, θ)er (θ) + vθ (r, θ)eθ (θ)] e) = er (er ∂r vr + eθ ∂r vθ ) + r−1eθ (er ∂θ vr + eθ ∂θ vθ + vr ∂ e +vθ ∂ | θ{z θ} | θ{z r} =eθ =−er = er er ∂r vr + er eθ ∂r vθ + eθ er r−1(∂θ vr − vθ ) + eθ eθ r−1(∂θ vθ + vr ) ∴     ∂vr 1 ∂vθ vr  1 ∂vθ 1 ∂vr vθ  Λ = er er + eθ eθ  + + − + 2 (er eθ + eθ er )  ∂r r ∂θ r ∂r r ∂θ r Drag force: Total force exerted by the fluid on the sphere is ˆ ˆ F D = (−P) · ds = (−pI − π) · ds, S S where the integral covers the whole surface of the sphere. At the surface of the sphere, r = a and ds = er a2 sin θ dθ dϕ so ˆ π ˆ 2π F D = −a2 dθ sin θ dϕ [(p + πrr ) er + πθr eθ ] . 0 Hydrodynamics – Spring 2010 0 117 The relevant components of π are, thus, ∂vr =0 ∂r   1 ∂v v 3U µ ∂v θ θ r πθr = −µ  + − = sin θ. r ∂θ ∂r r 2a πrr = −2µ Symmetry ⇒ non-zero component of force only in the z -direction. Therefore, 2 F D = −a ez 2 ˆ π dθ sin θ 0 = −2πa ez ˆ 0 ˆ 0 π 2π dϕ[(p + πrr ) |er {z· ez} +πrθ e ·e ] | θ {z z} =cos θ =− sin θ dθ sin θ(p cos θ − πrθ sin θ) = 6πU µa ez Stokes’s law Thus, a drag proportional to viscosity, velocity and circumference of the sphere is obtained. Stokes’s law is obtained from ∇2ω = 0, i.e, neglecting the non-linear term v · ∇v from N–S equation because of its assumed smallness relative to the linear viscous term ν∇2v for small Re. However, far from the sphere v ≈ U so v · ∇v ≈ U · ∇v = U · ∇(v − U ), and |v − U | ∼ U ar−1. Thus, |v · ∇v| ∼ U 2ar−2. On the other hand, 1 U 2 a2 1 a U aν |v · ∇v| ∼ ν|∇ v| = ν|∇ (v − U )| ∼ 3 = r Re r3 Re r 2 2 so the neglected term is larger than the viscous term at a/r . Re. This leads to a correction term in the drag force at Reynolds numbers not much smaller than unity. Hydrodynamics – Spring 2010 118 Flow at larger Reynolds numbers One could expect that at large Re, the viscocity term could be neglected and the flow around a solid body, like a cylinder, would be identical to that of an ideal flow. The real situation, found experimentally, is quite different. When Re & 20, two vortices appear behind the cylinder; when Re is further increased, the vortices on the two sides appear alternately, and a string of vortices, called Kármán vortex street (von Kármán 1911), appears in the wake of the cylinder. Finally, when Re & 104 the wake becomes turbulent. Hydrodynamics – Spring 2010 119 Drag at large Reynolds numbers In general, drag on a solid body can be written as FD = 21 CDAρU 2 where A is the reference area and CD is the drag coefficient of the body that for Stokes flow around a sphere are given by 24 12µ A = πa2; CD = = aρU Re Here the Reynolds number is taken to be Re = 2aU/ν . There is a positive correction for a steady laminar flow, as mentioned earlier. For non-steady wake flow, the drag coefficient increases even further. Hydrodynamics – Spring 2010 120 Boundary layers Even when increasing Re ≫ 1, the flow past an object does not become potential. Why? The reason is that the potential flow has to allow for large tangential velocities on the surface of the object, whereas in real fluids, v = 0 at solid boundaries. Next to the boundary, there should be a boundary layer where the velocity increases from zero to large values somewhat away from the surface. Thus, in the boundary layer, velocity increases rapidly because of the term ν∇2v and viscosity cannot be neglected however small it is. This was first realized by Ludwig Prandtl in 1905. Consider a horizontal flow in the xy plane above a semi-infinite horizontal plate at y = 0, x > 0. Without viscosity, a potential velocity field of v = V ex would be a solution. Thus, outside the boundary layer, this can be taken as the form of the velocity field. Consider the transport of vorticity generated at the tip of the plate. We have y δ V ∂ω = ∇ × (v × ω) + ν∇2ω ∂t ∂ 2ω ∂ω 2 +ν 2 = −v · ∇ω + ν∇ ω ≃ −V ∂x ∂y L Thus, denoting τ = x/V and assuming that vorticity is generated at τ = y = 0 0= 2 ∂ ω ∂ω ≃ν 2 ∂τ ∂y Hydrodynamics – Spring 2010 ⇒ω≃√ ω0 4πντ y2 − 4ντ e ⇒ δ(L) ≃ √ 2ντ = x v u u u t 2νL V 121 Boundary layer equations The N–S equation and the equation of continuity in planar flow read   2 1 ∂p ∂ 2 vx  ∂vx ∂vx  ∂ vx  + vy =− +ν 2 + vx ∂x ∂y ρ ∂x ∂x ∂y 2 vx  2 2  ∂vy 1 ∂p ∂ vy ∂ vy  ∂vy   + vy =− +ν + 2 2 ∂x ∂y ρ ∂y ∂x ∂y ∂vx ∂vy + =0 ∂x ∂y When applying N–S to boundary value problems, we make use of the fact that ∂y ∼ 1/δ , ∂x ∼ 1/x. Equation of continuity then implies that vx/x ∼ vy /δ . It is useful to render the equation into a dimensionless form. We use • a typical linear length scale, l, which makes x′ = x/l of unit order. Thus, δ ′ = δ/l is small. • the asymptotic flow velocity V as the scale of velocities making vx′ = vx/V of unit order and vy′ = vy /V ∼ δ ′ small. • the dynamic pressure ρV 2 as the scale of pressure. Thus, p′ = p/ρV 2 is of order Eu. r Finally, since δ ∼ νx/V , we note that ν ′ = Re−1 = ν/V l ∼ δ ′−2. Hydrodynamics – Spring 2010 122 We can now write down the equations in dimensionless form as   ′ ′ 1  ∂ 2vx′ ∂ 2vx′  ∂p′ ′ ∂vx ′ ∂vx   + =− ′+ vx ′ + vy ∂x ∂y ∂x Re ∂x′2 ∂y ′2 1 · 1 δ′ 1 · 1 Eu δ′ 1 δ′ 1 2 1 1 2 δ′   ′ ′ ∂p′ 1  ∂ 2vy′ ∂ 2vy′  ′ ∂vy ′ ∂vy   + vx ′ + vy ′ = − ′ + ∂x ∂y ∂y Re ∂x′2 ∂y ′2 1 · δ′ δ′ 1 ∂vx′ + ∂x′ · δ′ Eu δ′ ∂vy′ ∂y ′ 1 δ′ 1 δ′ δ′ 2 δ′ δ′ 1 δ′ 2 δ′ =0 Below the equations, the order of each term has been indicated. Keeping only the leading order terms in δ ′ we get from the y -component that ∂y′ p′ = 0 ⇒ p′ = p′(x′). Thus, ′ ′ ∂vx vx ′ ∂x + 1 ∂ 2vx′ dp′ =− ′+ ∂y dx Re ∂y ′2 ′ ′ ∂vx vy ∂vx′ ∂vy′ + =0 ∂x′ ∂y ′ The pressure gradient is determined from the potential part of the flow. Thus, −dp′/dx′ ≃ V ′ dV ′/dx′. Hydrodynamics – Spring 2010 123 One can then apply the so-called boundary-layer transformation, ȳ = y v̄y = ′ vy′ √ √ Re Re to get the boundary layer equations ′ ∂ 2vx′ ∂vx′ ′ dV = −V + v̄y + ∂x ∂ ȳ dx′ ∂ ȳ 2 ′ ′ ∂vx vx ′ ∂vx′ ∂ v̄y =0 + ∂x′ ∂ ȳ with the boundary conditions ȳ = 0 : ȳ → ∞ : vx′ = 0; v̄y = 0 vx′ = V ′(x). For the plate flow example, V ′ = 1. Thus, ∂ 2vx′ ∂vx′ = + v̄y ∂x ∂ ȳ ∂ ȳ 2 ′ ′ ∂vx vx ′ ∂vx′ ∂ v̄y =0 + ∂x′ ∂ ȳ ȳ = 0, x′ > 0 : with ȳ → ∞ : Hydrodynamics – Spring 2010 vx′ = 0; vx′ = 1 v̄y = 0 124 Solution of the plate flow BL problem Boundary layer problem, being an incompressible 2D N–S problem, can be most easily solved by introducing a streamfunction: ∂ψ ∂ψ vx′ = ; v̄y = − ′ . ∂ ȳ ∂x r r √ We deduced that δ ≃ 2νx/V = 2x′/Re, i.e., δ̄ ≃ 2x′. Thus, when solving the boundary layer equations, use of a similarity variable √ η = ȳ/ 2x′ will be made. It turns out that a streamfunction of form ψ= is useful. Thus, √ √ 2x′f (η) df ∂η df = ≡ f˙ ∂ ȳ dη dη  √ ∂η df f 1  ˙ ′ √ √ v̄y = − = ηf − f − 2x ′ ∂x dη 2x′ 2x′ vx′ = 2x′ with ... f + f f¨ = 0 with Hydrodynamics – Spring 2010 η=0: η→∞: f˙ = 0; f˙ = 1 f =0 Blasius equation 125 http://faculty.virginia.edu/ribando/modules/xls/ √ √ ′ Note: some sources (like the one above) use the similarity transformation η = ȳ/ x and ψ = x′f (η), which yields the Blasius equation in form ... 1 f + 2 f f¨ = 0; Hydrodynamics – Spring 2010 f (0) = f˙(0) = 0; f˙(∞) = 1 126 Aerodynamic lift We will, next, try to understand the aerodynamic lift, which allows birds, airplanes and helicopters to fly and contributes to the down-force of the race cars. For simplicity, we shall limit the discussion to airfoils. One can think of an airfoil as the cross-section of a translationally symmetric wing that has a much longer span than chord and a leading edge perpendicular to the flow direction. The flow around the wing is, thus, two-dimensional. The lift of an airfoil is based on the principle of Bernoulli: if the upper surface of the wing is longer than the lower surface and/or if the angle of attack is positive, the flow above the wing is faster than below. Thus, there is underpressure above the wing that generates the lift. We shall see that the lift of an airfoil is proportional to the circulation Γ around the airfoil, just as for the ideal flow around a cylinder considered in the exercises. But what produces the circulation around the airfoil? The answer is viscosity. The flow pattern around the airfoil without any circulation contains extremely large velocity gradients near the pointed end of the airfoil, which are relaxed by the generation of vorticity, i.e., circulation around the airfoil. Hydrodynamics – Spring 2010 viscous forces very large 127 Calculation of the lift on an airfoil While viscosity is responsible for the initial generatation of circulation around the airfoil, we expect its effect to be small outside the boundary layer around √ the airfoil. The thickness of a (laminar) boundary layer near the trailing edge of the airfoil is δ ∼ 5L/ Re. The kinematic viscosity of air near the ground is ν ∼ 10−5 m2s−1. If the speed of a model aircraft is V ∼ 10 ms−1 and the chord of the wing is L ∼ 0.1 m, we get Re ∼ 105 (still laminar). Thus, δ/L ∼ 0.2% and the velocity field around the airfoil is well represented as a potential flow outside the thin boundary layer attached to the airfoil. Consider the two-dimensional velocity field in the complex xy plane, i.e., take w(z) = vx − i vy , z = x + i y. Because the flow is potential, w is an analytic function. Thus, it can be expanded as a Laurent series around the origin (taken to be inside the airfoil) as ˛ 1 a−1 a−2 + 2 + ...; a−k = z k−1w(z) dz, w(z) = v∞ + z z 2π i C where the integral can be taken along any (counter-clockwise) closed curve C enclosing the origin, and use has been made of w → v∞ at |z| → ∞. Thus, ˛ ˛ ˛ ˛ 2πi a−1 = w dz = (vx − i vy ) (dx + i dy) = (vxdx + vy dy) − i (vy dx − vxdy). C Hydrodynamics – Spring 2010 C C C 128 Choosing now the curve to follow the streamline just above the thin boundary layer, we find that the second integral vanishes, since dx/vx = dy/vy along a streamline.5 Thus, 2πi a−1 = ˛ (vxdx + vy dy) = C ˛ C The net force on the airfoil per unit span can be obtained by integrating the force exerted by pressure, ˛ F = pn dl, where n is the inward normal of the surface. Taking θ = (dl, ex), we have n = − sin θ ex + cos θ ey . v · dl = Γ. y ey C dl x ex n Since p is not a function of the normal distance through the boundary layer, we can again take the integral along the streamline enclosing the airfoil and the boundary layer. Bernoulli’s principle states that p + 21 ρv 2 = constant along a streamline (or along any line in irrotational flow), so ˛ F = − 12 ρ v 2 n dl. C 5 There are two sections of C of length ∼ δ above the stagnation points, where the integration cannot be taken along a streamline, but as δ/L ≈ 0, these can be neglected. Hydrodynamics – Spring 2010 129 Introduce a complex force Z = Fy + i Fx = − 12 ρ = − 21 ρ ˛ = − 21 ρ ˛ C C ˛ C v 2 (ny + i nx) dl = − 12 ρ v 2 (cos θ dl − i sin θ dl) = − 21 ρ ww∗ dz ∗ = − 21 ρ ˛ C ˛ C ˛ C v 2 (cos θ − i sin θ) dl v 2 (dx − i dy) w(w dz + w∗ dz ∗ − w dz) and since w∗ dz ∗ − w dz = −2i ℑ(w dz) = −2i ℑ{(vx − i vy )(dx + i dy)} = −2i(vxdy − vy dx) vanishes along any streamline, we have ˛ Z = − 12 ρ Now, using the Laurent expansion of we get w2 dz. C v∞ Γ 2v∞a−1 2 2 + ... = v∞ + ... + w2 = (v∞ + a−1/z + ...)2 = v∞ + z πi z ˛ ˛ v∞ Γ 1 ⇒ = w2 dz ⇒ w2 dz = 2v∞Γ πi 2π i ˛ Z = − 21 ρ w2 dz = −ρv∞Γ, C i.e., Fx = 0 and Fy = −ρΓv∞. Hydrodynamics – Spring 2010 Kutta–Joukowski 130 Joukowski transform. Kutta condition Lift is Fy = −ρv∞Γ regardless of the airfoil shape. But what is the value of Γ? This can be deduced by conformal mapping of a cylinder to an airfoil shape. For example, by  Joukowski transform: z = e  1 ζ + . ζ −iα  r Cylinder centered at µ, radius a = (1 − µx)2 + µ2y ζ = µ + ae−iβ ; β ∈ [0, 2π] Thus, ζ = 1 when sin β = µy /a. This point maps to the cusp. The stagnation point is at sin(α+β) = −Γ/4πv∞a (exercise). Mapping it to the cusp gives µ = −0.12 + 0.08 i; α = 10◦ Hydrodynamics – Spring 2010 −1 Γ = −4πv∞a sin α + sin µy ! a Kutta condition 131 Accretion of mass Accretion of mass on compact massive objects is one of the most efficient energy release mechanisms. For example, the loss of potential energy of a particle of mass m falling through the gravitational field of a neutron star (from infinitely far to the surface) is GM⋆m GM⋆ = 2 mc2, a ca where the neutron star mass M⋆ ≃ M⊙ and the neutron star radius a ≃ 10 km. The factor GM⋆/c2a ≃ 0.15, so the gravitational energy available in such accreting material is more than what can be obtained from a fusion process converting hydrogen to heavier elements! If accreting matter around a massive compact object possesses angular momentum, it forms an accretion disk around the object. Accretion disks are important ingedients in many astrohysical objects such as X-ray binaries and active galactic nuclei. In general, accretion is a complicated phenomenon to model. We shall limit our discussion to only two forms of accretion, thin accretion disks (to be studied now) and spherical accretion (to be studied later together with stellar winds). Furthermore, we shall also neglect self-gravity in all discussion, thus assuming that the disk is much less massive than the accreting object. Hydrodynamics – Spring 2010 132 Accretion disks Consider matter orbiting an object of mass M in a nearly circular orbit of radius r. Thus, balancing the gravitational acceleration −GM er /r2 with the centrifugal acceleration vθ2er /r gives for the angular velocity Ω = vθ /r 1/2  GM Ω(r) =  3  r which is often referred to as Keplerian motion (for obvious reasons). As the angular velocity is a function of radius, the rotation of matter is not rigid. Thus, a velocity shear is produced. It is the viscous forces related to this shear that allows the dissipation of mechanical energy to occur. In the absence of viscosity, particles on circular orbits conserve their total mechanical energy ǫ = T + V = 21 r2Ω2 − GM GM =− = 12 V r 2r Particles losing mechanical energy cannot stay on fixed circular orbits. A slow loss of mechanical energy means that ṙ = −rǫ̇/ǫ < 0 as ǫ̇ < 0, ǫ < 0 and, thus, particles are slowly spiraling inwards while orbiting the massive object. Thus, a particle starting at infinity (ǫ = 0) and spiralling to the surface of the object (ǫ = − 12 GM/a) has lost half of the potential energy difference on the way in. This energy goes to heat and, ultimately, to radiation. As ǫ̇ = Fθ vθ and Fθ ∝ ν , we expect luminosity L = 21 GM ṁ/a ∝ ν , where ṁ ∝ ν is the accretion rate on the object. Thus, although there were no fixed solid bodies in the accretion disk and although the Reynolds number of the system is enormous, viscosity is still in a central role in the physics of accretion. Hydrodynamics – Spring 2010 133 The basic disk dynamics Consider a thin accretion disk in cylindrical coordinates with Ω k ez . Assume azimuthal symmetry and vz = 0. Thus, equation of continuity and the azimuthal component of the N–S equation are ∂ρ 1 ∂ + (rvr ρ) = 0 ∂t r ∂r   v v ∂v ∂v θ r θ θ  = {∇ · [2µ(Λ − 1 I tr Λ)]} , ρ + vr + θ 3 ∂t ∂r r where we have assumed that bulk viscocity vanishes. Note that because we are now considering the effects of viscocity at large scales, we cannot assume that µ is constant. Rather than writing down the viscous terms in terms of the components of Λ, we´ shall derive the term from first principles. First, however, integrate the equation over z writing Σ = dz ρ: ∂Σ 1 ∂ + (rvr Σ) = 0 ∂t r ∂r ˆ   ∂v v v ∂v θ r θ θ  = Σ dz {∇ · [2µ(Λ − 31 I tr Λ)]}θ , + vr + ∂t ∂r r where we have assumed that vθ and vr do not depend on z . Multiply the upper equation by rvθ and add to the lower multiplied by r to get 1∂ ∂ (Σr2Ω) + (rvr Σr2Ω) = G; ∂t r ∂r Hydrodynamics – Spring 2010 Ω= vθ r 134 Note that Σr2Ω is the surface density of angular momentum and vr Σr2Ω is its flux in the radial direction. As it does not depend on any other coordinate than r, the second term is clearly the divergence of the angular momentum flux. The equation, therefore, states the conservation of angular momentum, and G must be related to the viscous torque. Multiply the angular momentum transport equation by 2πr. Thus, ∂ ∂ (2πr Σr2Ω) + (vr 2πr Σr2Ω) = 2πr G ∂t ∂r As dr L(r, t) = 2πr dr Σr2Ω is the angular momentum carried by an annulus of width dr of the disk at r, we have   ∂L ∂(v L) r    = G(r + dr) − G(r), dr  + ∂t ∂r where G(r) is the viscous torque at radius r. It is the torque exerted on the material inside r by the material outside r through viscous forces. Thus, 2πr G = ∂G/∂r and so we need to find G. The viscous stress is τθr = µr ∂Ω/∂r so the viscous torque per unit area of a cylindrical surface is r τθr . The viscous torque is, thus, ˆ ˆ ∂Ω ∂Ω = 2πνΣr3 G(r, t) = r dθ dz µr2 ∂r ∂r and Hydrodynamics – Spring 2010   ∂ 1∂ 1∂  ∂Ω νΣr3  (Σr2Ω) + (rvr Σr2Ω) = ∂t r ∂r r ∂r ∂r 135 The basic equations governing the disk dynamics are, thus, ∂Σ 1 ∂ + (rvr Σ) = 0 ∂t r ∂r   ∂ 1∂ 1∂  3 ∂Ω  2 2 νΣr . (Σr Ω) + (rvr Σr Ω) = ∂t r ∂r r ∂r ∂r In order to obtain a dynamical theory based on these two equations, one of the unknowns Σ, vr and Ω needs to be prescibed. Considering the angular velocity to be nearly Keplerian, we can regard Σ and vr as the dynamical variables of the theory and solve the equations. As Ω ∝ r−3/2, we can write the second equation (LHS first) as  ∂Σ vr Σ 1 ∂ 3/2 ∂Σ 1 ∂ 1 ∂  3 1/2 = νΣr − + (r v Σ) = + (rv Σ) + r r 2 r3/2 ∂r ∂t r3/2 ∂r r ∂r 2r |∂t {z } =0  ∂  1/2 ⇒ rvr Σ = −3r νΣr ∂r    ∂Σ 3 ∂  1/2 ∂  r = νΣr1/2  = ∇ · (D∇Σ − V Σ), ∴ ∂t r ∂r ∂r 1/2 where D = 3ν and V = −3r−1/2∂r (r1/2ν) er . The diffusion term tries to spread the disk whereas convection term drives the accretion on the central mass. Note that the mass influx, −2πrvr Σ is positive as long as r1/2νΣ increases with with distance. The disk equation can be solved analytically, e.g., for a constant ν and Σ(r, 0) = (m/2πr0) δ(r − r0) (see Choudhuri). Hydrodynamics – Spring 2010 136 Near-Keplerian orbits Consider the two remaining components of the N–S equation: 1 ∂p GM ∂vr vθ2 − =− − 2 vr ∂r r ρ ∂r r 0=− 1 ∂p GM z − , ρ ∂z r3 where the radial √ component of the gravitational force components are approximated assuming z ≪ r, i.e., taking r2 + z 2 ≈ r, and where viscosity has been neglected as the other terms are much larger. The thin disk approximation immediately yields another simplification. Since the thickness of the is disk h, then |∂p/∂z| ∼ p/h so that from the z component of N–S we get GM h p ∼ ρh r3 GM h2 GM p 1 ∂p ⇒ ∼ 2 2 ≪ 2 ∼ ρ ∂r ρr r r r so that the pressure term can be neglected from the radial component of N–S. Once we again realize that vr ≪ vθ ,we see that the radial component yields the balance between centrifugal and gravitatinal force that led to the Keplerian angular velocities. To summarize, near-Keplerian motion is valid for thin disks, for which pressure forces can be neglected. If pressure terms are important, the disk is no longer thin and the analysis complicates considerably. Hydrodynamics – Spring 2010 137 Thin disk in steady state In steady state 1d (rvr Σ) = 0 r dr   1d 1 d ∂Ω νΣr 3  (rvr Σr2Ω) = r dr r dr ∂r Both equations integrate immediately to give rvr Σ = C1 Σr3Ωvr − νΣr3 dΩ = C2 dr where C1 and C2 are constants of integration. For a steady disk, the mass inflow rate ṁ = −2πr vr Σ. Thus, ṁ C1 = − . 2π The second integration constant is obtained at the inner radius of the disk, r = r⋆, where we impose the condition of rigid rotation, dΩ/dr = 0 and Ω = Ω⋆. Thus, C2 = C1r⋆2Ω⋆ = − Hydrodynamics – Spring 2010 ṁ 2 r Ω⋆. 2π ⋆ 138 In summary,  ṁ 2 ṁ  2 2 νΣ = Σr Ωvr + r⋆ Ω⋆ = r Ω⋆ − r Ω r ∂r 2π 2π ⋆ ṁ r⋆2Ω⋆ − r2Ω , νΣ = 2π 3 ∂Ω r ∂r 3 ∂Ω 3 i.e., ṁ ∝ ν , as expected. The viscous dissipation per unit volume is µr2(dΩ/dr)2. This gives the energy emitted from the unit area of the disk as ˆ   2 2 2 2 ṁ Ω − r Ω ∂Ω dΩ dΩ r dE ⋆ ⋆ 2 2  dz = νΣr   = = µr − dt dr dr 2π r ∂r Integrating this over the disk area gives the luminosity ∞ ˆ ∞ ṁ r⋆2Ω⋆ − r2Ω ∂Ω ∂Ω L= (r⋆2Ω⋆ − r2Ω) 2πr dr = ṁ dr 2π r ∂r ∂r r⋆ r⋆ ˆ ∞ ∂Ω r2Ω = ṁr⋆2Ω⋆(−Ω⋆) − ṁ dr ∂r r⋆   ˆ ∞ ˆ ∞ 1 2 2  − r Ω + rΩ2 dr rΩ2 dr = ṁ  = −ṁr⋆2Ω2⋆ + 12 ṁr⋆2Ω2⋆ + ṁ 2 ⋆ ⋆ ˆ r⋆ r⋆ which equals 12 GM ṁ/r⋆, as expected, if the Keplerian formula for Ω is used. Hydrodynamics – Spring 2010 139 The Keplerian formula does not apply in the region close to the inner radius of the disk, where the matter is rotating rigidly. Thus, one should write  1 2 2 − r Ω + L = ṁ  2 ⋆ ⋆ ∞ ˆ = − 21 ṁr⋆2Ω2⋆ + ṁ r r⋆ ˆ vθ2 ∞ r⋆   dr  ∂vr 1 ∂p GM  v + + 2 dr r ∂r ρ ∂r r     GM 1 2 GM  2 = −ṁ  12 r⋆2Ω2⋆ + 12 vr,⋆ = −ṁ − + 2 v⋆ + P⋆ + P⋆ − r⋆ r⋆ where ˆ P(r) = − r ∞ 1 ∂p dr = ρ ∂r ˆ p(r) p∞ dp = ∆w − ∆q, ρ(p) ∆w = w(r) − w∞ is the change in specific enthalpy, w = ǫ + p/ρ, and ∆q the deposited heat per unit mass. Here, we have used the first law of thermodynamics in form dw = dq + dp . ρ As we assumed a detailed balance between radiative losses and viscous heating, we should take ∆q = 0, so we get the luminosity as L = −ṁ ∆( 12 v 2 + w + Φ), where Φ = −GM/r is the gravitational potential. Thus, the gravitational potential energy is converted to kinetic energy, enthalpy and heat, which in turn is converted to radiation. Hydrodynamics – Spring 2010 140 Gas dynamics Hydrodynamics – Spring 2010 141 Compressible fluids Problems so far studied neglected the effects of compressibility of the fluid. Next, we will turn to the dynamics of compressible fluids. When comparing the Euler and N–S equations, we found that the second derivative in the N-S equation, however small in most regions of space, changed the nature of the differential equation completely allowing solutions with the no-slip boundary conditions to be obtained. In a similar manner, compressibility affects the fundamental nature of the solutions of fluid dynamical equations. Introducing compressibility gives a finite speed of propagation for mechanical disturbances in the fluid, i.e., the speed of sound cs. At low speeds (v ≪ cs) the effect of finite speed of information is negligible, and incompressibility is a good approximation. As long as v < cs, the fundamental nature of the solution does not change – the HD equations will be of the elliptic type. At high flow velocities (v > cs), however, information cannot propagate against the flow. This changes the HD equations from elliptic to hyperbolic type. The solutions are then fundamentally different from the incompressible solutions, as we shall see. Compressibility effects are most pronounced for problems involving dilute gases rather than liquids (with the important exception of sound propagation in liquids). Thus, we will consider the fluid as a perfect gas. Hydrodynamics – Spring 2010 142 Thermodynamics of a perfect gas The perfect gas law can be stated as p = RρT, where R = κB/m is the gas constant. The internal energy per unit mass is ǫ = cV T, where cV = 12 f R is the specific heat per unit mass. Here, f is the number of degrees of freedom (DoF) of the molecule, i.e., f = 3 (3 translational DoF) for monoatomic gas and f = 5 for a diatomic gas (3 translational and 2 rotational DoF). It is typical, instead of f , to use γ = (f + 2)/f = cp/cV , i.e., cV = R γ−1 ⇒ ǫ= RT . γ−1 Consider the the usual thermodynamic definition of entropy as δQ = T dS . Introduce then the extensive quantities as given per unit mass, δQ = ρ δq and dS = ρ ds, and write down the first law of thermodynamics in form     dT dρ dρ 1 − (γ − 1)  T ds = dǫ + p d   = dǫ − p 2 = cV T  ρ ρ T ρ     dρ dρ dp p dT dρ + − γ  = cV T  − γ  = cV T d ln γ = cV T  T ρ ρ p ρ ρ p ⇒ s = s0 + cV ln γ ρ Hydrodynamics – Spring 2010 143 Adiabatic flows For an adiabatic process, 0 = δQ = ρT ds giving ds =0 dt ⇒   d p =0 dt ργ which replaces the energy equation. The enthalpy per unit mass is defined through p γ w =ǫ+ = RT ρ γ−1   1 dp dp ⇒ dw = dǫ + p d   + = T ds + ρ ρ ρ so for adiabatic flows (ds = 0) dp dw = ρ ⇒ ˆ dp γ =w= RT. ρ γ−1 Thus, the principle of Bernoulli for adiabatic flows states that 1 2 2v + γ RT + Φ = constant γ−1 along a streamline, Φ being the potential of the conservative external force field, F = −∇Φ. Hydrodynamics – Spring 2010 144 Acoustic waves In order to concentrate on the effects of compressibility, we shall neglect transport phenomena (i.e., set µ = 0 = K ) for a while. Consider a homogeneous perfect gas at rest, and denote the quantities in this equilibrium (in the absence of external forces) by ρ = ρ0; p = p0; v = v 0 = 0. Suppose that pressure is perturbed to p = p0 + p1(x, t) so that the corresponding density and velocity are ρ = ρ0 + ρ1(x, t) and v = v 1(x, t). Assuming that the perturbations are adiabatic, we can relate   p dρ1 1 dp1 d p − γ γ+1 0 =  γ = γ dt ρ ρ dt ρ dt ⇒ dp1 = c2s dρ1 where c2s = γp/ρ. Assuming that the perturbations are small, ρ1 ≪ ρ0, p1 ≪ p0, we can replace the quantities in c2s by their equilibrium values, i.e., write γp0 c2s = , ρ0 which linearizes the equation and gives p1 = c2s ρ1. Using the same assumption in the equation of continuity, we get ∂ρ1 ∂ρ1 + ∇ · [(ρ0 + ρ1)v 1] ≈ + ρ0 ∇ · v 1 0= ∂t ∂t providing an equation between ρ1 and v 1. Hydrodynamics – Spring 2010 145 Next, apply the same technique of linearization to the Euler equation   ∂v 1 ∂v 1 ∂v 1 + (v 1 · ∇)v 1 + ∇p1 ≈ ρ0 + ∇p1 = ρ0 + c2s ∇ρ1. 0 = (ρ0 + ρ1)  ∂t ∂t ∂t Combining this with the linearised equation of contiunuity gives ∂ 2 ρ1 − c2s ∇2ρ1 = 0 2 ∂t which is the equation for acoustic waves showing that cs = v u u u t γp0 q = γRT ρ0 is the sound speed. Newton (1689) first derived the expression for sound speed assuming that the perturbations were isothermal, corresponding to γ = 1 in our equation. This expression gives cs ≈ 280 m/s for a gas at 0◦C, much lower than the experimental value. It was Laplace (1816) who realised that the perturbations should be regarded adiabatic and taking γ = 1.4, appropriate for air, gives cs ≈ 332 m/s, which agrees well with the experimental value. We note that assuming that p1 = (dp/dρ)0ρ1 gives c2s = (dp/dρ)0 which applies for all fluids. Liquids, which are hard to compress, yield a small ρ1 with a large p1 meaning that the sound speed in liquids is large. (For water, cs ≈ 1430 m/s.) Hydrodynamics – Spring 2010 146 Harmonic analysis Taking ρ1 = ρ̄1 exp{i(k · x − ωt)} and substituting this into the equation for acoustic waves gives (−ω 2 + c2s k 2)ρ̄1 exp{i(k · x − ωt)} = 0 ⇒ ω 2 = c2s k 2 as the dispersion relation for acoustic waves. The phase and group velocities of the acoustic waves are given by vp ≡ ωk k ; = c s k2 k vg ≡ ∂ω k = cs ∂k k The phase velocity is always perpendicular to and gives the velocity of the wave fronts (i.e., planes with constant values of the phase, φ = k · x − ωt). The group velocity gives the velocity of individual wave packets (and, thus, the propagation velocity of information and energy). Acoustic waves in a homogeneous medium are non-dispersive, i.e., have phase and group speeds, vp = cs = vg , independent of ω . Waves in more complicated media, like gravitationally stratified atmospheres, have more complicated dispersion relations and turn out to be dispersive. It follows from the Euler equation, −iωρ0v 1 + c2s ikρ1 = 0 that v 1 k k. Thus, acoustic waves are longitudinal waves. Hydrodynamics – Spring 2010 147 Non-linearity. Steepening into shock waves When deriving the dispersion relation of acoustic waves, we assumed that non-linear terms in the HD equations are vanishingly small. What happens if we relax this assumption? For simplicity, let us consider a 1-D situation, where a purely longitudinal wave propagates in the xdirection. Thus, ∂ ∂ρ + (ρv) = 0 ∂t ∂x ∂v ∂v 1 ∂p c2s ∂ρ +v =− =− ∂t ∂x ρ ∂x ρ ∂x where a polytropic relation between the pressure and the density is assumed and where c2s = dp/dρ. Clearly, there are several non-linear terms in this system of equations. However, assuming that the Mach number of the flow, v M≡ cs is large (cs ≈ 0), we obtain ∂v ∂v +v = 0, ∂t ∂x which already contains non-linearity but is much simpler to solve than the full set of equations. Let us, next, consider the solution of the simplified Euler equation, containing only the inertial force, in the xt-plane. Hydrodynamics – Spring 2010 148 Thus, ∂v ∂v +v = 0. ∂t ∂x Consider the values of v along the characteristic curves given by dx = v(x, t). dt Since dv ∂v ∂v dx ∂v ∂v = + = +v = 0, dt ∂t ∂x dt ∂t ∂x the value of the velocity v does not change along the characteristics, so the curves are actually straight lines in the xt-plane. The exact solution can be formally written as v = g(x − vt), where g(x) = v(x, 0) is the initial velocity profile. E.g.,          v0 , g(x) =  −v0,        x < −L x>L −v0x/L, |x| < L          v0 , ⇒ v(x, t) =  −v0,        x < v0t − L x > L − v0 t −v0x/(L − v0t), |x| < L − v0t Hydrodynamics – Spring 2010 149 At t = L/v0 the solution          v0 , x < v0t − L v(x, t) =  −v0,  x > L − v0 t       −v0x/(L − v0t), |x| < L − v0t has reached a state where it becomes discontinuous, i.e., a step from v = −v0 to v = v0 at x = 0. Clearly, the solution only applies for t < L/v0 since for larger values of time the function becomes multivalued. At times t < L/v0 the velocity profile is becoming steeper and steeper around x = 0. This phenomenon, caused by the term v ∂v/∂x and called non-linear steepening, leads to the formation of a shock wave for any large-amplitude acoustic wave packet. For a sinusoidal wave, i.e., g = −v0 sin kx, the region near x = 0 has g ≈ −v0kx. Thus, for |x| ≪ k −1 v ≈ −v0k(x − vt) ⇒v≈− v0kx 1 − v0kt so the time needed for shock formation through non-linear steepening can be estimated as t ≈ (kv0)−1, where k is the wavenumber and v0 is the amplitude of the wave. The physics of the shocks, once they have formed, can be analysed using the conservation laws for mass, momentum and energy. Hydrodynamics – Spring 2010 150 Structure of shock waves A shock wave is a region of small thickness over which the dynamical variables of the fluid change rapidly. Shocks are always related to compression, as steepening to shocks will only occur in regions with ∂v/∂x < 0. (A region of fluid with rarefaction, ∂v/∂x > 0, will get wider as a result of the inertial term.) At the limit of very thin transition, shocks can be treated as discontinuities in density, velocity and pressure. As we treat shocks as discontinuities, we can regard them as locally planar fronts with unit normal n. Observing the shock from the rest frame of the ambient fluid, we can take the velocity of the shock front as v s = vsn. Transforming to the coordinate system co-moving with the shock at this velocity, we can treat the flow near the shock as being independent of time. In this frame of reference, the gas is flowing in from the ambient medium to the shock at velocity −vsn. It is customary to choose one of the coordinate axes, say x, parallel to the shock normal so that the velocity of the gas through the shock is positive, i.e., ex = −n. Thus, vx = vs at x < 0. This region is called the upstream region of the shock. In the downstream region (i.e., x > 0) the velocity is also in the x direction because of symmetry. Hydrodynamics – Spring 2010 ex vs n vx = v s x 151 Notation and jump conditions It is customary to denote the values of the hydrodynamic quantities (ρ, v and p) measured in the upstream region by subscript p1, ρ 1 p2, ρ 2 1 and those measured in the downstream region by subscript 2. The ideal hydrodynamic equations in conservation form can v1 v2 be written as ∂ρ ∂ + (vρ) = 0 ∂t ∂x x ∂ ∂ (ρv) + (ρv 2 + p) = 0 ∂t ∂x ∂ 1 2 ∂ ( 2 ρv + ρǫ) + [( 12 v 2 + w)ρv] = 0. ∂t ∂x In the shock frame the system is in a steady state and the conservation laws can be integrated accross the discontiunuity, using wρ = γp/(γ − 1), to give v 1 ρ 1 = v 2 ρ2 ρ1v12 + p1 = ρ2v22 + p2 γ p1 1 2 γ p2 1 2 = . v v + + 2 1 γ − 1 ρ 1 2 2 γ − 1 ρ2 Regarding the upstream quantities as known, the downstream quantities can be solved from these jump conditions (or Rankine–Hugoniot conditions). Hydrodynamics – Spring 2010 152 Solution of the jump conditions The jump conditions can be solved to give the ratios v2/v1, ρ2/ρ1 and p2/p1 as a function of the Mach number v1 v1 =r M= cs1 γp1/ρ1 of the shock. The solution is (exercise) v 1 ρ2 (γ + 1)M2 = = v2 ρ1 2 + (γ − 1)M2 p2 2γM2 − (γ − 1) = p1 γ+1 Clearly, the upstream and downstream quantites agree if M = 1. For values M < 1, the compression ratio ρ2/ρ1 < 1, and because shocks are always compressive, the equation does not apply. Furthermore, the compression ratio of the shock increases monotonically as a function of the Mach number and approaches a finite limit at infinity: γ+1 ρ2 , → ρ1 γ−1 as M → ∞. For γ = 53 , γ+1 = 4. γ−1 The pressure ratio, on the other hand, increases without a limit as M → ∞. It is straightforward to show that the specific entropy s = κB ln(p/ργ ) increases over the shock as long as ρ2 > ρ1, consistent with the second law of thermodynamics. Hydrodynamics – Spring 2010 153 More shock properties The downstream temperature is p 2 ρ1 [2γM2 − (γ − 1)][(γ − 1)M2 + 2] p2/ρ2 = T1 = T1 T2 = T1 p1/ρ1 p 1 ρ2 (γ + 1)2M2 which also increases as a function of Mach number. Thus, shock waves heat the gas. For strong shocks ǫ2 ≈ 4 RT1 2γ(γ − 1) 2 1 2 M = v γ − 1 (γ + 1)2 (γ + 1)2 2 1 showing that strong shocks efficiently convert upstream kinetic energy to downstream thermal energy. The Mach number of the downstream flow is M22 ρ2v22 ρ1v12 ρ1 p1 (γ − 1)M2 + 2 v22 2 ρ1 p1 = =M = = 2 = cs2 γp2 γp1 ρ2 p2 ρ2 p2 2γM2 − (γ − 1) which yields M2 < 1 for all M > 1. Thus, shock waves are transitions from supersonic to subsonic flow. This property is important, as it shows that shock waves forming ahead of obstacles in a supersonic flow provide the means for information about the obstacle to propagate to the gas. The thickness of a shock is, of course, finite in reality. The steepening of an acoustic wave is limited by the transport processes (viscosity and heat conduction), which gives the shock thickness roughly as L ∼ ν2/∆v . Hydrodynamics – Spring 2010 154 Blast waves Consider a sudden release of energy within of a localized region of gas, e.g., a supernova explosion or a nuclear detonation in the atmosphere. We expect the explosion to generate a (spherical) shock wave traveling through the gas. The task is to find out how the radius of the shock evolves with time and how the hydrodynamic quantities behave inside the spherical shock. Suppose that the amount of energy released is E and that the density of the ambient medium is ρ1. Suppose further, that the temperature of the ambient medium is cold enough so that the pressure of the ambient medium is negligible compared to the pressure inside the blast wave. Thus, the ambient temperature can be assumed to be insignificant for the evolution of the blast wave. Let λ be a length scale giving determining radius of the blast wave at a given time t. Let us use dimensional analysis to determine λ. Thus, assume that λ = E a ρb1 tc. Dimensionally L = (M L2T −2)a (M L−3)b T c = M a+b L2a−3b T c−2a          a+b=0 ⇒ 2a − 3b = 1        c − 2a = 0          a= 1 5 ⇒ b = − 51        c= 2 5  1/5 2  Et   ∴λ= ρ1 Thus, the radius of the blast wave rS(t) ∼ (Et2/ρ1)1/5 within a factor of unit order. The expansion of a blast wave has been confirmed to follow this law by observations of nuclear explosions. Thus, it gives a valid model of the expansion of supernova remnants. Hydrodynamics – Spring 2010 155 Similarity variable. Fluid variables at the shock Now, consider the expansion of the gas behind the spherical blast wave. Let us consider the solution of the hydrodynamic equations under the assumption of self-similarity, i.e., that the hydrodynamic quantities behind the shock depend on a dimensionless similarity variable r ρ1 !1/5 ξ= =r . λ(t) Et2 Since the radius of the shock scales like λ, we assume that the position of the shock in terms of the similarity variable is ξ0, i.e., that   2 1/5 Et    . rS(t) = ξ0λ = ξ0  ρ1 Thus,   drS 2rS 2ξ0  E 1/5 vS(t) = = = . dt 5t 5 ρ1t3 The shock can be treated at the limit M → ∞. Thus, the values of the hydrodynamic variables immediately behind the shock front are vS ρ2 γ + 1 = = ρ1 γ − 1 v S − v 2 ⇒ v2 = 2 vS γ+1 2γM2 2 p2 = p1 = ρ1vS2 . γ+1 γ+1 Hydrodynamics – Spring 2010 156 Similarity solution Let us, then, try to find a solution behind the shock in form ρ(r, t) = ρ2ρ′(ξ) = ρ1 γ+1 ′ ρ (ξ) γ−1 ξ 4 r ′ r v(r, t) = v2 v ′(ξ) = v2 v ′(ξ) = v (ξ) ξ0 rS 5(γ + 1) t 2 !2 2 ξ 8ρ r r 1 ′ ′ p(r, t) = p2   p (ξ) = p2   p (ξ) = p′(ξ), ξ0 rS 25(γ + 1) t     where ρ′, v ′ and p′ are dimensionless functions of the similarity variable. Subsitute, next, these into to the ideal HD equations, ∂ρ 1 ∂ 2 + (r ρv) = 0 ∂t r2 ∂r ∂v ∂v 1 ∂p +v =− ∂t ∂r ρ ∂r   p ∂ ∂  + v  ln γ = 0 ∂t ∂r ρ Hydrodynamics – Spring 2010 157 d 2  ′ ′ dρ′ 3ρ v + ξ (ρ′v ′) = 0 + −ξ dξ γ + 1 dξ  ∴ ′   ′   ′  4 2(γ − 1) 1  ′ dp  2 dv  ′2 ′ dv  v + v ξ  = −  2p + ξ + −v ′ − ξ ′ 5 dξ 5(γ + 1) dξ 5(γ + 1) ρ dξ   d  p′  5(γ + 1) − 4v ′ ξ ln ′γ  = dξ ρ 2v ′ − (γ + 1) Thus, as r and t are eliminated from the system of equations, our trial form of the dependent variables works. The task is now to solve this system of first-order ordinary differential equations with the boundary conditions at the shock ρ′(ξ0) = v ′(ξ0) = p′(ξ0) = 1. The value of ξ0 can be found by requiring that the total energy of the blast wave remains constant: ˆ r   S 1 p  4πr 2 dr  ρv 2 + E= 2 γ−1 0 ˆ ξ0 32π ′ ′ ′2 4 (p + ρ v ) ξ dξ. ⇒ 1= 25(γ 2 − 1) 0 Solving the equations for the primed quantities gives the similarity solution for the blast wave that can be scaled to any values of E and ρ1. Obviously, the solution is not valid (1) at the earliest times, as it predicts v2, vS → ∞ with t → 0, and (2) at the latest times, when radiative cooling renders the assumption of constant energy invalid. Hydrodynamics – Spring 2010 158 Plot of the similarity solution Hydrodynamics – Spring 2010 159 One-dimensional gas flow Let us consider, next, steady-state solutions involving supersonic flows. Consider a one-dimensional, steady, adiabatic gas flow through a pipe with variable cross section. The governing equations are: 1 d (Avρ) = 0 A dx   dp d p 2 dρ = c = 0 ⇒ s dx ργ dx dx with c2s = γp ρ dv 1 dp c2s dρ v =− =− dx ρ dx ρ dx Thus, c2s dρ 1 dA 1 dv  = c2s  + − ρ dx A dx v dx  ⇒ (1 − M2) 1 dv 1 dA =− , v dx A dx  where M = v is the local Mach number. cs If the flow is subsonic, then dv/dx and dA/dx have opposite sign and decreasing cross-sectional area leads to an accelerating flow (as intuitively expected). However, if the flow is supersonic, increasing area yields (counterintuitively) an accelerating flow. Hydrodynamics – Spring 2010 160 de Laval nozzle We found that (1 − M2) 1 dA 1 dv =− . v dx A dx Consider a case, where the flow starts as subsonic and ends up supersonic. This requires that the sonic point M = 1 is reached at dA/dx = 0 and that the pipe has to first decrease and then increase in cross section. This kind of a pipe is called de Laval nozzle and it is used in rocket engines. The exit velocity of the gas is ve = v u u u u t  (γ−1)/γ  e     p 2γ RTi  1 − γ−1 pi Hydrodynamics – Spring 2010 . 161 Extragalactic jets Radio galaxies often appear with huge jets. These can be modeled as being channeled through de Laval nozzle -like geometry. (Blandford & Rees 1974) Hydrodynamics – Spring 2010 162 Hydrodynamics – Spring 2010 163 Spherical accretion and winds Let us, finally, treat the problem of mass flow onto or from a massive object in spherical symmetry, i.e., spherical accretion or spherical wind. These problems were first treated by Bondi (1952) and Parker (1958), respectively. While the problem of spherical accretion is mainly of academic interest, a spherical stellar wind is a qualitatively reasonable first approximation to the problem of the expansion of a nonconfined hot stellar atmosphere. Consider a steady flow in spherical geometry under the effect of grqvitational field of a massive central object. The governing equations are: 1 ∂ 2 (r ρv) = 0 r2 ∂r ⇒ 2 1 dρ 1 dv + + =0 r ρ dr v dr dv c2s dρ GM 1 dp GM 1 dv 2  GM v =− − 2 − 2 =− − 2 = c2s  + dr ρ dr r ρ dr r v dr r r   For simplicity, consider the isothermal expansion/accretion, where c2s = RT is a constant. (The case of polytropic expansion/accretion is left as an exercise.) Thus,   c2s  dv 2c2s GM v−  = − 2 . v dr r r   Clearly, v = cs can be realized only at r = rc ≡ Hydrodynamics – Spring 2010 GM 2c2s 164 The equation is easily integrated as r cs !2 cs !2 2GM + C, − ln = 4 ln + v v rc r c2s where C is a constant of integration. The solutions are depicted below. Solutions I and II are double valued and, hence, unphysical. Solution III would correspond to supersonic flow everywhere, and this does not correspond to either the wind or the accretion problem. Solution IV is everywhere subsonic. It would correspond to a stellar breeze in the outflow case. Transsonic solutions are obtained setting C = −3. Solution V is Parker’s solar (or stellar) wind and the other transsonic solution corresponds to Bondi’s spherical accretion. Hydrodynamics – Spring 2010 165 Waves and instabilities Hydrodynamics – Spring 2010 166 Perturbation analysis Consider a ball at rest on a horizontal point of a surface, as depicted below. All cases represent equilibrium states. In the leftmost case, the equilibrium is clearly unstable: even a slight perturbation of the position of the ball will lead to a destruction of the equilibrium and the ball will roll away from the hill top. In the case in the middle, the equilibrium is stable: if displaced from its equilibrium position and let to move, the ball will eventually return to the bottom of the valley. Before that, it will most likely overshoot to the other side of the equilibrium position and oscillate aroud the equilibrium position until frictional forces bring it back to rest. The equilibrium in the rightmost case is stable against small perturbations but unstable against large perturbations. Such a state is called conditionally stable or metastable. Obviously, it is not enough for a system to be an equilibrium solution of the equations of motion: in order for the system to represent a state to be found in Nature by any reasonable probability, the solution has to be stable against perturbations. The most robust equilibria are those that are stable against arbitrary perturbations. A stability analysis is also necessary in hydrodynamics to find out if any found mathematical equilibrium state is going to represent a physical equilibrium or not. Hydrodynamics – Spring 2010 167 Perturbation analysis for fluids We already introduced the linearization technique, when discussing the acoustic waves. The same technique can be applied in analysing the possible unstable growth of perturbations around any equilibrium solution of the HD equations. As for the acoustic waves, one can then apply Fourier analysis (i.e., write down the perturbed part of the solution as ∝ exp{i(k · x − ωt)}) and obtain, using the HD equations, a dispersion relation for the perturbations in form ω = ω(k). If the solution contains a positive imaginary part, i.e., ω(k) = ωr (k) + iωi(k) with ωi > 0 for some values of k then the perturbations at that particular wavenumber grow exponentially as ∝ exp{ωi(k)t} and the underlying equilibrium solution is unstable. If ωi < 0, then the perturbations are damped and the equilibrium is stable against perturbations at that wavenumber. In some cases we find that the equilibrium is stable or unstable, depending on the values of some parameters of the equilibrium (conditional stability). As an example, we mention the onset of convection in a vessel heated from below, once the temperature gradient exceeds a threshold value. Of course, one can justify the linearization of HD equations only if the perturbations are small. Thus, the linear perturbation analysis can only be applied to the initial evolution of unstable perturbations. Eventually, the perturbations grow and we may finally end up with an irregular state (in space and time) usually referred to as turbulence. Hydrodynamics – Spring 2010 168 Acoustic waves in a viscous fluid Let us consider a viscous Newtonian fluid, for which the complete N–S equation in the absence of external forces reads dv = −∇p + (ζ + 13 µ)∇(∇ · v) + µ∇2v. ρ dt Consider the evolution of adiabatic perturbations on top of a static equilibrium solution: ρ = ρ0 +ρ1(x, t), p = p0 + p1(x, t) and v = v 1(x, t). Linearizing the N–S equation and eliminating the pressure using p1 = c2s ρ1 and the density using the equation of continuity gives ∂v 1  µ 2 ∂v 1 ζ + 13 µ  ∂ 2v1 2 ∇ · ∇(∇ · v ) + = c + ∇ ∇ 1 s ∂t2 ρ0 ∂t ρ0 ∂t   where c2s = γp0/ρ0. Seeking for longitudinal solutions in form v = v1ei(kx−ωt)ex gives 2 ω = c2s k 2 ζ + 43 µ 2 k ω, i.e., −i ρ0 v u u u u u t  2 4 ζ + 43 µ 2 ζ + 3µ  k − i k ω = ±csk 1 −  2ρ0cs 2ρ0 Clearly, the equilibrium is stable (ωi < 0 at all k ). Acoustic waves also become dispersive in a viscous fluid and cease to propagate (ωr = 0) at k > 2ρ0cs/(ζ + 43 µ). At large wavelengths (low k ) the results derived from Euler equation are recovered. Hydrodynamics – Spring 2010 169 Convective instability Consider a perfect gas in hydrostatic equilibrium in a uniform gravitational field, g = −gez . Thus, p(z) and ρ(z) decrease upward with z . Assume that a blob of gas is displaced vertically. Let the original pressure and denisity of the blob (and the surrounding gas) be p and ρ, respectively. Let the external pressure and density at the new position be p′ and ρ′. Because acoustic waves transmit momentum in the fluid quite fast, it is reasonable to assume that the pressure inside the blob has equilibriated with the surroundings. However, heat transport via conduction is usually a slower process, so it is reasonable to assume that the density of the blob behaves adiabatically. Denote the density at the new position by ρ∗.   ρ∗ = ρ   ′ 1/γ p  p Thus, ρ∗ = ρ + On the other hand, dρ ρ′ = ρ + ∆z dz  ρ’, p’ ρ, p ρ, p dp ∆z. dz ρ dp ∆z γp dz ρ=p/RT =  and p′ = p + ρ *, p’   ρ dp ρ dT  ∆z ρ+ − p dz T dz   1 ρ dp ρ dT  ∴ ρ − ρ = − 1 −  ∆z + γ p dz T dz ′ Hydrodynamics – Spring 2010 ∗ 170 We got     1 ρ dp ρ dT  + ∆z ρ′ − ρ∗ = − 1 −  γ p dz T dz Clearly, if the difference is positive, the blob will be buoyant and the system is unstable. Thus, we derive the stability criterion of the atmosphere as   1 T dp dT < 1 −  dz γ p dz Schwarzschild criterion. If the temperature profile is steeper than this limit, then the atmosphere will be unstable to convection, i.e., up- and downward movement of gas. The Schwarzschild criterion is important to consider when discussing models of the internal structure of stars. An eq. of motion for the blob is 2 ∗ d ∆z ρ dt2 = −(ρ∗ − ρ′)g where we neglect any viscous forces arising from the movement of the blob relative to the surrounding gas. Thus, d2∆z + N 2∆z = 0, 2 dt where v     u u g dT 1 T dp   t − 1 −  N =u T dz γ p dz is the Väisälä–Brunt frequency. Clearly, if the atmosphere is stable, the blob ends up in oscillatory motion. A full analysis of this motion will reveal that the stably stratified atmosphere will support internal gravity waves in addition to ordinary acoustic waves. Hydrodynamics – Spring 2010 171 Rayleigh–Bénard convection Experimentally (Bénard 1900) it is found that a layer of liquid heated from below becomes unstable to convection, which occurs in cells with dimensions of the order of the depth of the layer, as the temperature gradient increases above a certain value. Let us, next, try to understand this theoretically. This requires us to analyze the problem using full, nearly incompressible N–S equation (Rayleigh 1916). Hydrodynamics – Spring 2010 172 Rayleigh–Bénard convection - set up Consider a layer of liquid between z = 0 and z = d, in a gravitational field g = −gez , heated from below. Assume that the liquid is nearly incompressible, i.e., that its density does not change in reaction to varying pressure, but decreases as a result of increasing temperature. The energy equation for a nearly incompressible fluid reads   ∂ǫ ρ  + v · ∇ǫ = ∇ · (K∇T ), ∂t where the internal energy is ǫ = cpT and cp is the specific heat. Thus, ∂T + v · ∇T = κ∇2T, ∂t where κ = K is the thermometric conductivity. ρcp Let Tb and Tt be the temperatures at the bottom and the top of the layer. Thus, the equilibrium solution is Tb − Tt . T0(z) = Tb − βz, where β = d The corresponding density is ρ0(z) = ρb(1 + αβz), where α is the coefficient of volume expansion with temperature. The equilibrium pressure satisfies the hydrostatic balance equation dp0 = −ρ0(z) g. dz Hydrodynamics – Spring 2010 173 Next, introduce perturbations around the equilibrium state. If T = T0 + T1, then ρ = ρ0 − ρbαT1. Write p = p0 + p1 and subsitute in the N–S equation to get   ∂v 1 + v 1 · ∇v 1 = −∇(p0 + p1) + (ρ0 − ρbαT1)g + µ∇2v 1. (ρ0 − ρbαT1)  ∂t Clearly, ∇p0 = ρ0g . We may also neglect the non-linear terms v 1 · ∇v 1 and ρbαT1 ∂v 1/∂t, and replace ρ0 by the constant ρb on the left-hand side, as it multiplies a small term (Boussinesq approximation). With these simplifactions, we get ρb ∂v 1 = −∇p1 − ρbαT1 g + µ∇2v 1 ∂t (LNS) Substituting the perturbed variables in the temperature equation and linearising we get ∂T1 = βv1z + κ∇2T1. ∂t (T) There are five variables (v 1, p1 and T1) and four equations in (LNS–T), but the system could be closed by using the kinematic constraint ∇ · v 1 = 0. However, it is also possible to eliminate extra variables by taking a curl of (LNS) twice and then considering the z component:   2 ∂ 2T1  ∂ 2  ∂ T1 4  + ν∇ v1z . (∇ v1z ) = αg  2 + ∂t ∂x ∂y 2 (V) Equations (T) and (V) form a closed system with v1z and T1 as the only unknowns. Hydrodynamics – Spring 2010 174 Rayleigh–Bénard convection – marginal stability We derived ∂T1 = βv1z + κ∇2T1 ∂t   2 2 ∂ T ∂ T ∂ 2 1 1 4  + ν∇ v1z .  + (∇ v1z ) = αg  ∂t ∂x2 ∂y 2 We can try to find solution of these equations in form v1z = W (z) exp(σt + ikxx + iky y) T1 = θ(z) exp(σt + ikxx + iky y) We expect the solution to be periodic in the xy plane and allow for a general time dependence (σ can, in principle, be complex). The functions W (z) and θ(z) then satisfy ordinary differential equations obtained by substituting the trial forms into the equations. Of course, if σr > 0, then the perturbation grows and if σr < 0 it decays. The limiting case is σr = 0, which is called marginally stable. In the following, we shall assume that σ = 0 to find a criterion for the onset of the instability. This gives   2  d θ = 0 βW + κ  2 − k 2 dz  2 2  d 2  W = 0, −αgk θ + ν  2 − k 2 dz Hydrodynamics – Spring 2010 where k 2 = kx2 + ky2. 175 Thus,  2 3 d 2  W = −αβgk W. − k2 2 dz Substituting z = z ′d and k = k ′/d, i.e., using the layer depth as the unit of distance, gives  νκ     2 3 d ′2 ′2  W = −Rk W, − k dz ′2 αβgd4 is the (dimensionless) Rayleigh number of the system. where R = κν For a correct physical solution of this sixth-order equation, one needs six boundary conditions. Four of them are immediately obtained as both W = 0 and θ = 0 (i.e., [(d/dz ′)2 − k ′2]2W = 0) have to be required both at z = 0 and z = d. The remaining two boundary conditions are obtained from the requirement that the tangential stresses at a free surface have to vanish, leading6 to d2W/dz ′2 = 0, and that the no-slip condition has to be maintained at a solid boundary, leading to dW/dz ′ = 0. The simplest solution is obtained, when both boundaries are considered as free surfaces (e.g., liquid floating on top of a heavier liquid), whence d2W/dz ′2 = 0 is taken to apply at both surfaces. Such boundary conditions are satisfied only by W (z ′) = W0 sin πnz ′, n ∈ N, with R and k ′ related by (n2π 2 + k ′2)3 . R= k′2 6 see Chandrasekhar’s Hydrodynamic and hydromagnetic stability, Dover 1981, p. 21–22. Hydrodynamics – Spring 2010 176 We obtained as a condition of marginal stability that (n2π 2 + k ′2)3 R= k′2 αgd3∆T . where R = κν Intuitively, an increase of the temperature difference yieds to instability, so values of R that fall above this value for given n and k correspond to an instability. Clearly, the most unstable situation corresponds to n = 1. Thus, the curve (π 2 + k ′2)3 R= k′2 separates the unstable parameter range (above the curve) from the stable range (below the curve) in the k ′R-plane. The separating curve has a minimum value at R Unstable π k ′ = kc′ = √ 2 4 with R = Rc = 27π ≈ 657.5. 4 For a Rayleigh number above this limit, there always exists an unstable range of wavenumbers. One can find that for two rigid boundaries, the critical Rayleigh number is 1707.8 and for one rigid and one free boundary, it is 1100.7. Most experiments yield Rc ≈ 1700. Hydrodynamics – Spring 2010 Rc Stable kc k 177 Perturbations at a two-fluid interface Consider a horizontal, planar interface of two fluids (which may be in uniform horizontal motion) in a vertical gravitational field. (Example: the surface of a lake or a pond.) Let us study the evolution of small perturbations of the bounding surface. We will consider the following simplifications of the problem: • Both fluids are assumed incompressible • Both fluids are assumed ideal • Surface tension is neglected As we assume incompressible and ideal fluids, Kelvin’s theorem of vorticity applies. Thus, assuming that the steady-state flow is irrotational also the perturbations of the fluid velocity have to remain such. We can, therefore, assume that the velocity on both sides of the interface is obtained from a scalar potential, v = −∇φ. Thus, v2 p ∂φ + ∇ = −∇ − ∇Φ, −∇ ∂t 2 ρ where Φ = gz is the gravitational potential. Integrating the equation gives − ∂φ 1 2 p + v + + Φ = F (t), ∂t 2 ρ where F (t) is an arbitrary function of time. Because of the assumption of irrotational velocity (at t = 0) this equation, resembling Bernoulli’s principle, applies throughout the two fluids. Hydrodynamics – Spring 2010 178 U’, ρ’ z ξ1 Consider a two-fluid interface (densities ρ and ρ′) with uniform flows in the x direction (speeds U and U ′), as depicted in the figure. This configuration satisfies the steady-state ideal HD equations, as long as pressure is continuous at the interface, i.e., x     p0(0) + gρ′z, z > 0 p0(z) =   p (0) + gρz, 0 U, ρ z<0 , where z = 0 has been chosen as the position of the interface in steady state. Consider small perturbations of the interface and denote its perturbed position by z = ξ1(x, t). The aim is to find out whether ξ1 grows with time, decays or oscillates. Velocity potential below the surface is φ = −U x + φ1(x, z, t), where the first term gives the uperturbed part. The perturbed part satisfies ∇ · v 1 = 0, i.e., ∇2φ1 = 0. Above the surface, φ′ = −U ′x + φ′1(x, z, t) with ∇2φ′1 = 0. Hydrodynamics – Spring 2010 179 For fluid elements at an infinitesimally small distance from the interface, we can write vz = dξ1/dt = ∂ξ1 ∂ξ1 + vx . ∂t ∂x Thus, to the linear order in small quantities, ∂ξ1 ∂φ′1 ∂ξ1 ∂ξ1 ∂φ1 ∂ξ1 = +U and − = + U′ − ∂z ∂t ∂x ∂z ∂t ∂x at z = 0. Consider the following Fourier components ξ1 = Aei(kx−ωt) φ1 = Cei(kx−ωt)+lz ′ φ′1 = C ′ei(kx−ωt)−l z Apply the Laplace equations to the velocity potentials to get −k 2 + l2 = 0 2 −k 2 + l′ = 0 which give l = l′ = |k|, where the signs have been chosen so that the perturbations vanish at infinity. Hydrodynamics – Spring 2010 180 Consider k > 0. Applying the boundary conditions at the interface then gives, −kC = i(kU − ω)A kC ′ = i(kU ′ − ω)A Dividing the equations gives kU − ω C = C ′ kU ′ − ω so we clearly need one more equation to find out the ratio of the amplitudes of the velocity potentials in order to arrive to a dispersion relation for the perturbations. This is provided by the requirement that pressure must be continuous accross the interface. − We have  − p = −ρ    ∂φ1 v  + ρF (t) + + gz  ∂t 2 p′ = −ρ′  −  2 ∂φ′1 ∂t + ′2  v  ′ ′ + gz   + ρ F (t) 2 Comparing with the steady-state solution (that needs to be recovered far from the interface where v 1 = 0), we get 2 ρF (t) = p0(0) + 21 ρU 2 and ρ′F ′(t) = p0(0) + 21 ρ′U ′ . Hydrodynamics – Spring 2010 181 Thus, applying p = p′ at the boundary, we have  − ρ 2  2  ∂φ′1  ′2 ′2 ∂φ1 v − U v −U  ′  = ρ  + + gξ1 + + gξ1  − ∂t 2 ∂t 2 and, using v 2 = (U − ∂φ1/∂x)2 + (−∂φ1/∂z)2 ≈ U 2 − 2U (∂φ1/∂x), this gives   ′ ′ ∂φ1 ∂φ1 ′  ∂φ1 ′ ∂φ1    −U + gξ1 = ρ − −U + gξ1 ρ − ∂t ∂x ∂t ∂x   at z = 0. Thus, ρ[−i(kU − ω)C + gA] = ρ′[−i(kU ′ − ω)C ′ + gA]. Substituting C and C ′ as a function of ω , k and A from above gives ρ(kU − ω)2 + ρ′(kU ′ − ω)2 = kg(ρ − ρ′) as the dispersion relation for the perturbations. Solving for the phase speed gives ′ ′ v u u u u t ω ρU + ρ U g ρ − ρ′ ρρ′(U − U ′)2 = ± − k ρ + ρ′ k ρ + ρ′ (ρ + ρ′)2 This general equation can be applied to several special cases of interest to derive the appropriate dispersion relations for several two-fluid-interface phenomena. Hydrodynamics – Spring 2010 182 Surface gravity waves Consider two fluids at rest with a lighter fluid lying above a heavier one. Thus, U = U ′ = 0, ρ > ρ′ and v u u u t ω g ρ − ρ′ =± . k k ρ + ρ′ Clearly, the perturbation of the surface moves as a wave (surface wave) with a phase velocity that depends on the wavenumber. Thus, surface waves are dispersive such that the longest waves have the highest phase speeds. If one considers a case where ρ′ ≪ ρ, like air on ρ’ < ρ top of water, the equation simplifies to t 1 > t0 ω/k r ω = ± gk. z x This expression neglects two effects that you are requested to take into account as exercises: ρ • Surface tension will somewhat increase the propagation speed. Intuitively this is clear, as the restoring force (trying to restore the “flatness“ of the surface) will be greater as a result of this extra tension acting on the separating surface. • A finite depth h of the lower fluid will modify the phase speed of waves with k −1 & h. This provides an upper limit to the phase speed. Hydrodynamics – Spring 2010 183 Rayleigh–Taylor instability Let us still consider fluids at rest. Now assume, however, that the heavier fluid is on top of the lighter one, i.e., ρ < ρ′. For a perfectly flat, horizontal separating surface, the configuration is still an equilibrium solution of the HD equations. We, however, know from experience that the heavier fluid should fall below the lighter one, so that the equilibrium with the heavy fluid on top should be unstable. This we also find from the dispersion relation: v u u u t ρ’ > ρ t 1 > t0 z x ρ v u u u u t ρ − ρ′ |ρ′ − ρ| ω = ± gk = ±i gk ρ + ρ′ ρ + ρ′ which corresponds to purely imaginary frequency, i.e., a purely growing (and a decaying), non-propagating perturbation. This instability is called Rayleigh–Taylor instability. We note that because of the equivalence principle, a situation corresponding to Rayleigh–Taylor instability can also be generated in an accelerating frame of reference. If there is no gravitational field but the interface is propagating in the −z direction with a decelerating speed in the inertial frame, we find a heavy fluid pushing against a light one, and the same instability analysis applies. Hydrodynamics – Spring 2010 184 The Crab nebula (and many other supernova remnants) provide us with examples of such a two-fluid interface. The ejected gas from the explosion tends to pile up in a thin layer just behind the interface between the ejecta and the surrounding (shocked) gas. As the interface is decelerating, this creates an effective outward g -field in the frame of the interface. We, thus, have a situation, where the lighter fluid is below the heavy on in this effective g -field and the situation is Rayleigh–Taylor unstable. Hydrodynamics – Spring 2010 185 Kelvin–Helmholtz instability Finally, let us consider the case with non-zero velocities. Let us restrict our attention to the Rayleigh– Taylor stable case (ρ > ρ′). Thus, in the coordinate system co-moving with the lower fluid v u u u u t ρ ∆U ρ − ρ′ ρρ′k 2∆U 2 ω= k ± gk − , ρ + ρ′ ρ + ρ′ (ρ + ρ′)2 ′ where ∆U = U ′ − U is the velocity of the upper fluid with respect to the lower fluid. (Note that this relation applies also if U and U ′ are not co-aligned, whence ∆U = |U ′ − U |.) It is clear that if ρ’ < ρ 2 ρρ′k∆U 2 > g(ρ2 − ρ′ ) ω is complex. The perturbation is both propagating at speed ωr /k = ρ′∆U/(ρ + ρ′) and growing exponentially at rate U’ − U t 1 > t0 z x ρ g ρ − ρ′ ρρ′∆U 2 − ωi = k (ρ + ρ′)2 k ρ + ρ′ This instability is called Kelvin–Helmholtz instability. Note that in the absence of gravitational field, an ideal two-fluid interface with a tangential velocity difference (i.e., a tangential discontinuity) is always Kelvin–Helmholtz unstable. v u u u u t Hydrodynamics – Spring 2010 186 Beyond the linear theory of waves and instabilities Our analysis of waves and instabilities has been restricted by the assumption that perturbation are small so that the equations can be linearised around the equilibrium solution. Only acoustic waves were treated with a simplified account of non-linearity, and shock waves, i.e., thin transition fronts were shown to be produced as a result of non-linear steepening. Dispersive waves (such as surface waves) will behave oppositely: their wave packets will be dispersed as the different Fourier components have different velocities. Sometimes a balance between dispersive effects and non-linear effects may be found. The resulting waves, called solitary waves, retain their shape during propagation. If solitary waves are, in addition, able to survive collisions they are called solitons. Consider a surface wave propagating in a shallow water of depth h ≪ k −1. In the linear limit, it has the disperison relation (exercise)   r 1 2 2 3  ω = gk tanh kh ≈ k gh 1 − k h ≡ c0k − σk with c0 = gh and σ = 16 c0h2 6 r r A wave fulfilling this dispersion relation would have to satisfy the differential equation ∂ 3v ∂v ∂v + c0 + σ 3 = 0. ∂t ∂x ∂x Adding the non-linear term v(∂v/∂x) and transforming to a coordinate system moving with the linear Hydrodynamics – Spring 2010 187 wave speed c0 [i.e., using X = x − c0t] gives ∂v ∂v ∂ 3v +v +σ = 0, ∂t ∂X ∂X 3 which is the celebrated Korteweg–de Vries equation that admits a soliton solution of form v= 3v0 . cosh2[ 21 (v0/σ)1/2(X − v0t)] Such soliton solutions are localized wave pulses that retain their shapes. Their amplitude, 3v0, propagation r r speed, (c0+) v0, and width, 2 σ/v0 = h 2c0/3v0, are all related to each other. More on solitons during the student seminar. In the non-linear regime of instabilities two kinds of behavior may be expected. The instability may be limited by non-linear effects so that the amplitude of the perturbation stops growing after a certain limit, often depending on the system parameter driving the instability. In this case a new stable laminar solution may be achieved, replacing the original equilibrium. Rayleigh–Benard convection at Rayleigh numbers close to the critical value provides an example. It is found that the amplitude of the convection velocity increases linearly with R − Rc so when the temperature gradient is increased, more and more energy is transported by convection. Another possibility is that no laminar state will be achieved but that the unstable system develops an irregular flow (in space and time) that can be analysed only in a statistical sense (e.g., flow past a body at the largest Reynolds number or convection at the largest Rayleigh numbers). These kind of systems are referred to as being turbulent. Hydrodynamics – Spring 2010 188 Turbulence Hydrodynamics – Spring 2010 189 The need for a statistical theory Let us consider the trajectories of a dynamical system in phase space around an unstable equiq2 librium point P . If the system is excactly at P at time t = 0 it will stay there. If, on the other hand, the system is even slightly displaced from this point, the distance from the point will start to grow exponentially and after a finite amount of time the system lies in a very different region of phase space. It will, therefore, be practically impossible to predict the state of the system q1 deterministically after a finite amount of time. Thus, for an unstable fluid system, we expect that the growth of perturbations may lead to a loss of q3 our prediction capability. One often encounters in practice systems in which fluid velocities seem to vary randomly in space and time. Such a state of a fluid is called turbulence. In addition to hydrodynamically unstable fluids, stable fluids stirred with random forces may develop turbulence. P Since it is impossible to develop a deterministic theory of turbulence, we wish to have a statistical description of its properties based on averages. Ensemble averiging would be most appropriate theoretically, observationally temporal averaging is most easily realised. For ergodic systems, the two are equivalent. We shall assume ergodicity and denote statistical averaging by an overbar. Hydrodynamics – Spring 2010 190 Correlation length The velocity in each point can be written as v = v̄ + v ′, where v̄ and v ′ are the statistical average and the fluctuating part. Obviously, v ′ = 0. Thus, a more complicated average of the fluctuating field is needed to describe its statistical properties. Maybe the simplest quantity to consider, proposed by Taylor (1935), is v ′(x, t) · v ′(x + r, t). At r = 0, this equals v ′2(x, t) which is proportional to the average kinetic energy density of the fluctuations. At r → ∞, the fluctuating velocities are uncorrelated, so we get lim v ′(x, t) · v ′(x + r, t) = v ′(x, t) · v ′(x + r, t) = 0 r→∞ Thus, v ′(x, t) · v ′(x + r, t) has a substantial value only if r . a finite distance. This distance is called the correlation length of the turbulence. Thus, correlation function like v ′(x, t) · v ′(x + r, t) contain the information on the strength and correlation length of the turbulence and are appropriate measures of its properties. Hydrodynamics – Spring 2010 191 Correlation tensor The velocity correlation tensor of the turbulent flow is defined as R(r; x, t) = v ′(x, t)v ′(x + r, t); i.e., Rij (r; x, t) = vi′ (x, t)vj′ (x + r, t) In general, it has nine components but symmetries can reduce the number of independent components drastically. From now on, we will not denote the depence on time explicitly if both velocities are evaluated at the same time. We already concluded that in case of unstable equilibria, it is impossible to prescribe the instantaneous values of the fluctuating fields. But what about the velocity correlations? A proper dynamic theory of turbulence would give the values of Rij (r; x) (and higher-order correlation function such as vi(x1)vj (x2)vk (x3)) from given initial and boundary conditions for the mean variables and statistics of the fluctuations. Such theory does not exist yet. In the following, we will consider incompressible, homogeneous and isotropic turbulence. The second assumption means that we assume the correlation tensor components not to depend on x. The third assumption means that the medium has no preferred direction, i.e., that the only vector the correlation tensor can depend on is r . It implies also that v̄ = 0. Homogeniety and incompressibility imply that Rij (r) = Rij (r; x) = Rij (r; x − r) = vi(x − r)vj (x) = Rji(−r) ∂Rij ∂vj (x + r) ∂Rji(r) = vi(x) = vi(x) (∇ · v)(x + r) = 0 ⇒ =0 ∂rj ∂rj ∂rj Hydrodynamics – Spring 2010 192 The form of the correlation tensor in isotropic turbulence can be derived as follows. Consider a Cartesian coordinate system K̃ , where one of the base vectors, say ẽ3, is directed along r = rẽ3. Denote the correlation function R̃33 = ṽ3(x)ṽ3(x + rẽ3) = 13 v 2f (r) with f (0) = 1. Because of symmetry, the correlation functions of the velocity components perpendicular to r have to agree, i.e., R̃ij = ṽi(x)ṽj (x + rẽ3) = 31 v 2g(r) with i = j = 1, 2 and g(0) = 1. Off-diagonal components of R̃ij vanish in this coordinate system as we may assume that the fluctuating velocity components are independent of each other. Thus,   g(r) 0 0  v 2  (R̃ij ) =  0 g(r) 0  .  3  0 0 f (r) Transforming to a Cartesian coordinate system with another orientation of axes ei, the tensor components transform as Rij = Mir MjsR̃rs, where Mij = ei · ẽj is an orthogonal tensor, i.e., M · MT = I ⇔ Mik Mjk = δij . Thus, Rij = 13 v 2[(Mi1Mj1 + Mi2Mj2)g(r) + Mi3Mj3f (r)] = 31 v 2{Mik Mjk g(r) + Mi3Mj3[f (r) − g(r)]} | Mik Mjk = δij = 13 v 2{δij g(r) + Mi3Mj3[f (r) − g(r)]} As Mi3 = ei · ẽ3 = ei · (r/r) = ri/r, we have Rij = 1 2 3v ) ri rj δij g(r) + 2 [f (r) − g(r)] r ( Hydrodynamics – Spring 2010 i.e., R = 1 2 3v rr ) g(r)I + [f (r) − g(r)] 2 . r ( 193 We can derive a connection between f (r) and g(r) using the incompressibility condition, i.e,     rirj  df dg  ∂r ∂ rirj ! ∂Rij 1 2  dg ∂r − = 3v  δij + 2 + [f (r) − g(r)] 0= ∂ri dr ∂ri r dr dr ∂ri ∂ri r2     !   dg ri df r r dg ∂ r r r k j i j i δki 2  δij + 2  −  + [f (r) − g(r)] = 13 v 2  dr r r dr dr r ∂ri r    !  dg rj  r r dg df ∂ r r r i i j k j 1 2   = 3v  + 2 − + [f (r) − g(r)]δki dr r r dr dr r ∂ri r2      df rj δ r + r δ r r r ik j k ij k j i  + [f (r) − g(r)]δki  − 2 = 31 v 2  dr r r2 r3 r            r i r j r i   (δki δik + 1)rj 1 2 df rj   − 2 + [f (r) − g(r)] = 3v    dr r r2 r3 r  = 31 v 2   df rj r [(I · I)kk + 1]rj  + [f (r) − g(r)]  − 2 dr r r2 r   j   2   = 13 v 2       rj df 2[f (r) − g(r)] +   dr  r  r f giving r df g(r) = f (r) + 2 dr    1.0   r df  rirj df  . − ∴ Rij (r) = 13 v 2 δij f (r) + 2 dr 2r dr  r Thus, the correlation tensor in incompressible, homogeneous, isotropic turbulence is fully determined by the longitudinal correlation function, f (r). It should look as sketched above. Hydrodynamics – Spring 2010 194 Turbulent diffusion Convection sets on in a fluid heated from below to ensure efficient transport of heat through the fluid. Under normal conditions, convection is orders of magnitude more efficient method of heat transport than conduction. We should ask, therefore, what is the role of the turbulent velocity fields in heat transport. More generally, it is of interest to study the effect of turbulent velocity fields to transport phenomena. Let us look at the problem of heat transport again. The energy equation in an incompressible fluid could be written in form ∂T + v · ∇T = κ∇2T, ∂t where κ = K is the thermometric conductivity. ρcp If κ ≈ 0 (its smallness will be discussed later), we find that temperature is a passively advected scalar in the fluid. This means that its value for each fluid element is fixed. Thus, if the fluid elements should diffuse as a result of turbulent velocity fields, we find that turbulence would provide a diffusion term in the energy equation even without heat conduction. Let us consider the displacements of the fluid elements, ξ(t) = ˆ t v L(t′) dt′, 0 where v L(t) = v(x0 + ξ(t), t) is the velocity of the fluid element, intially at x = x0, in the Lagrangian description and v(x, t) is the turbulent velocity field. (For simplicity, we’ll assume v̄ = 0.) Hydrodynamics – Spring 2010 195 Obviously, the mean displacement is zero, ξ = 0. However, ξ2 = ˆ t ˆ t vL 0 = 0 ′ dt (t′) dt′ ˆ t 0 · ˆ t vL 0 (t′′) dt′′ = t ˆ 0 dt′ t ˆ dt′′ v L(t′) · v L(t′′) 0 dt′′ v L(t′) · v L(t′′) is non-zero. Here, v L(t′) · v L(t′′) is a measure of correlation between the velocities at the positions of a fluid element at times t′ and t′′. If the turbulent velocity field is in a steady state, the value of the correlation may depend only on the time difference, τ = t′′ − t′. Thus, we may assume v L(t′) · v L(t′′) = v 2R(τ ). Clearly, R(0) = 1 and R(τ ) should have a substantial value only if τ is smaller than or of the order of a correlation time, τ . τcor. We may also take R(τ ) = R(−τ ) if we assume homogeniety of turbulence. Thus, ˆ t ˆ t ˆ t−t′ ˆ t ξ 2(t) = dt′ v 2 dτ R(τ ). dt′′ R(τ ) = dt′ v 2 0 0 −t′ 0 For t ≪ τcor, R(τ ) ≈ 1 and ξ 2 ≈ v 2t2. Now, consider times t ≫ τcor. Thus, for overwhelmingly most values of t′, the limits of the inner integral may be replaced with ±∞ and we get ξ 2(t) = ˆ 0 Hydrodynamics – Spring 2010 t ′ dt v 2 ˆ +∞ −∞ dτ R(τ ) = 2t v 2 ˆ +∞ dτ R(τ ). 0 196 We have found that for times t ≫ τcor ξ 2(t) = 6DTt, where DT = 1 3 v2 ˆ 0 +∞ dτ R(τ ) =: 13 v 2 τcor The linear time dependence of the mean square displacements is a signature of a diffusive process. Consider a diffusion process of test particles with density n, which all lie in the origin at time t = 0. Thus,   ∂ D ∂n ∂n = D∇2n = 2 r2  (spherical symmetry assumed) ∂t r ∂r ∂r ´∞ 2 r n(r, t) d3r 0 2 ´ ξ (t) = ∞ . 3r n(r, t) d 0 Multiplying the diffusion equation by r2 and integrating over all space gives ˆ ∞ ˆ ˆ     ∂ ∂ ∂n ∂n ∂ 2 2 3 2 2 3 r  d r = D r  dr 4πr r nd r = D ∂t ∂r ∂r ∂r ∂r 0 ˆ ˆ ∞ ˆ ∞ ∂n 4πr2n dr = 6D n d3r 4πr3 dr = 6D = −2D ∂r 0 0 ˆ ˆ ∴ r2n d3r = 6Dt n d3r i.e., ξ 2(t) = 6Dt We may, therefore, identify DT = 31 v 2τcor as the turbulent diffusion coefficient. Hydrodynamics – Spring 2010 197 Properties of turbulent diffusion Let us analyze turbulent diffusion more closely. In turbulent diffusion, diffusion coefficient does not depend on molecular properties or on quantity being advected (passively), but only on the statistics of the velocity fluctuations. same coefficient may be used as the effective thermometric conductivity (related to energy and kinematic viscosity (related to momentum transport). However, this statement does not in a general case, if the quantity being transported can react back on the fluid. the scalar Thus, the transport) quite hold Molecular diffusion coefficients in a dilute gas are all of the order D, ν, κ ∼ λhui √ where hui is the averge molecular speed in the rest frame of the fluid, λ = 1/( 2σn) is the mean free path of the molecules and σ and n are the collision cross-section and the density of the collision targets (other molecules of the gas). However, turbulent diffusion is often the dominant mode of diffusion, up to the extent that in a turbulent flow molecular diffusion may often be completely negligible. Thus, for air in room temperature, for example, λ ∼ 3 µm and hui ∼ 470 m s−1 giving D ∼ 10 cm2 s−1. The time it takes for heat to diffuse across a room of L = 10 m via molecular diffusion is, thus, L2 ∼ 1 day. τD ∼ D This would make the design of heating systems quite difficult. Fortunately, convection by turbulent velocity fields is a much faster process. Hydrodynamics – Spring 2010 198 The mean equations Let us, then, estimate the effect of turbulent diffusion on momentum transport. Consider an incompressible fluid obeying N–S equation   ∂ ∂vi  ∂  −p δij − ρvi vj + µ  (ρvi) = ρFi + ∂t ∂xj ∂xj Break up the quantities into mean part and fluctuating part, i.e., vi = v̄i + vi′ ; p = p̄ + p′; Fi = F̄i + Fi′. Substitute in N–S and average. All linear terms in fluctuating quantities vanish in averaging. The only non-zero contribution from fluctuations comes from vivj = (v̄i + vi′ )(v̄j + vj′ ) = v̄iv̄j + vi′ vj′ . Thus,   ∂  ∂ v̄i  ∂ ′ ′ −p̄ δij − ρv̄i v̄j − ρvi vj + µ , (ρv̄i) = ρF̄i + ∂t ∂xj ∂xj which was first derived by Reynolds (1895). This Reynolds equation is very similar to N–S equation, containing only the extra term ∝ vi′ vj′ . This term is known as the Reynolds stress. Hydrodynamics – Spring 2010 199 One option for evaluating the Reynolds stress is to use ∂vj′ ∂ ′ ′ ∂vi′ ′ ′ v + vi vv = ∂t i j ∂t j ∂t and ∂vi′ ∂vi ∂ v̄i = − ∂t ∂t ∂t and substitute the time derivatives from the N–S and Reynolds equations, respectively. This procedure, however, clearly leads to an equation that involves the three-correlations vi′ vj′ vk′ . If one tries to derive an equation for their evolution, that involves four-correlations etc. Thus, we again meet a closure problem: this chain of correlation function needs to be truncated by applying some sort of consitutive relation. One option is to apply  Thus, for constant DT  ∂ v̄i ∂ v̄j    vi′ vj′ = −DT  + ∂xj ∂xi   ∂ v̄i  ∂  ∂ −p̄ δij − ρv̄i v̄j + ρ(DT + ν) , (ρv̄i) = ρF̄i + ∂t ∂xj ∂xj and, thus, DT is the turbulent kinematic viscocity of the fluid. As already deduced, in practice ν ≪ DT usually holds. It has to be kept in mind that the adopted crude closure scheme does not always lead to qualitatively correct conlusions, but that the Reynolds stresses sometimes have to be modeled by other types of terms than diffusive. Hydrodynamics – Spring 2010 200 Energy spectrum of turbulent fluctuations Let us, then, return to the problem of developing a statistical description of turbulent fluctuations. We already introduced the correlation tensor Rij (r) = vi(x)vj (x + r) = v2 3   f (r) − g(r) rirj  ; 2 r g(r) δij +   with g(r) = f (r) + r df . 2 dr Let us, instead of Rij , consider its Fourier transform 1 Φij (k) = (2π)3 The inverse transform gives Rij (r) = ˆ ˆ Rij (r)eik·r d3r. Φij (k)e−ik·r d3k. The incompressibility condition ∂Rij /∂ri = ∂Rij /∂rj = 0 implies that kiΦij = kj Φij = 0. Symmetry considerations, similar to those applied for Rij , imply that the most general form of Φij in isotropic turbulence is Φij (k) = C(k)kikj + D(k)δij and on applying the incompressibility conditions we find that k 2C(k) = −D(k). Hydrodynamics – Spring 2010 201 Thus, Φij (k) is given by one scalar function E(k) as   E(k)  ki k j  Φij (k) = δ − . ij 4πk 2 k2 Consider, then, the kinetic energy of the turbulence (per unit mass) 1 1 1 2 v = R (0) = 2 2 ii 2 ˆ 1 Φii(k) d3k = 2 ˆ   ki ki  3 E(k)  dk δ − ii 2 4πk 2 | k {z } =3−1=2 = ˆ E(k) 3 dk⇒ 4πk 2 1 2 2v = ˆ ∞ E(k) dk. 0 Thus, E(k) is the energy spectrum of turbulence, E(k) dk giving the amount of turbulent kinetic energy in the wavenumber range [k, k + dk]. Specifying E(k) is enough for a complete description of Rij in isotropic homogeneous turbulence. In the following, we will present the celebrated Kolmogorov’s universal equilibrium theory for isotropic homogeneous turbulence, which obtains an expression for E(k) from dimensional arguments. Hydrodynamics – Spring 2010 202 Turbulent fluid in equilibrium Our next task is to derive an expression for the kinetic energy distribution E(k) corresponding to a dynamic equilibrium of the fluid equations. The procedure differs from that used for finding the steadystate distribution function of the molecules, because the fluid system is dissipative. Thus, when left alone, turbulence in a fluid typically decays. We shall assume that in a steady state, a turbulent fluid is constantly stirred, i.e., kinetic energy of turbulent fluctuations are constantly fed into the system. As turbulent kinetic energy is constantly dissipated by viscous forces, this implies that in steady state the rates of energy input and dissipation agree. We shall assume also that the stirring of the fluid is performed so that the turbulence produced is homogeneous and isotropic. Kolmogorov (1941) was the first to propose a theory for such a turbulence. Turbulent velocity field can be throught of as being made of eddies of different sizes. The input energy is typically fed into the largest scales. Kolmogorov’s idea was that the largest eddies can feed energy to smaller eddies and these eddies to even smaller ones, etc., producing a cascade of energy from the largest scales to the smallest, where the dissipation occurs. The cascade may be understood as a consequence of Kelvin’s theorem. Consider a vortex tube that extends from point Q to point P . As a result of the random velocity field (generated by stirring), the distance of the points tends to get larger. Because of incompressibility, the volume of a vortex tube stays constant; thus, its diameter will tend to decrease. Thus, energy cascades to smaller scales. Hydrodynamics – Spring 2010 P P Q Q 203 Turbulent cascade Eddies of size ℓ are expected to have a certain velocity v associated with them. The corresponding Reynolds number Reℓ = vℓ/ν is expected to be large for large eddies. Thus, we may assume that energy is fed into the system at rate ǫ per unit mass per unit time at the largest eddy size, L, and velocity, V , for which VL Re = ≫1 ν Energy is then cascaded to smaller and smaller eddies. For eddies of sufficiently small size, ℓd, we expect Red ∼ 1, i.e., ℓdvd ∼ ν and the energy in these eddies is dissipated by viscosity at rate ǫ per unit mass per unit time in order to maintain equilibrium. The eddies at intermediate scales merely transmit the energy from scales L to d. The intermediate eddies are characterised by their velocity, v , and scale, ℓ. Dimensionally, the only way to express ǫ in terms of v and ℓ is v3 ǫ∼ ⇒ v ∼ (ǫℓ)1/3. ℓ This applies down to the smalles scales so ǫ∼ As ǫ ∼ V 3/L, we have Hydrodynamics – Spring 2010 vd3 ℓd = vd3ℓ3d ℓ4d ∼  1/4 3 L  L4 ǫ  ∼ 3 ℓd ν vd4 ν ⇒ ∼ ℓ4d ν   ℓd ∼  1/4 3 3 L V  = 3  ν   3 1/4 ν ǫ   = Re3/4; , vd ∼ (ǫν)1/4. V ∼ Re1/4. vd 204 log E We can now construct the approximate form of the energy spactrum E(k). As the largest scale eddies have sizes ∼ L, there is a cutoff in the spectrum at wavenumbers smaller than kL ∼ 1/L. On the other hand, the smallest eddies have sizes ∼ d so there is also a cutoff in the spectrum at wavenumbers larger than kd ∼ 1/d. The range of wavenumbers E ~ k −5/3 kL ≪ k ≪ kd ∼ Re3/4kL is called the intertial range. Within this range, the energy spectrum is expected to depend only on ǫ and k , so dimensionally E(k) = Cǫ2/3k −5/3, kL ≪ k ≪ kd kL kd log k is the only possible form of the spectrum. Here, C is a dimensionless number called the Kolmogorov constant. This is the famous “5/3-law” of Kolmogorov. We already noted that v ∼ (ǫℓ)1/3 = (ǫ/k)1/3. Thus, E(k) dk ∼ v 2 ∼ (ǫ/k)2/3 giving the same form of the spectrum once we recognize that we should take dk ∼ k in this equation. Kolmogorov spectrum was not obtained above through a rigorous analysis. It has, however, been confirmed by many observations since the 1960s. Hydrodynamics – Spring 2010 205 Rotation and hydrodynamics Hydrodynamics – Spring 2010 206 Differential rotation Start rotating a bucket filled with water. Initially, the solid wall of the rotating bucket transmits angular momentum to the outer edge of the water but the central parts of the water inside the bucket are still at rest. The water, therefore, rotates differentially, i.e., with a non-constant angular velocity. As time goes on, viscous forces tend to even out the differences in the angular velocity and after a while, the water rotates with the same angular velocity everywhere. This kind of rotation is called rigid-body (or solid-body) rotation (although the body rotating is not rigid). Differentially rotating systems in astrophysics are ubiquitous. We already met thin accretion disks, which have negligible pressure gradients and, thus, a balance between gravity and centrifugal force leading to r Ω(r) = vθ /r ≈ |gr |/r The effect of viscosity is to produce a slow inflow of material rather than to produce rigid-body rotation. Inside a slowly rotating star, gravitational force and the pressure gradient nearly balance each other. The pressure can slightly adjust itself to counter the effect of (small) centrifugal force. Viscous forces could, thus, be expected to turn the rotation finally into a solid-body rotation. Molecular viscosity, however, is completely negligible inside a star and angular momentum is transported via turbulent convection. The effective turbulent viscosity inside a star may not be isotropic, and some remaining differential rotation as a result of this anisotropy may be expected. The Sun demonstrates a difference in angular velocities between the poles and the equator by & 10 percent. Hydrodynamics – Spring 2010 207 Stability of differential rotation Consider a fluid of uniform density rotating differentially around an axis of symmetry. Suppose that a ring of fluid at distance r0 from the axis rotating at velocity v0 is interchanged with a ring at distance r1 > r0 rotating at v1. The system is stable if the interchanged rings tend to move towards their original positions. Assume that the angular momentum is conserved. Thus, rv = constant and the displaced ring originally at r0 now rotates at velocity v0(r0/r1). The ring originally at r1 had a centripetal acceleration of v12/r1, which was provided by pressure gradient forces left after balancing gravity in the equilibrium. The ring moved to r1 now requires the centripetal force [v0(r0/r1)]2/r1 = v02r02/r13 to remain at its position, and if v12/r1 is greater than this, the system is stable. Thus, r02v02 v12 < ⇔ (r02Ω0)2 < (r12Ω1)2 3 r1 r1 is the condition for stability. This reads d 2 2 [(r Ω) ] > 0 Rayleigh’s criterion. dr Hydrodynamics – Spring 2010 208 N–S in a rotating frame of reference Consider a coordinate system rotating at constant angular velocity Ω = Ωe3 in the inertial frame. This rotating frame has the Cartesian base vectors E 1 = cos Ωt e1 + sin Ωt e2; E 2 = − sin Ωt e1 + cos Ωt e2; E 3 = e3, where ei are the fixed base vectors of the inertial frame. A position vector of a point with respect to an origin at the axis of rotation is x = xiei = XiE i = X. Using dE i/dt = Ω × E i and denoting V = ẊiE i and V̇ = ẌiE i we get dx d dXi dE i = (XiE i) = E i + Xi = ViE i + XiΩ × E i = V + Ω × X dt dt dt dt   d  dx  d2Xi dXi dE i d2E i dv = + Xi 2 = V̇ + 2Ω × V + Ω × (Ω × X) Ei + 2 = f= dt dt dt dt2 dt dt dt v= Here, V and V̇ are the velocity and acceleration as measured by an observer in the rotating frame and the incompressible N–S equation in the rotating frame of reference is, therefore, ∂V 1 + V · ∇V = V̇ = − ∇p + F ext + ν∇2V − 2Ω × V − Ω × (Ω × X). ∂t ρ ∂ Note that the del operator is given by ∇ = E i ∂X . i Hydrodynamics – Spring 2010 209 Centrifugal and Coriolis forces. Effective potential. Rossby number Typically, F ext is of gravitational origin and can given by F ext = −∇Φ. Then, the centrifugal force, −Ω × (Ω × X) = 21 ∇(|Ω × X|2), can be combined with the external force introducing an effective gravitational potential , Φeff = Φ − 12 |Ω × X|2. The effect of the Coriolis force, −2Ω × V , is more complicated but it is non-zero only if there are motions with respect to the rotating frame. Thus, the pressure in a static incompressible equilibrium in the rotating frame is obtained as p/ρ + Φeff = constant. Clearly, the effect of centrifugal force needs to be taken into account when it modifies the effective potential significantly. The importance of Coriolis force is determined by the Rossby number , U U 2/L = , ǫ= ΩU ΩL which compares the two terms V · ∇V and 2Ω × V in a dynamic situation, where they do not vanish. Here, U and L are the typical velocity and length scales of the flow, respectively. At large Rossby numbers, Coriolis force is unimportant, but for slow, large-scale flows (ǫ . 1), like those of the atmosphere and oceans, it plays a crucial role in the dynamics of the fluid. In geophysical fluid dynamics, one considers the dynamics of a thin spherical shell of fluid with small Rossby number. Hydrodynamics – Spring 2010 210 Geostrophic approximation From now on, we will denote the velocity and position in the rotating frame by v and r . Fluid flows associated with large scale circulations (atmospheres, oceans) are nearly horizontal. If the flows are changing only slowly and have large Rossby numbers, dv/dt ≈ 0. The Reynolds number for a global circulation is also very large so the viscosity term can be neglected. Furthermore, if the centrifugal correction to the effective potential is small (like for the Earth), then −∇Φeff ≈ −ger and the N–S equation becomes ∇p − ger − 2Ω × v = 0. − ρ Consider the horizontal (⊥ er ) and vertical (k er ) components of the equation. Gravity dominates over the Coriolis force for typical circulation velocities, so the vertical component of the Coriolis force may be neglected. Thus, the vertical component of the force balance states that 1 ∂p = −g. ρ ∂r Since gravity is absent in the horizontal direction, Coriolis force is the one balancing pressure gradients in the horizontal spherical surface. Thus, ∇hp = −2ρ(Ω × v)h. Thus, the horizontal pressure gradient is expected to be much smaller than the vertical gradient. These equations contitute the geostrophic approximation. Note that ∇hp ⊥ v ! A familiar example of the use of this approximation is the explanation of the circulation of air around a cyclone. Hydrodynamics – Spring 2010 211 Vorticity in the rotating frame. Bjerknes’s theorem From now on, consider ideal incompressible fluids. Thus, set ν = 0 and write ρ−1∇p = ∇(p/ρ). Making use of v · ∇v = 12 ∇v 2 − v × (∇ × v), and ω ≡ ∇ × v , we can write the Euler equation in the rotating frame in form   v2 1 ∂v p  − 2Ω × v. − v × ω = −∇  + + Φ − |Ω × v|2 ∂t ρ 2 2 On taking the curl of this equation, we again get rid of the “potential part” of the force, but the Coriolis term does not vanish. Thus, in a rotating coordinate system, ∂ω = ∇ × (v × ω − 2Ω × v) = ∇ × [v × (ω + 2Ω)]. ∂t Since Ω does not depend on time, this can be written as ∂Q = ∇ × [v × Q] ∂t with Q = ω + 2Ω. Thus, we can immediately conclude that Kelvin’s theorem on vorticity is generalized to ˆ d (ω + 2Ω) · ds = 0 Bjerknes’s theorem. dt S Hydrodynamics – Spring 2010 212 Taylor–Proudman theorem. Taylor column Consider a steady flow in the rotating frame of reference. Thus, ∇ × [v × (ω + 2Ω)] = 0. If the flow is slow (or the system is large) so that ω ≪ Ω, i.e., the Rossby number is small, then 0 = ∇ × (v × Ω) = Ω · ∇v − Ω(∇ · v) − v · ∇Ω + v(∇ · Ω) = Ω · ∇v for incompressible flow. This means that v does not change in the direction of Ω. Thus, slow steady flows are invariant in the direction of Ω, which is known as the Taylor–Proudman theorem. If a flow meets an obstacle in a rapidly rotating frame, a Taylor column is generated above it. The flow does not only flow around the obstacle but also around its extension in the direction of Ω. Hydrodynamics – Spring 2010 213 Self-gravitating rotating masses Consider a rigidly rotating fluid of constant density under self-gravity. It may be expected that the shape of the body differs from a sphere being flatter at the poles. In the rotating frame   p ∇  + Φ − 12 |Ω × r|2 = 0 ρ where Φ is the gravitational potential. This implies that p + Φ − 12 |Ω × r|2 = constant. ρ On the outer surface of the body, pressure is contant. Choose the z -axis co-aligned with the axis of rotation. Thus, Φ(x, y, z) − 21 Ω2(x2 + y 2) = constant on the outer surface. Knowing Φ would give us an equation for the outer surface. We want to establish that the shape of the body is an ellipsoid. For this purpose, we need the gravitational potential of an ellipsoid. For an ellipsoid of uniform density inside the bounding surface x2 y 2 z 2 + + =1 a2 b 2 c 2 the gravitational potential inside the surface is given by Φ = πGρ(α0x2 + β0y 2 + γ0z 2 − χ0), Hydrodynamics – Spring 2010 214 where α0 = abc ˆ ∞ dλ ; (a2 + λ)∆ γ0 = abc ˆ ∞ dλ ; (c2 + λ)∆ 0 0 β0 = abc ˆ ∞ dλ (b2 + λ)∆ χ0 = abc ˆ ∞ dλ , ∆ 0 0 r and ∆ = (a2 + λ)(b2 + λ)(c2 + λ). If the rotating mass of fluid takes the shape of an ellipsoid, then Φ = πGρ(α0x2 + β0y 2 + γ0z 2 − χ0), must satisfy Φ − 12 Ω2(x2 + y 2) = constant ⇒ (2πGρα0 − Ω2)x2 + (2πGρβ0 − Ω2)y 2 + 2πGργ0z 2 = constant on the surface given by Thus, x2 y 2 z 2 + + = 1. a2 b 2 c 2 a−2 ÷ b−2 ÷ c−2 = (2πGρα0 − Ω2) ÷ (2πGρβ0 − Ω2) ÷ 2πGργ0. If we can find solutions of this equation, then we would have established that rotating fluids take up ellipsoidal configurations. Hydrodynamics – Spring 2010 215 Maclaurin spheroids Intuitively, one expects a the rtating fluid to be flattened at the poles and symmetric with respect to the rotation axis. Thus, one adopts a = b > c, i.e., a spheroidal configuration with eccentricity c2 e = 1 − 2. a Such solutions of the problem are called Maclaurin spheroids in honor of Maclaurin, who demonstrated these solutions in 1740. For a = b, the integrals can be solved as √ 1 − e2 1 − e2 −1 α0 = β0 = sin e − e3 e2   −1 √ sin e 2 1 −  1 − e2 2  γ0 = 2  e e and substituting these to gives Hydrodynamics – Spring 2010 a−2 2πGρα0 − Ω2 1 − e = −2 = c 2πGργ0 √ 6 2 1 − e2 Ω2 2 −1 2 (3 − 2e ) sin e − (1 − e ). = πGρ e3 e2 216 (sketch) 0.5 (sketch) 1.0 Ω2 0.4 πGρ 0.3 L 0.8 0.6 0.2 0.4 0.1 0.2 0.2 0.4 0.6 0.8 1.0 e 0.2 0.4 0.6 0.8 1.0 e Ω2/πGρ reaches a maximum value of 0.499 at e ≈ 0.930. Thus, after the maximum, less rotation makes the system more flattened. But if, instead of√Ω, one looks at the angular momentum L = 52 M a4Ω of the spheroid of mass M = 43rπa2cρ = 43 πa3 1 − e2ρ, the situation looks different. L can be made dimensionless dividing it by GM 3(a2c)1/3, which is constant for all spheroids (of different e) with constant mass. Thus, L̄ = L r GM 3(a2c)1/3 √ 3 a = 5 c v !2/3 u u u u t √ v u u u u t Ω2 3 Ω2 −1/3 (1 − e) = πGρ 5 πGρ is a measure of the angular momentum of the body, and that increases without a limit as e → 1. The increase, of course, comes from the increase of the moment of inertia of the body. Hydrodynamics – Spring 2010 217 Jacobi ellipsoids Jacobi (1834) established that also ellipsoidal solution with three unequal axes are admitted under some circumstances. Such solutions are called Jacobi ellipsoids. It can be shown that for L̄ < 0.304 (i.e., e < 0.813), Maclaurin spheroids are the only possible ellipsoidal solutions. However, at L̄ > 0.304 a Jacobi ellipsoid with three unequal axes becomes possible. When Jacobi ellipsoids are admitted, we must ask which of the two possible solutions is more probable in Nature. For a given angular momentum, it can be shown that a Jacobi ellipsoid has less rotational energy than a Maclaurin spheroid. Thus, an internal dissipation mechanism (like viscosity) may make a Maclaurin spheroid unstable and relax to a Jacobi ellipsoid after some dissipation of rotational energy. Thus, the solution has spontaneously broken the symmetry which was present in the formulation of the problem. Are there any astronomical objects that would be examples of Jacobi ellipsoids? Some elliptic galaxies show tri-axial shape, so one might be tempted to call these stellar-dynamics analogues of Jacobi ellipsoids. Elliptical galaxies, however, have very little angular momentum, so the explanation is probably incorrect. Rather, the shape stems from the anisotropic velocity distribution of the stars inside the galaxy. Hydrodynamics – Spring 2010 218 Magnetohydrodynamics (MHD) Hydrodynamics – Spring 2010 219 Scope of magnetohydrodynamics So far we have considered neutral fluids, in which collisions are the main means of transporting momentum and energy between the molecules of the gas. The only force acting at large distance in neutral fluids is gravity. We introduced some fluid effects, where gravitational field acted as the restoring force, mainly involving surface phenomena (recall, however, also the internal gravity waves). For fluids made of charged particles, the situation is completely different. Besides gravitational fields, also the electromagnetic fields now affect the motion of the molecules. This inevitably leads to a more complicated dynamical theory than ordinary hydrodynamics, since fluids consisting of charged particles not only react to external electromagnetic forces but act as sources of the fields as well. Therefore, our former usual convention of treating body forces as external forces no longer apply. Magnetohydrodynamics is primarily a macroscopic dynamical theory of plasmas, which are quasineutral7 gases containing enough free charges to make collective electromagnetic phenomena important for their physical behaviour. The latter means that plasmas are dense enough so that a large number of particles are interacting with each other through Coulomb interactions at the same time. For a more quantitative definion of a plasma state, see the lectures of “Plasma physics”. For our purposes it is enough to understand that overwhelmingly largest part of the ordinary matter in the unverse is in the plasma state. (All stars are made of plasma; inside the orbit of Pluto, 0.1% of mass (mainly in Jupiter) and 10−15% of volume is in non-plasma state.) Besides plasmas MHD can also be applied to electrically conducting liquids like mercury. 7 Quasineutral ionized gases consist of equal amounts of positive and negative charges. Thus, the net charge of all the particles of a quasineutral gas is zero. Hydrodynamics – Spring 2010 220 Equation of continuity and momentum equation Equations governing the motion of the fluid in MHD are very similar to the HD equations. We will justify them here very heuristically and refer the student to the lectures of “Plasma physics” for a more rigorous derivation. The equation continuity (conservation of mass) ∂ρ + ∇ · (ρv) = 0 ∂t of course applies in MHD as well. The momentum equation has to be supplemented with forces that come from the interactions of the molecules with the electric (E ) and magnetic (B ) fields permeating the fluid. A charged particle (charge q ) feels the Lorentz force, q(E + v × B) so a large number N of charged particles (indexed by i) inside a volume ∆V feels the resultant force per unit volume of N 1 X qi(E + v i × B). f= ∆V i=1 P In a quasineutral situation, the sum i qi = 0, so the first term vanishes. The second term corresponds to f = j × B, where j = (∆V )−1 i qi v i P Hydrodynamics – Spring 2010 is the electric current density (SI unit A m−2) carried by the fluid. 221 Thus, the momentum equation of MHD reads ρ dv = ρF − ∇p + j × B + ρν∇2v, dt where all other body forces than the Lorentz force (e.g., gravity) have been included in F . The next task is to find out, how the current density and the magnetic field evolve. In MHD, one limits the theory to describe only “slow” phenomena, i.e., one takes Ampère’s law to apply. Thus, µ0j = ∇ × B, where µ0 is the permeability of vacuum. (In plasmas, all charges can be regarded as free.) Thus, the current density is immediately eliminated from the equation of motion and one gets ρ dv 1 = ρF − ∇p + (∇ × B) × B + ρν∇2v. dt µ0 Before turning to the evolution of the magnetic field, let us examine the Lorentz force somewhat more closely. First, use the vector identity   2 B  (∇ × B) × B = (B · ∇)B − ∇   . 2 Comparing the second term with the pressure gradient shows that the magnetic field introduces a pressure of B 2/2µ0. The first term represents a tension along the magnetic field lines. Hydrodynamics – Spring 2010 222 Magnetic force in tensorial form In fact, the magnetic body force can be written in a form of a surface force with the aid of the identity ∇ · (BB) = (∇ · B)B + (B · ∇)B. As magnetic field is always source free (∇ · B = 0), we get for the Lorentz force   2 BB  1 B   = −∇ · M. (∇ × B) × B = −∇ ·  I− µ0 2µ0 µ0 Thus, ∂ dvi = ρFi − (Pij + Mij ), dt ∂xj where Mij = (B 2/2µ0)δij − BiBj /µ0. The second term is a tension along the field lines, which can be seen by writing the tensor in matrix form in a (local) coordinate system with ez k B : ρ (Mij ) =        2 B /2µ0 0 0 0 B 2/2µ0 0 0 0 2 B /2µ0        −        0 0 0 0 0 0 0 0 2 B /µ0        The first term corresponds to an isotropic pressure (generating expansion) and the second term to a tension along the magnetic field (trying to keep the lines of force as short as possible). Hydrodynamics – Spring 2010 223 Faraday’s and Ohm’s laws. Induction equation The evolution of the magnetic field is governed by Faraday’s law ∂B = −∇ × E. ∂t Thus, although the electric field drops out from the equation of motion, it re-enters the problem through Faraday’s law. Electric field is related to the current density via the Ohm’s law. In the plasma rest frame, we write E ′ = ηj = η ∇ × B, µ0 where η is the resistivity of the fluid. The electric field transforms from the plasma rest frame to the laboratory frame as E = E′ − v × B giving η ∂B = ∇ × (v × B) + ∇2B, ∂t µ0 where resistivity has been assumed constant and the vector identity ∇ × (∇ × B) = ∇(∇ · B) − ∇2B along with ∇ · B = 0 has been used. This is the induction equation of MHD. Hydrodynamics – Spring 2010 224 MHD equations. Energy equation The MHD system of equations is ∂ρ + ∇ · (ρv) = 0 ∂t ∇p (∇ × B) × B ∂v + v · ∇v = F − + + ν∇2v ∂t ρ µ0 ρ ∂B η = ∇ × (v × B) + ∇2B ∂t µ0 appended with an appropriate equation of state or the energy equation. The MHD energy equation is   ∂ǫ ρ  + v · ∇ǫ = −p∇ · v + ∇ · (K∇T ) + ηj 2, ∂t where the last term is due to Ohmic heating of the fluid, i.e., due to dissipation of the electromagnetic fields. For highly conducting plasmas, the last two terms are often neglected, and the internal energy written as ǫ = cV T = p/ρ(γ − 1). Then, as for neutral fluids, one obtains   d p = 0. dt ργ Instead of an energy equation, it is also often customary to use a barotropic equation of state or simply assume the fluid to be incompressible. In those cases, of course, the energy equation becomes redundant. Hydrodynamics – Spring 2010 225 Magnetic Reynolds number Consider the MHD induction equation ∂B = ∇ × (v × B) + λ∇2B, ∂t where λ = η/µ0 is the magnetic diffusivity. Taking the ratio of the two terms on the right-hand side, i.e., V B/L LV = ReM = λB/L2 λ defines a dimensionless number known as the magnetic Reynolds number . The direct proportionality of ReM on L means that it is much larger for astrophysical plasmas compared to laboratory plasmas. For example, for a hydrogen plasma at T = 104 K the magnetic diffusivity is λ ≈ 103 m2 s−1. Thus, for a typical laboratory system with L ∼ 1 m and V ∼ 0.1 m s−1, ReM ∼ 10−4. Instead, for example, for a convection cell at the solar surface L ∼ 106 m and V ∼ 103 m s−1 we find ReM ∼ 106. Hydrodynamics – Spring 2010 226 Ideal MHD. Alfvén’s theorem on flux-freezing Thus, under typical conditions: Laboratory: Astrophysics: ∂B ≈ λ∇2B ∂t ∂B ≈ ∇ × (v × B) ∂t Of course, these are very crude simplifications: in many laboratory conditions ∇ × (v × B) cannot be neglected and in astrophysical plasmas λ∇2B is important in localized regions with large currents. The limit of λ → 0, i.e., ∂B = ∇ × (v × B) ∂t is known as the ideal MHD limit. Our general proof of Kelvin’s theorem on vorticity again becomes handy, because it shows immediately (since ∇ · B = 0) that ˆ d B · ds = 0 dt S for any material surface S following the fluid flow. This result was first pointed out by Alfvén (1942) and it is called Alfvén’s theorem on flux-freezing . It states that the flux inside a magnetic flux tube remains constant in an ideal MHD flow. Alfvén’s theorem is the magnetic analogue of Kelvin’s theorem on vorticity. Hydrodynamics – Spring 2010 227 Consequences of Alfvén’s theorem In astrophysical systems with high magnetic Reynolds number we can then imagine the magnetic field lines to be frozen-in to the plasma flow. Suppose we have a straight plasma column with homogeneous magnetic field. If the plasma column is bent or twisted, then the magnetic field lines become bent or twisted as well. In general, if two plasma elements are at some time connected with a magnetic field line, they remain forever connected with a magnetic field line in ideal MHD. Hydrodynamics – Spring 2010 228 Magnetic field amplification/generation by plasma motions Consider, first, the evolution of a magnetic field in a collapsing magnetized plasma cloud, i.e., a star forming from a cloud of interstellar plasma. As magnetic flux is frozen to the plasma, the magnetic field increases in inverse proportion to the cross-sectional area of the plasma cloud. Starting from typical values of interstellar density and magnetic field, the solar magnetic field would be of the order of ∼ 106 T in case no mechanism would be available for annihilating magnetic energy (actually, B ∼ 1 mT). Solar magnetic field, however, follows a 22-year activity quasi-cycle involving polarity changes. This is thought to occur because of the continuous re-generation of the field through a dynamo mechanism involving differential rotation and convection inside the Sun. Hydrodynamics – Spring 2010 229 MHD waves MHD contains a wealth of wave phenomena that modified with respect to or not observed at all in neutral fluids. By linearizing the MHD equations around a homogeneous static equilibrium with ρ = ρ0 +ρ1(x, t), v = v 1(x, t), p = p0 + p1(x, t) and B = B 0 + B 1(x, t), we get ∂ρ1 + ρ0 ∇ · v 1 = 0 ∂t ∂v 1 ∇p1 (∇ × B 1) × B0 =− + ∂t ρ0 µ 0 ρ0 ∇p1 = c2s ∇ρ1 ∂B 1 = ∇ × (v 1 × B 0) ∂t with c2s = γp0/ρ0 for an adiabatic fluid with zero viscosity and resistivity. Taking the time derivative of the equation of motion and using the other three equations to eliminate ρ1, p1 and B 1 we get the MHD wave equation, ∂ 2v1 = c2s ∇(∇ · v 1) + {∇ × [∇ × (v 1 × v A)]} × v A, 2 ∂t ⇒ B0 is the Alfvén velocity where v A = √ µ 0 ρ0 ω 2v 1 = c2s k(k · v 1) − {k × [k × (v 1 × v A)]} × v A Hydrodynamics – Spring 2010 MHD dispersion relation. 230 Applying the vector identity a × (b × c) = b(a · c) − c(a · b) (i.e., the “BACCAB formula”) three times, the dispersion relation can be written in form 2 0 = −ω 2v 1 + (c2s + vA )k(k · v 1) + (v A · k)[(v A · k)v 1 − k(v A · v 1) − v A(k · v 1)] 2 )kk − (v A · k)[kv A + v Ak]} · v 1 ≡ D · v 1. = {[(v A · k)2 − ω 2]I + (c2s + vA The non-trivial (v 1 6= 0) solutions of this equation are obtained as det D = 0, 2 where Dij = [(v A ·k)2−ω 2]δij +(c2s +vA )kikj −(v A ·k)(kivAj +vAikj ) is the (symmetric) dispersion tensor of the MHD fluid. Choosing the coordinate system so that v A = vAez and k = k(sin θex + cos θez ) we can write 0 = D · v1 =        2 −ω + 2 k 2 vA k 2c2s sin2 θ + 0 c2s k 2 sin θ cos θ 0 2 −ω + 2 2 vA k cos2 θ 0 c2s k 2 sin θ cos θ 0 −ω 2 + k 2c2s cos2 θ        v1x v1y v1z        Setting the determinant of this matrix to zero gives three solutions, i.e., three different MHD wave modes v 1 k ey ; 2 ω 2/k 2 = vA cos2 θ shear Alfvén wave v 1 ⊥ ey ; 2 2 2 + c2s ) + [ 41 (vA + c2s )2 − c2s vA cos2 θ]1/2 fast MHD wave ω 2/k 2 = 12 (vA v 1 ⊥ ey ; 2 2 2 + c2s ) − [ 41 (vA + c2s )2 − c2s vA ω 2/k 2 = 12 (vA cos2 θ]1/2 slow MHD wave Hydrodynamics – Spring 2010 231 Propagation parallel to the field Let us, first, investigate the MHD wave modes at the limit θ → 0, i.e., when the waves are propagating parallel to the field lines. Longitudinal waves (v 1 k k k v A) are simple sound waves with ω = ±csk. Either the fast (if cs > vA) or the slow (if vA > cs) MHD mode reduces to a sound wave at θ → 0. The dispersion relation of the shear Alfvén wave becomes ω = ±vAk, v1 ⊥ k which shows that these waves are transverse waves. The magnetic field of the wave is obtained as −iωB 1 = ik × (v 1 × B 0) = ikB0v 1 ⇒ B1 = − B0 v1 ⊥ B 0 vA showing that the magnetic field fluctuations are also transverse and bend the magnetic field lines. As vA = v u u u u t B 2/µ0 ρ0 = v u u u t magnetic tension , mass density shear Alfvén waves are a direct analogue of transverse waves propagating along a string. Either the fast (if vA > cs) or the slow (if vA < cs) MHD mode reduces to a shear-wave at θ → 0. Hydrodynamics – Spring 2010 232 Propagation perpendicular to the field If θ → π/2, then both the shear Alfén wave and the slow MHD wave have zero phase speed, i.e., these waves do not propagate perpendicular to the magnetic field. The dispersion relation of the fast mode is r 2 + c2 k, ω = ± vA s r 2 + c2 is called the magnetosonic speed. Correspondingly, this wave is referred to as where the speed vA s the magnetosonic wave. Magnetosonic waves compress the magnetic field, i.e., B1 = v1 B 0. vA The restoring force in this wave is the gradient of total pressure. This can be seen by writing 2 vA + c2s γp0 γB pB γp0 B02 + = + , = µ 0 ρ0 ρ0 ρ0 ρ0 where pB = B02/2µ0 is the magnetic pressure and γB = 2 is the adiabatic index corresponding to a twodimensional gas. This is natural, since the magnetic force j × B 0 is always in the plane perpendicular to the magnetic field (i.e., two-dimensional) and any motions of the fluid in the direction of the magnetic field leave the field unchanged. Hydrodynamics – Spring 2010 233 Magnetohydrostatics Static solutions of MHD, in the strict sense, would need to fulfill λ∇2B = ∂B/∂t = 0 in addition to fulfilling the force balance equation ∇p = ρF + j × B. However, it is customary to consider any configuration that fulfills the force balance equation a magnetohydrostatic equlibrium. This, of course, assumes implicitly that the transport of momentum via MHD waves has a much shorter times scale than the decay of the magnetic field via diffusion. In the absence of the external forces, the force balance equation reduces to ∇p = (∇ × B) × B . µ0 In astrophysics, we often deal with plasmas, for which the plasma beta, β= p 2µ0p = pB B2 is small. Thus, the force balance equation reduces to (∇ × B) × B = 0 ⇒ ∇ × B = α(x)B, where α(x) is a scalar. Magnetic fields fulfilling this equation are called force-free magnetic fields. Hydrodynamics – Spring 2010 234 As ∇ · (∇ × B) = 0 and ∇ · B = 0 identically, the force-free equation implies B · ∇α = 0, i.e., that α is constant along magnetic field lines. If α is a constant over all space, the equation ∇ × B = αB is linear. Linear force-free fields form an important category of MHD equilibrium solutions in solar and heliospheric physics, for example. The force-free magnetic field solutions include also the potential fields, for which α = 0. We note that a force-free field may be a valid solution of the magnetostatic problem even if the plasma beta is not small, if there is an external force field (like gravitational field) present in the system. In that case, the balance is obtained so that the hydrostatic balance and the force-free condition are satisfied simultaneously: ∇p = ρF 0 = (∇ × B) × B Finally, we note that in a general equilibrium in the absence of external forces, the pressure gradient has to be balanced by the magnetic force. Such solutions are considered, for example, when studying the magnetic confinement of plasma, which is important for fusion research. In this respect, also the analysis of the stability of the MHD equilibria is extremely important, and it turns out that many equilibrium solutions of plasmas confined by magnetic fields are actually unstable. More on MHD on a special course next fall. Welcome! Hydrodynamics – Spring 2010 235