Doru Caraeni CD-adapco, USA CFD Futures Conference, August 6-8, - - PowerPoint PPT Presentation

Doru Caraeni

CD-adapco, USA

CFD Futures Conference, August 6-8, 2012

Why I did Residual-based schemes research ?

• (1996) Leading the CFD/CAE group (Centrifugal Compressors) at

COMOTI Bucharest

• Challenge: to perform LES of turbulence inside high-PR CC
• Write a new CFD code (together with Aerospace Department at

“Polytechnica” Institute Bucharest) for industrial LES

• Found a few papers about early Residual-Distribution schemes
• Learned more about these scheme at a VKI Advanced CFD course
• Went to Lund Institute to learn LES of turbulence and develop a

(hopefully) best-in-class LES algorithm, for industry (Dec.1997)

Multidimen idimensional ional Upwind d Residual dual Distribut bution ion scheme me :

• Fluctuation-Splitting (RDS) proposed in 1986 by professor Phil Roe
• Developed by professors and students at Michigan University (Roe), VKI

(Deconinck), Bordeaux University (Abgrall), Polytechnica di Bari, Lund (Caraeni), Univ. of Leeds (Hubbard) etc.

• Compact matrix distribution schemes for steady Euler and Navier-Stokes

equations (E.van der Weide, H. Paillere), 1996.

• Second order RD scheme for LES of turbulence using a residual-

property preserving, dual time-step approach (Caraeni, 1999).

A short early history of MU-RDS

Multidimen idimensional ional Upwind d Residual dual Distribut bution ion scheme: me:

• Second order space-time RD scheme for unsteady simulations (2000, VKI)

(using space-time integration/residual-distribution to achieve accuracy)

• Third order RD scheme for steady inviscid flow simulations (2000, LTH)

(node gradient-recovery for quadratic solution representation)

• Third order RD scheme for the unsteady turbulent flow simulations (2001,

LTH). (node gradient-recovery and residual-property satisfying)

• Third order results with above gradient-recovery idea reported by Rad and

Nishikawa (2002, MU).

• High-order (>3) RD scheme for scalar transport equations (2002, BU & MU).

(sub-mesh reconstruction for high-order solution representation) “ Third-order non-oscillatory fluctuation schemes for steady scalar

conservation laws ” M. Hubbard, 2008.

A short early history of MU-RDS (cont.)

Ricchiuto, CEMRACS, 2012

A 3rd Order Residual-Distribution scheme for Navier-Stokes simulations

(Residual-property satisfying formulation)

(A Third Order Residual-Distribution Method for Steady/Unsteady Simulations: Formulation and Benchmarking, including LES, Caraeni, VKI, 2005)

             

j j j ij i j j t j ij i j j i t i j j t

T u h u e p u u u u

, , , , , , , , , , ,

) ( ) ( ) ( ) ( ) ( ) ( ) (         

V j j C j j t

F F U

, , ,

 

             h u p u u p u u p u u u F

j j j j j j j j C j

       

3 3 2 2 1 1

            

j j j j j j j V j

T u u u F

, 3 3 2 2 1 1 3 2 1

      

t

e u u u U ) , , , , (

3 2 1

     

t v j j c j j

U F F U

, , , ,

   



 

 

    dv U F F dv U

t j v j c j

] ) [(

, , ,

Jameson dual-time algorithm

High-order RD scheme

for Navier-Stokes equations

dv

T c j c T

F .

,



 

Convection flux cell-residual

k n T i T uns T v T c T T i k n i k n i

B U U

i

V

, 1 , , 1 1 , 1

) (

     



       

Update scheme for steady/unsteady simulations (Caraeni):

T i

B Upwind matrix residual Distribution coefficient (bounded)

dv

T v j v T

F .

,



 

Diffusion flux cell-residual

dv U

T t uns T

.

,



 

Unsteady term cell-residual







T i T i

I B

High-order RD scheme (cont.)

for Navier-Stokes equations

(conservativity)

] ) ( 2 1 4 1 [

1 , 



 

j j i LW T i

K K I

B



Distribution schemes (for preconditioned system):

Low Diffusion A (LDA) Lax-Wendroff (LW)

1 ,

) (

   



j j i LDA T i

K K B

1 1 1 , ,

3 1 3 1 3 1 3 1

      

    

i i i i i i i i i i i i j U j i

R R K R R K R R n F K   



 

4 1 i i c T

U K

Computes the convective cell residual with second

• rder accuracy (linear data)

High-order RD scheme (cont.)

for Navier-Stokes equations

High-order RD scheme (cont.)

How to construct a 3rd-order RDS (Ph.D. 2000, LTH):

• 1. Use (upwind or upwind-biased) uniformly bounded residual-distribution

coefficients (linearity/accuracy preserving RD scheme), and apply to total cell- residual

• 2. Compute the total cell residual (convective + diffusive + unsteady terms)

with the required accuracy:

• we used condition-1 + linear solution, second order accurate integration for

2nd order RDS

• we need to use condition-1 + use quadratic reconstruction, 3rd order

accurate integration for 3rd order RDS The idea is to use the same accuracy-preserving RD scheme, as for second

• rder schemes, but compute the total cell residual with 3rd order accuracy.

) ( 8 ) ( 2

1 1 1

, , ) ( i i i j i j i i mid

r r

Z Z Z Z Z

    

 



  

T c c T j c T

dS F dv F . .

,

                     

j j j j j j j j c

Z Z p Z Z p Z Z p Z Z Z Z F

4 3 3 2 2 1 1

  

Use parameter variable Z and assume a quadratic variation

• ver the tetrahedral cell.

) 2 ( 1

2 3 2 2 2 1 4

Z Z Z

Z Z p       

) , , , , 1 (

3 2 1

H u u u Z  

2

2 i p

u T c H  

High-order RD scheme (cont.)

Convection residual discretization, 3rd order.

Cell-residual in integral form:

Z,j computed with 2nd order accuracy

(multi-step algorithm)

High-order RD scheme (cont.)

Convection residual discretization, 3rd order.

  

 

  

4 1 , ,

} ) ( { .

i face i j i j c j T c c T

i

d n n F dS F 

                                 

     

i i i i i i

face j i j face j j i j face j j i j face j j i j face j i j face i j c j

d Z Z n d p Z Z n d p Z Z n d p Z Z n d Z Z n d n F          ) ( ) ( ) ( ) ( ) ( ) (

4 , 3 3 , 2 2 , 1 1 , , ,





i i k J

face j k face Z Z

d Z Z

I

 ) (

  





6 1 , ) ( ) (

1 1 1

} { ) (

i i face i i i j i k face j k

i i

d H H d Z Z

Z Z

 

k j HH face face k j

I A d H H

i i

, ) ( ) (

] [ 





k j HH

I

,

] [

Pre-computed matrix

ds F dv F

T v T v j v T

.

,

 



  

k k face v T v v T

n F ds F

k .

.

4 1 ) (

 

 

  

  

  

] ) ( ) [(

, , , i i i

u u u  

Assuming a quadratic variation of the Z variables over the cell, the

diffusive flux vector integral can be computed over the cell-face. Use the values of the Z variable and its gradients, defined in the nodes of the high-order FEM tetrahedral-cell.

High-order RD scheme (cont.)

Diffusion residual discretization, 3rd order.

High-order RD scheme (cont.)

Unsteady residual discretization , 3rd Order.

t U U U U

n n k n t

   

 

. 2 4 3

1 , 1 ,

  



  

T t T t v T

dv U Q dv U

9 ) ( , ) ( ,

. .

  

dv Q U

T t v T

.

) ( 9 ) ( ,

 



 

  

9 ,.., 4 ; 20 . 3 ,.., ; 05 . ; .

) (

         

    T T T

V I V I dv Q I

2nd order discretization in time, and 3rd order in space:

High-order RD scheme (cont.)

Monotone shock capturing

• 1. Shock detection or or
• 2. Blending between the high-order scheme and

a first order positive RD scheme (the N-scheme)

,..) , ( : . ). 1 (     P f where

N i LDA i T i

      



= 0 for a smooth flow = 1 (discontinuity detected)

) (

2 2 av

p p h    

 

  

 nodes N j # 1

) . (

2 2

V V      

• Uses a Multi-D Upwind Residual-Distribution scheme
• Formulated for fully unstructured grids (tetrahedrons),
• Compact scheme, highly efficient parallel algorithm.
• Implicit time integration (dual time-stepping algorithm).
• 3rd - order accuracy in space (using FEM integration)
• 2nd - order time discretization (BDF2 scheme)
• Acceleration techniques: preconditioning, point-

implicit relaxation, geometric multi-grid, etc.

Summary of this 3rd order RD algorithm

• Steady inviscid flows
• a. Sine-bump channel flow inlet Mach 0.5.
• Steady viscous flows
• c. Laminar flat-plate boundary layer, Reynolds 2000.
• Unsteady inviscid flows
• b. Vortex transport by uniform flow Mach 0.04.
• Shock capturing
• d. Shock vortex interaction.
• Large Eddy Simulation
• e. LES of turbulent channel flow.

Results

Inviscid sine bump channel flow:

• Inlet Mach number 0.5

Steady Euler

LDA 2nd LDA 3rd

Asymptotic accuracy Steady Euler 0.0001 0.001 0.01 0.1 0.1 1 h Error Classical Num3rd

RDS Solution on 32x8 grid

Maximum entropy production:

• 2nd order scheme E=2.3 e-4
• 3rd order scheme E=5.1 e-6
• Cell centered FVM(2ndO) E=4.2 e-4

Laminar viscous flow over a flat plate:

• Infinite Mach # 0.5, Re 2000
• Grid0 of 32x18 grid points

Laminar flat plate boundary layer Re=2000

2 4 6 8 10 12 14 16 18 0.5 1 1.5

U/Uinf y*sqrt(Uinf/miu/x)

Blasius-U/Uinf Num.3rdO Lev0 Num.3rdO Lev1 Num.3rdO Lev2 Num.3rdO Lev3

Laminar flat plate boundary layer Re=2000

2 4 6 8 10 12 14 16 18 0.01 0.02 0.03

V/Uinf y*sqrt(Uinf/miu/x)

Blasius-V/Uinf Num.3rdO Lev0 Num.3rdO Lev1 Num.3rdO Lev2 Num.3rdO Lev3

Steady Navier-Stokes

Asymptotic accuracy Laminar flat plate 0.0001 0.001 0.01 0.1 0.1 1

h Error Num3rd

Unsteady Euler

Vortex transport by Inviscid flow.

• Uniform flow Mach = 0.04

Third order results, grid 64x64

Vortex transport by Inviscid flow.

• Uniform flow Mach = 0.04

Unsteady Euler (cont.)

Vorticity X=0.05, T=5

• 10
• 5

5 10 15 20 25 30 0.025 0.05 0.075 0.1 Y OmegaZ

3rdO Sol. (32x32x2)

• Sol. (32x32x2) T=0

3rdO Sol. (64x64x4)

• Sol. (64x64x4) T=0

2ndO Sol. (32x32x2) 2ndO Sol. (64x64x4)

Vorticity Y=0.05 , T=5

• 10
• 5

5 10 15 20 25 30 0.025 0.05 0.075 0.1 X

OmegaZ

3rdO Sol. (32x32x2) Sol.at T=0 (32x32x2) 3rdO Sol. (64x64x4)

• Sol. (64x64x4) T=0

2ndO Sol.(32x32x2) 2ndO Sol. (64x64x4)

Asymptotic accuracy Vortical Euler flow

0.00001 0.0001 0.001 0.01 0.1 0.1 1

h Error Num3rd

Unsteady Euler (cont.)

Vortex transport by Inviscid flow.

• Uniform flow Mach = 0.04

Third order results, grid 64x64

T= 0 T = 12 periods T = 24 periods

Shock vortex interaction

Shock-vortex interaction:

• Steady shock in mid channel
• Vortex moves from left to right

Note: vortex preserving strength, before and after crossing shock

Shock vortex interaction

LES of turbulent channel flow

Turbulent channel flow:

• Reynolds # = 5400
• Ret= 344

0.2 0.4 0.6 0.8 1 1.2 1.4 0.25 0.5 0.75 1 y/h

U/ut

DNS(Kim et al.) 2ndO LES 3rdO LES 0.5 1 1.5 2 2.5 3 3.5 4 0.1 0.2 0.3 0.4 0.5

y/h u'rms/ut

• Exp. (Kreplin)

DNS (Kim et al.) 2ndO LES 3rdO LES 0.2 0.4 0.6 0.8 1 1.2 0.1 0.2 0.3 0.4 0.5

y/h v'rms/ut

Exp. DNS(Kim et al.) 2ndO LES 3rdO LES 0.2 0.4 0.6 0.8 1 1.2 1.4 0.1 0.2 0.3 0.4 0.5

y/h w'rms/ut

• Exp. (Kreplin)

DNS(Kim et al.) 2ndO LES 3rdO LES

LES Smagorinsky ut/Ub Uc/Ub 3rd O MDU 0.0640 | 1.164 2nd O MDU 0.0627 | 1.186

• DNS (Kim et al.) 0.0643 | 1.162

Multidimensional Residual-Distribution Solving for flow and “optimal” mesh

(Grids and solutions from Residual Minimization, Nishikawa, Rad, Roe, 2001)

Solving for flow and solution using RDS

Main ideas:

• Use multidimensional RDS to compute solution at vertices,
• There are 5-6 times less vertices than cells in the tetrahedral-

cells mesh …

• Use the extra “conditions” (cell-residual must be driven to

zero) to define mesh motion equations, using an LSQ approach,

• Algorithm computes an improved solution on a “optimized”

mesh, which minimizes the overall error in a specific norm.

(Nishikawa, 2001)

Solving for flow and solution using RDS

Flow over Joukowsky airfoil (known theoretical solution)

Original mesh Adapted mesh (Nishikawa, 2001)

Solving for flow and solution using RDS

Flow over Joukowsky airfoil (known theoretical solution)

Original mesh solution Cp, o Adapted mesh Cp, o (Nishikawa, 2001) Comparison with theoretical solution ----

• Resolves better real complex multidimensional physics (!)
• It is much more accurate that 2nd order Finite Volume method,
• It is capable of handling complex geometry (formulated tetrahedrons),
• Has a compact stencil algorithm, at every step (which leads to very

efficient parallelization),

• It is relatively to easy to extend to high order accuracy (at least from

2nd to 3rd order), and 3rd order results are significantly more accurate,

• Can be used to solve for flow and node location - using the

combined RDS/LSQ approach - for an optimal solution, on a given mesh topology.

Why using Multi-D Residual-Distribution schemes ?

Backup slides

High-order Residual Distribution Scheme for Scalar Transport Equations

• n Triangular Meshes

From “High-order fluctuations schemes on triangular meshes” R.Abgrall and Phil Roe, 2002

High-order RD scheme for scalar equations

) ( 3 2 2

III II I IV H

        

High-order RD scheme for scalar equations

High-order RD scheme for scalar equations

High-order RD scheme for scalar equations

From (Hubbard, J. Computational Physics 2007)

“Status of Multidimensional Residual Distribution Schemes and Applications in Aeronautics”, Deconinck et al. AIAA 2000-2328.

Space-time Residual Distribution schemes for unsteady simulations

“Construction of 2nd order monotone and stable residual distribution schemes: the unsteady case”, Abgrall et al. VKI 2002

Space-time RD for unsteady simulations

t t t x u t t t x u t x u

n n n n h

     

 

) ( ) ( ) , (

1 1







      

1

. ) (

.

n n

t t i i h t x

dt dx x u A t u

 

     

            

i n i n i i s n i n i t i i s t i t x

u u K t u u ) ( 2 ) ( 3

1 1 ,

) ( 2 ) ( 2 ) ( 3

1 1 1 1 , , n n i i n n i i n i n i N n i N n i

u u K t u u K t u u V           

      

   

           

             

i n j j i i n i n j j i i n

u K K u u K K u

1 1 1 1

“Upwind-in-time”

Space-time RD for unsteady simulations

s i i j t j i t i i LDA n i LDA n i

          



 

   ) 12 1 ( ) 6 1 (

1 , ,

1

) (

   



j j i i

K K



LDA space-time scheme: LDA+N space-time scheme:



   

          

j N n i t x LDA n i N n i B n i B n i

l l l

1 , , 1 , 1 , 1 , ,

) 1 (

“Status of Multidimensional Residual Distribution Schemes and Applications in Aeronautics” Deconinck et al. 2000

Space-time RD for unsteady simulations

Reflection of a planar shock from a ramp (density plot) Shock reflection on a forward facing step (density plot) Shock-vortex interaction From “Construction of 2nd order monotone and stable residual distribution schemes: the unsteady case”, Abgrall et al. 2002

Space-time RD for unsteady simulations

Convection of vortex

• Periodic BC’s
• One revolution simulated
• 2nd and 3rd order ST-RDS

compared

• Pressure contours displayed

From (Nadege Villedieu , VKI Ph.D., 2009)