Fluid Dynamics: Theory and Computation, Study notes of Fluid Dynamics

Download Fluid Dynamics: Theory and Computation and more Study notes Fluid Dynamics in PDF only on Docsity!

Fluid Dynamics:

Theory and Computation

Dan S. Henningson

Martin Berggren

August 24, 2005

• 1 Derivation of the Navier-Stokes equations
  • 1.1 Notation
  • 1.2 Kinematics
  • 1.3 Reynolds transport theorem
  • 1.4 Momentum equation
  • 1.5 Energy equation
  • 1.6 Navier-Stokes equations
  • 1.7 Incompressible Navier-Stokes equations
  • 1.8 Role of the pressure in incompressible flow
• 2 Flow physics
  • 2.1 Exact solutions
  • 2.2 Vorticity and streamfunction
  • 2.3 Potential flow
  • 2.4 Boundary layers
  • 2.5 Turbulent flow
• 3 Finite volume methods for incompressible flow
  • 3.1 Finite Volume method on arbitrary grids
  • 3.2 Finite-volume discretizations of 2D NS
  • 3.3 Summary of the equations
  • 3.4 Time dependent flows
  • 3.5 General iteration methods for steady flows
• 4 Finite element methods for incompressible flow
  • 4.1 FEM for an advection–diffusion problem
    • 4.1.1 Finite element approximation
    • 4.1.2 The algebraic problem. Assembly.
    • 4.1.3 An example
    • 4.1.4 Matrix properties and solvability
    • 4.1.5 Stability and accuracy
    • 4.1.6 Alternative Elements, 3D
  • 4.2 FEM for Navier–Stokes
    • 4.2.1 A variational form of the Navier–Stokes equations
    • 4.2.2 Finite-element approximations
    • 4.2.3 The algebraic problem in 2D
    • 4.2.4 Stability
    • 4.2.5 The LBB condition
    • 4.2.6 Mass conservation
    • 4.2.7 Choice of finite elements. Accuracy
• A Background material
  • A.1 Iterative solutions to linear systems
  • A.2 Cartesian tensor notation
    • A.2.1 Orthogonal transformation
    • A.2.2 Cartesian Tensors
    • A.2.3 Permutation tensor
    • A.2.4 Inner products, crossproducts and determinants
    • A.2.5 Second rank tensors
    • A.2.6 Tensor fields
    • A.2.7 Gauss & Stokes integral theorems
    • A.2.8 Archimedes principle
  • A.3 Curvilinear coordinates
• B Recitations 5C1214
  • B.1 Tensors and invariants
  • B.2 Euler and Lagrange coordinates
  • B.3 Reynolds transport theorem and stress tensor
  • B.4 Rankine vortex and dimensionless form
  • B.5 Exact solutions to Navier-Stokes
  • B.6 More exact solutions to Navier-Stokes
  • B.7 Axisymmetric flow and irrotational vortices
  • B.8 Vorticity equation, Bernoulli equation and streamfunction
  • B.9 Flow around a submarine and other potential flow problems
  • B.10 More potential flow
  • B.11 Boundary layers
  • B.12 More boundary layers
  • B.13 Introduction to turbulence
  • B.14 Old Exams
• C Study questions 5C1212

Preface

These lecture notes has evolved from a CFD course (5C1212) and a Fluid Mechanics course (5C1214) at

the department of Mechanics and the department of Numerical Analysis and Computer Science (NADA)

at KTH. Erik St˚alberg and Ori Levin has typed most of the L

A T E

Xformulas and has created the electronic

versions of most figures.

Stockholm, August 2004

Dan Henningson

Martin Berggren

In the latest version of the lecture notes study questions for the CFD course 5C1212 and recitation

material for the Fluid Mechanics course 5C1214 has been added.

Stockholm, August 2005

Dan Henningson

Chapter 1

Derivation of the Navier-Stokes

equations

1.1 Notation

The Navier-Stokes equations in vector notation has the following form

∂u

∂t

• (u · ∇) u = −

ρ

∇p + ν∇

2

u

∇ · u = 0

where the velocity components are defined

u = (u, v, w) = (u 1 , u 2 , u 3

the nabla operator is defined as

∂x

∂y

∂z

∂x 1

∂x 2

∂x 3

the Laplace operator is written as

2

=

2

∂x

2

2

∂y

2

2

∂z

2

and the following definitions are used

ν − kinematic viscosity

ρ − density

p − pressure

see figure 1.1 for a definition of the coordinate system and the velocity components.

The Cartesian tensor form of the equations can be written

∂u i

∂t

• u j

∂u i

∂x j

ρ

∂p

∂x i

• ν

2

u i

∂x j ∂x j

∂u i

∂x i

where the summation convention is used. This implies that a repeated index is summed over, from 1

to 3, as follows

u i u i = u 1 u 1

• u 2 u 2
• u 3 u 3

Thus the first component of the vector equation can be written out as

PSfrag replacements

x, x 1

y, x 2

z, x 3

u, u 1

v, u 2

w, u 3

Figure 1.1: Definition of coordinate system and velocity components

PSfrag replacements

x, x 1

y, x 2

z, x 3

u, u 1

v, u 2

w, u 3

x 1

x 2

x

0

i

r i

x

0

i

t

u i

x

0

i

t

P

t = 0

Figure 1.2: Particle path.

i = 1

∂u 1

∂t

• u 1

∂u 1

∂x 1

• u 2

∂u 1

∂x 2

• u 3

∂u 1

∂x 3

ρ

∂p

∂x 1

• ν

2 u 1

∂x 1

2

2 u 1

∂x 2

2

2 u 1

∂x 3

2

or

∂u

∂t

• u

∂u

∂x

• v

∂u

∂y

• w

∂u

∂z

ρ

∂p

∂x

• ν

2

u

∂x

2

2

u

∂y

2

2

u

∂z

2

∇

2 u

1.2 Kinematics

Lagrangian and Euler coordinates

Kinematics is the description of motion without regard to forces. We begin by considering the motion of a

fluid particle in Lagrangian coordinates, the coordinates familiar from classical mechanics.

Lagrange coordinates: every particle is marked and followed in flow. The independent variables are

x i

0

− initial position of fluid particle

t − time

where the particle path of P, see figure 1.2, is

r i

= r i

x i

0

,

t

and the velocity of the particle is the rate of change of the particle position, i.e.

PSfrag replacements

x, x 1

y, x 2

z, x 3

u, u 1

v, u 2

w, u 3

x 1

x 2

x

0

i

r i

x

0

i

t

u i

x

0

i

t

P

t = 0

x

u 1

u 2

u 1

Figure 1.3: Acceleration of fluid particles in converging river.

PSfrag replacements

x, x 1

y, x 2

z, x 3

u, u 1

v, u 2

w, u 3

x 1

x 2

x

0

i

r i

x

0

i

t

u i

x

0

i

t

P

t = 0

x

u 1

u 2

u 1

x 1

x 2

x 3

x

0

i

u i

d

t

dx

0

i

d

u i

d

t

ri

d

t

d

r

i

d

t

P

1

P

2

Figure 1.4: Relative motion of two nearly particles.

where we have used

t

( dr i

t

∂r i

∂x j

0

dx j

0

=

∂x j

0

∂r i

t

dx j

0

∂u i

∂x j

0

dx j

0

= du i

and where

∂u i

∂x j

0

− change of u i

with initial pos.

dx j

0 − difference in initial pos.

du i

− difference in velocity

We can transform the expression

t

( dr i

) = du i

in Lagrange coordinates to an equation for a material

line element in Euler coordinates

D

Dt

( dr i

) = du i

= {expand in Euler coordinates}

∂u i

∂x j

dr j

where

∂u i

∂x j

− change in velocity with spatial position

dr j − difference in spatial pos. of particles

Relative motion associated with invariant parts

We consider the relative motion du i

∂u i

∂x j

dr j

by dividing

∂u i

∂x j

in its invariant parts, i.e.

∂u i

∂x j

= ξ ij

• e ij

• e ij

e ij

where e ij

is the deformation rate tensor and

ξ ij

∂u i

∂u j

∂u j

∂x i

anti-symmetric part

e ij

∂u i

∂x j

∂u j

∂x i

∂u r

∂u r

δ ij

traceless part

e ij

∂u r

∂u r

δ ij isotropic part

The symmetric part of

∂u i

∂x j

, e ij

, describes the deformation and is considered in detail below, whereas the

anti-symmetric part can be written in terms of the vorticity ω k

and is associated with solid body rotation,

i.e. no deformation.

The anti-symmetric part can be written

We have dropped the quadratic terms since we are assuming that dt is small. dR

1

becomes

dR

1

= R

∂u 1

∂x 1

dt − R = R

∂u 1

∂x 1

dt + ...

− R

= R

∂u 1

∂x 1

dt = Re 11 dt

which implies that

dR

1

dt

= Re 11

Thus the deformation rate of side 1 depends on e 11 , both traceless and isotropic part of

∂u i

∂x j

Second, we consider the deformation of the angle between side 1 and side 2. This can be expressed as

cos

π

• dφ 12

r

1

· r

2

|r

1 | |r

2 |

= r

1

k

r

2

k

r

1

m

r

1

m

· r

2

n

r

2

n

− 1 / 2

δ k 1

∂u k

∂x 1

dt

δ k 2

∂u k

∂x 2

dt

∂u 1

∂x 1

dt

− 1 / 2

∂u 2

∂x 2

dt

− 1 / 2

∂u 1

∂x 2

dt +

∂u 2

∂x 1

dt

∂u 1

∂x 1

dt

∂u 2

∂x 2

dt

∂u 1

∂x 2

∂u 2

∂x 1

dt = 2 e 12 dt

where we have dropped quadratic terms. We use the trigonometric identity

cos

π

• dϕ 12

= cos

π

· cos dϕ 12

− sin

π

· sin dϕ 12

≈ − dϕ 12

which allow us to obtain the finial expression

dϕ 12

dt

= − 2 e 12

Thus the deformation rate of angle between side 1 and side 2 depends only on traceless part of

∂u i

∂x j

Third, we consider the deformation of the volume of the cube. This can be expressed as

dV =

r

1

r

2

r

3

− R

3

= R

3

∂u 1

∂x 1

dt

∂u 1

∂x 2

dt

∂u 1

∂x 3

dt

∂u 2

∂x 1

dt 1 +

∂u 2

∂x 2

dt

∂u 2

∂x 3

dt

∂u 3

∂x 1

dt

∂u 3

∂x 2

dt 1 +

∂u 3

∂x 3

dt

− R

3

= R

3

∂u 1

∂x 1

dt

∂u 2

∂x 2

dt

∂u 3

∂x 3

dt

− R

3

= R

3

∂u 1

∂x 1

∂u 2

∂x 2

∂u 3

∂x 3

dt

= R

3

∂u k

∂x k

dt

where we have again omitted quadratic terms. Thus we have

dV

dt

= R

3

e rr

and the deformation rate of volume of cube (or expansion rate) depends on isotropic part of

∂u i

∂x j

In summary, the motion of a fluid particle with velocity u i can be divided into the following invariant

parts

R

3

R

2

R

1

r

3

r

2

r

1

π

2

• dϕ 12

V (t)

S

(

t

)

V (t + ∆t) − V (t)

S (t + ∆t)

u

n

V (t + ∆t)

Figure 1.6: Volume moving with the fluid.

i) u i solid body translation

ii) ξ ij

1

2

∂u i

∂x j

∂u j

∂x i

1

2

kij ω k solid body rotation

iii) e ij

1

2

∂u i

∂x j

∂u j

∂x i

∂u r

∂u r

δ ij

volume constant deformation

iv) e ij

1

3

∂u r

∂u r

δ ij

volume expansion rate

1.3 Reynolds transport theorem

Volume integral following with the fluid

Consider the time derivative of a material volume integral, i.e. a volume integral where the volume is moving

with the fluid. We obtain the following expressions

D

Dt

V (t)

T

ij

dV = lim

∆t→ 0

∆t

V (t+∆t)

T

ij

(t + ∆t) dV −

V (t)

T

ij

(t) dV

= lim

∆t→ 0

∆t

V (t+∆t)

T

ij (t + ∆t) dV −

V (t)

T

ij (t + ∆t) dV

∆t

V (t)

T

ij (t + ∆t) dV −

V (t)

T

ij (t) dV

= lim

∆t→ 0

∆t

V (t+∆t)−V (t)

T

ij (t + ∆t) dV

V (t)

∂T

ij

∂t

dV

The volume in the first integral on the last line is represented in figure 1.6, where a volume element

describing the change in volume between V at time t and t + ∆t can be written as

u · n∆t = u k n k ∆t ⇒ dV = u k n k ∆t dS

V (t + ∆t) − V (t)

S (t + ∆t)

u

n

V (t + ∆t)

dS

∆

s

(

k

)

︷

︸︸

︷

n

I n

II

n

III(k)

R(n

I )

R(n

II

)

dl

Figure 1.7: Momentum balance for fluid element

R

r i

d

t

dr i

d

t

P

1

P

2

x 1

x 2

x 3

R

3

R

2

R

1

r

3

r

2

r

1

π

2

• dϕ 12

V (t)

S (t)

V (t + ∆t) − V (t)

S (t + ∆t)

u

n

V (t + ∆t)

n

R

Figure 1.8: Surface force and unit normal.

We can put the momentum conservation in integral form as follows

D

Dt

V (t)

ρu i dV =

V (t)

ρF i dV +

S(t)

R

i dS

Using Reynolds transport theorem this can be written

V (t)

ρ

Du i

Dt

dV =

V (t)

ρF i

dV +

S(t)

R

i

dS

which is Newtons second law written for a volume of fluid: mass · acceleration = sum of forces. To

proceed R i must be investigated so that the surface integral can be transformed to a volume integral. In

order to do that we have to define the stress tensor.

The stress tensor

Remove a fluid element and replace outside fluid by surface forces as in figure 1.7. Here R(n) is the surface

force per unit area on surface dS with normal n, see figure 1.8. Momentum conservation for the small fluid

particle leads to

ρ

Du i

Dt

dS dl = ρF i

dS dl + R i

n

I

j

dS + R i

n

II

j

dS +

k

R

i

n

III(k)

j

∆s

(k)

dl

Letting dl → 0 gives

0 = R

i

n

I

j

dS + R i

n

II

j

dS

Now n j = n

I

j

= −n

II

j

which leads to

R

i (n j

) = −R

i (−n j

implying that a surface force on one side of a surface is balanced by an equal an opposite surface force

at the other side of that surface. Note that it is a general principle that the terms proportional to the

volume of a small fluid particle approaches zero faster than the terms proportional to the surface area of

the particle. Thus the surface forces acting on a small fluid particle has to balance, irrespective of volume

forces or acceleration terms.

u

n

V (t + ∆t)

T

11

T 21

T

31

x 1

R(e 1

) = ( T

11

, T

21

, T

31

Figure 1.9: Definition of surface force components on a surface with a normal in the 1-direction.

P

1

P

2

x 1

x 2

x 3

R

3

R

2

R

1

r

3

r

2

r

1

π

2

• dϕ 12

V (t)

S (t)

V (t + ∆t) − V (t)

S (t + ∆t)

u

n

V (t + ∆t)

T

12

T 22

T

32

x 2

R(e 2

) = ( T

12

, T

22

, T

32

Figure 1.10: Definition of surface force components on a surface with a normal in the 2-direction.

We now divide the surface forces into components along the coordinate directions, as in figures 1.9 and

1.10, with corresponding definitions for the force components on a surface with a normal in the 3-direction.

Thus, T ij is the i-component of the surface force on a surface dS with a normal in the j-direction.

Consider a fluid particle with surfaces along the coordinate directions cut by a slanted surface, as in

figure 1.11. The areas of the surface elements are related by

dS j = e j · n dS = n j dS

where dS is the area of the slanted surface. Momentum balance require the surface forces to balance

in the element, we have

0 = R

1

dS − T 11

n 1

dS − T 12

n 2

dS − T 13

n 3

dS

which implies that the total surface force R i

can be written in terms of the components of the stress

tensor T ij

as

R

i

= T

i 1 n 1

+ T

i 2 n 2

+ T

i 3 n 3

= T

ij n j

PSfrag replacements

x, x 1

y, x 2

z, x 3

u, u 1

v, u 2

w, u 3

x 1

x 2

x

0

i

r i

x

0

i

t

u i

x

0

i

t

P

t = 0

x

u 1

u 2

u 1

x 1

x 2

x 3

x

0

i

u i d

t

dx

0

i

du i

d

t

r i

d

t

dr i

d

t

P

1

P

2

x 1

x 2

x 3

R

3

R

2

R

1

r

3

r

2

r

1

π

2

• dϕ 12

V (t)

S (t)

V (t + ∆t) − V (t)

S (t + ∆t)

u

n

V (t + ∆t)

e 1

e 2

e 3

R

R

1

n

dS 1

dS 2

dS 3

Figure 1.11: Surface force balance on a fluid particle with a slanted surface.

1.5 Energy equation

The energy equation is a mathematical statement which is based on the physical law that the rate of change

of energy in material particle = rate that energy is received by heat and work transfers by that particle.

We have the following definitions

ρ

[

e +

1

2

u i u i

]

dV energy of particle, with e the internal energy

ρu i

F

i dV work rate

︸ ︷︷ ︸

force · velocity

of F i on particle

u j

R

j

dS work rate of R j

on particle

n i q i dS heat loss from surface, with q i the heat flux vector, directed outward

Using Reynolds transport theorem we can put the energy conservation in integral form as

D

Dt

V (t)

ρ

[

e +

u i

u i

]

dV =

V (t)

ρF i

u i

dV +

S(t)

[n i

T

ij

u j

− n i

q i

] dS =

V (t)

[

ρF i

u i

∂x i

(T

ij

u j

− q i

]

dV

Compare the expression in classical mechanics, where the momentum equation is m u˙ = F and the

associated kinetic energy equation is

m

2

d

dt

(u · u) = F · u

work rate = force · velocity

( work = force · dist. )

From the integral energy equation we obtain the total energy equation by the observation that the

volume is arbitrary and thus that the integrand itself has to be zero. We have

ρ

D

Dt

e +

u i u i

= ρF i u i

∂x i

(pu i

∂x i

(τ ij u j

∂q i

∂x i

The mechanical energy equation is found by taking the dot product between the momentum equation

and u. We obtain

ρ

D

Dt

u i

u i

= ρF i

u i

− u i

∂p

∂x i

• u i

∂τ ij

∂x j

Thermal energy equation is then found by subtracting the mechanical energy equation from the total

energy equation, i.e.

ρ

De

Dt

= −p

∂u i

∂x i

• τ ij

∂u i

∂x j

∂q i

∂x i

The work of the surface forces divides into viscous and pressure work as follows

∂x i

(pu i ) = −p

∂u i

∂x i

− u i

∂p

∂x i

∂x i

(τ ij u j ) = τ ij

∂u i

∂x j

• u j

∂τ ij

∂x i

where following interpretations can be given to the thermal and the mechanical terms

© 1 : thermal terms

( force · deformation ): heat generated by compression and viscous dissipation

© 2 : mechanical terms

( velocity · force gradients ): gradients accelerate fluid and increase kinetic energy

The heat flux need to be related to the temperature gradients with Fouriers law

q i = −κ

∂T

∂x i

where κ = κ (T ) is the thermal conductivity. This allows us to write the thermal energy equation as

ρ

De

Dt

= −p

∂u i

∂x i

∂x i

κ

∂T

∂x i

where the positive definite dissipation function Φ is defined as

Φ = τ ij

∂u i

∂x j

= 2 μ

[

e ij

e ij

∂u k

∂x k

2

]

= 2 μ

e ij

∂u k

∂x k

δ ij

2

Alternative form of the thermal energy equation can be derived using the definition of the enthalpy

h = e +

p

ρ

We have

Dh

Dt

De

Dt

ρ

Dp

Dt

p

ρ

2

Dρ

Dt

︸︷︷︸

−

(

ρ

∂u i

∂x i

)

which gives the final result

ρ

Dh

Dt

Dp

Dt

∂x i

κ

∂T

∂x i

To close the system of equations we need a

i) thermodynamic equation e = e (T, p) simplest case: e = c v T or h = c p

T

ii) equation of state p = ρRT

where c v

and c p

are the specific heats at constant volume and temperature, respectively. At this time

we also define the ratio of the specific heats as

γ =

c p

c v

1.6 Navier-Stokes equations

The derivation is now completed and we are left with the Navier-Stokes equations. They are the equation

describing the conservation of mass, the equation describing the conservation of momentum and the equation

describing the conservation of energy. We have

Dρ

Dt

• ρ

∂u k

∂x k

ρ

Du i

Dt

∂p

∂x i

∂τ ij

∂x j

• ρF i

ρ

De

Dt

= −p

∂u i

∂x i

∂x i

κ

∂T

∂x i

where