adaptive mesh refinement - Semantic Scholar

European Congress on Computational Methods in Applied Sciences and Engineering ECCOMAS 2000 Barcelona, 11-14 September 2000

c ECCOMAS

ADAPTIVE MESH REFINEMENT: A WAVELET POINT OF VIEW Olivier Roussel?y and Marc P. Errera? ? Oce National

d'Etudes et de Recherches Aerospatiales (ONERA), DSNA/MSDH, BP 72, 92322 Ch^atillon Cedex, France E-mail: [email protected], [email protected], Web page: http://www.onera.fr/ yLaboratoire

d'Informatique pour la Mecanique et les Sciences de l'Ingenieur (LIMSI), Universite Paris-Sud, 91405 Orsay Cedex, France E-mail: [email protected], Web page: http://www.limsi.fr/

Key words: Adaptive mesh re nement, Wavelet, Multiresolution, Computational uid

dynamics

Abstract. In this paper, it is shown that adaptive mesh re nement (AMR) methods

applied to partial dierential equations in nite volume formulation are equivalent to a weak formulation of these equations on a particular basis of test-functions: the Haar wavelet basis. This equivalence enables us to give a new meaning to some re nement sensors used in AMR methods in order to predict the regions that need to be re ned. More precisely, it appears that local truncation error estimators and physical sensors such as gradients are no equivalent predictors, but both necessary. Moreover, new terms must be added to the usual re nement sensors. All of these terms are viewed in a 2D laminar Navier-Stokes test-case: the backward-facing step ow.

1

Olivier Roussel and Marc P. Errera

1 INTRODUCTION

Local mesh re nement methods appeared in the 1970s in order to locally improve solutions of partial dierential equations (PDE). Literature on these methods show many possible re nement strategies and sensors. Adaptive strategies can be divided into two classes: those which aim at re ning a unique mesh and those which aim at dispatching re nement levels into dierent subdomains or subscales. In the rst class, several methods have been developed in the framework of unstructured meshes (Batina et al [4], HarranKlotz [19], Lohner et al [24], Marcum et al [25], Powell and Ashford [29]). Re ning by adding grid points to a unique mesh is actually well-adapted to unstructured meshes. Other works propose such a strategy with structured meshes (Aftosmis [1], Powell and De Zeeuw [28], Schonfeld and Rudgyard [33]), but such a re ned mesh must necessarily be treated as a generalized unstructured mesh. In the second class, some works used an unstructured mesh (Braaten and Connell [10], Mavrilis [26]), but most are based on embedded structured meshes. Historically, the rst methods used topologically similar, lined up subdomains - the so-called Multi-Level Adaptative Technique (MLAT) (Brandt [11], McCormick et al [27], Angot et al [2]). Then, authors turn to topologically dierent, lined up subdomains - the so-called Adaptive Mesh Re nement (AMR) method (Berger [7], Berger and Colella [8], Quirk [30], Jouhaud and Borrel [22]) - but also to rotating subdomains (Berger [7], Benoit [6]). As for adaptive strategies, literature shows many dierent re nement sensors. They are usually divided into two classes: numerical sensors, which usually measure the local truncation error (Berger [7], Benoit [6]), and physical sensors, which are based on physical quantities of the ow (Hentschel and Hirschel [21], Coirier and Powell [13], Quirk [30]). At the same time, several works in dierent scienti c elds on scale-space decomposition have led to wavelet theory (see Daubechies [14]). In computational uid dynamics (CFD), these methods have been used to analyze or simulate turbulent ows (Farge [17], Schneider et al [32]). The wavelet approach was also used to re ne solutions of conservation laws in a nite volume formulation. This led to the so-called multiresolution schemes (Harten [20], Bihari et al [9], Chiavassa and Donat [12]). Thus, multiresolution schemes and adaptive mesh re nement methods reach the same goal, although the paths are dierent. In this paper, the comparison between the wavelet approach and adaptive mesh re nement enables us to give a new meaning to re nement sensors used in AMR methods. Particularly, numerical sensors and physical sensors do not belong to dierent classes.

2 ADAPTIVE MESH REFINEMENT AND WAVELET APPROACH 2.1 The Haar wavelet basis

The Haar wavelet [18] is the most elementary wavelet and historically the oldest (1910). It is de ned by two basic functions : the analyzing wavelet , usually called the mother wavelet, and the scaling function ', usually called the father wavelet. 2


The scaling function ' = '0;0 on [0; 1] is de ned as (

x 2 [0; 1] '(x) = 10 ifelsewhere The analyzing wavelet = 0;0 on [0; 1] is de ned as 8 > > > > > > > < (x) = > > > > > > > :

h

?1 if x 2 0; 12

(1) h

1

h i if x 2 21 ; 1

0

elsewhere

(2)

The translations and dilatations of the mother wavelet constitute an orthonormal basis of L2([0; 1]) (Haar [18], Daubechies [14]).

2j x ? i

j;i (x) = 2 2 j

(3)

where j > 0 and 0 i < 2j . ψ ϕ

1

00

0

1 ψ

00

1 0

1

−1 2

x

x

ψ1 1

ψ1 0

0

1

x

- 2

Figure 1: Haar wavelet basis,primary elements on [0; 1]

3


Thus, each function f 2 L2([0; 1]) can be expressed as

f (x) = c0;0 '0;0(x) + where c0;0 =< f; '0;0 >[0;1] and dj;i =< f;

j ?1 1 2X X

j =0 i=0

dj;i

j;i (x)

(4)

j;i >[0;1].

Note that the Haar mother wavelet does not have good time-frequency localization: its Fourier transform ^() decays like jj?1 for ! 1. Nevertheless, its wavelet basis is an ecient tool for approximating discontinuous L2 functions.

2.2 Finite volume formulation

Let us consider a conservation law in one space dimension 8 @Q + @ F (Q; @ Q) = 0 > > > < x @t @x > > > : Q(x; 0) = Q (x) 0

(5)

Q(x; t) is a vector of q components, and F denotes the total ux of Q. It can be split into an inviscid contribution Fi and a viscous contribution Fv following F (Q; @xQ) = Fi (Q) + Fv (@xQ) We consider a coarse grid G = f(xk = k x); k = 0; : : : ; Nxg of a computational domain

= [0; L]. Therefore, we have Nx x = L. For each cell !k = [xk ; xk+1], 0 k Nx ? 1, we de ne a characteristic function k as k (x) =

(

1 if x 2 !k 0 elsewhere

(6)

The nite volume formulation of Eq. (5) on the grid G is in fact a classical weak formulation of this same equation on the basis of test-functions k . Z @Q Z @F k k (E ) @t (x; t) (x) dx + @x (x; t) k (x) dx = 0 (7)

ie (E k )

Z

!k

@Q dx + Z @F dx = 0 @t ! @x k

4

(8)


2.3 Adaptive mesh re nement

We then decide to re ne a cell wk on one scale. We split this cell and create two re ned cells !1k;0 and !1k;1. In the same manner, we de ne two characteristic functions k1;0 and k1;1 on each cell. The nite volume formulation on each subcell is (i = 0; 1) Z @F Z @Q k k dx + @x k1;i dx = 0 (9) (E1;i)

@t 1;i where k1;i(x) = xk

(

1 if x 2 !1k;i 0 elsewhere

(10)

xk+1

ωk

(j=0) ω k1 0

ω k1 1 (j=1)

ω k2 0

ω k2 1

ω k2 2

ω k2 3 (j=2)

Figure 2: Subcells of !

k

For a given scale j , we can de ne the same nite volume formulation in each subcell (0 i < 2j ) Z @Q Z @F k k k dx = 0 (Ej;i) dx + (11) j;i

@t

@x j;i where

(

x 2 !j;ik kj;i(x) = 10 ifelsewhere

(12)

This system of equations is redundant. There are in fact more equations than degrees of freedom, given that Z @Q Z @Q Z @Q k dx = k dx + k dx (13)

@t j +1;2i

@t j +1;2i+1

@t j;i 5


2.4 Wavelet approach

Let us de ne a Haar wavelet basis f'k0;0 ; j;ik j j > 0; 0 j < 2j g in each cell !k . 8 > k (x) = p 1 if x 2 ! k ; 0 elsewhere > > ' > 0;0 > > xk < > > > > > > :

k j;i (x)

k 2j x?xxk

= p2 k x j

2

?i

!

(14)

where j > 0, 0 i < 2j and xk = xk+1 ? xk . We can easily show that 8 > > > > > > < > > > > > > :

k 'k0;0 = p k x

= p2 k kj+1;2i+1 ? kj+1;2i x j

2

k j;i

(15)

Therefore, applying a similar linear combination to the system of equations (Ej;ik ), 8 ~ k p 1 Ek > > > E < > > > : E ~j;ik

xk

p2 2 j

xk

h

Ejk+1;2i+1 ? Ejk+1;2i

i

(16)

we get the non-redundant system 8 > > > >
Z @F Z @Q > > > k dx + k dx = 0 k) : (E~j;i

@x j;i

@t j;i (E~ k )

(17)

Thus, the adaptive mesh re nement method applied to a conservation law in a nite volume formulation is equivalent to writing this same equation in a classical weak formulation, but replacing the characteristic functions on each cell by a Haar wavelet basis (Roussel [31]). Moreover, in the rst case, we obtain more equations than degrees of freedom whereas, in the second case, the number of equations is equal to the number of degrees of freedom. Nevertheless, the nite volume formulation is incomplete and approximate numerical schemes have to be added. Therefore, computations made on two dierent scales and reported on the coarsest grid are not equal. 6


3 CONSEQUENCES ON REFINEMENT SENSORS In this section, denotes the standard L2 dot product on . 8f; g 2 L2 ( ); < f; g >=

Z

3.1 Computational improvement

f (x) g(x) dx

Let us now consider a cell !j;ik and its two subcells !jk+1;2i and !jk+1;2i+1. At time tn = n t, Qkj;i(x) denotes the restriction of the approximate solution Q on !j;ik , and Q kj;i denotes its mean value. Therefore, we have

Qkj;i(x) = Q kj;i kj;i(x)

(18)

On the j + 1-th scale, the approximate solution becomes

Q kj;i(x) = Q kj+1;2i kj+1;2i(x) + Q kj+1;2i+1 kj+1;2i+1(x)

(19)

Given that kj;i = kj+1;2i + kj+1;2i+1, we can decompose Q kj;i(x) into its mean value and the rest, following 8 Q kj+1;2i + Q kj+1;2i+1 Q kj+1;2i ? Q kj+1;2i+1 > > k > + = Q > > < j +1;2i 2 2 > > k k k k > > > kj+1;2i+1 = Qj+1;2i + Qj+1;2i+1 ? Qj+1;2i ? Qj+1;2i+1 : Q 2 2 we get Q k + Q k Q k ? Q k Q kj;i(x) = j+1;2i 2 j+1;2i+1 kj;i(x) + j+1;2i+12 j+1;2i kj+1;2i+1 ? kj+1;2i (20)

Given that

< Q;

Z

k j;i > =

Q

k j;i dx

p

k k Q ? Q k = x j+1;2i+1 j+1;2i 2 22 j

we deduce from Eqs. (15) and (20) that Q k + Q k Q kj;i(x) = j+1;2i 2 j+1;2i+1 kj;i(x) + < Q; 7

k j;i >

k j;i

(21)

(22)


The improvement on the estimation of Q between the j -th and the j + 1-th scale is

D Qkj;i = Q kj;i(x) ? Qkj;i(x) = kj;i kj;i + < Q;

k j;i >

k j;i

(23)

where

Q kj+1;2i + Q kj+1;2i+1 k ? Qj;i 2 Thus, the improvement between two scales can be decomposed into a corrective term on the j -th scale and an additive term which is nothing more than the next term in the development of Q in the Haar orthonormal wavelet basis. Therefore, we can consider that the computation is suciently accurate locally on a given scale if the improvement is small. In the other case, the computation needs to be improved and one more scale is necessary. Let us denote kj;i =

2 +1 = p2 k < Q; j;ik >= Q kj+1;2i+1 ? Q kj+1;2i (24) x The term kj;i is nothing more than the local truncation error. The term j;ik is proportional to the mean gradient in the x-direction. Thus, both local truncation error and mean gradients must be small to ensure that the computation is suciently accurate. In most papers that deal with local mesh re nement methods [23], local truncation error and gradients are considered as dierent kinds of sensors. The rst is often called a numerical sensor, whereas the second is considered as a physical sensor. They are in fact complementary, and both are necessary to determine the regions that should be re ned.

j;ik

j

3.2 Non-dimensional sensors

In most AMR methods, the sensors are divided by their maximum on the ow eld [23]. The main advantage is to show values between 0 and 1. Unfortunately, such a choice does not allow the evolution of the sensor with the scale number to be shown. Actually, there is always a point where the sensor value is 1. Moreover, if the sensor value is low but homogeneous in the ow eld, its value divided by the maximum will be close to 1 almost everywhere. In this article, the terms kj;i and j;ik are chosen because their dimension is the same as Q. We can therefore compute non-dimensional values of the sensors following 8 > > > > > > < > > > > > > :

kj;i

kj;i =Q max

j;ik k j;i = Q max 8

(25)


We can then chose a unique sensibility " for all components of Q and specify that the computation in !j;ik is acurate enough if

jkj;i j < " and j kj;ij < " 3.3 Extension to two dimensions

(26)

Having developed a one-dimensional orthonormal basis, we would like to use these functions as building blocks in higher dimensions. One way of doing so is to take the tensor product of two one-dimensional bases [14] and to de ne k k k j;ix ;iy (x; y ) = j;ix (x) j;iy (y )

(27)

Considering 'k0;0 = ?k 1;0, the resulting functions constitute an orthonormal wavelet basis (j ?1 ; 0 i < 2j if j 0, i = 0 elsewhere). Observe that here the scale parameter j simultaneously controls the dilatation in x and in y. We now consider the domain = [0; Lx] [0; Ly ], a two-dimensional cell !k of , 0 k (Nx ? 1)(Ny ? 1), and its subcells on the j -th scale !j;ik ;i , where j > 0 and 0 ix;y < 2j . The improvement between the estimations of Q on the j -th and the j +1-th makes four terms appear x y

8 > > > > > > > > > > > > > < > > > > > > > > > > > > > :

Q kj+1;2i ;2i + Q kj+1;2i +1;2i + Q kj+1;2i ;2i +1 + Q kj+1;2i +1;2i +1 k ? Qj;i 4 j;ik = Q kj+1;2i +1;2i +1 ? Q kj+1;2i ;2i +1 + Q kj+1;2i +1;2i ? Q kj+1;2i ;2i +1 kj;i =

x

x

y

x

y

y

x

x

x

y

y

x

y

y

x

y

j;ik = Q kj+1;2i +1;2i +1 ? Q kj+1;2i +1;2i + Q kj+1;2i ;2i +1 ? Q kj+1;2i ;2i +1 j;ik = Q kj+1;2i +1;2i +1 ? Q kj+1;2i +1;2i ? Q kj+1;2i ;2i +1 + Q kj+1;2i ;2i +1 x

y

x

y

x

y

x

y

x

y

x

y

x

y

x

y

(28)

kj;i is the local truncation error ; j;ik is proportional to the mean value of @x Q ; j;ik is proportional to the mean value of @y Q ; j;ik is proportional to the mean value of @xy2 Q. Our interest will focus on the last term. Actually, it is not a standard sensor like the other ones. In the numerical results, we will re ne the mesh according to the values of this particular sensor. 9


4 NUMERICAL SCHEME

A computer code, known as the MSD code, has been developed since the mid 1980s. This code can compute turbulent and reactive ows about realistic aerospace con gurations and is widely used in a great variety of scienti c and engineering problems.

4.1 Governing equations

The governing equations are the two-dimensional compressible Navier-Stokes equations. The system of equations in curvilinear coordinates (; ) in the strong conservation law form can be written as @ Q + @ (F + F ) + @ (G + G ) = S (29) @ J @ i v @ i v where Q = (; u; v; E )t J is the inverse cell volume, and E are the density and total energy per unit volume, respectively. (u; v) denotes the velocity eld. The ux vectors Fi and Gi denote the inviscid uxes, Fv and Gv denote the viscous ones.

4.2 Numerical methodology

Details of the schemes implemented in the MSD code have been extensively reported (see e.g. Dutoya and Errera [15]). Only a brief review is given here. Using an implicit time discretization and after linearizing the uxes, the discretized equations may be written in the following form :

fNUMERICS g Q = PHY SICS

(30)

where NUMERICS represents the implicit operator containing the ux Jacobians, PHY SICS is the residual, and Q denotes the change in the solution vector at each time step. Accurate approximations are made on the governing equations in the right-hand side. The role of the left-hand side is to convey the solution in a stable manner. Eq. (29) is used to compute the residual with an upwind ux-dierence splitting method, whereas viscous uxes are evaluated with a second-order accurate centered scheme. A brief description of the implementation of the scheme is presented. The state of the

ow at each interface is described by two vectors of variables on each side, QL and QR. The superscripts L, R stand for the values from the left and right sides of the cell face of the computational volume. They are obtained by upwind extrapolation to the cell surface with a MUSCL approach using the variables stored at the cell centers in conjunction with a limiter. The order of accuracy for the upwinding is determined by the accuracy of the extrapolated left and right states. 10


The variable can be the conservative vector Q or the primitive variables. In the latter case, the numerical procedure restores Q at the end of the computation. The interface variables are obtained with a hybrid upwind scheme up to third-order of accuracy. In the nite-volume procedure, the dierences among all the numerical schemes lie mainly in the de nition of the numerical ux evaluated at the cell interface. In this study, the ux is computed in the following way i h i h Fi+ 21 = 12 F QR + F QL + 21 F QR ? F QL (31) where = sign (A) and A = @F @Q Introducing the diagonalizing matrices T and T ?1 and the diagonal matrix of eigenvalues A, is de ned as

= T sign () T ?1

(32)

This expression requires the computation of all the right and left eigenvectors. It is much more ecient to make use of the formal decomposition of the Jacobian (Dutoya and Errera [16]), which represents an interesting and powerful property, so that the Jacobian matrix may be rewritten in the following way :

(33) A = u n I + c r+ l + ? r ? l ? where c denotes the speed of sound, I the identity matrix and r+ , l+, r? , l? are the right and left eigenvectors of A corresponding to + = un + c and ? = un ? c, respectively.

The expressions of the eigenvectors are given in [16]. They are constructed as a sum of a thermodynamic vector and a kinematic vector. It can be shown that any operator of the Jacobian matrix can be written as i h i h (A) = 0 + + ? 0 r+ l+ + ? ? 0 r? l?

(34)

Expressions (33) and (34) neither depend on the number of variables, nor on the thermodynamic properties of the uid. Eq. (34) shows that the Jacobian matrix and all derived operators can be written with only three eigenvalues and two tensor products. From a practical point of view, all functions of the Jacobian can be directly expressed, and the explicit computation of the eigenvectors corresponding to 0 = un can be avoided. Consequently, the implementation of most numerical schemes, such as ux-vector splitting,

ux-dierence splitting or any other type using upwind operators, can now be made in a straightforward manner. 11


Conditions at the various boudaries of the ow domain are implemented in an implicit manner. At the in ow boundary, one quantity is extrapolated from the interior and all the others are prescribed. At the exit, only one quantity, a constant pressure, is speci ed, while the others are extrapolated from the interior. The implicit boundary conditions are applied by appropriate modi cation of the block matrices. To further increase the rate of convergence to steady state, the method uses local time steps.

4.3 Implicit method

The time integration is obtained by an implicit factored method (Beam and Warming [5]) . The left-hand-side of Eq. (29) is split into three one-dimensional operators. This operation introduces an O (t2) factorization error, which aects the convergence properties of the method. However, use of an ADI scheme results in a stable scheme for small to moderate time steps. Moreover, the ADI procedure involves solving a sequence of easily invertible tri-diagonal systems. This is eciently accomplished with little computational storage. The main disadvantage of the factorization technique lies in their relatively slow rate of convergence to a steady state.

4.4 Adaptive algorithm

Adaptive algorithm is made through directional cell division. This methodology is

exible and ecient for resolving ner length scales. It avoids grid redistributions and minimizes computational cost since cells are only added in zones of the ow eld that require further enhancement. The construction of new cells takes advantage of the coordinate directions available and is one-dimensional in each coordinate direction. Moreover, this technique requires a ow solver which is able to divide an arbitrary collection of cells without unstructured data storage. The MSD code employs a multidomain approach and oers great exiblity in handling complex geometries. Each geometry is meshed separately with structured grids and then overset onto a main grid in an arbitrary manner. A grid is de ned as a simple collection of structured (hexahedral) cells not necessarily linked together. In the general organization establishing the structure of the relationships among the grids, all the re ned cells of the given level of re nement form a new domain. In this paper, a coarse grid is generated in the primarily ow region and another ner grid is overset on the major grid by doubling the grid density in each direction. This con guration allows the direct computation of the wavelet coecients. Two types of comunication are possible : 1. the solutions in each grid are completely independent 2. the ow variables between the grids are exchanged via a source term in the coarser grid. This is not a direct injection of ne-grid quantities onto the coarse grid. The communication is achieved thanks to a source term in the coarse grid. 12


In the interest of brevity, the implementation of the zonal boudary conditions is not discussed here. Details may be found in [15]. The question of conservative vs. nonconservative interpolation is beyond the scope of the present study. It will certainly be a topic of future research.

5 NUMERICAL RESULTS

The laminar ow over a backward-facing step is an interesting test-case for several reasons:

the velocity eld shows at least one zone of recirculation ; the laminar state enables us to study the re nement quality without any interaction

with a turbulence model ; literature shows several experimental results.

5.1 Geometry and notations

Geometry and notations are given in Fig. 3. Let us denote:

h and l the height and length of the upwind channel ; H and L the height and length after enlargement ; S =H ?h ; u0 the upstream ow velocity. x’B

x"B B

u max

h

H S

y

O

A

C xA

x

l

L

x’C

x"C

Figure 3: Geometry and notations

The values of these parameters are close to the ones used by Armaly et al [3] in their experimental results: 13

Olivier Roussel and Marc P. Errera 8 > > > > > > < > > > > > > :

h H l L u0

= = = = =

5 10?3 m 10 10?3 m 0:2 m 0:5 m 10 m:s?1

(35)

The uid is air, considered as a perfect monospecies gas with the following thermodynamic properties: 8 > > > < > > > :

= M = Cp = Pr =

1:177 kg:m?3 29 10?3 kg 1003:4 J:kg?1 :K ?1 0:7

(36)

As the upstream ow is very slow, the pressure forces are very weak. This is apt to damage the computational accuracy. It was therefore decided to use a moderate velocity and increase the viscosity to reach the desired Reynolds number. This one is de ned by a characteristic velocity u~ = 32 umax and a characteristic length D = 2h. h = 400 Re = 43 umax (37) Following Armaly et al [3], the backward-facing step ow at such a Reynolds number is incompressible, laminar and two-dimensional. It shows a zone of recirculation in the A-region and a second one in the B-region slightly appears. 0.1 0.05 0

-0.1

0

0.1

0.2

0.3

0.4

0.5

-0.1

0

0.1

0.2

0.3

0.4

0.5

0.01

0.005

0

Figure 4: Mesh with correct scale (top) and stretched 10 times in y (bottom)

14


The coarse mesh is a 70 20 grid. A coarser 35 10 grid is generated to compute the sensors on two grids (Fig. 4).

5.2 Sensor plotting

The sensors k0;0, k0;0, k0;0 and k0;0 are quadri-dimensional vectors like Q. We will focus on the second component, which corresponds to u. The main ow is in the x direction. Moreover, there are regions where the ow is almost parallel, particularly before the step and in the outgoing ow region. The second component of the four sensors is plotted in Fig. 5. 0.01

0.005

0

-0.1

0

0.1

0.2

0.3

0.4

0.5

0.01

0.005

0

0.1 0.01 0.001 0.0001

-0.1

0

0.1

0.2

0.3

0.4

0.5

-0.1

0

0.1

0.2

0.3

0.4

0.5

-0.1

0

0.1

0.2

0.3

0.4

0.5

0.01

0.005

0 0.01

0.005

0

Figure 5: Sensors on the rst scale: 0 0 ( rst), 0 0 (second), 0 0 (third), 0 0 (fourth) k

k ;

;

k

;

k ;

We can notice that the local truncation error is maximum in the incoming channel region and around the step. The gradient in the x-direction is low, except around the upwind boundary region and in a narrow region around the step. The gradient in the y-direction shows that boundary and mixing layers should be re ned. The last term is high in a large region around the upwind boundary and after the step. The colors follow 15


a logarithmic law. Actually, the sensibility can be set between 10?1 and 10?4 . We can also notice that the maximum of all these sensors corresponds to the gradient in the y-direction, which is a not surprising. Nevertheless, we decide to choose the fourth one as the re nement sensor. The sensiblity " is 10?2 .

5.3 Re ned solution 0.01

0.005

0

-0.1

0

0.1

0.2

0.3

0.4

0.5

-0.1

0

0.1

0.2

0.3

0.4

0.5

0.01

0.005

0

Figure 6: Fine mesh (top) and global re ned mesh (bottom)

0.01

0.005

0

-0.1

0

0.1

0.2

0.3

0.4

0.5

-0.1

0

0.1

0.2

0.3

0.4

0.5

U 14 9 4 -1

0.01

0.005

0

Figure 7: Velocity eld u and streamlines: coarse mesh (top), global re ned mesh (bottom)

A new mesh is generated and added to the coarse mesh (Fig. 6). The velocity eld in the x-direction and the streamlines are plotted in Fig. 7. We can notice that the B-region starts to appear in the re ned mesh, whereas it was invisible in the coarse one. 16


As expected, the recirculation in the A-region is higher in the re ned mesh than in the coarse one. Thus, the choice of a sensor is not without consequences. Each sensor indicates dierent regions that need to be re ned. Therefore, the numerical results strongly depend on this choice in an adaptive mesh re nement procedure.

5.4 Velocity pro les 0.01

0.005

0 -5

0

5

10

15

20

Figure 8: Velocity pro le U(y) at x = 0:04 m. 35 10 (green), 70 20 (blue), global re ned mesh (red)

Velocity pro les are extracted from 35 10 and 70 20 coarse grids and from the global re ned mesh at x = 0:04 m. In this pro le, the recirculation in the A-region is characterized by a negative velocity in the lowest part of the ow (Armaly [3]). We notice that this negative velocity is easily seen only in case of the re ned mesh (Fig. 8). At y = 0 and y = 0:01 m, the velocity should be strictly equal to zero but is not due to interpolation errors in the post-treatment.

6 CONCLUSION

In this study, the aim was to show the links between adative mesh re nement methods and the wavelet approach. In the choice for a re nement sensor, for example, most authors work with heuristic sensors. In this article, this choice is motivated by the estimation of the computational improvement between two scales, which is equal to a corrective term on the low scale - the local truncation error - and the next term in the development of the function in the Haar wavelet basis. Computations performed for the backward-facing step ow at low Reynolds number 17


show that all the sensors are not equivalent. The choice of a sensor determines the regions that will be re ned.

ACKNOWLEDGEMENTS

The authors would like to thank Gilles Chaineray, whose assistance and contribution have been most valuable. This work is supported by ONERA general resources contract.

REFERENCES References

[1] Aftosmis M. J., \A second-order TVD method for the solution of the 3D Euler and Navier-Stokes equations on adaptively re ned meshes", 13th ICNMFD, lecture notes in physics, Vol. 414, 1991, pp. 235-239. [2] Angot P., Caltagirone J. P., Khadra K., \Une methode adaptative de ranement local: la correction du ux a l'interface", C. R. Acad. Sci. Paris, t. 315, serie I, 1992, pp. 739-745. [3] Armaly B. F., Durst F., Pereira J. C. F., Scho nung B., \Experimental and theoretical investigation of backward-facing step ow", J. Fluid Mech., Vol. 127, 1983, pp. 473-496. [4] Batina J. T., Rausch R. D., Yang H. T. Y., \Spatial adaptation of unstructured meshes for unsteady aerodynamic ow computations", AIAA Journal, Vol. 30, No. 5, 1992, pp. 1243-1251. [5] Beam R. M., Warming R. F., \An implicit nite-dierence algorithm for hyperbolic systems in conservation-law form", J. Comp. Phys., Vol. 22, 1976, pp. 87-110. [6] Benoit C., \Methode d'adaptation de maillages au moyen d'algorithmes genetiques pour le calcul d'ecoulements compressibles", Ph. D. Thesis, Ecole Nationale Superieure des Arts et Metiers, Paris, 1999. [7] Berger M. J., \Adaptive mesh re nement for partial dierential equations", Ph. D. Thesis, Standford, 1982. [8] Berger M. J., Colella P., \Local adaptive mesh re nement for shock hydrodynamics", J. Comput. Phys., Vol. 82, 1989, pp. 67-84. [9] Bihari B. L., Harten A., \Multiresolution schemes for the numerical solution of 2-D conservation laws I", SIAM J. Sci. Comput., Vol. 18, No. 2, 1997, pp. 315-354. [10] Braaten M. E., Connell S. D., \A 3D unstructured adaptive multigrid scheme for the Navier-Stokes equations", AIAA Paper 95-1728, 1995. [11] Brandt A., \Multi-level adaptive solutions to boundary value problems", Math. Comp., Vol. 31, 1977, pp. 333-390.

18


[12] Chiavassa G., Donat R., \Numerical experiments with point value multiresolution for 2D compressible ows", GrAN Report 99-02, Departament de Matematica Aplicada, Universitat de Valencia, 1999. [13] Coirier W. J., Powell K. G., \An accuracy assessment of Cartesian-mesh approaches for the Euler equations", AIAA Paper 95-3335, 1995. [14] Daubechies I., \Ten lectures on wavelets", CBMS lecture notes series, SIAM, 1991. [15] Dutoya D., Errera M. P., \Le code Mathilda. Modeles physiques, reseau de calcul et methode numerique", ONERA RT 42/3473, 1991. [16] Dutoya D., Errera M. P., \Une decomposition formelle du jacobien des equations d'Euler. Application a des schemas numeriques decentres", Rech. Aerosp., Vol. 1, 1992, pp. 25-35. [17] Farge M., "Wavelet transforms and their applications to turbulence", Annu. Rev. Fluid. Mech., Vol. 24, 1992, pp. 395-457. [18] Haar A., \Zur Theorie der orthogonalen Funktionensysteme", Math. Ann., Vol. 69, 1910, pp. 331-371. [19] Harran-Klotz P., \Maillages auto-adaptatifs et approximation des systemes hyperboliques separables. Application aux equations d'Euler tridimensionnelles", Ph. D. Thesis, Ecole Nationale Superieure de l'Aeronautique et de l'Espace, Toulouse, 1991. [20] Harten A., \Multiresolution algorithm for the numerical solution of hyperbolic conservation laws", Comm. Pure Appl. Math., Vol. 48, 1995, pp. 1305-1342. [21] Hentschel E., Hirschel H., \Self adaptive ow computations on structured grids", 1st ECCOMAS CFD Conference, 1994. [22] Jouhaud J. C., Borrel M., \A hierarchical adaptive mesh re nement method: application to 2D ows", 3rd ECCOMAS CFD Conference, Paris, 1996. [23] Jouhaud J. C., \Methode d'adaptation de maillages structures par enrichissement. Application a la resolution numerique des equations d'Euler et de Navier-Stokes", Ph. D. Thesis, Universite de Bordeaux I, 1997. [24] Lo hner L., \An adaptive nite element scheme for transient problems in CFD", Comput. Meth. Appl. Eng., Vol. 61, 1987, pp. 323-328. [25] Marcum D. L., Weatherhill N. P., Marchant M. J., Beaven F., \Adaptive unstructured grid generation for viscous ow applications", AIAA Paper 95-1726, 1995. [26] Mavrilis D. J., \Accurate multigrid solution of the Euler equations on unstructured and adaptive meshes", NASA Contractor Report 181679, ICASE Report No. 88-40, 1988.

19


[27] McCormick S. F., \Multilevel adaptive methods for elliptic eigenproblems: a two-level convergence theory", SIAM J. Num. Anal., Vol. 31, No. 6, 1994, pp. 1731-1745. [28] Powell K. G., De Zeeuw D., \An adaptively-re ned cartesian mesh solver for the Euler Equations", AIAA Paper 91-1542, 1991. [29] Powell K. G., Ashford G. A., \Adaptive unstructured triangular mesh generation and

ow solvers for Navier-Stokes equations at high Reynolds numbers", AIAA Paper 95-1724, 1995. [30] Quirk J. J., \An adaptive grid algorithm for computational shock hydrodynamics", Ph. D. Thesis, Cran eld Institute of Technology, College of Aeronautics, 1991. [31] Roussel O., \Ranement local de maillage par une methode de formulation faible en ondelettes", ONERA RT 1/6181 DSNA/N, 1999. [32] Schneider K., Kevlahan N.K.-R., Farge M., \Comparison of an adaptive wavelet method and nonlinearly ltered pseudospectral methods for two-dimensional turbulence", Theoret. Comput. Fluid Dynamics, Vol. 9, 1997, pp. 191-206. [33] Scho nfeld T., Rudgyard M., \A cell-vertex approach to local mesh re nement for the 3D Euler equations", AIAA Paper 94-0318, 1994.

20