Physics 526 Notes: Fluid Dynamics: Jean Eilek
Physics 526 Notes: Fluid Dynamics: Jean Eilek
Physics 526 Notes: Fluid Dynamics: Jean Eilek
Jean Eilek
Physics Department, New Mexico Tech
Socorro, NM 87801, U.S.A. and
[email protected]
Our goal in this course is to explore the basic D. Work in a Rotating Frame 3
physics of fluid systems, both with and without the ef- E. Dimensional Analysis 4
fects of magnetic fields. Fluid physics by itself – hydro- F. Appendix: When can we use hydrodynamics? 5
dynamics – applies to any neutral fluid (think about wa- 1. hard sphere collisions 5
ter, or molasses, or the earth’s atmosphere), whether or 2. plasmas: the coulomb cross section 6
not there are magnetic fields around. “Hydro” or “HD” 3. collisionless plasmas 7
has broad applications – smooth flows, turbulent flows,
shocks, sound waves, instabilities. If the fluid is ionized 2. VISCOSITY AND LAMINAR FLOW: BASICS 8
(think of liquid sodium, the higher reaches of the earth’s A. One-dimensional Laminar Flows 8
atmosphere, or just about any astrophysical system you B. Steady Flow: Cartesian Applications 9
can bring to mind), it can (and almost certainly will) 1. flow between parallel plates 9
carry a current. The interaction of the current with a 2. flow in an open channel 9
magnetic field (external or self-generated) modifies all 3. hele shaw flow 9
of the phenomena above and adds some new ones. This C. Steady Flow: Cylindrical Applications 10
is described by magnetohydrodynamics (MHD). 1. pipe flow 10
Traditionally, HD and MHD are treated separately, 2. circulating flow 10
but they have a lot in common, and I have learned a lot D. Viscous Stresses, Generally 11
by comparing and contrasting HD and MHD phenom- 1. do it physically first 11
ena. I’m hoping to share that with you as we go through 2. then do it formally 12
this course. E. The Navier-Stokes Equation (in Cartesian) 13
A word of caution: you should especially note units F. Appendix: Navier-Stokes in other coordinates 13
and dimensions. These notes are in cgs. That makes 1. Cylindrical coordinates 13
very little difference for “rocks” (analyses that involve 2. Spherical Polar Coordinates 14
mass, length, time); but it makes a big difference for
electrodynamics and MHD. The E and B fields, as well 3. POTENTIAL FLOW 16
as the fundamental charge, have different dimensions in A. Setups and Assumptions 16
cgs than in SI; and the coupling constants in Maxwell’s 1. velocity potential 16
equations are different. Both systems appear in the lit- 2. stream function 16
erature; while there is some trend towards SI, many im- 3. other coordinates 16
portant references are in cgs. I am going to follow my B. Two-dimensional Planar Problems 17
experience and preference and use cgs; I’ll also give the 1. sources and sinks, flow past half body 17
SI versions of most critical equations. 2. flow past a cylinder 18
C. Axisymmetric 3D Problems 19
point source and stream flow 19
Contents flow around a sphere 19
a line source 19
1. BASIC CONCEPTS AND TOOLS 1 D. d’Alembert’s “Paradox” 20
A. Kinematics: How to Describe a Flow 1
stream function 1 4. VORTICITY 22
B. Mass Conservation: The Continuity Equation 1 A. Vortex Kinematics 22
lagrangian derivative 2 B. Vortex Dynamics 22
C. Momentum Conservation: Euler’s Equation 2 C. Conservation of Circulation: Kelvin’s Theorem 23
1. control volumes 3 D. The Magnus Effect 24
2. bernoulli’s relation 3 E. Vortex lines and their behavior 25
ii
A; let the surface haveRan outward normal, n̂. The mass In the second expression, we have used Gauss’s law for
within this volume is V ρdV , if ρ is the mass density. tensors (noting that ρvv is a second-rank tensor).
The net rate of change of this mass is Now, the net rate of change of momentum in the vol-
ume must be equal to the net force exerted on the vol-
d
Z
ρdV ; ume. We consider external forces which act throughout
dt V the volume (“body” forces, such as gravity, electomag-
netism, bouyancy, radiation pressure if the fluid is opti-
if there are no sources or sinks of matter, this quantity cally thin; we let g be the net force per mass), and also
must equal zero. Now, there are two ways this integral the force exerted on the surface by the fluid outside V .
can change with time. (i) there can be intrinsic variation The net force on the volume V is, then,
of ρ, ∂ρ/∂t 6= 0; or (ii) there can be flow into or out of Z Z Z Z
the volume, at a rate ρv · n̂ per surface area. The sum ρgdV − pn̂dA = ρgdV − ∇pdV (1.8)
of (i) and (ii) must balance out to zero: V A V V
d
Z Z
∂ρ
Z where we have again used vector identities in the last
ρdV = dV + ρv · n̂dA = 0 (1.3) step. This balance, (1.7) and (1.8), can then be written
dt V V ∂t A
in differential form,
R
But
R the surface integral can be written as A ρv· n̂dA = ∂
(ρv) + ∇ · (ρvv) = −∇p + ρg (1.9)
V ∇ · (ρv)dV . Since V is arbitrary, we can set the ∂t
integrand to zero, and we get the differential form of This is our basic force equation, known as Euler’s equa-
this basic equation: tion. Note that the vv term is a tensor.1
∂ρ It is conventional to simplify this, by expanding the
+ ∇ · (ρv) = 0 (1.4) derivatives on the left hand side and using (1.4); we get
∂t
∂v
This is, of course, the continuity equation, applied to ρ + ρ(v · ∇)v = −∇p + ρg (1.10)
mass conservation. ∂t
Another version of Euler’s equation holds if there is no
LAGRANGIAN DERIVATIVE body force. To get to this, write the pressure in terms of
the pressure tensor (for an isotropic gas pressure):
Note, the two terms in the above expression describe
the two “intrinsic” ways in which the mass in the ele-
p 0 0
mental volume can change. We collect them as
P = 0 p 0 = pδij (1.11)
D ∂
= +v·∇ (1.5) 0 0 p
Dt ∂t
and with this, the continuity equation becomes and note that the divergance of a tensor is found (in
Cartesian) by
Dρ
+ ρ∇ · v = 0 (1.6) ∂Pij
Dt (∇ · P)i = (1.12)
∂xj
C. Momentum Conservation: Euler’s Equation with the Einstein summation convention assumed.
Consider again our surface A, enclosing With this notation, the Euler equation in the case of
R volume V .
The momentum within this surface is V ρvdV . The no body force can be written,
net rate of change of this quantity again must reflect in- ∂
trinsic (∂/∂t 6= 0) variation and advection (flow across (ρv) + ∇ · (ρvv + P) = 0 (1.13)
∂t
the surface). Thus, we write the net rate of change of
momentum as
1
In general, a “vector product” tensor ab expands out as
∂
Z Z
(ρv)dV + (ρv)v · n̂dA
2 3
a 1 b 1 a1 b 2 a 1 b 3
V ∂t A ab = 4a2 b1 a2 b2 a2 b3 5
(1.7)
∂
Z Z
a 3 b 1 a3 b 2 a 3 b 3
= (ρv)dV + ∇ · (ρvv)dV
V ∂t V
3
This, and (1.4), are conservative forms of the basic But now: the right hand side of (1.16) is normal to both
equations. the local flow field (that is normal to streamlines) and
to the local vorticity ω~ . Thus, we have one form of
1. CONTROL VOLUMES Bernoulli’s relation: in inviscid, steady flow, the term
We argued above that equating (1.7) and (1.8) will give in brackets has zero gradient in the direction of the local
an integral form of the momentum balance equation. velocity field. Thus, we have one version of Bernoulli’s
This is true if the volume V and surface A are fixed (not law:
1 2 dp
Z
themselves moving). It can also be useful to consider v + + Φg = constant along streamline
a moving volume (such as in a rocket problem). To 2 ρ
specify, let V ∗ (t) be the instantaneous control volume, (1.17)
and A∗ (t) be its area. Let b be the local velocity of the Further, in an adiabatic gas, p ∝ ργ if γ is the adi-
boundary. Our integral equations become atabic index (the ratio of specific heats). The second
term simplifies, so that Bernoulli’s relation for an invis-
d
Z Z
ρdV + ρ(v − b) · n̂dA = 0 cid adiabatic gas is
dt V ∗ (t) A∗ (t) 1 2 γ p
v + + Φg = constant along streamline
and 2 γ−1ρ
(1.18)
d
Z Z
ρvdV + ρv(v − b) · n̂dA Alternatively, in an incompressible fluid, ρ is constant,
dt V ∗ (t) A∗ (t) and the second term in (1.16) becomes simply p/ρ.
Thus, for an incompressible fluid, Bernoulli’s relation
Z Z
= ρgdV − pn̂dA is
V ∗ (t) A∗ (t)
1 2 p
Exercise for the student: can you derive or justify these v + + Φg = constant along streamline (1.19)
2 ρ
relations?
D. Work in a Rotating Frame
2. BERNOULLI ’ S RELATION
Euler’s equation effectively expresses force balance. It
It might be comforting to prove that we can extract must therefore be modified in a rotating system, where
Bernoulli’s relationship from what we have so far. we expect Coriolis and centrifugal forces. Recall your
Start with Euler’s equation, in the form (1.10). But basic mechanics. We want to transform vectors from
now, note two useful facts. The first is that if the fluid is our (inertial) frame to a frame rotating at Ω . Let r be
barotropic – that is if pressure is a function of density the usual vector from the coordinate origin to the obser-
only and vice-versa (as in an adiabatic gas), we have vation point (polar coordinates), and let R be a vector
from the rotation axis, perpendicular, to the observation
1 dp
Z
∇p = ∇ (1.14) point (cylindrical coordinates). The time derivative of a
ρ ρ general vector, P, transforms as
(this can be verified using the chain rule; take p =
µ ¶ µ ¶
dP dP
F (ρ), F being some function, and go from there). = +Ω ×P
dt i dt r
Thus, this term is a perfect differential. The second use-
ful fact is that Thus, velocities and accelerations transform as
µ ¶
1 2 vi = vr + Ω × r ;
v · ∇v = −v × ω + ∇ v (1.15)
2
Ω × vr + Ω × (Ω
ai = ar + 2Ω Ω × r)
(this is easiest to verify by expanding out in Cartesian
coordinates). Thus, this term is also a perfect differen- Ω × vr − Ω 2 R
= ar + 2Ω
tial. The first term on the right hand side is written in
terms of ω = ∇ × v, the local vorticity (more in chap- From this, we can write the incompressible Euler equa-
ter 4 on this). Specify the force to gravity, which can be tion in a rotating frame (now dropping the r subscripts):
expressed in terms of a potential: g = ∇Φg . If we then Dv 1
consider steady flow, we can rewrite (1.10) as = − ∇p + g + Ω2 R − 2Ω Ω×v (1.20)
Dt ρ
· ¸
1 2 dp
Z
The two new terms are the centrifugal and Coriolis
∇ v + + Φg = v × ω (1.16)
2 ρ forces.
4
VL
Re = (1.24) 2
The SI expression is β = 2µo p/B 2 .
ν
5
2. PLASMAS : THE COULOMB CROSS SECTION (clearly, we cannnot integrate from bmin = 0 to bmax =
∞, since the integral in (1.38) would diverge). bmin is
There is another type of encounter which is important in usually taken to be the distance corresponding to maxi-
plasmas: a long-range encounter between two particles mum energy transfer,
which feel a 1/r2 Coulomb force. The particles never
directly collide with each other; but the long-range scat- e2
tering allows exchange of energy and momentum, and bmin ≃ (1.41)
2me v 2
thus makes the system act like a fluid.
• Start with a single encounter, in which particle A bmax is less straightforward. A common choice is the
(an electron, say) scatters on particle B (a proton, say; Debye shielding length (the scale over which an extra
with mp >> me , we can assume the proton stays at charge causes charge separation in a plasma):
rest. Let the incoming particle have velocity v and mass ¢1/2
bmax ≃ λD = kB T /4πne2
¡
me , and let in come in at impact parameter b. We can (1.42)
solve this problem exactly, from classical mechanics, (kB is the Boltzmann constant, and T is the temper-
and find the deflection angle, θ, and the resultant veloc- ature). Thus, the best choice of ln Λ clearly depends
ity and momentum changes, ∆p = m∆v. Here, we on the exact situation one is considering. Luckily, for
will approximate this analysis. our purposes, this is only a logarithmic uncertainty, and
R • The net impulse on the electron will be ∆p = will not be critical for most of our calculations. The
F(t)dt, integrated over the collision. Now, the force choices above, with typical astrophysical parameters,
is strong only when the two particles are close. Since give ln Λ ≃ 10 − 20, in almost any diffuse-matter set-
they are close for a period of time ∆t ≃ 2b/v, we can ting.
approximate F ≃ e2 /b2 and ∆p ≃ 2F b/v. (Since we • Numerically, for a thermal plasma with 12 me v 2 ≃
know the net deflection is perpendicular to the initial kB T , the Coulomb cross section becomes,
direction of motion, we can also drop the vector nota-
tion). This gives us the net energy gain per collision, ln Λ
σc ≃ 7 × 10−13 2 cm2 (1.43)
T4
(∆p)2 2e4
∆E = ≃
2me me b2 v 2 where T4 = T /104 K; so that
3. COLLISIONLESS PLASMAS
2. VISCOSITY AND LAMINAR FLOW: BA- Beware: you must be careful when working with
SICS viscosity. The definitions, symbols and even units
of viscosity vary from author to author. I pre-
Let’s start by looking at some simple fluid flow so- fer the choice, that ν∇v has the dimensions of a
lutions. We’ll consider steady-state, incompressible force per area per gram: and that ρν∇v is a force
flows. Some very simple solutions can be found when per area. Thus, the dimensions of ν are [L]2 /[T].
viscosity is important to the flow. Because viscosous Many authors (including Faber) include the den-
stresses can, however, get very complicated, we’ll be- sity; he uses η (dimensions [M]/[L][T]) which is
gin with simple examples, then see the general stuff. the same as our ρν. Still other authors use µ where
Faber uses η.
A. One-dimensional Laminar Flows So: putting this back into vector form, we’ve discov-
ered that our governing equation for this simple system
Think about a long, 1D channel such as in Figure 2.1. should be1
The flow is driven by a pressure gradient, p2 < p1 ; and
we assume the velocity must go to zero at the walls (this ∇p = ρν∇2 v
is a no-slip boundary). Because the flow is incompress-
ible – has the same density everywhere – the velocity But now ... looking back to equation (1.10), our re-
profile must be independent of x, in order to enable sult here suggests that we should add a new term to de-
mass conservation (the same flow rate across any point scribe the effects of viscous stresses. Taking this as true
within the channel). for the moment (we’ll derive it more generally below),
we have the basic equation governing viscous flow:
y
dx v
(p1) dy (p2 ) ∂v
+ ρ(v · ∇)v + ∇p = ρg + ρν∇2 v
2b
x
ρ (2.2)
∂t
Figure 2.1. Defining the geometry for this problem. The
flow is driven by a pressure gradient; p2 < p1 . This driver,
Compare this back to our first version of the momentum
combined with the no-slip condition assumed at the walls equation, namely (1.10), which we called Euler’s equa-
and the action of viscous stresses within the flow, leads to tion. This version of the force equation, with viscosity
the quadratic velocity profile shown. included, is called the Navier-Stokes equation.
In the rest of this chapter, we look at smooth (time-
Now, think about a little differential volume dxdy, steady) flows in which viscosity plays the dominant
as shown in the figure. It feels a net force to the right, role. This situation is often called laminar flow. Faber
due to dp/dx; what balances this force to keep the flow describes laminar flow in a more restricted sense, im-
from accelerating? This must be viscosity: friction on plying that “the fluid can be treated as an assembly
the little volume due to the shear in the flow. We as- of laminae of uniform thickness, whose boundaries re-
sume the shear stresses are linearly proportional to the main fixed as the flow moves between them”. The more
velocity gradient, transverse to the surface (this is called general, and more usual, definition, uses laminar flow
a Newtonian flow). The little volume will feel friction to mean smooth, non-turbulent flow. In most applica-
forces on its top and bottom surfaces; the net force on tions, we note that a real fluid does not slip freely past
the volume will be the difference between the forces a boundary (as was the implicit assumption in potential
on the top and bottom. If any one surface feels a force flow theory), but rather sticks: we will generally use
ρν(∂v/dy), the net force (per length in the z direction) no-slip boundary conditions.
on our element is
1
" What does ∇2 a mean, if a is a vector? For Cartesian it’s simple:
¶ #
∂2v
µ ¶ µ
∂v ∂v
∇ 2 a = ∇ 2 a x , ∇ 2 a y , ∇ 2 az
` ´
ρν − ρν dx = ρν dx
∂y y+dy ∂y y dy 2
(2.1) where
∂2 ∂2 ∂2
„ «
In this last I’ve assumed dy is small, of course, and also ∇2 = , ,
∂x2 ∂y 2 ∂z 2
that ν is independent of y. I’ve also written the constant is the usual Laplacian operator. For cylindrical and spherical co-
of proportionality as ρν, and the constant ν is what I ordinates, ∇2 a is a more complicated form – check the vector
like to call the coefficient of viscosity; I’ll return to this references I’ll put up on the web, or your favorite vectors-in-
shortly. funny-coordinates reference book.
9
(strictly, the mean over the z-direction): to get (2.5), we argue here that each side of 2.13 must
be constant, giving a flow profile
Q b2 dp
V = =− (2.10)
b 12ρν dx r2 dp
v(r) = + A ln r + B (2.14)
4ρν dz
(the minus sign just notes that V points opposite to the
pressure gradient). But I can write the right hand side where A, B are again integration constants, set by the
as the gradient of a scalar: boundary conditions.
µ 2 ¶
b
V =∇ p (2.11)
12ρν
That is, the horizontal velocity can be treated in terms
of a potential, φ = −(b2 p/12ρν). This shows that all
the potential flow apparatus of chapter 3 can be used
here; and, conversely, that flows where this works are
good illustrations for two-dimensional potential flow Figure 2.3. The circular Poiseuille flow solution,
solutions. showing the velocity field and also the stress field (τ ; our σ)
From Kundu figure 9.5.
But then: what does it take to justify this? We must
be able to ignore the viscous forces in the x and y direc-
But now, A = 0 to keep v finite as r → 0; and a
tions, and to ignore the inertial force in the x direction.
no-slip boundary at r = a gives
That is, we need
∂vx ∂ 2 vx ∂ 2 vx ∂ 2 vx r2 − a2 dp
vx ≪ν 2 ; ≪ v(r) = (2.15)
∂x ∂z ∂y 2 ∂z 2 4ρν dz
But now: the second condition simply requires that that is, another centered parabolic profile. The mass
L ≫ b; while the first requires flux here is
Z a
V2 νV b2 1 πa4 dp
≪ 2 ; ≪ (2.12) Q=ρ 2πvrdr = − (2.16)
L b L2 Re 0 8ν dz
Thus, the Reynolds number determines the required di- This is called Poiseuille’s Law, or the Hagen-Poiseuille
mensions of the apparatus. law.
C. Steady Flow: Cylindrical Applications 2. CIRCULATING FLOW
We can repeat the exercise for cylindrical geometry; Another case is flow between two concentric, rotating
differences from the planar cases highlight the effects cylinders, called circular Couette flow.
of geometry. The flow is now taken to be in the z direc-
tion, as is customary in cylindrical coordinates.
1. PIPE FLOW
tion of motion are the shear force points along the interface between the
vφ2 two fluid parcels). The coefficient ν is the viscosity or
1 dp
= shear viscosity.
r ρ dr
· ¸ (2.17)
d 1 d
0=ν (rvφ )
dr r dr
The general solution for the azimuthal velocity is
B
vφ (r) = Ar + (2.18)
r
and the pressure gradient, required to offset centrifugal
force, is
B 2
µ ¶
dp ρ
= Ar + (2.19)
dr r r
Using boundary conditions vφ (R1 ) = Ω1 R1 and
vφ (R2 ) = Ω2 R2 gives the specific values for A and
B.
This has some interesing limits. One has Ω = Ω2 Figure 2.5. The two-dimensional geometry for the stress
and R1 = 0, that is a rotating cylindrical tank. This has tensor. pi are the normal forces on the fluid square, and sij
are the shear forces. From Faber, Figure 1.2
vφ = Ωr (2.20)
showing that the fluid inside goes into solid body ro- Life is a bit more complicated, however. Consider
tation. Alternatively, if the outer cylinder is at infin- a two-dimensional Cartesian system, where s12 is the
ity, with Ω2 = 0, and the inner one has R = R1 and force in direction 2, due to adjacent fluid in direction 1;
Ω = Ω1 , then the solution is and s21 is the force in direction 1, due to adjacent fluid
in direction 2. These two forces must be equal: s12 =
ΩR2 s21 ; if not we would have a net torque on the parcel of
vφ = (2.21)
r fluid. Thus, we must symmetrize our expression for the
We’ll see this again in chapter 4 when we work with the shear force:
µ ¶
irrotational vortex. ∂v1 ∂v2
s12 = s21 = ρν + = s3 (2.22)
∂x2 ∂x1
D. Viscous Stresses, Generally
where the last equality simply defines s3 , to shorten the
OK, we’ve avoided this long enough .. we need to write notation.
down the more general form of the viscous stresses. There are also normal forces on the parcel: due to
This will be a tensor. We start by determining the sur- the fluid pressure, and also due to the fluid deformation.
face forces on a piece of fluid, due to deformations Faber approaches this by considering a frame rotation,
(arising from velocity gradients) of the adjacent bits of into the 45◦ degree primed frame. From force balance,
fluid. we find that
1. DO IT PHYSICALLY FIRST 1
s′3 = (p1 − p2 )
2 (2.23)
Following Faber here. Consider a planar flow in the x̂ 1
′
direction, with transverse velocity gradients (along ŷ). p1 = (p1 + p2 − 2s3 )
2
We expect the shear stress to be linear: ν x̂dvx /dy is the
force between adjacent streams in the fluid.2 (Note that and also, from a Taylor expansion, we find that
µ ′
∂v2′
¶
∂v1
s3 = ρν −
2
This is a BIG assumption – it’s traditional, but not obvious to ∂x′1 ∂x′2
me that it should always hold. Look ahead to the effects of a from which we get
magnetic field and j × B forces, for a possible counter-example.
Luckily, however, we rarely need to worry about viscous stresses ∂v1 ∂v2 ∂v3
in MHD applications ....
p1 + 2ρν = p2 + 2ρν = p3 + 2ρν (2.24)
∂x1 ∂x2 ∂x3
12
But now, taking the mean (isotropic) pressure to be and we note that the Ωij term is 1/2 of the usual vor-
1 ticity, Ω = ∇ × v. Also, note that the trace of Dij
p= (p1 + p2 + p3 ) (2.25) is
3
we end up with an expression for the normal force: ∂v1 ∂v2 ∂v3
Dmm = + + =∇·v (2.30)
2
µ
∂v1 ∂v2 ∂v3
¶ ∂x1 ∂x2 ∂x3
p1 = p − ρν 2 − − (2.26)
3 ∂x1 ∂x2 ∂x3 Now, we argue that Σij should depend only on the sym-
metric part of the velocity gradient, Dij (otherwise,
So, generalizing to three dimensions, we have all 9 simple rotation, such as solid body, would lead to a spu-
components (6 independent) of the surface forces act- rious stress term). Thus, requiring linearity, we have in
ing on this parcel, which we collect as the tensor ~~σ : general
p1 s12 s13 Σij = αijmn Dmn (2.31)
σij = − s21 p2 s23 (2.27)
which involves no fewer than 81 separate α’s. However,
s31 s32 p3 we can make arguments about symmetry and isotropy
to reduce this to the most general form useful for us:
(the minus sign has been added for consistency with the µ ¶
usual treatment). 1
Σij = 2ρν Dij − δij Dmm + ρνb δij Dmm (2.32)
3
References
I’ve mostly followed Kundu and Faber here, as well as Thompson (whose discussion of viscosity and stress
tensors I like). For the non-Cartesian forms of the NS equation, Pozrikidis is one good reference.
∂2f
µ ¶ µ ¶
1 ∂ ∂f 1 ∂ ∂f 1
∇2 f = r2 + 2 sin θ + 2 2
r2 ∂r ∂r r sin θ ∂θ ∂θ r sin θ ∂φ2
Cylindrical Coordinates, Axisymmetric This describes radial flow away from a source (or in-
The velocity is connected to the potential by wards to a sink if C < 0). Note that ∇ · v 6= 0 at the
origin, and that the potential is undefined there. Inte-
∂φ ∂φ grating the mass flux around a source or sink shows us
vr = ; vz = (3.12)
∂r ∂z that ρC describes the amount of mass per time being
added at the source, per length in the z direction (recall
The incompressibility condition is that this is a 2-D problem in 3-D).
∂vz 1 ∂ Another basic building block is that of a constant
+ (rur ) = 0 (3.13) flow of speed U in the x direction, represented in Carte-
∂z r ∂r
sian by
From this we can define a stream functionn:
φ = Ux ; v = ∇φ = U x̂ (3.20)
1 ∂ψ 1 ∂ψ
vz = ; vr = − (3.14)
r ∂r r ∂z 1. SOURCES AND SINKS , FLOW PAST HALF BODY
Spherical Polar Coordinates, Axisymmetric We can combine these two building blocks (and work
The velocity is connected to the potential by in plane polar):
∂φ 1 ∂φ C
vr = ; vθ = (3.15) φ = U r cos θ + ln r
∂r r ∂θ 2π
C (3.21)
The incompressibility condition is vr = U cos θ +
2πr
1 ∂ ¡ 2 ¢ 1 ∂
r vr + (vθ sin θ) = 0 (3.16) vθ = −U sin θ
r ∂r sin θ ∂θ
From this we get a stream function: to get what is called flow past a half body. We see why,
by noting that the corresponding stream function is
1 ∂ψ 1 ∂ψ
vr = ; vθ = − (3.17) C
r2 sin θ ∂θ r sin θ ∂r ψ = U r sin θ + θ. (3.22)
2π
B. Two-dimensional Planar Problems
Consider Laplace’s equation in two dimensions, in
plane polar coordinates (r, θ).
1 ∂2φ
µ ¶
2 1 ∂ ∂φ
∇ φ= r + 2 2 =0
r ∂r ∂r r ∂θ
We know the possible general solutions:
φ = constant ; φ ∝ ln r ; φ ∝ rn cos(nθ);
φ ∝ rn sin(nθ) ; φ ∝ r−n cos(nθ) ; φ ∝ r−n sin(nθ); Figure 3.1. Potential flow past a two-dimensional “half
body” (represented mathematically by a source function).
(3.18) The boundary streamline is labelled by ψ = m/2 (that’s
C/2 in our notation); the shape of this streamline defines
the shape of the half body. S shows the stagnation point in
And, of course, any superposition of these is a solution. front. From Kundu figure 6.9
So, we proceed by finding useful superpositions . . .
One important building block will be a simple The flow clearly has a stagnation point: if it’s at po-
source or sink: sition (a, π), then v(a, π) = 0 requires a = C/(2πU ).
Now, the value of the stream function passing through
Cθ C this point – its label – is
ψ= ; φ= ln r
2π 2π (3.19)
C C C
vr = ; vθ = 0 ψs = U a sin π + π=
2πr 2π 2
18
C C
U r sin θ + θ= (3.23)
2π 2
This stagnation streamline separates the flow into an
upstream region (v → U x̂ at ∞), and the region in-
side the streamline, a region with a smooth nose, called
a half body. We can, therefore, put a physical body with
this shape into a flow, and keep the same external solu-
tion.
Once we have the velocity solution, we can find the
drag on the half body – the net force the fluid exerts on
Figure 3.2. A doublet: a source and a sink which
the body (or vice versa). As long as we ignore friction approach each other, keeping the product of their strength C
(i.e., viscosity), the net force on the body is just the and separation ǫ constant. The streamlines become circles,
(vector) integral of the pressure over the surface of the whose centers lie on the y-axis and which are tangent to the
body. That is, for this simple case, we expect the drag x-axis at the origin. From Kundu figure 6.8
force to be
Now: combine this with a uniform stream flow, φ =
Z
FD = − pn̂dS (3.24) U x and ψ = U y: combining and putting in polars, with
S
a2 = µ/U , gives
if n̂ is the normal to the surface S, and the integral is
a2 a2
µ ¶ µ ¶
taken over the surface of the body. Bernoulli’s relation
for an incompressible fluid (1.19) can be written as φ=U r+ cos θ ; ψ=U r− sin θ
r r
a2 a2
µ ¶ µ ¶
1 1
p + ρv 2 = p∞ + ρU 2 (3.25) vr = U 1 − 2 cos θ ;
r
vθ = −U 1 + 2 sin θ
r
2 2
where we have evaluated the constant at “infinity”, that (3.27)
is far upstream. Work this out: the (normalized) excess
pressure, But, now, this shows that ψ(a, θ) = 0: the streamline
ψ = 0 is a circular cylinder of radius a. Thus, flows
Cp = 2(p − p∞ )/ρU 2 inside this cylinder have no influence on flows outside;
we can again replace this doublet with a physical cylin-
turns out to be positive at the front of the body, zero at der, and find the same flow outside.
θ = 113◦ , and negative along the sides. You can show
from this, by setting up the integrals, that the net pres-
sure, integrated over the surface of the body, is zero:
there is no drag in this system.
2. FLOW PAST A CYLINDER
(r, θ)
θ α
a Figure 3.5. Illustrating the nature of the flow field around
0 ξ
an arbitrarily shaped body. From Currie Figure 9.1.
Figure 3.4. The geometry of a line source, which extends
from ξ = 0 to ξ = a.
surface of the body, and then integrate it over the sur-
ξ = a (the end of the source) to the observation point. face. Or we can try a different approach.
The net stream function is then P Q
Z a Z a
k k
ψ(r, θ) = dψ = − cos α(ξ)dξ = (r−r1 )
0 4π 0 4π
0000000
1111111
(3.39) 1111111
0000000
0000000
1111111
0000000
1111111
0000000
1111111 u(y)
where the last step is obtained by straightforward geo- U
o
0000000
1111111
0000000
1111111
0000000
1111111
0000000
1111111
metrical identities and changes
R of the integration vari- 0000000
1111111
0000000
1111111
0000000
1111111
p 0000000
1111111
able. [Exercise: what is φ = dφ?] o 0000000
1111111
D. d’Alembert’s “Paradox” y
S R
We have now seen several examples of flow past a solid x
body, described by potential flow. In each case we Figure 3.6. Momenteum balance of flow past an arbitrary
found that the net force exerted by the flow on the body body. The upstream flow has constant speed Uo and
is zero. Thus, there is no drag in a potential flow around pressure po . The downstream velocity may be a function of
a body. This result turns out to hold for any body, sym- y; in this sketch it is lower immediately behind the body (as
metric or not. due to a wake). Note, the body need not be spherical, that
was just easy for me to draw... Following Kundu figure 4.7.
Does this sound wrong? We know experimentally
that a drag force does exist when a body is immersed in Consider in figure 3.6 the region of fluid within the
a fluid flow. This apparent contradiction, between the rectangle SRQP; the net drag on the body within this
theory and experiment, became known as d’Alembert’s region is just the rate of change of momenteum flux:
paradox. As you might suspect, the solution of the para- FD = Ṗout , i.e. the net outflow of x-momentum from
dox is that the theory is incomplete (no, we don’t say the region. The outflow through PS and QR is easy to
“wrong”). In particular, we have ignored two important calculate:
pieces of the physics here. One, is that these solutions Z L
can have nonzero velocity at a solid surface: check Fig- PS
ure 3.1 or 3.3. This is naive: any real fluid has a finite, Ṗ =− Uo ρUo dy = −LρUo2 ;
o
if small, viscosity. But viscosity will force to a no-slip Z L (3.40)
boundary condition. This leads to a viscous boundary QR 2
Ṗ =ρ u(y) dy
layer, which will modify the solutions close to the sur- 0
face, and add tangential stress which contributes some
Now: there must be mass outflow, Ṁ , through the
drag. Two, the flow behind the body may very well be
sides, SR and PQ (this is necessary to conserve mass
turbulent – as this viscous boundary separates from the
if u(y) < Uo , as shown). Mass conservation requires
body and connects into the wake. We will see later that
this also contributes to a drag force. Both of these ideas Z L
PQ SR
are illustrated in Figure 3.5. LρUo = Ṁ + Ṁ + ρ u(y)dy (3.41)
O
We can also carry out a general drag calculation, for
a finite body, that illustrates this argument. We could But if the sides are far enough away from the rock that
work out the pressure in the fluid, everywhere on the the x velocity ≃ Uo there, then the rate of outflow of
21
x-momentum is
Z L
Ṗ P Q + Ṗ SR = ρUo (Uo − u(y))dy (3.42)
0
How some have been ill-posed, some sin- where Stokes’ theorem is used in the second equality.
gular, As we noted above, K can be taken as a vector, with
Some poisoned by their self induction, direction given by the sense of the flow, using the right
some core size killed, hand rule. We can think of K as a “flux” of vorticity;
Some haunted by the mathematics they it measures the amount of ω passing through some sur-
have involved. face bounded by dl.
Let’s look at some examples. Solid body rotation is
All murderous.
described by vθ = Ωr;
To proceed, take the curl of (4.5). Using the vector re- Proof #1 Kundu proves this as follows. Start with Eu-
lation v · ∇v = (∇ × v) × v + 12 ∇v 2 , we get the ler’s equation, (1.10) (note we aren’t worried about vis-
dynamical equation for vorticity: cosity right now):
ω
∂ω 1 Dv 1
ω × v) = ν∇2ω + 2 ∇ρ × ∇p (4.6)
+ ∇ × (ω = − ∇p + ∇Φ (4.11)
∂t ρ Dt ρ
• The behavior of vortex lines can be very in- Figure 4.4. A vortex ring approching a wall. From
triguing. They move under their own power. A vor- Kundu figure 5.14
tex line will move in a transverse flow, in a direction
U × K. (Recall the Magnus effect). A bent vortex line • Vortex rings are formed when a vortex line con-
will move in a direction K × n̂, if n̂ is the direction nects back on itself. A vortex ring will move itself along
of its radius of curvature R; Batchelor shows its speed through the background fluid, at a speed ∼ K/4πb if
∼ K/4π. b is the radius of the ring. (Any two opposite points
• Images are useful in calculating the behavior of around the ring act like paired, anti-parallel K lines.
vortex lines near boundaries (remember your images Alternatively, this is a case of well-defined curvature of
in electrostatics?). A wall, such as in Figure 4.3, is a a single vortex line.) One vortex ring will spread out as
streamline, and must have zero normal velocity. This it approaches a wall. Two rings will move around and
means we can add an imaginary image vortex, behind through each other – Faber shows this in his figures.
the wall and paired with the real one, to understand the • Vortex lines can reconnect. Viscosity breaks the
motion of the vortex line. conservation of circulation, and allows adjacent vor-
26
tex lines or rings to reconnect, thereby changing their and Kundu (for the math); but there’s also good ma-
topology. terial in Shivamoggi, not to mention a whole book by
Saffman (quite mathematical), and several useful ar-
F. Generation of Vorticity ticles in Annual Review of Fluid Mechanics over the
We argued above that circulation is conserved; this was years.
Kelvin’s theorem. How, then, can a fluid which starts ir-
rotational, ever gain vorticity? The answer, once again,
is the combination of viscosity and boundary effects.
To illustrate, I combine a figure from Kundu with the
discussion from Batchelor.
References
Vorticity is a large, interesting and complex topic.
I’ve mostly followed Faber and Tritton (for the words)
27
5. LAMINAR FLOW: MORE APPLICA- rock is not very disturbed by the rock’s motion (that is,
TIONS the flow will vary with height, above the rock). Now
consider the same thing but with a rotating water tank.
We continue the previous discussion with more exam- The Taylor-Proudman result says that the flow at all
ples and applications from here and there. We start heights above the rock must be the same. This produces
with a simple discussion of flow in a (quasi) two- a Taylor column: the water above the rock forms a ver-
dimensional, rotating system – all sorts of fun things tical column which is fixed relative to the rock; water
happen when you’re rotating. We then look at two more not over the rock flows around this column, as though
mathematically involved examples of viscous flows. it were a solid object.
A. Geostrophic Flow ATMOSPHERIC CIRCULATION
What happens in a rotating system? Because our impor- One striking example of our first point is the at-
tant application is to terrestrial (geophysical) flows, we mospheric circulation around a pressure extremum:
can work in the incompressible limit. Referring back to the circulation is counterclockwise (seen from above)
chapter 1, the general equation of motion comes from around a low, and clockwise around a high – as shown
(1.20), but now has viscosity added, as per Chapter 2: in the figure.
Dv 1
= − ∇p + g + Ω2 R − 2Ω
Ω × v + ν∇2 v (5.1)
Dt ρ
pressure High pressure
We consider the limit which is interesting for terrestrial Low
coriolis coriolis
in this geometry would have to follow contours of con- The boundary conditions are: v = 0 at z = 0 (no-slip
stant depth h. Consider, then, a small radial displace- at the solid surface), and v = U x̂ as z → ∞. This
ment of a column of the fluid, to a region with h + δh. system solves to
This violates the T-P theorem; the column will try to h i
return to its initial position. But it has some inertia, so vx = U 1 − e−z/δ cos(z/δ)
it will overshoot. We thus have an oscillation: this is (5.4)
one version of a Rossby wave. vy = U e−z/δ sin(z/δ)
Rossby waves exist in the atmosphere, and in the
ocean; they are found to be low frequency ω < Ω with δ 2 = ν/Ω. Thus, the velocity vector rotates
and long wavelength. The restoring force for terres- through a spiral, while it increases from zero at the
ground to ∼ U x̂ at z ∼> δ.
trial waves is due to variations in the depth of the atmo-
sphere parallel to Ω (here meaning the earth’s rotation). • Free boundary. Our driver here is the wind stress
In the atmosphere they are easily observed as the large- at the ocean surface, again in the x direction; call it τ x̂.
scale meanders of the mid-latitude jet stream. In the We can ignore pressure gradients in this case; thus our
ocean they are harder to find, being very small ampli- equations are, simply,
tude (several cm) and very long wavelength (hundreds
of km), but now detected in satellite data. d 2 vx
−2Ωvy = ν
dz 2 (5.5)
d 2 vy
Ω 2Ωvx = ν 2
dz
α Take z = 0 at the ocean surface, and negative below.
Our boundary conditions now are v → 0, z → −∞,
and dvy /dz = 0, dvx /dz = τ /ρν, z = 0 (refer back to
h
the definitions of stress in chapter 4). These solve to
³ z π´
vx = Vo ez/δ cos − +
Figure 5.2. The geometry for a simple version of Rossby δ 4 (5.6)
waves. The curved arrow marks the deviation of the upper z/δ
³ z π´
plate from parallel. Following Tritton Figure 16.13. vy = −Vo e sin − +
δ 4
√
THE EKMAN LAYER where Vo = τ /ρ 2Ων, and δ is the same as above.
Again, the velocity vector rotates through a spiral; note
Coriolis forces also do interesting things at boundaries. that it is at a 45◦ angle to ~τ at the ocean surface.
The two applications relevant to geophysical flows are Both cases have similar consequences. A prevail-
solid boundaries (such as land) and free boundaries ing wind or flow results in a transverse flow within the
(such as the ocean). In each case, the combination of Ekman boundary layer (this is Ekman transport). Mass
a driving wind (think of geostrophic flow) and viscos- conservation then forces a net vertical motion (why?)
ity within the boundary layer results in a transverse flow into or out of the layer (this is Ekman pumping). In
within that layer. In this section I present an overview; the atmosphere this is related to updrafts in low pres-
details will appear in the homework. sure regions (thus storm formation); in the ocean this
• Solid boundary. At some high altitude, let the is related to phenomena such as upwelling of cold sub-
wind velocity be U x̂. Thus, (5.2) reduces at high alti- surface water in regions with a steady prevailing wind
tudes to direction.
1 dp
2ΩU = − B. Viscous Flow: Time-Dependent Problems
ρ dy
Viscosity is dissipative: we expect a viscous flow to
(note dp/dx = 0, right?). Within the boundary layer
decay with time. In this section we consider solutions
we must keep the viscous terms. Here, (5.1) becomes
of the basic equations, (4.5) or (4.6), but here omitting
d 2 vx gravity and the advective terms. Thus we want solu-
−2Ωvy = ν tions of
dz 2
(5.3)
d 2 vy 1 dp ∂v
2Ωvx = ν 2 −
dz ρ dy ρ + ∇p = ρν∇2 v (5.7)
∂t
29
1
·
1 ∂2ψ 1 ∂
µ
1 ∂ψ
¶¸ (Note that we couldn’t use Bernoulli’s law here; why?)
ωφ = + 2 (5.19) Thus, the pressure is highest at the forward stagnation
r sin θ ∂r2 r ∂θ sin θ ∂θ
point, and lowest (negative in fact) at the rear stagna-
Putting this into (5.17), we get a fourth order equation tion point. For the rest of the stress tensor, we use the
for ψ: velocity solution, from (5.23), to get
3a3
µ ¶
·
∂2 sin θ ∂
µ
1 ∂
¶¸2 ∂vr 3a
+ 2 ψ=0 (5.20) Σrr = 2ν = 2νU cos θ − 4
∂r 2 r ∂θ sin θ ∂θ ∂r 2r2 2r
3νU a3
· ¸
∂ vθ
³ ´ 1 ∂vr
This is now something that can be solved. Our bound- Σrθ = ν r + =− sin θ
∂r r r ∂θ 2r4
ary conditions are
(5.26)
ψ(a, θ) = 0 Putting this into the expression for FD , multiplying by
4πa2 and integrating over θ, gives the drag
dψ/dr(a, θ) = 0 (5.21)
1 FD = 6πρνaU (5.27)
ψ(∞, θ) = U r2 sin2 θ
2
3
Compare our earlier discussion in Chapter 3, in and around
2
For once the name makes sense! (3.24).
31
References
I’ve mostly followed Kundu, who does the more
mathematical examples. He also has a good chapter
on geophysical fluid dynamics. Tritton is also good for
geostrophic flows – as always his discussions are help-
ful (not just pure math).
32
Putting this into (6.5) and differentiating with r, the ba- Numerical solutions of this equation are called poly-
sic equation becomes tropes, or Lane-Emden solutions, for polytropic index
n = 1/(γ − 1). These are physically well-behaved, in
d r2 dρ
µ ¶
G that the density is finite at the center and falls faster than
=− 4πr2 ρ (6.7)
dr ρ dr RT 1/r2 at large r.
The solutions of this are less straightforward. First, we 3. REALITY: NONISOTHERMAL ATMOSPHERES
note that a basic scale length appears:
Just a note of caution here. Both of the preceding ex-
9RT 1/2
µ ¶
amples assumed a constant temperature throughout the
ao = (6.8)
4πGρo atmosphere. This naive assumption allowed us to find
(nearly) analytic solutions. More realistically, however,
if ρo is some characteristic density, say that at the ori- we expect the temperature to be a function of position,
gin. (the numerical constants appear to simplify things controlled by the thermodynamics of the system. For
later on.) We expect physical solutions of (6.7) to in- earth’s atmosphere, T (z) is determined by the com-
volve lengths scaled to ao . bined actions of solar radiation coming in at the top,
Consider an isothermal gas; these are the simplest all the chemistry and energy transfer effects within the
solutions. One solution of (6.7) is the simple power atmosphere, and radiative losses back into space. For a
law, ρ ∝ 1/r2 . The divergence as r → 0 keeps this star, T (r) is determined by nuclear energy generation
from being an interesting solution, however. More in- in the star’s core, radiative energy transfer within the
teresting physical solutions must be found numerically, star, and eventual radiative losses into space.
starting with ρ(r = 0) = ρo and working out. These
solutions do, indeed, have a turnover at r ≃ ao ; at large 4. ADIABATIC ATMOSPHERE
radii they do approach ρ ∝ 1/r2 . (These solutions also
Another analytic model of an atmosphere assumes the
have problems, for the total mass M (r) diverges ∝ ln r
gas is adiabatic: that is, as it compresses or expands
as r → ∞. Physically, the problem is that the gas can-
(while satisfying HSEq), it heats or cools accordingly.
not be maintained isothermal everywhere.)
This is useful in analyzing the stability of an atmo-
sphere to convection (which we’ll do immediately be-
low).
To develop this idea, go back to basic HSEq, (6.1):
dp/dz = −ρg. But now, assume the gas is adiabatic,
so that
µ ¶(γ−1)/γ µ ¶1/γ
T p ρ p
= ; = (6.10)
To po ρo po
We can use these to relate dT /dz to dp/dz for an adi-
abatic atmosphere:
1 dT γ − 1 1 dp 1 dρ 1 1 dp
= ; = (6.11)
T dz γ p dz ρ dz γ p dz
Figure 6.1. Density (and projected density, Σ) solution
for a self-gravitating isothermal sphere. The dotted line is From here, assuming p = ρRT (ideal gas) and using
the asymptotic ρ ∝ 1/r2 power-law solution. From Binney cp − cv = R (refer back to §6.1) we get the condition
& Tremaine figure 6.7 for an adiabatic atmosphere:
Self-gravitating spheres can also be modelled as adi- dTad g (γ − 1) mg
abatic gases. Take the equation of state as p = Kργ . =− =− (6.12)
dz cp γ kB
The hydrostatic equation, (6.5), and the differential
form (6.7), become where cp is the isobaric specific heat, and m is the mass
per particle. This clearly gives a linear temperature
dρ dΦ drop with altitude z; g/cp ≃ 10◦ C/km for typical atmo-
Kγργ−2 = −ρ ;
dr µdr spheric conditions. The vertical temperature gradient is
¶ (6.9)
1 d 2 dΦ often caled Γ = dT /dz; and is sometimes called the
r = 4πGρ
r2 dr dr “lapse rate”.
34
measures the rate of change of a quantity, due to its mo- cluding viscosity, equation (2.2); multiply it by vi and
tion through a region in which the flow field changes. interchange dummy indices, to get an expression for the
We get, time-change of kinetic energy:
µ ¶ µ ¶
D 1 2 D p D v2
ρ e+ v +ρ = ρg (6.33) ρ = σik,i uk + ρgk vk (6.35)
Dt 2 Dt ρ Dt 2
Now, using g = ∇Φ, Bernoulli’s relationship becomes Subtract this from the full energy equation (also written
in Cartesian), and symmetrize the tensor term (in the
1 p 1
e + v 2 + + Φ = h + v 2 + Φ = constant (6.34) second step), to get two forms of the energy equation
2 ρ 2 including viscosity:
which does, indeed, recover the Bernoulli relation that
De De dvk
we derived above from momentum conservation. We ρ = σik vk,i = σik Dik ; ρ= −p +Σik Dik
recall, again, that this law holds along any one stream- Dt Dt dxk
(6.36)
line in the flow.
This brings back the (symmetrized) deformation tensor
in Cartesian:
µ ¶
1 ∂vi ∂vj
References Dij = + (6.37)
2 ∂xj ∂xi
I mostly follow Thompson for the basic energy
conservation and dissipation analysis. The isothermal and introduces what I’ll call the dissipation function:
sphere discussion follows Binney & Tremaine (Galac- D = Σij Dij . Referring back to (2.37) and (2.38),
tic Dynamics); the terrestrial atmosphere discussion it’s normal to assume νb = 2ν/3, which simplifies the
leans on Kundu and Tritton. stress and dissipation tensors. In this limit, still Carte-
sian, the energy equation becomes
µ ¶2
∂e ∂e ∂vk ∂vk
ρ + ρvi +p = ρν (6.38)
∂t ∂xi ∂xk ∂xk
E. Appendix: Viscous Dissipation
Well, that was so much fun, let’s do it again: in Carte- To close, I write out explicitly the dissipation func-
sian to be explicit. Start with our force equation, in- tions for all 3 coordinate systems.
Cartesian:
2ρν £
(D11 − D22 )2 + (D22 − D33 )2 + (D33 − D11 )2
¤
D=
3 (6.39)
+ ρνb (D11 + D22 + D33 )2
¡ 2 2 2
¢
+ 4ρν D12 + D13 + D23
where, as above,
µ ¶
1 ∂vi ∂vj
Dij = + (6.40)
2 ∂xj ∂xi
Cylindrical:
¡ 2 2 2 2 2 2
+ ρ(νb − 2ν/3)(∇ · v)2
¢
D = 2ρν Drr + Dθθ + Dzz + 2Drθ + 2Dθz + 2Dzr (6.41)
where
µ ¶
∂vr 1 ∂vθ vr 1 1 ∂vr ∂ vθ
Drr = ; Dθθ = + ; Drθ = +r
∂r r ∂θ r 2 r ∂θ ∂r r
µ ¶ µ ¶ (6.42)
∂vz 1 ∂vθ 1 ∂vz 1 ∂vr ∂vz
Dzz = ; Dθz = + ; Drz = +
∂z 2 dz r ∂θ 2 ∂z ∂r
38
Spherical is even better (recall that θ is the polar angle and φ is the azimuthal angle):
¡ 2 2 2 2 2 2
+ ρ(νb − 2ν/3)(∇ · v)2
¢
D = 2ρν Drr + Dθθ + Dφφ + 2Drθ + 2Dθφ + 2Dφr (6.43)
where
∂vr 1 ∂vφ vr vθ cot θ 1 ∂vθ vr
Drr = ; Dφφ = + + ; Dθθ = +
∂r r sin θ ∂φ r r r ∂θ r
µ ¶ µ ¶
1 1 ∂vr ∂ vφ 1 1 ∂vr ∂ vθ
Drφ = +r ; Drθ = +r (6.44)
2 r sin θ ∂φ ∂r r 2 r ∂θ ∂r r
µ ¶
1 sin θ ∂ vφ 1 ∂vθ
Dθφ = +
2 r ∂θ sin θ r sin θ ∂φ
We can also demonstrate that cs is the speed of a travel- SOUND WAVES : A FORMAL APPROACH
ling wave. Let some perturbation (δρ, δp, δT ) be mov-
We can also use a more mathematical approach, deriv-
ing at some cs . Ahead of the wave the fluid has v = 0;
ing a formal wave equation by linearizing our two basic
behind the wave the fluid has δv, in the same direction
equations, continuity & momentum:
as the wave motion.
∂ρ
+ ∇ · (ρv) = 0 (7.5)
∂t
(assuming no body forces)
∂v
ρ + ρ(v · ∇)v = −∇p (7.6)
∂t
Our question is, at what speed does a small, compres-
sional disturbance in the gas travel? We start with a
uniform, static unperturbed state, described by ρo , po
and vo = 0. We add small perturbations, ρ1 , p1 and v1 .
If we put these into the mass and momentum equations,
(7.5) and (7.6), noting that the “o” terms are constant in
space and in time, we get
∂ρ1
Figure 7.1. Propagation of a sound wave: (a) into a still + v1 · ∇ρo + (ρo + ρ1 )∇ · v = 0
∂t
fluid; (b) a stationary wave. From Kundu figure 15.1. Note,
u here is the same as v in the text. and
∂v1
Mass conservation at the wave front, in a frame moving (ρo + ρ1 ) + (ρo + ρ1 )v1 · ∇v1 = −∇p1
with the wave front, gives ∂t
ρcs = (ρ + δρ)(cs − δv) In these equations, we now drop terms which are sec-
ond order in the perturbed quantities, and we write
and to lowest order in things small, this gives ∇p = (∂p/∂ρ)∇ρ. This gives us
δρ ∂ρ1
δv ≃ cs (7.2) + ρo ∇ · v 1 = 0 (7.7)
ρ ∂t
40
From these we immediately find the dispersion rela- Here, we have introduced the Mach number, M =
tion: v/cs .
Thus, we have verified that we have non-dispersive We have just demonstrated that the sound speed is the
waves propagating at speed cs . This again verifies our speed at which information propagates; it is also critical
guess in (7.1). to the dynamics of the flow.
The derivative in (7.1) is usually taken assuming the There are important difference between subsonic
disturbance is adiabatic, so that p/ργ is constant, and and supersonic flows. Subsonic flows can be thought of
c2s = γp/ρ. This corresponds, physically, to the heat- as quasi-hydrostatic. That is, the flow field is strongly
ing and cooling times for the perturbation to be long influenced by pressure gradients which are determined
compared to the wave travel time. An alternative choice by conditions a long distance away (such as at bound-
which is sometimes used, is to consider isothermal per- aries).
turbations: T = constant – which is the limit in which Supersonic flows, however, are quasi-ballistic. This
the local cooling and heating times are short compared distinction is due to the fact that information travels at a
to the wave travel time. In this case, c2s = p/ρ. finite speed, the speed of a simple sound wave. Both 1D
41
and 3D cartoons, in Figures 7.1 and 7.2, can illustrate • hypersonic: M ∼ > 3, say. Qualitatively the same as
this point. Pressure gradients have only a limited range supersonic flow, but some interesting new effects – such
of influence, and conditions far away have little or no as strong heating and ionization of boundary layers –
effect on a solution locally. We’ll see that supersonic come into play. Shocks can be analyzed in the strong-
flows can (and usually do) contain discontinuous jumps shock limit (chapter 9).
in the flow properties (shocks). They can violate our
subsonic intuition, for instance a supersonic flow in a C. Weak Waves and Causality
diverging channel will accelerate (as we’ll see below). How can we use the signal speed – the sound speed – to
perturbation understand a flow? One way to approach this, following
Undisturbed flow Currie, is to consider “weak waves”. Specialize to a
cs cs 1D system, and let co be the undisturbed value of cs .
reverse waves forward waves Equations (7.7) and (7.8) become
Figure 7.2. Physical illustration of simple waves. The
∂ρ1 ∂v1 ∂v1 ∂ρ1
information that the flow has been “whacked” at point a, + ρo = 0 ; ρo + c2s =0
propagates by simple sound waves, moving at speed cs ∂t ∂x ∂t ∂x
relative to the fluid in the pipe. Following Thompson figure In the previous derivation, cs can be a function of den-
8.6
sity, and thus can vary within the wave. Here, we sim-
plify by taking cs to be a constant, co . But now: be-
cause ρo is assumed constant, and because vo = 0, we
ct
s
can write
vt vt
∂ρ ∂v ∂v ∂ρ
cst + ρo = 0 ; ρo + c2o =0 (7.18)
∂t ∂x ∂t ∂x
M>1
M=0 M<1
Compare these to the originals, (7.7) and (7.8): we have
Figure 7.3. Mach’s construction for the propagation of a quietly removed the bothersome, nonlinear terms. This
disturbance. Consider a point source of sound (Thompson
suggests a bumblebee) in a moving medium. If the source is valid only as long as we continue to assume a small
and flow are stationary, the sound propagates spherically perturbation – a “weak wave”. We can now divide the
from the source. If the source/flow are moving subsonically. first of (7.18) by ρo , the second by ρo co , and add and
the motion only distorts the spherical wavefronts. If, subtract, to get
however, the motion is supersonic. all disturbances are µ ¶ µ ¶
confined to a Mach cone; an observer located outside of this ∂ v ρ ∂ v ρ
+ + co + =0 ;
cone does not receive any information about the bee. The ∂t co ρo ∂x co ρo
opening angle of the cone is called the Mach angle: µ ¶ µ ¶
sin µ = 1/M. ∂ v ρ ∂ v ρ
− − co − =0
∂t co ρo ∂x co ρo (7.19)
Some authors use the following classifications,
These now have the form of a total (Lagrangian) deriva-
based on the Mach number M = v/cs (I follow Kundu,
tive (compare equation 1.5). Thus, we have an impor-
for instance):
tant result:
• incompressible: M ∼ < 0.3 everywhere; can ignore µ ¶
density variations due to pressure changes. v ρ
+ = constant on x − co t = constant
< M < 1 everywhere. Need to be co ρ o
• subsonic: 0.3 ∼ ∼ (7.20)
careful with density fluctuation (now above the 10%
and
level), but no shock waves in the flow. µ ¶
< M < 1.2, say; M is “around v ρ
• transonic: 0.8 ∼ ∼ − = constant on x + co t = constant
unity”. Shock waves appear, increasing drag. These are co ρ o
the hardest flows to analyze analytically, as the nonlin- (7.21)
ear terms in the governing equations are important, but An alternate form of (7.20 and 7.21) can be found with
the simplifying effects of M > ∞ flow can’t be used a bit of algebra. Using (7.7) and ρ1 ≪ ρo , it’s easy to
yet. express ρ/ρo in terms of p/po and γ, to write
< M < 3, typically. Shock waves J+
µ ¶
• supersonic: 1 ∼ ∼ v 1 p
= + = constant
are usually present. Analysis is easier because informa- c0 co γ p o (7.22)
tion propagates along well-defined directions in (x, t)
space, called characteristics. on x − co t = constant
42
and the low-pressure gas. The information that this has hap-
pened can only travel at the sound speed, co . Thus, a
J−
µ ¶
v 1 p compression wave will move to the right, and an ex-
= − = constant
c0 co γ p o (7.23) pansion wave to the left, both at co . The problem: what
is the velocity and pressure of the gas everywhere, as a
on x + co t = constant
function of time?
The interpretation is simple for (7.20) or (7.22). The Figure 7.4 shows an (x, t) diagram for this system,
lines x = ±co t are the loci of forward and backward with the two wave paths as solid lines. Gas ahead of
propagating sound waves. Our results say that the quan- the forward wave must be undisturbed: it has p =
tities on the left of equations (7.21) and (7.23) are con- po , v = 0. Similarly, gas behind the reverse wave has
stant along these trajectories. p = p1 , v = 0. What of the middle region? Con-
sider an observation point P (xp , tp ). Two wave lines
Note that the above derivation is only valid for weak intersect P : one forward wave which originates from
shocks, working in a frame where the fluid speed is (x < xp , t = 0), and a reverse wave which originates
much less than the sound speed. However, these condi- from (x > xp , t = 0). But now: we know v = 0 ev-
tions may easily be relaxed. In a moving fluid, one finds erywhere at t = 0, so we can evaluate the constants in
that the characteristics are (not surprisingly) subject to (7.22, 7.23) for each of these wave lines. For the reverse
a Galilean transformation by speed v. The invariants, one, we have J − /c0 = −1/γ; and for the forward one,
and associated characteristics, become we have J + /c0 = p1 /γpo (verify this for yourself!).
µ Z
∂ρ
¶ Thus, at point P we know that
+
J = v + cs = constant
ρ (7.24) v 1 p 1 p1 v 1 p 1
+ = and − =−
on x − (cs + v)t = constant co γ p o γ po co γ p o γ
2. PISTON PROBLEM U
(piston)
Consider another 1D tube filled with gas; at t = 0 a
piston is moved into the gas, at velocity U . How does t
(v=U)
which starts from the piston at t = 0. Also again, con-
sider an observation point P (x, t). It connects the the
forward x-axis by a reverse wave; we can evaluate the
constant J − /c0 , because we know v = 0, p = po at
the x-axis. The point P also connects to the piston by a
forward wave. We cannot evaluate J + /c0 yet, however; (v=0, p=p0)
x
we know the velocity vp = U , but not the pressure, pp , Figure 7.5. Illustrating the geometry of the piston
at the piston. We thus need a third wave, another re- problem, in the x,t plane. The two solid lines are the locus
verse one which connects the piston to the x-axis. We of the piston (assumed to be close to vertical; U ≪ co ), and
thus have three algebraic equations: the forward wave starting from the origin. The dotted lines
are the two waves which intersect the observation point P ,
U 1 pp 1 and the third (reverse) wave used to connect the piston to
− =− the x > 0 axis. Following Currie Figure 11.5.
co γ p o γ
v 1 p U 1 pp
+ = + (7.27)
co γ p o co γ p o
v 1 p 1
− =−
co γ p o γ
This now solves the system: we can find pp , and then
solve for conditions in the gas between the piston and
the forward wave:
p U
v=U ; =γ +1 (7.28)
po co
3. WAVES AT BOUNDARIES
References
Good references here include Thompson, Currie,
Kundu and Faber.
44
We must note, however, that M = 1 is not nec- We also apply energy conservation, assuming an adia-
essarily reached at the throat. Other combinations of batic shock:
initial/boundary conditions allow other flows. Equation
(8.13) also allows smooth solutions where the velocity γ p1 1 2 γ p2 1 2
+ v1 = + v (8.16)
v is an extremum: everywhere subsonic flow can have a γ − 1 ρ1 2 γ − 1 ρ2 2 2
maximum of v at the throat, and everywhere supersonic
flow can have a minimum. Figure 8.5 illustrates this. Now, these equations can be used to express the three
post-shock quantities, ρ2 , v2 , and p2 , in terms of their
pre-shock counterparts. Working through the algebra
(see chapter 9) gives the normal shock jump conditions:
ρ1 γ−1 1 2
= + 2
ρ2 γ+1 M γ+1
p2 2γM2 − (γ − 1)
= (8.17)
p1 γ+1
v2 γ−1 1 2
= +
v1 γ + 1 M2 γ + 1
Figure 8.5. More convergent-divergent nozzles, in which Once you have these, the jumps in T and cS can also be
the velocity does not reach M = 1 at the throat. From derived from the usual adiabatic and ideal gas relations.
Kundu figure 15.7
In a nutshell: a shock decelerates the flow (v2 < v1 ).
The density must increase, due to mass conservation
(ρ2 > ρ1 ). The “lost” kinetic energy goes to heat: T2 >
3. NORMAL SHOCKS T1 and p2 > p1 .
Most flows that go between subsonic and supersonic And a caveat: equations (8.17) hold for the rest
aren’t so smooth, however (in fact it takes very special frame of the shock. In many applications,however, the
conditions to achieve the smooth, de Laval transition of shock moves through the fluid. If this is the case, you’ll
Figure 8.4). More typically, mismatched variations in need to remember to do a Galilean transform to the
channel area and/or external pressure lead to shocks in shock-rest frame before applying the jump conditions.
the flow. Chapter 9 is devoted to shock properties; but 4. 1D FLOWS WITH INTERNAL SHOCKS
we also need some short discussion here.
In brief, a shock is a discontinuous transition in the Now, let’s return to 1D flows, and consider what hap-
flow. If you force a supersonic flow to change suddenly pens if smooth flow is not reached.
(compared to the sound travel time), a shock forms at The pressure jump maintained across the nozzle can
the right place and strength to accomodate the change. determine the flow possibilities. From (8.12), the sign
Here we’ll look at normal shocks, which means of dp is opposite to the sign of dv; accelerating flow
shocks normal to the flow.2 Put yourself in a frame has a pressure drop, and vice versa. Consider, say, a
in which the shock is at rest; let “1” be conditions up- nozzle connecting an input region of pressure po , and
stream (flow coming into the shock) and “2” be condi- an exit at pe . What happens? (I follow Anderson in this
tions downstream. discussion; refer to Figure 8.6.)
We apply mass and momentum conservation across (a) For pe = po , there is no flow of course. Let pe
the shock. These two conditions become drop a bit (say to pa in the figure): the flow is com-
pletely subsonic, with highest M and lowest p at the
ρ 1 v1 = ρ 2 v 2 (8.14) throat.
(b) At some specific value of pe , the flow reaches
and sonic at the throat. It remains subsonic everywhere
else, however (pb is not small enough to allow a smooth
sonic transition); this is called choked flow.
ρ1 v12 + p1 = ρ2 v22 + p2 (8.15)
(c) As the exit pressure is further reduced, a region
of supersonic flow occurs downstream of the throat.
2
No, there are no “abnormal” shocks, sorry. Because it can’t continue so all the way to the exit, a
47
normal shock forms in the diverging nozzle; the flow flow. In particular, a smooth transition from subsonic to
slows down to subsonic at the shock. supersonic flow is possible if the gas stays hot enough
(d) For some specific value of pe , the shock is lo- (extended heating sources are required).
cated exactly at the exit. The fully isentropic, subsonic- The basic solution is due to Parker. Consider a
supersonic flow pattern now exists throughout the entire steady, spherical outflow. Mass conservation in this
duct, except at the exit. case is ρvr2 =constant; or,
1 dρ 1 dv 2
+ + =0 (8.18)
ρ dr v dr r
while the momentum equation becomes in this case
(noting that gravity from the central star is important),
dv dp GM
ρv + = −ρ 2 (8.19)
dr dr r
Writing dp/dr = c2s dρ/dr, these two equations com-
bine to give the basic wind equation,
c2s dv 2c2 GM
µ ¶
v− = s − 2 (8.20)
v dr r r
This does not have analytic solutions over the whole
range of r. It must be solved numerically; examples are
shown in Figure 8.7. However, as with our analysis of
channel flow, we can learn a lot by simple inspection of
(8.20).
S2 11111111111111
00000000000000
Whether or not a particular flow satisfies this condition C
00000000000000
11111111111111
00000000000000
11111111111111
00000000000000
11111111111111
00000000000000
11111111111111
depends on the starting conditions, such as with what 00000000000000
11111111111111
00000000000000
11111111111111
00000000000000
11111111111111
S1 00000000000000
11111111111111
velocity and temperature it left the stellar surface, and 00000000000000
11111111111111
00000000000000
11111111111111
00000000000000
11111111111111
00000000000000
11111111111111
also what the boundary conditions at large distances 00000000000000
11111111111111
00000000000000
11111111111111
00000000000000
11111111111111
00000000000000
11111111111111
a 00000000000000
11111111111111
are. If it does not start in such a way to satisfy this con- 00000000000000
11111111111111
00000000000000
11111111111111
R(t)
00000000000000
11111111111111
b 00000000000000
11111111111111
dition, it either stays subsonic (corresponding to finite c
d
00000000000000
11111111111111
00000000000000
11111111111111
00000000000000
11111111111111
00000000000000
11111111111111
pressure at infinity), or cannot establish a steady flow. 00000000000000
11111111111111
00000000000000
11111111111111
(ii) The solution beyond the sonic point depends on Figure 8.8. Cartoon of the structure of a stellar wind,
and its interaction with the ISM. Left: the shock structure
the temperature structure of the wind. The only solu- within the wind. Right: The outer shell of dense,
tions with dv/dr > 0 for r > rs are those for which snowplowed ISM. From Dyson & Williams figures 7.3 and
c2s (r) drops off more slowly than 1/r; it is only these 7.4.
for which the right-hand side stays positive. In the case
of an isothermal wind, with c2s = constant, (8.20) can
S1 be the inner shock; region “b” be the wind-gas which
be solved in the limit r >> rs :
has been through the shock; C be the contact surface be-
v 2 (r) ≃ 4c2s ln r = constant (8.22) tween the wind and the ISM; region “c” be the shocked
ISM; and S2 be the outer shock (moving into the ISM).
Thus, the wind will be supersonic, by a factor of a We expect S1 to be an adiabatic shock (since the wind
few, as r → ∞. The question of how the solar wind is probably hot and low density, and thus will have a
manages to stay nearly isothermal is not solved; it is long cooling time); region “b” will contain hot, shocked
probably due to energy transport by some sort of waves 2
3 mvwind
wind, with Tb ∼ 16 kB ∼ several × 107 K (noting
(MHD or plasma waves, for instance) which are gener-
1
ated in the photosphere and damped somewhere far out that m = 2 mp is the mean mass per particle if region
in the wind. “b” is fully ionized). The outer shock will probably be
isothermal, since the ISM is denser and cooler than the
(iii) Inside the sonic point, the gravity term will
wind. Thus, the shocked ISM will be in a thin shell,
dominate the right hand side of (8.20). Thus, solutions
containing all of the original ISM that lay between S2
with dv/dr > 0 , and v 2 < c2s , will obey
and the star.
dv GM We can’t say anything about the details of the shock
v ≃− 2
dr r transitions until we know more about shocks ... and
This equation looks as if gravity is driving the wind out! that’s next.
This unlikely-looking result comes from the fact that
the flow is nearly subsonic in this region; therefore, the
dp/dr term in (8.19) – which actually drives the wind
out – is nearly equal to the gravity term. References
What happens at the outer edge of such a wind? Re- The basic 1D flow material can be found in sev-
member that the wind is not flowing into vacuum, but eral places; I’ve followed Thompson, Kundu and An-
rather into some external medium (call it the ISM, for derson, mostly. The more specialized spherical-wind
interstellar medium: this theory has been developed to applications are found in, for instance, Frank, King &
describe the solar wind). The pressure in the wind is Raine, Accretion Power in Astrophysics; and Dyson &
dropping with radius (because the density is dropping, Williams, Physics of the Interstellar Medium.
right?). When the wind pressure is close to the am-
bient pressure, the wind must slow down. In fact, we
expect a two-shock structure. At this outer boundary
of the wind, we expect some sort of shock transition,
since the wind is supersonic. Past this shock, the hot,
shocked wind-gas will expand into the ISM (at about
its own sound speed, to start); as long as this expansion
is supersonic relative to the ISM, the expanding hot gas
will drive a “snowplowed” shell of ISM, and a second
shock, out into the ISM.
A cartoon of this region, at some point in time,
would be that in Figure 8.8. Let region “a” be the wind;
49
9. SHOCKS IN FLUID FLOW is the basic momentum flux. One useful consequence
of this (resembling Bernoulli’s equation) comes from
Next, assume we have formed a shock, and consider dotting (again) it with n̂:
how it affects the fluid moving through it. We idealize
the shock as an infintesimally thin discontinuity in the p + ρ(v · n̂)2 = 0
££ ¤¤
(9.6)
flow. This allows us to use conservation laws to deter-
mine the “jumps” in macroscopic quantities across the We can also dot (9.5) with t̂, and use (9.4), to show
shock (but of course it does not allow us to say anything
about the fluid flow within the shock).
££ ¤¤
(v · t̂) = 0 (9.7)
A. Jump conditions that is, the component of velocity in the shock plane
To start, recall our three basic equations – conservation does not change.
of mass, momentum, energy – written in conservative • Energy flux (recall ρe + p = γp/(γ − 1):)
form. Here, we ignore body forces, viscosity, and heat- ··µ ¶ ¸¸
ing/cooling losses. The basic equations, then, are: 1 2 γ
ρv + p v · n̂ = 0 (9.8)
2 γ−1
∂ρ
+ ∇ · (ρv) = 0 (9.1)
∂t B. Normal shocks
(from 1.4);
We first consider normal shocks, for which v k n̂. In
∂ this case, the system above reduces to the following.
(ρv) + ∇ · (ρvv + P) = 0 (9.2) Continuity gives
∂t
from (1.13), recalling P is the pressure tensor; and ρ 1 v1 = ρ 2 v2 (9.9)
µ ¶ · µ ¶¸
∂ 1 1
ρe + ρv 2 + ∇ · v ρe + p + ρv 2 =0 Momentum conservation gives
∂t 2 2
(9.3) ρ1 v12 + p1 = ρ2 v22 + p2 (9.10)
(from 6.28). Comment: we are ignoring viscosity (it
makes this job so much easier ...). You remember that The energy jump condition, factoring out ρv from each
viscosity depends on the second derivative of the ve- side, can be written in the case of an adiabatic shock,
locity field. It follows that viscosity is critical to the
internal structure of the shock; but not to the jump con- γ p1 1 2 γ p2 1 2
+ v = + v (9.11)
ditions. γ − 1 ρ1 2 1 γ − 1 ρ2 2 2
We use these conservation laws to find jump con-
ditions. We work in a frame moving with the shock, Now, these equations can be used to express the three
and assume a steady flow in that frame, so ∂/∂t = 0. post-shock quantities, ρ2 , v2 and p2 , in terms of their
(Other frames require a Galilean transform to this shock pre-shock counterparts. One way to do this is to take
frame). General notation: the incoming (upstream) the density jump as the basic variable: X = ρ2 /ρ1 .
quantities are labelled “1”, the outgoing (downstream) One can then eliminate the ratios p2 /p1 , and cs2 /cs1 ,
quantites are labelled “2”, and [[A]] = A2 − A1 is the from the equations, deriving a quadratic in X:
jump in A across a boundary. Let n̂ be the vector nor-
X 2 2 + M2 (γ − 1) −2X 1 + M2 γ +M2 (γ + 1) = 0
£ ¤ ¡ ¢
mal to the boundary, t̂ be a unit vector tangential to the
boundary, and work in a frame in which the shock is and the positive solution of this gives the interesting so-
stationary. Our conservation laws become, then, lution. Writing out the ratios, we get the jump condi-
• Mass flux: tions:
[[ρv · n̂]] = 0 (9.4) ρ1 γ−1 1 2
= + 2
ρ2 γ+1 M γ+1
which describes continuity of mass flux.
p2 2γM2 − (γ − 1)
• Momentum flux (note P is a trace-only tensor if = (9.12)
we ignore viscosity:) p1 γ+1
v2 γ−1 1 2
= +
[[ρv(v · n̂) + pn̂]] = 0 (9.5) v1 γ+1 M γ+1 2
50
(Jumps in T and cS can also be derived). We can also STRONG SHOCK LIMIT, NORMAL SHOCKS
find a condition for the change in the Mach number:
When M >> 1, these equations simplify, to
C. Oblique shocks
Now, consider oblique shocks; in which the upstream
velocity v makes an angle β < π/2 with the shock sur-
face. We know, from (9.7), that the tangential velocity
component u = v cos β does not change in the shock.
Thus, the normal upstream velocity, w = v sin β, con-
trols the jump conditions. It follows from this that an
oblique shock exists only if
M1 sin β ≥ 1
[[ρu]] = 0 ; [[v]] = 0 ;
p + ρw2 = 0
££ ¤¤
(9.15)
·· ¸¸
1 2
h+ w =0
2
Figure 9.1. Downstream conditions (solutions of 9.12, It is also useful to define the deflection angle: δ =
9.13) as a function of the shock Mach number. Top, for − [[β]] = β − β2 . This is the amount by which the
γ = 1.4 (useful for air); bottom, for γ = 5/3 (useful for a flow bends towards the shock (or away from the shock
monatomic gas, such as ionized hydrogen). The differences normal). Geometry gives us, directly,
in γ are apparent in the
u2 tan β2 tan(β − δ)
= =
u1 tan β tan β
51
mach number,
w
2
(γ − 1)M2 sin2 β + 2
M22 sin2 (β − δ) = (9.18)
w1 δ 2γM2 sin2 β + (1 − γ)
β Note that both M2 > 1 and M2 < 1 are possible.
1. TWO POSSIBLE DEFLECTIONS
[[w]] tan δ
− = (9.16)
v cos2 β (1 + tan β tan δ)
We can use the jump conditions to find the usual Figure 9.3. Oblique shock solutions. The dashed lines
are the strong shock branch, the solid lines are the weak
post-shock quantities. If we start by specifying β, we branch, and the heavy dotted line shows the maximum
can find relations analogous to (9.12 & 9.13) deflection angle δmax . From Kundu figure 15.15. Note, σ in
this figure is β in our notation.
ρ2 (γ + 1)M2 sin2 β
= This contains important results. First, the incoming
ρ1 (γ − 1)M2 sin2 β + 2
angle, β, corresponding to a given deflection, δ, is dou-
p2 2γM2 sin2 β γ − 1
= − (9.17) ble valued.1 That is: for any specified deflection angle,
p1 γ+1 γ+1
u1 (γ + 1)M sin2 β
2
= 1
This analysis gives us analytic expressions, (9.17, 9.19 and 9.20),
u2 (γ − 1)M2 sin2 β + 2 for the post-shock conditions as functions of M and β. Al-
ternatively, in shock tables, one may find M and δ specified;
(remembering that M refers to the upstream value, then there are two sets of post-shock conditions, corresponding
M1 .) We can also find an expression for the post-shock to strong and weak shocks.
52
there are two possible shock angles (relative to the in- 2. HIGH MACH NUMBER LIMIT
coming flow). A given δ can come from (a) a strong
shock, corresponding to high β values and giving sub- Finally, as one would expect, oblique shocks also have
sonic postshock flow; and also a weak shock, from low a M ≫ 1, “strong-shock” limit. The limits are just
β values and, generally, supersonic postshock flow. The the same as for a normal shock, (9.14), if the normal
flow will generally find the weak shock solution if it velocity w is substituted for the full velocity v. The
can. strong-shock limit of the deflection angle is still double-
valued, and can be found from the high-M solutions in
strong
Figure 9.3.
shock
M D. The Weak Shock Limit
1
weak
shock The weak shock limit describes shocks at oblique an-
gles (close to the Mach angle) to the upstream flow, and
θ thus shocks in which the velocity and pressure jumps
Figure 9.4. Weak and strong shocks. The geometry (θ) of across the shock are small. These “mach wave” shocks
the wall forces the flow to bend through θ; the mach number are important in bending supersonic flows.
M1 (and flow properties) determine whether forced bend For weak shocks, we can treat [[p]] /p1 as a small
will happen through a weak shock, a strong shock, or parameter. (Notation: Faber calls this σ, and does not
neither.
necessarily assume it is small. I’m following Thomp-
son here, in the weak-shock case.) To avoid notation
A second important result from this analysis is that confusion, here we let
there is a maximum deflection angle possible for any [[p]]
incoming flow. This connects to the the existence of Π=
γp1
detached shocks – that’s what happens when the flow
can’t bend enough at an obstacle to form an attached (the γ factor simplifies the algebra later), and we as-
shock. Figure 9.5 illustrates this. sume Π ≪ 1. With this, we can expand the jump con-
ditions in Π, as
attached [[ρ]] [[u]]
≃Π; ≃ −Π
shock ρ2 cs1
M [[cs ]] γ−1
1 ≃ Π (9.21)
cs1 2
γ+1 γ+1
M1n ≃ 1 + Π; M2n ≃ 1 − Π
4 4
Thus, both the incoming and downstream Mach num-
bers are close to unity. Seen in the lab frame, the shock
is propogating at nearly the sound speed (as one would
detached expect for a low-amplitude wave). We can also find the
shock deflection and the shock angle, for weak oblique shocks
M
1 ¡ 2 ¢1/2
M −1
δ≃ 2
Π
µ M ¶ (9.22)
γ+1 Π
β − µM ≃ 1 −
4 (M2 − 1)
Thus, the flow deflection is small, and the shock angle
Figure 9.5. Attached vs. detached shocks. Top, the (relative to the incoming flow) becomes the Mach an-
opening angle of the wedge is less than the maximum
deflection angle of the flow; an attached shock forms. gle, in the weak shock limit.
Bottom, the opening angle of the wedge is larger than the
PRANDTL - MEYER FUNCTION
maximum possible deflection of the flow; a detached shock
forms.
The relation between the deflection, δ, and the Mach
number, M, turns out to be a perfect differential (as we
53
will show). The integral form of this function will be This can be written in terms of M: using
useful later.
Rearranging the weak-shock relations, we find to dM dw dcs cs dcs
lowest order in [[u]], = = ; wdw + =0
M w cs Γ−1
¡ 2 ¢1/2
M −1 [[u]]
δ≃− (9.23) where the factor Γ = (γ + 1)/2 for a perfect gas. This
M 2 cs1 gives
This is useful now. Since only the normal velocity
changes at the shock, we have w12 − w22 = u21 − u22 ; dw dM/M
and to first order in the jumps, w1 [[w]] ≃ u1 [[u]]. Thus, = (9.26)
w 1 + (Γ − 1)M2
with u1 /v ≃ sin µM = 1/M,
In the limit of infinitesimal strength, the oblique shock This can be integrated formally. Taking the constant of
becomes a Mach wave, and with δ → dδ, [[w]] → dw, integration so that δ(1) = 0, the function describing the
¢1/2 dw deflection angle, in the weak shock limit, as a function
dδ ≃ ± M2 − 1
¡
(9.25) of total upstream mach number is
w
References
I’m taking this mostly from Thompson, Kundu and
Faber.
54
to downstream, and convert ordered to random energy, Thus, this function F is simply a constant times the lo-
just as interparticle collisions do in a shock in a neutral cal sound speed. Now, the pressure derivatives can be
gas. The width of a collisionless shock is typically a written in terms of F , as
fundamental plasma scale, such as the ion gyroradius or
the length derived from the plasma frequency, ∼ cs /ωp . ∂F ∂F ∂p 1 ∂p
= =
∂t ∂p ∂t ρcs ∂t
B. The Method of Characteristics ∂F ∂F ∂p 1 ∂p
= =
In chapter 7 we introduced weak waves – waves of such ∂x ∂p ∂x ρcs ∂x
small amplitude that the local conditions (density, pres-
so that (10.5) can be written as
sure) are not significantly disturbed. These waves prop-
agate at a constant speed, co , thus their loci in an (x, t) ·
∂ ∂
¸
diagram are straight lines. ± (v ± cs ) (v ± F ) = 0 (10.8)
∂t ∂x
Our goal in this section is to extend this approach to
the general case, rather than limiting ourselves to small- At this point, each equation (the “+” and “-” ones) con-
amplitude or subsonic flows. To do this, we first find the tain only one derivative operator:
equations that describe the characteristics. Start with
the 1D momentum and continuity equations (in a uni- D+ ∂ ∂
= + (v + cs ) ;
form channel, no A terms, and with no external body Dt ∂t ∂x (10.9)
forces): D− ∂ ∂
= − (v − cs )
∂v ∂v 1 ∂p Dt ∂t ∂x
+v + =0 (10.2)
∂t ∂x ρ ∂x so that our two equations are, finally,
and
D+ D−
∂ρ ∂ρ ∂v (v + F ) = 0 ; (v − F ) = 0 (10.10)
+v +ρ =0 (10.3) Dt Dt
∂t ∂x ∂x
These equations turn out to be quite useful. We must
But this latter can be written, using ρ = ρ(p) and dp = interpret them physically, as follows. From (10.9) and
c2s dρ, as (10.10), we see that the quantities
1 ∂p v ∂p ∂v
+ + cs =0 (10.4) 2
ρcs ∂t ρcs ∂x ∂x J+ = v + F = v + cs
γ−1
(10.11)
Now, these last two can be added and subracted to form − 2
J =v−F =v− cs
two new equations, γ−1
µ ¶
∂v 1 ∂p ∂v 1 ∂p are constant along the lines described (and labelled) by
± + (v ± cs ) ± = 0 (10.5)
∂t ρcs ∂t ∂x ρcs ∂x
dx
As Thompson notes, it appears doubtful that this is any C+ : = v + cs ;
dt (10.12)
simplification, but it is nonetheless. Define a new funci- dx
−
ton F = F (p, s) (where s is the specific entropy), C : = v − cs
dt
which satisfies
Z p These lines are called characteristic lines. They de-
dp
F = (10.6) scribe the paths of signals (sound waves) travelling for-
po ρcs ward or backwards in a flow, at the local sound speed.
where po is some (useful) reference state, and the inte- This, the quantity v +F is constant for an observer trav-
gration is carried out at constant entropy. For an adia- elling at velocity v + cs , or on any plane perpendicular
batic perfect gas, the F integral can be done straight- to the x axis which moves at v + cs (that is, with a for-
forwardly: ward, or positive, sound wave in a flow with velocity v.
this is the plus characteristic, or the C + wave). Simi-
dp dρ 2dcs 2
Z Z Z
larly, the quantity v − F is is constant for an observer
F = = cs = = cs
ρcs ρ γ−1 γ−1 travelling at velocity v − cs (that is, a backward or neg-
(10.7) ative sound wave; the minus characteristic or the C −
56
Piston advance can be treated in the same way, tak- ary region ahead of the wave. On these characteristics,
ing dX/dt > 0. The difference is that a shock forms J − = Jo− = constant everywhere, so that
ahead of the piston. The analysis above carries through 2 2
just the same – we again note that v+cs = cs,o +dX/dt J− = v − cs = − cs,o ; (10.15)
is the inverse slope of the C + lines. But now, with γ−1 γ−1
dX/dt > 0, the slope decreases moving along the path and thus
of the accelerating piston: the characteristic lines must 2
cross. As we will see next, a shock forms at the point v= (cs − cs,o ) (10.16)
γ−1
where they first intersect.
Thus, fluid within compression regions (cs > cs,o )
moves in the +x direction, the direction of travel of
the wave. In addition, any portion of the C + wave must
move with speed
γ+1 2
v + cs = cs − cs,o (10.17)
γ−1 γ−1
Both of these connections verify the monotonic in-
crease of wave speed with wave amplitude (cs − cs,o )
which gives rise to the distortion of the wave into a
shock. We can visualize this from characteristics also.
The slope of each characteristic, in a (t, x) diagram,
is the inverse of the local sound speed – so regions
Figure 10.6. The geometry and characteristics of higher overdensity have higher cs , and thus lower
associated with smooth piston advance. From Thompson slopes. Thus, characteristics from some region will in-
figure 8.21. tersect at some later time – shown as point s in the fig-
ure.
2. CONNECTION TO SHOCK FORMATION
3. TRAFFIC SHOCKS describes the car density at t = 0. But each point in the
(x, t) plane can be mapped back, along its local char-
Finally, let’s consider a different approach: characteris- acteristic, to a starting point xo . Thus, if things are
tics and shocks in traffic flow. single-valued, the solution can be found directly from
ρ(x, t) = f [xo (x, t)]. (If things are multivalued, that
is characteristics from more than one xo pass through a
given (x, t), then we find that a shock develops. Exam-
ples of both of these are given in the homework.
∂ρ ∂q
+ =0 (10.18)
∂t ∂x
But this is, of course, just the continuity equation. Now,
traffic-flow people simplify this by assuming the flow
can be taken as a function only of the density: q =
q(ρ).1 This allows the basic equation to be written
∂ρ ∂ρ
+ c(ρ) =0 (10.19)
∂t ∂x
if c(ρ) = ∂q/∂ρ.
Now, we know how to treat the equation (10.19):
the solution ρ is constant along the characteristic lines
dx
= c(ρ) (10.20)
dt
Why? Consider the case when the “total derivative” is
zero:
Dρ ∂ρ dx ∂ρ
= + =0
Dt ∂t dt dx
But this is just equation (10.18) if the dx/dt term is
given by (10.20). Thus, (10.19) shows that ρ is constant
along the “path” in (x, t) space given by (10.20). That
is, the lines on which ρ is constant, are simple waves,
straight lines in the (x, t) plane.
Typically, the problem is set up in terms of initial
conditions: ρ(x, 0) = f (x), that is some function f (x)
1
A typical choice assumes the velocity is v(ρ) = V (1 − ρ/ρo ),
ρ < ρo . That is, a linear relation between density and speed, up
to some maximum density ρo at which gridlock sets in and no
further motion occurs.
59
11. TWO-DIMENSIONAL STEADY FLOW But now, we note a significant difference for M > 1
and M < 1 flow. For subsonic flow, M < ∞, (11.5)
We have studied transonic, one-dimensional flows, both is essentially Laplace’s equation, (3.3); so we can solve
time-steady and not. We now move to 2D, steady flows. this linear problem with our potential flow techniques
There are two mathematical approaches that can be use- of chapter 2. We also note that such flows are totally
ful: linearized potential flow, and characteristics. governed by their boundary conditions.1
A. The Nature of Steady, two-dimensional flow. However, if M > 1, the nature of equation (11.5)
changes. In particular, it becomes a wave equation; and
As we have noted already, subsonic and supersonic has as the general solution,
flows behave very differently. A formal way to demon-
strate this is through an extension of potential-flow φ(x, y) = F (x − βy) + G(x + βy) (11.6)
techniques to a general, compressible flow. Consider a
steady, irrotationa, inviscidl flow. Because ∇ × v = 0, where β 2 = M2 − 1, and F , G are any function. This
we choose a velocity potential, v = ∇φ. The continu- is getting back into the regime of characteristics – the
ity and Euler equations for steady flow are flow is determined only by those regions nearby from
1 which information can propagate.
∇ · (ρv) = 0 ; (v · ∇)v + ∇p = 0 (11.1)
ρ We can restate this using some results from formal
PDE theory. Equation (11.5) is a second order¡ PDE; its¢
Now, we note the useful fact, still for steady flow nature depends depends on the sign of the 1 − M2
(whence?): term. If M < 1, the equation is elliptic – and does not
1 have any real characteristics. If M > 1, the equation is
v · ∇p = −c2s ∇ · v hyperbolic, and does have real characteristics. We can
ρ
exploit those characteristics in the solution of steady,
We can combine these results to get a second order two-dimensional problems; which is the point of this
equation in the velocity potential: chapter.
1 A final note: this breakdown into elliptic vs. hy-
∇φ · ∇(∇φ)2 = c2s ∇2 φ
£ ¤
(11.2) perbolic does not depend on the linearization; Schreier
2
demonstrates that the full potential system, (11.2) and
To be more explicit, we can pick Cartesian coordinates, (11.3) also changes type at M = ∞.
and write
∂2φ 2 2 ∂2φ 2 2 ∂φ ∂φ ∂ 2 φ B. Signal Propagation in Flows
2
(cs −φx )+ 2 (cs −φy ) = 2 (11.3)
∂x ∂y ∂x ∂y ∂x∂y If you don’t care for such a mathematical approach, we
These equations (11.2) and (11.3) are general, but not can try a more physical one. Consider the Mach con-
very useful as they are horribly nonlinear in φ. struction (due to Ernst Mach in an 1887 paper on super-
We can simplify it, however, by linearizing. In this sonic projectiles). Let some object emit a steady signal
section I describe one particular version of this: two- while moving at some speed v through a fluid. (This
dimensional Cartesian flow. That is, consider a simple, signal might, for instance, be the simple fact that the ob-
unperturbed flow – for instance, a uniform flow U in ject is moving and thus disturbing the local fluid). The
the x̂ direction. Add a thin body, such as an airfoil, signal propagates at cs relative to the fluid. If the object
nearly aligned with the x-axis. This will perturb the is moving subsonically, the signals “converge” ahead
flow slightly; we can specify the velocity potential to of the object, and “diverge” behind it: this leads to the
describe the perturbed flow. Going through the alge- familiar Doppler effect. Note, all of space eventually
bra (check Thompson, or Schreier, for the details), the receives the signals. The situation is different however,
linearized potential equation becomes if the object is moving supersonically. In this case the
spherical wavefronts emitted as the object moves, de-
¢ ∂2φ 2
2∂ φ fine an (upstream) cone of influence: fluid inside this
U 2 − c2s
¡
= c s (11.4) cone can receive information about the motion, while
∂x2 ∂y 2
Or, in terms of the Mach number,
1
Remember all those problems you solved in electrostatics, for
¢ ∂2φ ∂2φ the electric potential, where the boundary conditions determined
1 − M2
¡
+ 2 =0 (11.5)
∂x2 ∂y everything?
60
fluid outside cannot. The opening angle of this cone is corner. Thus then tend to lean “forward”. Minus char-
acteristics come from ahead of a given point in the flow;
cs 1 at any point the intersection of plus and minus char-
sin µM = = (11.7)
v M acteristics determines θ and M (the latter through the
value of δ(M)). In particular, the flow angle and speed
which is called the Mach angle. The surface bounding
are constant along any given characteristic. Thus, the
the range of influence of the object is called the Mach
flow must follow the bending of the Mach lines. This
surface or characteristic surface; in two dimensions it
geometry forces the streamlines to diverge, so the flow
becomes the Mach lines, or simply characteristics.
must expand. Expanding a supersonic flow accelerates
it – just as in chapter 7.
By the way: where does the energy for the acceler-
ation come from?
C. How Does Supersonic Flow Turn a Corner? Figure 11.2. Flow turning corners, with streamlines and
plus (forward) characteristics (Mach lines) shown. The left
Now, we illustrate the use of characteristics in steady, diagram shows flow into a bending wall. The streamlines
supersonic flow with one example: flow around a cor- come closer together as they follow the wall around, so the
ner. flow must compress and decelerate. A shock forms where
the characteristics intersect. The right picture shows the
Before we hit the math, however, let’s think about opposite case, flow around a corner. The characteristics
the physics. Consider supersonic flow along a wall, and diverge, as does the flow; this expansion results in an
let the wall have a corner. Remember that information acceleration. From Kundu figure 15.18
about the corner cannot propagate very effectively up-
stream, so the incoming fluid cannot “know” very easily
about the corner. Just how, then, does the flow turn the
corner? Now, consider the converse case: flow “into” a wall
To be specific, think about flow around a wall which which bends, as in Figure 11.1a. The geometry forces
turns a corner away from the flow, as in Figure 11.2b. the flow to be compressive, and thus (still being super-
Since we’re working with supersonic flow, we can think sonic) it must slow down. Here, however, the corner-
back to the channel flow of Chapter 7, in particular flow turning takes place through a set of shocks; it cannot
into an expanding channel. The flow must expand, and be smooth. We can see this by drawing in characteris-
being supersonic it must also speed up. This can be tics again. The nature and angle of the characteristics
seen in terms of characteristics (also called Mach lines follows the same rules as above; the angle of the wall
in this context). In the next section we will show that sets the angles of the Mach lines which arise from it.
the Riemann invariants in this situation are given in But now, it’s apparent that these lines must intersect: a
terms of two angles: θ, the angle of the flow relative shock will form. (Formally: think about characteristics
to some axis (say the wall before it turns), and δ(M, intersecting. If “adjacent” lines carry different values
the Prandtl-Meyer angle defined in chapter 9. The in- of θ ± δ, the values of θ and/or δ will not be uniquely
variants are θ ± δ(M). Our boundary (“initial”) condi- defined at the intersection. That is not good. The fluid
tion must be along the wall: the flow (streamline) must will respond by setting up a shock – an infinitely thin
parallel the wall, at the wall. Thus, the characteristics (well, nearly) jump between “front” and “back” prop-
must be straight lines which originate from the wall. erties.
Now, the plus (forward) characteristics have values of By the way: when the flow decelerates, where does
θ which go more and more negative, moving along the the energy go?
61
But, ds/dm+ = cos µM and dn/∂m+ = sin µM ; so everywhere. The other characteristics have
that
µ ¶ θ − δ(M) − 2θ − δ1 on m−
dF ∂F ∂F
= cos µM + tan µM
dm+ ∂s ∂n where δ1 = δ(M1 ). From this, we get
A similar result obtains for derivatives along m− .
Thus, the two characteristic equations for steady, two- on m− : θ = constant and
dimensional flow can be written, (11.17)
δ(M) = δ(M)1 = θ = constant
d sin µM sin θ
(θ + δ(M)) = +
dm+ r (11.15) Thus, each m− characteristic is a straight line, with
d sin µM sin θ
(θ − δ(M)) = − constant δ(M) and θ values everywhere along it. There
dm− r will then be a first intersting m− characteristic, on
Now, there are more general ways to treat characteristic which θ = 0◦ and δ(M) = δ1 = δ(M1 ); and a last
problms, if the derivative of some function along a char- intersting one, on which θ2 = −∆, and δ2 = δ1 + ∆;
acteristic is not zero.2 For now, let’s switch back from the final M value can be found from this, through the
cylindrical to plane geometry: that means take r → ∞ Prandtl-Meyer function.
in (11.15). This gives
θ + δ(M) = constant on m+
(11.16)
θ − δ(M) = constant on m − References
Mostly from Schreier and Thompson, here.
and these equations are good analogs of (10.10). Thus,
we have our characteristics (m± ) and our invariants
[θ ± δ(M)].
2
We didn’t bother with this in chapter 10, feel free to bug a math-
ematician if you’re interested.
63
12. SIMILARITY SOLUTIONS We see from this that the only constant dimensional
quantities which enter the problem are E and ρo ; these
Similarity solutions can be a powerful method when along with r and t consitutue the independent variables
they apply. Consider an example: a point explosion of the problem. (The dependent variables are v, p, ρ).
going off in a uniform atmosphere (such as a supernova
We want to put the basic equations into a dimension-
explosion in a uniform interstellar medium). We expect
less form. The standard approach in similarity solutions
a spherical shell will be driven out into the medium.
is to search for a dimensionless variable, η, which rules
How does this shell expand with time (what is the law
the problem. In general we want
for its size, R(t)?), and what is the gas structure within
the shell? η = (constant)rλ t−µ (12.2)
We can skip ahead if we realize that there are only
two scaling parameters in the problem: the total energy with λ, µ and the constant to be determined from the
E, and the ambient density ρo . We expect the shell ra- specific problem. Next, we set up scalings of the other
dius must involve only these two parameters, and the physical quantities:
time t since the explosion. If we further guess that the
expansion goes as a power law right (this turns out to v = rt−1 ṽ(η) ; p = r−1 t−2 p̃(η) ;ρ = r−3 ρ̃(η)
be right; if we had made a bad guess our later analy- (12.3)
sis would show up the inconsistency), there is only one Derivatives of all the quantities must also be trans-
way to form a distance out of energy, density and time. formed. For instance,
That is, R(t) ∝ (E/ρ)1/5 t2/5 . This must be the form
of the dymamical law for the outer shell. We can skip ∂ µη d ∂ λη d
→− ; →− (12.4)
further ahead: the internal structure of the expanding ∂t t dη ∂r r dη
bubble must also be described in terms of this scaling
law. For instance, r must appear only in combination (Think: why?). When this is done, the PDE’s from the
as r/R(t); velocity only as v/Ṙ(t); and so on. basic equations are transformed into ODE’s. If there
are no independent physical scales in the problem (for
Putting this scaling into the basic dynamical equa- instance, does the ambient density have a length scale?
tions looks rather horrible, but actually leads to a very If so then we can’t solve the problem this way), then
useful simplification of the analysis. This is best done these ODE’s can be solved once and for all time, giving
by example, and I present three in what follows. You us the solution to the problem.
will see that the art of similarity solutions lies in under-
standing the problem, physically, beforehand, so that Returning to the blast wave problem, the arguments
you can make a good guess as to the scaling laws and in the introduction give us the dimensionless combina-
similarity variables. As often in solving physics prob- tion of the basic variables:
lems, hindsight is a great help. ³ ρ ´1/5 r
o
η= (12.5)
A. Blast Waves: the Sedov-Taylor Solution E t2/5
and this is the similarity variable for this problem. The
We start with our initial example, a spherical blast wave
outer shock must then be labelled by some value of η,
resulting from a point explosion.
say ηo ; so that the shock position as a function of time
Dump a quantity E of energy, instantaneously, into is
a perfect-gas atmosphere with density ρo . The energy
µ ¶1/5
appears in the gas as kinetic and internal energy, with a E
strong psherical shock propagating away from the ori- Rs (t) = ηo t2/5 (12.6)
ρo
gin. Let Vs is the shock velocity, and ρs , ps and vs are
the conditions just behind the shock. The jump condi- That is, we expect the radius of the outer shell to be
tions (9.12) for a strong shock become close to our dimensional analysis guess, (12.6), but we
can’t be sure of the exact coefficient yet. From (12.6),
2
ps = ρo Vs2 the shell/shock velocity it
γ+1
γ+1 2 Rs 2
µ
E
¶1/5
ρs = ρo (12.1) Vs (t) = = ηo t−3/5 (12.7)
γ−1 5 t 5 ρo
2
vs = Vs
γ+1 We will find that ηo ≃ 1.
64
We can go further and determine the gas structure at the shock front. The gas flow behind the shock is, as
inside the outer shock. We work with dimensionless usual, described by the basic three equations:
forms of the basic variables:
5(γ + 4) u
v̂ =
4 r/t ∂v ∂v 1 ∂p
+v + =0
25(γ + 1) p ∂t ∂r ρ ∂r
p̂ = (12.8)
8 ρo r2 /t5 ∂ρ ∂v ∂ρ 2ρv
+ρ +v + =0 (12.10)
γ−1 ρ ∂t µ ∂r ∂r¶ r
ρ̂ = ∂ ∂ p
γ + 1 ρo +v ln γ = 0
∂t ∂r ρ
and each of these are only functions of η. (Check back
with equations 12.3). The boundary conditions of the
problem are
After a decent bit of algebra, these can be written in
v̂(ηo ) = 1 ; p̂(ηo ) = 1 ; ρ̂(ηo ) = 1 . (12.9) terms of the “hatted” variables:
dv̂ dp̂ 1
ρ̂ (2v̂ − γ − 1) + (γ − 1) = [ρ̂v̂(5γ + 5 − 4v̂) − 4(γ − 1)p̂]
dη dη 2η
µ ¶
γ + 1 1 dρ̂ dv̂ 3v̂
v̂ − + =− (12.11)
2 ρ̂ dη dη η
d p̂ 1 5(γ + 1) − 4v̂
ln γ =
dη ρ̂ η 2v̂ − (γ + 1)
sion gets denser and becomes a dense, thin shell. Past (the corner) to the point (r, ψ) under consideration. The
this point the expansion is governed by momentum con- irrotational equation becomes
servation rather than energy conservation. . . which is
another topic that we won’t address here. dθ 1 dv
+ cot(ψ − θ) = 0 (12.15)
dψ v dψ
B. Prandtl-Meyer flow, revisited
and, using the useful connection,
In Chapter 11 we slogged through the characteristic-
based solution to this problem; now, let’s repeat the
µ ¶ µ ¶
dρ dh ∂ρ dv ∂ρ ρv dv
analysis using similarity methods. Here, we do not have = =v =− 2
dψ dψ ∂h s dψ ∂h s c dψ
an easy identification with a physical scale (there is no (12.16)
“radius of the shell” to consider). Instead, we guess that (where we used the chain rule; energy equation; and
our fundamental variables must be angles (or, equiva- Gibbs relation, dh = T ds + vdp, to find (∂ρ/∂h|s =
lently, the two cartesian coordinates (x, y) must only ρ/c2 ). With this, the continuity equation becomes
appear as the ratio x/y). Figure 12.2 has the geometry,
in today’s notation. dθ ¡ 2 ¢ 1 dv
cot(ψ − θ) + M −1 =0 (12.17)
dψ v dψ
C. Gravitational Collapse: the Shu Solution These are the heart of this problem (compare the sets
12.11, or 12.15 and 12.17). Shu presents numerical so-
A nice third example of this is the gravitational collapse lutions, which are illustrated in Figure 12.3.
of a self-gravitating cloud – such as a protostar. Frank
Shu has developed this area extensively, and I take this
(quite basic) analysis from his book.
What is the physical picture? Let’s start with a
static, isothermal sphere; take a = (kT /m)1/2 as the
sound speed. For this argument, let’s consider the sim-
ple version with ρ(r) ∝ 1/r2 , and ignore the unphys-
ical divergeance at the origin. Imagine some perturba-
tion, at t = 0, that makes the inside of the clour (pro-
tostar) slightly overdense. This will lead to an “inside-
out” collapse of the protostar. As the expansion starts,
for r ∼ 0, we do not expect the outer regions to know
about the expansion at all; an expansion wave will move
outwards, at a speed a, as the collapse proceeds. This
suggests a similarity variable, x = r/at. If we as-
sume there is no important core (ignore the core dis-
cussed in (6.8) and following), then there is no physical
scale length in the problem, and we can apply similarity
methods. Figure 12.3. Self-similar solutions for Shu’s collapse
Pick our basic equations, then. It is convenient to model. The solid line describes the density; the dotted line
the infall velocity; and the dashed line, the mass function.
rewrite the continuity equation (1.4) in terms of mass From table 18.1 of Shu.
shells, taking M (r) as the mass within radius r:
∂M ∂M ∂M Applications and extensions of this may be dis-
= 4πr2 ρ ; +v =0 (12.23) cussed in class or in the homework (depending on
∂r ∂t ∂r
Jean’s mood and time).
For isothermal flow, the force equation is
∂v ∂v a2 ∂ρ GM
+v =− − (12.24)
∂t ∂r ρ ∂r r
References
We convert to similarity variables:
The Sedov solution is the classical example here; I
α(x) followed Thompson and Shore, but it’s also in many
ρ(r, t) = other references. I took the Prandtl-Meyer solution
4πGt2
a3 t from Thompson, and the Gravitational Collapse from
M (r, t) = m(x) (12.25) Shu’s book.
G
v(r, t) = av(x)
Alternatively, using vector algebra and defining the The first term describes intrinsic changes in B, while
quantity η = c2 /4πσ, this becomes the second describes advection of the field through the
boundary. Now, Stokes’ theorem rewrites this as
∂B
= ∇ × (v × B) + η∇2 B
Z µ ¶
(13.7) dΦB ∂B
∂t = − ∇ × (v × B) · dS (13.12)
dt S ∂t
which is often called the induction equation. Two com- Now, (13.7) allows us to write this, generally, as
ments on this equation.
dΦB
Z
• Compare this to equation (4.8), which describes = η∇2 B · dS (13.13)
the evolution of vorticity, in the barotropic limit. They dt S
are just the same: and so we should expect the magnetic Thus, as η → 0 (the perfect conductivity limit, an
field and the vorticity to behave similarly in hydrody- “ideal” plasma), ΦB becomes a constant of the mo-
namic situations. One example is flux freezing, which tion. That is, the magnetic flux through a loop which
we derive immediately below: it’s the MHD analog of is “attached to the plasma”, stays constant as that loop
Kelvin’s theorem. Another example is the existence of is stretched (or compressed) by the flow field. The
magnetic flux ropes: a field subject to fluid forces tends mean magnetic field in that loop therefore decreases (or
to bunch into linear, high-field regions (think of a tor- grows), in proportion to the change of loop area. This
nado or a vortex line). is Flux Freezing.2
• Equation (13.7) has two important limits, as fol- Magnetic Field Lines? (A discussion from
lows. When we worked with fluids, we noted that the Schmidt...) The concept of a magnetic line of
behavior of a system is sensitive to the ratio of the ad- force is an abstraction. IN general no identity can
vective to dissipative terms. We expressed this ratio in be attached to these lines (they cannot be labelled
terms of the Reynolds number. The induction equation in a varying field), nor can we speak of “motion”
admits a similar separation of effects. We take the mag- of field lines. In a perfect conductor, however, the
netic Reynolds number concept of field lines becomes meaningful, due to
flux freezing.
LV Consider a material line in the fliud (say a chain
Rm = (13.8) of labelled droplets, or particles painted pink), de-
η
fined by intersecting material surfaaces. Choose
where L, V are characteristic length and velocity scales these surfaces everywhere tangential to B at t =
of the problem. This quantity again measures the radio 0. The flux through both surfaces is therefore zero
of inertial/advection, to dissipative, terms. to start, and their intersection defines a field line at
that point. Flux freezing guarantees that these sur-
1. IDEAL LIMIT: FLUX FREEZING faces continue to ΦB = 0 at any later time. Thus,
their intersection continues to define a field line,
An imporant limit in MHD is the limit when Rm ≫ 1, in fact the same field line – it has become identi-
called ideal MHD . This is the case of a highly con- fiable; labelling the material (painting it pink) has
ducting plasma. Here, the induction equationk (13.7), labelled the field line, and the local fluic velocity
reduces to v(x, t) is also the velocity of that section of the
field line. The field line is attached to – “frozen
∂B into” – the fluid.
≃ ∇ × (v × B) (13.9)
∂t 2. RESISTIVE LIMIT: FLUX ANNIHILATION
while Ohm’s law becomes In a fluid with finite conductivity, flux freezing no
v longer holds. We can explore this by going to the other
E+ ×B≃0 (13.10) limiting case, when Rm ≪ 1. This is diffusive limit.
c
If we simply ignore the advection term, equation (13.7)
Consider, now, a closed curve (C) bounding a surface becomes
(S) which is moving withR the plasma. The magnetic ∂B
flux through S is ΦB = B · dS; and its rate of change = η∇2 B (13.14)
∂t
can be written
dΦB ∂B
Z I
2
Compare Kelvin’s Theorem, from chapter 4; you will see that the
= · dS − B · (v × dl) (13.11)
dt S ∂t C
math is identical, as is the nature and interperation of the result.
69
This describes the effect of Ohmic dissipation on the 1. ENERGETICS OF THE E AND B FIELDS
magnetic field; note that it is a standard diffusion equa-
tion. To remind yourself of this...go to Jackson.3 Keep all of
Diffusion of field lines? We know how solutions the terms – the full Maxwell’s equations. Consider the
to (18.4) behave: an initial field will decay on a quantity E · j: it expands as
timescale ∼ L2 /η. Some authors discuss this in ∂E
terms of field line “diffusion” or “slippage” out 4πj · E = cE · (∇ × B) − E ·
of the fluid. Remember that the density of field ∂t
lines is related to the strength of the field; so a But now, use the vector identity ∇·(E×B) = B·(∇×
lower density of field lines, with time, should cor- E) − E · (∇ × B), and Maxwell again, to get
respond to field lines “diffusing” out of the field.
In particular, when η is finite, field lines are no c ∂E ∂B
−j · E = ∇ · (E × B) + E · +B·
longer tied to parcels of the plasma; some authors 4π ∂t ∂t
talk of field lines “moving through” the plasma in Now: we identify the the energy density in the fields, as
dissipative regions. uE = E 2 /8π; and uB = B 2 /8π. And we integrate the
C. Fluid Equations: Lorentz force equation above over some volume V , with a surface S,
giving
The effect of the B field on the force equation is Z Z
d
straightforward. We simply add the Lorentz force to − j · EdV = (uE + uB ) dV
the momentum equation: V V dt
(13.17)
c
Z
Dv j + (E × B) · dS
ρ = −∇p + × B (13.15) S 4π
Dt c
Note, I have ignored viscosity here, as well as any ex- Thus: the first term is the rate of change of field energy
ternal forces (such as gravity). Now: expand out the in V . The second term – involving the Poynting flux,
Lorentz force as S = (c/4π)(E×B) – is the flow of EM energy through
j 1 the surface. The last term, then, is the work done by
×B= (∇ × B) × B the fields on the sources in the volume (and it can have
c 4π (13.16)
1 1 either sign ... thinking about driving vs. dissipation .. I
2
= − ∇B + (B · ∇)B think). This can also be written in differential form:
8π 4π
This is an important breakdown of the Lorentz force; d
(uE + uB ) + ∇ · S = −j · E (13.18)
it demonstrates that the field exerts a magnetic ten- dt
sion and a magnetic pressure on the fluid. The first which looks much like any other conservation equation.
term in (13.16) represents the gradient of a scalar pres-
sure, pB = B 2 /8π. It appears in the momentum 2. ENERGETICS OF THE FLUID
equation parallel to the fluid pressure....you can think
Schmidt has a nice approach: I’ve redone it in cgs,
of trying to compress a magnetic field, parallel to it-
keeping all the terms. To start, dot v into the force
self, with the field resisting the compression (“fight-
equation (13.15):
ing back”). The second term in (13.16) is non-zero
only if the field varies parallel to itself. A simple il- Dv 1
ρv · + v · ∇p = v · j × B
lustration is a curved field line. The curvature means Dt · c (13.19)
there is a current flowing along the field line; the j × B
¸
1 1 ∂E
force points inwards (relative to the curvature). Thus, = v · (∇ × B) − ×B
4π c ∂t
curved field lines “want to straighten out”...Some au-
thors combine both effects by describing magnetic field where I’ve used the full Maxwell on the right. With
lines as “elastic bands within the fluid”, which resist some algebra, the terms on the LHS can be written:
being stretched: either pushed together, or pulled trans- ·
∂v
¸
verse to their length. ρv· + (v · ∇)v
∂t
µ ¶ µ ¶ (13.20)
D. Fluid Equations: Energetics ∂ 1 2 1 2
= ρv + ∇ · ρv v
∂t 2 2
We need two items here: first recall the energetics of
the fields, by themselves; then their effects on the ener-
3
getics of the fluid. Preferably the 1st edition, where he works in in cgs.
70
and just the same as the first part of the course. Now, con-
sider the last term, using vector algebra to reorganize
1 ∂p γ each term:
v · ∇p = + ∇ · (pv) (13.21)
γ − 1 ∂t γ−1
where we’ve used p ∝ ργ in this last; otherwise these
are both general, for compressible fluids. So far, this is
· ¸
1 1 ∂E 1 1 ∂E
v · (∇ × B) − × B = − (v × B) · (∇ × B) + (v × B) · (13.22)
4π c ∂t 4π 4πc ∂t
Now ... use v × B/c = j/σ − E; and E · (∇ × B) = B · ∇ × E − ∇ · (E × B); expand out and collect terms, to get
µ ¶
v 1 2 1 ∂E ∂B c
· (j × B) = j − E· +B· − ∇ · (E × B) (13.23)
c σ 4π ∂t ∂t 4π
This is still fully general; it combines the EM and fluid tween collisions, tcoll , or its inverse the collision rate,
terms in a conservative form... We can now follow the νcoll = 1/tcoll , as follows. Consider a free electron, in
usual MHD approach and drop E 2 in the ∂/∂t term (as a plasma, subjected to an external electric field E. The
Schmidt does). Also note that Parker writes the Poynt- net force on the particle can be estimated,
ing flux, in the MHD case, as S = B × (v × B)/4π =
v⊥ B 2 /4π, explicitly assuming ideal and an induction- ∆p
Fnet ≃ eE − (13.25)
only E field. ∆t
where ∆p/∆t is the mean rate of momentum change
per collision. But if the charges have a net drift velocity
vD , we can estimate ∆p/∆t ∼ me vD /tcoll ; then, in a
References steady state we have Fnet ≃ 0, so that
Once again, much of this basic material is “just from
eE ≃ nvνcoll (13.26)
me”. Possibly useful references might be Priest (who
has a very good MHD introduction), also Woods and and the drift velocity must be vD = eEtcoll /me . Next,
Davidson. I followed Schmidt for the energetics analy- we can use this in the (static) Ohm’s law, to relate the
sis, and cannot for the life of me remember where I first conductivity to the drift velocity:
found the Generalized Ohm’s law arguments (probably
Chen’s plasma book? ... sorry). j = ne evD = σE (13.27)
plasma, as in chapter 1 (check §1.F). Alternatively, the • microscopic. Alternatively, consider single parti-
collisions may be with microturbulence in the plasma cle motion. The particles undergo gyromotion around
... in which case the collision rate usually has to be put B; they undergo an E × B drift; and they suffer col-
in by hand (rather than calculated from first principles). lisions which disrupt these ordered motions. To illus-
trate (I take this from Park), let B = (0, 0, B) and
2. CROSS - FIELD CONDUCTIVITY E = (Ex , 0, Ez ). The equations of motion, in the ab-
Consider a general situation in which E and B exist; sence of collisions, are
there will be Ek and E⊥ components, relative to B. We
can anticipate the three terms in the Generalized Ohm’s dvx
law, as follows: (1) collisional conductivity, giving a m = qEx + qBvy
dt
field-aligned current proportional to Ek , as usual; (2) dvy
quasi-collisional conductivity, giving a cross-field cur- m = qBvx
dt
rent proportionao to E⊥ ; and (3) a transverse current, dvx
perpendicular to both E and B, based on single-particle m = qEz
dt
E × B drift.4 Note that this drift is independent of par-
ticle charge and mass .. but we get a net current be-
cause the collision frequencies are not. There are two with solutions
approaches, macro and micro.
• macroscopic. Now: just consider one species, vyo + qEx
vx = vxo cos Ωt + sin Ωt
still, and extend (??eqn13.26) as qB
vyo + qEx qEx
³ v ´
vy = cos Ωt − vxo sin Ωt − (13.33)
e E + × B + mvνcoll = 0 (13.29) qB qB
c
qEz
Retain the definition j = nev; so that vz = vzo + t
m
mνcoll 1
j + j × B = −eE (13.30)
ne nc (This shows the combination of gyromotion and E ×
Bdrift). Now, include collisions. Let νcoll be the colli-
Because we have cross products, things get complicated
sion frequency, and assume collisions occur randomly
.. choose the form,
in time. The probability of a collision is e−νcoll t n and
j = σo Ek + σ⊥ E⊥ + σH n̂ × E (13.31) the probability that a particle will escape a collision in
time t, t + dt is given by νcoll e−νcoll t dt. Inbetween col-
Put this into (13.30), do the algebra and collect the lisions the trajectories of (13.33) are followed. Thus,
terms5 We find, the drift speeds are
ne2
σo = Z
mνcoll hvx i = A νcoll e−νcoll t (Ex /B) sin Ωtdt
ν2
σ⊥ = σo 2 coll 2 (13.32) Ex νcoll Ω
νcoll + Ω = 2 + Ω2
B νcoll
νcoll Ω
σH = σo 2 Z
νcoll + Ω2 hvy i = A νcoll e−νcoll t (Ex /B)(cos Ωt − 1)dt
where Ω is the gyrofrequency, as usual. These are the Ex Ω2
collisional, Pederson and Hall conductivities. =− 2
B νcoll + Ω2
Z
4
Quick physics here: a single particle undergoes gyromotion hvz i = A νcoll e−νcoll t (Ex /B) sin Ωtdt
around the local magnetic field, at a frequency Ω = eB/mc. If
there is also a perpendicular E field, the center of the gyro-orbit
Ez
= (13.34)
shifts during one orbit. Following the particle’s trajectory, it’s Bνcoll
easy to show that the orbit-center moves at a steady drift speed,
∝ E × B.
if the normalizing factor A−1 = νcoll e−νcoll t dt. Fi-
R
5
Remember the definitions: ωp2 = 4πne2 /m, and Ω = eB/mc;
that’s how the frequencies ωp and Ω come into the expression. nally, use this in the definition of current density, j =
72
· ¸ · ¸
∂ 1 2 1 1 2 1 γ 1 1
ρv + p+ B + ∇ · ρv 2 v + pv + E × B = − j 2 (13.40)
∂t 2 γ−1 2µo 2 γ−1 µo σ
14. SIMPLE MHD EQUILIBRIA dynamics, as the fluid equations do not provide for in-
trinsic static equilibria (unless a gravitational field is
How many ways can you construct a steady magnetic included). The plasma equations do, in principle, pro-
field? When you have done that, can it confine a vide for self-confinement: if the plasma pressure can
plasma? There are several common applications. just balance the Lorentz forces from the fields. Most
commonly, flows are ignored, as is resistivity. The gen-
A. Potential Fields
eral condition for equilibrium is, then,
Now and then we can simplify by considering static
magnetic fields in vacuum – that is fields arising from j
× B = ∇p (14.5)
a current which is confined to some finite spatial re- c
gion. A typical example is the magnetic field of the
sun, above the solar surface. This is, of course, can be solved, subject to the con-
straints
Go back to Maxwell: in a current-free region, the
magnetic field satisfies ∇×B = 0. That means it can be c
∇·B=0; ∇·j=0; j= ∇×B (14.6)
found from a scalar potential, Φm . Because ∇ · B = 0 4π
always, the potential (and field) satisfy
(The second relation holds in steady state, right?). A
∇2 Φm = 0 ; B = −∇Φm (14.1) system which satisfies (14.5) also obeys
so that we have a simple, radial pressure balance. For One possible equilibrium (illustrated in the figure) is
instance, one possible equilibrium (illustrated in the fig-
ure) is Bo r
Bφ (r) = (14.12)
1 + r2 /a2
−r 2 /a2
p(r) = po e
h i (14.10) (exercise for the student: what are the correspond-
2 2
Bz (r) = Bo 1 − βo e−r /a ing pressure and current density profiles?) This type
of pinch can be self-confining.R Note that the cur-
r
with βo = 8πpo /Bo2 . rent within radius r is I(r) = 0 2πjz rdr, and from
Such a system is effectively confined by the box Maxwell the B field is Bφ (r) = 2I(r)/rc. Using these
(which carries the φ-currents). We will see that it is and the pressure balance condition (14.11),
stable; however losses out the ends are severe, so that it Z a
isn’t particularly useful for laboratory confinement. 1 2
prdr = I
2 a
(14.13)
0 4πc
z
3
2 B
z
2.5
0
0 1 2 3
r
2
plasma confines itself by carrying just the right net cur- Figure 14.2. Above, the geometry of a linear pinch: the
rent Io . The basic relation is field is azimuthal and the current is axial. Below, qualitative
à ! solutions for a linear pinch. Following Freidberg figures
d Bφ2 Bφ2 15.4, 15.5
p+ =− (14.11)
dr 8π 4πr
75
3. GENERAL SCREW PINCH The simplest approach here is to pick a simple form
for α, usually constant, and look for analytic solutions.
A combination of these two cases is clearly possible, Note one simple consequence: taking the divergeance
in which both Bz and Bφ are non-zero. The general of (14.16), we have
equilibrium is, then, (14.8). One example, which may
be considered in more detail in the homework, comes
∇ · (αB) = (B · ∇)α = 0 (14.17)
from Freidberg: consider
2I r so that α must be constant along field lines; that is, B
Bφ = (14.14) lies on surfaces of constant α; these are magnetic flux
c r + a2
2
surfaces. If we take α = constant everywhere, in fact,
where I is the plasma current and a is a scale length. (14.16) becomes
Let the pressure and Bz field be related by
¡ 2
∇ + α2 B = 0
¢
(14.18)
Bz2 (r) = Bo2 − λp(r) (14.15)
where Bo is the externally applied field and λ is a con- Now: the mathematically astute among you will rec-
stant, determined experimentally. ognized immediately that this is the vector Helmholtz
equation. (The rest of you – like me – can refer, for
A mixed pinch of this type turns out to be more sta-
instance, to Morse & Feshbach chapter 13). It has well-
ble than the simple Z pinch, due to its axial field. How-
understood, if complicated, solutions. Here, I quote a
ever it retains the problem of losses out the ends. The
couple of useful solutions without deriving them.
next step might be to bend the ends of the screw pinch
together, to form a toroid. The simplest form of that is 1. CYLINDRICAL GEOMETRY
a tokamak; many, many trees have been lost to suppor
the books and papers analyzing this geometry. It turns The system becomes (following Priest)
out still to suffer serious instability problems. You can
go one step further by twisting the cylinder, about its
à !
d Bφ2 + Bz2 Bφ2
axis, before connecting the ends...that gives you a stel- + =0
larator. I have no intention of pursuing the specialized dr 8π 4πr
(14.19)
geometry that these machines require, in these notes. dBz
− = αBφ
dr
C. Force-Free Fields
But this has Bessel function solutions:
Another track considers solutions to the MHD force
equation when both gravity and plasma pressure are
negligible. In this case, we must have Bφ (r) = Bo J1 (αr) ; Bz (r) = Bo Jo (αr) (14.20)
X
Br = Cn n(n + 1)r−3/2 Jn+(1/2) (αr)Pn (cos θ)
n
h i dP
n
X
Bθ = Cn −nr−3/2 Jn+(1/2) (αr) + αr−1/2 Jn−(1/2) (αr) (14.22)
n
dθ
X dPn
Bφ = − Cn αr−1/2 Jn+(1/2) (αr)
n
dθ
These solutions describe nested, toroidal flux surfaces; conditions, to keep things well behaved. For cylindri-
the field lines are again helices on the flux surfaces. cal solutions, a physical boundary at some outer radius
makes sense (say from external pressure or a wall). For
3. NON - LINEAR FIELDS the spherical solution just discussed, we can also im-
A third interesting case has uniformly twisted solutions, pose a wall (as in the lab), or can also assume the α
and does not assume α = constant. In cylindrical ge- constant is non-zero only inside of some R. This re-
ometry, a field line is traced in (φ, z) by quires that the usual magnetic boundary conditions be
met across r = R, and keeps the currents (sources)
rdφ dz confined to a finite region.
= (14.23)
Bφ Bz
D. Gravitational Equilibrium I
Thus, the twist of a field line between 0 and z is
Z z In Chapter 6 we addressed the question of hydrostatic
Bφ
Z
equilibrium: the static support of a fluid, by its own
ϕ(z) = dφ = dz (14.24)
0 rBz pressure, in a gravitational field. When we add a mag-
netic field, the static balance becomes much more com-
Now, a field whose line-twist is independent of r must plex:
satisfy
Bφ j
= br (14.25) ∇p = ρg + ×B (14.28)
Bz c
Requiring that this field also be force-free, that is has Recall the gravitational acceleration can be written
j × B = 0, gives in terms of a potential, g = −∇Φg . The three-
dimensional nature of the Lorentz force, j × B, makes
dBz Bφ d this a much more complicated problem if we treat it
Bz + (rBφ ) = 0 (14.26)
dr r dr fully generally. Let’s avoid that and find some simple
Combining this with b = constant gives this solution approaches.
1. “ CONVECTIVE ” INSTABILITY
This analysis can be extended by considering adi- buoyancy force will overcome the tension, leading to
abatic changes of the parcel density and temperature, instability. Comparing these two forces, we find insta-
as in chapter 6; we may explore this in class or in the bility occurs if
homework.
B2
2.SCALE OF UNSTABLE PERTURBATIONS : g∆ρ > ; λ > 2H (14.37)
4πλ
PARKER I NSTABILITY
We can make the same argument for the magnetic
Now, let the magnetic structure bend as it rises. The bouyant (in)stability of a flat flux sheet (as in the planar
gas is free to flow along the lines, and thus will accu- model of the galaxy, above). Our static solution gives a
mulate in the “valleys”; the “tops” of the bent field lines new scale height:
will thus be lighter than the surrounding gas, and will
rise due to bouyant forces. This instability turns out to p ρ B2
be important for long wavelenghts: only large perturba- = = 2 = e−z/Λ (14.38)
po ρo Bo
tions are unstable.
where Λ = c2s + vA 2 /2 /g, in direct extension of the
¡ ¢
To see this, first consider an isolated flux tube, such
as might give rise to a sunspot. Let the external gas have nonmagnetized scale height from above (check back to
a density scale height H = kB T /mg (refer back to Chapter 6). Let the flux sheet be raised some small ∆z,
chapter 6 for hydrostatic equilibrium with no magnetic with the gas sliding down into the “valleys” again. If the
field). If the magnetic field is confined in the flux tube, perturbation has horizontal scale λ, its effective curva-
and the intial state is in pressure balance, we have the ture radius is, from simple geometry,
internal/external balance,
λ2 2
r≃
B2 16 ∆z
pi + = pe (14.35)
8π The density at the top of the perturbation is now (take
“o” to be at the height of the undisturbed sheet)
µ ¶
−∆z/H ∆z
ρi (∆z) ≃ ρo e ≃ ρo 1 − (14.39)
H
(Note there is no magnetic pressure supporting the gas
inside the bent flux sheet). The change in external den-
sity does know about magnetic pressure however:
µ ¶
−∆z/Λ ∆z
ρ3 (∆z) ≃ ρo e ≃ ρo 1 − (14.40)
Λ
If we again compare bouyancy to the restoring force of
magnetic tension, we find the instability condition for
Figure 14.4. The geometry of a sub-surface flux tube
before it erupts from the sun due to bouyancy, and its
this case:
possible post-eruption state; from Tajima &Shibata figure 16Λ2
3.17. λ2 > (14.41)
(1 + 1/β)2
Assuming the gas inside is at the same temperature,
it must be at lower density than the outside. Thus leads if β = 8πp/B 2 is the ratio of gas to magnetic pressure.
to a bouyancy force,
B2
Fbuoy = g (ρe − ρi ) = g∆ρ = (14.36)
8πH References
Say, now, that the flux tube is bent upwards, locally, The general cylindrical MHD equilibria are dis-
with a radius of curvature λ. If λ is short, magnetic cussed in several books; Priest and Bateman are good
tension will pull the tube back towards its initial posi- sources. Moffatt discusses the force-free solutions. I
tion, giving a stable system. If, however, λ is long, the took the Magnetic Buoyancy discussion from Shu.
79
For our geometry, consider the magnetic flux through value of Ψ at the footpoint of the line (for instance
a surface taken in the (x, z) plane (think of a unit the surface of the sun). This result also shows us that
length in y to make a reasonable mental picture, so surfaces of constant pressure nearly coincide with field
da = dydz → dz). This is just lines. For scales ≪ H(z), the pressure is nearly con-
Z z stant, and hydrostatic balance becomes j × B ≃ c∇p.
∂Ψ Dotting this with B or with j gives
ΦB = dz = Ψ(z) − Ψ(0) (15.6)
o ∂z
B · ∇p ≃ 0 ; j · ∇p ≃ 0 .
Thus: if we take Ψ(0) = 0, then Ψ(z) measures the
magnetic flux contained “within” (below in this case) Thus, the magnetic field and current lie nearly in sur-
the surface of constant Ψ. faces of constant pressure. Deviations occur only on
scales ∼> H(z) (for an example see the homework).
2. APPLY THEM : SOLAR MAGNETIC ARCHES
Fact 3. we can consider the full force equation as a
Now: put this formalism into the magnetostatic equa- system to be solved for the function Ψ(z) (and thus the
tion, (15.1). We note several useful facts. field structure). The gravitational and pressure terms
Fact 1. The associated current is easy to find: can be combined as
j = (c/4π)(∇2 Ψ)ŷ ρ∇Φg + ∇p = −e−Φg /cs ∇q
2
(15.11)
You can find this in your favorite EM book, or derive it
where, in the last step, we have introduced the usual
directly from ∇ × (∇ × A).2 From (15.2) and (15.1),
gravitational potential, given by g = −∇Φg , and have
the components of the Lorentz force are 2
· ¸ defined the function q = peΦg /cs . Combining this with
j 1 ¡ 2 ¢ ∂Ψ ¡ 2 ¢ ∂Ψ (15.8), we can write the full force equation, (15.1), as
×B= ∇ Ψ , 0, ∇ Ψ (15.7)
c 4π ∂x ∂z
1 2
− (∇2 Ψ)∇Ψ = e−Φg /cs ∇q (15.12)
which can also be written, 4π
j 1 ¡ 2 ¢ This balance requires that ∇Ψ and ∇q be parallel to
×B= ∇ Ψ ∇Ψ (15.8)
c 4π one another; so that surfaces of constant Ψ must also
be surfaces of constant q. It follows that we can pick
Details here: The ∇2 Ψ term is a scalar function, which q = q(Ψ), i.e. that q can be written as a function of Ψ.
multiplies the vector ∇Ψ. We have taken an impor- Given this, our basic equation turns into
tant step: we have reduced our system to one scalar
unknown, Ψ, rather than the vector unknowns B and j. 2 dq
Fact 2. What is the pressure distribution? Consider ∇2 Ψ = −4πe−Φg /cs (15.13)
dΨ
some field line, locally B, at a direction θ to the vertical.
The component of (15.1) parallel to B is This is an important result: we have reduced the system
from a general vector equation, (15.1), to a quasilinear
dp dp PDE in Ψ(x, z). If we can pick q(Ψ) (for instance, con-
0=− − ρg cos θ; = −ρg (15.9)
ds dz stant or linear or some other simple function), we can
if ds = dz/ cos θ is the differential distance along the find solutions for Ψ(x, z) by standard PDE methods.
field line. The second equality just converts from ds An example of this will appear in class or in the
to dz, and shows that normal hydrostatic equilibrium homework.
obtains along the field line. This can be written,
Rz
B. Flux Functions II: Grad-Shafranov Equation
− dz/H(z)
p(z) = po (Ψ)e o (15.10)
Another version of flux functions is applied to plasma
if H(z) = kB T (z)/mg is the local scale height.3
I confinement problems which can be reduced to two in-
have noted the explicit dependence of po on the local dependent variables (such as an axisymmetric cylinder
problem). Again, the general equilibrium condition, in
the absence of gravity, is
2
BIG NOTE: ∇2 Ψ = ∂ 2 Ψ/∂x2 + ∂ 2 Ψ/∂z 2 is a scalar, with
simple form in Cartesian. j
3
Compare to the scale height in Chapter 6 for a uniform T . × B = ∇p (15.14)
c
81
Once again, it follows that j · ∇p = 0 and ∇p · B = 0. One function, ψ, is enough if the field has no φ com-
That is, B and j are normal to ∇p; each of B and j lie ponent. For the more general case we need a second
in constant-p surfaces. These are called magnetic sur- scalar function. Consider the total current J, within
faces or flux surfaces; they contain the field lines. Most surface S which is planar and lies perpendicular to the
interesting equilibria consist of sets of nested magnetic z axis (refer back to Figure 15.2 for such a surface):
surfaces. Once again, we will label these surfaces with Z Z 2π Z r
flux functions, ψ.4 We will have ψ = constant in a flux
J= j · da = jz rdφdr (15.17)
surface; and we can write p = p(ψ), so that p(ψ) is S 0 0
constant on the surface given by ψ = constant.
But now, since rjz = ∂(Bφ r)/4π∂r in axisymmetry,
z we can write
c
J = Bφ r (15.18)
2
B This is our second surface function, J(ψ), the axial cur-
rent crossing through the surface S. It pairs with p in
the analysis, both functions being constant on a ψ sur-
face.
r 2. APPLY THEM : THE G - S EQUATION
c ∂2ψ
· µ ¶¸
d 1 ∂ψ c
jφ = − 2
+r = − ∆∗ ψ
4πr ∂z ∂r r ∂r 4π
Figure 15.2. Geometry for the axisymmetric cylindrical (15.19)
problem. The magnetic flux through the circle,
R where we have followed convention and defined
perpendicular to the z-axis, is ΦB = B · da = Bz 2πrdr.
Using (15.16), this becomes ΦB (r) = 2πrAφ = 2πψ(r, z) 1 ∂2ψ
· µ
d 1 ∂ψ
¶¸
∗
(you should derive this!). See the text for more details. ∆ ψ= +r (15.20)
r ∂z 2 ∂r r ∂r
1. DEFINE THE FLUX FUNCTIONS to shorten the notation. From this we can write the gen-
eral B field in terms of two scalars J and ψ:
We need to introduce two scalar functions. The first is µ ¶
the φ- component of the vector potential, which can be 1 ∂ψ 2J 1 ∂ψ
related to ψ, the poloidal flux function B= − , , (15.21)
r ∂z cr r ∂r
ψ(r, φ) = Aφ r (15.15) Similarly, the total current density is
(Figure 15.2 shows the basic geometry for the simpler, 1
µ
1 ∂J c ∗ 1 ∂J
¶
axisymmetric case. As with the Cartesian case, the j= − , ∆ ψ, (15.22)
scalar function ψ(r) labels the flux within a circle of 2π r ∂z 2 r ∂r
radius r. It follows, as before, that the poloidal field
components can be written These surface functions, J and ψ, can be used to
µ ¶ reduce the equilibrium condition, (15.14), to one DE
ψ rather than three. As we noted abive, because both the
Bpol = ∇ × φ̂ φ ⇒ pressure p and the axial current J are surface functions,
r
(15.16) we can write p = p(ψ), and J = J(ψ). Our most
1 ∂ψ 1 ∂ψ
Br = − ; Bz = useful relation is the r-component of (15.11):
r ∂z r ∂r
∂p 1
4
= (jφ Bz − jz Bφ )
Yes, the notation has changed slightly ... its conventional to ∂r c
use (lower case) ψ for flux functions in non-cartesian problems
(which tend to be laboratory plasmas, in cylindrical or toroidal But now, since p = p(ψ), we can write ∂p/∂r =
geometries). (dp/dψ)(∂ψ/∂r). We can then cancel the ∂ψ/∂r
82
terms, insert the full expression for jφ , and get one form assuming an analytic form for p(ψ) and J(ψ).5 Most
of the G-S equation: realistic p(ψ) and J(ψ) choices requires numerical so-
lution of the equation; this is the current approach in
1 ∂J 2
· 2 µ ¶¸
∂p 1 ∂ ψ ∂ 1 ∂ψ the literature. However, some analytic solutions were
+ + +r =0
∂ψ 2πr2 c ∂ψ 4πr2 ∂z 2 ∂r r ∂r found early on, and are still useful for gaining insight.
(15.23) You will see an example in class or in the homework.
This is the simplest interesting version of the G-S equa-
tion, outside of that in Cartesian coordinates, that I C. Helicity and Taylor relaxation
know. (Most are set up for toroidal coordinate systems,
Now, a different track ... return to the force-free fields
to be relevant in the lab, and are truly horrible ..) If
which we discussed in §14.3. There, we presented some
p(ψ) and J(ψ) are known functions of ψ, then this is a
simple solutions; in this section we discuss why they
PDE in the scalar flux function, ψ. Once we solve this
might be interesting.
for ψ, everything else we need is known (check 15.19,
15.21, 15.22). This subject was motived by laboratory plasmas,
which exhibit self-organization. These plasmas, in tin
3. EXAMPLE : SIMPLE PINCHES cans, tend to evolve spontaneously towards a small
number of preferred states, which are independent of
To solve the G-S equation, one chooses functional the initial conditions of the system. The structure of
forms for dp/dψ, and dJ/dψ, then finds solutions of these preferred states depends only on a few global pa-
(15.23). rameters, such as total magnetic flux, and on the geom-
As a very simple example, let’s assume J(ψ) = 0, etry of the system (say, cylindrical or toroidal). Exam-
and look for a solution independent of z (think of ples are the reversed field pinch (RFP), the spheromak,
plasma confinement in a simple pinch geometry). The and the tokamak (each of which have particular stable
G-S equation becomes field configurations, that I’m not going to describe in
µ ¶ detail here). The process by which the plasma evolves
∂p 1 ∂ 1 ∂ψ from an arbitrary initial condition to its long-lived fi-
+ =0 (15.24)
∂ψ 4πr ∂r r ∂r nal state is called plasma relaxation. Taylor (1974, and
later) suggested that this is a variational process – in
But, we know Bz = (1/r)(dψ/dr); so this equation is which the system evolves toward a state of minimum
just energy, but constrained by one or more invariants. To
describe this we need two mathematical tools.
Bz2
µ ¶
d
p+ =0 (15.25) 1. HELICITY
dr 8π
First, we must define helicity. For a magnetic system,
(does that look familiar?). A slightly more complicated
the helicity is defined by the integral,
solution is the screw pinch. We again assume the equi-
librium is independent of z, but retain J(ψ) 6= 0. If we Z
multiply the G-S equation by (4πr2 )(dψ/dr), we get H = B · AdV (15.28)
let A → A′ = A + ∇χ, where χ is some arbitrary To work with ∂A/∂t, we need to make a choice of
function. The effect of this on the helicity is gauge. Recall the general expression,
Z Z 1 ∂A
E+ = −∇Φ
H′ − H = ∇χ · BdV = χB · n̂dS (15.29) c ∂t
v S
if Φ is the scalar electric potential. Note, further, that
Thus, H is uniquely defined only if B · n̂ = 0 (that is, E + v × B = 0 in an ideal system. If we work now in
if there is no component of B across the boundary; n̂ is a gauge where ∇Φ = 0, we can show H is constant.
the surface normal vector.) • method 1. One way to proceed (following Priest)
is in terms of conditions on A. In an ideal system, we
• How can we understand helicity? Moffatt & have
Tsinober (1992) give two useful illustrations. For
one illustration, consider a circular magnetic flux ∂B ∂A
= ∇ × (v × B) ; = v × B . (15.31)
tube, in which each field line is a circle running ∂t ∂t
parallel to the axis of the tube. The helicity of
this state is (clearly?) zero. Let the magnetic flux The first relation in (15.31) comes from the basic induc-
of this tube be ΦB . Imagine now that we cut the tion equation, for an ideal (η → 0) system. The second
tube at any section, twist it through an angle 2πγ, comes from “uncurling” B in the first, and choosing a
and reconnect. Each field line in the tube is now gauge with Φ = 0. From this, (15.30) can be written,
a torus knot; and any pair of field lines is now Z µ ¶
∂H ∂A
Z
linked. This new system has finite helicity. =− A· · n̂dS + 2 B · (v × B)dV
• Or, consider a second, more complicated exam- ∂t S ∂t V
ple. Two magnetic flux tubes, topologically linked (15.32)
through each other. The integral in (15.28) simpli- But the second term is clearly zero. Thus we find that
fies to two line integrals, around the axes of each H is invariant in systems for which A is held constant
tube: on the boundary.
I I • method 2. An alternative proof comes from
H = κ1 A · dl + κ2 A · dl Biskamp, who works with conditions on B and v. Still
C1 C2 using the fact that B · E = 0, we rewrite (15.30) as
important conservation laws. First, helicity is an invari- From this last, we find that H is invariant if B · n̂ =
ant in an ideal, closed MHD system. Second, the min- v · n̂ = 0; that is, if neither B nor v connect across the
imum magnetic energy available to a system at a fixed boundary.
(invariant) value of the helicity corresponds to that of a
force-free field. Here are the proofs for each of these. Both of these methods show that H is constant for
an ideal, closed system (highly conductive and not con-
2. I NVARIANCE OF THE HELICITY nected to the outside universe).
is minimized subject to the constraint H = constant, in which the Bessel-function solutions (14.20) are seen.
the resultant field is force free. It has also been suggested to occur in solar flares (which
We prove this using a Lagrange multiplier tech- involved a sudden, dramatic release of energy), and
nique. That is, we consider small perturbations to a even in radio jets (where it has been used to predict the
magnetic system, and search for the state which satis- magnetic field structure). We should note, however, that
fies applications to astrophysics, while very tempting, are
hampered by the lack of rigid boundaries. The deriva-
δ (W − αH) = 0 (15.35) tions above made heavy use of the tin-can boundary that
we have in the lab.
In particular, following Priest, consider small perturba-
tions which ahve δA = 0 on S (the surface; again a How does this work? How can magnetic energy be
fixed A thereon), and δB = ∇ × δA. The change in dissipated while magnetic helicity is not? There is still
W is then (linearizing and subtracting αH), a fair amount of discussion about this; the sense of the
Z literature is that this works because magnetic helicity
decays much more slowly than magnetic energy does.
8πδW = 2B · δB − α (δA · B + A · δB) dV Following Bellan, think about the dimensions of mag-
(15.36) netic energy vs. magnetic helicity: W ∼ B 2 L3 , while7
Algebra allows us to rewrite this as H ∼ ABL3 ∼ B 3 L4 , where L is the characteristic lin-
Z ear dimension of the system (or the scale of the Fourier
8πδW = ∇ · (−2B × δA + αA × δA) dV component which holds most of the power). Now, let
Z B be decreased (for instance by resistive dissipation);
+ 2 (∇ × B − αB) · δAdV (15.37) the total energy decays ∝ L3 , while the helicity decays
∝ L4 . Thus, the smaller the scale, the bigger is the ratio
of energy loss to helicity loss.
But the first term again can be written as a surface inte-
We can do this a bit more formally,following Or-
gral, which vanishes due to δA = 0 on S. We are left,
lanti & Schnack. Magnetic energy decays due to
then, with the second term:
resistivity, as
1
Z
δW = (∇ × B − αB) · δAdV (15.38) dW
Z
4π = −η j2 dv (15.40)
dt
But this says that we have a minimum energy state if
the field is force-free: while helicity decays as
dH
Z
δW = 0 if ∇ × B = αB (15.39)
= −2η j · BdV (15.41)
dt
Thus: we have shown that the state of minimum
energy, energy, given a specified value of helicty, is a (to get this last, we started with the definition,
force-free state. (15.28), used dA/dt = −cE, and Ohm’s law).
But now, if we
P think in terms of Fourier compo-
4. TAYLOR RELAXATION nents, B = k bk eikr for a set of wavenumbers
k, then we have
This result has an interesting application in the lab (and
presumably elsewhere). Lab plasmas are observed to dW X
→ −η k 2 b2k (15.42)
“relax” spontaneously. If a plasma starts in some ini- dt
k
tial state (specified by the experimenters, say), it goes
through a turbulent state (which develops due to insta- and
bilities in the initial configuration), and ends up in a
final, longer-lived configuration. Taylor hypothesized dH X
→ −2η kb2k (15.43)
that this corresponds to the plasma dissipating mag- dt
k
netic energy (through self-generated turbulence), while
maintaining its initial value of helicity (as determined Thus, again, at a given k, the energy decay ∝ k 2 ,
by the initial configuration). Its final state should then while the helicity decay ∝ k.
be a force-free state, as in (14.16). This is indeed ob-
served, most famously in an experiment called the Re- 7
verse Field Pinch, which has cylindrical geometry and Remember how A connects to B.
85
References
I’ve taken the flux-function material and solutions
from a mixture of Priest, Bateman and the plasma lit-
erature. The helicity and relaxation discussion comes
partly from me, and partly from references such as Mof-
fatt and Biskamp, as well as the more specialized books
by Bellan, and Orlanti & Schnack.
86
16. MHD EFFECTS IN FLUID FLOWS Integrate this over a finite box (without losses through
the sides):
We now have the basics: we know, in principle, how
d 1 2 1
Z Z
B fields interact with fluid flows. They exert magnetic ρv dV = − j 2 dV (16.3)
tension and pressure on the fluid; the fluid, in return, dt 2 σ
modifies and/or determines the local B field, through Thus, we find that the energy lost from the flow, due
inductive effects (that is, the v × B EMF creates local to the Lorentz force, is just balanced by ohmic losses –
currents). In this chapter we visit some applications. which of course appear as heat in the flow: (Compare
the more general energy equation, at the end of Chapter
A. Magnetic damping and stirring 13 – this is a simple application of the same physics.)
Magnetic fields can have some unexpected effects on 2. M AGNETIC STIRRING .
fluid flow – if you set things up right, they can either
accelerate or decelerate the flow. Here are a couple of Magnetic fields can also induce motion in a fluid. One
brief examples; we may discuss them more in class. simple example is a magnetic field rotated (at some an-
gular frequency Ω) around a fluid that is initially sta-
1. M AGNETIC DAMPING . tionary (think about external electromagnets rotating
around the cylinder).
Magnetic fields can decelerate a flow. For a concrete
example, let’s say you try to drive a flow (maybe a jet)
across a pre-existing B field. Qualitataively, we know Ω
that the flow will initially try to “stretch” the field lines;
the resultant magnetic tension will resist – and decel-
erate – the flow. More quantitatively, the EMF caused
by the flow, v × B, will generate a transverse current j; B
this will give a “backwards” j × B force which will try
to decelerate the flow.
R
This integrates to give the total current, I = jy dz:
v
µ ¶ H o
IBo ∂p Ha Haz
= 2l − 2ρνA sinh (16.14)
c ∂x l l z
B
H
(Note, this latter is an implicit solution for the velocity y
But now, we ignore gradients in the magnetic pres- The above is a perfectly good solution...but Cravens
sure (they are small), and estimate the second term as notes that simple dimensional analysis may not be good
(B · ∇)B ≃ (By Bo /H)ŷ. We thus have two force enough when By is small (and presumably numerical
89
methods may have trouble). We can also gain insight j⊥o H = j⊥1 H, or j⊥0 = j⊥1 (because we’ve cho-
by repeating the analysis in terms of the currents and sen the two slabs to have the same thickness). We can
electric fields induced by the motion. We expect two therefore combine the results above, as
cross-field currents, Pederson and Hall... The Pederson µ ¶
current is the interesting one.1 The Lorentz force in the 1 Dvo dp
v1 = − ρo + (16.21)
upper slab acts to decelerate it; that in the lower slab ρ1 νcoll Dt dy
accelerates the slab:
which allows, again, the possibility of the two drivers,
Dvo /Dt and dp/dy.
top : j × B ≃ −j⊥o Bo ŷ
To finish this, we still need to find j⊥o . Consider
bottom : j × B ≃ j⊥1 Bo ŷ (16.19) the top slab: if it is collisionless, we expect E =
−v × Bo = −vo Bo x̂ in the top slab. But also, be-
cause this is a steady-state, one-dimensional problem,
The two momentum equations become we have ∇ × E = 0, thus E must be the same in the
bottom slab (as well as in the intervening space). It fol-
Dvo dp lows that the Bo lines “map” the field from one slab
ρo + = −j⊥o Bo
Dt dy to the other: the magnetic field lines are equipotentials.
(16.20)
Dv1 We can, therefore, use Ohm’s law (for the Pederson cur-
ρ1 = j⊥1 Bo − ρ1 νcoll v1
Dt rent) to find j⊥o = j⊥1 = σ⊥ vo Bo ..thus finishing the
problem.
We can again find an equilibrium solution for the lower
slab, which has v1 = j⊥1 Bo /ρ1 νcoll . 3. ENERGETICS AND EQUIVALENT CIRCUIT
1
the Hall current is there, in principle infinite in extent, so without
gradients, and not interesting for this problem.
90
17. WAVES AND SHOCKS IN MHD where c2s = γkB To /m is the square of the (unper-
turbed) sound speed.
Nonmagnetized flows have one important characteris-
These equations (17.3) can be combined, algebra-
tic signal speed – that’s the sound speed, as we saw in
ed, and eventually made into one equation for the dis-
Chapter 7. When we add a B field, things get more
turbance velocity, v1 :
interesting ... there is more than one possible “signal-
carrying” wave. ∂ 2 v1
= c2s ∇(∇ · v1 )
A. MHD Waves ∂t2 (17.4)
Bo
In Chapter 7 we found that a density perturbation prop- = [∇ × (∇ × (v1 × Bo ))] ×
4πρo
agates at a speed cs = (∂p/∂ρ)1/2 . This is a criti-
cal quantity in fluid dynamics – the speed at which in- If we now assume a plane-wave solution,
formation can propagate in a gas. For a magnetized
plasma, the situation is (not surprisingly) more com- v1 (r, t) = v1 ei(k·r−ωt) (17.5)
plex. There are two characteristic waves, Alfven and
our linearized equation reduces to
magnetosonic. To start, we carry out a linear analysis
to discover what types of perturbations propagate, and
at what speeds. ω 2 v1 =c2s k(k · v1 )
Bo (17.6)
1. BASIC STRUCTURE : LINEAR ANALYSIS + [k × (k × (v1 × Bo ))] ×
4πρo
Start with a set of basic, unperturbed equations:
Now, finally, we’re in position to analyze basic wave
Dρ modes for this system.
+ ρ∇ · v = 0
Dt As a prelude, note that if Bo = 0, (17.6 reduces to
Dv 1
ρ = −∇p + (∇ × B) × B
Dt µ4π ¶ ω 2 v1 = c2s k(k · v1 ) (17.7)
D p
=0 Dotting this with k, and assuming k · v1 6= 0 (i.e.the
Dt ργ
∂B perturbed velocity has some component parallel to the
= ∇ × (v × B) wavevector), this gives
∂t
∇·B=0 (17.1) ω = cs k (17.8)
Linearize these in the usual way. That is, take for the dispersion relation. This should make you feel
good – we have recovered sound waves.
ρ = ρ o + ρ1 ; v = v1
2. ALFVEN WAVES
p = po + p1 ; B = Bo + B1 (17.2)
These are waves in which the magnetic field dominates.
where the unperturbed (“zero” subscript) state is taken It exerts the restoring force; fluctuations in the plasma
as uniform and homogeneous. (The more general case density and pressure are either exactly zero, or unim-
of background structure, won’t be addressed here). portant. We can start by guesstimating the likely wave
With this, we get speed. Recall that waves in an elastic wire propagate
∂ρ1 due to the tension; as the field lines in a plasma exert a
+ (v1 · ∇)ρo + ρo (∇ · v1 ) = 0 tension Bo2 /4π, one might expect a wave speed
∂t
∂v1 1 Bo
ρo = −∇p1 + (∇ × B1 ) × Bo vA = (17.9)
∂t 4π µ ¶ (4πρo )1/2
∂p1 2 ∂ρ1
+ (v1 · ∇)po − cs + (v1 · ∇)ρo = 0
∂t ∂t This is the Alfven speed, and it is, indeed, a useful scal-
∂B1 ing speed for waves in a magnetized plasma. We can
= ∇ × (v1 × Bo ) also note directly that ∇ · B = 0 ⇒ k · B1 = 0; so
∂t
(17.3) that the magnetic field perturbation must be normal to
∇ · B1 = 0 the wavevector.
91
Following Priest, we now specify to the magnetized- angle terms explicit letting θ be the angle between k
only limit, that is, dropping the po terms from (17.3). and Bo ), and letting b̂o be the unit vector along Bo , we
This is called the cold-plasma limit, and means we’re have the Alfvenic wave dispersion relation:
ignoring any internal-pressure effects (wo we no longer
have sound waves). Doing this, rewriting to make the
ω2 2 2
h i
2 v 1 = k cos θv 1 − (k · v 1 )k cos θ b̂o + (k · v 1 ) − k cos θ(b̂o · v 1 ) k (17.10)
vA
Dotting this with bˆo gives bˆo · v1 = 0; thus, waves in an Alfvenic disturbance described by
this limit just have the perturbed velocity normal to the
ambient magnetic field. Dotting this with k gives B1
v1 = − ; |Bo + B1 | = constant
¡ 2 (4πρo )1/2
ω − k 2 vA2
¢
k · v1 = 0 (17.11) (17.14)
satisfies the full equations, 17.1), without the need to
This has two separate solutions, as follows.
linearize. One important consequence of these waves,
B0 Priest notes, is that finite-amplitude waves do not tend
mag tension
v B1 to steepen, and so dissipate must less readily than other
1
wave modes.
B
mag tension
• compressional alfven waves The second solution
mag tension
to (17.11) is
Figure 17.1. Schematic of (shear) Alfven waves; the
perturbed B1 and v1 terms are perpendicular to the ω = kvA (17.15)
background field Bo . Following Cravens Figure 4.16.
which describes compressional Alfven waves. For these
• (shear) alfven waves. This is the case of an in- waves, the linearized equations (17.3) show that v1 is
compressible perturbation: normal to Bo , and lies in the (k, Bo ) plane. It therefore
∇ · v1 = 0 ; k · v1 = 0 (17.12) has components both along and transverse to k, and so
gives rise to density and pressure fluctuations.
which is, of course, one solution to (17.11). These are,
thus, transverse waves (particle motion transverse to k, 3. MAGNETOSONIC WAVES
B1 ⊥ Bo ). Using this, (17.10) gives When the gas pressure is dynamically comparable
ω = kvA cos θ (17.13) to the magnetic field, the wave nature is different.
We might expect the wave speed to be a mixture of
for Alfven waves (sometimes called shear Alfven compressive effects (through cs ) and magnetic effects
waves). These waves have a phase speed vA cos θ, (through vA ).
which agrees with our simple argument for waves prop-
B
agating exactly along Bo . Putting this solution back 1
k
into the linearized equations (17.3), the first two show v1
ω2 2 2
h i
ˆo + (1 + c2 /v 2 )(k · v1 ) − k cos θ(bˆo · v1 ) k
2 v 1 = k cos θv 1 − (k · v 1 )k cos θ b s A (17.16)
vA
Dotting this with k and bo , in turn, gives two useful B. MHD Shocks; Jump Conditions
relations:
In Chapter 9, we found that density perturbations –
cos θ(bˆo · v1 )
¡ 2
−ω + k 2 c2s + k 2 vA
2
(k · v1 ) = k 3 vA
2
¢
sound waves – will steepen and can develop into
(17.17) shocks. Just as there are more than one signal-carrying
and type of MHD wave, there are more than one type of
MHD shock.
k cos θc2s (k · v1 ) = ω 2 (bˆo · v1 ) (17.18) Also in chapter 9, we applied conservation laws
to determine the jump conditions at hydrodynamic
If k · v1 = 0, we recover the Alfven wave solution shocks. We will follow the same path here, and we will
(17.13). If this isn’t zero, we can combine these two re- find that there is more than one soltuion to the jump
lations to find the dispersion relation for magnetosonic conditions – more than one type of shock. In particular,
waves: magnetosonic waves are compressive, and can steepen
into shocks; the two MS modes, fast and slow, connect
ω 4 − ω 2 k 2 c2s + vA
2
+ c2s vA
2 4
k cos2 θ = 0 (17.19)
¡ ¢
to two types of MHD shocks. Alfven waves, not being
For forward-pointing (ω/k > 0) waves, there are two compressive, do not steepen into shocks (although there
distinct solutions: are formal discontinuities that one can find, associated
with Alfven waves.)
ω2 1¡ 2 2
¢ 1¡ 4 4
¢1/2
2
= cs + vA ± cs + vA − 2c2s vA
2
cos 2θ We work in a frame in which the shock is at rest; this
k 2 2 makes everything steady state. We also ignore dissipa-
(17.20)
tion; it matters within the shock, but here we idealize to
The minus sign in (17.20) gives what is known as an infintely thin jump. We use to the notation of Chap-
the slow mode; the plus sign gives the fast mode. The ter 9, that is [[A]] is the jump in A across a boundary.
Alfven speed lies inbetween these two wave speeds; The unit vectors are n̂, normal to the boundary, and t̂,
thus the Alfven wave is occasionally called the inter- tangential to it. We start with the basic equations, and
mediate mode. for each write down the jump across the shock.
These two magnetosonic modes may be thought of • Maxwell. We know ∇ · B = 0; this
(as Priest notes) as a sound wave, modified by the mag-
netic field, and a compressional Alfven wave, modified [[B · n̂]] = 0 ; Bn = constant (17.21)
by the plasma pressure. If the B field becomes small
(vA → 0), the slow mode disappears, and the fast wave
is the continuity of the normal component, Bn , of the
becomes a sound wave. If the gas pressure becomes
B field. (You recall that the jump in Bt can be finite if
small (cs → 0), the slow mode disappears again, and
there is a surface current).
the fast wave becomes a compressional Alfven wave.
• Mass flux. We have, in a steady state, ∇ · (ρv) =
4. VALIDITY OF MHD WAVE THEORY 0, so that
An important limitation on this MHD wave theory is
[[ρv · n̂]] = 0 ; ρvn = constant (17.22)
that it assumes simple, collective particle motion. At
higher frequencies, the single-particle motions in mixed
is the continuity of mass flux.
magnetic and electric fields cannot keep up with the
driving wave. The wave modes become more com- • Induction. In steady state, ∇ × (v × B) = 0.
plicated as a result, and require a plasma-physical ap- Across the boundary, this gives
proach (following individual particles, or at least indi-
vidual species). The common usage is that we can treat [[(v × B)t ]] = 0 ; [[vn Bt ]] = Bn [[vt ]] (17.23)
waves by MHD methods for wave frequencies below
Ωi = eBo /mi c, the ion gyrofrequency. where we’ve used Bn =constant in the last.
93
let the upstream B field be scaled by b = B 2 /4πρc2s = Solving (17.33) shows us that a perpendicular, mag-
vA2 /c2 . The various equations can be written as, netized shock is compressive only if the upstream flow
s
(relative to the shock) is above the fast magnetosonic
ρ2 speed. That is, X ≥ 1, if M2 ≥ 1 + vA 2 /c2 ; the up-
=X s
ρ1 stream flow must satisfy v12 ≥ c2s1 + vA1
2 . Compared to
3
ρvx = constant (17.34)
B2 1 2 B2
µ ¶ µ ¶
B·v
vx p + − Bx + vx ρe + ρv + = constant (17.39)
8π 4π 2 8π
95
These equations can, in principle, be solved for the [[vx ]] − m [[ν]] = 0 (17.42)
general shock jump conditions. But they are long and
complex, so most authors simplify by working in a use-
ful reference frame.
1
2. FIND A USEFUL REFERENCE FRAME n [[vx ]] + [[p]] + hBy i [[By ]] = 0 (17.43)
4π
At this point, it is very helpful to transform to a frame
moving along the shock face, at a speed
B1y 1
v1y = v1x (17.40) n [[vy ]] − Bx [[By ]] = 0 (17.44)
B1x 4π
µ ¶ µ ¶
2 2 5
2 sin θ1 − ǫ χ − 2 ǫ sin θ1 − 1 + S1 χ − (ǫ + 2S1 sin θ1 ) − 0 (17.48)
3 6
µ ¶
• method # 1. which can be solved directly for χ, given p2 5ǫ 1
=1+ χ− ǫ (17.50)
choices of S1 and θ1 (that is, choices for the upstream p1 3S1 2
magnetic field ratio and the flow angle). Slow shocks
correspond to χ < 0, and fast shocks to χ > 0. The and
more usual ratios can be found from the χ solution:
2
vx1 sin θ1 ǫ
2 = 1+ 2
(χ + sin θ1 ) + (χ + sin θ1 )
bx1 cos θ1 cos2 θ1
(17.51)
ρ2 χ + sin θ1 Woods presents a long discussion of the details of this
=1+ǫ (17.49)
ρ1 1 + χ sin θ1 system and some solutions.
96
Alternatively, the conservation laws, (17.41)-(17.46), Whichever solution method one uses, three solutions of
can be written as these equations can be found.
One has X = 1, is called an intermediate or Alfven
ρvx = constant
wave, and isn’t of great interst as it isn’t compressive.1
Bx = constant In the intermediate wave, the tangential field compo-
à ! nent is simply reversed in sign, so the overall field ro-
By2 Bx (17.52) tates without changing magnitude; hence the name.
ρvv = p+ n̂ − B = constant
8π 4π The other two solutions are compressive, with X >
1, and are slow and fast shocks. The slow shock has
1 2
(v + vy2 ) + h = constant v12 ≤ vA1
2 , and also has B < B . That is, the magnetic
2 1
2 x field bends towards the shock normal in a slow shock.
(recalling that h = γp/(γ − 1)ρ is the enthalpy). Similarly, in the slow shock the parallel flow is slowed
To solve these, we can (following Priest) again use down, v2y < v1y . Conversely, the fast shock has v12 ≥
XVA1 2 , and so B > B . Thus, the field bends away
the compression ratio, X = ρ2 /ρ1 . Keeping our defi- 2 1
nitions, c2s = γp/ρ, and vA2 = B 2 /4πρ (all evaluated from the shock normal in a fast shock. The parallel flow
upstream), the jump conditions become is speeded up in the fast shock, v2y > v1y .
v2x 1
=
v1x X
v2y v 2 − vA1
2
= 21 2
v1y v1 − XvA1
ρ2
=X (17.53)
ρ1
(a) (b) (c) (d) (e)
B2y v 2 − vA1
2
= 21 2 X
Figure 17.5. Left, the changes in magnetic field
B1y v1 − XvA1 direction that are caused by the three types of oblique MHD
shock waves; (a), slow mode shock, in which B bends
(γ − 1)Xv12 X2 − 1
µ ¶
p2
=X+ toward the shock normal; (b) intermediate wave, in which B
p1 2c2s1 X2 “flips”; and (c) fast shock, in which B bends away from the
shock normal. From Priest figure 5.6. Right, the two special
where θ is the angle between the upstream magnetic cases of switch-off (d) and switch-on (e) shocks. From
field and the shock normal. The compression ratio X is Priest figure 5.7.)
now the solution of
¡ 2 2
¢ There are two particular special cases of these. The
v1 − XVA1 F1 (X, θ)
(17.54) slow solution becomes a switch-off shock when v1 =
1 2 2 2 vA1 . In this limit, the tangential field vanishes behind
+ vA1 v1 sin θXF2 (X) = 0
2 the shock, B2y = 0. The fast solution becomes a
where switch-on shock when v1 = X 1/2 vA1 . In this limit, the
upstream field is parallel to the shock normal: B1y = 0.
F1 (X, θ) = Xc2s1
(17.55)
1
+ v12 cos θ [X(γ − 1) − (γ + 1)]
2
and References
The two best references I know for MHD waves and
F2 (X) = (γ + X(2 − γ)) v12
(17.56) shocks are Priest and Woods. I’ve followed both of
− 2
XvA1 ((γ + 1) − X(γ − 1)) them here.
18. MAGNETIC RECONNECTION has evolved around it. In this chapter I’ll try to give an
overview of this culture, both the simple treatments of-
Magnetic reconnection is an innocent-sounding name ten found in the literature, and an impression of what
for an important, and initially unintuitive, process. It might realy happen.
can be defined as the breaking, and topological rear- Space examples are often quoted – the interaction
rangement, of B field lines. That means, it violates of the solar wind with the earth’s magnetosphere, or
what we’ve learned about magnetic flux conservation – the rapid energy release that creates a solar flare. But
which says that the B lines are “frozen into” the plasma. reconnection is also important in laboratory plasmas,
Magnetic flux is conserved only when resistive dissipa- where it underlies the “Taylor relaxation” process we
tion can be ignored; reconnection is an effect of local- talked about in Chapter 15; and it is critical to the exis-
ized, intense dissipation regions. It results in magnetic tence of the dynamos that are responsible for the earth’s
energy being converted to plasma energy (both heat and and sun’s magnetic fields.
bulk motions).
Because reconnection allows field lines to be rear- A. Advection vs. dissipation of magnetic field
ranged, and to “move through” a plasma, it can “con- To set the stage, think back to the induction equation,
nect” regions of plasma which had previously been
magnetically isolated. To be specific: think of space ∂B
plasma, where the electrical conductivity is very high = ∇ × (v × B) + η∇2 B (18.1)
∂t
(all those free charges can move where they want to),
and therefore the flux freezing limit is a good one. This The first term on the right-hand side describes inductive
means that the plasma particles can freely move and effects of the plasma flow on the B field; the second
mix along field lines, but cannot do so across the field. term describes resistive diffusion or dissipation. The
A particle always remains tied to the same field line as ratio of the two terms is given by the magnetic Reynolds
it moves with the plasma flow. Think, then, about two number,
initially separate plasma regions which come into con-
tact with one another.1 It follows that the two plasmas Lv
Rm = (18.2)
will not mix; rather, a thin boundary layer will form η
between them, which keeps them quite separate. (The
location of the boundary layer is, of course, determined where the magnetic diffusivity η = c2 /4πσ and σ is the
by pressure balance and large-scale fluid dynamics.) In electrical conductivity. It’s worth noting, here, that we
general, the magnetic fields on either side of the bound- can “uncurl” (18.1), to get
ary layer will have different strengths and orientations,
v η
so that a current sheet must exist within the boundary E= ×B+ ∇×B (18.3)
layer. c c
Thus, flux freezing leads inevitably to the predic- – and remember from Maxwell that ∂B/∂t = −c∇ ×
tion that in plasma systems space will be divided into E; so there is also an E field in the system when we
separate cells, which contain magnetized plasma from have either advection (v × B) or gradients (∇ × B).
different sources, separated from each other by thin cur- Now, envision a region of adjoining, magnetically
rent sheets. Based on what we’ve seen up to now, the isolated regions, as in Figure 18.1, and think about how
plasmas from these separate domains should be able to it might evolve.
mix only very slowly, as resistivity dissipates the cur-
rent and magnetic flux. Looking ahed to (18.5), this
should happen on the (usually very slow) timescale,
B
τdiss ∼ L2 /η. But: the answer is much more inter-
esting. A mixture of dynamics and resistivity allows x j
the topology of the magnetic field to change, and the
plasmas to mix, on a much shorter timescale. This pro- B
cess is called magnetic reconnection, and quite a culture
Figure 18.1. The setting for a reconnection event – the
1
To be specific, think of the solar wind where it encounters the oppositely directed B field separates plasma in the regions
earth’s magnetosphere. Both of these are diffuse, magnetized, above and below the current sheet (j, located at x ≃ 0).
flux-frozen plasmas.
98
1. RESISTIVE LIMIT: FIELD DECAY ON AXIS a parameter of the problem .. and so people usually
prefer to estimate an Alfven time:
First, consider the highly resistive case, when Rm ≪ 1.
We have simple diffusion, L
τA ≃ (18.8)
∂B vA
= η∇2 B (18.4)
∂t where vA is evaluated at large distances from the cur-
As we saw earlier, we expect the field close to the x = 0 rent sheet.
line to decay, and the current sheet to broaden, with We might expect that τA would be the best we could
time (the first is also described as B lines “diffusing do to annihilate magnetic flux – simply drive the field
through the plasma”), This basic equation has well- into the dissipation region. We’ll see below that things
known solutions; I’m storing one in the Appendix to aren’t this easy – that τA underestimates the timescales
this chapter. But we don’t need the full solution to (or overestimates the reconnection rate). The inflow
guesstimate the timescale for this process. If the x- must be regulated by the dissipation rate ... so the final
extent of the system is initially ∼ L, then from (18.4) answer for the reconnection timesscale will be some-
we can directly estimate the diffusion time: where between τA and τdif f .
L2 B. Steady, 2D Reconnection
τdif f ≃ (18.5)
η
Our next task, then, is to combine advection and diffu-
This limit destroys magnetic energy, and heats the sion. When we keep both terms in (18.1), we lose the
plasma – but it doesn’t, by itself, change the magnetic chance at simple analytic solutions. Instead, we’ll work
topology. with dimensional and scaling arguments2
Alternatively, consider a highly conductive case, when To understand how resistivity can “break” and “recon-
Rm ≫ 1. Now, the first term in (18.1) dominates; nect” field lines, think about the geometry in Figure
18.2. We know resistivity is important in regions of
∂B high current density - such as the central region (around
= ∇ × (v × B) (18.6)
∂t OP). If we set up this geometry and waited awhile, the
and we talk about magnetic flux freezing. In the context central magnetic field would decay as the current layer
of Figure 18.1, envision the plasma flowing towards the supporting them is dissipated. This would deplete the
current sheet, with some velocity Vo . Because of flux magnetic pressure in this region. The plasma above
freezing, the B lines will be carried with the plasma, and below this region would be pushed inwards by its
in towards x = 0. Thus, B will grow, close to the x- own pressure, bringing a in a fresh supply of field (and
axis (the field lines must bunch up as they’re carried in- plasma).
wards). I’ve also put one analytic solution to this prob-
A V D
lem in the Appendix to this chapter. Once again, we i
E F
can also estimate a (dynamic) timescale for this pro-
cess, without knowing the details:
B C
L
τdyn ≃ (18.7) V
o O P V
o
Vo
B’ C’
Its also worth noting that this process can’t continue
forever – at some point the magnetic energy will build E’
up enough to exert a back pressure and shut down the Vi F’
flow ... unless this energy is dissipated as heat which A’ D’
can escape the system. That means that we have to in- Figure 18.2. Illustrating a sinple reconnection geometry;
see text for discussion. Following Choudhuri figure 15.2.
corporate both advection and dissipation ... which leads
us to reconnectoin.
But first, one more timescale. We’ll find below that 2
This is the true state of the field, folks .. it’s either these scaling
the inflow velocity Vo is something to be solved for, not arguments, or numerical simulations.
99
Now, look at this in more detail (I’m following annihilating magnetic field and energy (where does it
Choudhuri’s discussion here). The field lines ABCD go?).
and A’B’C’D’ move inwards, with velocity vin . Even- • Next, use the induction equation. In a steady state,
tually the BC and B’C’ parts of the field lines decay it is ∇ × (v × B) = η∇2 B, which gives (by dimen-
away. The AB part of the field line is moved to EO, sional/scaling analysis)n
and the A’B’ part of that field line is moved to E’O.
Thus, these “fragments” of two original field lines now vin B ηB η
make up one new field line, EOE’. And similarly, the ≃ 2 ; vin ≃ (18.11)
l l l
parts CD and C’D’ eventually make up a new field line,
FPF’. Thus, “cutting and pasting” of field lines (other- This is saus that we can keep the diffusion region steady
wise known as reconnection) takes place in the central if the rate at which flux is brought in is equal to the rate
region. And, of course, there must be plasma flow away at which it is annihilated. We can also note that there
from the region – sideways in this cartoon – to conserve must be an E field, as shown in the figure, to maintain
mass. the current sheet. From Maxwell, we get
Bo2
v ≃ pmax − po (18.13)
out 8π
l
v Along and within the current sheet (call that the ŷ di-
L out
rection), there is no v × B force, so only the pressure
gradient accelerates the flow. We have then,
v
in 2
∂vy ∂p vout
Figure 18.3. Geometry of Sweet-Parker reconnection. ρvy ≃− ; ≃ pmax − po (18.14)
The current sheet (grey shaded area) has transverse width l ∂y ∂y 2
and lateral extent L The input velocity is vin and the output
velocity is vout . The induction E ∝ v × B is normal to the Thus, the outflow speed must be
system (into/out of the paper), as is the current in the current
sheet. Quantities far away from the current sheet are 2 2 Bo2
labelled with subscript o. Following Cowley figure 5.5. vout ≃ vA = (18.15)
4πρ
Our analysis is strictly two-dimensional and steady OK: we can combine these three parts to get our
state. I’ll assume the flow is incompressible – that it main results, namely, the inflow velocity ahd thickness
stays at constant density (which turns out to be a good of the dissipation layer:
approximation if vin and vout are both subsonic. There
are three important steps here. 2 vA η ηL
vin ≃ ; l2 ≃ (18.16)
• First, mass conservation requires L vA
vin L = vout l (18.9) To connect with the literature, we can also rewrite our
results in terms of Rm:
From this, directly, conservation of magnetic flux gives
2
vAi B2 L2
Bo vin = Bout vout (18.10)
2
vin = 2
; Bout = o ; l2 = (18.17)
Rm Rm Rm
Thus: the outflow speed must be significantly faster Thus, when Rm ≪ 1 – which we expect for the highly
than the inflow speed (because we’re assuming l ≪ L). conductive plasmas in our examples – reconnection is
The outflow field is much smaller than the inflow; we’re indeed slow.
100
solutions assume that the diffusion region will adjust electric fields (again on microscopic scales; they won’t
itself to respond to the external boundary (driving) con- affect the macroscopic dynamics that we’re studying in
ditions, in doing so making L∗ ∼ l and thus making the this course).
reconnection rate nearly independent of η.
Plasma physics details: the most common origin
Do Petschek-type mechanisms work? It looks at of microturbulence is streaming of charges rela-
present as if the Petschek mechanism is not obtained tive to the background plasma. If the streaming
in real situations. Several numerical simulations (e.g., speed is high enough – typically
Biskamp 1986, Ugai 1995) have looked for steady p greater than the
plasma thermal speed, vth ∝ kT /m – the sys-
reconnection solutions, and found that simple, two-
tem is unstable, and some of the streaming kinetic
dimensional systems will reach the Sweet-Parker so-
energy is turned into microturbulence. The most
lution for low inflow rates. At higher rates, however,
common way to set up relative streaming is with
the simulations find a transition to unsteady behavior,
a high current density. As above, j = nevD ; a
rather than a transition to the Petschek mode. In addi-
higher j leads to a higher drift velocity vD , and
tion, recent lab work which was configured explicitly
when vD > vth streaming instabilities occur. But
to measure reconnection physics (cf.review by Kulsrud,
the current density is controlled by external con-
1998) seems to support the Sweet-Parker model.
ditions – for instance by the scale over which B
2. COMPRESSIBLE FLOW changes significantly, in the simple reconnection
models. Thus, thin current sheets ⇒ high current
Another suggestion is that the plasma in the dissipa- density ⇒ microturbulence.
tion region (DR) is compressed. If this is the case, then Once microturtubulence is created, stochastic electric
(18.9) is replaced by vin ρo L ≃ vout ρDR vout . This will fields will scatter the plasma particles; if the plasma tur-
increase vin relative to vout – its a plausible idea, but bulence is strong enough (and very weak levels will do)
I’m not aware of any lab tests yet. This will be investi- this turbulent scattering will be much more important
gated in the homework. than single-particle interactions. Because this is a non-
linear process (at what amplitude does the turbulence
3. ANOMALOUS DIFFUSION
saturate?), we can’t write down collision times from
This is my personal favorite, and is a good example of first principles. One common approach (way to get
how plasma microphysics can be important in macro- around this) argues that we can get an upper limit to the
scopic situations. That is: the electrical resistivity σ, resistivity (lower limit to τcoll is to set τcoll = 2π/ωp ,
which sets the size of the η term, is traditionally cal- if ωp2 = 4πne2 /m is the square of the plasma fre-
culated assuming Coulomb collisions are the mecha- quency (in cgs units). This is used in (18.20) to estimate
nism for energy transfer and dissipation. (Check the anomalous conductivity or anomalous resistivity.
appendix of Chapter 1 for Coulomb scattering, and the
D. Other approaches
back of Chapter 13 for conductivity discussions). The
critical point is that electrical conductivity is regulated Reconnection is a very active field these days. I’m just
by the collision rate of the current-carrying electrons: putting a few short comments here; once again, these
notes would be very long if I tried to review the entire
j nevD ne2 field.
σ= = ≃ τcoll (18.20)
E E m
1. SPONTANEOUS RECONNECTION
where τcoll is the electron collision time.
Now, Coulomb collisions are the “hands-off” The simple models assumed a steady state, with some
particle-particle collisions that result from the long- mass inflow to balance the flux annihilation. Thus they
range nature of the electrical force. They must happen are “driven” in a general sense. But the plasma will find
in any ionized plasma – but they are slow (at least in its own way; oppositely directed magnetic fields across
low-density plasmas, such as space or the lab), and thus a thin current sheet will spontaneously reconnect. This
not very effective at dissipating a current. Are Coulomb is called a tearing instability; and typically results in
collisions the only dissipation mechanism? The an- formation of magnetic islands in the current sheet. Fig-
swer is, almost certainly not. If the plasma is turbu- ure 18.6 illustrates this behavior.
lent on microscopic scales (which is very likely to be This instability can be studied analytically; the math
the case), and if the turbulence involves charge separa- is long and horrendous, and I’ll put it in a later chapter.
tion (ditto), we expect the plasma to contain fluctuating It can also be studied with numerical simulations; an
102
example of that will also be in that later chapter. field, in a reconnection event. Flares have motivated
much of the work on steady-state models; however as
an outsider I suspect time-dependent, patchy, localized
reconnection must be taking place.
• Reconnection at the magnetopause, where the so-
lar wind hits the earth’s magnetic field. Spacecraft ob-
servations suggest this is very patchy, localized and
time-dependent. Picture, for instance, magnetic flux
tubes being carried along in the solar wind; and let one
of them impact the magnetopause, which also has its
field bunched into ropes. This process can allow solar
wind plasma, and field, to penetrate into the magneto-
sphere.
There are many numerical simulations of dynamic
Figure 18.6. Effect of a perpendicular plasma reconnection in the literature; I’ll bring some to class.
displacement in a neutral sheet, leading to field compression
for the ideal case (η = 0; top figure) and to field disruption 4. THREE - DIMENSIONAL RECONNECTION
and island formation for the resistive case (η 6= 0; lower
figure). Following Biskamp Figure 4.4. Another observation: only rarely can a reconnection
event be well described by a two-dimensional analy-
sis. My last example in fact assumed this – because
2. DRIVEN RECONNECTION
the intersection of two magnetic flux ropes is clearly a
This is a different approach. It’s worth remembering three-dimensional process. Going to 3D is challenging,
that the 2D models above are all “passive”: put two mis- and work is only starting here (helped significantly by
aligned B fields together and wait to see how quickly recent advances in resistive MHD codes). If I can I’ll
they reconnect. But nature does not always work this bring some recent images to class.
way. One can envision a situation in which the two
anti-parallel magnetic structures are driven together, by
large-scale flows in the system (one example of this is
MHD turbulence, in which different parts of the plasma References
move in random directions, at a speed set by the ener-
Some of this is “just from me”, from my reading
getics of the turbulence). In this case one (this one at
in the field. I’ve taken the basic Sweet-Parker and
least) expects that the inflow speed (vin ) in our nota-
Petschek material from Priest; and the two formal so-
tion above) will be set by the large-scale flow (say the
lutions in the Appendix from Priest & Forbes.
turbulent speed). How can the reconnection site adjust
to this? Some authors (e.g. Parker) suggest that the
internal structure of the reconnection layer – its thick-
ness, density or resistivity (set by microscale turbulence
therein) – will adjust as necessary. I personally tend to
agree with them. E. Appendix
Here’s another personal impression: who says steady- Refer back to Figure 18.1, so that we have a 1D prob-
state reconnection is relevant to any natural situation? lem, B = (0, B(x, t), 0); and take the figure as an ini-
That is: the arguments above show that steady re- tial condition
connection models must be forced, sometimes rather
severely, to connect with what we think is occuring B = Bo , x > 0 ; B = −Bo , x < 0 (18.21)
in nature. But examples of non-steady reconnection
events are easy to find: The time evolution of the system is governed by
• Reconnection in solar flares. Flares are seriously
transient events; the large amount of energy released is ∂B ∂2B
believed to be due to very fast annihilation of magnetic =η 2 (18.22)
∂t ∂x
103
There are scads of different ways to solve this equation. 2. ADVECTION - ONLY SOLUTION
One uses a Green function:
Z Refer again to Figure 18.1, but this time specify a ve-
B(x, t) = G(x − x′ , t)B(x′ , 0)dx′ (18.23) locity field which describes inflow towards x = 0, but
also allows flow out the sides (to conserve mass):
and for this initial condition, the Green function is x y
vx = −Vo ; vy = Vo (18.28)
a a
1 ′ 2
G(x − x′ , t) = e−(x−x ) /4ηt (18.24) This satisfies ∇ · v = 0 (incompressible); and has
(4πηt)1/2 streamlines xy = constant.4 As an initial condition,
From this we get the desired solution, in terms of an pick a B field which increases going away from x = 0,
error function: and changes sign at x = 0:
µ ¶ x
2Bo x B(x, t = 0) = Bo cos (18.29)
B(x, t) = √ erf √ a
π 4ηt
Z x/√4ηt (18.25) Figure 18.8 shows this geometry. If we ignore diffu-
2Bo −u2 sion, our governing equation is
= √ e du
π 0
∂B x ∂B Vo
− Vo = B (18.30)
∂t a ∂x a
B(x,t)
B y
o
−Bo
This has the expected behavior, as illustrated in Fig- Figure 18.8. Solution for the magnetic field as a function
of time, being advected into a current layer at x = 0. The
ure 18.7: the B field decays with time, near the cur- dotted lines are velocity streamilines; the magnetic field
rent sheet, which itself spreads with time. One can also builds up close to x = 0 as time goes on. Solution details
show that the total current stays constant: are in the text. Following Priest & Forbes figure 3.3.
c ∂B cBo
Z
j= ; J = jdx = (18.26) Priest and Forbes present a cute solution to this.
4π ∂x 2π Consider the characteristic lines (in the (x, t) plane)
and that the magnetic energy decays with time: x(t) = x∗ e−Vo t/a (18.31)
duB B ∂B B ∂2B (clearly x∗ is the initial value of x, on a given velocity
Z Z
= dx = η dx line, at t = 0). Along these lines, B obeys
dt 4π ∂t 4π ∂x2
∂B 2 dB ∂B dx ∂B
Z µ ¶
η
=− dx (18.27) = −
4π ∂x dt ∂t dt ∂x (18.32)
∂B x ∂B Vo
4π 2
Z µ ¶ Z 2
η j = − Vo = B
=− j dx = − dx ∂t a ∂x a
4π c σ
(where we’ve integrated by parts in the second line, and 4
Check back to Chapter 1 for the definition of streamlines, and
used basic Maxwell in the third). how to find them from the velocity field.
104
(the first line is just chain rule, and the second line no-
tices connects back to eqn 18.30). But now it’s easy:
along the characteristic lines, (18.31), the solution is
³x ´
∗
B(x, t) = Bo eVo t/a cos (18.33)
a
All we have left to do is write this in terms of (x, t)
everywhere:
³x ´
B(x, t) = Bo eVo t/a cos eVo t/a (18.34)
a
This shows that the field accumulates – grows in
strength – near the x = 0 axis; in fact it concentrates
there (closer and closer to x = 0 as time goes on.
However, this clearly isn’t a physical solution, be-
cause if B grows exponentially, magnetic pressure will
react back on the flow.
105
of a light fluid into a heavy one (as in an explosion in a (??) represents a perturbation with grows exponentially
dense atmosphere) is also subject to this instability. with time: the system is unstable. Let k 2 = kx2 +ky2 , and
denote z-derivatives by Df = df /dz (In this and the
2. THE MATH next section only; Chandra’s notation). These equations
Our basic approach here is modal analysis. That is: first then become
we linearize, assuming small perturbations, so that we
can drop terms which are second order in these pertur- ikx δp = −iωρu
bations. We have done this before – in deriving sound
waves and MHD waves, for instance. Our goal is to find iky δp = −iωρv
the time evolution of these perturbations: do they grow
Dδp = −iωρw − gδρ (19.3)
with time (which is an instability), or simply oscillate
(which says we have wave solutions, but not instabili-
ikx u + iky v + Dw = 0
ties). In order to do this, we do a modal analysis. We
assume that any perturbation can be treated as a sum of iωδρ = wDρ
Fourier components (which is legal if the defining equa-
tions are linear in the perturbations – right??). We can These can be used to show
then focus on the evolution of one such Fourier compo-
nent – for instance (19.2), below. k 2 δp = −iωρDw (19.4)
My analysis here follows Chandrasekhar, Hydrody-
namic and Hydromagnetic Stability. He works through and
a modal analysis, and we will pick up on some of his g
results. Recall our formal derivation of sound waves, Dδp = −iωρw + wDρ (19.5)
iω
in Chapter 4: we carried through a linear analysis, ex-
panding the variables v, p, ρ in (unperturbed) + (small These two now combine – eliminating δp between them
perturbation) terms. Our basic equations are those – as
of mass conservation, incompressibility of the unper-
gk 2
turbed state, and momentum conservation. Geome- D (ρDw) = wDρ + ρwk 2 (19.6)
try: take g = gẑ; assume the unperturbed state had ω2
v = 0, and only has gradients on the ẑ direction. and let But now, this is an equation for w(ω, k), because Dρ
the perturbed velocity components be δv = (u, v, w). is specified in the problem. Picking w = 0 on the (dis-
The linearized equations – for the perturbations δρ, δp, tant) boundaries of the system, this becomes a Sturmian
keeping only lowest order “small” terms, are: characteristic value problem for ω 2 . One can show,
∂u ∂ then, that ω 2 > 0 if Dρ < 0 everywhere; and con-
ρ = − δp versely, ω 2 < 0 if Dρ > 0 somewhere in the fluid.
∂t ∂x That is: a density structure which decreases with ẑ is
∂v ∂
ρ = − δp stable; while a density structure which increases with ẑ
∂t ∂y is unstable.
∂w ∂ The next step is to find out the growth rates, and
ρ = − δp − gδρ (19.1)
∂t ∂z whether all wavevectors are unstable. To do this, we
∂u ∂v ∂w simplify to the case of two uniform fluids, separated by
+ + =0
∂x ∂y ∂z a horizontal interface at z = 0. Now, away from this
∂ dρ boundary, (19.6) becomes
δρ = −w
∂t dz
D2 w = k2 w (19.7)
Now, assume each perturbed quantity has the form (a
“mode”) which has simple exponential solutions. Picking the
ones which stay well-behaved at z → ±∞, we have
δf (x, y, z, t) = δf (z)ei(kx x+ky y+ωt) (19.2) away from the boundary
We will generally treat k as given, and want to find the
nature of ω. In particular: if ω is real, (19.2) simply w(z) = Ae−kz ; z>0
represents waves or oscillations – so that the system is (19.8)
kz
stable. However if ω has a negative imaginary part, then w(z) = Ae ; z<0
107
To connect solutions through the boundary, we need (Why is it proportional to the perturbed velocity?) The
jump conditions. Define perturbed momentum equations become
∆s f = f (zs + 0) − f (zs − 0) iρωu = −ikx δp
that is, the jump in f from the below some surface zs to B
above that surface. Here zs = 0. Applied to (??), we iρωv = (ikx by − iky bx ) − iky δp
4π
find the jumps in p and ρ are related by B g
iρωw = (ikx bz − iDbx ) − iky δp + wDρ
g 4π iω
∆o δp = ws ∆ o ρ (19.9)
ω (19.14)
Manipulation of the system (??) also gives
Now: using (19.11) to convert from b components to v
2
k ∆o δp = −∆o (−iωρ) (19.10) components, we can show
kx2 B 2
µ ¶
Finally, these combine to give iρω + (ikx v − iky u) = 0 (19.15)
iω 4π
k2
∆o (ρDw) = gw∆o ρ (19.11) which requires kx v = ky u. Using this result, and fur-
ω2
ther algebra, we find
But in this system, ∆o ρ = ρ2 − ρ1 (upper = lower val-
ues); and Dw = ±kw (- for upper, + for lower). The iρωDw = −k 2 δp (19.16)
result (19.11) thus becomes
and
µ ¶
2 ρ2 − ρ 1 B 2 kx2 ¡ 2 g
ω = −gk (19.12) iρωw − D − k 2
¢
w = −Dδp + wDρ
ρ2 + ρ1 2
4πk iω iω
(19.17)
This is our result: ω 2 < 0 – the system is unstable – And now, we can look for the new results. First,
if ρ2 > ρ1 . This recovers what we guess from potential combine (19.16) and (19.17) to get
energy arguments, at the start. We also find that all
wavenumbers are unstable 2
√ (that is, have ω < 0), with B 2 kx2 ¡ 2 gk 2
D (ρDw)−ρwk 2 − 2
¢
D − k w = wDρ
a growth rate |ω| ∝ k. This last result is true for 4πk 2 ω 2 ω2
small perturbations; however when the instability starts (19.18)
to grow, larger-scales dominate the system. which is the generalization of (19.4) for the non-
magnetic case.
This analysis could also have included surface ten-
sion at the interface: when included, one finds that sur- From this, first, we note that if kx = 0, which cor-
face tension stabilitizies high wavenumbers (small spa- responds to propagation perpendicular to the field, (??)
tial scales). just reduces to (19.6): cross-field perturbations have no
effect on the RT instability. The ones with kx 6= 0,
3. THE MAGNETIZED CASE however, do change the results. We can repeat the non-
magnetized analysis . . . the jump condition (19.11)
What happens if we add a magnetic field? We expect generalizes to
little effect from a vertical field, as it doesn’t hamper
vertical fluid motions. This turns out to be the case – B 2 kx2 k2
the stability conditions are not affected. It does have ∆o (ρDw) − ∆ o (Dw) = gw∆o ρ (19.19)
4πω 2 ω2
some effect on the growth rate, however (note that a
Again taking the case of two uniform fluids separated
vertical perturbed velocity must also involve horizon-
by an interface, the result (19.12) becomes
tal, cross-field flows – by the continuity equation). The
more interesting case is a horizontal magnetic field: we 2B 2 kx2
µ ¶
2 ρ2 − ρ 1
expect this to affect the stability as well as the growth ω = −gk + (19.20)
ρ2 + ρ1 4π (ρ2 + ρ1 )
rates, and it does.
In this section I still follow Chandra. Consider a Thus: a horizontal magnetic field stabilizes high fre-
uniform B = B x̂, with perturbed field quency perturbations. Instability (ω 2 < 0) requires
kx 2B 2 kx2
δB = (bx , by , bz ) = Bv (19.13) gk (ρ2 − ρ1 ) > (19.21)
ω 4π
108
Finally, we note that this analysis assumed an in- The basic perturbed equations are
compressible fluid. That simplifies the algebra but isn’t
∂u ∂u dU ∂
required for the physics. Extending to the compressible ρ + ρU + ρw = − δp
case (as done by Shivamoggi, Theory of Hydromagnetic ∂t ∂x dz ∂x
Stability), one finds that compressibilty has a stabiliz- ∂v ∂v ∂
ρ + ρU = − δp
ing effect, in that it reduces the growth rates in both the ∂t ∂x ∂y
magnetized and unmagnetized cases, compared to their ∂w ∂w ∂
ρ + ρU = − δp (19.22)
incompressible analogs. ∂t ∂x ∂z
∂u ∂v ∂w
C. The Kelvin Helmholtz Instability + + =0
∂x ∂y ∂z
Well, that was so much fun, let’s do it again. The other ∂ ∂ dρ
δρ + U δp = −w
well-known instability of an interface is the Kelvin Hel- ∂t ∂x dz
moltz Instability. We again assume the perturbation has the form, (19.2).
The system (??) becomes
1. THE PHYSICS
ikx δp = −iνρu − ρwDU
This instability arises at a velocity shear. Consider two
fluids, or two pieces of the same fluid, in relative mo-
iky δp = −iνρv
tion. For instance, think of a horizontal interface; fluid
below is at rest, and fluid above is moving parallel to Dδp = −iνρw (19.23)
the interface. (This might describe a situation in the at-
mosphere, with high-velocity winds going past a low- ikx u + iky v + Dw = 0
velocity cloud layer). When the interface is displaced
upwards, it forces the overlying flow to deviate over iνδρ = −wDρ
the perturbation. This means the pressure there drops
(Bernoulli, remember?). The boundary thus feels a lift. where we have defined the “shifted” frequency, ν =
But of course the situation is mirror symmetric across ω + kx U . From these we can find
the interface; the downward displaced boundary feels
a downward “lift”. As the boundary follows these per- iρνDw − iρkx wDU = −k 2 δp (19.24)
turbations, it is sheared by the flow into which it pene- and
trates, and rollw up into a vortex sheet.
ik 2 Dδp = ρk 2 νw (19.25)
(compare equations 19.4 and 19.5). These combine
v
2 (eliminate δp) to give
D [ρνDw − ρkx wDU ] − k 2 ρνw = 0 (19.26)
and from this, our useful jump condition , at an inter-
v
1 face zs , is
Figure 19.2. The origin of the Kelvin-Helmholtz ∆s [ρνDw − ρkx wDU ] = 0 (19.27)
instability. Shear at an interface between two fluid layers in
relative motion generates a vortex sheet; pressure (what quantities must be held constant across an inter-
differences (think of the Bernoulli effect) lead to instability, face? how did we throw out some of the terms in ???)
and eventual mixing of the two fluids.
We again apply this to two uniform fluids. Above
and below the boundary, we again have (19.7), as all
2. THE MATH
else is constant. But now, the ratio w/(ω + kx U ) must
be continuous across the boundary (why?). Thus, (19.8)
We again follow Chandrasekhar. The methods are sim- is replaced by
ilar to those presented for the RT instability, but the ge-
ometry is slightly more complicated. We add a hori- w(z) = Aν2 e−kz ; z>0
zontal flow U = ux̂ to the unperturbed system, but we (19.28)
kz
now ignore gravity. w(z) = Aν1 e ; z<0
109
(where the subscripts 2,1 refer to (top),(bottom) flow U; we don’t expect cross-field flows to differ from the
speeds). Our characteristic equation is simpler here. unmagnetized case.
Taking ρ1 = ρ2 = ρ, we can derive The modal equations become
(ω + kx U2 )2 + (ω + kx U1 )2 = 0 (19.29)
But this has only complex roots (if kx 6= 0). The ikx δp = −iνρu − ρwDU
quadratic solves to
B
1 i iky δp = −iνρv + (ikx by − iky bx )
ω = − kx (U1 + U2 ) ± kx (U1 − U2 ) (19.30) 4π
2 2 B
Dδp = −iνρw + (ikx bz − Dbx )
4π (19.31)
Thus, any such perturbation is unstable (ω has an
imaginary part), no matter how small the difference ikx u + iky v + Dw = 0
U1 − U2 may be.
iνδB = ikx Bv + bz DU
As Chandrasekhar quotes Helmholtz: “Every per-
fect geometrically sharp edge by which a fluid iνδρ = −wDρ
flows must tear it asunder and establish a surface
of separation, however slowly the rest of the fluid
may move”.
From these, after a fair bit of algebra, and again using
Equation (19.30) shows that higher k’s grow the the “shifted” frequency ν = ω+kx U , we find (compare
fastest when the perturbation is small. When it becomes 19.16 and 19.17)
finite, however, this analysis breaks down; experiment
shows that a characteristic and large scale will dominate
the instability.
ρνDw − ρkx wDU = ik 2 δp (19.32)
3. THE MAGNETIZED CASE
2 k2 w 2
· µ ¶ ¸ µ ¶
2 2 B Dw 3B DU
ik Dδp = ρk νw + kx2 D − − kx D w (19.33)
4π ν ν 4π ν2
These two combine as
2 k2 w 2
· µ ¶ ¸ µ ¶
B Dw 3B DU
D (ρνw − ρkx wDU ) + kx2 D − − kx D w (19.34)
4π ν ν 4π ν2
4B 2
(U1 − U2 )2 ≤ (19.38)
4πρ
References
As in the text, I’ve followed Chandrasekhar, Hy-
drodynamic and Hydromagnetic Stability for the hydro,
and Shivamoggi, Theory of Hydromagnetic Stability),
for the MHD.
111
W
20. IDEAL MHD INSTABILITIES W
(perturbed position) and the change in the potential energy, due to this per-
11
00
00
11 turbation, is
ξ
r
11
00
00 (equilibrium position)
11 1
Z
δW = − δx · FdV (20.5)
2
r
o
B12
Z · ¸
1 j 2
δW = − δx⊥ · × B1 + γp|∇ · δx| + (δx⊥ · ∇p)∇ · δx⊥ dV (20.6)
2 4π c
An alternative version of this which explicitly uses the field line (inverse) curvature, ~κ = b · ∇b if b is the unit
vector along B, is
Z " 2
1 B1⊥ B 2
δW = + |∇ · δx⊥ + 2δx⊥ · ~κ|2
2 4π 4π
# (20.8)
jk
+ γp|∇ · δx|2 − 2(δx⊥ · ∇p)(~κ · δx⊥ ) − (x⊥ × b) · B1 dV
c
This last form is the so-called “intuitive form”1 ; each some equilibrium configuration (as in Chapter 17, say),
term allows a simple physical interpretation. The B1⊥ and asking whether some perturbation δx can be found
term represents the energy needed to bend field lines. for which δW < 0. If this is the case, then we know
The second term is the energy needed to compress the our initial configuration is unstable. (The converse isn’t
magnetic field. The γp term, the third term, is the en- necessarily useful: if we try some perturbation and find
ergy needed to compress the plasma. Each of these con- δW > 0, this doesn’t prove stability unless all pertur-
tributions is stabilizing. The last two terms can be pos- bations have been considered).
itive or negative, and thus can drive instabilities. The
first of these depends on ∇p ∼ j⊥ × B, while the sec- B. Apply: Pinch Instabilities
ond depends on the parallel current jk . Thus, a vacuum
field is stable (but not very interesting), while either 1. THETA PINCH
perpendicular or parallel currents are potential sources
of instability. We first look at stability of the θ pinch: which is simply
Thus: the energy method is applied starting with described (cf. §17.3) by
1
oh yeah? Bz2 B2
p(r) + = o (20.9)
8π 8π
113
πL a
Z
δW = δW (r)dr (20.13) Again evaluating δW (r) (to go in 20.12), one finds
4π 0
m2 Bφ2 ¡ 2
where the integrand can be written (with ko2 = k 2 + 2 (20.18)
¢
δW (r) = ξr + ξz
m2 /r2 ) r2
· µ ¶ ¸2
¶2 2 d ξr 8π dr 2
+ Bφ r + ikrξz + ξ
µ
δW (r) im d
= ko ξφ − (rξr ) dr r r dr r
Bz2 ko r2 dr
(20.14)
k2
µ
d
¶2 Now, this expression can be rearranged so that ξz ap-
2 2
+ 2 2 (rξr ) + k ξr pears only in quadrature. Thus, as with ξφ above, the
ko r dr
energy perturbation is minimized here if
Now, the first term is the only one containing ξφ : the
ikr2
µ ¶
perturbed δW can therefore be minimized by choosing d ξr
ξz = 2 (20.19)
im d m + k 2 r2 dr r
ξφ = (rξr ) (20.15)
m2 2 2
+ k r dr When this is chosen, the energy perturbation can be
But then, with this choice, the integrand in (20.13) for written
δW is positive definite for any choice of ξr . Thus, the µ ¶ 2
minimum of δW is positive for k > 0, and approaches dp ξr
δW (r) = 8πr + m2 Bφ2
zero as k → 0: this system is stable for finite wave- dr r2
lengths and approaches marginal stability for very long (20.20)
m2 r2 Bφ2
µ µ ¶¶2
d ξr
wavelengths. + 2
m + k 2 r2 dr r
What is the underlying physics? This equilibrium
has no parallel currents, so current-driven modes (such Still hunting for a minimum of δW , we note that k →
as the pinch) can’t be excited. In addition, as the field ∞ (long wavelengths) kills the last term. Thus, our final
lines are straight, pressure-driven modes (such as the expression to be analyzed for the energy perturbation is
kink) can be excited. Any perturbation to this equi-
πL a
Z µ ¶ 2
librium either bends or compresses the field lines, and dp 2 2 ξr
δW = 8πr + m Bφ rdr (20.21)
both are stabilizing influences. 4π 0 dr r2
114
And now: the magnitude and sign of dp/dr control the requires a particular mix of axial and azimuthal fields.
result. In particular, if In particular, it can be shown that the axial magnetic
field must be made strong enough, and the plasma col-
2 2
dp m Bφ umn fat enough, that no part of the field between the
2r + <0 (20.22) plasma and the wall closes on itself over the length of
dr 4π
the plasma column. This is the Kruskal-Shafranov sta-
anywhere, then one can find a δx function localized in bility criterion, which can be written
this region which will make the full integral negative.
Thus, this equilibrium is unstable against m 6= 0 per- rBz (r) > LBφ (r) (20.25)
turbations. This is the classic kink instabilitey.
If we use the starting equilibrium relation, (20.15), if L is the length of the plasma column. Thus, this im-
the stability condition (20.21) can be rewritten in two poses a severe limit on the current than can be driven
ways: through the plasma.
1 d ¡ 2¢
rBφ < m2 − 1
Bφ2 dr
(20.23)
2r2 d Bφ References
µ ¶
< m2 − 4
Bφ2 dr r I’ve mostly followed a mixture of Priest & Bateman,
here.
Now, for typical Z pinches, Bφ /r is a decreasing func-
tion of r. For such profiles, (20.22) predicts stability for
m ≥ 2. Also, at large radii Bφ ∼ 1/r corresponding
to a vacuum field; but near the origin, in the current-
carrying region, Bφ ∼ r. Thus, a Z pinch is always
instable to the m = 1 mode. Figure 20.7 shows the
physical picture of this instability. Physically, the kink
instability may be attributed to a localized twist in the
column, pushing the field lines together on the inside
edge of the kink, and pulling them apart on the outside
edge. This sets up an internal magnetic pressure gradi-
ent which acts to enhance the kink.
A similar analysis can be done for m = 0 pertur-
bations; it turns out that one can’t assume incompress-
ibility in this case, so it must be treated separately. One
finds an analogous result: the system is unstable if
r dp 2γBφ2 /4π
− > (20.24)
p dr γp + Bφ2 /4π
21. RESISTIVE MHD INSTABILITIES 21.2. Now, consider a finite resistivity. This will allow
the field to diffuse through the fluid (or, vice versa, will
In the previous section we considered ideal MHD in- allow the fluid to move across the field lines). Under the
stabilities – those for which resistivity is not important. right conditions (mainly long wavelength disturbances;
Another very important type of MHD behavior involves see below), the new, reconnected state has lower energy
the effects of resistivity. Resistivity allows plasma to than the initial one. Thus the system is unstable.
move across magnetic field lines. This violates mag-
In this geometry, the unstable tearing mode causes
netic flux conservation (what we called “frozen in field”
magnetic surfaces close to the tearing layer disrupt and
or “flux freezing”). Unlike the ideal case, resistive ef-
reconnect, forming a chain of filaments, as illustrated
fects allow changes in the topology of the field lines.
in Figure 21.2. These filaments are also called “mag-
The most important application of this is in field line
netic islands”, but note this really only describes a two-
reconnection, which we treat in the next chapter.
dimensional slice of the system.
In this chapter we consider resistive instabilities, in
which spontaneous field line readjustment occurs. Re-
sistive instabilities are most relevant in a current sheet.
That is, a thin region, generally formed at the interface
of two different plasmas, whre B varies rapidly. Sev-
eral such instabilities are known; here we consider the
tearing mode, which is fundamental to magnetic recon-
nection.
B. Tearing Mode: the Math they overlap (formally this is an asymptotic analysis),
and this matching allows us to solve for γ and δ. One
I mostly follow Biskamp in these notes; note however useful quantity in this analysis is the jump in the varia-
that the original derivation, given by Furth, Kileen & tion scale of ψ1 at the tearing layer. Its standard, if ugly,
Rosenbluth (1963) is quite different on the surface, and notation is
is often quoted. Both sources do agree on the answer,
1 £ ′
∆′ = lim ψ1 (δ) − ψ1′ (−δ)
¤
however. (21.6)
δ→0 ψ1
We start in the usual place, with the equations of
motion (10.2) and induction (10.4). This instability is Here and in the rest of this section, primes on ψ1 and
well described in the incompressible limit; to shorten φ1 mean derivatives with respect to x; the prime on ∆′
the notation we set ρ = 1, If we linearize the basic is standard notation but does not mean it is a derivative
equations, and ignore both pressure effects and fluid (mutter, mutter . . .). We will consider cases where
viscosity, we get ψ1 varies only slowly with x near the tearing plane, so
that ψ1 . ≪ ψ1 /δ, but where the second derivative ψ1′′ ≃
∂ ψ)1 ∆′ /δ, has a jump at x = 0.
v1 = (∇ × B1 ) × Bo + (∇ × Bo ) × B1
∂t (21.3) To do the actual solution, we work first in the inner,
∂ 2 resistive region (where the dissipation terms are most
B1 = ∇ × (v1 × Bo ) + η∇ B1
∂t important, so we can ignore the inductive terms). We
Assuming incompressibility allows us to introduce a then work in the outer, nearly-ideal region (where we
velocity stream function, φ. Working in two dimen- can ignore the dissipation terms). Finally, we require
sions allows us to use a magnetic flux function, ψ. In that the two solutions match well, which gives us our
this coordinate system, we thus have final answer for the growth rate (or timescale).
• Inner, Resistive Region Now: we pick
v = ẑ × ∇φ ; B = ẑ × ∇ψ (21.4) a wavenumber for the perturbation, ψ1 (x, y) =
ψ1 (x)eiky , and we assume the unperturbed ψ is con-
(Strictly speaking, these are the perpendicular veloc- stant across this region. Keeping the resistive terms,
ity and field components; in the details of the analysis the system (21.5) becomes
we’re ignoring variations in the y direction). Resistive
instabilities are driven mainly by the parallel current: γφ′′1 = ixkB ′ ψ1′′ − ikj ′ ψ1
j ≃ jk ≃ jz . Therefore we also neglect j⊥ , which (21.7)
means Bz ≃ constant. Also, as we’re interested in sta- γψ1 = ixkB ′ φ1 + ηψ ′′
bility and the growth rate, we write ∂f /∂t = γf (so
that γ real and positive corresponds to instability) and It turns out that j ′ term does not affect the resulst much,
make solving for γ the primary goal of our analysis. and it is usually ignored in this analysis. Doing that, we
With this, the equations (21.3) become solve the system by (1) choosing parity: pick ψ1 (x)
even and φ1 (x) odd; (2) approximate derivatives: φ′′1 ≃
γ∇2⊥ φ1 = B · ∇j1 + B1 · ∇j −φ1 /δ 2 , ψ1′′ ≃ −∆′ φ1 /δ. Putting these in and doing
(21.5) some algebra gives us an intermediate solution for the
γψ1 = B · ∇φ1 + η∇2⊥ ψ1 inner region:
Our plan now is (1) pick a geometry for the unper- γ ≃ η 3/5 (∆′ )4/5 (kB ′ )2/5
(21.8)
turbed field; (2) solve these equations for the spatial ′ 2/5 ′ 1/5 ′ −2/5
behavior of ψ1 and φ1 , and (3) from these results deter- δ ≃ δ∆ /γ ≃ η (∆ ) (kB )
mine γ. Remember that dissipation is a second-order Biskamp’s comments here: this argument assumed
derivative term; it will be important only in small re- ∆′ > 0. A full stability analysis (cf. FKR 1963) finds
gions where things change rapidly. In this problem that that instability requires ∆′ > 0, and in fact shows that
means close to the current sheet. The traditional way to −∆′ is the perturbation energy, δW (so that δW < 0
solve the system (21.5) is in two parts. Close to the cur- gives instability). Thus our choice of ∆′ > 0 here is
rent sheet, we keep all the physics, but use simplified OK.
geometry to get an answer. We call the width of this
•Outer, Ideal Region. In this region, we ignore
inner, resistive region δ. Far from the current sheet, we
second derivatives, so the basic equation becomes
use simplified physics (that is, we ignore the resistiv-
ity). We then require the two solutions to match where B · ∇j1 + B1 · ∇j ≃ 0 (21.9)
117
Putting in the geometry of Figure 21.1 explicitly, this The critical point is that this instability leads to a
becomes topological change in the field lines; it is not at all
simple dissipative diffusion. The field lines reconnect
j ′ (x)
µ ¶
′′
ψ1 − k +2
=0 (21.10) across the initial neutral sheet. Some authors discuss
B(x) the tearing mode as a change from one possible equilib-
rium state (say, an ideal one) to another (call it a resis-
Now we must pick a form for B(x). The usual choice tive one). If the new, reconnected state has a lower mag-
for this problem is B(x) = tanh(x/a), which intro- netic energy, then the instability can go spontaneously;
duces a as the width of the current layer. We pick the reduced magnetic energy will appear as heat.
boundary conditions ψ1 → 0 as |x| → ∞. The solution
of (21.10) is now analytic: (a) (b)
µ ¯ x ¯¶
1 10 10
′ 2 1
∆ = − ka
a ka 0 0
-5 -5
Recalling that a full analysis shows instability only
when ∆′ > 0, we see from this that only large scales -10 -10
final result, the growth rate and the physical scale of the (c) (d)
1 1
γ≃ (21.12)
(ka) τ τ 2/5
2/5 3/5 -5 -5
D A
-10 -10
a τA
δ∼ (21.13) Figure 21.3. A 2D simulation illustrating the evolution
(ka)3/5 τD of a reconnecting tearing mode. Solid lines are magnetic
field lines (lines of constant flux function ψ); dotted lines
Thus, we find that such a current sheet is always un- show the velocity field. The initial state (top left) and three
stable to tearing mode on long scales (low k), and that later times are shown; the formation of magnetic islands is
the growth time (τtear ≃ 1/γ) is, indeed, inbetween the apparent. From Garasi, PhD thesis (NMSU), 2002.
Alfven (dynamic) and resistive (diffusion) timescales.
δ (x)
Turbulence is an old topic which remains fresh today. y, v
It was studied as long ago as the 16th century, when
Leonardo da Vinci studied turbulence generated by ob-
x, u
stacles placed in a water flow. It is still being studied
today, when such issues as the role of small-scale vor- Figure 22.1. Flow past a thin plate. An initially uniform
tex ropes, the effect of magnetic fields in MHD turbu- flow comes in from the left; when it encounters the plate,
lence, or “quantum turbulence” in superfluids, are top- the no-slip requirement at the plate surface causes a viscous
ics of current research. In fact, many of the basic ques- boundary layer to develop. The width of the layer, δ(x),
tions are still unanswered today. Because turbulence grows with x in this situation.
is fundamentally nonlinear, analytic solutions are hard
to come by; traditional work relies heavily on scaling a solution of the form u = uo g(y/δ), if uo is the up-
laws, while large numerical simulations are critical to stream/incoming velocity, g is some unknown function,
modern work. and δ can be a function of x. Because this is incom-
To set the stage, I know of no better opening, than pressible, we can use a stream function ψ: u = ∂ψ/∂y;
to paraphrase Tennekes & Lumley. v = −∂ψ/∂x. We therefore look for a ψ solution with
the form ψ(x, y) = uo δf (y/δ); f is another ast-yet-
Most flows occuring in nature and in engineer-
unknown function. If we carry out the algebra and put
ing applications are turbulent. These include
the boundary layer in the earth’s atmosphere; jet these into the basic equations (22.1), we find (a) a new
streams in the upper troposphere; and cumulus ODE; and (b) a consistency condition on δ. The results
clouds which are in turbulent motion. Subsur- are
face ocean currents are turbulent. Stellar atmo-
spheres are turbulent, as is the gaseous interstel- δ(x) ∝ (νx/uo )1/2 ; f f ′′ + f ′′′ = 0 (22.2)
lar medium. Boundary layers on aircraft wings The first is our desired result: the shape of the bound-
are turbulent. Most combustion processes involve
ary layer (the constant of proportionality depends on
turbulence. The flow of natural gas and oil in
just how δ is defined, for instance is it the point where
pipelines is turbulent, as is water flowing in rivers
and canals. THe wakes of ships, cars and aircraft u(y) = 0.9uo ? 0.99uo ?). The ODE has to be solved
are in turbulent motion. In fluid dynamics lami- numerically; the primes represent differentiation with
nar flow is the exception, not the rule: one must respect to the similarity variable, η = y/δ(x).
have small dimensions and high viscosities to en- Now, assume we have the full solution for the
counter laminar flow. boundary layer (rather than just the sketch above). Sub-
ject it to the same type of analytic instability analy-
A. The transition to turbulence sis as we did in chapter 19. That is, assume the per-
turbed (x, y) velocities, and also the perturbed pres-
The instabilities discussed in the last chapter can de-
sure, go as f (y)ei(kx+ωt) ; and determine whether the
velop into full-fledged turbulence. To explore this, let’s
frequency has any imaginary part. If it does, the system
choose a particular setting: consider a boundary layer in
is unstable. The key parameter here is the viscosity,
an initially laminar flow. In order to do stability analy-
ν – which is folded into the boundary-layer Reynolds
sis, we need a mathematical description of the boundary
number, Re = uo δ/ν. Figure 22.2 shows a typical re-
layer. Here’s a sketch of one.
sult: any values of (k, Re) inside the curve are unsta-
A standard representation of such a boundary layer ble, while values of (k, Re) outside the curve are stable.
is the Blasius solution. (I follow Tritton, 2nd edition, Thus: there is a minimum Reynolds number for insta-
chapters 8, 11). Assume incompressible flow, and let u bility – for lower values viscosity stabilizes the flow.
b, the x-velocity, v be the y-velocity. The basic equa- Above this value, there is a finite range of scales which
tions are are unstable: very small and very large scales remain
∂u ∂u ∂2u ∂u ∂v stable.
u +v =ν 2 ; + =0 (22.1) What happens next? From this type of analysis, as
∂x ∂u ∂x ∂x ∂y
well as from experiment, we expect that all flows be-
To solve these, we do dimensional analysis (more for- come turbulent at high enough Re. Linear stability the-
mally, look for a similarity solution). We want to find ory can tell you when it starts .. but the waves pre-
119
δv(l)2 , between two points separated by a distance l, stretching – which increases the magnitude of the vor-
obeys δv(l)2 ∝ l2/3 . Put into wave-number space,1 ticity (wny?), but also reduces the cross-section of the
this becomes what is now known as the Kolmogorov vortex tube.
law: v 2 (k) ∝ k −5/3 .
• Finite energy dissipation: in an experiment, vary
viscosity while keeping everything else the same: find
energy dissipation per unit mass behaves in a way con-
sistent with a finite limit.. Following Frisch, think about
drag. Consider a fluid moving past a rock (or a car mov-
ing through air). The drag force F = (1/2)Cd ρAv 2 ,
(you have seen this, right?) if A = L2 is the projected Figure 22.3. Schematic representation of the evolution
of a “marked blob” of fluid within turbulent motion. Its
area of the rock/car, and CD is the coefficient of drag. shape becomes more and more distorted by the velocity
Experiments show that CD is only a very slow function fluctuations, with smaller and smaller scales appearing, as
of Re. Thus, the work done W = F v; and the energy time goes on. Eventually the scales become small enough
dissipated per unit mass is that diffusion matters, and the marked fluid mixes with its
surroundings.
W 1 v3
ε= 3
= CD (22.3)
ρL 2 L
3. THE KOLMOGOROV SCALING ARGUMENTS
Thus, ε does not depend (very much) on ν, and thus
will have a finite limit as ν → 0. We can talk about the main properties and scaling laws
for homogeneous, isotropic turbulence, following Kol-
2. EDDIES AND THE ENERGY CASCADE mogorov’s analysis, without needing the details of the
statistics. Kolmogorov argued that properties of the
One critical fact about fluid turbulence is that energy is
flow are determined by the scale l, and the energy rate
transferred from large scales to smaller scales. Picture,
ε, only. (motivated by the observations, above). Then,
say, initially large eddies, on the order of the driving
from this, K. argued (or “derived”), the important law:
scale of the turbulence. This could be the width of a
boundary layer, or the diameter of a pipe. Now, con- δv(l)3 = (4/5)εl (22.4)
sider slightly smaller eddies. Due to vortex stretching,
they are strained by the velocity field of the largest ones,
What is the rate of energy transfer in the cascade?
thus growing in strength; they extract energy from the
Let vt = hvi be the mean turbulent speed, and let
larger ones. This continues on to still smaller scales.
λt = hλi be the largest (driving) scale.2 The time for
Thus: the turbulent kinetic energy cascades down from
these largest eddies to “turn over” must be τt ∼ vt /λt .
large to small eddies in a series of small steps. This
Observations find that the large eddies transfer much
process is essentially inviscid, since the vortex stretch-
of their energy to smaller ones, in one or two τt times.
ing mechanism arises from the nonlinear terms in the
Thus the energy cascade rate, which also must become
equations of motion.
the energy dissipation rate, is
Following Tritton: the figure illustrates the evolu-
tion of a blob of fluid, due to the actions of the locally vt3
turbulent velocity field. Showing (1) the process of re- ε≃ (22.5)
λt
peated instability; each stage gives rise to motions of
greater complexity and smaller scales, than the previous Kolmogorov argued that scales in this range should
stages; (2) energy is extracted from larger scales (or the not be affected either by kt = 2π/λt (as eddies at k are
mean flow), to smaller (refer ahead to mean-flow equa- driven only by their immediate neighbors in k-space),
tions and Re stress); (3) vortex stretching is part of the nor by kd = 2π/λd (basically for the same reason; en-
process. The random nature of turbulent motions is dif- ergy is moving downward in k-spacd). Thus, letting
fusive, as two particles that happen to be close, at a give
time, find themselves much further apart soon after. If
2
these two particles are on the same vortex line, we get A formal version of this is the Taylor microscale, defined by
some authors as λ2T = [v 2 ]/[∇v]2 , and by others as λ−2 T =
d2 C(r)/dr2 |0 , that is the curvature of the correlation function at
1
How? change variables, to k ≃ 2π/l; and by energy conserva- zero lag. Either definition has the same content: this is the scale
tion, note that we must have v 2 (l)dl = v 2 (k)dk. where most of the turbulent energy is concentrated.
122
vl be the velocity typical of scale l, the energy flow Finally, it’s worth noting that we can also find the
rate ε = vl3 /l must be independent of the scale l (or dissipation scale by noting that the local energy trans-
k = 2π/l). What power spectrum W (k) is consistent fer rate, vl3 /l, must equal the energy dissipation rate,
with this picture? The constancy of ε tells us that νvl2 /l2 , on this scale. Equating these two recovers
(22.8) nicely.
µ ¶1/3
1/3 2πε Below kt the power falls off, due to a lack of driv-
vl ≃ (εl) ; vk ≃ (22.6) ing, and the fact that turbulent energy only cascades for-
k
ward in hydrodynamic turbulence. Above kd the power
But also, the power “at k” is vk2 ≃ kW (k). Thus we spectrum is usually taken to fall off more steeply, ad
find the equilibrium spectrum, W (k) ∝ k −3 . The data agree with this very well (lots
of people have looked at this).
W (k) ∝ ε2/3 k −5/3 (22.7)
4. WHAT IF THE FLUID IS MAGNETIZED ?
W(k)
−5/3 (direction of
k
energy flow
in the cascade)
D. Appendix: fun facts from Fourier transforms
1
Z
3
Check: can you show this? You need to define Re in terms of the ṽ(k) = v(x)eik·x dx (22.9)
turbulent quantities vt and λt . (2π)3/2
123
4
Note, in going from ṽ(k) to P (k) we have lost information on
the relative phase of the kth mode. This is assumed not to be im-
portant in the usual homogeneous turbulence models, but phase
information does matter in some applications, such as intermit-
tency.
5
Some authors use the structure function:
In the last chapter we talked about the basic nature of This analysis can be repeated for the other basic equa-
turbulence: what it’s like and how it’s described statis- tions; I’m not going to write down all the steps, how-
tically. In this chapter I carry on with more “classical ever. The momentum equation starts as
turbulence”, namely the effect of turbulent stresses on
the mean-field flow; and say a little about “modern tur- ∂
(V + v)+(V + v) · ∇(V + v) =
bulence”, including how 2D turbulence is different, and ∂t
(23.5)
1
turbulence on small scales (which includes vorticity & − ∇(P + p) + ν∇2 (V + v)
intermittency). ρ
A. Mean Field Equations The mean of this equation becomes – written out in
Cartesian
In many systems, we can treat the large-scale flow as ∂Vi ∂Vi 1 ∂P
steady, or at least slowly varying, with the turbulence as +Vj =−
∂t ∂xj ρ ∂xi
a rapidly fluctuating additive term. That is: let V be the à ! (23.6)
mean velocity, and v the turbulent one, so that the net ∂ 2 Vi ∂ 2 Vi ∂
velocity field is V+v. Treat the pressure field similarly, +ν + − hvi vj i
∂x2i ∂x2j ∂xj
P + p; ditto for temperature. (As long as we stay in
the incompressible limit – which is where most turbu- Thus: the turbulent field contributes a net reaction
lence analysis stays – there are no density fluctuations, back on the mean field. This last term is called the
right?) If these sums are put into the basic dynam- Reynolds stress. Even though hvi i = 0 for each com-
ical equations, we can (borrowing terminology from ponent i, the mean of their product does not necessar-
MHD turbulence), isolate the dynamical equations for ily vanish: hvi vj i =
6 0. It turns out that this is so for
the “mean field” quantities, and find how the turbulence anisotropic turbulence – and real turbulence is com-
affects the mean flow. We will find, for instance, that monly anisotropic.1
the mean-flow momentum equation contains what are Note that the Reynolds stress arisews from averag-
called Reynolds stresses: non-zero terms involving sec- ing the nonlinear convective term in the full Navier-
ond moments of the fluctuating velocity field. Thus, Stokes equation; it is not truly a stress term, more like
the mean flow and turbulence are intimately connected, the mean momentum flux due to the turbulent fluctua-
with the one affecting the other. tions. Nonetheless, it is often talked about in terms of a
1. THE CONTINUITY EQUATION
turbulent viscosity. Because this extra force term in the
basic momentum equation is the gradient of a tensor, it
For instance, take the incompressible continuity equa- is generally combined with the viscous stress tensor, as
tion: µ ¶
∂Vi ∂Vj
∇ · (V + v) = ∇ · V + ∇ · v = 0 (23.1) σij = −P δij + ρν + − ρhvi vj i (23.7)
∂xj ∂xi
Now, take means: that means ensemble averages, or We can go further, and write the turbulent contribution
time averages for the case of stationary flows. This to the viscous stress tensor as
gives µ ¶
dVi dVj
∇ · hVi + ∇ · hvi = 0 (23.2) σij,turb = −ρhvi vj i = ρνt + (23.8)
∂xj ∂xi
But now: first, we have hVi = V; that is the “mean This last step is an assumption which turns out to be
flow” velocity. And, we assume the turbulent fluctua- quite well justified. That is, the Reynolds stress is lin-
tions have zero mean: hvi = 0. Thus, (23.2) becomes early proportional to the mean flow shears. But then,
∇·V =0 (23.3)
1
(recovering incompressibility of our mean state). And, This is the case because, microscopically, viscosity mixes the rel-
subtracting (23.2) from (23.1) gives ative phases of each vi component; they do not stay in phase, and
thus the mean of their product is non-zero. Alternatively: think
∇·v =0 (23.4) about how Brownian motion, on the microscopic level, trans-
ports momentum components laterally.. this gives us the usual
showing that our turbulent field is also incompressible macro stress tensor. Here, the turbulent fluctuations provide a
(by itself). mean transport of momentum...
125
we can see from this that the turbulent viscosity coef- the constancy of τ we get τ = ydPw /dx + τw and
ficient νt ∼ vt L, if L is the local gradient scale. This τw is another constant, the mean viscous shear stress
connects back to our dimensional argument, above. at the wall y = 0. We find its value by noting that
the flow must be symmetric about the midplane; thus
3. EXAMPLE : 2D CHANNEL FLOW τw = −DdPw /dx. Using this in the definition of τ , we
I follow Mathieu & Scott here. Consider a channel get
in the x-direction, with y the transverse coordinate. ³ y´ dVx
(Once again, compare this to the laminar flows in chap- τw 1 − + ρhvx vy i = ν (23.11)
D dy
ter 2). Assume the mean flow independent of z, and
that Vz = 0. Note, this does not imply that the flow This is almost the answer for the mean flow .... as we’re
is fully 2D, indeed the turbulence will be 3D. Rather, turbulent, we know Re ≫ 0, so that the RHS is small
we’re assusming that it’s 2D in the mean, and all sta- except very close to the walls. Thus, away from the
tistical properties are unchanged under reflection about boundary layer, −ρhvx vy i ≃ τw (1 − y/D), that is lin-
the z = 0 plane. As before, put the overall mean flow in ear behavior. And, if this balance holds, we expect Vx
x direction, and let it be driven by some dP/dx. Also to be independent of y; and “measurements find Vx is
as before, all quantites except the pressure depend only approximately constant in this region”. Near the wall,
on y; we extend this here to include Re stresses. we expect steep gradients in Vx , so that viscosity mat-
111111111111111111111111111111111111111111
000000000000000000000000000000000000000000
000000000000000000000000000000000000000000
111111111111111111111111111111111111111111
ters...and a more complicated behavior of the Reynolds
stress.
y
dhvy2 i and noting that δij ∂Vi /∂xj = ∂Vi /∂xi = 0, we have a
∂P
+ρ =0 useful form:
∂y dy
(23.9) µ ¶ µ ¶
d 2 Vx ∂P dhvx vy i D 1 2 ∂ P Vj
ν 2 = +ρ V = − + 2νVi Eij − hvi vj iVi
dy ∂x dy Dt 2 i ∂xj ρ
∂Vi (23.14)
From the first of these, we get P (x, y) = Pw (x) − − 2νEij Eij + hvi vj i
∂xj
ρhvy2 i: where Pw (x) is the pressure at the wall (where
vy = 0 by no-slip); and the second term is the diagonal Here, the first 3 terms (inside the brackets) are advective
contribution of the Re stress to the net pressure. Putting energy transport; the term with ν is viscous dissipation;
this into the second equation in (23.9), we get and the last term is energy lost (in a decelerating flow)
· ¸ to drive the turbulence.
dPw dτ d ∂Vx
= = ν − ρhvx vy i (23.10)
dx dy dy ∂y 5. WHAT ABOUT THE TURBULENT TERMS ?
The total mean stress is defined as τ = −ρhvx vy i + In these notes I’ve said nothing about the dynamics of
νdVx /dy. But now: the left hand side of (23.10) is the turbulent terms (p, v, etc). That was on purpose
only a function of x, while the right hand side is only .... refer back to the scale separation we used in de-
a function of y; thus, each side is a constant. The left riving the mean-field equations (23.3), (23.6), (23.14).
hand side is just the driving pressure gradient. From We could equally well have separated out the D/Dt
126
terms involving the turbulent fields. Remember, now, scales: one author described it as the buildup of coher-
that flow equations are nonlinear (the v · ∇v term, for ent vortices, in which nonlinear distortions nearly van-
instance). When one isolates the dynamical equations ish, and which continue to grow by vortex coalesence.
for the fluctuations, one finds higher-order terms are in- One can use arguments similar to those in §22.3 to
volved. For instance, the equation for Dhvi i/Dt in- predict the turbulent spectrum. Let ko be the driving
volves ∂hvi vk i/∂xk ; the equation for Dhvi vj i/Dt in- wavenumber, and assume it is well between the sys-
volves ∂hvi vj vk i/∂xk ; and so on. Some additional as- tem scale ks , and the dissipation scale kd . Energy cas-
sumption is always needed to close a system like this; cades to lower wavenumbers; the arguments of §22.3
and pursuing that would take us too far afield. still apply, so that a power spectrum W (k) ∝ k −5/3
should obtain for kmin < k < ko , where kmin cor-
B. Two-dimensional Turbulence
responds to some largest scale reached by the reverse
. . . differs significantly from turbulence in three- cascade. Above ko , the enstrophy cascade rules. De-
dimensions. The most striking difference is that tur- scribing enstrophy and turnover time at scale l as ωl ∼
bulent power in 2D cascades both forward (to higher vl2 /l2 , τl ∼ l/vl , we get the enstrophy power to be
wavenumber) and backwards (to lower wavenumber). ǫl ∼ Ωl /τl ∼ vl3 /l3 . But this last must be indepen-
This is still a topic of active discussion in the literature. dent of l (repeating the Kolmogorov argument); from
I mostly follow Biskamp’s discussion in these notes. this we need vl ∝ l, so that kW (k) ∝ vk2 ∝ k −2 . This
should describe the cascade up to the dissipation range.
Recall 3D turbulence: we saw that turbulent power
Thus, for 2D fluid turbulence driven at ko , we expect
undergoes only a forward cascade. We also recall that
three-dimensional fluid flow has two significant invari-
ants:
1 W (k) ∝ k −5/3 ;
Z
kmin < k < ko
E= v 2 dV, energy (23.15) (23.19)
2 −3
W (k) ∝ k ; ko < k < kd
and
1
Z
HV = v · ∇ × vdV, velocity helicity (23.16)
2 This is supported by observations. As in 3D turbulence,
(Invariance of the first should be obvious; invariance of the high-k part of the cascade can be steady-state. En-
the second may be proved in the homework.2 ) Both en- ergy input to the system at ko will cascade forward,
ergy and velocity helicity are conserved in the mode- reaching the dissipation scale kdiss where it goes into
mode interactions which set up the cascade. Such hear. On the contrary, however, the low-k part of the
small-scale interactions drive both E and HV forward, cascade cannot be stable. There is no low-k dissipation
to smaller wavenumbers. The second invariant in 2D to balance the energy input. One would expect kmin to
turbulence, however, is different. The two invariants in decrease with time, until the system size ks is reached;
2D are at and after this point energy will continue to accumu-
late at the largest scales allowed in the system.
1
Z
E= v 2 dV, energy (23.17) What sets the cascade direction? I have not found
2 any clear answer in the literature. On small scales,
and the nature of energy, enstrophy and helicity transfer in
mode-mode interactions (say, k1 +k2 → k3 ) must have
1
Z
Ω= ω 2 dV, enstrophy (23.18) a preference for forward or reverse transfer. The details
2 of these processes seem not to be obvious, however.
Biskamp reports on work in the literature addressing
(The enstrophy is a fancy name that turbulence types
the (thermodynamic) equilibrium distribution of E(k),
like to use for the mean square vorticity.) In 2D tur-
HV (k), and Ω(k). It appears that the equilibrium dis-
bulence, wave-wave interactions drive Ω forward, but
tributions of both E(k) and HV (k) are both weighted
drives E to smaller wavenumbers – in an inverse cas-
towards high k’s in 3D; while in 2D the distribution of
cade. This cascade drives power to larger and larger
Ω(k) is weighted to high k, but that of E(k) is weighted
to low k. The inference, then, is that the cascade di-
2
The methods are similar to the proof that magnetic helicity is rection is set by the tendency of the system to move
invariant, back in chapter 15. towards a statistical equibrium state.
127
C. Small scales and intermittency time/space scales. This turns out not to be so: on
small enough scales (high enough frequencies, think of
Much of the action these days seems to be “what’s hap- a high-pass filter), you find intermittent “bursts” of en-
pening on small scales”, which means close to the Kol- ergy. Connect this to statistics: a self-similar random
mogorov dissipation scale. This is potentially a big signal, which is the same on different scales, is a white
topic; I’m storing only a brief overview here. noise signal, and has Gaussian statistics. When you go
intermittent, you have a greater chance of getting large
1. INTERMITTENCY bursts, and less chance of getting low-amplitude signal
.. so the probability distribution function (PDF) flat-
The term intermittency is used in two different ways — tens. Several authors show PDF plots with tails much
which drove me crazy when I was trying to learn this flatter than Gaussian. This becomes conspicuous only
field. on scales comparable to, or smaller than, the dissipa-
Older: macroscopic intermittency. Turbulence tion scale. Thus, it is characteristic of the dissipation
observed on large scales (comparable to the system range, and does not imply breakdown of the entire K41
size, well above the dissipation range) is intermittent: analysis.
only a fraction of the volume is filled with turbulent
spots or eddies at any instant. Sit at one point in a tur-
bulent region .. you will find periods of high frequency 3. THE ROLE OF VORTEX FILAMENTS
References
The mean-field material is “traditional”, and can be
found in various books which treat turbulence math-
ematically. Tennekes & Lumley, or Hinze, are good
sources.
The newer turbulence material, especially on inter-
mittency, I’ve taken from “here and there”. The An-
nual Review of Fluid Mech has several useful reviews
(including the Moffat homage-to-Batchelor article from
2003).
129
acquired. It is now established that solar wind turbu- modifications in the Kraichnan and GS models. Note
lence is anisotropic; that’s good, because the wind def- that this allows a steady state spectrum W (k): the for-
initely has an ordered background B. Unfortunately, ward cascade transfers energy from the input/driving
however, newer data don’t clearly support the IK spec- scale, to the dissipation scale, at which scale it goes
trum – they suggest a steeper spectrum, maybe closer to heat. Because dissipation is more important at small
to k −5/3 . I personally suspect that we just don’t know, scales (high k’s), we expect always to be able to find a
yet. steady state.
Moving to a larger playing field, there is a much- For 3D-HD turbulence, we only have a forward cas-
quoted paper (Armstrong etal ApJ 1995) which at- cade. But in 2D HD turbulence, and in 2D and 3D
tempted to pull together several different types of mea- MHD turbulence, ond also finds an inverse cascade:
surement – all indirect, of necessity – of electron den- power flows to larger scales (smaller wavenumbers).
sity fluctuations in the interstellar medium (ISM). I Inverse cascades move power to large scales, where dis-
have to point out that the ISM is very definitely a mag- sipation is always small. Thus, they do not reach steady
netized plasma, but also very definitely highly inho- states. They tend to evolve to lower and lower k’s, until
mogeneous – different experiments probe different re- they reach the scale of the system, and then accumulate
gions and scales of the galaxy. This paper patched all power on those scales.
the measurements together and concluded that the ISM To summarize:
obeys a Kolmogorov-type scaling, W (k) ∝ k −5/3 ,
• 3D HD: all cascades are direct. This is the basic
over 12(!!) orders of magnitude in k. I remain person-
Kolmogorov picture.
ally skeptical of this result, because of the difficulties
and uncertainties in combining these disparate types of • 2D HD: enstrophy has a forward cascade, but total
measurement. energy has an inverse cascade. This seems to connect
to the coherence and large-scale eddies seen in 2D tur-
C. Cascades and Related Things bulent flows.
• 3D MHD: energy and cross helicity have forward
In addition to MHD effects on the inertial range, there cascades; while helicity has an inverse cascade. This
are several other important differences. has dramatic consequences for maintaining large-scale
1. IDEAL INVARIANTS
ordered magnetic fields; this is the turbulent dynamo
which we see in the next section.
How to HD and MHD turbulence compare, in 2D and • 2D MHD: energy and cross helicity have forward
3D? What are the important invariants? cascades; while the mean vector potential has an inverse
We saw in Chapter 23 that magnetic helicity is in- cascade.
variant in MHD flows; as long as dissipation is small
(for scales above the dissipation range), total energy is 3. SELF - ORGANIZATION IN MHD
as well. This analysis can be carried out for HD and
MHD turbulence. One finds the following invariants. Some interesting properties of MHD turbulence are
• 3D HD: total energy E = 12 v 2 dV ; velocity he-
R connected to these invariants.
licity HV = 12 v · ∇ × vdV • Force-free fields. We have already seen that
R
turbulent relaxation evolves a plasma towards a self-
• 2D HD: total energy E; and total vorticity, also
~ = 1 ω 2 dV (recall ω organized state, in which j k B. This is due to the
called enstrophy, Ω 2 ~ = ∇ × v).
invariance of the magnetic helicity (that is, to its much
1
R • 3D MHD: total energy E; magnetic1 Rhelicity HB = slower rate of decay than the magnetic energy shows).
A · BdV ; and cross helicity, K = 2 v · BdV
2 • Dynamic alignment. There is evidence (in
• 2D MHD: total energy E; cross helicity
R K; and the the solar wind, backed up by dynamical models) that
mean square magnetic potential, A = 21 ψ 2 dV strong MHD turbulence shows a second form of self-
organization. It tends to evolve to a state in which
2. CASCADE DIRECTIONS
±v k B. That is, it shows dynamic alignment of the
In 3D HD turbulence we encountered a forward cas- velocity and magnetic fields. This is due to the invari-
cade (also called a direct cascade). That is, the tur- ance of the cross-helicity; a formal minimization of E,
bulent shear stresses transfer energy to smaller scales subject to K being constant, derives this condition.
(higher wavenumbers). The original Kolmogorov pic- • Energy equipartition. We should also note that
ture involved a forward cascade, of course, as does its there is evidence that MHD turbulence reaches a state
132
of approximate equipartition between kinetic and mag- sun/earth/galaxy formed — because we know the re-
netic energies: hvi2 ∼ hB 2 i/4π. This is a common sistivity of the plasmas in question, and thus we know
result in numerical simulations how long it would take a primordial current to dissipate.
Such calculations predict that primordial fields would
D. MHD Dynamos long ago have died away; but we know that stars, plan-
ets, and galaxies are still magnetized.3
Now, a larger question: where do magnetic fields come
from? In the lab, the answer is easy: “currents”. In Thus, we still must ask, “what supports astrophysi-
magnetic solids, the currents are those of well-ordered cal B fields? The answer is still, “currents”; but what
electrons spins in ferromagnetism. More typically, cur- drives astrophysical currents? We can’t expect a device
rents in the lab — and their consequent B fields — such as in Figure 24.1 exists inside a planet, or star, or
come from obvious things like batteries and wires. The whatnot ... so we need to find a way to drive fluid mo-
issue is then, what drives the currents? My dictionary tions which can maintain the B fields we observe. This
defines a dynamo as question gets us into what’s called dynamo theory.
As with vorticity, the question of fluid dynamos
“a device for converting mechanical en- could occupy a full course on its own. All we can do,
ergy into electrical energy, usually by ex- once again, is a brief introduction. To set the stage,
pending the mechanical energy in produc- I paraphrase Roberts & Soward (1994, Ann Rev Astr
ing a periodic motion of a conductor and Ap):
a surrounding magnetic field”.
We know dynamos exist; one can buy them in
shops.4 What, one may then ask, is all the fuss
A simple lab version of this is called the unipolar about? To answer, we must make an impor-
dynamo, in Figure 24.1. This involves a conducting tant distinction between commercial dynamos and
disk, threaded by a B field, which rotates about its axis. naturally occuring ones. In order to ensure that
This induces a radial E field, v×B/c, and thus a poten- the induced currents don’t short-circuit, an engi-
tial drop between the axis and the edge of the disk. If neer makes sure that the geometry of the machine
you hook up wires in the right way you’ll have a current is an asymmetric and multiply-connected region.
– and this current will create its own magnetic field. By contrast, the conducting fluid in as astrophys-
ica body usually occupies a symemtric, simply-
connected domain – such as a sphere (the earth
or sun) or a simple disk (the galaxy). It is by no
means obvious that a dynamo can operate in such
a simple system. Than, then, is the point of this
field [by which they mean dynamo theory].
1. COWLING ’ S THEOREM
ηk 2 vo2 ∂b
hv × bi = (Bo · k) ẑ (24.12) =∇ × (V × bv × B)
(η 2 k 4 + ω 2 ) ∂t (24.16)
2
+ ∇ × G + η∇ b
This verifies our assumption, that hv×bi belongs in the
large-scale induction equation, (24.9). Thus, this EMF where
acts back on the large-scale field, in fact acts as a source
(driving) term on that field. G = v × b − hv × bi (24.17)
135
Now: the standard analysis argues that the ∇ × G – which connects back to Taylor relaxation, which we
term can be dropped from (24.16) (due to being second- saw in Chapter 15: a plasma will spontaneously relax
order small), as can the dissipative term η∇2 b (assum- to a mininum-energy state, which is force-free. The re-
ing we are still far from the disspation lengths). This laxation proceeds through self-generated turbulence –
means that we can formally integrate (24.16) to get thus the plasma finds its own solution of ∇ × (αB) +
Z η∇2 B = 0.
b ≃ ∇ × [v × B] dt′ (24.18)
G. Astrophysical dynamos in the lab
(where v′ is evaluated at t′ , and unprimed means eval- Finally, a few words about trying to do this in the lab.
uated at t. From this, Everything above is still pure theory — it would be
good to verify directly that an αω dynamo (rotation plus
turbulence, as in the sun), or an α2 dynamo (pure tur-
Z
ε ≃ − dt′ hv × ∇ × v′ × B i
£ ¡ ¢¤
(24.19)
bulence) can really make a large-scale ordered B field.
Several groups are working on this, including NMT’s
But now, we can note that ε is linear in B and its curl: very own Stirling Colgate. The experiments use liquid
metal — usually liquid sodium — in some sort of ro-
ε ≃ αB + β∇ × B (24.20) tating system (the ω in an αω dynamo), and try various
ways to induce turbulence (the α) in the flow. The last I
Working out the integrals in detail, one finds
heard, no one had sucessfully made their dynamo work
1 τ — but I think the field’s progressing. Check the Feb
Z
α= dt′ hv · ∇ × v′ i = HV 2006 issue of Physics Today if you’d like more details.
3 3
(24.21)
1 τ
Z
β= dt′ hv · v′ i = hv2 i
3 3
if τ is the velocity correlation time. Thus: the mean References
field equation now becomes For the IK/GS turbulent cascade discussion, I’ve
∂B followed two excellent recent review papers: Diamond
= ∇ × (V × Bo ) + ∇ × (αB) + (η + β)∇2 B etal, A Tutorial on basic concepts in MHD turbu-
∂t
(24.22) lence and turbulent processes (available through the
Thus: the β term simply adds to the dissipation – an- physics.ucsd web site, I can’t find the publication ref-
other meeting with turbulent dissipation. The α term, erence); and Schekochihin & Cowley, Turbulence and
however, acts as a source term: one can show that it magnetic fields in astrophysical plasmas, to appear
leads to a growth of magnetic flux. in MHD: historical evolution and trends (Molokov,
Thus, we have a big result: helical turbulence in a Moreau & Moffatt, eds, 2006). These papers have ref-
plasma can amplify the large-scale field. That’s a tur- erences to the original papers.
bulent dynamo. I’ve also used Biskamp’s book for discussions on
If this works – if ε = αB — then we can see two intermittency, inverse cascades & invariants, etc.
useful astrophysical consequences. One is balancing The general dynamo material is in lots of books;
ohmic losses, as in the sun or the earth — and (in prin- I like Priest’s Solar MHD, and also Moffett’s original
ciple) accounting for the occasional field reversals in Mean Field Electrodynamics.
each body. The second is “growing” the B field in the
first place. To see this, note that (9.14) allows solutions
B ∝ eαt , if α is constant in time. That’s a growing B
field, with growth time ∼ L/α (some large-scale length
scale L). We might expect that a small seed field would
grow exponentially until some other physics (dissipa-
tion? back reaction on the driving fluid?) comes into
play.
In addition, we note that when V = 0 (no large-
scale flows), (24.14) has particularly simple steady-
state solutions: B k j. That is just a force-free field