Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

Solvers, schemes SIMPLEx, upwind

REDUCED lecture…

Computer Fluid Dynamics E181107CFD7r

2181106

FINITE VOLUME METHODCFD7r

x

x

y

z

zSouthWest

Top

E

North

Bottom

x

y

FINITE CONTROL VOLUME x = Fluid element dx

Finite size Infinitely small size

[ ( ) ] 0

( ) 0

CV

A CV

S dV

Gauss theorem

n dA SdV

FVM diffusion problem 2DCFD7r

, ,....

(without sources)

P P E E W W N N S S

e e w wE W

PE WP

P E W N S

a a a a a S

A Aa a

x x

a a a a a

( ) ( ) 0Sx x y y

W P e E

N

S

A

Boundary conditions

( )

|

| ( ( ) )

A

A A

A e

A

qx

k Ax

First kind (Dirichlet BC)

Second kind (Neumann BC)

Third kind (Newton’s BC)

Example: insulation q=0

Example: Finite thermal resistance (heat transfer

coefficient )

( ) 0A CV

n dA SdV

... 0E Pe e w w

PE CV

A A SdVx

FVM convection diffusion 1DCFD7r

CV - control volume

A – surface of CV

( ) [ ( ) ]

( ) ( )

CV CV

A A CV

u dV S dV

n u dA n dA S dV

Gauss theorem

1D case

A W P E B

w e

( ) ee e e

PE

F u Dx

( ) ww w w

WP

F u Dx

( ) ( )e w e Ew P w Pe WF F D D S F-mass flux D-diffusion conductance

FVM convection diffusion 1DCFD7r

Relative importance of convective transport is characterised by Peclet number of cell (of control volume, because Pe depends upon size of cell)

2

3

Re [ . ]

Re Pr [ ]/

Re [ ]

p p

m ABAB

u x F kgPe Pa s

D m s

u x k W kg K kgPe

k c c m K J m s

u x kg m kgPe Sc D

D m s m s

Prandtl

Schmidt

Binary diffusion coef

Thermal conductivity F

PeD

FVM convection diffusion 1DCFD7r

Methods differ in the way how the unknown transported values at the control

volume faces (w , e) are calculated (different interpolation techniques)

Result can be always expressed in the form

Properties of resulting schemes should be evaluated

Conservativeness

Boundedness (positivity of coefficients anb)

Transportivity (schemes should depend upon the Peclet number of cell)

P P nb nbneighbours

a a

FVM central scheme 1DCFD7r

Conservativeness (yes)

Boundedness (positivity of coefficients)

Transportivness (no)

2 2

P P W W E E

w eW w E e

a a a

F Fa D a D

1 1( ) ( )

2 2e P E w P W

| |2e

ee

FPe

D

Very fine mesh is necessary

FVM upwind 1st order 1DCFD7r

Conservativeness (yes)

Boundedness (positivity of coefficients) YES for any Pe

Transportivness (YES)

max ,0 max ,0P P W W E E

W w w E e e

a a a

a D F a D F

for flow direction to right 0 else

for 0 else e P e e E

w P w w W

F

F

But at a prize of decreased order of accuracy (only 1st

order!!)

FVM Upwind FLUENTCFD7r

= V A

dV dA

0

1f

fA

f

ff AV

1

)(2

110 f

Upwind of the first order

Upwind of the second order0f s

0 f

FVM exponential FLUENTCFD7r

For Pe<<1 the exponential profile is reduced to linear profile

0

1 0

( ) 1=

1

rPe

L

Pe

r e

e

Exact solution of )()(rr

vr r

LvPe r

r

0

1

Projection of velocity to the direction r

Called power law model in Fluent

There are two basic errors caused by approximation of convective terms

False diffusion (or numerical diffusion) – artificial smearing of jumps, discontinuities. Neglected terms in the Taylor series expansion of first derivatives (convective terms) represent in fact the terms with second derivatives (diffusion terms). So when using low order formula for convection terms (for example upwind of the first order), their inaccuracy is manifested by the artificial increase of diffusive terms (e.g. by increase of viscosity). The false diffusion effect is important first of all in the case when the flow direction is not aligned with mesh, see example

Dispersion (or aliasing) – causing overshoots, artificial oscillations. Dispersion means that different components of Fourier expansion (of numerical solution) move with different velocities, for example shorter wavelengths move slower than the velocity of flow.

FVM dispersion / diffusionCFD7r

0 - boundary condition, all

values are zero in the right triangle (without diffusion)

1Numerical diffusion

Numerical dispersion

u=v

FVM Navier Stokes equationsCFD7r

Beckman

Different methods are used for numerical solution of Navier Stokes equations at

Compressible flows (finite speed of sound, existence of velocity discontinuities at shocks)

Incompressible flows (infinite speed of sound)

Explicit methods (e.g. MacCormax) but also implicit methods (Beam-Warming scheme) are used. These methods will not be discussed here.

Vorticity and stream function methods (pressure is eliminated from NS equations). These methods are suitable especially for 2D flows.

Primitive variable approach (formulation in terms of velocities and pressure). Pressure represents a continuity equation constraint. Only these methods (pressure corrections methods SIMPLE, SIMPLER, SIMPLEC and PISO) using staggered and collocated grid will be discussed in this lecture.

FVM Navier Stokes equationsCFD7r

2

2

( ) ( )

( ) ( )

0

u

v

u uv p u uS

x y x x x y y

uv v p v vS

x y y x x y y

u v

x y

Example: Steady state and 2D (velocities u,v and pressure p are unknown functions)

This is equation

for u

This is equation

for v

This is equation for

p? But pressure p is not in the continuity equation! Solution of this problem is in SIMPLE methods, described later

FVM checkerboard patternCFD7r

10

10

10

10

10

10

10

10 10

1000

0

0

0

0

0

00 0

0

Other problem is called checkerboard pattern

This nonuniform distribution of pressure has no effect upon NS equations

The problem appears as soon as the u,v,p values are concentrated in the same control volume

W P E

2

( ) ( ) u

vS

x y x x x

u u up

y

u

y

10 10| 0

2 2E W

P

p pp

x x x

FVM staggered gridCFD7r

10

10

10

10

10

10

10

10 10

1000

0

0

0

0

0

00 0

0

Different control volumes for different equations

Control volume for momentum balance

in y-direction

Control volume for momentum balance

in x-direction

Control volume for continuity equation,

temperature, concentrations

2

0 10

( ) ( ) u

p

u u u

x x

vS

x y x x x

p u

y y

FVM staggered gridCFD7r

Location of velocities and pressure in staggered grid

I-2 I-1 I I+1 I+2i-1 i i+1 i+2

J+1

J

J-1

J-2

j+2

j+1

j

j-1

pu

v

Pressure I,J

U i,J

V I,j

Control volume for continuity equation,

temperature, concentrations

Control volume for momentum balance

in x-direction

Control volume for momentum

balance in y-direction

FVM momentum xCFD7r

I-2 I-1 I I+1 I+2i-1 i i+1 i+2

J+1

J

J-1

J-2

j+2

j+1

j

j-1

pu

Pressure I,J

U i,J

V I,j

p

Control volume for momentum balance

in x-direction

, , , 1, ,( )i J i J nb nb i J I J I J iJneighbours

a u a u A p p b

2

( ) ( ) u

vS

x y x x x

u u up

y

u

y

FVM SIMPLE step 1 (velocities)CFD7r

, , 1 ,

* * * *, , ,( )

I j nb I J I JI j nb I j I jneighbours

a v va A p p b

Input: Approximation of pressure p*

, 1, ,

* * * *, , ( )

i J nb I J I Ji J nb i J iJneighbours

a u Au a p p b

FVM SIMPLE step 2 (pressure correction)CFD7r

* * *' ' 'p p p u u u v v v

' ' 'IJ IJ nb nb IJneighbours

a p a p b

, 1, ,

* * * *, , ( )

i J nb I J I Ji J nb i J iJneighbours

a u Au a p p b

, , , 1, ,( )i J i J nb nb i J I J I J iJ

neighbours

a u a u A p p b “true”

solution

, 1, ,

' ' ' ', , ( )

i J nb I J I Ji J nb i Jneighbours

a au u A p p

subtract Neglect!!

!

, 1, , , 1 ,

' ' ' ' ' ', , ,( ) ( )

i J I J I J I J I Ji J I j I jd p p d p pu v

Substitute u*+u’ and v*+v’ into continuity equationSolve equations for pressure corrections

FVM SIMPLE step 3 (update p,u,v)CFD7r

*, , ,

*, , , 1, ,

*, , , , 1 ,

'

( ' ' )

( ' ' )

I J I J I J

i J i J i J I J I J

I j I j I j I J I J

p p p

u u d p p

v v d p p

Continue with the improved pressures to the step 1

Only these new pressures are necessary for next iteration

FVM SIMPLER (SIMPLE Revised)CFD7r

Idea: calculate pressure distribution directly from the Poisson’s equation

2

2

2

( )

( )

( , )

uu p u

uu p u

p f u v

This technique is called Rhie-Chow interpolation.

FVM Rhie Chow (Fluent)CFD7r

Rhie C.M., Chow W.L.: Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA Journal, Vol.21, No.11, 1983

The way how to avoid staggering (difficult implementation in unstructured meshes) was suggested in the paper

FVM Rhie Chow (Fluent)CFD7r

Rhie C.M., Chow W.L.: Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA Journal, Vol.21, No.11, 1983

0

1

1

y

v

x

u

fy

p

y

vv

x

vu

fx

p

y

uv

x

uu

y

x

ynvnPvPvPvP

xnunPuPuPuP

fvASy

pBSv

fuASx

pBSu

*

**

**

*

|

|P Ee

)||(**

**EPuEPEuPu x

p

x

pBuuSS

**

**

* ** * * * * * * ***

|

|

1| ( ( ))

2 2 2 2 2

P Pu u P

E Eu u E

Pu Eu E WE P P E E P EE Pe eu u e u u

pu S B

x

pu S B

x

S S p pp p u u p p p ppu S B B B

x x x x x

How to interpolate ue velocity at a control volume face

Fluent - solversCFD7r

, , cp,…

u…

v…

w…

p… (SIMPLE…)

T,k,,……

Segregated solver

u

u uu

v vut

w wu

E Eu

, , cp,…

k,,……

Coupled solver

Nonstationary problem always solved by implicit Euler

Nonstationary problem solved by implicit Euler (CFL~5), explicit Euler (CFL~1) or Runge Kutta

Algebraic eq. by multigrid AMG Algebraic eq. by multigrid AMG, or FAS (Full Aproximation Storage)

CFL Courant Friedrichs Lewy criterion

Cx

tu

FVM solversCFD7r

Dalí

FVM solversCFD7r

TDMA – Gauss elimination tridiagonal matrix

(obvious choice for 1D problems, but suitable for 2D and 3D problems too, iteratively along the x,y,z grid lines – method of alternating directions ADI)

PDMA – pentadiagonal matrix (suitable for QUICK), Fortran version

CGM – Conjugated Gradient Method (iterative method: each

iteration calculates increment of vector in the direction of gradient of minimised function – square of residual vector)

Multigrid method (Fluent)

Solution of linear algebraic equation system Ax=b

Solver TDMA tridiagonal system MATLABCFD7r

function x = TDMAsolver(a,b,c,d)%a, b, c, and d are the column vectors for the compressed tridiagonal matrixn = length(b); % n is the number of rows % Modify the first-row coefficientsc(1) = c(1) / b(1); % Division by zero risk.d(1) = d(1) / b(1); % Division by zero would imply a singular matrix. for i = 2:n id = 1 / (b(i) - c(i-1) * a(i)); % Division by zero risk. c(i) = c(i)* id; % Last value calculated is redundant. d(i) = (d(i) - d(i-1) * a(i)) * id;end % Now back substitute.x(n) = d(n);for i = n-1:-1:1 x(i) = d(i) - c(i) * x(i + 1);endend

Forward step (elimination entries bellow diagonal)

1 1 1, 2,...,i i i i i i ia x b x c x d i n

Solver TDMA applied to 2D problems (ADI)CFD7r

x

y

Y-implicit step (repeated application of TDMA for 10 equations, proceeds from

left to right) X-implicit step (proceeds bottom up)

Solver MULTIGRIDCFD7r

Rough grid (small wave numbers) Fine grid (large wave numbers details)

Interpolation from rough to fine grid

Averaging from fine to rough grid

Solution on a rough grid takes into account very quickly long waves (distant boundaries etc), that is refined on a finer grid.

Unsteady flowsCFD7r

Discretisation like in the finite difference methods discussed previously

n

n+1

∆x

∆tnPT

nWT

nET

1nPT

n-1

n

n+1

∆x

∆tnPT

1nWT

1nET

1nPT

n-1

n

n+1

∆x

∆tnPT

nWT

nET

1nPT

n-1

1nWT

1nET

Example: Temperature field

EXPLICIT scheme IMPLICIT scheme CRANK NICHOLSON scheme

2

2

xt

a

2

2

T Ta

t x

First order accuracy in time First order accuracy in time Second order accuracy in time

Stable only if Unconditionally stable and bounded Unconditionally stable, but bounded 2x

ta

only if

FVM Boundary conditionsCFD7r

Boundary conditions are specified at faces of cells (not at grid points)

INLET specified u,v,w,T,k, (but not pressure!)

OUTLET nothing is specified (p is calculated from continuity eq.)

PRESSURE only pressure is specified (not velocities)

WALL zero velocities, kP, P calculated from the law of wallPyP

Recommended combinations for several outlets

INLET

PRESSURE

PRESSURE

INLET

OUTLET

OUTLET

PRESSURE

PRESSURE

PRESSURE

Forbidden combinations for several outlets

PRESSURE

PRESSURE

OUTLET

there is no unique solution because no flowrate is specified

Stream function-vorticityCFD7r

Ghirlandaio

Introducing stream function and vorticity instead of primitive variables (velocities and pressure) eliminates the problem with the checkerboard pattern. Continuity equation is exactly satisfied.

However the technique is used mostly for 2D problems (it is sufficient to solve only 2 equations), because in 3D problems 3 components of vorticity and 3 components of vector potential must be solved (comparing with only 4 equations when using primitive variable formulation of NS equations).

Stream function-vorticityCFD7r

The best method for solution of 2D problems of incompressible flows

vy

p

y

vv

x

vu

ux

p

y

uv

x

uu

2

2

1

1

Pressure elimination

2( ) ( ) ( )u v u v u v

u vx y x y y x y x

Introduction vorticity (z-component of vorticity vector)

,u vy x

wvuzyx

kji

u

2u vx y

Velocities expressed in terms of scalar stream function

Vorticity expressed in terms of stream function

2

0

y

v

x

u

y

u

x

v

Stream function-vorticityCFD7r

Three equations for u,v,p are replaced by two equations for , . Advantage: continuity equation is exactly satisfied. Vorticity transport is parabolic equation, Poisson equaion for stream function is eliptic.

| , |w w w wu vy x

2u vx y

How to do it, is shown in the next slide on example of flow in a channel

2

Lines of constant are streamlines (at wall =const). The no-slip condition at wall (prescribed velocities) must be reflected by boundary condition for vorticity

Stream function-vorticityCFD7r

( | , | )w w w wu vy x

x

y

=0

H

w dyyu0

)(

H

Increases from 0 to w for example =uy

1 2 N-1 N

1

2

M

N-1=N

y

u

x

v

Stream function

Vorticity

Zero vorticity at axis (follows from definition)

Prescribed velocity profile at inlet u1(y) and v1(y)=0

y

u

x

12

2

1

221

2

12

21

21

12 2

1

2

1x

xx

xx

x

2 1 11 2

2( ) u

x y

Vorticity at wall (prescribed velocity U, in this case zero)

)(2

12yU

y MMMM

Example: cavity flow 1/4CFD7r

x

y

H

1 2 N-1 N

1

2

M

UM

Cavity flow. Lid of cavity is moving with a constant velocity U. There is no in-out flow, therefore stream function (proportional to flowrate) is zero at boundary, at all walls of cavity

)(2

12yU

y MMMM

Vorticity at moving lid

Vorticity at fixed walls 1 1 2 1 1 1 22 2 2

2 2 2( ) ( ) ( )N N Nx x y

CFD7r

FVM applied to vorticity transport equation (upwind of the first order). Equidistant grid is assumed =x=y

x

y

W P E

UM

N

S

))0,max(

())0,max(

())0,max(

())0,max(

()4||||

(22222

P

NP

SP

EP

WPP

P

vvuuvu

Eliptic equation for stream function

2 2

4 1( )P W E S N P

This systém of algebraic equation (including boundary conditions) is to be solved iteratively (for example by alternating direction method ADI described previously)

2

SNPu

2W E

Pv

Example: cavity flow 2/4

CFD7r Example: cavity flow 3/4Matlab% tok v dutině. metoda vířivosti a pokutové funkice. Upwind. Metoda% střídavých směrů. h=0.01; um=0.0001; visc=1e-6;% parametry site, a casovy krok relax=0.5; n=21; d=h/(n-1); d2=d^2; niter=50; psi=zeros(n,n); omega=zeros(n,n); u=zeros(n,n); v=zeros(n,n); a=zeros(n,1); b=ones(n,1); c=zeros(n,1); r=zeros(n,1); for iter=1:niter for i=1:n u(i,n)=um; end for i=2:n-1 for j=2:n-1 u(i,j)=(psi(i,j+1)-psi(i,j-1))/(2*d); v(i,j)=(psi(i-1,j)-psi(i+1,j))/(2*d); end end

%stream function x-implicitfor j=2:n-1 for i=2:n-1 a(i)=-1/d2;b(i)=4/d2;c(i)=-1/d2; r(i)=(psi(i,j-1)+psi(i,j+1))/d2+ omega(i,j); end r(1)=0;r(n)=0; ps=tridag(a,b,c,r,n); for i=2:n-1 psi(i,j)=(1-relax)*psi(i,j)+relax*ps(i); endend%stream function y-implicitfor i=2:n-1 for j=2:n-1 a(j)=-1/d2;b(j)=4/d2;c(j)=-1/d2;r(j)=(psi(i-1,j)+psi(i+1,j))/d2+ omega(i,j); end r(1)=0;r(n)=0; ps=tridag(a,b,c,r,n); for j=2:n-1 psi(i,j)=(1-relax)*psi(i,j)+relax*ps(j); endend% vorticity boundary conditionsfor i=1:n omega(i,n)=relax*omega(i,n)+(1-relax)*2/d2* (psi(i,n)-psi(i,n-1)-um*d); omega(i,1)=relax*omega(i,1)+(1-relax)*2* (psi(i,1)-psi(i,2))/d2; omega(1,i)=relax*omega(1,i)+(1-relax)*2*(psi(1,i)-psi(2,i))/d2; omega(n,i)=relax*omega(n,i)+(1-relax)*2*(psi(n,i)-psi(n-1,i))/d2;end

%vorticity x-implicitfor j=2:n-1 for i=2:n-1 up=u(i,j);vp=v(i,j); a(i)=-visc/d2-max(up,0)/d; b(i)=4*visc/d2+(abs(up)+abs(vp))/d; c(i)=-visc/d2-max(-up,0)/d; r(i)=omega(i,j-1)*(visc/d2+max(vp,0)/d)+ omega(i,j+1)*(visc/d2+max(-vp,0)/d); end r(1)=omega(1,j);r(n)=omega(n,j); ps=tridag(a,b,c,r,n); for i=2:n-1 omega(i,j)=(1-relax)*omega(i,j)+relax*ps(i); endend%vorticity y-implicitfor i=2:n-1 for j=2:n-1 up=u(i,j);vp=v(i,j); a(j)=-visc/d2-max(vp,0)/d; b(i)=4*visc/d2+(abs(up)+abs(vp))/d; c(j)=-visc/d2-max(-vp,0)/d; r(j)=omega(i-,j)*(visc/d2+max(up,0)/d)+ omega(i+1,j)* (visc/d2+max(-up,0)/d); end r(1)=omega(i,1);r(n)=omega(i,n); ps=tridag(a,b,c,r,n); for j=2:n-1 omega(i,j)=(1-relax)*omega(i,j)+relax*ps(j); endend end

CFD7r Example: cavity flow 4/4

5 10 15 20 25 30

5

10

15

20

25

30

-9

-8

-7

-6

-5

-4

-3

-2

-1

0x 10

-8

5 10 15 20 25 30

5

10

15

20

25

30

-1

0

1

2

3

4

5

6

7

8

9x 10

-5

u

5 10 15 20 25 30

5

10

15

20

25

30

-0.4

-0.35

-0.3

-0.25

-0.2

-0.15

-0.1

-0.05

0

0.05

0.1

5 10 15 20 25 30

5

10

15

20

25

30

-3

-2

-1

0

1

2

3x 10

-5

v

Re=1

10 20 30 40 50 60 70

10

20

30

40

50

60

70

-10

-9

-8

-7

-6

-5

-4

-3

-2

-1

0x 10

-6

Re=1000

CFD7r Example: cavity flow FLUENT

Re=1

u v