Hydrodynamical equations
Derivation and simple solutions
Jiˇr´ı Krtiˇcka
Masaryk University
Derivation of hydrodynamical equations
Boltzmann equation
Particle distribution function F(t, x, ξ) gives the number of particles in
the element of the phase space dx dξ = dx1 dx2 dx3 dξ1 dξ2 dξ3 with
coordinates x and momenta ξ as
F(t, x, ξ) dx dξ.
The time evolution of the particle distribution function under the
influence of external force f acting on partice with mass m and taking
into account particle collisions is
∂F
∂t
+
ξh
m
∂F
∂xh
+ fh
∂F
∂ξh
=
dF
dt coll
,
which is the Boltzmann equation. Here used the Einstein summation
convention for index h.
1
Boltzmann equation
Using the Poisson bracket
{H, F} =
∂H
∂xh
∂F
∂ξh
−
∂H
∂ξh
∂F
∂xh
,
the Boltzmann equation for the system that obeys the Hamilton equation
can be rewritten as
∂F
∂t
− {H, F} =
dF
dt coll
.
For stationary collisionless system the distribution function depends on
the particle energy only,
{H, F} = 0.
2
Momentum equations
The Boltzmann equation can be solved numerically to derive the particle
distribution function. However, for most of practical applications, the
distribution function is very close to the Maxwelian distribution expressed
at given location in the frame comoving with the fluid. In such a case,
just mean quantities are of real importance for the description of the
flow. These are moments of the Boltzmann equation
m F dξ = ρ, (0th moment, flow density),
1
m
ξF dξ = v, (1st moment, flow velocity).
These can be derived by multiplying the Boltzmann equation by m and
ξ/m and by integrating. However, the equation for n-th moment contains
n + 1-th moment. Consequently, we shall close the equations somehow to
avoid obtainig infinite set of equations. This is done for the equation for
the 2nd moment using thermodynamical relations for pressure.
3
The continuity equation
Multiplicating the Boltzmann equation by particle mass m and
integrating over the velocity space
m
∂F
∂t
dξ
1
+ m
ξh
m
∂F
∂xh
dξ
2
+ m fh
∂F
∂ξh
dξ
3
= m
dF
dt coll
dξ
4
1 = m
∂
∂t
F dξ = m
∂n
∂t
=
∂ρ
∂t
,
2 =
∂
∂xh
ξhF dξ = m
∂
∂xh
(nvh) =
∂ (ρvh)
∂xh
,
3 = fh[F]∞
−∞ dξ′
= 0 is fh does not depend on ξ,
4 = 0 for conserved quantity (m),
where
• n = F dξ is number density of particles,
• ρ = mn is the density,
• vh = 1
N ξhF dξ is the mean speed.
4
The continuity equation
This gives
∂ρ
∂t
+
∂ (ρvh)
∂xh
= 0,
or
∂ρ
∂t
+ ∇ · (ρv) = 0,
which is the continuity equation.
5
The continuity equation: interpretation
Integration over volume fixed in space gives
−
V
∂ρ
∂t
dV =
V
∇ · (ρv) dV ,
or, using the Stokes theorem
−
d
dt V
ρ dV =
∂V
ρv dS,
which is the expression of the law of conservation of mass.
V
6
The continuity equation: Lagrangian picture
Introducing the Lagrangian derivative, describing the time change of any
quantity q(t, x) following a moving fluid particle,
Dq(t, x)
Dt
=
∂q(t, x)
∂t
+
∂q(t, x)
∂xh
∂xh
∂t
,
D
Dt
≡
∂
∂t
+ v · ∇
the continuity equation can be rewritten as
Dρ
Dt
+ ρ∇ · v = 0,
which for incompressible fluid (ρ = const.) is
∇ · v = 0.
7
Equation of motion
Multiplicating the Boltzmann equation by ξi and integrating
ξi
∂F
∂t
dξ
1
+ ξi
ξh
m
∂F
∂xh
dξ
2
+ ξi fh
∂F
∂ξh
dξ
3
= ξi
dF
dt coll
dξ
4
1 =
∂
∂t
ξi F dξ = m
∂
∂t
(nvi ) =
∂ (ρvi )
∂t
,
2 =
1
m
∂
∂xh
ξi ξhF dξ = m
∂
∂xh
(ci + vi )(ch + vh)F dξ =
m
∂
∂xh
vi vh F dξ + vh ci F dξ + vi chF dξ + ci chF dξ =
∂
∂xh
(mnvi vh + 0 + 0 + phi ) =
∂
∂xh
(ρvi vh + phi ) ,
3 = h fh[ξi F]∞
−∞ dξ′
− h δihfhF dξ = −nfi = −ρgi ,
4 = 0 for conserved quantity (ξ), where
• ch = ξh/m − vh is the thermal speed,
• phi = m ci chF dξ is the pressure tensor, phi = pδhi ,
• gi = fi /m is force per unit of mass (acceleration). 8
Equation of motion
This gives
∂ (ρvi )
∂t
+
∂
∂xh
(ρvi vh + p δhi )
Πik
= ρgi ,
which is, after differencing and using the continuity equation,
ρ
∂vi
∂t
+ ρvh
∂vi
∂xh
= −
∂p
∂xi
+ ρgi ,
where Πik is the momentum flux density tensor, or
ρ
∂v
∂t
+ ρv · ∇v = −∇p + ρg,
the momentum equation. Introducing the Lagrangian derivative the
momentum equation has a form of Newton’s second law
ρ
Dv
Dt
= −∇p + ρg.
9
Energy equation
Multiplicating the Boltzmann equation by ξi ξj /m and integrating
1
m
ξi ξj
∂F
∂t
dξ
1
+
1
m2
ξi ξj ξh
∂F
∂xh
dξ
2
+ ξi ξj
fh
m
∂F
∂ξh
dξ
3
=
1
m
ξi ξj
dF
dt coll
dξ
4
1 =
1
m
∂
∂t
ξi ξj F dξ = m
∂
∂t
(ci + vi )(cj + vj )F dξ =
∂
∂t
(ρvi vj + pij) ,
2 =
1
m2
∂
∂xh
ξi ξj ξhF dξ = m
∂
∂xh
(ci + vi )(cj + vj )(ch + vh)F dξ =
∂
∂xh
(ρvi vj vh + vhpij + vi phj + vj phi ) ,
3 =
0, terms with h = i and h = j (direct integration),
−fi nvj − fj nvi , terms with h = i or h = j (per-partes),
4 = 0 when contraction is performed, where
• phij = chci cj F dξ/m is phij = 0 when neglecting viscosity.
10
Energy equation
After the contraction and multiplication by 1
2 we derive
∂
∂t
1
2
ρv2
+
3
2
p +
∂
∂xh
1
2
ρvhv2
+
5
2
pvh − ρvi gi = 0,
or, introducing the specific energy ρǫ = 3
2 p,
∂
∂t
ρǫ +
ρv2
2
+ ∇ · ρv ǫ +
v2
2
+ pv − ρvg = 0,
which is the energy equation.
11
Energy equation: some manipulations
Multiplication of momentum equation by vi and summation gives
ρvi
∂vi
∂t
+ ρvi vh
∂vi
∂xh
= −vi
∂p
∂xi
+ vi ρgi ,
or
∂
∂t
1
2
ρv2
−
1
2
v2 ∂ρ
∂t
+
∂
∂xh
1
2
ρv2
vh −
1
2
v2
=0 (continuity equation)
∂ (ρvh)
∂xh
= −vi
∂p
∂xi
+vi ρgi .
Substracting this from the energy equation yields equation for the
internal energy
∂ (ρǫ)
∂t
+ ∇ · (ρǫv) = −p∇ · v,
which can be rewritten using the continuity equation as
ρ
Dǫ
Dt
= −p∇ · v.
12
Energy equation: second law of thermodynamics
The conservation of entropy for isentropic flow requires that
Ds
Dt
= 0,
which for the specific entropy of ideal gas s = cV ln(pvκ
) + const.
= cV ln(pρ−κ
) + const. is (using p = 2
3 ρǫ)
∂ (ρǫρ−κ
)
∂t
+ v · ∇ ρǫρ−κ
= 0.
Derivating and multiplying by ρκ
we arrive at
∂ (ρǫ)
∂t
+ ∇ · (ρǫv) − ρǫ∇ · v − κǫ
∂ρ
∂t
− κǫv · ∇ρ = 0.
Eliminating the last two terms using the equation of continuity and
noting that κ − 1 = 2
3 for ideal gas we derive the equation for the
internal energy once again
∂ (ρǫ)
∂t
+ ∇ · (ρǫv) = −p∇ · v.
13
Many faces of the beast
Collecting the nuggets: the hydrodynamical equations
∂ρ
∂t
+ ∇ · (ρv) = 0
ρ
∂v
∂t
+ ρv · ∇v = −∇p + ρg,
∂
∂t
ρǫ +
ρv2
2
+ ∇ · ρv ǫ +
v2
2
+ pv = ρvg
- system of nonlinear first-order partial differential equations
- unknowns ρ, v, p, and ǫ (+equation of state)
- initial and boundary conditions crucial
- inviscid flow, no magnetic field
- some special analytic solutions, general solution only numerically
- stationary solutions are important (∂/∂t = 0, but v = 0)
14
The hydrodynamical equation in planar symmetry
In a planar symmetry the hydrodynamic quantities do not depend on x
and y coordinates, there is no flow in x and y directions (v = v(z)z)
and the hydrodynamical equations are
∂ρ
∂t
+
∂
∂z
(ρv) = 0,
∂v
∂t
+ v
∂v
∂z
= −
1
ρ
∂p
∂z
+ g.
15
The hydrodynamical equations in spherical coordinates
In spherical coordinate system, the components of the velocity vector are
v = (vr , vθ, vφ) and the components of force are g = (gr , gθ, gφ). The
equation of continuity is
∂ρ
∂t
+
1
r2
∂
∂r
(r2
ρvr ) +
1
r sin θ
∂
∂θ
(sin θρvθ) +
1
r sin θ
∂
∂φ
(ρvφ) = 0
and the components of equation of motion take the form of
∂vr
∂t
+ vr
∂vr
∂r
+
vθ
r
∂vr
∂θ
+
vφ
r sin θ
∂vr
∂φ
−
v2
θ + v2
φ
r
= −
1
ρ
∂p
∂r
+ gr ,
∂vθ
∂t
+ vr
∂vθ
∂r
+
vθ
r
∂vθ
∂θ
+
vφ
r sin θ
∂vθ
∂φ
+
vr vθ
r
−
v2
φ cot θ
r
= −
1
rρ
∂p
∂θ
+ gθ,
∂vφ
∂t
+vr
∂vφ
∂r
+
vθ
r
∂vφ
∂θ
+
vφ
r sin θ
∂vφ
∂φ
+
vr vφ
r
+
vθvφ cot θ
r
= −
1
rρ sin θ
∂p
∂φ
+gφ.
16
The hydrodynamical equations in spherical symmetry
In a spherical symmetry the hydrodynamic quantities do not depend on θ
and φ coordinates, there is no flow in θ and φ directions (v = v(r)r) and
the hydrodynamical equations are (v ≡ vr )
∂ρ
∂t
+
1
r2
∂
∂r
(r2
ρv) = 0,
∂v
∂t
+ v
∂v
∂r
= −
1
ρ
∂p
∂r
+ g.
17
The hydrodynamical equations in cylindrical coordinates
In cylindrical coordinate system, the components of the velocity vector
are v = (vR , vφ, vz ) and the components of force are g = (gR , gφ, gz).
The equation of continuity is
∂ρ
∂t
+
1
R
∂
∂R
(RρvR ) +
1
R
∂
∂φ
(ρvφ) +
∂
∂z
(ρvz ) = 0
and the components of equation of motion take the form of
∂vR
∂t
+ vR
∂vR
∂R
+
vφ
R
∂vR
∂φ
+ vz
∂vR
∂z
−
v2
φ
R
= −
1
ρ
∂p
∂R
+ gR,
∂vφ
∂t
+ vR
∂vφ
∂R
+
vφ
R
∂vφ
∂φ
+ vz
∂vφ
∂z
+
vR vφ
R
= −
1
ρR
∂p
∂φ
+ gφ,
∂vz
∂t
+ vR
∂vz
∂R
+
vφ
R
∂vz
∂φ
+ vz
∂vz
∂z
= −
1
ρ
∂p
∂z
+ gz.
18
Hydrostatic equilibrium
Static case: ∂/∂t = 0, v = 0
In a static case the equation of continuity is fulfilled identically and the
momentum equation leads to
∇p = ρg
the equation of hydrostatic equilibrium. The energy equation Qrad = 0
gives the radiative equilibrium equation.
19
Atmosphere in hydrostatic equilibrium
The equation of hydrostatic equilibrium in homogeneous gravitational
field directed along the z-axis is
dp
dz
= −ρg,
which, using the ideal gas equation of state p = ρkT/(µmH), leads to
d (ρT)
dz
= −
µgmH
k
ρ.
In isothermal atmosphere T = const. this has the solution
ρ = ρ0e−z/H
, H =
kT
µmHg
,
where H is the atmospheric scale-height. For z → ∞ is ρ → 0, as it
should be.
20
Density scale height: sweeping the whole Universe
Because many astrophysical object are close to the hydrostatic
equilibrium, the density scale height is one of the most important
charasteristics length scales in astrophysics. Besides others, it determines
the density scale of atmosphere.
T [K] g/g⊙ H
neutron star 106
1010
3 mm
white dwarf 104
103
100 m
Earth 300 4 × 10−2
10 km
Sun 6000 1 200 km
Supergiant 104
10−3
1 R⊙
galaxy cluster 107
10−13
105
pc
H =
kT
µmHg
= 30 km
T
103 K
g
g⊙
−1
21
Reverting the equation: temperature estimate
Reverting the equation to estimate the temperature, we get
T
103 K
=
H
30 km
g
g⊙
.
Assuming that the coronal plasma is in hydrostatic equilibrium,
H ≈ 0.1 R⊙ gives an estimate of the coronal temperature of 2 × 106
K.
22
The mass dependence
The density scale height depends on mass of particles
H =
kT
µmHg
.
This explains variations of He abundance in closed loops observed using
HERSCHEL experiment (Moses et al. 2020, Ofman et al. 2024). 23
Atmosphere in hydrostatic equilibrium: spherical symmetry
The equation of hydrostatic equilibrium in spherically symmetric
isothermal case is
dp
dr
= −ρg,
which, with g = GM/r2
, has the solution
ρ = ρ0 exp
µmHGM
kT
1
r
.
There are two problems with this solution applied for gas spheres. For
r → 0 the equation is not applicable, because one should insert
M = M(r): Lane-Emden equation. Moreover, for r → ∞ is
ρ → ρ0 = ρISM. Solution: Bonnor-Ebert spheres with external pressure.
Matter may escape from the regions, where the thermal speed is higher
than the escape speed: atmospheric escape: loss of planetary
atmospheres, solar-type (coronal) winds.
24
Lane-Emden equation
Consider a spherical mass in equilibrium. The hydrostaic equilibrium
equation is
dp
dr
= −ρg = −
ρGM(r)
r2
.
The polytropic relation p = Cρ1+1/n
with the definition of mass inside
radius r, which is M(r) = 4π
r
0 ρr′2
dr′
, gives after differentiation
1
r2
d
dr
r2
ρ
d
dr
ρ1+1/n
= −
4πG
C
ρ.
Introducing new variables θ and ξ via ρ = λθn
and ξ = r/α, where λ is
arbitrary dimensional constant and
α =
C(1 + n)
4πGλ1+1/n
we arrive at Lane-Emden equation
1
ξ2
d
dξ
ξ2 dθ
dξ
= −θn
.
25
Hydrostatic atmospheres with radiative force
The equation of hydrostatic equilibrium in spherically symmetric
atmosphere in radiative equlibrium is
dp
dr
= −ρg + ρgrad,
with the radiative force grad = 1
c κρFν dν = κρL
4πcr2 . The temperature is
given by the energy transport equation
dT
dr
= −
3
4aT3
κρL
4πcr2
.
This can be rewritten in terms of prad = (a/3)T4
as
dprad
dr
=
4aT3
3
dT
dr
= −
κρL
4πcr2
= −ρgrad.
Therefore, the equation of hydrostatic equilibrium is
dptot
dr
=
d(p + prad)
dr
= −ρg.
26
Atmospheres close to the Eddington limit
Dividing the last two equations
dptot
dprad
=
dp
dprad
+ 1 =
g
grad
=
1
Γ
, or
dp
dprad
=
1
Γ
− 1,
where the generalized Eddington factor is
Γ =
κρL
4πcGM
.
The derivative dp/dprad is a function of Γ only. Moreover, the the point
at which the envelope solution crosses the Eddington limit Γ = 1 needs
to be an extremum in p (Gr¨afener et al. 2012).
27
Envelope inflation close to the Eddington limit
6 R
4 R
2 R
7
6
5
4
3
8 7 6 5 4 3
log10 (Prad /dyn cm-2
)
log10(P/dyncm-2
)
1.0
0.5
0.0
-0.5
-1.0
Pgas
Prad
Ptot
Logarithm of the Eddingtor factor Γ
(colors) in the prad − p plane with
Iglesias & Rogers (1996) opacities.
Black arrows denote slopes
dp/dprad.
The numerical solution for 23 M⊙ star almost precisely follows a path
with Γ = 1 and crosses the Eddington limit at the lowest gas pressure
(corresponding to the Fe-opacity peak). The gas density increases
outwards leading to density inversion. This explains why WR and LBV
stars have extended envelopes.
(Gr¨afener et al. 2012) 28
Eddington limit
The Eddington parameter due to pure electron scattering can be
evaluated as
Γ =
σT
ne(r)
ρ(r) L
4πcGM
, in scaled quantities, Γ ≈ 10−5 L
1 L⊙
M
1 M⊙
−1
.
Therefore, our Sun is very far from the Eddington limit Γ = 1, while hot
O stars with L ≈ 106
L⊙ and M ≈ 10 M⊙ are very close to the Eddington
limit Γ ≈ 0.1.
From the Eddington limit Γ = 1 one can estimate the minimum mass of
radiating object that can be gravitationally bound. From the above
equation
MEdd ≈ 10−5
M⊙
L
1 L⊙
.
The observed bolometric flux from quasars of the order of 1013
L⊙ yields
the Eddington mass of 108
M⊙ providing a hint that the objects are
powered by a supermassive black hole.
29
Suggested reading
G. K. Batchelor: An Introduction to Fluid Dynamics
D. Mihalas & B. W. Mihalas: Foundations of Radiation Hydrodynamics
F. H. Shu: The physics of astrophysics: II. Hydrodynamics
A. Feldmeier: Theoretical Fluid Dynamics
30