Flow Dynamics, Lecture Notes- Physics 3, Study notes of Physics

Download Flow Dynamics, Lecture Notes- Physics 3 and more Study notes Physics in PDF only on Docsity!

BASIC CONCEPTS

IN COMPUTATIONAL FLUID

DYNAMICS

finite differences

time-marching

stability

false diffusion

false convection

w

iggles; cell Reynolds number

CFD is the art of simulating numerically the motion of fluid. Thekey difficulty is that in reality fluid motion takes place in aCONTINUUM described by continuous DIFFERENTIALEQUATIONS. A computer can only handle a limited quantity ofnumbers (2Mb stores 500,000 32bit numbers) and perform basicarithmetical operations. The art of CFD is to construct anappropriate DISCRETE REPRESENTATION of the continuum.This construction is governed by:

•PHYSICAL BEHAVIOUR of the continuum equations –conservation; the possible presence of characteristics or wave-like behaviour; whether the problem is “initial value” or“boundary value”.•NUMERICAL BEHAVIOUR of the discrete representation –accuracy; stability; false diffusion false convection.There are four main methodologies adopted in constructing thenumerical models

• finite difference
• finite volume

-spectral

• finite element

The main ideas can be most easily demonstrated using the finite difference approach.

COMPUTATIONAL FLUID DYNAMICS

It is clear that both differences are “consistent” in the sense that as

X

0 the desired

differential is recovered. However, for practical computations

X is very much finite (of

order 10

chord lengths!) and the two difference equations have very different properties.

The leading error term in the expansion is retained and called the “truncation error”; one-sided dif

f

erences are “first order accurate” as the truncation error is O(

X) whereas centred

differences are “second order accurate” as the error term is O(

X

2

Thus,•

ONE-SIDED DIFFERENCE

CENTRED DIFFERENCE

DESIRED

TRUNCATION ERROR

−

2 2

1 i

i

X

f

X 2

f X

X

f

f

−

3 3 2 1 i 1 i

X

f

X 6

f X

X

f

f

A one-sided difference replacement for a sense derivative is Note how the lack of symmetry in the one-sided difference has caused the truncation errorto deteriorate Appropriate difference operators can be derived in a similar way for any requireddifferential. For example, a centred difference replacement for a second derivative is

3 3

2 2

2

2 i

1 i

i

X

f

X

X

f

X

f

f 2

f

−

−

4 4 2 2 2 2

1 i

i

1 i

X

f

X

X

f

X

f

f 2

f

−

If we use this “dif

f

erence molecule” and choose a centred representation for

2

u/

X

2

and a

one-sided representation for

u

t we have

Time level

Space location

2

n 1 i i 1 i

n i

1 n i

X

U

U

U

t

U

U

υ

−

This expression is “first order accurate in time and second order accurate in space”. The newvalue of U

i^

can be found EXPLICITLY as all other references to U

i^

are at the old time level.

All we have to do is choose a value of

t and “Time-March” the solution forward: see Fig 1.

i = 21

X

5 9

X

t^2

=

Δ

Δ ν

/ x (^0510)

X t^

2

=

Δ

ν

u

i = 1

(^511)

X

t^2

=

Δ

Δ ν

2 (^2) X u

u t

∂ ∂ ν =

∂ ∂

(^

n) 1 i i 1 i 2

ni

1 ni^

u

u 2

u

X

t

u

u

−

−

Δ

Δ ν

=

0

BOUNDARY

u(1) = 0

21 space nodesCONDITIONS u(21) = 0

i i + 1

Δ

X

5

X t^ 2

=

Δ

ν (^1020)

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

x

Fig. 1: 1D Heat Flow

X

MOVIE

MOVIE

10

We imagine the finite difference solution to be composed of two parts, the exact solution,

and

a perturbation,

ε

i

We substitute this into the finite difference equations and follow the growth or decay of

ε

This approach leads to strong physical insight into how the numerical model mimics theunderlying physics of the governing differential equation.For example: introduce a point disturbance at point i

DISCRETE PERTURBATION STABILITY ANALYSIS

i U

i

i

i^

U

U

ε

ε i

[

]

[

]

Δ υ − = ε ∴

ε

Δ υ + ε = ⇒

ε + = + − Δ

υ

−

2

1 n i

2

1 n i

i

i

i

n

1 i

n i

n

1 i

2

n i

1 n i

X

t

U

X

t

U

U U ; U U 2 U X

t

U

U

For stability,

X

t

e . i

X

t

Thus

u

2

2

1

n i

ν

Δ ν − ≤ − ≤ ε +

Note: if we require no overshoot (ie no non-physical negative values of U) we musthave

X

t^2

ν

Substituting the “wiggle” disturbance into the finite difference equation gives: This shows that for a given

ν

and

X, the stable choice of

t is restricted; Fig.l

shows solution attempts either side of the stability boundary. This stability bound is very restrictive since if we half the mesh spacing (to checkaccuracy for example) then

t is reduced by a factor of four and so we must perform

four times as many time marching steps, each costing twice as much, to get to thesame physical time – an increase in cost of factor 8!

X

t

Hence

1 X / t 4 1 1

U

stability

for

X / t 4 1 U. e

i

X

t

U

2

(^12) n i

2

1 n i

2

1 n i

ν

• ≤ Δ Δ ν − ≤ − ∴

ε

Δ Δ ν − = ε

ε + ε − − ε Δ

Δ ν = ε − −

14

VON NEUMANN STABILITY ANALYSIS In this common analysis we

represent

the variation of the computed variable by a Fourier

Series:

SHORTEST WAVELENGTH

π

=

λ k

x

X

LONGEST WAVELENGTH

λ

k

n i

U

ˆ^ kXj

k

n

n i^

e

U

U

U

n

is the Fourier coeff. @ time n

k = wave number = 2

π

/wavelength

X = i

λ

X with I the space increment

ˆj

16

The stability analysis is linear so we consider one harmonic

i ˆj

n

kXj ˆ

n

n i^

e

U

e

U

u

θ

with

θ

= k

X, the “phase angle” and substitute into the difference equation:

(^

)^

(^

)

(

)

(

)

cos

X

t

e

e

X

t

U

U

e e 2 e U X

t

e

U

e

U

2

ˆj

ˆj

2

n

1 n

1 i ˆj

i ˆj

1 i ˆj

n

2

i ˆj

n

i ˆj

1 n

θ

Δ ν + = λ ∴

Δ ν + = λ = ∴

ν

θ −

θ

− θ θ + θ θ θ +

AMPLIFICATION FACTOR

For stability :

π ≤ θ ≤ ≤ λ

for

for

all

harmonics

Hence,

before

as

X

t^2

ν

17

If we are only interested in steady solutions and we regard time marching as simply aniterative strategy we might reasonably ask what is the best value of

t to choose. It is not

always the case that the bigger

t the faster the “convergence rate” (i.e. the fewer the steps

needed to attain a steady solution). We assess the convergence rate by examining themagnitude of the von Neumann

λ

which, if less than unity, gives the rate at which errors

associated with a given harmonic are damped. CONVERGENCE RATE

cos

X

t

have

we

above

From

2

θ

ν

λ

π

θ

λ^ 1. 0.

SHORTEST WAVELENGTH ERRORS

2

X

t 4 1 X x 2

x

K

Δ ν − = λ ∴ = Δ ⋅ Δ

π

θ

π

-thus for a typical choice of

X

t^2

λ

ν

π

This means that “wiggles” arevigorously damped!

2D EXAMPLE^ Consider the 2D heat conduction equation^ A numerical representation is

ν

2 2

2 2

Y

u

X

u

u t

n

2

1 ij

ij

1 ij

2

j 1 i

ij

j 1 i

n ij

1 n ij

Y

u

u 2

u

X

u

u 2

u

t

u

u

ν

=

−

−

i+1j

ij-

i-1j

ij+

t Δ

Y Δ

X Δ

n + 1 n

ij

t

Y

X

We analyse for stability using the 2D form of the von-Neumann method:

(^

) η

ξ

j

i ˆj

n

n i^

e

U

u

no

wave

wise

Y

k

Y

k

no

wave

wise

X

k

X

k

Y j Y , X i X , 1

ˆj

Y

Y

X

X

η

ξ

Substituting this trial function: The restrictive stability bounds of “explicit” time-marching can, under some circumstances,be relieved by “implicit” techniques, i.e. by updating the variables simultaneously

(

)

η

ξ

Δ ν ⋅ + = λ = ∴

− + Δ ν = − +

η − η ξ − ξ + ξ +

η

cos

X Y

cos

X

t

U

U

Y

e

e

X

e

e

e

tU

e

U

e

U

2

2

n

1

n

ˆj 2

ˆj

ˆj 2

ˆj

...j ˆ

n

...j ˆ

n

j

i ˆj

1

n

This can be analysed for stability etc… In particular, if

X=

Y,

⏐λ⏐≤

for

which is twice as restrictive as in 1D!

X

t^2

ν IMPLICIT TIME-MARCHING

i+

i-

i

n+

Δ

t

Δ

X

n

1 n 1 i i 1 i 2

ni

1

ni

2 2

u

u 2

u

X

t

u

u

X

u

u t

−

ν

ν

=

t

X