Second-order-accurate explicit fluctuation splitting schemes for unsteady problems

Computers & Fluids 38 (2009) 1384–1393

Contents lists available at ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

Second-order-accurate explicit fluctuation splitting schemes for unsteady problems G. Rossiello, P. De Palma, G. Pascazio, M. Napolitano * Dipartimento di Ingegneria Meccanica e Gestionale, Centro di Eccellenza in Meccanica Computazionale, Politecnico di Bari, Via Re David 200, 70125 Bari, Italy

a r t i c l e

i n f o

Article history: Received 7 May 2007 Accepted 23 January 2008 Available online 12 April 2008

a b s t r a c t This paper is concerned with the issue of obtaining explicit fluctuation splitting schemes which achieve second-order accuracy in both space and time on an arbitrary unstructured triangular mesh. A theoretical analysis demonstrates that, for a linear reconstruction of the solution, mass lumping does not diminish the accuracy of the scheme provided that a Galerkin space discretization is employed. Thus, two explicit fluctuation splitting schemes are devised which are second-order accurate in both space and time, namely, the well known Lax–Wendroff scheme and a Lax–Wendroff-type scheme using a three-pointbackward discretization of the time derivative. A thorough mesh-refinement study verifies the theoretical order of accuracy of the two schemes on meshes with increasing levels of nonuniformity. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction The residual distribution or fluctuation splitting (FS) method has been developed at first for computing compressible steady flows [1–6]. The FS approach is based on a tessellation of the computational domain and on a continuous reconstruction of the solution over linear cell-vertex (triangular/tetrahedral) elements. The method is based on three fundamental steps: (i) evaluating the fluctuation, namely, the flux balance over each computational cell; (ii) distributing the residual contributions (signals) among the vertices of the cell using suitable coefficients; (iii) updating the solution at each node by summing up all contributions from the elements sharing that node. More recently, the FS method has been extended to the case of unsteady hyperbolic problems [7–12]. The first and most logical choice was the FS Lax–Wendroff scheme [7], which however was never proved to be second-order-accurate on a general triangular mesh. Thus several approaches have been pursued by different research groups to design second-order-accurate FS schemes. In particular März and Degrez have shown that a Runge–Kutta (or method-of-lines) approach is unsuitable to develop second-order-accurate FS schemes; thus they have reformulated the FS spatial distribution as a Petrov–Galerkin Finite Element method and have shown that a consistent mass matrix can be obtained, leading to an implicit scheme which is second-order-accurate in both space and time [8]. More recently Csík et al. [9] proposed an implicit approach in the framework of space–time residual distribution schemes, by introducing linear tetrahedral elements in two dimensions, which again involves a consistent mass matrix. In order to render such an approach less cumbersome and applicable to three dimensions, Mezine et al. [10] have replaced the tetrahedral space–time elements with hexahedral ones, * Corresponding author. E-mail address: [email protected] (M. Napolitano). 0045-7930/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2008.01.021

which use the same spatial mesh at different time levels. In his doctoral thesis [11], Caraeni proposed an alternative approach based upon a second-order-accurate finite-difference time discretization of FS schemes, leading to a different consistent mass matrix. Recently, the authors derived the sufficient conditions to obtain consistent mass matrices, and proposed an alternative implicit scheme based on a second-order-accurate finite-difference time discretization [12]. Then, they provided a stable third-orderaccurate FS scheme, in which the space accuracy is achieved by reconstructing the gradient of the solution at the vertices of each element and the time one is obtained by an implicit scheme employing a four-point-backward finite-difference formula [13]; and the sufficient conditions for an FS scheme to be ðr þ 1Þthorder-accurate in the case of unsteady problems were derived for the case of a finite-difference discretization of the time derivative as well as for that of a space–time approach [13]. The common feature of such unsteady flow FS schemes enjoying an order of accuracy greater than one is the use of a non-diagonal consistent mass matrix for discretizing the physical time derivative. On the other hand, to the authors knowledge, a rigorous and complete analysis of the sufficient conditions required by an explicit FS scheme (using a lumped mass matrix) to achieve second-order accuracy in both space and time has never been performed and is given here, as described in the following. After briefly recalling the FS methodology for both steady and unsteady problems, a truncation error analysis provides the sufficient conditions for an FS scheme to be second-order-accurate in both space and time, leading to a class of mass matrices satisfying such conditions. Then, it is shown that mass lumping preserves second-order accuracy provided that a Galerkin space discretization is employed. Finally, two explicit FS schemes are derived from implicit ones with lumped mass matrix, by Taylor expanding the spatial residual, which are shown to achieve second-order accuracy in both space and time: the Lax–Wendroff (LW) scheme, and a LW-type

1385

G. Rossiello et al. / Computers & Fluids 38 (2009) 1384–1393

scheme using a three-point-backward approximation of the time derivative. It is indeed a fortunate coincidence, as well as an additional pleasure for the authors, that such a finding appears in honour of Prof. Lerat, who has contributed significantly to the development of LW-type schemes. 2. The FS method In this section, the fundamentals of FS schemes are provided with reference to a scalar conservation equation on a 2D spatial domain ou þ $  FðuÞ ¼ 0; ot with u : H ! R; H ¼ X  ½0; þ1½; X # R2 ;

F ¼ ðf ðuÞ; gðuÞÞs : ð1Þ

For linear advection, the flux vector can be expressed as F ¼ ku, where k ¼ ða; bÞs is the velocity. For a divergence-free k, Eq. (1) governs the advection of the scalar quantity u. Eq. (1) is discretized in space and time: the spatial computational domain, X, is divided into cell-vertex triangular elements, T, with linear dimension h; the generic node and time-level are labeled i and n, respectively, Dt being the time step. The integral and weak forms of Eq. (1) read:  Z  ou þ $  F dX 8C # X ð2Þ ot C

þ1 0

Z

  ou x þ $  F dX dt; ot R2

ð3Þ

respectively, x being a continuous function with compact support in R2  ½0; 1½. Either such formulation can be taken as the starting point for developing an FS method, according to the following three steps: 1. Approximate the sought solution, u, with a polynomial reconstruction, uh , based on nodal values. 2. Set up a system of equations by integrating Eq. (2) on a dual cell or, equivalently, Eq. (3). 3. Solve the system and update the solution at the nodes. Given a triangular element T: i, j, and k indicate its three vertices, numbered counterclockwise, and ni is the vector normal to the edge opposite to node i, pointing towards the interior of T and having the edge’s length: ni ¼ ðyj  yk Þ^i  ðxj  xk Þ^j; and j ¼ modði þ 1; 3Þ;

for

The median dual cell, C i , associated with the node i is sketched in Fig. 1; its area, Si , is equal to one-third of the sum of the areas of the elements sharing node i: Si ¼

X ST : 3 T3i

The spatial variation of the solution is assumed to belong to the functional space of linear polynomials on linear (P1) elements, i.e.:

and Z

Fig. 1. Median dual cell of node i.

i ¼ 1; 2; 3;

k ¼ modði þ 2; 3Þ;

where ^i and ^j are the Cartesian unit vectors. By definition, one has P3 i¼1 ni ¼ 0, so that, defining the inflow parameter as:

uh ðx; y; tÞ ¼

3 X

N j ðx; yÞuj ðtÞ

8x ¼ ðx; yÞ 2 T;

ð5Þ

j¼1

with local base and gradient on the element equal to: N i ðxÞ ¼ 1 þ $N i ðxÞ ¼

ðx  xi Þ  ni ; 2ST

ni 2ST ;

ð6Þ

so that $uh ¼

3 1 X uj nj ; 2ST j¼1

ð7Þ

and 3 X o h ouj ðtÞ u ðx; y; tÞ ¼ N j ðx; yÞ: ot ot j¼1

The flux balance over a triangle T is called fluctuation or residual; in the case of an advection equation and linear reconstruction, it takes the following form: Z /T ¼ $  F h dS ¼ k  $uh ST ; ð8Þ T

ki ¼

1 k  ni ; 2

ð4Þ

k being the averaged advection velocity over the element, one has P3 i¼1 ki ¼ 0. A given edge is said to be outflow (respectively, inflow) when the inflow parameter relative to the opposite node is smaller (respectively, larger) than zero. The area ST of a triangle T can be computed as: ST ¼

1 kni  nj k 2

for any i; j ¼ 1; 2; 3 with i–j:

where, for non-linear advection, k is an averaged value of the velocity which maintains conservation at the discrete level. As a consequence of the linear reconstruction, using Eqs. (4) and (7), the fluctuation over a triangle can be computed exactly, as: /T ¼

3 X

kj uj :

ð9Þ

j¼1

In order to evaluate the residual at any given node (i.e., at each dual cell around node i), each signal /Ti , namely, the portion of the ele-

1386

G. Rossiello et al. / Computers & Fluids 38 (2009) 1384–1393

ment fluctuation /T which is sent from the element to node i, is computed as: /Ti ¼ bTi /T ;

ð10Þ

bTi being the distribution coefficients. By definition, the following conservation constraint holds: 3 X

/Tj ¼ /T )

j¼1

3 X

bTj ¼ 1:

ð11Þ

integrals of both the unsteady and steady terms are employed for computing the fluctuations and bounded distribution functions are used. In general, the signals need to be computed implicitly, as shown in [12], which analyzes the conditions which make a mass-matrix consistent and thus the scheme second-order-accurate, and defines a class of consistent mass matrices. Here, such an analysis is recalled as the starting point to obtain second-order-accurate explicit schemes, by a suitable mass lumping.

j¼1

The final equation, for the steady case, reads: X T /i ¼ 0:

3.1. Truncation error analysis ð12Þ

T3i

An explicit update derived by a two-level explicit Euler integration scheme is used, uinþ1

Dt X T ¼ uni  / ; Si T3i i

ð13Þ

as the simplest smoother to drive the spatial residual to zero and thus to obtain the steady state solution. In such a case, the scheme is uniquely defined by its distribution coefficients, namely, its signals. Distributing the fluctuation using bounded coefficients, bTi , a Linearity Preserving ðLPÞ scheme is obtained, namely, a scheme which preserves an initial exact linear solution, thus being second-order-accurate in space for homogeneous advection equations [2]. Any such scheme can be recast into the form of a Petrov–Galerkin finite element method, with the weight function xi obtained by adding an elemental term to the classical Galerkin shape function, Ni ,   1 xTi ðx; yÞ ¼ N Ti ðx; yÞ þ bTi  ; ð14Þ 3 as: /Ti ¼ ¼

Z T

Z

T

xTi k  $uh dX N Ti k  $uh dX þ

 Z  1 bTi  k  $uh dX ¼ bTi /T : 3 T

ð15Þ

In particular, the Lax–Wendroff [2] and UCV [14] schemes, of interest here, are given as:   1 ki Dt 2 ST bT;LW ¼ ; with Dt ¼ þ min i 3 2ST 3 T3i maxj2T jkj j and bT;UCV ¼ i

1 2 ki : þ P 3 3 3j¼1 jkj j

3. Second-order accuracy for unsteady problems When the unsteady problem is considered, time accuracy becomes a fundamental issue. It has been recognized in [8] that a Runge–Kutta (or method-of-lines) approach is unsuitable to develop second-order-accurate FS schemes. This is related to the consistent discretization of the unsteady fluctuation defined by the signal: Z ouh dX: ð16Þ xi wTi ¼ ot T It has been shown in [13] that, for the case of the unsteady advec2 tion equation, the truncation error is Oðh Þ if: (i) a second-orderaccurate space reconstruction is used for the numerical flux and the solution; (ii) a second-order-accurate finite-difference discreti3 zation in time is employed for ðouh =otÞi ; and iii) /Ti þ wTi ¼ Oðh Þ. The last condition is fulfilled if suitable approximations of the space

Following a space–time approach, the conservation law, Eq. (1), can be written as: ou þ $  FðuÞ ¼ $t  GðuÞ ¼ LðuÞ ¼ 0; ot

ð17Þ

where $t is the space–time divergence and G ¼ ðu; FÞs . The corresponding weak formulation reads: Z

xi Lh ðuh ðuj ÞÞ dH ¼

Z 0

H

þ1

Z

xi Lh ðuh ðuj ÞÞ dX dt ¼

X

X

UEi ðuj Þ ¼ 0;

E3i

ð18Þ where xi indicates a generic weight function, E and UEi indicate the space–time elements and signals, respectively, and h is the linear measure of E. Furthermore, uh is a polynomial representation based on the computed nodal values, uj , with j 2 E. In [13], it has been shown that sufficient conditions for the truncation error of the rþ1 scheme to be Oðh Þ are: (i) the use of an ðr þ 1Þth-order-accurate rþdþ1 h Þ, d being the number of approximation G for G; (ii) UEi ¼ Oðh space dimensions. The latter condition on the signals is fulfilled if a suitable approximation of the space–time integrals is employed for computing the fluctuation and bounded distribution functions are used. Alternatively, using a finite-difference approximation in time based on a multi-step scheme, the weak formulation of Eq. (17) reads: Z xi X

 h  X ou ðuj Þ UTi ðuj Þ ¼ 0; þ $  F h ðuh ðuj ÞÞ dX ¼ ot T3i

ð19Þ

where both terms in the integral must be computed at the same time level, leading to an implicit scheme. In this case, the conditions rþ1 for the truncation error to be Oðh Þ are [13]: (i) an ðr þ 1Þth-orderaccurate space reconstruction is used for F h and uh ; (ii) an ðr þ 1Þthorder-accurate finite-difference discretization in time is employed rþd for ðou=otÞi ; and (iii) UTi ¼ Oðh Þ. The last condition on the signals is fulfilled if suitable approximations of the space integrals of both the unsteady and steady terms are employed for computing the fluctuations, and bounded distribution functions are used. Integrating Eqs. (18) or (19), an expression for the signals is obtained, involving in general a mass matrix, m, and thus leading to implicit schemes. Only this second approach will be considered here, for brevity. In order to obtain a second-order-accurate scheme, linear reconstructions for uh and F h are employed. Concerning the signals, one has:  h  ou þ $  F h ðuh Þ dX ot T  h  Z Z ou ou ¼  xi xi ð$  F h ðuh ÞÞ  $  FðuÞÞÞ dX dX þ ot ot T T

UTi ¼

Z

xi

3

¼ Oðh Þ;

ð20Þ

provided that the time derivative is approximated with secondorder accuracy, namely,

1387

G. Rossiello et al. / Computers & Fluids 38 (2009) 1384–1393

ouh ou 2  ¼ OðDt2 Þ ¼ Oðh Þ; ot ot

ð21Þ

T

ouh dX ¼ ot ¼

Z xi

3 X

Nj

T

j¼1

3 X

ouj ; mTij ot

j¼1

3 X ouj ouj dX ¼ ot ot j¼1

xi N j dX

T

T

xi dX ¼ bTi ST :

ð32Þ

3 X

mTij ¼

j¼1

3 X

cij bTi ST ;

ð33Þ

j¼1

so that the following constraints hold: 3 X

cij ¼ 1;

for

i ¼ 1; 2; 3:

ð34Þ

3 X

ð23Þ

3 Z X j¼1

xj

T

ouh dX ¼ ot

Z T

3 ouh ST X ouj dX ¼ : ot 3 j¼1 ot

ð35Þ

ð24Þ 3 Z X i¼1

mTij

ouj þ bTi /T ; ot

ð25Þ

whereas the semi-discrete equation at point i can be written as: " #   3 X X X T T T T ouj ð26Þ þ bi / ¼ 0: Ui ¼ Ui ¼ mij ot T3i T3i j¼1 The integral at the right hand side of Eq. (24) can be evaluated as Z Z mTij ¼ xi N j dX ¼ N j ðxH xi dX; ð27Þ i Þ T

T

xi

3 X 3 3 X 3 X ouh ouj X ouj dX ¼ ¼ : mTij cij bTi ST ot ot ot i¼1 j¼1 i¼1 j¼1

ð36Þ

Comparing the right hand sides of the last two equations, the following conservation constraints are obtained: 3 X i¼1

cij bTi ¼

1 ; 3

for

j ¼ 1; 2; 3:

ð37Þ

The second-order-accurate implicit scheme proposed in [12] fulfills the constraints (34) and (37), and is based on a three-level discretization of the time derivative and a dual-time-stepping technique [15]. Each step in the fictitious time can be written as:

T

xH i

uinþ1;kþ1 ¼ uinþ1;k 

2 T.

Proof 3.1. Since N j ðxÞ is bounded for x 2 T: xi ðxÞ inf N j ðxÞ 6 xi ðxÞN j ðxÞ 6 xi ðxÞ sup N j ðxÞ x2T

8x 2 T;

ð28Þ

x2T

thus, inf N j ðxÞ

due to the linear reconstruction employed for u, one has:

On the other hand, using Eqs. (23) and (31),

j¼1

where

xj ¼ 1;

j¼1

are the components of the element mass matrix. Thus, the general form of the signals for an implicit scheme reads:

Z

xi ðxÞ dX 6 T

Z

N j ðxÞxi ðxÞ dX T Z 6 sup N j ðxÞ xi ðxÞ dX x2T

ð29Þ

T

which, being N j continuous in T, implies that Z Z 9 xH such that N j ðxÞxi ðxÞ dX ¼ N j ðxH xi ðxÞ dX; i 2 T i Þ T

ð30Þ

T

which ends the proof. h Recalling the definition of bTi , Eq. (22), and indicating N j ðxH i Þ ¼ cij , one has: mTij ¼ cij bTi ST :

ð31Þ

The constraints to be satisfied by the mass-matrix coefficients for the scheme to be consistent and conservative are now recalled [12]. Since 3 X

Z

Moreover, using Eq. (31), one has

3 X

T

T

x2T

j¼1

xi N j dX ¼

Furthermore, since for conservation

Z

where Z mTij ¼ xi N j dX

UTi ¼

3 Z X

j¼1

and xi

mTij ¼

j¼1

and the weight functions, xi , are bounded. Starting from Eq. (20), recalling that a linear reconstruction for F h is employed, so that $  F h ðuh Þ is constant over the element T, and using the linear Lagrangian shape functions, N j ; j ¼ 1; 2; 3, one has: Z Z 1 xi $  F h ðuh Þ dX ¼ bTi /T with bTi ¼ xi dX ð22Þ ST T T

Z

3 X

N j ¼ 1;

j¼1

it follows that

Ds nþ1;k U Si i

where the time derivative at the nodes is evaluated as:   3unþ1  4uni þ uin1 ou : ¼ i 2Dt ot i

ð38Þ

ð39Þ

The choice of the mass matrix was performed so as to guarantee that, for an upwind space discretization scheme, also the unsteady term is distributed only between the downstream nodes, see [12], for details, leading to the following form: 2 T 3 b1 ð2  bT1 Þ bT1 ð1  bT2 Þ bT1 ð1  bT3 Þ ST 6 T T T T T T T 7 ð40Þ mij ¼ 4 b2 ð1  b1 Þ b2 ð2  b2 Þ b2 ð1  b3 Þ 5: 3 bT3 ð1  bT1 Þ bT3 ð1  bT2 Þ bT3 ð2  bT3 Þ We now investigate whether a suitable mass lumping can render such schemes explicit, while maintaining their second-order accuracy in both space and time. Mass lumping is obtained by taking cij ¼ dij , dij being the Kronecker delta function. Such a choice satisfies the constraints (34) and the corresponding diagonal mass matrix reads: mTij ¼ dij bTi ST :

ð41Þ

Therefore, in order to satisfy constraints (37), a unique choice of the distribution coefficients exists, namely, bTi ¼ 1=3, for i ¼ 1; 2; 3. Such coefficients are obtained using a Galerkin scheme, namely, employing xi ¼ N i . In conclusion, for linear elements, when using a Galerkin space integration, the mass-lumping procedure maintains the second-order accuracy of the method.

1388

G. Rossiello et al. / Computers & Fluids 38 (2009) 1384–1393

butions of the internal edges cancel out; using Eq. (6), the second term is discretized as:

3.2. Second-order-accurate FS schemes with mass lumping 3.2.1. Two-time-level schemes The general space–time conservation statement at a generic node i, during a single time step Dt, in which Galerkin weight functions are employed, leads to the following discrete equation: Uinþ1 ¼

XZ T3i

tþDt t

Z T

 h  ou þ $  F h ðuh Þ dX dt ¼ 0 Ni ot

ð42Þ

and then, integrating in time,  Z X Z Uinþ1 ¼ N i ðuh;nþ1  uh;n Þ dX þ Dt n N i $  F h;nþ1=2 ðuh Þ dX T3i

T

Dt 2

n Z  k  ni Dt ki T;n k  $uh dX ¼ / : 2ST 2 ST T

ð52Þ

Therefore, the final scheme reads: X 1 Dt ki  T;n / ¼ 0; Uinþ1 ¼ Si ðuinþ1  uni Þ þ Dt þ 3 2 ST T3i

ð53Þ

which is the explicit Lax–Wendroff scheme, named here FSLW.

T

¼ 0: ð43Þ The first term represents the unsteady residual, which, as shown in the section above, can be integrated by mass lumping without loosing second-order accuracy, obtaining XZ N i ðuh;nþ1  uh;n Þ dX ¼ Si ðuinþ1  uni Þ: ð44Þ T3i

0

0.25

0.5

0.75

1

1

1

0.75

0.75

0.5

0.5

0.25

0.25

T

The second term can be discretized linearizing the flux function in time Z 1 N i $  F h;nþ1=2 ðuh Þ dX ¼ /T;nþ1=2 : ð45Þ 3 T Thus, a Crank-Nicolson implicit Galerkin scheme is obtained: Uinþ1

¼

Si ðunþ1 i



uni Þ

1 X T;nþ1=2 þ Dt / ¼ 0; 3 T3i

ð46Þ

with /T;nþ1=2

0

 nþ1  u þ un ; ¼/ 2

0 0

ð47Þ

/ðunþ1 Þ þ /ðun Þ : 2

$  F ¼ k  $u; with constant advection k, for simplicity, one has:  n oðk  $uÞ Dt ðk  $uÞnþ1=2 ¼ ðk  $uÞn þ þ OðDt2 Þ ot 2 Dt ¼ ðk  $uÞn  ½k  $ðk  $uÞn þ OðDt 2 Þ: 2

and Dt 2

T

1

0.25

0.5

0.75

1

1

1

0.75

0.75

0.5

0.5

0.25

0.25

Z

Dt N i k  ½k  $un n d‘ 2 oT Z Dt ½k  $N i k  $un dX; þ 2 T

N i ½k  $ðk  $uÞn dX ¼ 

0

ð49Þ

Then, integrating in space the two terms at the right hand side, one has: Z 1 N i ðk  $uÞn dX ¼ /T;n ; ð50Þ 3 T



0.75

ð48Þ

In order to obtain an explicit scheme, a Taylor expansion of the convective (implicit) part of the residual is performed. Considering the case

Z

0.5

Fig. 2. Coarsest mesh A (h1 ¼ 3:108  102 ).

or, alternatively, /T;nþ1=2 ¼

0.25

ð51Þ

respectively. The first terms of Eq. (51) sum to zero when the signals to node i are collected, since N i has compact support and the contri-

0

0 0

0.25

0.5

0.75

1 2

Fig. 3. Coarsest mesh B (h1 ¼ 3:534  10 ).

1389

G. Rossiello et al. / Computers & Fluids 38 (2009) 1384–1393

0

0.25

0.5

0.75

which can be considered a LW-type scheme using a three-pointbackward time derivative, named here FSLW2.

1

1

1

4. Results

0.75

0.75

0.5

0.5

0.25

0.25

0

0 0

0.25

0.5

0.75

1 2

Fig. 4. Coarsest mesh C (h1 ¼ 3:125  10 ).

3.2.2. Three-time-level schemes Using a finite-difference approximation in time based on three levels, Eq. (39), the weak formulation of the problem reads:  h nþ1 XZ ou h h þ $  F ¼ N ðu Þ dX Unþ1 i i ot T T3i  nþ1  3ui  4uni þ uin1 1 X T;nþ1 þ / ¼ 0: ð54Þ ¼ Si 2Dt 3 T3i Starting from such a scheme, an explicit scheme can be derived by a Taylor series expansion of the convective term in Eq. (54), namely:   nþ1  X 3ui  4uni þ uin1 1 ki /T;n ¼ 0; þ þ Dt 2Dt ST 3 T3i

ð55Þ

-0.5

0

-1

-0.5

-1.5

-1

-2

-1.5

LOG [L∞]

LOG [L1]

Unþ1 ¼ Si i

In this section, the accuracy of the schemes has been verified numerically by solving two unsteady problems on the same periodic domain ½0; 1  ½0; 1 up to t ¼ 1. In order to assess the performance of the efficient explicit schemes, FSLW and FSLW2, the very accurate but costly implicit one of [13], with consistent mass matrix and UCV distribution coefficients (CN-UCV), has also been considered for comparison. It is noteworthy that the two explicit schemes directly derive by the corresponding implicit ones, given by Eqs. (46) and (54). Thus, such schemes have also been used to solve the test cases presented here. However, they have proven to be as accurate as their explicit counterparts in spite of their implicit nature, which renders their (dual-time-stepping) solution procedure much dearer. Therefore, such results have been omitted here, for the sake of brevity. In order to provide an undisputable numerical validation of the order of accuracy of the schemes, three sets of meshes have been employed for discretizing the computational domain, called A, B, and C. Each coarsest mesh has 32 cells along each edge of the computational domain: the coarsest mesh A is a regular unstructured mesh obtained using a Delaunay triangulation, see Fig. 2; the coarsest mesh B is an unstructured mesh, obtained by a triangulation with a non uniform discretization of the edges of the computational square, which is made up of regular triangles having largely different sizes, see Fig. 3; the coarsest mesh C is a considerably skewed mesh obtained by a perturbed structured module, shown in Fig. 4. Starting from each coarsest mesh, four refined meshes were obtained by mapping it onto each portion of the computational domain obtained by dividing it into 4, 16, 64, and 256 pffiffiffiffiffiffi identical squares, (hi ¼ h1 =4i for i ¼ 1; . . . ; 4), with h1 ¼ 1= N T , N T being the total number of triangles of the corresponding coarsest mesh. It is noteworthy that such a refinement does not preserve the shape of each cell, locally, thus providing a very severe test for the order of accuracy of the schemes. For all schemes, the time step has been computed as

-2.5

-3

-2

-2.5

-3

-3.5 FS LW FS LW2 CN UCV

-4

FS LW FS LW2 CN UCV

-3.5

-4

-4.5 -3

-2.5

-2

LOG [h]

-1.5

-3

-2.5

-2

-1.5

LOG [h]

Fig. 5. Linear advection: L1 norms (left) and L1 norms (right) of the errors for mesh A.

G. Rossiello et al. / Computers & Fluids 38 (2009) 1384–1393

-0.5

0

-1

-0.5

-1.5

-1

-2

-1.5

LOG [L∞]

LOG [L1]

1390

-2.5

-3

-2

-2.5

-3

-3.5 FS LW FS LW2 CN UCV

-4

FS LW FS LW2 CN UCV

-3.5

-4

-4.5 -3

-2.5

-2

-3

-1.5

-2.5

-2

-1.5

LOG [h]

LOG [h]

-0.5

0

-1

-0.5

-1.5

-1

-2

-1.5

LOG [L∞]

LOG [L1]

Fig. 6. Linear advection: L1 norms (left) and L1 norms (right) of the errors for mesh B.

-2.5

-3

-2

-2.5

-3

-3.5 FS LW FS LW2 CN UCV

-4

FS LW FS LW2 CN UCV

-3.5

-4

-4.5 -3

-2.5

-2

-3

-1.5

-2.5

-2

-1.5

LOG [h]

LOG [h]

Fig. 7. Linear advection: L1 norms (left) and L1 norms (right) of the errors for mesh C.

Dt ¼

  2 ST ; min 3 T3i maxj2T jkj j

which guarantees the stability for both the FSLW and FSLW2 schemes. The accuracy of the schemes has been verified at first by computing the unsteady linear advection of the double-sine-shaped function, 2

The logarithms of the L1 and L1 norms of the errors are reported in Figs. 5–7, which confirm the predicted order of accuracy of the schemes. A second test case has then been considered, namely, the circular advection of the double-sine-shaped function, ( uðx; yÞ ¼

2

uðx; yÞ ¼ sin ð2pxÞ sin ð2pyÞ; with k ¼ ð1; 2Þ.

with

2

2

sin ð2prÞ sin ð2hÞ

for r 6 0:5

0

for r > 0:5

;

1391

-0.5

0

-1

-0.5

-1.5

-1

-2

-1.5

LOG [L∞]

LOG [L1]

G. Rossiello et al. / Computers & Fluids 38 (2009) 1384–1393

-2.5

-3

-2

-2.5

-3

-3.5 FS LW FS LW2 CN UCV

-4

FS LW FS LW2 CN UCV

-3.5

-4

-4.5 -3

-2.5

-2

-3

-1.5

-2.5

-2

-1.5

LOG [h]

LOG [h]

-0.5

0

-1

-0.5

-1.5

-1

-2

-1.5

LOG [L∞]

LOG [L1]

Fig. 8. Circular advection: L1 norms (left) and L1 norms (right) of the errors for mesh A.

-2.5

-3

-2

-2.5

-3.5

-3 FS LW FS LW2 CN UCV

-4

-4.5

FS LW FS LW2 CN UCV

-3.5

-4

-3.5

-3

-2.5

-2

-1.5

LOG [h]

-3.5

-3

-2.5

-2

-1.5

LOG [h]

Fig. 9. Circular advection: L1 norms (left) and L1 norms (right) of the errors for mesh B.

r 2 ¼ ðx  0:5Þ2 þ ðy  0:5Þ2 ;

1

h ¼ sin

  y  0:5 ; r

and k ¼ ð2py; 2pxÞ. The logarithms of the L1 and L1 norms of the errors are reported in Figs. 8–10, which again confirm the second-order accuracy of the schemes. For such a test case, two further mesh refinements ði ¼ 5; 6Þ were performed on the mesh type B for the FSLW and 2 FSLW2, to verify that all error norms tend to zero as h , for h ! 0.

Finally, Figs. 11–13 provide the contour levels for the circular advection test case at t ¼ 1, obtained using mesh A with h2 ¼ h1 =4 by the FSLW, FSLW2, and CN-UCV schemes, respectively. The FSLW and FSLW2 schemes provide similar solutions, which are affected by a lagging phase error. The CN-UCV scheme solution coincides within plotting accuracy with the exact one, except for the small oscillations, which are expected for any second-orderaccurate linear scheme.

G. Rossiello et al. / Computers & Fluids 38 (2009) 1384–1393

-0.5

0

-1

-0.5

-1.5

-1

-2

-1.5

LOG [L∞]

LOG [L1]

1392

-2.5

-3

-2.5

-3

-3.5 FS LW FS LW2 CN UCV

-4

FS LW FS LW2 CN UCV

-3.5

-4

-4.5 -3

-2

-2.5

-2

-3

-1.5

-2.5

-2

-1.5

LOG [h]

LOG [h]

Fig. 10. Circular advection: L1 norms (left) and L1 norms (right) of the errors for mesh C.

0

0.25

0.5

0.75

1

1

0 1

0.25

0.5

0.75

1

1

1

0.75

0.75

0.75

0.75

0.5

0.5

0.5

0.5

0.25

0.25

0.25

0.25

0

0 0

0.25

0.5

0.75

1

0

0 0

0.25

0.5

0.75

1

Fig. 11. Circular advection: contours levels computed using the FSLW scheme and mesh A with h2 ¼ h1 =4 (Du ¼ 0:1, minimum plotted level is zero).

Fig. 12. Circular advection: contours levels computed using the FSLW2 scheme and mesh A with h2 ¼ h1 =4 (Du ¼ 0:1, minimum plotted level is zero).

From all of the results shown in Figs. 5–13 it clearly appears that the implicit scheme of [12] provides a smaller error in all cases. However it is much more expensive than the two explicit schemes (see [12] for details). In the present calculations, using a dual-time-stepping procedure and reducing the L1 norm of the residual to 1010 at each physical time step, it has a CPU cost about 20 times greater than those of the explicit ones. Needless to say, such a comparison is somewhat unfair to the implicit scheme, the convergence criterion for the pseudo-time convergence being very demanding.

5. Conclusions This paper addresses the issue of obtaining explicit (masslumped) fluctuation splitting schemes for unsteady flow problems, which achieve second-order accuracy in both space and time on arbitrary unstructured meshes. The general sufficient conditions required by a scheme to be second-order-accurate in space and time are provided and fluctuation splitting Lax–Wendroff schemes are shown to satisfy such conditions in spite of their lumped mass discretization of the time derivative. A very thorough mesh refine-

G. Rossiello et al. / Computers & Fluids 38 (2009) 1384–1393

0

0.25

0.5

0.75

References

1

1

1

0.75

0.75

0.5

0.5

0.25

0.25

0

0 0

0.25

0.5

0.75

1393

1

Fig. 13. Circular advection: contours levels computed using the CN-UCV scheme and mesh A with h2 ¼ h1 =4 (Du ¼ 0:1, minimum plotted level is zero).

ment study, using three sets of unstructured meshes, provides a numerical validation of such a result, which the authors are very pleased to dedicate to their friend, Prof. Lerat. Acknowledgement This work has been supported by MiUR and Politecnico di Bari, Grant Cofinlab 2000.

[1] Struijs R, Deconinck H, Roe P. Fluctuation splitting schemes for the 2D Euler equations. In: VKI LS 1991-01; 1991. [2] Paillere H. Multidimensional upwind residual distribution schemes for the Euler and Navier–Stokes equations on unstructured grids, PhD thesis. von Karman Institute; 1995. [3] van der Weide E, Deconinck H, Issmann E, Degrez G. Fluctuation splitting schemes for multidimensional convection problems: an alternative to finite volume and finite element methods. Comput Mech 1999;2(23):199–208. [4] De Palma P, Pascazio G, Napolitano M. A hybrid fluctuation splitting scheme for two-dimensional compressible steady flows. In: Hafez MM, Chattot JJ, editors. Innovative methods for numerical solution of partial differential equations. New Jersey: World Scientific Publishing; 2002. p. 305–33. [5] Catalano LA, De Palma P, Pascazio G, Napolitano M. Cell-vertex adaptive Euler method for cascade flows. AIAA J 1995;33(12):2299–304. [6] Bonfiglioli A, De Palma P, Pascazio G, Napolitano M. An implicit fluctuation splitting scheme for turbomachinery flows. ASME J Turbomach 2005;127:395–401. [7] Hubbard ME, Roe PL. Compact high resolution algorithms for time dependent advection problems on unstructured grids. Int J Numer Meth Fluid 2000;33:711–36. [8] März J, Degrez G. Improving time accuracy for residual distribution schemes. VKI-DC report 1996-17; 1996. [9] Csík A, Ricchiuto M, Deconinck H, Poedts S. Space–time residual distribution schemes for hyperbolic conservation laws. AIAA Paper 2001-2617. [10] Mezine M, Ricchiuto M, Abgrall R, Deconinck H. Monotone and stable residual distribution schemes on prismatic space–time elements for unsteady conservation laws. In: von Karman Institute Lecture Series 2003-05; 2003. [11] Caraeni DA. Development of a multidimensional residual distribution solver for large eddy simulation of industrial turbulent flows, PhD thesis. Lund Institute of Technology; 2000. [12] De Palma P, Pascazio G, Rossiello G, Napolitano M. A second-order-accurate monotone implicit fluctuation splitting scheme for unsteady problems. J Comput Phys 2005;208:1–33. [13] Rossiello G, De Palma P, Pascazio G, Napolitano M. Third-order-accurate fluctuation splitting schemes for unsteady hyperbolic problems. J Comput Phys 2007;222:332–52. [14] Giles M, Anderson W, Roberts T. Upwind control volume, a new upwind approach. AIAA Paper 90-0104. [15] Jameson A. Time dependent calculations using multigrid with applications to unsteady flows past airfoils and wings. AIAA Paper 91-1596.