Smoother Scheme Oren Peles and Eli Turkel Department of Applied - - PowerPoint PPT Presentation

Improvements of Unsteady Simulations for Compressible Navier Stokes Based on a RK/Implicit Smoother Scheme

Oren Peles and Eli Turkel

In memoriam of Prof. Saul Abarbanel. August 20-24, 2018

Department of Applied Mathematics, Tel-Aviv University

• Background and motivation
• Governing equations
• Two-fluids model
• RK/Implicit smoother with source terms
• Dual time stepping approach with backward differencing
• Dual time stepping approach with DIRK
• Adaptive time steps
• Lagrangian trajectories of burning droplets
• Examples
• Conclusions and future work

Lecture outline

Motivation - Leidenfrost effect of droplets moving on a ratchet -Leidenfrost maze*

layer keeping vapor , produces an insulating boiling point in which a liquid, in near contact with a mass significantly hotter than the liquid's phenomenon is a physical effect Le Leid idenfrost The : cooking is most commonly seen when is

• rapidly. Because of this 'repulsive force', a droplet hovers over the surface rather than making physical contact with it. Th

boiling that liquid from across the pan and takes longer to evaporate skitters point, the water Leidenfrost is at or above the temperature

• ne sprinkles drops of water in a pan to gauge its temperature: if the pan's

than in a pan below the temperature of the Leidenfrost point (but still above boiling temperature).

*Cheng, C., Guy, M., Narduzzo, A., & Takashina, K. (2015). The Leidenfrost Maze. European Journal of Physics, 36(3), 035004

• We wish to simulate the maze.
• Problem statement – very low Mach number, unsteady, two-

phase flow (surrounding air and vapours from the droplets) with moving finite size droplets.

Background

• We wish to accelerate the computation of time dependent problems in the

framework of a dual-time stepping based code with an explicit scheme for pseudo time and an implicit scheme for the physical time.

• Physical models – low Mach compressible flows with perfect gas, two-fluids

mixing model, inviscid and viscous flows.

• Acceleration methods within the pseudo-time marching: multigrid, local

time steps, implicit residual smoothing

• For the pseudo time stepping a low-storage RK/Implicit Smoother method

is used.

• We consider the schemes for the physical time – Backward differencing

(BDF) and Diagonally Implicit Runge-Kutta (DIRK).

• Adaptive time steps based on DIRK scheme.

Governing equations

The Navier-Stokes equations for the conservative variables 𝑅 = 𝜍, 𝜍𝑣, 𝜍𝑤, 𝜍𝑥, 𝜍𝐹 𝑈 in conservation form are:

𝜖𝑅 𝜖𝑢 + 𝛼𝐺 𝑅 = 𝛼𝐺

𝑤 𝑅 + 𝑇

𝐺 is the inviscid flux 𝐺 = 𝜍𝑣, 𝜍𝑣𝑣 + 𝑞𝑦 , 𝜍𝑣𝑣 + 𝑞𝑧 , 𝜍𝑣w + 𝑞𝑨 , 𝜍𝑣 𝐹 + 𝑞 𝑈 𝐺

𝑤 is the viscous flux

𝐺

𝑤 = 0, 𝜐

𝑦, 𝜐 𝑧, 𝜐 𝑨, 𝑣𝜐 𝑦 + v𝜐 𝑧 + 𝑥𝜐 𝑨 + 𝑙𝛼𝑈

𝑈

𝜍 is the gas density, p is the pressure, u, v and w are the velocity vector components and E is the internal energy

𝜍𝐹 = 𝑞 𝛿 − 1 + 1 2 𝜍 𝑣2 + 𝑤2 + 𝑥2

𝜐 𝑦, 𝜐 𝑧, 𝜐 𝑨 are the stress vectors, k is the thermal conductivity and T is the temperature

Two-Fluids model

• Simple modelling of steady and unsteady mixing of two perfect gases

Model governing equations

         

 

 

                                                                                                               

           

                                h h Sc u z u y u x T H u t E z p w u t u w y p v u t u v x p u u t u u Sc u t Sc u t

z y x z y x

                    Pr ˆ ˆ ˆ

• The continuity equation is replaced by two continuity equations for the

two gases composing the mixture

• The modified Navier-Stokes equations are given by

     

   w y w y W    

1

Energy and Total enthalpy

p E H u p E          

2

2 1 1 

v p v p

C C w y w y w C w y w y w C                                

             

1 1 1 1 1 1

Heat capacities Total density and mean molecular weight Equation of state

𝑞 = 𝜍𝑆𝑈 𝑋

             

                                                                 1 1 2 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 / / / 1 / / / 1 / / / 1

2 2

U c w w p U c w w p w v u w w v v u u W U

p b b p a a

                   

Transformation Jacobian for primitive to conservative variables

                                                1 1 2 1 1 1 1 2 1 1 1 1 1 1

2 2

U c w w p U c w w p w v u w w v v u u U W

v b b v a a

          

Conservative to primitive variables transformation Jacobian

RK/Implicit smoother with source terms

 

       

faces all S

ds n F V S s d n F V S t Q   1 1

We solve the system in conservation form: Applying the Gauss theorem for an arbitrary control volume and a discretization in time yields the basic scheme:  

Q S F t Q     

The low-storage Runge-Kutta time marching scheme is:

The residual of the k-th step is:

For accelerating the calculations using large CFL numbers the residual is replaced by a spatially smoothed residual Δ𝑅

𝑅(0) = 𝑅𝑜 𝑅(𝑙+1) = 𝑅(0) − 𝛽𝑙+1Δ𝜐𝑆𝑙 𝑙 = 0 … 𝑞-1 𝑅n+1 = 𝑅(𝑞) 𝛽𝑙 are the RK coefficients

𝑆k = 1 𝑊 𝐺

all faces

⋅ 𝑜 ⅆ𝑡 − 𝑇𝑙

Following Rossow 2006, Swanson et al. 2007, we start with the spatially discretized equation: Linearizing F and S in time we obtain:

• the flux Jacobian

Q S  

• the source Jacobian

1      



S ds n F V t Q

faces all



k faces all

R Q Q S t ds n A V t I                 





Q F A   

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

  

 A A A

where

 

A A A  



2 1

Finally, the implicit smoothing scheme is given by:

𝐽 + 𝜁 Δ𝜐 𝑊 𝐵+ ⋅ 𝑜ⅆ𝑇

𝑏𝑚𝑚 𝑔𝑏𝑑𝑓𝑡

− Δ𝜐 𝜖𝑇 𝜖𝑅 Δ𝑅𝑚𝑝𝑑𝑏𝑚 = 𝑆𝑙 − 𝜁 Δ𝜐 𝑊 𝐵− ⋅ 𝑜 ⋅ Δ𝑅𝑂𝐶ⅆ𝑇

𝑏𝑚𝑚 𝑔𝑏𝑑𝑓𝑡

ε is a relaxation factor. The implicit smoothing scheme is solved iteratively, using symmetric Gauss- Seidel method or red-black iterations (RB is useful for parallel computing).

Dual time stepping approach with backward differencing

The system we solve in the pseudo time 𝜐 is 𝜖𝑅 𝜖𝜐 + 𝜖𝑅 𝜖𝑢 + 𝛼𝐺 𝑅 = 𝛼𝐺

𝑤 𝑅

For each physical time step, we solve the NS equation for a steady state . The derivative with respect to the physical time is a source term and approximated via backward differencing scheme: 𝑇 = − 𝜖𝑅 𝜖𝑢 ≈ − 3𝑅𝑜+1 − 4𝑅𝑜⋅ + 𝑅𝑜−1 2Δ𝑢

Since 𝑅𝑜+1 is unknown, we replace it by the best approximation 𝑅(𝑙+1) which is the solution in the next pseudo-time RK stage. 𝑅(𝑙+1) = 𝑅(0) − 𝛽𝑙+1Δ𝜐 𝑆k 1 + 𝛽𝑙+1 3Δ𝜐 2Δ𝑢 𝑇 ≈ − 3𝑅(𝑙+1) − 4𝑅𝑜⋅ + 𝑅𝑜−1 2Δ𝑢 We have Identify 1 + 𝛽𝑙+1

3Δ𝜐 2Δ𝑢 −1

as a point-implicit smoothing operator and inspired by this formulation, the modified RK/implicit smoother operator for time dependent problem is:

We note that this formulation is not the standard one* in which the original residual is given by 𝑆 k = 1 𝑊 𝐺

all faces

⋅ 𝑜 ⅆ𝑡 − 3𝑅(𝑙) − 4𝑅𝑜⋅ + 𝑅𝑜−1 2Δ𝑢 The original smoother operator is 𝐽 + 𝜁 Δ𝜐 𝑊 𝐵+ ⋅ 𝑜ⅆ𝑇

𝑏𝑚𝑚 𝑔𝑏𝑑𝑓𝑡

− α𝑙 3 2 Δ𝜐 Δ𝑢 + Δ𝜐 𝜖𝑇 𝜖𝑋 Δ𝑅 𝑚𝑝𝑑𝑏𝑚 = 𝑆𝑙 − 𝜁 Δ𝜐 𝑊 𝐵− ⋅ 𝑜 ⋅ Δ𝑅 𝑂𝐶ⅆ𝑇

𝑏𝑚𝑚 𝑔𝑏𝑑𝑓𝑡

𝐽 + 𝜁 Δ𝜐 𝑊 𝐵+ ⋅ 𝑜ⅆ𝑇

𝑏𝑚𝑚 𝑔𝑏𝑑𝑓𝑡

− 3 2 Δ𝜐 Δ𝑢 + Δ𝜐 𝜖𝑇 𝜖𝑋 Δ𝑅 𝑚𝑝𝑑𝑏𝑚 = 𝑆𝑙 − 𝜁 Δ𝜐 𝑊 𝐵− ⋅ 𝑜 ⋅ Δ𝑅 𝑂𝐶ⅆ𝑇

𝑏𝑚𝑚 𝑔𝑏𝑑𝑓𝑡

*Vatsa, V., & Turkel, E. (2003). Choice of variables and preconditioning for time dependent problems. In 16th AIAA Computational Fluid Dynamics Conference (p. 3692)

𝑆k = 1 𝑊 𝐺

all faces

⋅ 𝑜 ⅆ𝑡 − 3𝑅(0) − 4𝑅𝑜⋅ + 𝑅𝑜−1 2Δ𝑢 where the modified residual is

Shu-Osher problem – comparison of Convergence between old and new method

Standard method New method

Why DIRK?

• DIRK is a good choice for a high order scheme.
• Similar to BDF, DIRK requires only one forward time step.
• High order scheme allows larger time steps.
• Internal stages cost CPU time (good for GPUs).
• Allows varying time steps and hence adaptive time steps.

 Good for sharply varying in time transient problems.  Analysis for the time step is not needed.

Dual time stepping approach with Diagonally Implicit Runge- Kutta (DIRK)

We wish to solve the semi-discrete equation of the form 𝜖𝑣 𝜖𝑢 = −𝑆 𝑣 The DIRK time marching scheme is 𝑣(𝑘) = 𝑣𝑜 − Δ𝑢 𝑏𝑗𝑘𝑆 𝑗

𝑘 𝑗=1

𝑘 = 1 … 𝑡 𝑣(𝑜+1) = 𝑣𝑜 − Δ𝑢 𝑐

𝑘𝑆 𝑘 𝑡 𝑘=1

Where 𝑆 𝑘 = 𝑆 𝑣 𝑘 , 𝑣 𝑜 = 𝑣 𝑢𝑜 , 𝑣 𝑘 is the solution of the j-th RK step and also 𝑐j = 𝑏js ⟹ 𝑣𝑜+1 = 𝑣 𝑡 𝑏jj ≠ 0 ∀𝑘 > 1 𝑏11 = 0 ⟹ 𝑣 1 = 𝑣𝑜

We choose the DIRK3 scheme. The non-zero coefficients are given by

21 22 33 32 22 22 31 32 33 1 2 3

.788

1 3 2 6 1 12 (2 1) 1 For the lower order embedded scheme 5 3 ˆ 12 12 9 3 ˆ 12 12 2 2 3 ˆ 12 12

j js

a a a a a a a a a b a b b b                

𝑣(𝑘) = 𝑣𝑜 − Δ𝑢 𝑏𝑗𝑘𝑆 𝑗

𝑘 𝑗=1

Rewrite as 𝑆∗ 𝑣(𝑘) ≡ 𝑆 𝑘 − 𝑣𝑜 − 𝑣 𝑘 𝑏𝑘𝑘Δ𝑢 − 1 𝑏𝑘𝑘 𝑏𝑗𝑘𝑆 𝑗

𝑘−1 𝑗=1

= 0 (∗) Using a low-storage RK/Implicit smoother method we solve 𝜖𝑣 𝑘 𝜖𝜐 = 𝑆∗ 𝑣 𝑘 for 𝑣 𝑘 in pseudo-time until 𝑆∗ 𝑣 𝑘 is sufficiently small. The residual smoother operator is 𝐽 + 𝜁 Δ𝜐 𝑊 𝐵+ ⋅ 𝑜ⅆ𝑇

𝑏𝑚𝑚 𝑔𝑏𝑑𝑓𝑡

− 1 𝑏𝑘𝑘 Δ𝜐 Δ𝑢 Δ𝑅 𝑚𝑝𝑑𝑏𝑚 = 𝑆𝑙 − 𝜁 Δ𝜐 𝑊 𝐵− ⋅ 𝑜 ⋅ Δ𝑅 𝑂𝐶ⅆ𝑇

𝑏𝑚𝑚 𝑔𝑏𝑑𝑓𝑡

Adaptive time steps

• For a given DIRK scheme of order p, the embedded scheme is a scheme of

lower order, p-1 with the same 𝑏𝑗𝑘 coefficients and 𝑐 𝑗 instead of 𝑐𝑗

• Let 𝑣𝑜+1 be the solution of the given DIRK scheme after one time step of

size Δ𝑢 and let 𝑣 𝑜+1 be the solution of the embedded scheme, both with initial condition 𝑣𝑜, we have:

𝑣 𝑜+1 = 𝑣𝑜 + 𝑑 ∆𝑢𝑞−1 𝑣𝑜+1 = 𝑣𝑜 + 𝑑 ∆𝑢𝑞

We want to bound 𝑣

𝑜+1 − 𝑣𝑜+1 relative to the norm of the solution 𝑣𝑜 so that

𝑣 𝑜+1 − 𝑣𝑜+1 = 𝑣 𝑜 − 𝑣𝑜 ∙ Δ𝑢𝑜+1 Δ𝑢𝑜

𝑞−1

≤ 𝑈𝑃𝑀 ∙ 𝑣𝑜 Δ𝑢𝑜+1 ≤ Δ𝑢𝑜 ∙ 𝑈𝑃𝑀 ∙ 𝑣𝑜 𝑣 𝑜 − 𝑣𝑜

1 𝑞−1

The TOL parameter determines the required accuracy. A large TOL enables large time steps with low accuracy while smaller TOL yields a more accurate solution with smaller time steps (so an increased CPU time).

Lagrangian trajectories of burning droplets

• Droplets position and velocity vectors are calculated in a Lagrangian

manner with

ⅆ𝑦𝑞 ⅆ𝑢 = 𝑤𝑞 and ⅆ𝑤𝑞 ⅆ𝑢 = 𝑔𝐸 𝜐𝑣 𝑤𝑕 − 𝑤𝑞

Where 𝑔

𝐸 is the drag coefficient, 𝜐𝑣 is the droplet relaxation time

and 𝑤𝑕 is the velocity vector of the gas.

• A source term is added to the gas flow equation to account the

mass flux burning rate 𝑛

• f the droplets
• Energy equation

ⅆℎ ⅆ𝑢 = 𝑛

⋅ 𝐷𝑞 ⋅ 𝐻 𝑦 , 𝑢

• Continuity equation

ⅆ𝜍 ⅆ𝑢 = 𝑛

⋅ 𝐻 𝑦 , 𝑢

𝐻 𝑦 , 𝑢 = 1 𝜏3 2𝜌 3 exp − 𝑦 − 𝑦 𝑞 𝑢

2

2𝜏2 Where 𝐷𝑞 is the heat capacity of the burning gases escaping from the droplet and 𝐻 𝑦 , 𝑢 is a moving Green function which is taking in to account the finite size of the droplets cantered in its temporal position with a Gaussian width equal the droplet’s radii - 𝜏 = 3𝑛 4𝜌𝜍

1/3

 Riemann problem with perfect gas  Two-fluids Riemann problem  Shock-Sine wave interaction (Shu-Osher problem)  Shock wave/bubble interaction  Leidenfrost maze

Results

Sod’s shock tube - Riemann problem with perfect gas

One-dimensional problem with initial conditions:

𝑞 = 1 𝑄𝑇𝐽 𝑏𝑜ⅆ 𝑈 = 231: 11𝐿 𝑝𝑜 𝑢ℎ𝑓 𝑠𝑗𝑕ℎ𝑢 𝑡𝑗ⅆ𝑓 𝑝𝑔 𝑢ℎ𝑓 𝑢𝑣𝑐𝑓 𝑞 = 10 𝑄𝑇𝐽 𝑏𝑜ⅆ 𝑈 = 288: 89𝐿 𝑝𝑜 𝑢ℎ𝑓 𝑚𝑓𝑔𝑢 𝑡𝑗ⅆ𝑓

• Typical time step for BDF is 1e-6
• Time steps for Tol=1e-3 - Δ𝑢 ≈ 3 × 10−6
• Time steps for Tol=1e-2 - Δ𝑢 ≈ 1.3 × 10−5

Initial conditions

PL=68940 Pa ρL=1 kg/m3 Cp/Cv=1.1 ML=15 gr/mole PR=6894 Pa ρR=0.125 kg/m3 Cp/Cv=1.4 MR=29 gr/mole

Two-fluids Riemann problem

In the next slides we present the solution at t = 200ms

Velocity

Density Pressure

Shock-Sine wave interaction (Shu-Osher problem)

𝜍 = 3:857143; 𝑣 = 2:629369; P=10.33333 for x ≤ −0.4 𝜍 = 1 + 0.2 sin 5𝑦; 𝑣 = 0; 𝑄 = 1 for x > −0.4

One-dimensional problem with initial conditions: BDF - Δ𝑢 = 10−3 Tol = 2e-3 Δ𝑢 ≈ 2x10−3 Tol = 1e-2 Δ𝑢 ≈ 4x10−3

Shock-wave-bubble interaction (using two-fluids model)

Air, Pre-shock conditions Air, Post-shock conditions

Helium, Bubble

Shockwave direction M=1.22 Schematic description Initial conditions

Time steps evolution

• Typical time step for BDF is 1e-4
• For this problem the total CPU

time gain was approximately a factor of four.

Comparing results between methods

Mach number Density

The maze with the ratchet directions

Mach number contour map and streamlines

Conclusions and fu futu ture work

• We developed a method using dual time stepping with the combination of a

low-storage RK/Implicit residual smoother for pseudo time marching and DIRK for the physical time

• The parameter (TOL) controls the accuracy and also the CPU time
• The method requires less CPU time for problems with time scales which are

sharply varying

• DIRK based dual time stepping still requires a convergence and stability

analysis