a P f.
(8.22)
Now using the basic kinematic relation (19), we arrive at the following Euler equation:18
Euler
v
equation
v v P f .
(8.23)
t
Generally, this equation has to be solved together with the continuity equation (21) and the equation of state of the particular fluid, = ( P ). However, as we have already discussed, in many 16 If the reader still needs a reminder, see MA Eq. (12.1).
17 See, e.g., EM Sec. 4.1, QM Sec. 1.4, and SM Sec. 5.6.
18 It was derived in 1755 by the same Leonhard Euler whose name has already been (reverently) mentioned several times in this course.
Chapter 8
Page 9 of 30
CM: Classical Mechanics
situations the compressibility of water and other important liquids is very low and may be ignored, so that may be treated as a given constant. Moreover, in many cases the bulk forces f are conservative and may be represented as a gradient of a certain potential function u(r) – the potential energy per unit volume:
f
u ;
(8.24)
for example, for a uniform, vertical gravity field, u(r) = gy, where y is referred to some (arbitrary) horizontal level. In this case, the right-hand side of Eq. (23) becomes -( P + u). For these cases, it is beneficial to recast the left-hand of that equation as well, using the following well-known identity of vector algebra19
2
v
v v
v v
.
(8.25)
2
As a result, the Euler equation takes the following form:
2
v
v
v v P u
0.
(8.26)
t
2
In a stationary flow, the first term of this equation vanishes. If the second term, describing fluid’s vorticity, is zero as well, then Eq. (26) has the first integral of motion,
Bernoulli
2
P u v const ,
(8.27) equation
2
called the Bernoulli equation.20 Numerous examples of the application of Eq. (27) to simple problems of stationary flow in pipes, both with and without the Earth gravity field, should be well known to the readers from their undergraduate courses, so I hope I can skip their discussion without much harm.
In the general case, an ideal fluid may have vorticity, so that Eq. (27) is not always valid.
Moreover, due to the absence of viscosity in an ideal fluid, the vorticity, once created, does not decrease along the so-called streamline – the fluid particle’s trajectory, to which the velocity is tangential at every point.21 Mathematically, this fact22 is expressed by the following Kelvin theorem: (v) dA = const along any small contiguous group of streamlines crossing an elementary area dA.23
However, in many important cases, the vorticity is negligible. For example, even if the vorticity exists in some part of the fluid volume (say, induced by local turbulence, see Sec. 6 below), it may decay due to the fluid’s viscosity, to be discussed in Sec. 5, well before it reaches the region of our interest. (If this viscosity is sufficiently small, its effects on the fluid’s flow in the region of interest are 19 It readily follows, for example, from MA Eq. (11.6) with g = f = v.
20 Named after Daniel Bernoulli (1700-1782), not to be confused with Jacob Bernoulli or one of several Johanns of the same famous Bernoulli family, which gave the world so many famous mathematicians and scientists.
21 Perhaps the most spectacular manifestation of the vorticity conservation is the famous toroidal vortex rings (see, e.g., a nice photo and a movie at https://en.wikipedia.org/wiki/Vortex_ring), predicted in 1858 by H. von Helmholtz, and then demonstrated by P. Tait in a series of spectacular experiments with smoke in the air. The persistence of such a ring, once created, is only limited by the fluid’s viscosity – see the next section.
22 This theorem was first formulated (verbally) by Hermann von Helmholtz.
23 Its proof may be found, e.g., in Sec. 8 of L. Landau and E. Lifshitz, Fluid Mechanics, 2nd ed., Butterworth-Heinemann, 1987.
Chapter 8
Page 10 of 30
CM: Classical Mechanics
negligible, i.e. the ideal-fluid approximation is still acceptable.) Another important case is when a solid body of an arbitrary shape is embedded into an ideal fluid whose flow is uniform (meaning, by definition, that v(r, t) = v0 = const) at large distances,24 its vorticity is zero everywhere. Indeed, since
v = 0 at the uniform flow, the vorticity is zero at distant points of any streamline, and according to the Kelvin theorem, should equal zero everywhere.
In such cases, the velocity distribution, as any curl-free vector field, may be represented as a gradient of some effective potential function,
v
.
(8.28)
Such potential flow may be described by a simple differential equation. Indeed, the continuity equation (21) for a steady flow of an incompressible fluid is reduced to v = 0. Plugging Eq. (28) into this relation, we get the scalar Laplace equation,
2
0 ,
(8.29)
which should be solved with appropriate boundary conditions. For example, the fluid flow may be limited by solid bodies, inside which the fluid cannot penetrate. Then the fluid velocity v at the solid body boundaries should not have a normal component; according to Eq. (28), this means
0 .
(8.30)
surfaces
n
On the other hand, if at large distances the fluid flow is known, e.g., uniform, then:
v const,
r
at .
(8.31)
0
As the reader may already know (for example, from a course on electrodynamics25), the Laplace equation (29) is analytically solvable in several simple (symmetric) but important situations. Let us consider, for example, the case of a round cylinder, with radius R, immersed into a flow with the initial velocity v0 perpendicular to the cylinder’s axis (Fig. 8). For this problem, it is natural to use the cylindrical coordinates, with the z-axis coinciding with the cylinder’s axis. In this case, the velocity distribution is obviously independent of z, so that we may simplify the general expression of the Laplace operator in cylindrical coordinates26 by taking / z = 0. As a result, Eq. (29) is reduced to27
1
1
2
,
0
at .
R
(8.32)
2
2
The general solution of this equation may be obtained using the variable separation method, similar to that used in Sec. 6.5 – see Eq. (6.67). The result is28
24 This case is very important, because the motion of a solid body, with a constant velocity u, in the otherwise stationary fluid, gives exactly the same problem (with v0 = -u), in a reference frame bound to the body.
25 See, e.g., EM Secs. 2.3-2.8.
26 See, e.g., MA Eq. (10.3).
27 Let me hope that the letter , used here to denote the magnitude of the 2D radius vector = { x, y}, will not be confused with the fluid’s density – which does not participate in this boundary problem.
28 See, e.g., EM Eq. (2.112). Note that the most general solution of Eq. (32) also includes a term proportional to
, but in our geometry, this term should be zero for such a single-valued function as the velocity potential.
Chapter 8
Page 11 of 30
CM: Classical Mechanics
a b ln
n
c
n s
n a b n
(8.33)
0
0
cos
sin
n
n
n
n
,
n 1
where the coefficients an and bn have to be found from the boundary conditions (30) and (31). Choosing the x-axis to be parallel to the vector v0 (Fig. 8a), so that x = cos , we may spell out these conditions in the following form:
,
0
at R,
(8.34)
v cos , at R,
(8.35)
0
0
where 0 is an arbitrary constant, which does not affect the velocity distribution and may be taken for zero. The condition (35) is incompatible with any term of the sum (33) except the term with n = 1 (with s 1 = 0 and c 1 a 1 = – v 0), so that Eq. (33) is reduced to
c b
1 1
v
cos.
(8.36)
0
Now, plugging this solution into Eq. (34), we get c 1 b 1 = – v 0 R 2, so that, finally, 2
R
v
cos.
(8.37a)
0
y
(a)
(b)
2v
r
0
v0
v
0
0
x
2v0
Fig. 8.8. The flow of an ideal, incompressible fluid around a cylinder: (a) equipotential surfaces and (b) streamlines.
Figure 8a shows the surfaces of constant velocity potential given by Eq. (37a). To find the fluid velocity, it is easier to rewrite that result in the Cartesian coordinates x = cos, y = sin: 2
2
R
R
v x 1
v x 1
.
(8.37b)
0
2
0
2
2
x y
Chapter 8
Page 12 of 30
CM: Classical Mechanics
From here, we may readily calculate the Cartesian components of the fluid’s velocity: 29
2
2
2
y x
R
v
v 1
2
x
R
v
0
x
2
2
x y
1
cos 2
,
2
0
2
(8.38)
xy
R
2
2
2
v
v R
v
y
0
y
2 2
x y
sin 2 .
2
0
2
These expressions show that the maximum fluid’s speed is achieved at the transverse diameter’s ends (
= R, = /2), where v = 2 v 0, while at the longitudinal diameter’s ends ( = R, = 0, ), the velocity vanishes – the so-called stagnation points.
Now the pressure distribution may be calculated by plugging Eqs. (38) into the Bernoulli equation (27) with u(r) = 0. The result shows that the pressure reaches its maximum at the stagnation points, while at the ends of the transverse diameter x = 0, where the velocity is largest, it is lower by 2 v 2
0 . Note that the distributions of both the velocity and the pressure are symmetric with respect to the transverse axis x = 0, so that the fluid flow does not create any net drag force in its direction. It may be shown that this result, which stems from the conservation of the mechanical energy of an ideal fluid, remains valid for a solid body of arbitrary shape moving inside an infinite volume of an ideal fluid – the so-called D’Alembert paradox. However, if a body moves near an ideal fluid’s surface, its energy may be transformed into that of the surface waves, and the drag becomes possible.
Speaking about the surface waves: the description of such waves in a gravity field30 is one more classical problem of the ideal fluid dynamics.31 Let us consider an open surface of an ideal liquid of density in a uniform gravity field f = g = - gn y – see Fig. 9.
y
q( y ,
0 z, t)
P P 0
k
A
v( y ,
0 z, t)
0
z
A
Fig. 8.9. Small surface wave on a deep
heavy liquid. Dashed lines show particle
trajectories.
(For
clarity,
the
displacement amplitude A is strongly
g
exaggerated.)
If the wave amplitude A is sufficiently small, we may neglect the nonlinear term (v)v A 2 in the Euler equation (23) in comparison with the first term, v/ t, which is linear in A. For a wave with 29 Figure 8b shows the flow streamlines. They may be found by the integration of the obvious equation dy/ dx =
vy( x, y)/ vx( x, y). For our simple problem, this may be done analytically, giving y(1 – R 2/2) = const, where the constant is specific for each streamline.
30 The alternative, historic term “gravity waves” for this phenomenon may nowadays lead to confusion with the relativistic effect of gravity waves – which may propagate in free space.
31 It was solved by Sir George Biddell Airy (1801-19892), of the Airy functions’ fame. (He was also a prominent astronomer and, in particular, established Greenwich as the prime meridian.)
Chapter 8
Page 13 of 30
CM: Classical Mechanics
frequency and wave number k, the particle’s velocity v = dq/ dt is of the order of A, so that this approximation is legitimate if 2 A >> k( A)2, i.e. when
kA ,
1
(8.39)
i.e. when the wave’s amplitude A is much smaller than its wavelength = 2/ k. Due to this assumption, we may neglect the liquid vorticity effects, and (for an incompressible fluid) again use the Laplace equation (29) for the wave’s analysis. Looking for its solution in the natural form of a sinusoidal wave, uniform in one of the horizontal directions ( x),
i( kz
)
Re
t
( y) e
,
(8.40)
we get a very simple equation
2
d
2
k 0 ,
(8.41)
2
dy
with an exponential solution (properly decaying at y -), = A exp{ ky}, so that Eq. (40) becomes
Re ekyei( kz t) eky cos
,
(8.42)
A
A
kz t
where the last form is valid if A is real – which may be always arranged by a proper selection of the origins of z and/or t. Note that the rate k of the wave’s decay in the vertical direction is exactly equal to the wave number of its propagation in the horizontal direction – along the fluid’s surface. Because of that, the trajectories of fluid particles are exactly circular – see Fig. 9. Indeed, using Eqs. (28) and (42) to calculate velocity components,
v ,
0
v
k eky cos
(8.43)
x
y
A
kz t,
v
k
eky sin
z
A
kz t,
y
z
we see that vy and vz, at the same height y, have equal real amplitudes, and are phase-shifted by /2. This result becomes even more clear if we use the velocity definition v = dq/ dt to integrate Eqs. (43) over time to recover the particle displacement law q( t). Due to the strong inequality (39), the integration may be done at fixed y and z:
ky
ky
k
q q e sin kz t
,
q q e cos kz t
,
q
with
.
(8.44)
y
A
z
A
A
A
Note that the phase of oscillations of vz coincides with that of qy. This means, in particular, that at the wave’s “crest”, particles are moving in the direction of the wave’s propagation – see arrows in Fig. 9.
It is remarkable that all this picture follows from the Laplace equation alone! The “only”
remaining feature to calculate is the dispersion law ( k), and for that, we need to combine Eq. (42) with what remains, in our linear approximation, of the Euler equation (23). In this approximation, and with the bulk force potential u = gy, the equation is reduced to
P gy
0 .
(8.45)
t
Chapter 8
Page 14 of 30
CM: Classical Mechanics
This equality means that the function in the parentheses is constant in space; at the surface, and at negligible surface tension, it should be equal to the pressure P 0 above the surface (say, the atmospheric pressure), which we assume to be constant. This means that on the surface, the contributions to P that come from the first and the third term in Eq. (45) have to compensate for each other. Let us take the average surface position for y = 0; then the surface with waves is described by the relation y( z, t) = qy( y, z, t) – see Fig. 9. Due to the strong relation (39), we can use Eqs. (42) and (44) with y = 0, so that the above compensation condition yields
k
sin
.
(8.46)
A
kz t g sin
A
kz t 0
This condition is identically satisfied on the whole surface (and for any A) as soon as Surface
waves’
2
gk ,
(8.47)
dispersion
This equality is the dispersion relation we were looking for. Looking at this very simple result (which includes just one constant, g), note, first of all, that it does not involve the fluid’s density. This is not too surprising, because due to the weak equivalence principle, particle masses always drop out from the solutions of problems involving gravitational forces alone. Second, the dispersion law (47) is strongly nonlinear, and in particular, does not have an acoustic wave limit at all. This means that the surface wave propagation is strongly dispersive, with both the phase velocity u ph / k = g/ and the group velocity u gr d/ dk = g/2 u ph/2 diverging at 0.32
This divergence is an artifact of our assumption of the infinitely deep liquid. A rather straightforward generalization of the above calculations to a layer of a finite thickness h, using the additional boundary condition vy y=- h = 0, yields a more general dispersion relation:33
2
gk tanh kh.
(8.48)
It shows that relatively long waves, with >> h, i.e. with kh << 1, propagate without dispersion (i.e.
have / k = const u), with the following velocity:
u gh1/ 2 .
(8.49)
For the Earth’s oceans, this velocity is rather high, close to 250 m/s (!) for the average ocean depth h 5
km. This result explains, in particular, the very fast propagation of tsunami waves.
In the opposite limit of very short waves (large k), Eq. (47) also does not give a good description of typical experimental data, due to surface tension effects – see Sec. 2 above. Using Eq. (13), it is easy (and hence also left for the reader’s exercise) to show that their account leads (at kh >> 1) to the following modification of Eq. (47):
2
3
k
gk
.
(8.50)
32 Here, unlike in Chapters 6 and 7, the wave velocity is denoted by the letter u to avoid any chance of its confusion with the velocity v (43) of the liquid’s particles.
33 This calculation (left for the reader’s exercise), shows also that at finite h, the particle trajectories are elliptical rather than circular, becoming more and more stretched in the wave propagation direction near the bottom of the layer.
Chapter 8
Page 15 of 30
CM: Classical Mechanics
According to this formula, the surface tension is important at wavelengths smaller than the capillary constant a c given by Eq. (14). Much shorter waves, for that Eq. (50) yields k 3/2, are called the capillary waves – or just “ripples”.
8.5. Dynamics: Viscous fluids
The viscosity of many fluids, at not overly high velocities, may be described surprisingly well by adding, to the static stress tensor (2), additional elements proportional to the velocity v dq/ dt:
~
P
(v) .
(8.51)
jj '
jj'
jj'
In view of our experience with the Hooke’s law (7.32) expressing a stress tensor proportional to particle displacements q, we may expect a similar expression with the replacement q v = dq/ dt:
1
1
~
2 e
(8.52a)
jj'
Tr
jj'
jj
jj'
'
e 3
Tre ,
3
3
where ejj’ are the elements of the symmetrized strain derivative tensor:
ds
1 v
v
e
jj' j
j' .
(8.52b)
jj'
dt
2 r
r
j'
j
Experiment confirms that Eq. (52) gives a good description of the viscosity effects in a broad range of isotropic fluids. The coefficient is called either the shear viscosity, or the dynamic viscosity, or just viscosity, while is called the second (or bulk) viscosity.
In the most frequent case of virtually incompressible fluids, Tr (e) = d[Tr (s)]/ dt = ( dV/ dt)/ V = 0, so that the term proportional to vanishes, and is the only important viscosity parameter.34 Table 1
shows the approximate values of the viscosity, together with the mass density , for several representative fluids.
Table 8.1. Important parameters of several representative fluids (approximate values)
Fluid (all at 300 K, until indicated otherwise)
(mPas)
(kg/m3)
Glasses
1021–1024
2,200–2,500
Earth magmas (at 800 to 1,400 K)
104–1014
2,200–2,800
Machine oils (SAE 10W – 40 W)
65-320
900
Water
0.89
1,000
Mercury
1.53
13,530
Liquid helium 4 (at 4.2K, 105 Pa)
0.019
130
Air (at 105 Pa)
0.018
1.3
34 Probably the most important effect we miss by neglecting is the attenuation of the (longitudinal) acoustic waves, into which the second viscosity makes a major contribution – whose (rather straightforward) analysis is left for the reader’s exercise.
Chapter 8
Page 16 of 30
Essential Graduate Physics
CM: Classical Mechanics
One can see that may vary in very broad limits; the extreme cases of fluids are glasses (which, somewhat counter-intuitively, are not stable solids even at room temperature, but rather may “flow”, though extremely slowly, until they eventually crystallize) and liquid helium (whose viscosity is of the order of that of gases,35 despite its much higher density).
Incorporating the additional elements of jj’ to the equation (23) of fluid motion, absolutely similarly to how it was done at the derivation of Eq. (7.107) of the elasticity theory, and with the account of Eq. (19), we arrive at the famous Navier-Stokes equation:36
Navier-
v
2
Stokes
v v P f v v.
(8.53)
equation
t
3
The apparent simplicity of this equation should not mask an enormous range of phenomena, notably including turbulence (see the next section) that are described by it, and the complexity of its solutions even for some simple geometries. In most problems interesting for practice, the only option is to use numerical methods, but due to a large number of parameters (, , , plus geometrical parameters of the involved bodies, plus the distribution of bulk forces f, plus boundary conditions), this way is strongly plagued by the curse of dimensionality that was discussed in the end of Sec. 5.8.
Let us see how the Navier-Stokes equation works, on several simple examples. As the simplest case, let us consider the so-called Couette flow of an incompressible fluid layer between two wide, horizontal plates (Fig. 10), caused by their mutual sliding with a constant relative velocity v0.
v
y
0
d
v( y)
Fig. 8.10. The simplest problem of
the viscous fluid flow.
0
z
Let us assume a laminar (vorticity-free) fluid flow. (As will be discussed in the next section, this assumption is only valid within certain limits.) Then we may use the evident symmetry of the problem, to take, in the coordinate frame shown in Fig. 10, v = nz v( y). Let the bulk forces be vertical, f = n yf, so they do not give an additional drive to the fluid flow. Then for the stationary flow (v/ t = 0), the vertical, y- component of the Navier-Stokes equation is reduced to the static Pascal equation (6), showing that the pressure distribution is not affected by the plate (and fluid) motion. In the horizontal, z-
component of the equation, only one term, 2 v, survives, so that for the only Cartesian component of the fluid’s velocity we get the 1D Laplace equation
2
d v
.
0
(8.54)
2
dy
35 Actually, at even lower temperatures (for He 4, at T < T 2.17 K), helium becomes a superfluid, i.e. loses its viscosity completely, as a result of the Bose-Einstein condensation – see, e.g., SM Sec. 3.4.
36 Named after Claude-Louis Navier (1785-1836) who had suggested the equation, and Sir George Gabriel Stokes (1819-1903) who has demonstrated its relevance by solving the equation for several key situations.
Chapter 8
Page 17 of 30
CM: Classical Mechanics
In contrast to the ideal fluid (see, e.g., Fig. 8b), the relative velocity of a viscous fluid and a solid wall it flows by should approach zero at the wall,37 so that Eq. (54) should be solved with boundary conditions
,
0
at y ,
0
v
(8.55)
v , at y d.
0
Using the evident solution to this boundary problem, v( y) = ( y/ d) v 0, illustrated by the arrows in Fig. 10, we can now calculate the horizontal drag force acting on a unit area of each plate. For the bottom plate, F
v
v
z
0
.
(8.56)
zy y0
y0
A
y
d
y
(For the top plate, the derivative v/ y has the same value, but the sign of dAy has to be changed to reflect the direction of the outer normal to the solid surface so that we get a similar force but with the negative sign.) The well-known result (56) is often used, in undergraduate physics courses, for a definition of the dynamic viscosity , and indeed shows its meaning very well.38
As the next, slightly less trivial example let us consider the so-called Poiseuille problem:39
finding the relation between the constant external pressure gradient - P/ z applied along a round pipe with internal radius R (Fig. 11) and the so-called discharge Q – defined as the mass of fluid flowing through the pipe’s cross-section in unit time.
higher
R
lower
pressure
v()
pressure
0
z
Fig. 8.11. The Poiseuille problem.
Again assuming a laminar flow, we can involve the problem’s uniformity along the z-axis and its axial symmetry to infer that v = n zv(), and P = - z + f(, ) + const (where = {, } is again the 2D
radius vector rather than the fluid density) so that the Navier-Stokes equation (53) for an incompressible fluid (with v = 0) is reduced to the following 2D Poisson equation:
2 v .
(8.57)
2
After spelling out the 2D Laplace operator in polar coordinates for our axially-symmetric case / = 0, Eq. (57) becomes a simple ordinary differential equation,
37 This is essentially an additional experimental fact, but may be understood as follows. The tangential component of the velocity should be continuous at the interface between two viscous fluids, in order to avoid infinite stress –
see Eq. (52), and solid may be considered as an ultimate case of fluid, with infinite viscosity.
38 The very notion of viscosity was introduced (by nobody other than the same Sir Isaac Newton) via a formula similar to Eq. (56), so that any effect resulting in a drag force proportional to velocity is frequently called the Newtonian viscosity.
39 It was solved by G. Stokes in 1845 to explain the experimental results obtained by Gotthilf Hagen in 1839 and (independently) by Jean Poiseuille in 1840-41.
Chapter 8
Page 18 of 30
Essential Graduate Physics
CM: Classical Mechanics
1 d
dv
(
) ,
(8.58)
d
d
which has to be solved on the segment 0 R, with the following boundary conditions: v ,
0
at R,
dv
(8.59)
,
0
at 0.
d
(The latter condition is required by the axial symmetry.) A straightforward double integration yields:
v
2 2
R ,
(8.60)
4
so that the (easy) integration of the mass flow density over the cross-section of the pipe, R
2
Q vd r 2
( 2
2
R ) d ,
(8.61)
4
A
0
immediately gives us the so-called Poiseuille (or “Hagen-Poiseuille”) law for the fluid discharge: Poiseuille
law
4
Q
R .
(8.62)
8
The most prominent (and practically important) feature of this result is the very strong dependence of the discharge on the pipe’s radius.
Of course, the 2D Poisson equation (57) is so readily solvable not for each cross-section shape.
For example, consider a very simple, square-shaped cross-section with side a (Fig. 12).
y
a / 2
Fig. 8.12. Application of the finite-difference
a
method with a very coarse mesh (with step h
= a/2) to the problem of viscous fluid flow in
2
a pipe with a square cross-section.
0
x
In this case, it is natural to use the Cartesian coordinates aligned with the cross-section’s sides, so that Eq. (57) becomes
2
2
v v
0
for
const,
x, y a,
(8.63)
2
2
x
y
and has to be solved with boundary conditions
v ,
0
at x, y ,
0 a .
(8.64)
Chapter 8
Page 19 of 30
CM: Classical Mechanics
For this boundary problem, analytical methods such as the variable separation lead to answers in the form of infinite sums (series), which ultimately require computers anyway – at least for their plotting and comprehension. Let me use this pretext to discuss how explicitly numerical methods may be used for such problems – or for any partial differential equations involving the Laplace operator. The simplest of them is the finite-difference method40 in which the function to be calculated, f(r), is represented by its values f(r1), f(r2), … in discrete points of a rectangular grid (frequently called mesh) of the corresponding dimensionality – Fig. 13.
(a)
(b)
f
f
f
( f / r )
j
f
/ r
j
f
f
f
f
h
Fig. 8.13. The idea of the finite-
f
difference method in (a) one and
h
(b) two dimensions.
r
h
h
j
In Sec. 5.7, we have already discussed how to use such a grid to approximate the first derivative of the function – see Eq. (5.97). Its extension to the second derivative is straightforward – see Fig. 13a:41
2
f
f
1 f
f
1 f f
f f
f f 2 f
.
(8.65)
2
2
r
r
r
h
r
r
h
h
h
h
j
j
j
j
j
The relative error of this approximation is of the order of h 24 f/ r 4 j, quite acceptable in many cases. As a result, the left-hand side of Eq. (63), treated on a square mesh with step h (Fig. 13b), may be approximated with the so-called five-point scheme:
2
2
v v
v v 2 v
v v 2 v
v v v v 4 v
.
(8.66)
2
2
2
2
2
x
y
h
h
h
(The generalization to the seven- point scheme, appropriate for 3D problems, is straightforward.) Let us apply this scheme to the pipe with the square cross-section, using an extremely coarse mesh with step h
= a/2, shown in Fig. 12. In this case, the fluid velocity v should equal zero at the walls, i.e. at all points of the five-point scheme except for the central point (in which the velocity obviously reaches its maximum), so that Eqs. (63) and (66) yield42
40 For more details see, e.g., R. Leveque, Finite Difference Methods for Ordinary and Partial Differential Equations, SIAM, 2007.
41 As a reminder, at the beginning of Sec. 6.4 we have already discussed the reciprocal transition – from a similar sum to the second derivative in the continuous limit ( h 0).
42 Note that value (67) of v max is exactly the same as given by the analytical formula (60) for the round cross-section with the radius R = a/2. This is not an occasional coincidence. The velocity distribution given by (60) is a quadratic function of both x and y. For such functions, with all derivatives higher than 2 f/ r 2 j equal to zero, equation (66) is exact rather than approximate.
Chapter 8
Page 20 of 30
Essential Graduate Physics
CM: Classical Mechanics
0 0 0 0 4 v
1 a
max
2
,
i.e. v
(8.67)
2
( a / )
2
max
16
This result for the maximal velocity is only ~20% different from the exact value. Using a slightly finer mesh with h = a/4, which gives a readily solvable system of three linear equations for three different velocity values (the exercise left for the reader), brings us within just a couple of percent from the exact result. So numerical methods may be practically more efficient than the “analytical” ones, even if the only available tool is a calculator app on your smartphone rather than an advanced computer.
Of course, many practical problems of fluid dynamics do require high-performance computing, especially in conditions of turbulence with its complex, irregular spatial-temporal structure – see the next section). In such cases, the finite-difference approach discussed above may become unsatisfactory, because it implies the same accuracy of the derivative approximation through the whole area of interest.
A more powerful (but also much more complex for implementation) approach is the finite-element method in which the discrete-point mesh is based on triangles with unequal sides, and is (in most cases, automatically) generated from the system geometry, giving many more mesh points at the location(s) of the highest gradients of the calculated function (Fig. 14), and hence a better calculation accuracy for the same total number of points. Unfortunately, I do not have time for going into the details of that method, so the interested reader is referred to the special literature on this subject.43
Fig. 8.14. A typical finite-element
mesh generated automatically for a
system with relatively complex
geometry – a round cylindrical shell
inside another one, with mutually
perpendicular axes. (Adapted from
the
original
by
I.
Zureks,
https://commons.wikimedia.org/w/in
dex.php?curid=2358783, under the
CC license BY-SA 3.0.)
Before proceeding to our next topic, let me mention one more important problem that is analytically solvable using the Navier-Stokes equation: a slow motion of a solid sphere of radius R, with a constant velocity v0, through an incompressible viscous fluid – or equivalently, a slow flow of the fluid (uniform at large distances) around an immobile sphere. In the limit v 0, the second term on the left-hand side of Eq. (53) is negligible (just as at the surface wave analysis in Sec. 3), the equation takes the form
2
P v ,
0
for R r ,
(8.68)
43 I can recommend, e.g., C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, Dover, 2009, or T. Hughes, The Finite Element Method, Dover, 2000.
Chapter 8
Page 21 of 30
CM: Classical Mechanics
and should be complemented with the incompressibility condition v = 0 and the boundary conditions v 0, at r R,
(8.69)
v v , at r .
0
In spherical coordinates, with the polar axis directed along the vector v0, this boundary problem has the axial symmetry (so that v/ = 0 and v = 0), and allows the following analytical solution:
3
3
R
R
3
3
R
R
3 v R
v v cos 1
,
0
v v
P
.
(8.70)
r
0
2
sin 1
,
cos
2 r
2
0
r
4 r
4 2
r
2 2
r
Now calculating the tensor elements (52b) at r = R, using them to find the stress tensor elements from Eq. (52a), and integrating the elementary forces (7.18) over the surface of the sphere, it is straightforward to obtain the famous Stokes formula for the drag force acting on the sphere:44
F 6 Rv .
(8.71) Stokes
0
formula
For water drops with a 1-micron diameter, usually taken for the border between aerosols and droplets, descending in the ambient-condition air under their own weight, it predicts an equilibrium velocity v of close to 0.1 meter per hour, with the further scaling v R 2.45 (Note, however, that at R below ~10 m, corrections due to air molecule discreteness become noticeable.)
For what follows in the next section, it is convenient to recast this result into the following form: 24
C
,
(8.72)
d
e
R
where C d is the drag coefficient defined as
F
C
,
(8.73)
d
2
v
A/ 2
0
with A R 2 being the sphere’s cross-section “as seen by the incident fluid flow”, and Re is the so-called Reynolds number.46 In the general case, the number is defined as
vl
e
R
,
(8.74) Reynolds
number
where l is the linear-size scale of the problem, and v is its velocity scale. (In the particular case of Eq.
(72) for the sphere, l is identified with the sphere’s diameter D = 2 R, and v with v 0). The physical sense of these two definitions will be discussed in the next section.
44 This formula has played an important role in the first precise (better than 1%) calculation of the fundamental electric charge e by R. Millikan and H. Fletcher from their famous oil drop experiments in 1909-1913.
45 These numbers are of key importance not only for the contagious disease transmission analysis, but also for many other fields including atmospheric physics. For example, for an average water droplet in clouds, with R ~ 10
m, Eq. (71) (even with a due account of the slightly lower air viscosity at typical cloud heights) yields the equilibrium descent velocity of the order of 10 m/hr, substantiating the correct answer to the popular high-school question, “Why clouds do not fall?” (The answer is: each water droplet does descend, but so slowly that it has ample time to evaporate at the lower surface of the cloud, so that the cloud as a whole may maintain its height.) 46 This notion was introduced in 1851 by the same G. Stokes but eventually named after O. Reynolds who popularized it three decades later.
Chapter 8
Page 22 of 30
Essential Graduate Physics
CM: Classical Mechanics
8.6. Turbulence
As Fig. 15 shows, the Stokes result (71)-(72) is only valid at Re << 1, while for larger values of the Reynolds number, i.e. at higher velocities v 0, the drag force is larger. This very fact is not quite surprising, because at the derivation of the Stokes’ result, the nonlinear term (v)v in the Navier-Stokes equation (53), which scales as v 2, was neglected in comparison with the linear terms, scaling as v. What is more surprising is that the function C d( Re) exhibits such a complicated behavior over many orders of the velocity’s magnitude, giving a hint that the fluid flow at large Reynolds numbers should be also very complicated. Indeed, the reason for this complexity is a gradual development of very intricate, time-dependent fluid patterns, called turbulence, rich with vortices – for example, see Fig. 16. These vortices are especially pronounced in the region behind the moving body (the so-called wake), while the region before the body remains almost unperturbed. As Fig. 15 indicates, the turbulence exhibits rather different behaviors at various velocities (i.e. values of Re), and sometimes changes rather abruptly – see, for example, the significant drag’s drop at Re 5105.
C d
Re
Fig. 8.15. The drag coefficient for a sphere and a thin round disk as functions of the Reynolds number.
Adapted from F. Eisner, Das Widerstandsproblem, Proc. 3rd Int. Cong. on Appl. Mech., Stockholm, 1931.
In order to understand the conditions of this phenomenon, let us estimate the scale of various terms of the Navier-Stokes equation (53) for a generic body with characteristic size l, moving in an otherwise static incompressible fluid, with velocity v. In this case, the time scale of possible non-stationary phenomena is given by the ratio l/ v,47 so that we arrive at the following estimates: 47 The time scale of phenomena in externally-driven systems may be different; for example, for forced oscillations with frequency , it may be the oscillation period T 2/. For such problems, the ratio S ( l/ v)/ T, commonly called either the Strouhal number or the reduced frequency, serves as another dimensionless constant.
Chapter 8
Page 23 of 30
Essential Graduate Physics
CM: Classical Mechanics
v
v
2
v
f
v
Equation term:
t
(8.75)
2
2
v
v
v
g
2
Order of magnitude:
l
l
l
(I have skipped the term P because as we have seen in the previous section, in typical fluid flow problems, it balances the viscosity term, and hence is of the same order of magnitude.)
uniform
fluid flow
solid
sphere
Fig. 8.16. A snapshot of the turbulent tail ( wake) behind a sphere moving in a fluid with a high Reynolds number, showing the so-called von Kármán vortex street. Adapted from the original (actually, a very nice animation, http://www.mcef.ep.usp.br/staff/jmeneg/cesareo/vort2.gif) by Cesareo de La Rosa Siqueira, as a copyright-free material, available at https://commons.wikimedia.org/w/index.php?curid=87351.
Estimates (75) show that the relative importance of the terms may be characterized by two dimensionless ratios.48 The first of them is the so-called Froude number 49
2
v
/
2
l
v
F
,
(8.76)
g
lg
which characterizes the relative importance of the gravity – or, upon appropriate modification, of other bulk forces. In most practical problems (with the important exception of surface waves, see Sec. 4
above), F >> 1 so that the gravity effects may be neglected.
Much more important is another ratio, the Reynolds number (74), which may be rewritten as 2
vl
v
/ l
e
R
,
(8.77)
2
v / l
and hence is a measure of the relative importance of the fluid particle’s inertia in comparison with the viscosity effects.50 So again, it is natural that for a sphere, the role of the vorticity-creating term (v)v 48 For substantially compressible fluids (e.g., gases), the most important additional dimensionless parameter is the Mach number M v/ v l, where v l = ( K/)1/2 is the velocity of the longitudinal sound – which is, as we already know from Chapter 7, the only wave mode possible in an infinite fluid. Especially significant for practice are supersonic effects (including the shock wave in the form of the famous Mach cone with half-angle M = sin-1 M-1) that arise at M > 1. For a more thorough discussion of these issues, I have to refer the reader to more specialized texts – either Chapter IX of the Landau-Lifshitz volume cited above or Chapter 15 in I. Cohen and P. Kundu,
Fluid Mechanics, 4th ed., Academic Press, 2007 – which is generally a good book on the subject.
49 Named after William Froude (1810-1879), one of the applied hydrodynamics pioneers.
50 Note that the “dynamic” viscosity participates in this number (and many other problems of fluid dynamics) only in the combination /, which thereby has deserved a special name of kinematic viscosity.
Chapter 8
Page 24 of 30
CM: Classical Mechanics
becomes noticeable already at Re ~ 1 – see Fig. 15. What is very counter-intuitive is the onset of turbulence in systems where the laminar (turbulence-free) flow is formally an exact solution to the Navier-Stokes equation for any Re. For example, at Re > Re t 2,100 (with l 2 R and v v max) the laminar flow in a round pipe, described by Eq. (60), becomes unstable, and the resulting turbulence decreases the fluid discharge Q in comparison with the Poiseuille law (62). Even more strikingly, the critical value of Re is rather insensitive to the pipe wall roughness and does not diverge even in the limit of perfectly smooth walls.
Since Re >> 1 in many real-life situations, turbulence is very important for practice. (Indeed, the values of and for water listed in Table 1 imply that even for a few-meter-sized object, such as a human body or a small boat, Re > 1,000 at any speed above just ~1 mm/s.) However, despite nearly a century of intensive research, there is no general, quantitative analytical theory of this phenomenon, and most results are still obtained either by rather approximate analytical treatments, or by the numerical solution of the Navier-Stokes equations using the approaches discussed in the previous section, or in experiments (e.g., on scaled models51 in wind tunnels). A rare exception is the relatively recent theoretical result by S. Orszag (1971) for the turbulence threshold in a flow of an incompressible fluid through a gap of thickness t between two parallel plane walls (see Fig. 10): Re t 5,772 (for l t/2 and v
v max). However, even for this simplest geometry, the analytical theory still cannot predict the turbulence patterns at Re > Re t. Only certain general, semi-quantitative features of turbulence may be understood from simple arguments.
For example, Fig. 15 shows that within a very broad range of Reynolds numbers, from ~102 to
~3105, C d of a thin round disk perpendicular to the incident flow, C d is very close to 1.1 for any Re > 103, and that of a sphere is not too far away. The approximate equality C d 1, meaning that the drag force F is close to v 2
0 A/2, may be understood (in the picture where the object is moved by an external force F with the velocity v 0 through a fluid that was initially at rest) as the equality of the force-delivered power Fv
2
0 and the fluid’s kinetic energy ( v 0 /2) V created in volume V = v 0 A in unit time. This relation would be exact if the object gave its velocity v 0 to each and every fluid particle its cross-section runs into, for example by dragging all such particles behind itself. In reality, much of this kinetic energy goes into vortices, where the particle velocity may differ from v 0, so that the equality C d 1 is only approximate.
Another important general effect is that at very high values of Re, fluid flow at the leading surface of solid objects forms a thin, highly turbulent boundary layer that matches the zero relative velocity of the fluid at the surface with its substantial velocity in the outer region, which is almost free of turbulence and many cases, of other viscosity effects. This fact, clearly visible in Fig. 16, enables semi-quantitative analyses of several effects, for example, the so-called Magnus lift force 52 Fl exerted (on top of the usual drag force Fd) on rotating objects, and directed across the fluid flow – see Fig. 17.
An even more important application of this concept is an approximate analysis of the forces exerted on non-rotating airfoils (such as aircraft wings) with special cross-sections forming sharp angles at their back ends. Such a shape minimizes the airfoil’s contacts with the vortex street it creates in its 51 The crucial condition of correct modeling is the equality of the Reynolds numbers (74) (and if relevant, also of the Froude numbers and/or the Mach numbers) of the object of interest and its model.
52 Named after G. Magnus, who studied this effect in detail in 1852, though it had been described much earlier (in 1672) by I. Newton, and by B. Robins after him (in 1742).
Chapter 8
Page 25 of 30
CM: Classical Mechanics
wake, and allows the thin boundary layer to extend over virtually all of its surface, enhancing the lift force.
v0
Fl
Fd
Fig. 8.17. The Magnus effect.
Unfortunately, due to the time/space restrictions, for a more detailed discussion of these results I have to refer the reader to more specialized literature,53 and will conclude this chapter with a brief discussion of just one issue: can turbulence be explained by a single mechanism? (In other words, can it be reduced, at least on a semi-quantitative level, to a set of simpler phenomena that are commonly considered “well understood”?) Apparently the answer is no,54 though nonlinear dynamics of simpler systems may provide some useful insights.
In the middle of the last century, the most popular qualitative explanation of turbulence had been the formation of an “energy cascade” that would transfer the energy from the regular fluid flow to a hierarchy of vortices of various sizes.55 With our background, it is easier to retell that story in the time-domain language (with the velocity v serving as the conversion factor), using the fact that in a rotating vortex, each Cartesian component of a particle’s radius vector oscillates in time, so that to some extent the vortex plays the role of an oscillatory motion mode.
Let us consider the passage of a solid body between two, initially close, small parts of the fluid.
The body pushes them apart, but after its passage, these partial volumes are free to return to their initial positions. However, the dominance of inertia effects at motion with Re >> 1 means that the volumes continue to oscillate for a while about those equilibrium positions. (Since elementary volumes of an incompressible fluid cannot merge, these oscillations take the form of rotating vortices – see Fig. 16
again.)
Now, from Sec. 5.8 we know that intensive oscillations in a system with the quadratic nonlinearity, in this case, provided by the convective term (v)v, are equivalent, for small perturbations, to the oscillation of the system’s parameters at the corresponding frequency. On the other hand, as was briefly discussed in Sec. 6.7, in a system with two oscillatory degrees of freedom, a periodic parameter change with frequency p may lead to the non-degenerate parametric excitation (“down-conversion”) of oscillations with frequencies 1,2 satisfying the relation 1 + 2 = p.
Moreover, the spectrum of oscillations in such a system also has higher combinational frequencies such as (p + 1), thus pushing the oscillation energy up the frequency scale (“up-conversion”). In the presence of other oscillatory modes, these oscillations may in turn produce, via the same nonlinearity, even higher frequencies, etc. In a fluid, the spectrum of these “oscillatory modes” (actually, vortex 53 See, e.g., P. Davidson, Turbulence, Oxford U. Press, 2004.
54 The following famous quote is attributed to Werner Heisenberg on his deathbed: “When I meet God, I will ask him two questions: Why relativity? And why turbulence? I think he will have an answer for the first question.”
Though probably inaccurate, this story reflects rather well the frustration of the fundamental physics community, renown for their reductionist mentality, with the enormous complexity of phenomena that obey simple (e.g., the Navier-Stokes) equations, i.e. from their point of view, do not describe any new physics.
55 This picture was suggested in 1922 by Lewis F. Richardson.
Chapter 8
Page 26 of 30
Essential Graduate Physics
CM: Classical Mechanics
structures) is essentially continuous, so that the above arguments make very plausible a sequential transfer of the energy from the moving body to a broad range of oscillatory modes – whose frequency spectrum is limited from above by the energy dissipation due to the fluid’s viscosity. When excited, these modes interact (in particular, mutually phase-lock) via the system’s nonlinearity, creating the complex motion we call turbulence.
Though not having much quantitative predictive power, such handwaving explanations, which are essentially based on the excitation of a large number of effective degrees of freedom, had been dominating the turbulence reviews until the mid-1960s. At that point, the discovery (or rather re-discovery) of quasi-random motions in classical dynamic systems with just a few degrees of freedom altered the discussion substantially. Since this phenomenon, called the deterministic chaos, extends well beyond the fluid dynamics, I will devote to it a separate (albeit short) next chapter, and in its end will briefly return to the discussion of turbulence.
8.7. Exercise problems
8.1. For a mirror-symmetric but otherwise arbitrary shape of a ship’s hull, derive an explicit expression for the height of metacenter M – see Fig. 3. Spell out this expression for a rectangular hull.
8.2. Find the stationary shape of the open surface of an
incompressible, heavy fluid in a container rotated about its vertical axis with
a constant angular velocity – see the figure on the right.
g
8.3. In the first order in the so-called flattening f ( R e – R p)/ R p << 1 of the Earth (where R e and R p are, respectively, its equatorial and polar radii), calculate it within a simple model in that our planet is a uniformly-rotating nearly-spherical fluid ball, whose gravity field is dominated by a relatively small spherical core. Compare your result with the experimental value of f, and interpret the difference.
Hint: You may use experimental values R e 6,378 km, R p 6,357 km, and g 9.807 m/s2.
8.4.* Use two different approaches to calculate the stationary shape of the
surface of an incompressible fluid of density near a vertical plane wall, in a uniform gravity field – see the figure on the right. In particular, find the height h g
h
of liquid’s rise at the wall surface as a function of the contact angle
?
c.
2 R
8.5.* A soap film with surface tension is stretched between two
similar, coaxial, thin, round rings of radius R, separated by distance d – see the figure on the right. Neglecting gravity, calculate the equilibrium shape of the d film, and the external force needed for keeping the rings at this distance.
Chapter 8
Page 27 of 30
Essential Graduate Physics
CM: Classical Mechanics
8.6. A solid sphere of radius R is kept in a steady, vorticity-free flow of an ideal incompressible fluid, with velocity v 0. Find the spatial distribution of velocity and pressure, and in particular their extreme values. Compare the results with those obtained in Sec. 4 for a round cylinder.
8.7.* Solve the same problem for a long and thin solid strip of width 2 w, with its plane normal to the unperturbed fluid flow.
Hint: You may like to use the so-called elliptic coordinates {, } defined by their relations with the Cartesian coordinates { x, y}:
x C cosh
cos , y C sinh sin ,
0
with ,
,
where C is a constant; in these coordinates,
1
2
2
2
.
2
2
2
2
2
C (cosh cos )
8.8. A small source, located at distance d from a plain wall of a container
filled with an ideal, incompressible fluid of density , injects additional fluid isotropically, at a constant mass current (“discharge”) Q dM/ dt – see the figure d
on the right. Calculate the fluid’s velocity distribution, and its pressure on the wall, created by the flow.
Q
Hint: Recall the charge image method in electrostatics,56 and contemplate
its possible analog.
8.9. Calculate the average kinetic, potential, and full energies (per unit area) of a traveling sinusoidal wave, of a small amplitude qA, on the horizontal surface of an ideal, incompressible, deep fluid of density , in a uniform gravity field g.
8.10. Calculate the average power (per unit width of the wave’s front) carried by the surface wave discussed in the previous problem, and relate the result to the wave’s energy.
8.11. Derive Eq. (48) for the surface waves on a finite-thickness layer of a heavy liquid.
8.12. The utmost simplicity of Eq. (49) for the velocity of waves on a relatively shallow ( h << ) layer of an ideal incompressible liquid implies that they may be described using a very simple physical picture. Develop such a picture, and verify that it yields the same expression for the velocity.
8.13. Use the solution of the previous problem to calculate the energy and power of the shallow-layer waves, and use the result to explain the high tides on some ocean shores, using two models: (i) the water depth h decreases gradually toward the shore, and
(ii) h decreases sharply, at some distance l from the shore – as it does on the ocean shelf border.
8.14.* Derive a 2D differential equation describing the propagation of relatively long ( >> h) waves on the surface of a broad, plane layer of thickness h, of an ideal, incompressible fluid, and use it 56 See, e.g., EM Secs. 2.9, 3.3, and 4.3.
Chapter 8
Page 28 of 30
Essential Graduate Physics
CM: Classical Mechanics
to calculate the longest standing wave modes and frequencies in a layer covering a spherical planet of radius R >> h.
Hint: The second task requires some familiarity with the basic properties of spherical harmonics.57
8.15. Calculate the velocity distribution and the dispersion relation of the waves propagating along the horizontal interface of two ideal, incompressible fluids of different densities.
8.16. Derive Eq. (50) for the capillary waves (“ripples”).
8.17. Use the finite-difference approximation for the Laplace operator, with the mesh step h =
a/4, to find the maximum velocity and total mass flow Q of a viscous, incompressible fluid through a long pipe with a square-shaped cross-section of side a. Compare the results with those described in Sec.
5 for the same problem with the mesh step h = a/2, and for a pipe with the circular cross-section of the same area.
x
h
8.18. A layer, of thickness h, of a heavy, viscous,
incompressible fluid flows down a long and wide inclined plane, 0
under its own weight – see the figure on the right. Find the
v( x) ?
stationary velocity distribution profile, and the total fluid discharge
,
(per unit width.)
g
8.19. An external force moves two coaxial round disks of radius R, with an incompressible viscous fluid in the gap between them, toward each other with a constant velocity u. Calculate the applied force in the limit when the gap’s thickness t is already much smaller than R.
8.20. Calculate the drag torque exerted on a unit length of a solid round cylinder of radius R that rotates about its axis, with angular velocity , inside an incompressible fluid with viscosity , kept static far from the cylinder.
8.21. Solve a similar problem for a sphere of radius R, rotating about one of its principal axes.
8.22. Calculate the tangential force (per unit area) exerted by an incompressible fluid, with density and viscosity , on a broad solid plane placed over its surface and forced to oscillate, along the surface, with amplitude a and frequency .
8.23. A massive barge, with a flat bottom of area
V t
A, floats in shallow water, with clearance h << A 1/2 – y
?
see the figure on the right. Analyze the time dependence
of the barge’s velocity V( t), and the water velocity
M
,
profile, after the barge’s engine has been turned off.
x
h
0
Discuss the limits of large and small values of the
dimensionless parameter M/ Ah.
57 See, e.g., EM Sec. 2.8 and/or QM Sec. 3.6.
Chapter 8
Page 29 of 30
CM: Classical Mechanics
8.24.* Derive a general expression for mechanical energy loss rate in a viscous incompressible fluid that obeys the Navier-Stokes equation, and use this expression to calculate the attenuation coefficient of the surface waves, assuming that the viscosity is small. (Quantify this condition).
8.25. Use the Navier-Stokes equation to calculate the coefficient of attenuation of a plane, sinusoidal acoustic wave.
8.26.* Use two different approaches for a semi-quantitative calculation of the Magnus lift force Fl exerted by an incompressible fluid of density on a round cylinder of radius R, with its axis normal to the fluid’s velocity v0, and rotating about the axis with an angular velocity – see Fig. 17. Discuss the relation of the results.
Chapter 8
Page 30 of 30
CM: Classical Mechanics
Chapter 9. Deterministic Chaos
This chapter gives a very brief review of chaotic phenomena in deterministic maps and dynamic systems with and without dissipation, and an even shorter discussion of the possible role of chaos in fluid turbulence.
9.1. Chaos in maps
The possibility of quasi-random dynamics of deterministic systems with a few degrees of freedom (nowadays called the deterministic chaos – or just “chaos”) had been noticed before the 20th century,1 but became broadly recognized only after the publication of a 1963 paper by theoretical meteorologist Edward Lorenz. In that work, he examined numerical solutions of the following system of three nonlinear, ordinary differential equations,
q a ( q q ),
1
1
2
1
Lorenz
system
q a q q q q ,
(9.1)
2
2 1
2
1 3
q q q a q ,
3
1 2
3 3
as a rudimentary model of heat transfer through a horizontal layer of fluid separating two solid plates.
(Experiment shows that if the bottom plate is kept hotter than the top one, the fluid may exhibit turbulent convection.) He has found that within a certain range of the constants a 1,2,3, the solution to Eq.
(1) follows complex, unpredictable, non-repeating trajectories in the 3D q-space. Moreover, the functions qj( t) (where j = 1, 2,3) are so sensitive to initial conditions qj(0) that at sufficiently large times t, solutions corresponding to slightly different initial conditions become completely different.
Very soon it was realized that such behavior is typical for even simpler mathematical objects called maps, so that I will start my discussion of chaos from these objects. A 1D map is essentially a rule for finding the next number qn+1 of a discrete sequence numbered by the integer index n, in the simplest cases using only its last known value qn. The most famous example is the so-called logistic map:2
Logistic
q
f ( q ) rq 1
( q ).
(9.2)
map
n 1
n
n
n
The basic properties of this map may be understood using its (hopefully, self-explanatory) graphical representation shown in Fig. 1.3 One can readily see that at r < 1 (Fig. 1a) the logistic map sequence rapidly converges to the trivial fixed point q(0) = 0 because each next value of q is less than the previous one. However, if r is increased above 1 (as in the example shown in Fig. 1b), the fixed point 1 It may be traced back at least to an 1892 paper by the same Jules Henri Poincaré who was already reverently mentioned in Chapter 5. Citing it: “…it may happen that small differences in the initial conditions produce very great ones in the final phenomena. […] Prediction becomes impossible.”
2 Its chaotic properties were first discussed in 1976 by Robert May, though the map itself is one of the simple ecological models repeatedly discussed much earlier, and may be traced back at least to the 1838 work by Pierre François Verhulst.
3 Since the maximum value of the function f( q), achieved at q = ½, equals r/4, the mapping may be limited to segment x = [0, 1], if the parameter r is between 0 and 4. Since all interesting properties of the map, including chaos, may be found within these limits, I will discuss only this range of r.
© K. Likharev
Essential Graduate Physics
CM: Classical Mechanics
q(0) becomes unstable. Indeed, at qn << 1, the map yields qn+1 rqn, so that at r > 1, the values qn grow with each iteration. Instead of the unstable point q(0) = 0, in the range 1 < r < r 1, where r 1 3, the map has a stable fixed point q(1) that may be found by plugging this value into both parts of Eq. (2):
)
1
(
q
f )1(
q
)
1
(
rq 1
(
)
1
(
q ),
(9.3)
giving q(1) = 1 – 1/ r – see the leftmost branch of the plot shown in Fig. 2.
y
(a)
y
(b)
y x
y x
q
r 1
2
r 1
q
y f ( x) rx 1
( x)
1
q 2
y f ( x)
q 3
0
q
q
)
1
(
1
1
x
0
q
q
q
q
1
x
2
1
2
3
Fig. 9.1. Graphical analysis of the logistic map for: (a) r < 1 and (b) r > 1.
r 3
r
r
1
2
(2)
q
Fig. 9.2. The fixed points and
)
1
(
1
chaotic regions of the logistic
q
q
1
r
(2)
q
map. Adapted, under the CCO
1.0 Universal Public Domain
Dedication, from the original
by Jordan Pierce, available at
http://en.wikipedia.org/wiki/Lo
gistic_map. (A very nice live
simulation of the map is also
available on this website.)
r
However, at r > r 1 = 3, the fixed point q(1) also becomes unstable. To prove that, let us take q q )1
(
q~ , assume that the deviation q~ from the fixed point q(1) is small, and linearize the map (2) n
n
n
in q~ – just as we repeatedly did for differential equations earlier in this course. The result is n
df
q~
.
(9.4)
(1 q
~
r 1
(
2 q )1
(
q~
)
(2 r q~
)
n 1
)
n
n
n
q q
dq
Chapter 9
Page 2 of 14
CM: Classical Mechanics
It shows that at 0 < 2 – r < 1, i.e. at 1 < r < 2, the deviations q~ decrease monotonically. At –1 < 2 – r n
< 0, i.e. in the range 2 < r < 3, the deviations’ sign alternates, but their magnitude still decreases – as in a stable focus, see Sec. 5.6. However, at –1 < 2 – r, i.e. r > r 1 3, the deviations grow by magnitude, while still changing their sign, at each step. Since Eq. (2) has no other fixed points, this means that at n
, the values qn do not converge to one point; rather, within the range r 1 < r < r 2, they approach a limit cycle of alternation of two points, q (2)
(2)
+
and q- , which satisfy the following system of algebraic
equations:
(2)
q
f q
q
f q
.
(9.5)
(2)
(2)
,
(2)
These points are also plotted in Fig. 2, as functions of the parameter r. What has happened at the point r 1
= 3 is called the period-doubling bifurcation.
The story repeats at r = r 2 1 + 6 3.45, where the system goes from the 2-point limit cycle to a 4-point cycle, then at r = r 3 3.54, where the limit cycle begins to consist of 8 alternating points, etc.
Most remarkably, the period-doubling bifurcation points rn, at that the number of points in the limit cycle doubles from 2 n- 1 points to 2 n points, become closer and closer to each other. Numerical simulations show that at n , these points obey the following asymptotic behavior: Feigenbaum
C
bifurcation
r r
,
where r
...,
5699
.
3
...
6692
.
4
(9.6)
n
n
sequence
The parameter is called the Feigenbaum constant; for other maps, and some dynamic systems (see the next section), period-doubling sequences follow a similar law, but with different values of .
More important for us, however, is what happens at r > r. Numerous numerical experiments, repeated with increasing precision,4 have confirmed that here the system is disordered, with no reproducible limit cycle, though (as Fig. 2 shows) at r r, all sequential values qn are still confined to a few narrow regions.5 However, as parameter r is increased well beyond r, these regions broaden and merge. This is the so-called deep chaos, with no apparent order at all.6
The most important feature of the chaos (in this and any other system) is the exponential divergence of trajectories. For a 1D map, this means that even if the initial conditions q 1 in two map implementations differ by a very small amount q 1, the difference qn between the corresponding sequences qn is growing, on average, exponentially with n. Such exponents may be used to characterize chaos. Indeed, an evident generalization of the linearized Eq. (4) to an arbitrary point qn is df
Δ q
e Δ q ,
e
.
(9.7)
n 1
n
n
n
q qn
dq
4 The reader should remember that just like the usual (“nature”) experiments, numerical experiments also have limited accuracy, due to unavoidable rounding errors.
5 The geometry of these regions is essentially fractal, i.e. has a dimensionality intermediate between 0 (which any final set of geometric points would have) and 1 (pertinent to a 1D continuum). An extensive discussion of fractal geometries and their relation to the deterministic chaos may be found, e.g., in the book by B. Mandelbrot, The Fractal Geometry of Nature, W. H. Freeman, 1983.
6 This does not mean that chaos’ depth is always a monotonic function of r. As Fig. 2 shows, within certain intervals of this parameter, the chaotic behavior suddenly disappears, being replaced, typically, with a few-point limit cycle, just to resume on the other side of the interval. Sometimes (but not always!) the “route to chaos” on the borders of these intervals follows the same Feigenbaum sequence of period-doubling bifurcations.
Chapter 9
Page 3 of 14
Essential Graduate Physics
CM: Classical Mechanics
Let us assume that q 1 is so small that N first values qn are relatively close to each other. Then using Eq.
(7) iteratively for these steps, we get
N
N
q
Δ
q
Δ
Δ q
e ,
that
so
ln
ln e .
(9.8)
N
1
N
n
n
Δ
1
q
n
1
n1
Numerical experiments show that in most chaotic regimes, at N such a sum fluctuates about an average, which grows as N, with the parameter
N
1
lim
Lyapunov
q
lim
ln
,
(9.9)
Δ
0
N
en
exponent
1
N n1
called the Lyapunov exponent,7 being independent of the initial conditions. The bottom panel in Fig. 3
shows as a function of the parameter r for the logistic map (2). (Its top panel shows the same pattern as Fig. 2, which is reproduced here just for the sake of comparison.)
q
Fig. 9.3. The Lyapunov exponent
0
for the logistic map. Adapted, with
permission, from the monograph by
Schuster and Just (cited below). ©
Wiley-VCH Verlag GmbH & Co.
KGaA.
Note that at r < r, is negative, indicating the sequence’s stability, besides the points r 1, r 2, …
where would become positive if the limit cycle changes (bifurcations) had not brought it back into the negative territory. However, at r > r, becomes positive, returning to negative values only in limited intervals of stable limit cycles. It is evident that in numerical experiments (which dominate the studies of deterministic chaos) the Lyapunov exponent may be used as a good measure of the chaos’ depth.8
7 After Alexandr Mikhailovich Lyapunov (1857-1918), famous for his studies of the stability of dynamic systems.
8 N-dimensional maps that relate N-dimensional vectors rather than scalars, may be characterized by N Lyapunov exponents rather than one. For chaotic behavior, it is sufficient for just one of them to become positive. For such systems, another measure of chaos, the Kolmogorov entropy, may be more relevant. This measure and its relation with the Lyapunov exponents are discussed, for example, in SM Sec. 2.2.
Chapter 9
Page 4 of 14
CM: Classical Mechanics
Despite the abundance of results published for particular maps,9 and several interesting observations (like the already discussed existence of the Feigenbaum bifurcation sequences), to the best of my knowledge, nobody can yet predict the patterns like those shown in Fig. 2 and 3 from just studying the mapping rule itself, i.e. without carrying out actual numerical experiments. Unfortunately, the understanding of deterministic chaos in other systems is not much better.
9.2. Chaos in dynamic systems
Proceeding to the discussion of chaos in dynamic systems , it is more natural, with our background, to illustrate this discussion not with the Lorenz equations, but with the system of equations describing a dissipative pendulum driven by a sinusoidal external force, which was repeatedly discussed in Chapter 5. Introducing two new variables, the normalized momentum p ( dq/ dt)/0, and the external force’s full phase t, we may rewrite Eq. (5.42), describing the pendulum, in a form similar to Eq.
(1), i.e. as a system of three first-order ordinary differential equations:
q p
0
,
p
sin q 2 p ( f / ) cos ,
