Comparison of discontinuous Galerkin time integration schemes for the solution of flow problems with deformable domains

Applied Mathematics and Computation xxx (2015) xxx–xxx

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Comparison of discontinuous Galerkin time integration schemes for the solution of flow problems with deformable domains Ondrˇej Bublík ⇑, Jan Vimmr, Alena Jonášová European Centre of Excellence, NTIS – New Technologies for the Information Society, Faculty of Applied Sciences, University of West Bohemia, Univerzitní 8, CZ-306 14 Pilsen, Czech Republic

a r t i c l e

i n f o

Keywords: Implicit method Local time stepping method Discontinuous Galerkin finite element method Compressible Navier–Stokes equations ALE method

a b s t r a c t The computational efficiency of numerical solvers is known to strongly depend on the chosen time integration scheme. Thus, when solving viscous flow problems on time-varying domains, an efficient and reliable solver is one of the prerequisites for a successful solution. For this reason, we provide a comparison of two time integration schemes in terms of numerical error and CPU time usage. The contribution of the paper lies in the mutual comparison of two different approaches for time integration of the compressible Navier–Stokes equations in arbitrary Lagrangian– Eulerian (ALE) description that are spatially discretised by the discontinuous Galerkin finite element method. The computational ability of implicit Crank–Nicolson scheme and explicit three-step Runge–Kutta scheme, computational performance of which is improved by the local time-stepping technique, are tested on two flow problems: the propagation of an isentropic vortex, which has a known analytical solution, and the viscous flow around an oscillating NACA 0012 airfoil. All numerical simulations are carried out on unstructured triangular meshes by using a continuous mapping between the reference and time-varying domains. Ó 2015 Elsevier Inc. All rights reserved.

1. Introduction During the last decade, the discontinuous Galerkin (DG) finite element method, first proposed by Reed and Hill [1], has become very popular for the solution of various problems from the field of fluid mechanics [2–4], electrodynamics or electromagnetism [5–7] and plasma physics [8]. The reasons behind the popularity can be attributed to the method’s ability to achieve high order of spatial accuracy combined with low artificial damping, robustness and overall stability, thus, making it an ideal method for the simulation of laminar and turbulent flows with complex vortical structures. Although able to produce stable and high-order accurate solutions on fully unstructured meshes, the DG method is associated with a large number of unknowns, which with increasing problem size can substantially increase the computational demands. For this reason, the choice of a time integration method is of crucial importance [9,10]. Most often, the time discretisation is carried out with an explicit time integration scheme that represents the simplest choice. It is, however, well known that the computational efficiency of explicit schemes is directly affected by the maximal ⇑ Corresponding author. E-mail address: [email protected] (O. Bublík). http://dx.doi.org/10.1016/j.amc.2015.03.063 0096-3003/Ó 2015 Elsevier Inc. All rights reserved.

Please cite this article in press as: O. Bublík et al., Comparison of discontinuous Galerkin time integration schemes for the solution of flow problems with deformable domains, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.03.063

2

O. Bublík et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx

size of the global time-step restricted by the CFL condition of stability. This restriction becomes particularly notable in computational meshes with large variations in element size that are very common in computational practice. For example, when a local mesh refinement is needed due to the presence of sharp corners in the geometry or due to the presence of shock waves. In this case, the smallest element can significantly reduce the size of the global time step and, thus, hamper the overall computational efficiency. To overcome the disadvantage, a local time-stepping (LTS) technique can be used, see, e.g., [7]. As the name suggests, the LTS method utilises the principle of local time steps that are computed for each element independently. In this way, a massive reduction of CPU time can be achieved without the need to significantly alter the computational algorithm developed for an explicit scheme. Another way, how to overcame the restriction imposed by the CFL condition, is to introduce an implicit time integration method with a large interval of stability. Although this approach may seem straightforward, it is a fact that the application of implicit schemes gives rise to several problems. Among them, it is possible to mention the need to solve a system of nonlinear algebraic equations during every time iteration, which in large size problems may negatively affect the CPU time. It is also not insignificant that implicit schemes are generally difficult to program and cannot be easily parallelised as the explicit ones. Moreover, meshes with elements of highly disparate size are known to lead to ill-conditioned problems with implicit schemes. By taking into consideration all the advantages and disadvantages of explicit and implicit schemes, a question naturally arises, whether the modification of explicit schemes in the form of the LTS method can provide higher or at least comparable computational performance than their implicit counterparts. For unsteady 3D viscous flow simulations in non-deforming domains, an answer to this question can be found in [9], where the authors compared various Jacobian free implicit methods with selected explicit ones. On the basis of the presented results, the authors concluded that explicit schemes can achieve approximately the same computational performance as the implicit ones when the LTS technique is applied, as well. By following up the work of Birken et al. [9], we try to find an answer to the question above by solving unsteady flow problems in time-varying domains, where the numerical solution is also dependent on the mapping procedure between the fixed reference configuration and the real domain. To determine the impact of time integration on CPU time usage and numerical error, we consider two well-known explicit and implicit time integration schemes – the explicit three-step Runge–Kutta method, performance of which is improved by the local time-stepping method, and the implicit Crank– Nicolson method. The structure of the paper is outlined as follows: First, the equations governing the viscous flow in a time-varying domain are introduced and discretised in the sense of the DG method. Next, the time integration with the explicit LTS and implicit schemes is carried out, followed by a description of the mapping procedure necessary for the solution of flow problems on moving domains. In the next section, the computational ability of the developed explicit and implicit solvers is first verified for an inviscid model problem consisting of a moving isentropic vortex [11] and then used for the numerical simulation of viscous flow around an oscillating NACA 0012 airfoil, motion of which is prescribed by a time-dependent function. Finally, we summarise our work, give some concluding remarks and provide a brief outline of future work. 2. Mathematical model and DG discretisation 2.1. Governing equations The mathematical model describing the compressible viscous flow in a two-dimensional time-varying computational domain consists of the non-linear conservative system of the Navier–Stokes equations in ALE form. The equations written in the compact vector form are given as

 X 2  2 D Aw X @ @w @ v ¼ þ f s ðwÞ  U s f ðw; rwÞ; Dt @xs @xs @xs s s¼1 s¼1 where

DA Dt

ð1Þ

denotes the so-called ALE time derivative, ½x1 ; x2  ¼ x 2 XðtÞ  R2 ; t 2 ½0; T ; wðx; tÞ ¼ ½.; .u1 ; .u2 ; ET is the vector

. is the density, u ¼ ½u1 ; u2 T is the velocity vector, E is the total energy per unit volume, T f s ðwÞ ¼ ½.us ; .us u1 þ ds1 p; .us u2 þ ds2 p; ðE þ pÞus T and f vs ðwÞ ¼ ½0; s1s ; s2s ; u1 s1s þ u2 s2s þ k@T=@xs  are the inviscid and viscous fluxes, respectively. The stress tensor sij is defined as

of conservative variables,

sij ¼ l

  @ui @uj 2 @uk þ  dij ; @xj @xi 3 @xk

i; j ¼ 1; 2;

ð2Þ

where l is the dynamic viscosity and dij is the Kronecker delta. By U s we denote the components of the domain velocity. The system of governing equations (1) is completed by the constitutive relations for ideal gas and temperature

" p ¼ ðj  1Þ E 

# 2 1 X . u2s ; 2 s¼1

where p is the pressure,

k

  @T j l @ p ; ¼ @xs j  1 Pr @xs .

ð3Þ

j ¼ 1:4 is the adiabatic index and Pr ¼ 0:72 is the Prandtl number.

Please cite this article in press as: O. Bublík et al., Comparison of discontinuous Galerkin time integration schemes for the solution of flow problems with deformable domains, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.03.063

3

O. Bublík et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx

The system of equations (1) is equipped with the initial condition

wðx; 0Þ ¼ w0 ðxÞ;

x 2 Xð0Þ

ð4Þ

and the boundary conditions:  At the inlet, the density prescribed.

.i , velocity vector ui , zero traction

 At the outlet, the static pressure po , zero traction

P2

s

s¼1 ks ns

P2

s

s¼1 ks ns

¼ 0; k ¼ 1; 2 and zero normal derivative

@T @n

¼ 0 are

@T ¼ 0; k ¼ 1; 2 and zero normal derivative @n ¼ 0 are prescribed.

@T ¼ 0 is prescribed and the velocity components us ; s ¼ 1; 2 are set to correspond to  On the wall, zero normal derivative @n the velocity of the moving domain boundary.

2.2. Discretisation The spatial discretisation of Eq. (1) is carried out with the discontinuous Galerkin finite element method (DGFEM). The viscous fluxes f vs ðwÞ are discretised in dependence on the chosen time integration scheme. Namely, in the case of an explicit scheme, we apply the easily implementable local discontinuous Galerkin (LDG) scheme as proposed by Cockburn and Shu [12]. By contrast, the implicit schemes make it difficult to implement the LDG scheme, thus, compelling us to use the incomplete interior penalty (IIP) scheme, see [13–18]. Let us divide the time-varying computational domain X  R2 into a set of non-overlapping triangles (control elements) Xk S with boundary @ Xk , i.e., it holds k Xk ¼ X. Let h ¼ maxk fdiamðXk Þg, the solution is considered on the space Sh ¼ ½Sh 4 , where

n o Sh ¼ g 2 L2 ðXÞ : gjXk 2 P q ðXk Þ; gjXnXk ¼ 0; 8Xk :

ð5Þ

Here, P q ðXk Þ; q P 0 denotes the space of polynomials of degree 6 q on the element Xk . The basis functions of the space P q ðXk Þ are chosen as a system of orthogonal polynomials [19,20] defined on the reference triangular element, Fig. 1, as

Pij ðn; gÞ ¼ p0i



 i 2n 1g pj2iþ1 ðgÞ; 1g 2

ð6Þ

where

 n 1 d  a 2 n ð1  bÞ ðb  1Þ a n 2 n!ð1  bÞ db n

pan ðbÞ ¼

is the nth Jacobi polynomial. For the LDG scheme [12,2], Eq. (1) is rewritten to a system of first order partial differential equations using an auxiliary variable for the gradient

h ¼ rw:

ð7Þ

By scalar multiplication of Eqs. (1) and (7) by a test function v , integrating over a control element Xk and using the Green’s theorem, we get the following integral identities:

Fig. 1. Reference triangle according to [19,20].

Please cite this article in press as: O. Bublík et al., Comparison of discontinuous Galerkin time integration schemes for the solution of flow problems with deformable domains, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.03.063

4

O. Bublík et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx

d dt

I 2 X @v ðf s ðwÞ  wU s Þ  dX þ ðF ðwþ ; w ; nÞ  wU n Þ  v dl @xs Xk Xk s¼1 @ Xk Z X I X I 2 2 @v f vs ðw; hÞ  dX þ hf vs ðw; hÞins  v dl  c ½w  v dl; ¼ @xs Xk s¼1 @ Xk s¼1 @ Xk

Z

Z

w  v dX 

hs  v dX þ

Z

Xk

Z

w

Xk

@v dX  @xs

I

hwins  v dl ¼ 0;

s ¼ 1; 2;

ð8Þ

ð9Þ

@ Xk

P where n ¼ ½n1 ; n2 T is the outward unit vector normal to the boundary @ Xk ; U n ¼ 2s¼1 U s ns is the normal domain velocity, c is þ  1 the penalty parameter, hwi ¼ 2 ðw þ w Þ denotes the mean value of quantity w; ½w ¼ ðwþ  w Þ represents the jump of w and

F ðwþ ; w ; nÞ ¼ f n ðwÞ ¼

2 X

f n ðwþ Þ þ f n ðw Þ 1  kmax ðwþ  w Þ; 2 2

f s ðwÞns

s¼1

pffiffiffiffiffiffiffiffiffiffiffiffi is the Lax–Friedrichs flux [21] with kmax ¼ jU n j þ a, where a ¼ jp=. is the local speed of sound. If a part of the control element boundary @ Xk lies on the domain boundary @ XðtÞ, it is necessary to adjust the inviscid and viscous fluxes F and f vs ; s ¼ 1; 2 and to ensure the satisfaction of boundary conditions for w on the basis of the boundary conditions defined on @ XðtÞ. The implementation is carried out in the sense of [2]. The volume and curve integrals in identities (8) and (9) are computed with the help of the Gauss integration rules of appropriate order. The solution w and the test functions v are considered to be in the space Sh defined by Eq. (5). The elements wm of the vector w are taken as a linear combination of basis functions

wm ¼

N X wim ui ;

m ¼ 1; . . . ; 4:

ð10Þ

i¼1

The substitution of wm and

v

into Eqs. (8) and (9) yields the following system of ordinary differential equations:

dðMWÞ ^ s Þ; ¼ RðW; U dt

ð11Þ

where MðtÞ is a time-dependent diagonal matrix with mass matrices on its diagonal, W is the vector of unknown coefficients ^ s is the domain velocity vector defined in nodal points of the triangular mesh. of linear combinations wim and U 2.2.1. Explicit local time-stepping method For the explicit time integration of Eq. (11), the two-step Runge–Kutta scheme of second order accuracy is employed

 1   ^ nÞ ; W1 ¼ Mnþ1 Mn Wn þ DtRðWn ; U s   1 1 1 1 ^ nþ1 Þ ; Wnþ1 ¼ Mnþ1 Mn Wn þ Mnþ1 W1 þ DtRðW1 ; U s 2 2 2

ð12Þ

where Mn ¼ Mðtn Þ and Mnþ1 ¼ Mðt n þ DtÞ. As already mentioned in Introduction, the global time step size of the explicit Runge–Kutta scheme (12) is restricted by the CFL condition of stability

8 > > < CFL Dt 6 min  k > jk jþjk > 4 : xkDl yk j þ Re k

9 > > = 1  ; ð2q þ 1Þ> > 1 ; 2

Dlk

where kxk ¼ ju1 j þ a and kyk ¼ ju2 j þ a are the maximum wave velocities in the x- and y-directions, respectively, Re is the Reynolds number, Dlk is the diameter of control element Xk and p is the order of polynomial approximation. To increase the computational efficiency and to overcome the time step size restriction, we employ the local time-stepping (LTS) technique by computing local time steps Dtk on every element Xk of the computational mesh independently. Let wsk be the solution vector on element Xk at time t s ; s ¼ 1; 2; . . . In a triangular mesh, the element Xk is adjacent to three other elements Xk1 ; Xk2 ; Xk3 , each with corresponding solution vectors wsk11 ; wsk22 ; wsk33 . Due to the locality of the DG discretisation, the communication between the elements is performed through the numerical fluxes and penalisation terms on edges shared by these elements. Because of this, the solution vectors on the adjacent elements at the time t s are necessary for the evaluation of numerical fluxes on the three element edges. Let us assume that the following condition is satisfied during the whole computation

tsi 1 < ts < t si ;

i ¼ 1; 2; 3:

ð13Þ

Please cite this article in press as: O. Bublík et al., Comparison of discontinuous Galerkin time integration schemes for the solution of flow problems with deformable domains, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.03.063

O. Bublík et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx

5

Then, the solution vectors wski ; i ¼ 1; 2; 3 on the adjacent elements Xki can be interpolated with second-order accuracy at time ts by using solution vectors at times tsi and t si 1 as s

s 1

wski ¼ wi;1 wkii þ wi;2 wkii ;

ð14Þ

where wi;j ; j ¼ 1; 2 are the Lagrange interpolation polynomials defined as

wi;1 ¼

ts  tsi 1 ; t si  t si 1

wi;2 ¼

t s  t si : tsi 1  t si

The principle of the LTS algorithm based on the condition given by Eq. (13) is illustrated on a one-dimensional example in Fig. 2 for two consequent iteration steps. The figure shows three different elements Xk1 ; Xk and Xk2 , each with their own time step and maximum time level t s1 ; t s and ts2 , respectively. In the first example, Fig. 2(a), the element Xk is to be advanced in time due to ts < t s2 < ts1 . For this purpose, the solution vectors wski ; i ¼ 1; 2 of the two adjacent elements Xk1 and Xk2 are computed by means of the interpolation formula (14) and used for the solution update on the element Xk and the computation of a new time step. The following situation is shown in Fig. 2(b), where the comparison of the maximum time levels indicates the element Xk2 as the one that needs to be advanced in time. The computational algorithm of the explicit LTS method is implemented in Java programming language by employing the principles of object-oriented programming. 2.2.2. Implicit time integration method The implicit time integration of Eq. (11) is carried out using the incomplete interior penalty (IIP) discontinuous Galerkin scheme. After multiplying the system of Navier–Stokes equations (1) by a test function v and integrating them over a control element Xk , we apply the Green’s theorem and rewrite the viscous fluxes under the following linear form [13]

f vs ðw; rwÞ ¼

2 X @w K s;r ðwÞ ; @xr r¼1

s ¼ 1; 2

yielding the integral identity

d dt

Z

I 2 X @v ðf s ðwÞ  w U s Þ  dX þ ðF ðwþ ; w ; nÞ  w U n Þ  v dl @xs Xk Xk s¼1 @ Xk !  Z X I X I 2 2 2 X 2  X @w @v @w  ns  v dl þ cW K s;r ðwÞ dX þ K s;r ðwÞ ½w  v dl; ¼ @xr @xr @xs Xk s¼1 @ Xk s¼1 r¼1 @ Xk r¼1 w

v dX 

Z

ð15Þ

where cW is the penalty parameter. The implementation of boundary conditions is similar as in the case of the explicit method, following the papers [2,22]. Using similar principles as in the case of the LDG method, we arrive to a system of ordinary differential equations (11) time integrated with the second-order accurate implicit Crank–Nicolson scheme as

! ^ nÞ 1 nþ1 1 @RðWn ; U 1 n n 1 s ^ n Þ: Wnþ1 ¼ M  M W þ RðWn ; U s 2 @W 2 Dt Dt

ð16Þ

This discretisation approach leads to a system of linear equations Ax ¼ b that for the purpose of the present study is solved using the direct solver implemented in the UMFPACK library.

Fig. 2. Examples of the LTS iteration process.

Please cite this article in press as: O. Bublík et al., Comparison of discontinuous Galerkin time integration schemes for the solution of flow problems with deformable domains, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.03.063

6

O. Bublík et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx

Fig. 3. Example of mesh deformation as a result of airfoil movement. Here, the blending parameters are chosen as r 1 ¼ 0:5 and r 2 ¼ 2 with the centre at ½xc ; yc  ¼ ½0:5; 0.

3. Mesh deformation To carry out numerical simulations in deformable domains, it is necessary to select a suitable procedure to define a mapping between a fixed reference configuration and the real time-varying domain. Among the various known techniques described, e.g., in [11], we employ a continuous mapping in combination with a polynomial blending function. The principle of this simple and effective technique is illustrated in Fig. 3, where a triangular mesh around an airfoil is deformed as a result of airfoil rotation and movement. In accordance with [11,22,23], we first introduce a blending domain defined by two fixed circles centred at the airfoil, see Fig. 3(a). The first circle with radius r1 includes the moving airfoil and its boundary, while the second circle with radius r 2 closely adjoins the boundary of the computational domain. In this way, we obtain an annulus (r 1 6 r 6 r 2 ) that defines the part of the computational domain, mesh of which will be deformed in accordance with the movement of the airfoil. As the second step of the mapping procedure, the initial coordinates x0 of each mesh node are transformed into new coordinates y, Fig. 3(b), by using the following formula:

y ¼ p þ TðaÞx0 ;

ð17Þ

where p is the translation vector and TðaÞ is the rotation matrix. Finally, a new deformed mesh is generated according to the motion of the inner circle containing the airfoil, Fig. 3(c). The coordinates x of the new mesh points are computed by applying the blending function as

x ¼ bðrÞx0 þ ð1  bðrÞÞy;

ð18Þ

where the blending values are given as

8 0 if r 6 r 1 ; > > < if r P r 2 ; bðrÞ ¼ 1  > > : f rr1 if r1 < r < r2 : r2 r1

ð19Þ

For the purpose of this paper, the blending polynomial f occurring in Eq. (19) is chosen as a cubic polynomial in the form f ðxÞ ¼ 3x2  2x3 [11]. 4. Numerical results In this section, the developed computational algorithms based on the arbitrary Lagrangian–Eulerian (ALE) formulation of the compressible Navier–Stokes equations are used for the numerical solution of inviscid and viscous flow problems on moving domains. Note that all computations are carried out for linear polynomials (10) for N ¼ 3, regardless of the type of time integration method employed. 4.1. Isentropic vortex To verify the explicit LTS Runge–Kutta and implicit Crank–Nicolson solvers and to assess their convergence on time-varying domains, we consider an inviscid flow problem consisting of a compressible vortex in a square domain, Fig. 4 (top left). The two-dimensional model of isentropic vortex convection, described in detail in [24], is chosen because of the existence of exact solution [25]. At the beginning of the simulation ðt ¼ 0Þ, the isentropic vortex, also known as the Euler vortex [26], is located at the point ½x0 ; y0 , Fig. 4 (top right), and moves with the free-stream at an angle b with respect to the x-axis. In accordance with Please cite this article in press as: O. Bublík et al., Comparison of discontinuous Galerkin time integration schemes for the solution of flow problems with deformable domains, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.03.063

O. Bublík et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx

7

Fig. 4. Computational domain layout of the vortex benchmark (top left), initial state of the isentropic vortex (top right) and mesh deformation at the times t1 ¼ 0:25T (bottom left) and t2 ¼ 0:75T (bottom right).

the papers [11,25,26], the initial values of velocity vector components u1 ; u2 , density ., pressure p and total energy E are given as follows

   eðy  y0 Þ f ðx; yÞ ; u1 ¼ u1 cos b  exp 2p r c 2    eðx  x0 Þ f ðx; yÞ ; exp u2 ¼ u1 sin b þ 2p r c 2 1 !j1 e2 ðj  1ÞM 21 . ¼ .1 1  expðf ðx; yÞÞ ; 8p2 j !j1 e2 ðj  1ÞM21 expðf ðx; yÞÞ ; p ¼ p1 1  8p2 E¼

ð20Þ

p 1 þ . u21 þ u22 ; j1 2

where u1 ; .1 ; p1 and M 1 are the free-stream velocity, density, pressure and Mach number, j ¼ cp =cV is the adiabatic index and the function f ðx; yÞ is defined as

f ðx; yÞ ¼

 1 1  ðx  x0 Þ2  ðy  y0 Þ2 : 2 rc

The parameters rc and e occurring in equations (20) measure the size and strength of the moving vortex, respectively. For the computation of the free-stream values p1 ; .1 ; u1 , the following formulas are used

 j j  1 2 1j p1 ¼ p0i 1 þ M1 ; 2  11j j1 2 .1 ¼ .0i 1 þ M1 ; 2 sffiffiffiffiffiffiffiffiffiffi jp1 ; u1 ¼ M1

ð21Þ

.1

where the values of total pressure p0i and total density

.0i are chosen equal to 1.

Please cite this article in press as: O. Bublík et al., Comparison of discontinuous Galerkin time integration schemes for the solution of flow problems with deformable domains, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.03.063

8

O. Bublík et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx

As mentioned in [11,25], the exact solution of the isentropic vortex problem at ðx; y; tÞ is steady in the frame of reference moving with the free-stream. In other words, the analytical solution has the same shape as the initial conditions (20), with the exception that it propagates in the direction given by the angle b. The mesh deformation, which is shown in Fig. 4 (bottom) at two selected time instants, is prescribed in the form of time-dependent functions for the mesh point coordinates

xi;j ðtÞ ¼ xi;j ð0Þ þ A sin

      2p t 2p 2p sin xi;j ð0Þ sin yi;j ð0Þ ; t0 L L 

yi;j ðtÞ ¼ yi;j ð0Þ þ A sin

     2pt 2p 2p sin xi;j ð0Þ sin yi;j ð0Þ ; t0 L L

where L is the length of the square domain and A is the amplitude of mesh deformation. For the purpose of this paper, we consider a square domain of length L ¼ 20 and amplitude A ¼ 2 with free-stream boundary conditions and the initial position of the vortex centre set at the point ½x0 ; y0  ¼ ½5; 5, Fig. 4 (top right). For the simulation pffiffiffi parameters b ¼ p=4; M 1 ¼ 0:3; j ¼ 1:4; t 0 ¼ 2L=ð4 u1 Þ; r c ¼ 1:5 and e ¼ 0:8, the computations are performed for three qualitatively different meshes up to the time T ¼ 10=u1 when the vortex moves to the point ½xðTÞ; yðTÞ ¼ ½15; 15, i.e., to the opposite corner of the square domain. To assess the computational performance of our two solvers, we introduce the numerical error for the linear polynomial basis functions (10) in the following form:

error ¼

4 X

kwm 

wem kL2

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Z 4 X 2 ¼ ðwm  wem Þ dX;

m¼1

m¼1

ð22Þ

X

where k:kL2 denotes the L2 -norm and wm ; wem are the components of the vector of conservative variables w denoting the numerical and exact solutions, respectively. The error values computed with the explicit LTS and implicit methods using various CFL numbers and computational meshes are listed in Table 1. To approach the convergence analysis adequately, the numerical error and the corresponding CPU times are presented on four different mesh sizes h ¼ 0:5; 1; 2; 4, where h ¼ maxk fdiamðXk Þg. The table future shows convergence rate using the half-step method

r¼



1 log10 ð2Þ

log10

 errorðh ¼ 1Þ ; errorðh ¼ 0:5Þ

ð23Þ

where errorðh ¼ 1Þ and errorðh ¼ 0:5Þ are the numerical errors computed for meshes of element size 1 and 0.5, respectively. The results summarised in Table 1 indicate that the convergence rate r of both the explicit LTS and implicit methods is strongly dependent on the chosen CFL number and the element size h. For example, it can be noted that with increasing CFL values the convergence rate decreases and is lower than the formal second order of accuracy. To compare the efficiency of both the explicit and implicit methods, it is necessary to choose appropriate CFL numbers, which will guarantee that both methods will lead to approximately the same computational errors. For this purpose, we utilise the graphs shown in Fig. 5, displaying the dependence of the numerical error on the resulting CPU time. On the basis of the graphs, an optimal CFL number can be determined for a pre-selected numerical error. For example, if a numerical error of 0.06 is chosen, then Fig. 5 indicates that the optimal choice of CFL number is 0.385 in the case of the explicit LTS method and 2 in the case of the implicit method. Further from Table 1, it can be noted that the CPU times for the element size h ¼ 0:5 and the aforementioned CFL numbers are 96.6 s and 1102.4 s for the explicit LTS and implicit methods, respectively. By means of this test example, it can be demonstrated that for inviscid flow simulation, the second order two step explicit LTS Runge– Kutta method is about ten times faster than the implicit Crank–Nicolson method.

Table 1 The numerical error, CPU times and convergence rates for linear polynomial basis functions computed with the explicit LTS and implicit methods. Element size

h¼4

h¼2

h¼1

h ¼ 0:5

Method

CFL

Error

CPU (s)

Error

CPU (s)

Error

CPU (s)

Error

CPU (s)

r

Explicit LTS

0.07 0.385 0.7

0.1517 0.1517 0.1518

5.3 1.7 1.2

0.0961 0.1002 0.1061

16.2 3.8 2.8

0.0541 0.0700 0.0922

80.2 16.2 10.0

0.0194 0.0564 0.0929

508 96.6 52.6

1.4802 0.3102 0.0100

Implicit

1 2 5 10 50 100

0.1564 0.1566 0.1576 0.1580 0.1580 0.1587

39.4 22.3 12.6 8.8 5.1 3.7

0.1056 0.1096 0.1159 0.1212 0.1237 0.1247

102.1 57.4 27.3 19.1 9.5 7.1

0.0699 0.0791 0.0958 0.1093 0.1249 0.1229

408.3 209.4 93.1 55.1 21.9 17.4

0.0361 0.0489 0.0716 0.0925 0.1224 0.1234

2189.5 1102.4 462.7 251.9 80.7 59.1

0.9533 0.6938 0.4201 0.2408 0.0292 0.0059

Please cite this article in press as: O. Bublík et al., Comparison of discontinuous Galerkin time integration schemes for the solution of flow problems with deformable domains, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.03.063

O. Bublík et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx

9

Fig. 5. Numerical error as a function of CPU time for various CFL numbers – explicit LTS method (left) and implicit method (right).

Fig. 6. Position parameters of the airfoil (CG = centre of gravity).

Fig. 7. Computational domain layout of the airfoil problem.

4.2. Oscillating NACA 0012 airfoil To assess the accuracy and efficiency of the two DG time integration schemes in a more complex setting, we consider a viscous flow problem consisting of an oscillating NACA 0012 airfoil of length L ¼ 1. The oscillation of the airfoil is described by two superimposed movements, as illustrated in Fig. 6: 1. vertical motion of the centre of gravity (CG):

xðtÞ ¼ 0:5;

yðtÞ ¼ 0:3 sin t;

Please cite this article in press as: O. Bublík et al., Comparison of discontinuous Galerkin time integration schemes for the solution of flow problems with deformable domains, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.03.063

10

O. Bublík et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx

2. clockwise rotation around CG:

#ðtÞ ¼



1 p cos t 1  et : 6

The oscillating airfoil is set in a dimensionless rectangular domain of size 25-by-20 that is discretised into 14 344 triangular elements, see Fig. 8 for a detailed view at the near-airfoil region. The numerical solution of the dimensionless compressible Navier–Stokes equations in ALE form is carried out for the free-stream Mach number M 1 ¼ 0:3, the Reynolds number

Fig. 8. A detailed view at the structure of the computational mesh around the airfoil.

Fig. 9. Detailed view at the distribution of Mach number around the oscillating airfoil at time T ¼ 4p.

Table 2 Comparison of CPU times after reaching the simulation time T ¼ 4p. Time integration method

CPU time (s)

Explicit LTS (CFL ¼ 0:385) Implicit (CFL ¼ 2)

2028 110 035

Please cite this article in press as: O. Bublík et al., Comparison of discontinuous Galerkin time integration schemes for the solution of flow problems with deformable domains, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.03.063

O. Bublík et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx

11

Fig. 10. The y-component of force acting on the airfoil surface as a function of time computed with the explicit LTS and implicit methods.

Re ¼ 1000 and the zero angle of attack. The boundary conditions are implemented in a similar way as in the vortex benchmark (compare Fig. 4 (top left) and Fig. 7). The values of total pressure p0i and total density .0i are chosen equal to 1, whereas the static pressure po ¼ p1 at the outlet is computed using Eq. (21). With regard to the results obtained in the vortex benchmark discussed above, the CFL numbers for the current airfoil problem are chosen as 0.385 and 2 for the explicit LTS and implicit methods, respectively. This choice should guarantee that both methods will achieve a similar numerical error, making their comparison possible. CPU times corresponding to the two considered time integration methods after reaching the simulation time T ¼ 4p are listed in Table 2. The distribution of the Mach number at the aforementioned time instant is shown in Fig. 9, while Fig. 10 displays the time evolution of the y-component of force acting on the airfoil surface. On the basis of the relatively similar isolines and forces computed by both methods, it can be concluded that the results indicate approximately the same numerical error, supporting so our previous observation. By comparing the CPU times from Table 2, it can be noted that in the case of a more complex geometry and laminar viscous flow, the second order two step explicit LTS Runge–Kutta method is approximately 50 times faster than the implicit Crank–Nicolson method. 5. Conclusion In this work, the compressible Navier–Stokes equations in ALE form were solved by means of the discontinuous Galerkin (DG) finite element method formulated for unstructured triangular meshes. For the mapping between a fixed reference configuration and the real time-varying domain, we have used the continuous approach based on the application of a polynomial blending function that compared to other methods does not require the solution of additional equations and is, thus, computationally effective. To advance the numerical solution in time, the explicit Runge–Kutta scheme improved in the sense of the local time-stepping (LTS) technique and the implicit Crank–Nicolson scheme were employed and their computational performance analysed in terms of numerical error and CPU time usage. First, the proposed DG solvers were verified against the analytical solution of the Euler vortex in a time-varying domain and then applied for the numerical simulation of viscous flow around an oscillating NACA 0012 airfoil. On the basis of the vortex benchmark, optimal CFL numbers for the explicit LTS and implicit methods were determined in such a way as to ensure a similar numerical error. Then, using this setting of CFL numbers, the airfoil problem was solved and analysed. The obtained results confirmed that both time integration methods were able to achieve qualitatively comparable flow fields. A comparison of CPU time usage further revealed that the explicit LTS method is more efficient than the implicit one for both inviscid and laminar viscous flows. Finally, the present study showed that the explicit method improved by the LTS technique can be seen as a simple method compared to the implicit one, and as such can be taken as an alternative to the implicit method in terms of computational efficiency. The conclusions drawn in this study may be also applicable for the solution of coupled fluid–structure interaction problems [27,28]. In the future, the computational efficiency of explicit LTS schemes will be further studied for turbulent flow problems. Acknowledgements This work was supported by the European Regional Development Fund (ERDF), project ‘‘NTIS – New Technologies for the Information Society’’, European Centre of Excellence, CZ.1.05/1.1.00/02.0090 and by the project TE01020068 ‘‘Centre of research and experimental development of reliable energy production’’ of the Technology Agency of the Czech Republic. Please cite this article in press as: O. Bublík et al., Comparison of discontinuous Galerkin time integration schemes for the solution of flow problems with deformable domains, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.03.063

12

O. Bublík et al. / Applied Mathematics and Computation xxx (2015) xxx–xxx

References [1] W.H. Reed, T.R. Hill, Triangular mesh methods for the neutron transport equations, Technical Report LA-UR-73-479, Los Alamos Scientific Laboratory, 1973. [2] B. Cockburn, C.-W. Shu, Runge–Kutta discontinuous Galerkin methods for convection-dominated problems, J. Sci. Comput. 16 (2001) 173–261. [3] F. Bassi, S. Rebay, G. Mariotti, S. Pedinotti, M. Savini, A high-order accurate discontinuous finite element method for inviscid and viscous turbomachinery flows, in: Second European Conference on Turbomachinery, Fluid Dynamics and Thermodynamics, p. 99108. [4] M. Feistauer, V. Kucˇera, J. Prokopová, Discontinuous Galerkin solution of compressible flow in time-dependent domains, Math. Comput. Simul. 80 (2010) 1612–1623. [5] J. Hesthaven, T. Warburton, Nodal high-order methods on unstructured grids. I. Time-domain solution of Maxwell’s equations, J. Comput. Phys. 181 (2002) 186–221. [6] J. Hesthaven, T. Warburton, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, Springer, 2008. [7] L.D. Angulo, J. Alvarez, F. Teixeira, M. Pantoja, S. Garcia, Causal-path local time-stepping in the discontinuous Galerkin method for Maxwell’s equations, J. Comput. Phys. 256 (2014) 678–695. [8] L.D. Angulo, J. Alvarez, F. Teixeira, M. Pantoja, S. Garcia, A discontinuous Galerkin method for ideal two-fluid plasma equations, Commun. Comput. Phys. 9 (2011) 240–268. [9] P. Birken, G. Gassner, M. Haas, C.-D. Munz, Efficient time integration for discontinuous Galerkin method for the unsteady 3D Navier–Stokes equations, in: European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2012), pp. 4334–4353. [10] F. Renac, S. Gérald, C. Marmignon, F. Coquel, Fast time implicit–explicit discontinuous Galerkin method for the compressible Navier–Stokes equations, J. Comput. Phys. 251 (2013) 272–291. [11] P.-O. Persson, J. Bonet, J. Peraire, Discontinuous Galerkin solution of the Navier–Stokes equations on deformable domains, Comput. Methods Appl. Mech. 198 (2009) 1585–1595. [12] B. Cockburn, C.-W. Shu, The local discontinuous Galerkin method for time-dependent convection–diffusion systems, SIAM J. Numer. Anal. 35 (1998) 2440–2463. [13] V. Dolejší, M. Holík, J. Hozman, Efficient solution strategy for the semi-implicit discontinuous Galerkin discretization of the Navier–Stokes equations, J. Comput. Phys. 230 (2011) 4176–4200. [14] D.N. Arnold, An interior penalty finite element method with discontinuous elements, SIAM J. Numer. Anal. 19 (1982) 742–760. [15] R. Hartmann, P. Houston, Symmetric interior penalty DG methods for the compressible Navier–Stokes equations I: Method formulation, Int. J. Numer. Anal. Model. 3 (2006) 1–20. [16] R. Hartmann, P. Houston, Symmetric interior penalty DG methods for the compressible Navier–Stokes equations II: Goal-oriented a posteriori error estimation, Int. J. Numer. Anal. Model. 3 (2006) 141–162. [17] J. Douglas, T. Dupont, Interior penalty procedures for elliptic and parabolic Galerkin methods, Comput. Methods Appl. Sci., Lect. Notes Phys. 58 (1976) 207–216. [18] M.F. Wheeler, An elliptic collocation-finite element method with interior penalties, SIAM J. Numer. Anal. 15 (1978) 152–161. [19] S. Sherwin, G.E. Karniadakis, A new triangular and tetrahedral basis for high-order (hp) finite element methods, Int. J. Numer. Methods Eng. 38 (1995) 3775–3802. [20] M. Dubiner, Spectral methods on triangles and other domains, J. Sci. Comput. 6 (1991) 345–390. [21] E.F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer, 1999. ˇ esenek, M. Feistauer, A. Kosík, DGFEM for the analysis of airfoil vibrations induced by compressible flow, ZAMM – J. Appl. Math. Mech. 93 (2013) [22] J. C 387–402. [23] M. Feistauer, J. Horácˇek, M. Ruzˇicˇka, P. Svácˇek, Interaction of a flexibly supported airfoil and a channel flow, Eng. Mech. 15 (2008) 57–77. [24] G. Erlebacher, M. Hussaini, C.-W. Shu, Interaction of a shock with a longitudinal vortex, J. Fluids Mech. 337 (1997) 129–153. [25] E. Petrini, G. Efraimsson, J. Nordström, A numerical study of the introduction and propagation of a 2-D vortex, Technical Report FFA TN 1998-66, The Aeronautical Research Institute of Sweden, 1998. [26] C. Yee, N.D. Sandham, M.J. Djomehri, Low dissipative high order shock-capturing methods using characteristic-based filters, J. Comput. Phys. 150 (1999) 199–238. [27] M. Feistauer, J. Hasnedlová-Prokopová, J. Horácˇek, A. Kosík, V. Kucˇera, DGFEM for dynamical systems describing interaction of compressible fluid and structures, J. Comput. Appl. Math. 254 (2013) 17–30. [28] B. Froehle, P.-O. Persson, A high-order discontinuous Galerkin method for fluid–structure interaction with efficient implicit–explicit time stepping, J. Comput. Phys. 272 (2014) 455–470.

Please cite this article in press as: O. Bublík et al., Comparison of discontinuous Galerkin time integration schemes for the solution of flow problems with deformable domains, Appl. Math. Comput. (2015), http://dx.doi.org/10.1016/j.amc.2015.03.063