Implicit Smoothing with Runge-Kutta Schemes for Navier-Stokes CFD - - PowerPoint PPT Presentation

Rapid Convergence using Implicit Smoothing with Runge-Kutta Schemes for Navier-Stokes

CFD Futures Conference NIA, Hampton, VA August 6-8, 2012

Eli Turkel

Introduction

• Efficiency is a challenge with increasing
• Problem Size
• Flow Complexities
• Improvements in computational efficiency
• Multigrid
• Local Time Stepping
• Implicit Residual Smoothing

 Explicit time marching Runge-Kutta scheme  Space centered or upwind schemes

RK/implicit smoother scheme

RK/Implicit smoother scheme

finite volume + source terms

Extension to turbulent flow Extension to chemical reactions Time dependent - dual time step

Navier-Stokes Equations

 

T n

e w v u Q        ...

2 1



N-S equations in conservative form:

S z H H y G G x F F t Q

V V V

               ) ( ) ( ) (

• F , G and H are the inviscid fluxes
• Fv , Gv and Hv are the viscous fluxes

The Navier-Stokes (N-S) equations for the conservative variables:

k-ω SST Turbulent model

• In turbulent flow, the viscosity is the sum of the laminar

and the turbulent viscosity

• In addition to N-S equation we solve the k-ω SST

turbulence model equations -      

j j j T j j T k j k

x x k F x x S u t x k x k P k u t k                                                                

 

1 1 2

2 1 2 2 *

 

( )

1 1 2

max ,

t T

a k a SF ρ µ ρν ω = =

The turbulent viscosity μt is given by

• k and ω are the model variables,
• S is the magnitude of the vorticity,
• F1,2 are dump functions,
• ν, νT are the laminar and turbulent viscosity parameters

10

Real gas equation of state





 

n i i

RT W RT P

1

 



   n i i i

w W

1 1

1  

The equation of state is with the mean molecular weight

Chemical reactions source terms

E C C B A

f r f r

k k k k

, 2 , 2 , 1 , 1

2   

The source term vector, S, describes the rate of change of species k:

N 1,2,..., k dt d

reactions all i k

   



 i ik k k k k

q w        

where ωk and wk are the molar concentration and molecular weight of the species. νik are the stochiometric coefficients of the species k in the reaction i. qi are the rate of progress variables given by -

i r i f i

q q q

, , 







spieces k k i f i f

ki

k q

' , , 







spieces k k i r i r

ki

k q

' ' , , 



qf,i and qr,i are defined by:

;

The forward reaction rate kf and the reverse rate kr are empirically known functions of the temperature. The forward constant is given by an Arrhenius expression of the type:

RT E i i f

i ie

T A k

/ , −

=

β

• Ai, βi and Ei are the Arrhenius constants:
• Ai is the rate constant
• βi is the temperature exponent
• Ei is the activation energy.

The source term for the temperature is

∑

− =

k k k v

h c T ρ ρ   1 ∑       =

∈ i K k ik

RT P K k

atm p i r υ ,

      ∆ − ∆ = RT H R S K p exp

RK/Standard Scheme: Three Components

• The 3 components are: RK scheme, implicit residual smoothing, multigrid
• The qth stage of the RK component can be written as

W(q) = W(0) − α1∆t R(W(q−1))

where R is the vector residual function, ∆t is the time step, and the RK coefficients αq are [0.25, 0.1667, 0.375, 0.5, 1.0]. The residual function R is given by

R = R(W(q)) = 1

V

 LcW(q) +

q

• r=0

γqr LvW(r) +

q

• r=0

γqr LdW(r)

  ,

with the operators Lc, Lv, and Ld for convection, viscous diffusion and numerical dissipation.

RK/Implicit Scheme

• The change in the solution on the qth stage

δW(q) = W(q) − W(0) = −αq ∆t V L W(q−1) =

R(W(q−1)),

L is the complete difference operator.

• If we apply implicit residual smoothing, then

Li δW = δW(q), where Li is an implicit operator.

• Approximately inverting the implicit operator Li

δW = −αq ∆t V P L W(q−1), P is a preconditioner defined by the approximate inverse L−1

i

.

Transforming the equations to primitive variables, the flux Jacobian is written as:

  

 A A A

where

 

A A A  



2 1

local all faces NB all faces

t S I A ndS t Q V Q t Q A Q ndS V  

 

                   

 

Finally, the implicit smoothing scheme is given by:

Re Im

• 10
• 8
• 6
• 4
• 2
• 3
• 2
• 1

1 2 3

wgts: [1 1 1] wgts: [1 .5 .5]

(a)

Re Im

• 10
• 8
• 6
• 4
• 2
• 3
• 2
• 1

1 2 3

wgts: [1 1 1] wgts: [1 .5 .5]

(b) Figure 1: Fourier footprints of RK(3,3) scheme with two preconditioners for all modes with high-frequency components (64 × 64, M = 0.5, α = 0◦, CFL = 103, AR = 5, Re = ∞). (a) Without entropy fix, (b) with entropy fix.

0.35 . 2 5 0.2 0.1 0.35 0.3 0.4 0.15 0.15 0.1 0.1 0.35 0.25 0.1 0.2 0.2 0.35 0.2 0.4 0.6 0.8 0.8

θx θy

• 3
• 2
• 1

1 2 3

• 3
• 2
• 1

1 2 3

(a)

0.2 0.2 0.2 0.25 0.25 0.2 0.15 0.15 0.2 0.6 0.8 0.8 0.4

θx θy

• 3
• 2
• 1

1 2 3

• 3
• 2
• 1

1 2 3

(b) Figure 2: Effect of dissipation weights on damping behavior of RKI(3, 3) scheme (64×64, M = 0.5, α = 45◦, CFL = 103, AR = 1, Re = ∞, 2 SGS). (a) wgts: [1, 1, 1], (b) wgts: [1, 0.5, 0.5].

0.25 . 2 0.15 0.25 . 2 0.25 . 9 0.15 . 4

θx θy

• 3
• 2
• 1

1 2 3

• 3
• 2
• 1

1 2 3

(a)

0.2 0.2 0.2 0.2 0.15 0.15 . 2 5 . 3 5 0.9

θx θy

• 3
• 2
• 1

1 2 3

• 3
• 2
• 1

1 2 3

(b) Figure 3: Damping behavior of RKI(3, 3) scheme with variation in Re and AR (64 × 64, M = 0.5, α = 45◦, 2 SGS). (a) Re = 102, AR = 10, CFL = 103, (b) Re = 106, AR = 103, CFL = 104.

Extension to turbulent flow

             ~ ~ ~ k Q

                2

* *

k Q S

We solve two sets of equations for the smoothed residuals independently: a) For the Eulerian primitive variables b) For the turbulent model variables – k and ω

In the k- ω set we use the Jacobian of the source term:

Extension to chemical reactions

For a real gas, we insert part of the entries of the Jacobian source terms into the system for the Eulerian variables’ residuals. For the {ω,T, u} primitive variables, the Jacobian matrix we use is

{ }

                        ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = T T T T T u T RJ        , ,

3 3 3 2 2 2 1 1 1

ω ω ω ω ω ω ω ω ω ω

For {ρ, P, u} variables this matrix becomes { }

                          ∂ ∂         ∂ ∂ − ∂ ∂         ∂ ∂ − ∂ ∂         ∂ ∂ − ∂ ∂ ∂ ∂ ∂ ∂ − ∂ ∂ ∂ ∂ − ∂ ∂ − ∂ ∂ ∂ ∂ − ∂ ∂ − ∂ ∂ ∂ ∂ − ∂ ∂ ∂ ∂ − ∂ ∂ − ∂ ∂ − ∂ ∂ = T P P T T P P T w RT T P P T w RT T P P T w RT T R W w TW T TW T w w TW T w w T R W w TW T w w TW T TW T w w T R W w TW T w w TW T w w TW T u P RJ                       , ,

3 3 3 2 2 2 1 1 1 3 3 3 3 3 3 2 3 3 1 3 2 2 2 3 2 2 2 2 2 1 2 1 1 1 3 1 1 2 1 1 1 1

ω ω ω ω ω ω ω ρ ρ ω ω ω ρ ω ρ ω ω ρ ρ ω ρ ω ω ω ρ ω ω ρ ρ ω ρ ω ρ ω ω ω ρ

where

T P T T T RT T P

i

∂ ∂ + ∂ ∂ = ∂ ∂

∑

   ω

Determination of temperature from internal energy

For a given internal energy e0, we want to determine the temperature. We solve the equation f(T) = e(T) – e0 = 0 iteratively using the Newton-Rapson method

   

n n n n

T f T f T T '

1

 



 

v

c dT T de 

     

n v n n n v n

T c e T e T T c T

1

  



Time dependent - Dual time step

1 1

3 4 2

n n n

Q Q Q Q t t

 

     

For steady state calculations we use pseudo time approach. For time dependent calculations, in order to be able to use all of the acceleration methods we use a dual time step. An approximation to the physical time derivative now appears as a source term in the right hand side of the N-S equations.

Q Q F t        

t is the physical time and τ is the pseudo time. We approximate the physical time derivative using a backward difference scheme.

1

3 4 2

k n n k

Q Q Q R F t



     

Since we do not know Qn+1 we approximate it with Qk+1 After some rearrangement of the above equation we have

is the smoothed residual and R is

k k

R

t k k t k k k k

Q R Q Q

    

    

 

   

2 3 2 3 1

1 ~

Preliminary results

Convergence Histories for Three-Stage Schemes

Cycles Log(||Res||2)

50 100 150 200 250

• 14
• 12
• 10
• 8
• 6
• 4
• 2

Case 1 Case 9 Case 10 RAE 2822, RK3/Implicit, Matrix Dissip. 320 x 64 CFL = 1000

RAE 2822: Cases 1, 9, 10

!"#$"%&'($)(*+ !"#$"%&'($)(*+ ,-(%$"%.)& ,-(%$"%.)&/ / 012+- 012+-

3,(%*(-* 3,(%*(-* 456& 456&7 7 8&2+9+2$ 8&2+9+2$ :;& :;&7 7 8&$,(<+$=&)+%,-(2&*.>>+-+%)+$ 8&$,(<+$=&)+%,-(2&*.>>+-+%)+$ '4?@A '4?@A BC1&DEF2.).,&-+$.*1(2&3E"",C.%< BC1&DEF2.).,&-+$.*1(2&3E"",C.%< 'GH&@&IJKJ&$+)"%*$ 'GH&@&IJKJ&$+)"%*$ 5D:3 5D:3 456& 456&7 7 8&2+9+2$ 8&2+9+2$ :;& :;&7 7 A&$,(<+$=&)+%,-(2&*.>>+-+%)+$ A&$,(<+$=&)+%,-(2&*.>>+-+%)+$ '4?@ILLL '4?@ILLL M&363&$N++F$& M&363&$N++F$&!@LKOP 'GH&@&8KP&$+)"%*$ 'GH&@&8KP&$+)"%*$

!

x Y

0.1 0.2 0.3

• 0.05

0.05 0.1 0.15

M 2.1 1.7 1.3 0.9 0.5 0.1

(a)

x Y

0.1 0.2 0.3

• 0.05

0.05 0.1 0.15

M 0.1 0.08 0.06 0.04 0.02

(b) Figure 4: Rocket Motor (a) full range of Mach contours in the chamber and nozzle (b) zoom on the Mach values inside the motor where the sound speed is around 1000 m/s and the Mach values are low.

cycles log(err)

500 1000 10

• 7

10

• 6

10

• 5

10

• 4

10

• 3

10

• 2

10

• 1

10 10

1

implicit smoothing RK/implicit smoothing

Figure 5: Convergence rate for rocket motor for original and improved schemes

!"#$%&#'()*#&+$,&$+"-(Solver

DLR-F4 ! SA Turbulence model Re = 3.000.000 - M = 0.75

!

!"#$%&'()*+,,,+,,,'-.//0'1'23'1'!456/.'78.-9094:'1%;<

! )*'*%69=$784-0 ! >"?'=9@.&'(AB'@9:'C48'*'48D.80'-4:E.8F.:-.'GB,'-H-/.0I+

Inviscid computations on a ONERA M6 with a tetrahedral mesh

Geometry: ONERA M6

Flow conditions:

M = 0.8395

Re = 11,720,000  = 3.06° inviscid

Mesh

Euler mesh Made of 611,856 tetrahedral cells Generated using Gmsh

Mesh

5

Code performances: Lift and Drag

Decay of density residual for solvers with/without convergence acceleration

6

.

Code performances: Lift and Drag

Evolution of lift and drag

7

Code performances: Lift and Drag

Distribution of pressure coefficient on the wing surface, with Mach iso-countours at different spanwise locations

Convergence History for Various CFL Numbers

Log (err) Different CFL # are used for the N-S eqs. and the turb. Eqs.

• N-S CFL = 1x105; Turb. CFL = 2x104
• N-S CFL = 1000; Turb. CFL = 1000
• N-S CFL = 100; Turb. CFL = 100

Laminar reactive flow – Rapid expansion diffuser

Contours maps for the temperature (A), Mach number (B), H2O mass fraction (C) and OH radical (D) A B C D

Convergence history for CFL 20 without multigrid acceleration

Turbulent reactive flow – rocket motor plume

A rocket motor plume exiting from the motor nozzle into a low Mach number free stream flow is calculated. The plume boundary conditions are defined on the nozzle throat where the flow velocity is sonic and the species mass fractions are

• defined. The species used for this problem

are: H, O, OH, H2, O2, CO, CO2, H2O, HCL and N2.

Convergence history: The fluid CFL is 100,000 and turbulent CFL is 200

• Total density
• Turbulent parameter ω

Turbulent kinetic energy Temperature

CO2 mass fraction (A) and OH radical (B) A B

32

34

For CFL=100, we need on the average 5-10 subiterations per physical time step. For CFL=1000 this rises to 6-11 subiterations per physical time step with a total CPU increase of about 20-30%. For CFL=100, we typically need 7-15 subiterations while for CFL=1000 we need 8-17 subiterations. In summary: the preconditioned dual time step code needs

• nly about 10% of the CPU of the original code.

3

10 TOL





4

10 TOL





Conclusions

• centered, upwind schemes
• structured, unstructured grids
• large time step for 2 equation

turbulence model

• faster convergence chemical reactions
• faster convergence in subiterations of

dual time step

• more robust solutions

• Other turbulence models
• Multigrid in turbulence equations
• LES
• Higher Order accurate schemes

(DGS, Spectral Volume/Difference, WENO)

• Non RK algorithms (ADI,GS,Krylov)