An eﬃcient semi-implicit time integration method for extra large eddy simulations

(1)

An efficient semi-implicit time integration method for extra large eddy

simulations

B. Steerneman

X Y

Z

Grind around cylinder

SIRK-3A semi-implicit

SIRK-3A explicit

Institute for Mathematics and Computing Science

(2)

(3)

Master thesis

An efficient semi-implicit time integration method for extra large eddy

simulations

B. Steerneman

Supervisor:

A.E.P. Veldman

University of Groningen

Institute for Mathematics and Computing Science P.O. Box 800

9700 AV Groningen

The Netherlands October 2007

(4)

(5)

List of Figures

2.1 Stability regions of SIRK-3A, SIRK-3C and LSSIRK-3C . . . 13

2.2 Dispersion and dissipation of various methods: time discretisation only . . . 16

2.3 Dispersion and dissipation of various methods: time and space discretisation compared by varying k . . . 16

2.4 Dispersion and dissipation of various methods: time and space discretisation compared by varying the number of grid cells per wave length N_x . . . 16

2.5 Dispersion and dissipation of an implicit and explicit direction of B3: time and space discretisation compared by varying k . . . 17

2.6 Dispersion and dissipation of an implicit and explicit direction of B3: time and space discretisation compared by the number of gridcells per wave length . . 18

2.7 Dispersion and dissipation of an implicit and explicit direction of SIRK-3A: time and space discretisation compared by varying k . . . 18

2.8 Dispersion and dissipation of an implicit and explicit direction of SIRK-3A: time and space discretisation compared by the number of gridcells per wave length . . . 19

3.1 Convection-diffusion equation: time steady solution S(y) . . . 27

3.2 Convection-diffusion equation: grid . . . 29

3.3 Convection-diffusion equation: exact solution of the problem . . . 31

3.4 Convection-diffusion equation: convergences for the various (LS)SIRK methods 32 3.5 Burgers’ Equation: solution with ν = 1 . . . 34

3.6 Burgers’ Equation: convergence of the solution with second order spacial discretisation and its exact Jacobian using various (LS)SIRK methods . . . 35

3.7 Burgers’ Equation: convergence of the solution using fourth order spacial discretisation and second order Jacobian solved using various (LS)SIRK methods 36 3.8 Euler 1D: problem specification . . . 38

3.9 Euler 1D: close-up of grid cell Ωi . . . 38

3.10 Euler 1D tube: the test tube . . . 44

3.11 Euler 1D tube: time dependent solution . . . 44 3

(8)

3.12 Euler 1D tube: steady solution . . . 45

3.13 Euler 1D tube: test results with low CFL-number . . . 45

3.14 Euler 1D tube: test results with low CFL and residual tolerance . . . 45

3.15 Euler 1D tube: test results with high CFL-number . . . 46

3.16 Euler 1D tube: test results with high CFL and residual tolerance . . . 47

5.1 New Euler 1D Tube: results with ENSOLV with CFL = 2 . . . 54

5.2 New Euler 1D Tube: results with MATLAB with CFL = 2 . . . 55

5.3 New Euler 1D Tube: results with ENSOLV with CFL = 20 . . . 55

5.4 New Euler 1D Tube: convergence results with MATLAB and ENSOLV . . . 56

5.5 New Euler 1D tube: convergence results ENSOLV using the original ENSOLV Riemann boundary conditions . . . 56

5.6 Laminar cylinder: grid . . . 58

5.7 Laminar cylinder: the lift and drag coefficient up to the restart point . . . . 59

5.8 Laminar cylinder: the vortex street at the restart point . . . 59

5.9 Laminar cylinder: convergence results ENSOLV . . . 60

5.10 Laminar cylinder: the final point with B3 and 30 pseudo time steps . . . 62

5.11 Laminar cylinder: the final point after the restart with SIRK-3A and 5 Newton steps . . . 62

5.12 Laminar cylinder: the final point after the restart with SIRK-3A and 8 Newton steps . . . 63

(9)

Chapter 1 Introduction

CFD

The use of Computational Fluid Dynamics (CFD) is a common factor in the technical design industry these days. For example, the design of rockets, airplanes and Formula 1 cars would not be the same without CFD. When only the averaged drag coefficient or the averaged forces are needed, time independent simulations are sufficient. However, for a lot of applications knowing only the averaged values is not sufficient. One can think of the unsteady flow that comes from the central core of a space launcher which creates unsteady forces on its nozzle during take-off. Knowing only the averaged values makes it impossible to create a cost effective space launcher which is strong enough to withstand those unsteady forces. For those problems time dependent computations are necessary.

X-LES

For time dependent calculations of the dynamic flow, the extra-large eddy simulation (X-LES) model is often used. X-LES consists of a combination of the Reynolds-averaged Navier-Stokes (RANS) and large eddy simulation (LES) equations. RANS is used near the surface and LES is used in the other parts of the flow. X-LES simulations solve more details than RANS, without the cost of a full LES simulations.

Time dependent simulations always have restrictions in the time step size. For explicit calculations those restrictions are needed for the method to be stable, for implicit methods those restrictions are for accuracy considerations. Stability restrictions depend on the size of grid cells: the smaller the grid cell, the smaller the required time step size is needed.

On the surface the velocity is zero. Just above the surface (for example the wing of an aircraft) the velocity is high, which creates a very steep gradient in the velocity normal to the surface.

The area where this phenomenon is found is called the boundary layer. To get a proper representation of the boundary layer many grid cells are needed in the normal direction to the surface. This results in very fine grid cells in the normal direction, resulting in a very severe time step restricton. However, in the boundary layer the tangential direction does not have the steep gradient in velocity that the normal direction has. Therefore less grid cells are needed in the tangential direction, which makes the stability condition for the tangential direction less restrictive than the stability condition for the normal direction. For an explicit

5

(10)

time integration method, due to the restrictive stability conditions , the timestep has to be very small, which makes the computations expensive. Implicit calculations would require fewer time steps, but for large problems are very expensive, at least in terms of memory requirements.

SIRK

In the past few years, many semi-implicit methods were developed, for example the SIRK (semi-implicit Runge Kutta) methods, developed by X. Zhong et al (Ref [11]). This method enables separation of the computation of the direction normal to the surface and the direction tangential to the surface. The normal direction can be treated implicitly, while the tangential directions are treated explicitly. This removes the very restrictive stability conditions in the time step size due to the fine grid cells in the direction normal to the surface, but does not require a very expensive implicit solver at each time step.

The main purpose of this report is to investigate the useability of such a SIRK method in the RANS region of the X-LES model. This will be done using within the restrictions posed by requirements given in the next section.

ENSOLV

At the National Aerospace Laboratorium (NLR) an X-LES based CFD method has been developed. It is part of ENSOLV, the block-structured flow solver of the NLR. Therefore also a short introduction how X-LES is contained in ENSOLV is presented, which gives a few extra requirements for the SIRK methods.

Using X-LES, the RANS region is threated implicitly and the LES region is threaded explicitly. The current implicit time integration method is three-point backward (B3). This method is second order accurate. The set of implicit equations arriving from the time discretisation is solved in an iterative manner using a pseudo timestepping technique. The pseudo timestepping is performed by a semi-implicit scheme that is very similar to the SIRK methods, but uses an approximated Jacobian. The approximated Jacobian is a tridiagonal matrix, which is inverted with the Thomas algorithm. Inverting the approximated Jacobian is much faster than solving the complete Jacobian. It would be very favourable when this Jacobian could be used to solve the implicit term in the SIRK method. Since the current method is second order accurate, it is required that the new method is also second order accurate.

To get more simplicity in the coupling of the RANS and LES regions, it would be desirable when the semi-implicit could be used explicitly in the LES regions. Therefore the dissipation and dispersion of the time integration method needs to be considered also.

Requirements

The new method based on SIRK has to meet the following requirements:

• The method needs to be at least second order accurate in time,

• The dissipation and dispersion have to be comparable to the time integration methods now used (B3 for the implicit part and Runge-Kutta 4 for the explicit part),

(11)

7

• The approximated Jacobian used in the current line-implicit scheme should also work as Jacobian for the SIRK methods,

• The implicit solver has to be an order faster than the method with pseudo timesteps, now used for X-LES.

The requirements will be investigated extensively in this report.

Overview

After this introductory chapter there is a chapter about the theoretical aspect of the semi- implicit methods. Chapter 3 contains a few model problems solved with a semi-implicit method and in the end of that chapter a choice has been made for a certain method to be implemented. Chapter 4 is about the implementation of that method in ENSOLV. Chapter 5 contains some standard test cases with the newly implemented method in ENSOLV.

(12)

(13)

Chapter 2 Semi implicit time integration schemes

2.1 Semi implicit Runge Kutta methods

A theory about semi-implicit Runge-Kutta methods (SIRK) and its low-storage variants (LSSIRK) is introduced by Zhong respectively in reference [11] and reference [10]. In this section his argumentation is partly followed. As stated in the introduction many problems in CFD have a separable stiff and non-stiff part in their differential equations, which are treated separately. The equation ^du_dt = ¯f (u) is considered an already spatially discretised set of au- tonomous differential equations. The right hand side ¯f will be split in ¯f (u) = f (u) + g(u), where f (u) contains the non-stiff terms and g(u) contains the stiff terms. First some analysis on the SIRK schemes is done. The LSSIRK schemes are a special class of SIRK schemes with some extra restrictions on the coefficients. The most general way to write a r-stage SIRK scheme is,

uⁿ⁺¹ = uⁿ+ Xr j=1

ω_jk_j (2.1)

ki = δt



f (uⁿ+ Xi−1 j=1

bijkj) + g(uⁿ+ Xi−1 j=1

cijkj + diki)



 i ∈ {1, . . . , r}.

This scheme is called a method-A scheme. Method-A does not prescribe anything for the treatment of the implicit term. In method-C¹ the implicit term is linearised,

uⁿ⁺¹ = uⁿ+ Xr j=1

ω_jk_j (2.2)



I − δtdiJ(uⁿ+ Xi−1 j=1

c_ijk_j)



 ki = δt



f(uⁿ+ Xi−1 j=1

b_ijk_j) + g(uⁿ+ Xi−1 j=1

c_ijk_j)



 i ∈ {1, . . . , r}.

In these equations, J represents the Jacobian of g.

1it is called method-C and not method-B to be consistent with the notation in most of the literature (eg.

reference [11]), method-B is left out because it is not useful in this case

9

(14)

2.1.1 Accuracy conditions

First the accuracy conditions for the parameters of both methods are considered. This is done by comparing the Taylor series of f (u) and g(u) and the Taylor series of the methods.

Since third order accuracy is necessary, 3-stage SIRK schemes are considered and the Taylor series need to be determined up to the third order. This leads to the following conditions:

for first order accuracy for both method-A and method-C:

ω₁+ ω₂+ ω₃ = 1,

for second order accuracy for both method-A and method-C:

ω₂b₂₁+ ω₃(b₃₁+ b₃₂) = 1 2 ω₁d₁+ ω₂(d₁+ c₂₁) + ω₃(d₃+ c₃₂+ c31) = 1 2, for third order accuracy for both method-A and method-C:

ω₂(b²₂₁) + ω₃(b₃₂+ b₃₁) = 1 3 ω3b21b32 = 1 6 ω2(b21d2+ b21d1) + ω3(d1b31 + d2b32+ c21b32+ b21c32+ d3b31+ d3b32) = 1 3 ω₁d²₁+ ω₂(d²₂+ d₂c₂₁+ d₁c₂₁) + ω₃(d₁c₃₂+ d₂c₃₂+ c₂₁c₃₂+ d₃c₃₂+ d₃c₃₁+ d²₃) = 1 6, an extra condition for method-A third order accuracy:

ω₁d²₁+ ω₂(c₂₁+ d²₂) + ω₃(c₃₁+ c₃₂+ d₃)² = 1 3, and an extra condition for method-C third order accuracy

ω₂(c²₂₁+ 2d₂c₂₁) + ω₃(c₃₁+ c₃₂)²+ ω₃(2d₃(c₃₁+ c₃₂)) = 1 3.

When analysing these conditions, a few remarks must be made. At first, when only a second order accuracy is required, the same conditions can be used. Set ω₃ = 0 and only the conditions for first and second order are considered, and a 2-stage, second order accuracy SIRK method is created. As those conditions are the same for method-A and method-C, the same parameters can be used.

When the first and second order conditions are investigated more thoroughly, it can be seen that the parameters that represent the explicit part (bij) and those that represent the implicit part (c_ij and d_i) are in different equations. Thus, if only second order accuracy is required, an arbitrary second-order accurate explicit and an arbitrary second-order diagonally-implicit Runge-Kutta scheme can be used to created a SIRK-2A or SIRK-2C method as long as they have the same values for ω₁ and ω₂. When higher order Runge-Kutta schemes are used for

(15)

2.1. SEMI IMPLICIT RUNGE KUTTA METHODS 11 this, the result will still be second order, because then the coupling of the implicit and explicit part comes into play. Before the (LS)SIRK methods were introduced by Zhong, the coupling always was a problem, and often only second order accurate methods were used. A well- known example of such a method is the Adams-Bashforth/Crank-Nicolson (ABCN) method, which uses Adams-Bashforth for the explicit part and Crank-Nicolson for the implicit part.

2.1.2 Stability conditions

In addition to the accuracy conditions, a stability condition is needed to be sure that the stiff part is stable for all values of the time step δt. To obtain such an accuracy condition, a standard technique is used. Consider the following linear model equation,

du

dt = f (u) + g(u) = λ_fu + λ_gu. (2.3) Substitution of (2.3) in any of the SIRK methods gives the following amplification factor or stability function,

γ = uⁿ⁺¹ uⁿ = 1 +

Xr j=1

ω_jk_j (2.4)

k_j = δtλ_f(1 +Pi−1

j=1b_ijk_j) + δtλ_g(1 +Pi−1 j=1c_ijk_j)

1 − diδtλ_g . (2.5)

A problem is stable when for all combinations of eigenvalues λ_f and λ_g of the problem,

|γ(λf, λ_g)| < 1. In reference [11], Zhong introduced a stability condition which ensures that large values of λ_g in the left complex plane do not decrease the stability region for λ_f. Zhong demands that when

Re(δtλ_g) → −∞, γ(δtλ_f, δtλ_g) → 0. (2.6) In reference [11], Zhong proves that condition (2.6) is satisfied when

1 + Xr j=1

ω_jβ_j = 0 (2.7)

holds with

β_i = −1 d_i



1 + Xi−1 j=1

c_ijβ_j



 .

2.1.3 Low storage semi implicit method

Since ENSOLV usually works with large grids, it would be nice to have a low storage method.

A normal r-stage Runge-Kutta scheme occupies (r + 1) times the grid in the memory, once for uⁿ and r times for the Runge-Kutta stages. A low-storage scheme occupies only two times the grid in the memory: once for the r-th stage and once for the continuously updated Runge-Kutta stage. On the otherhand, low-storage schemes tend to have larger errors than

(16)

non-low-storage schemes. This will be evaluated during the model problems. The method-A variant (LSSIRK-rA) looks like:

k¯_j = a_jk¯_j−1+ δt

f (u_j−1) + g(u_j−1+ ¯c_jk¯_j−1+ c_jk_j) u_j = u_j−1+ b_j¯k_j

and the method-C variant (LSSIRK-rC) has the following form:

[I − δtcjJ(u_j−1+ ¯c_j)] ¯k_j = δt

f (u_j−1) + g(u_j−1+ ¯c_j¯k_j−1)

+ a_j[I − δtcjJ(u_j−1+ ¯c_j)] ¯k_j−1 uj = uj−1+ bj¯kj.

The variable u₀ represents uⁿ in the normal variant and u_r represents uⁿ⁺¹. Of course J represents the Jacobian of g again. The product of J and ¯k_j−1 is strictly not necessary, but gives better stability regions later on. When the LSSIRK methods are written as SIRK methods, their conversion parameters can be determined easily. Knowing them, the same stability and accuracy conditions can be used. For a three stage scheme the parameters should be defined as follows:

for both methods

ω₁ = b₁+ b₂a₂+ b₃a₃a₂ b₂₁ = b₁ d₁ = c₁ ω₂ = b₂+ b₃a₃ b₃₁ = b₁+ b₂a₂ d₂ = c₂

ω3 = b3 b32 = b2 d3 = c3

for method-A

c₂₁ = b₁+ ¯c₂+ c₂a₂

c₃₁ = b₁+ b₂a₂+ ¯c₃a₂+ a₂c₃a₂ c₃₂ = b₂+ ¯c₃+ c₃a₃

and for method-C

c21 = b1+ ¯c2

c₃₁ = b₁+ b₂a₂+ ¯c₃a₂ c₃₂ = b₂+ ¯c₃.

For a two-stage scheme, the parameters which have a subscript three can just be omitted.

As can be seen, because of the choice of an extra J in each Runge-Kutta stage, also the parameters for the second step already are different for method-A and method-C.

2.1.4 Parameters for various (LS)SIRK methods

In this section an overview is presented of the sets of parameters found by Zhong et al. for the method A and method C versions of the third order (LS)SIRK schemes.

(17)

2.1. SEMI IMPLICIT RUNGE KUTTA METHODS 13 SIRK−3A

0.4 0.6

0.6

0.8 0.8

0.8

1

Re(λ_f h) Im(λ f h)

−30 −2 −1 0

0.5 1 1.5 2

SIRK−3C

0.6 0.4 0.6

0.6

0.8

1

−30 −2 −1 0

0.5 1 1.5 2

LSSIRK−3C

0.6 0.4 0.6

0.60.8 0.8

0.8

1

1 1

−30 −2 −1 0

0.5 1 1.5 2

Figure 2.1: Stability regions of SIRK-3A, SIRK-3C and LSSIRK-3C SIRK-3A

w₁ = ¹₈ b₂₁ = ⁸₇ d₁ = ³₄ c₂₁ = ⁵⁵⁸⁹₆₅₂₄ w₂ = ¹₈ b₃₁ = ₂₅₂⁷¹ d₂ = ₂₃₃⁷⁵ c₃₁ = ₂₆₀₉₆⁷⁶⁹¹ w₃ = ³₄ b₃₂ = ₃₆⁷ d₃ = ₁₆₈⁶⁵ c₃₂ = −²⁶³³⁵₇₈₂₈₈

SIRK-3C

ω₁ = ¹₈ b₂₁ = ⁸₇ d₁ = 0.7970967740096232 c₂₁ = 1.058925354610082 ω₂ = ¹₈ b₃₁ = ₂₅₂⁷¹ d₂ = 0.5913813968007854 c₃₁ = ¹₂

ω₃ = ³₄ b₃₂ = ₃₆⁷ d₃ = 0.1347052663841181 c₃₂ = −0.3759391872875334

LSSIRK-3C

b₁ = ¹₄ a₂ = −¹₄ ¯c₂ = −1.143097033946135 c₁ = 2.267596813284564 b₂ = ²₉ a₃ = −²⁹₂₇ ¯c₃ = −2.031219208388789 c₂ = 2.685297589634163

b3 = 3 c3 = 2.309749357551431

The LSSIRK-3A method is not included, although it was tested during the stability tests.

The parameters Zhong gave, did not statisfy the stability condition (2.7). The larger the absolute value of λ_g became, the larger the absolute value of γ grew.

2.1.5 Stability regions

The stability regions of the considered methods are presented in this section. The stability region plots are based on the theory in section 2.1.2. For plotting these stability regions, the maximum absolute value of γ is computed for each λ_f over λ_g in the left complex plane. The results are shown in figure 2.1.

All three methods have an acceptable stability region. They all contain a part of the imaginary axis, which is important because the eigenvalues of high Reynolds number flows as appear in ENSOLV lie close to the imaginary axis. LSSIRK-3C seems to be the best, because it includes the largest part of the imaginary axis and even permits small positive values of λ_f near the imaginary axis.

(18)

2.1.6 Dissipation and dispersion

It is not only important that a method is stable, as discussed in section 2.1.2, a method should also have an acceptable dissipation (non physical damping) and dispersion (the phase error).

Dissipation and dispersion is usually defined as described below. When γ is the amplification factor of the method and γ_ex the exact amplification factor of the one dimensional version of (2.3), the quotient of both is considered in the following manner:

γ

γ_ex = |r|e^iφ,

where r is the dissipation and φ the dispersion. The exact solution for a problem has dissipation r one and the dispersion φ zero. The dissipation should always be less or equal than one, otherwise the problem is not stable. It is possible to create dissipation and dispersion figures with this definition, not considering the dissipation and the dispersion created by a spatial discretisation. For central convection discretisations, for example, it is usual to look at the dispersion and dissipation on the upper imaginary axis, so λ_f, λ_g ∈ [0, ∞]i.

Spatial discretisation

If the spatial discretisation is considered , some extra work needs to be done. Consider the convection differential equation

∂u

∂t +∂cu

∂x = 0, c > 0 (constant), x ∈ [0, L], t > 0 (2.8) u(x, 0) = U e^ikx, k = 2π

l , u(0, t) = u(L, t)

with L the wave length and k the wave number of the initial solution, which has the solution u(x, t) = U eîk(x−ct)= U eî(kx−ωt)= bu(t)eîkx,

with bu(t) = U e⁻^iωt and ω = ck. The time of one period will be written as T = ^2π_ω = _c^l. The spatial discretisation is defined arbitrary as long as it fits in the following semi-discretised form of the differential equation,

x_j = jδx, j = 0 . . . N, δx = L N duj

dt = c

δx X

m

α_mu_j+l (2.9)

Fourier analysis

The discrete solution for u_j = u(jδx, t) can be written as (analogous to the exact solution), u_j = bue^ikx^j = bue^iθj

(19)

2.1. SEMI IMPLICIT RUNGE KUTTA METHODS 15 with θ = kδx = ^2πδx_l = ^ωδx_c . θ is the wave number in the computational domain. Substituting the latter into equation (2.9) and dividing by e^iθj gives:

dbu dt = c

δx X

m

α_me^iθmbu = c

δxz(θ)bu = ωz(θ) θ bu with z(θ) =P

mα_me^iθm. The latter equation can be solved with a time integration method and divided by the exact solution for (2.8) to determine the dispersion and the dissipation of the spatial and time integration combined.

When the dissipation and dispersion rate of both the implicit and explicit part of a (LS)SIRK needs to be determined, one spatial dimension in (2.8) is not sufficient. It is fairly easy to add an extra spatial dimension y in (2.8),

∂u

∂t +∂c_xu

∂x + ∂c_yu

∂y = 0, c_x, c_y > 0 (constant), x, y ∈ [0, Lx] × [0, Ly], t > 0 u(x, y, 0) = U e^i(k^x^x+k^y^y), kx= 2π

l_x, ky = 2π l_y u(0, 0, t) = u(L_x, L_y, t)

which has the solution

u(x, t) = U eî(k^x^(x−c^x^t)+k^y^(y−c^y^t)) = U eî(k^x^x−ω^x^t+k^y^y−ω^y^t)= bu(t)eî(k^x^x+k^y^y),

with bu(t) = e⁻^i(ω^x^+ω^y^)t, ωx = cxkx and ωy = cyky. After doing the same steps as done with the one dimensional version the following equation for bu is found

dbu dt =

c_x

δxz(θ_x) + c_y δyz(θ_y)

bu =

ω_xz(θ_x)

θ_x + ω_yz(θ_y) θ_y

bu.

Comparing the spatial discretised version with the exact solution

When the discretised versions are compared, the CFL number ^cδt_δx is kept constant and θ is varied. This can be done in two ways. The first way is to fix the mesh (δx is constant) and vary the wave number k. Then e.g. after one time step δt the quotient of the amplification factors of the discretised version and the exact solution can be considered as a function of the wave number k.

The second way is to vary the mesh and fix the wave number. In this way the number of grid cells to accurately capture a wave length can be easily determined. Some more work is needed for this one. Since the CFL number is fixed and δx is varying, also δt has to vary. Let Nx

be the varying number of grid cells per wave length, N_x = _δx^l . Now the computational wave number can be written as θ = _N^2π

x. N_t = _δt^T represents the number of time steps needed for one period T. Knowing that c = _T^l, the CFL number can be written as CF L = ^δt/T_δx/l = ^N_N^x

t. Consider the amplification factors after a fixed time span T , which has in the discretised version a varying number of time steps Nt= _{CF L}^N^x = _{CF Lθ}^2π depending on the wave number in the computational domain θ.

(20)

0 0.2 0.4 0.99

0.992 0.994 0.996 0.998 1

iλ_f*δ t

Dissaption rate |r|

0 0.2 0.4

0 0.01 0.02 0.03 0.04

iλ_fδ t

Phase error φ

Real

Runge−Kutta 4 B3

(LS)SIRK−3A/C

Real

Runge−Kutta 4 B3

(LS)SIRK−3A/C

Figure 2.2: Dispersion and dissipation of various methods: time discretisation only

0 0.5 1 1.5

0.8 0.85 0.9 0.95 1

θ

Dissaption rate |r|

0 0.5 1 1.5

0 0.1 0.2 0.3 0.4 0.5 0.6

θ

Phase error φ

Real

Runge−Kutta 4 B3

(LS)SIRK−3A/C Real

Runge−Kutta 4 B3

(LS)SIRK−3A/C

Figure 2.3: Dispersion and dissipation of various methods: time and space discretisation compared by varying k

0 20 40

0.5 0.6 0.7 0.8 0.9 1

N

Dissaption rate |r|

0 20 40

0 2 4 6 8 10

N

Phase error φ

Real

Runge−Kutta 4 B3

(LS)SIRK−3A/C Real

Runge−Kutta 4 B3

(LS)SIRK−3A/C

Figure 2.4: Dispersion and dissipation of various methods: time and space discretisation compared by varying the number of grid cells per wave length N_x

(21)

2.1. SEMI IMPLICIT RUNGE KUTTA METHODS 17

b3: dissipation rate |r| (cfl_x,cfl_y) = (100,0.3)

0.9996 0.9997 0.9998 0.9999

0.9999

θ_x θy

0 0.2 0.4 0.6 0.8 1

x 10⁻³ 0

0.05 0.1 0.15 0.2 0.25 0.3 0.35

0.9994 0.9995 0.9996 0.9997 0.9998 0.9999

1 b3: phase error φ (cfl_x,cfl

y) = (0.3,100)

0.0005 0.0005

0.001 0.001

0.0015 0.0015

0.002 0.00250.003

θ_x

θy

0 0.2 0.4 0.6 0.8 1

x 10⁻³ 0

0.05 0.1 0.15 0.2 0.25 0.3 0.35

0 0.5 1 1.5 2 2.5 3 3.5 x 104 ⁻³

Figure 2.5: Dispersion and dissipation of an implicit and explicit direction of B3: time and space discretisation compared by varying k

2.1.7 Dissipation and dispersion figures

Explicit direction only

When a (LS)SIRK method is used in an X-LES model, it is likely that the explicit LES areas are also solved with a LS(SIRK) method. This is convenient because of the coupling of the explicit LES and semi-implicit RANS areas. Therefore it is necessary to compare the dissipation and dispersion of the explicit part of the (LS)SIRK methods with the methods now used in ENSOLV’s X-LES model: second order backward (B3) and Runge-Kutta 4.

Figures 2.2 to 2.4 show the comparisons described in the previous section using CF L = 1.

Figure 2.2 represents the dispersion and dissipation of the one dimensional version of (2.3).

Figures 2.2 and 2.3 respresent the convection equation (2.8) with a fourth order central spatial discretisation: the middle one with a varying wave number, the lower one with a varying number of grid cells per wave length.

The first remark that has to be made is that the results for LSSIRK-3C, SIRK-3C and SIRK- 3A are the same. This is because the explicit parts of those methods use exactly the same coefficients. It was expected that the results of the third order (LS)SIRK methods should be somewhere between the second order B3 results and the fourth order Runge-Kutta 4 results.

The (LS)SIRK results are much closer to the Runge-Kutta four results than to the B3 results, which is what was hoped for. In figure 2.4 can be seen, that per wave length twice the number of time steps are needed for the (LS)SIRK methods to get the same dispersion and dissipation on the same grid. Since CF L = 1, N in the figures also represents the number of time steps per wave length. With these figures it is expected that, from a dispersion and dissipation point of view, the (LS)SIRK methods are appropriate for the explicit LES part of a X-LES model.

(22)

b3: dissipation rate |r| (cfl_x,cfl_y) = (100,0.3)

0.450.50.50.550.550.60.60.650.650.70.70.750.750.80.8 0.85

0.85

0.9

0.90.9

0.95

0.950.95

0.95 0.95

Number of grid cells per wave length in implicit direction

Number of grid cells per wave length in explicit direction ^0.5 ¹ ^1.5 ² ^2.5 ³

x 10⁴ 10

20 30 40 50 60 70 80 90 100

0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9

0.95 b3: phase error φ (cfl_x,cfl_y) = (0.3,100)

11

1 1

22

2 2

3 3

4 4

5 6 5

x 10⁴ 10

20 30 40 50 60 70 80 90 100

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6

Figure 2.6: Dispersion and dissipation of an implicit and explicit direction of B3: time and space discretisation compared by the number of gridcells per wave length

sirk3a: dissipation rate |r| (cfl_x,cfl_y) = (100,0.3)

1 1

1

1 1

θ_x θy

0 0.2 0.4 0.6 0.8 1

x 10⁻³ 0

0.05 0.1 0.15 0.2 0.25 0.3 0.35

1 1 1 1 1 1 1

sirk3a: phase error φ (cfl_x,cfl_y) = (0.3,100)

1e−05 2e−05 3e−05

4e−05 5e−05 6e−05 7e−05 8e−05 9e−05

θ_x

θy

0 0.2 0.4 0.6 0.8 1

x 10⁻³ 0

0.05 0.1 0.15 0.2 0.25 0.3 0.35

0 1 2 3 4 5 6 7 8 x 109 ⁻⁵

Figure 2.7: Dispersion and dissipation of an implicit and explicit direction of SIRK-3A: time and space discretisation compared by varying k

(23)

2.1. SEMI IMPLICIT RUNGE KUTTA METHODS 19

sirk3a: dissipation rate |r| (cfl x,cfl

y) = (100,0.3)

0.750.750.80.80.850.850.90.90.950.95

1

1 1

1

x 10⁴ 10

20 30 40 50 60 70 80 90 100

0.7 0.75 0.8 0.85 0.9 0.95

1 sirk3a: phase error φ (cfl_x,cfl

y) = (0.3,100)

1 1

2 2

34 34

5 6 5

x 10⁴ 10

20 30 40 50 60 70 80 90 100

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6

Figure 2.8: Dispersion and dissipation of an implicit and explicit direction of SIRK-3A: time and space discretisation compared by the number of gridcells per wave length

Implicit and explicit direction

Ofcourse, also the dispersive and dissipative behaviour of the (LS)SIRK methods have to be considered in the RANS, semi-implicit, areas of the problem. Therefore, also some two dimensional plots are made considering both the implicit and explicit direction of the method.

To simulate the high aspect ratios of the grid cells typically found in the RANS area, the CFL number in the implicit direction is more than three hunderd times larger than the CFL number in the explicit direction. Just as in the previous figures the methods are compared with the methods already present in ENSOLV. Because of the high CFL number, it is not possible to compare the results with explicit Runge-Kutta methods, so only the second order backward (B3) is considered. To avoid too much figures in the report, only time and space discretisation together are considered. First with varying wave number k and secondly with a varying number of grid cells per wave length. Since the results of all (LS)SIRK methods are almost the same, only the results of SIRK-3A will be shown. In the figures a CFL number of 0.3 in the explicit direction is used, this is because the central spatial discretisation has pure imaginary eigenvalues and is not stable for pure imaginary eigenvalues greater than about 0.35 . In practice the eigenvalues have small negative parts due to the artificial dissipation and the diffusive term, so CFL numbers up to 1.5 can be used. More about this can be found in section 2.1.2.

Figures 2.5 to 2.8 show the results. Both B3 and SIRK-3A are almost symmetric in the line θ_x

= θy in the dissipation figures. Knowing this, the same conclusions can be made as with the one dimensional dissipation figures. The dissipation of the SIRK-3A method is significantly smaller than the dissipation of the B3 method. Also the dispersion of the SIRK-3A scheme is smaller than the dispersion of the B3 scheme. From a dissipation and dispersion point of view, (LS)SIRK methods are very suitable for the RANS areas.

(24)

2.2 The implicit term

To be able to use an A variant (as described in equation (2.1)) of a (LS)SIRK method, a method to solve the implicit term needs to be determined. There are a number of possibilities.

Three are considered here.

2.2.1 Exact solution

The first one is to determine an exact solution of the implicit equation. This can be very expensive, because in ENSOLV and a few of the test cases this is a set of non-linear equations.

This method can be used for example to determine a reference solution to which the other methods can be compared. In practice for model problems in MATLAB, a builtin non-linear equation solver can be used.

2.2.2 Linearised

The second one, a quite simple one, is to use the method used in the C variants as described in equation (2.2). A difference between the A and C variant is that the coefficients in the A variants are choosen in such a way that an exact solution of the implicit term is expected, while in the C variants the coefficients are choosen such that only this linearisation is considered. With non-linear problems it is expected that this treatment of the implicit term is not sufficient for the A variants.

2.2.3 The Newton method for solving a system of non-linear equations The Newton method is described more extensively. The Newton method for system of non- linear equations is a generalisation of the Newton method for the one-dimensional case for solving h(x) = 0. The one-dimensional case can be described by the following fixed-point method:

x^(l)= x^(l−1)− φ(x^(l−1))h(x^(l−1)),

where φ(x) = 1/h⁰(x). In for example reference [1] this method is derived and there is also proven that this method has quadractic convergence to h(x) = 0 when the initial estimation x⁽⁰⁾ is close enough to the real root.

The multi-dimensional case for solving H(x) = 0, where x and H is a vector now, can be constructed by replacing h⁰(x) with the Jacobian of H(x) denoted as K(x). Which results in:

x^(l)= x^(l−1)− K⁻¹(x^(l−1))H(x^(l−1)).

Again a proof of quadratic convergence when the initial estimation x⁽⁰⁾ is close enough to the real root is given in reference [1].

(25)

2.2. THE IMPLICIT TERM 21

The Newton methods for (LS)SIRK-rA

Now H(x) and K(x) will be determined for a Runge-Kutta stage in a (LS)SIRK-rA method.

First the equation to be solved for each Runge-Kutta stage, equation (2.1), is recalled:

ki = δt



f (uⁿ+ Xi−1 j=1

bijkj) + g(uⁿ+ Xi−1 j=1

cijkj+ diki)



 .

So each Runge-Kutta stage x := k_i is the variable to be solved, which leads to

H(x) = x − δt



f (uⁿ+ Xi−1 j=1

b_ijk_j) + g(uⁿ+ Xi−1 j=1

c_ijk_j+ d_ix)



 ,

where all the other variables and functions are known. The Jacobian matrix of H(x) can be determined now:

K(x) = I − δtdiJ(uⁿ+ Xi−1 j=1

c_ijk_j + d_ix),

where I is the N -dimensional identity matrix, J the Jacobian matrix of g, N the dimension of the vector x and n the indicator of the previous time step.

Two things additionally need to be considered. First what initial estimate x₍₀₎ is used. Since ki represents a part of the difference between two succesive time steps and the time steps are small, it will be assumed that the difference is not so big and x⁽⁰⁾ = 0 is close enough to the real difference between two succesive time steps.

Furthermore a stopping criterion is needed. Two kinds of stopping criteria are used in the model problems. The first one is to stop the Newton method after a fixed number of steps.

The second one is to stop the Newton method at step l when kH(x^(l))k or kx^(l)− x^(l−1)k is smaller than a certain tolerance . The drawback of using the easier kx^(l)− x^(l−1)k ≤ is that when two succesive Newton steps are close to each other but not to the root of H(x) it also stops. That is why kH(x^(l))k ≤ or a combination of both is used.

In order to be able to compare different time step sizes with the same base tolerance ^∗, is scaled with the time step size (δt)^s, where s is the order of convergence of the time integration method. So the residue of the Newton method stays scales with the error introduced by the time integration method. When kxk is not of order one, can also be scaled with kxk.

Resulting in

= ^∗(δt)^skx^k−1k. (2.10)

A major drawback of the Newton method is that, unlike for example the bisection method, it does not always converge to the nearest root. On the other hand, in the derivation of the Newton method it can be seen that the Jacobian can be approximated in stead of determined exactly as presumably needed in the C variant methods. In the test cases in the next chapter it will be determined whether the Newton method with x⁽⁰⁾ = 0 works as non-linear solver for the implicit part of (LS)SIRK-rA methods.

(26)

2.3 Computational efforts

In this section the theorical computational effort for ENSOLV for the current method and for the different (LS)SIRK methods are determined. In the next chapter, when for example the number of Newton steps in the Newton method are known, a comparison in speed can be made. The comparison is based on the number of residual computations that are made in the implicit block and explicit block and the number of approximated Jacobian inversions that are made in the implicit block. The computation effort in time of an implicit block, an explicit block and a Jacobian inversion are represented by respectively R_i,R_e and I_i.

2.3.1 Current line-implicit method with coupling

In the current scheme each physical time step N_p pseudo time steps with N_sⁱ Runge-Kutta stages are done in the implicit block, in each of the Runge-Kutta stages one residual is computed and one approximated Jacobian inversion is done. In the explicit part N_s^e Runge- Kutta stages are done, each of those steps one residual is computed. This leads to a total of

Ccurrent = N^pN_sⁱ(R_i+ I_i) + N_s^eR_e computations per physical time step.

2.3.2 (LS)SIRK methods

The (LS)SIRK methods are divided in two groups. First SIRK-3A with the implicit term linearised, SIRK-3C and (LS)SIRK-3C, they all require one approximated Jacobian inversion and one residual computation for the implicit block. The same method is used for the explicit block, but then only the residual computation is needed. Each physical timestep N_s^s Runge- Kutta stages are taken, which gives a total of

Csirk = N^s^s(R_i+ R_e+ I_i)

computations per physical time step. Secondly, when the Newton method is used for the implicit term each Runge-Kutta stage N_n Newton steps are needed to compute the implicit term. Which gives a total of

Csirknewton = N^s^s(Nn(Ri+ Ii) + Re) computations per physical time step.

2.3.3 Comparison

It is hard to compare the current line-implicit method and the (LS)SIRK methods. The time step sizes needed for a reasonable answer can be different as was shown in the dispersion and dissipation figures. For different problems the computational ratio between the implicit and explicit block are not the same and the number of Newton steps needed for a Newton based SIRK method is not known. A few example cases are compared.

(27)

2.4. DISCUSSION 23 Totally implicit with the same grid size

Now R_e = 0, which leads to Csirknewton = NⁿCsirk. So using a non-Newton (LS)SIRK method is N_n times cheaper than using a SIRK-3A with Newton. Also R_i and I_i are the same for the current method and the (LS)SIRK methods. Which gives the following ratio,

γimp = Ccurrent

Csirknewton = N_pN_sⁱ

N_nN_s^s. (2.11)

Usually N_sⁱ = 5, N_p = 100 and since in this report only 3 stage (LS)SIRK methods are considered N_s^s= 3. If the approximated Jacobian inversion are accurate enough for the non Newton methods (Nn = 1), γimp ≈ 167, so those methods would be 167 times faster than the current line-implicit method. Since N_n is not known yet, approximations for SIRK-3A with Newton will be made in the next chapter.

Coupling test case

In reference [7] Scheijbeler, did a test run with the current line implicit method. The per- formance results of that test case can be found on page 87 of reference [7]. Those results are used now to compare the theorical speed of the (LS)SIRK methods and the current line implicit method with coupling. From Scheijbelers results can be derived that for the current method and that test case R_i + I_i ≈ Re. Scheijbeler used N_p = 50 pseudo time steps and N_sⁱ = 5 and N_s^e = 4 Runge-Kutta stages. Again it is assumed that N_s^s = 3. Using the same grid size for the (LS)SIRK methods, now the following ratio is obtained,

γcoup = Ccurrent

Csirknewton = N_pN_sⁱ+ N_s^e

(N_n+ 1)N_s^s ≈ 85

N_n+ 1. (2.12)

When the non-Newton methods work, so Nn = 1, γcoup ≈ 43. The estimates made for SIRK-3A with the Newton method are made in the next chapter.

2.4 Discussion

Three methods discussed in this chapter are worth trying in the test cases: LSSIRK-3C, SIRK-3A and SIRK-3C. Low storage methods are preferred above non low storage methods, because they use less memory. The C variant methods are most likely quicker than the A variant methods because per Runge-Kutta only one time the Jacobian has to be inverted, while in the A methods a non-linear system of equations has to be solved. In section 2.2 a few methods for treating this implicit term are introduced, where solving it exact is an option for the simple test problems, but not for the more complicated problems typically solved by ENSOLV, because the expense in terms of computer time. Only linearising the term and solving that, will most likely not give third order results for non-linear problems. So the Newton method is in practice the only option left from the three.

In terms of memory use and computer time LSSIRK-3C is the most interesting method, providing it performs well with an approximated Jacobian. Followed by SIRK-3C with the same requirement with respect to the Jacobian. Most likely SIRK-3A in combination with

(28)

the Newton method will work for approximated Jacobians, at the cost of an iterative process each Runge-Kutta stage.

In the next chapter it will be evaluated,

• whether the methods converge third order,

• how sensitive the Newton method and the C variants are for approximating Jacobians (as used in ENSOLV)

• and which parameters and how many steps are needed for the Newton method to converge.

(29)

Chapter 3 Testing (LS)SIRK methods on model problems

3.1 Model problem: Convection-diffusion equation

For analysing semi-implicit time integration methods a model is needed that,

• has an exact solution, that can be determined,

• can be split in a stiff part in one direction and a non stiff part in another direction,

• has some sort of ‘difficulty’ at one of the boundaries,

• is a time dependent problem,

and since it is a model problem, it has to be as simple as possible.

A suitable problem is the convection-diffusion equation

∂u

∂t + R = 0 with

R = div(a u − µ grad u)

and a a vector which contains the convective velocities. When the incompressibility condition (div a = 0) is assumed, the equation can be written with R in its more common form

R = a grad u − div(µ grad u).

The two dimensional version of this equation for u(x, y, t) on the region x, y ∈ [0, 1] × [0, 1]

and t ∈ [0, ∞] can be writen as

∂u(x, y, t)

∂t = µ∂²u(x, y, t)

∂y² − a∂u(x, y, t)

∂x − b∂u(x, y, t)

∂y (3.1)

25

(30)

with the following initial and boundary conditions

u(x, 0, t) = 0 (3.2)

u(x, 1, t) = 1 (3.3)

u(0, y, t) = u(1, y, t) (3.4)

u(x, y, 0) = f (x, y) (3.5)

The y-direction is considered to be the stiff one. The x-direction is the non-stiff one. To make the x-direction even less stiff and to avoid any numerical restrictions, the diffusion constant in that direction is zero.

3.1.1 Analytical solution

The problem is not homogeneous, which makes it more difficult to solve. On physical grounds it is assumed that the solution will converge in time to a time steady solution,

S(y) = lim

t→∞u(x, y, t),

which is also constant in the x-direction. The solution can be written as u(x, y, t) = S(y) + v(x, y, t). Under those assumptions the boundary conditions for u(x, y, t) can be rewritten for v(x, y, t) and will become

v(x, 0, t) = v(x, 1, t) = 0 (3.6)

v(0, y, t) = v(1, y, t) (3.7)

v(x, y, 0) = f (x, y) − S(y) (3.8)

The time steady solution

First the solution for S(y) will be determined. Since S(y) only depends on y and is time independent, the equations simplify to

b∂S

∂y − µ∂²S

∂y² = 0

with S(0) = 0 and S(1) = 1. The solution of this problem is given for example in reference [8] and is

S(y) = 1 − e⁻^µ^b ^y 1 − e⁻^µ^b

.

Figure 3.1 shows that for _µ^b 1 the solution has a small boundary layer at y = 1, which is desirable for the model problem.

An eﬃcient semi-implicit time integration method for extra large eddy simulations

An efficient semi-implicit time integration method for extra large eddy

simulations

B. Steerneman

Institute for Mathematics and Computing Science

Master thesis

An efficient semi-implicit time integration method for extra large eddy

simulations

B. Steerneman

Supervisor:

A.E.P. Veldman

University of Groningen

Institute for Mathematics and Computing Science P.O. Box 800

9700 AV Groningen

The Netherlands October 2007

Contents

List of Figures

Chapter 1

Introduction

Chapter 2

Semi implicit time integration schemes

2.1 Semi implicit Runge Kutta methods

2.2 The implicit term

2.3 Computational efforts

2.4 Discussion

Chapter 3

Testing (LS)SIRK methods on model problems

3.1 Model problem: Convection-diffusion equation