Seyed Amin Ghaffari
email: ghaffari@l3m.univ-mrs.fr

Kai Schneider
email: kschneid@cmi.univ-mrs.fr

Keywords: Incompressible flows, Multiresolution analysis, Vorticity stream-function formulation, Volume penalization method, Central finite differences

Oh who is sitting at the eternal extreme of the existence, help me to rescue from the cemetery of

Développement d'une méthode multirésolution adaptative pour étudier le comportement des écoulements tourbillonnaires près des parois solides

Résumé L'objectif de cette recherche est de développer1 une méthode multi-échelle adaptative en temps et en espace plus efficace que les méthodes actuelles pour résoudre des équations aux dérivées partielles de type hyperbolique intervenant en mécanique des fluides. La nouvelle méthode, basée sur une discrétisation en différences finies d'ordre deux et une analyse multi-échelle, permet de réduire significativement le nombre de points nécessaires. La grille se raffine automatiquement dans les régions avec un fort gradient. La méthode est appliquée à l'équation de Burgers puis prolongé aux équations de Navier-Stokes bidimensionnelles. Pour étudier le comportement des écoulements tourbillonnaires près des parois solides, le problème du collision de dipôle sur paroi droit a été utilisé comme référence et ensuite le collision sur des parois courbes avec la méthode de pénalisation a été considéré. La stratégie d'adaptation est basée sur la transformée en ondelettes et le seuillage des coefficients. Pour intégrer les équations dans le temps, on utilise des méthodes de Runge-Kutta d'ordre différents, avec un pas de temps fixe ou adaptatif. Les résultats obtenus montrent que le temps de calcul peux être réduit considérablement avec cette méthode en conservant la précision.

Mots clés: Ecoulements incompressibles, Analyses multirésolutions , Formulation fonction de courant et vorticité, Méthode de pénalisation en volume, Différence fini centrée Chapter 1 Introduction

Objectives and historical background

The aim of present investigation is to develop a reliable self-adaptive direct numerical method to study the near wall behavior of unsteady incompressible flows. The governing equations are the Navier-Stokes equations with proper initial and boundary conditions. Conventional methods for spatial discretization of the PDEs (e.g., finite differences, finite volumes and finite elements) have limited order of accuracy especially near boundaries, but they are more flexible in dealing with complex geometries over a suitable grid. On the other hand standard spectral methods which are widespreadly used in direct numerical simulation of turbulence are limited to Cartesian grids. One can recognize the poor spectral localization (good spatial localization/resolution) of the former methods while good spectral localization (poor spatial localization/resolution) of the latter methods [START_REF] Powers | Manifold methods for energetic materials[END_REF]. Therefor for flows in/around complex geometries, the use of a grid capable to resolve all the scales will be limited to low-Reynolds flows. In general the small scales in a turbulent flow are not limited to near wall regions and can also move in the flow field. A high-order method with a prescribed grid spacing cannot resolve all the scales and thus the interactions present in the flow unless with the use of a fine grid in regions with high gradient. The limitation of mentioned methods for problems with widely disparate spatial scales, has encouraged the researcher to use alternative methods with limited accuracy but good spatial localization in regions with high gradient of flow variables. As a result recently there has been increasing interest in self adaptive numerical methods for solving time-dependent PDEs, see for example [START_REF] Siddiqui | Wavelet optimized finite difference method with non-static re-gridding[END_REF]. Adaptive methods can be divide into r-type (a fixed number of grid points are redistributed), h-type (regriding is performed occasionally) and p-type (the degree of the polynomial representing the solution is locally increased) each with their own advantages and disadvantages as detailed in literatures. Among different methods for grid adaptation h-type refinement proved to be more advantageous. By considering some advantages of the h-type refinement we explain some details of this methods. Existing numerical methods based on h-type refinement fall into two classes: error indicator based, where the grid is refined to resolve gradients of a physically relevant quantity, and error control based, where the error is estimated and the grid is refined to ensure this error is less than a prescribed tolerance. The error-indicating strategy does not control the error directly, but instead controls the mesh coarsening and refinement. The error-estimating strategy minimizes the error as measured in an appropriate norm, which leads to an optimal distribution of the grid and thus a more reliable solution.

Among different error-estimating adaptation strategies (which most of them belong to the finite element family) wavelet-based numerical methods have proved to be an efficient tool in developing adaptive numerical methods which control the global (usually L 2 ) approximation error. Wavelet transforms allow to estimate the local regularity of solutions of the given PDE, with a very efficient algorithm, and thus can define auto-adaptive discretization with local mesh refinement. Liandrat and Tchamitchian [START_REF] Liandrat | Resolution of the 1D regularized burgers equation using a spatial wavelet approximation[END_REF] proposed the first wavelet-based adaptive method for numerical simulation of PDEs. The currently existing wavelet-based algorithms can be classified as pure wavelet methods and wavelet optimized grid methods. Pure wavelet methods, divided to Galerkin and collocation schemes, employ wavelets directly for discretization of the governing equations. The wavelet transform facilitates the effective sparse representation of the functions and pseudo-differential operators (and their inverse) by filtering of wavelet coefficients of the function and of the matrix representing the operators. The characterization of function spaces in terms of wavelet coefficients and the corresponding norm equivalences leads to diagonal preconditioning of operators (compression) in wavelet space. In practice similar to pseudo-spectral methods the evaluation of nonlinear term in function space is not efficient, usually by performing an inverse wavelet transform in each time step, the solution is reconstructed on a locally refined grid, and the nonlinear term is then evaluated point-wise in physical space. For a comprehensive review of wavelet methods applied to the computational fluid dynamics see [START_REF] Schneider | Wavelet methods in computational fluid dynamics[END_REF]. On the other hand, wavelet optimized grid (WOG) methods combine classical discretizations of considered equations (e.g., finite differences or finite volumes) with wavelets, which are used to define the adaptive grid. The goal of the wavelet-based adaptation of the grid is to obtain the best approximation of the flow variables on a near optimal grid. This approximation has a one-to-one correspondence between the wavelet expansion coefficients and the grid points. Thus, nonlinear filtering of wavelet coefficients of the interested flow variables automatically coarsens the computational grid. Since with the use of wavelet transform and thresholding of the wavelet coefficients, a given functions can be reconstructed with a prescribed accuracy, thus adaptive methods based on wavelets provide global error control for numerical solution of PDEs. See [START_REF] Roussel | Développement d'un algorithme multirésolution adaptatif tridimensionnel pour les équations aux dérivées partielles paraboliques. Application l'étude des instabilités thermo-diffusives de flamme[END_REF] and [START_REF] Hejazialhosseini | High order finite volume methods on wavelet-adapted grids with local time-stepping on multicore architectures for the simulation of shock-bubble interactions[END_REF] where a finite volume discretization of governing equations combined with cell-averaged interpolating wavelet transform for grid adaptation. In the present work the Navier-Stokes equations are discretized with a classical second order central finite differences, thus similar to WOG methods the role of the wavelet transform is the adaptation of the grid and the fast interpolation of flow variables at new unknown points. Wavelet optimized finite difference (WOFD) methods based on fast wavelet transforms are very efficient (especially for problems with a large magnitude of scales, e.g., turbulence) since the computational cost scales linearly with the number of wavelet coefficients retained in the approximation [START_REF] Alam | Simultaneous space-time adaptive wavelet solution of nonlinear parabolic differential equations[END_REF].

In the procedure of solving the incompressible Navier-Stokes equations a difficulty arises because of the lack of a time-derivative term in the continuity equation, which limits the straightforward use of time-marching numerical methods designed for coupled equations. Another problem is the procedure of updating the pressure, which requires special treatment. Most numerical methods like projection, marker and cell and fractional step methods involve solving a decoupled Poisson equation for the pressure at each time step. The common application of this method consist of two steps: The first step is to solve for an intermediate velocity field using the momentum equations (Burgers-like), in which a fraction of the pressure-gradient term from the previous time step can be involved or it can be excluded entirely. In the second step, the pseudo pressure is computed via intermediate velocity field, which maps (Leray projector) the intermediate velocities onto a divergence-free (selonoidal) velocity field. This method is first order in time, originally proposed by Chorin [START_REF] Chorin | Numerical solution of Navier-Stokes equations[END_REF] in 1968 and developed by Kim and Moin [START_REF] Kim | Application of a fractional-step method to incompressible Navier-Stokes equations[END_REF] to second order accuracy. For a comprehensive review of projection methods see [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]. Another approach is the artificial compressibility method (ACM), which relaxes the ellipticity of the pressure field by adding a pseudotime derivative of pressure into the continuity equation, hence coupling continuity and momentum equations. Then iterations are performed in pseudo-time until attaining a divergence free velocity field. The method was originally designed for steady problems by Chorin [START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF] in 1967 and developed by Rogers and Kwak [START_REF] Rogers | Steady and unsteady solutions of the incompressible Navier-Stokes equations[END_REF] to second order accuracy in time, with the use of dual time stepping strategy. For a comprehensive review of ACM methods see [START_REF] Kwak | Computational challenges of viscous incompressible flows[END_REF]. However in two-dimensional problems vorticity and stream-function formulation would be the best for some purposes by eliminating the pressure gradient entirely and reducing the flow variables to one scalar value.

In the present work the concept of adaptive multiresolution method will be applied to the vorticity and stream-function formulation. However the developed concepts are also applicable to the primitive variable formulation. See [START_REF] Durbin | An approach to local refinement of structure grids[END_REF] and [START_REF] Ding | A stencil adaptive algorithm for finite difference solution of incompressible viscous flows[END_REF] for a finite difference discretization of governing equations using gradient approach for grid adaptation. To the best of our knowledge this investigation is the first attempt to develop an adaptive finite difference method based on the interpolating wavelet transform applied to the vorticity and stream-function formulation. A second-order central finite difference method with symmetric stencil over an adaptive Cartesian grid is used for spatial discretization of the equations. Finite difference method represents a suitable combination with the multiresolution analysis based on the Harten's pointvalue wavelet transform. The chosen discretization avoids adding excessive numerical dissipation which is usually introduced by upwind methods, thus has a reasonable accuracy in space over an adapted grid. The concept of symmetric stencils will lead to intermediate (hung) points, that their values can be interpolated accurately via inverse wavelet transform. Another advantage of the vorticity and stream-function formulation is that it allows the use of high-order time stepping rather easily, e.g., fourth-order Runge-Kutta method, to minimize the truncation error due to discrete time integration which is in nature an accumulative error during simulations [START_REF] Lomax | Fundamentals of Computational Fluid Dynamics[END_REF].

In the present investigation the unexpectedly difficult, dipole-wall collision is chosen as a benchmark computation. The problem is studied experimentally with the use of PIV technique by Naguib and Koochesfahani [START_REF] Naguib | On wall-pressure sources associated with the unsteady separation in a vortex-ring/wall interaction[END_REF]. Numerical simulations of the dipole-wall collision were first performed by Orlandi [START_REF] Orlandi | Vortex dipole rebound from a wall[END_REF] in two dimensions, who revealed that the creation of vorticity at the wall drives a vigorous rebound of the dipolar vortex. Each dipole half forms a new pair with a secondary vortex consisting of boundary layer vorticity. The new dipoles travel along circular trajectories away from the wall. The production of the secondary vorticity at the wall during the dipole-wall collision was later studied by Coutsias and Lynov [START_REF] Coutsias | Fundamental interactions of vortical structures with boundary layers in two-dimensional flows[END_REF] and by Clercx and van Heijst [START_REF] Clercx | Dissipation of kinetic energy in twodimensional bounded flows[END_REF]. The Reynolds number in the reported simulations is limited to Re = 5000 as accurate representation of the small-scale structures near the wall requires very high resolution. Walker [START_REF] Walker | The boundary layer due to a rectilinear vortex[END_REF] and Peridier et al. [START_REF] Peridier | Vortex-induced boundary layer separation. Part 1. The unsteady limit problem Re → ∞[END_REF][START_REF] Peridier | Vortex-induced boundary layer separation. Part 2. Unsteady interacting boundary-layer theory[END_REF] have investigated the behavior of the boundary layer at a no-slip wall affected by a nearby single (point) vortex. This setup permits the use of boundary-layer equations, i.e., a reduced version of the Navier-Stokes equations, on a specially tailored grid to obtain results for large Reynolds numbers. More recently Obabko and Cassel [START_REF] Obabko | Navier-Stokes solutions of unsteady separation induced by a vortex[END_REF] reported results on this problem based on the full Navier-Stokes equations [START_REF] Kramer | Vorticity dynamics of a dipole colliding with a no-slip wall[END_REF].

After validation of the developed adaptive multiresolution method with the results of previous studies of a dipole collision with a straight wall the extension to collision with curved walls via the volume penalization method will be presented. Nowadays, there exist several techniques to introduce solid curved boundaries in the numerical simulation of different types of flows. Among those, some of them use finite difference methods by mapping the equations to general curvilinear coordinates [START_REF] Rogers | Steady and unsteady solutions of the incompressible Navier-Stokes equations[END_REF], others use full unstructured grids with finite volume or element discretization of the governing equations, we will call them body-fitted grids. All of them have their own advantages and difficulties. Another approach which became more popular during recent years are immersed boundary methods for imposing complex geometries while using the advantages of Cartesian grids. They are easy-to-implement and more efficient than classical approaches such as body-fitted grids, in particular for moving and/or flexible geometries [START_REF] Kolomenskiy | Simulation numérique direct des écoulement tourbillonaires autour des obstacles mobiles[END_REF]. The name immersed boundary comes from the fact that the physical solid boundary does not always conform with the computational grid points at boundaries, instead it is inside of the solution domain. The volume penalization method originally introduced by Arquis and Caltagirone [START_REF] Arquis | Sur les conditions hydrodynamiques au voisinage d'une interface milieu fluid -milieux poreux: application à la convection naturelle[END_REF] for flows in porous media, does posses this property; therefore it is an example of immersed boundary methods. This technique is physically motivated, since it is based on the idea of modeling solids as porous media whose porosity coefficient tends to zero. In addition, it is mathematically justified, due to the convergence property proven rigorously by Angot et al. [START_REF] Angot | A penalization method to take into account obstacles in viscous flows[END_REF].

As a starting point in the present work we take the two-dimensional vorticity streamfunction solver developed in [START_REF] Sabetghadam | Implementation of vortex stretching into the two-dimensional Navier-Stokes equations via arbitrary external straining[END_REF] for a uniform grid and the adaptive one-dimensional Burgers solver developed by author at M2P2 (Marseilles) during the last summer. The manuscript is organized as follows: In the following a short introduction to the discrete wavelet transforms and the idea of point selection by filtering of the wavelet coefficients will be presented. Then in chapter two some theoretical background of the adaptive multiresolution method will be introduced. After that the results for the one-dimensional Burgers equation and the dipole collision with either straight or curved walls will be reported. Finally, conclusions and perspectives will be discussed in chapter four.

Wavelet transform

Like Fourier transforms, wavelet transforms can be viewed as a change of basis in function space, e.g., from time to frequency or from space to wavenumber. The basis functions of the Fourier transform are trigonometric functions, i.e., sines and cosines. For the wavelet transform there is a large choice of basis functions which are called mother wavelets. For the Fourier transform the decay of the Fourier coefficients depends on the global regularity of the decomposed function. In wavelet transforms some of the new basis functions are more complicated than trigonometric functions but what makes them interesting is that, unlike the Fourier basis, wavelets are well localized functions in space; simultaneously, like them, they are quit localized in function space or more precisely they have a characteristic scale. So they can unfold signals or fields into both space (or time) and scale, and possibly directions in dimensions higher than one. The continuous wavelet transform has been discovered by Grossmann and Morlet who published the first paper on wavelets in 1984 [START_REF] Goupillaud | Cycle-octave and related transforms in seismic signal analysis[END_REF]. The orthogonal wavelet transform has been discovered by Lemarié and Meyer (1986). Then, Daubechies (1988) found orthogonal bases made of compactly supported wavelets, and Mallat (1989) designed the fast wavelet transform (FWT) algorithm. More details can be found in [START_REF] Schneider | Wavelets: Mathematical Theory[END_REF]. The second generation wavelet transform (SGWT) was introduced independently by Sweldens [START_REF] Sweldens | The lifting scheme: A construction of second generation wavelets[END_REF], Donoho [START_REF] Donoho | Interpolating wavelet transforms[END_REF] and Harten [START_REF] Harten | Multiresolution Representation of data: A general framework[END_REF] in the early 1990s. It is a wavelet transform where the filters (or even the represented wavelets) are not designed explicitly, but rather the transform consists of calculating the differences between an interpolated value (with different accuracies) of the function at a given point from the corresponding neighbor points in lower levels, and its value.

The main advantage of this approach is that it is applicable to either periodic or non-periodic functions.

Biorthogonal wavelet transform

In a discrete multiresolution analysis framework, data are represented at different scale levels, and the main tools are appropriate transformations, e.g., wavelet transform, relating the information f J at the finest scale level J to the lower ones, and vice versa. After a transformation the output contains the information f 0 on the coarsest level, and d j that keeps the details between a scale level j and the next upper level j + 1. Schematically, we have

f J ↔ f M R J = (f 0 , d 0 , ..., d J-1 ) (1.1)
To illustrate this methodology, we consider the case of discretization by Harten's point values [START_REF] Bihari | Multiresolution schemes for the numerical solution of 2-D conservation laws[END_REF] for uniform grids, which is well adapted for finite difference methods, versus Harten's cell average method which is more suitable for finite volume methods. By considering in a unit interval the hierarchy of uniform dyadic grids will obtain from

X j = {x j,i ∈ R : x j,i = i2 -j , i = 0, • • • , 2 j }, j = 0, • • • , J (1.2)
with spacing 2 -j , where j is the level and i represents the position. The number of points must always be odd (N = 2 J + 1) to have a point in the middle. As indicated in Fig. 1.1, in multiresolution analysis to go from X j to a more refined grid X j+1 , Figure 1.1: Uniform dyadic grids with N = 2 J + 1 = 17 points in finest level (J = 4), including lower levels down to (j = 0).

we add to X j new midpoints between the old points, dividing by two the step size.

The direct and inverse multiresolution (MR) transforms are usually associated to bi-orthogonal multiresolution representation of functional spaces of the form

V J = V 0 ⊕ J ∑ j=0 W j (1.3)
where V J and W J are associated function spaces. Thus a given function, f (x) can be represented with

f (x) = 2 J ∑ i=0 f 0,i Φ 0,i (x) + J ∑ j=0 2 J ∑ i=0 d j,i Ψ j,i (x) (1.4)
where the orthonormal basis are scaling functions Φ j,i and wavelets Ψ j,i defined as

Ψ j,i (x) = 2 j/2 Ψ(2 j x -i), j, i ∈ Z Φ j,i (x) = 2 j/2 Φ(2 j x -i), j, i ∈ Z
At a given scale (j) the scaling function Φ j,i is orthonormal with respect to its translation by discrete steps i2 -j but not with respect to its dilates,

⟨Φ j,i , Φ j,k ⟩ = δ ik (1.5)
The wavelets Ψ j,i are orthonormal with respect to their translates by discrete steps i2 -j and their dilates by discrete steps 2 -j corresponding to octaves [START_REF] Schneider | Wavelets: Mathematical Theory[END_REF].

⟨Ψ j,i , Ψ j ′ ,i ′ ⟩ = δ jj ′ δ ii ′ (1.6)
where δ denotes the Kronecker symbol and the inner (scalar) product is defined by

⟨f, g⟩ = ∫ R f (x)g * (x)dx (1.7)
where * denotes the complex conjugate in the case of complexed valued wavelets. The scaling function coefficients, f j,i , essentially encode the smooth part of the function, while the wavelet coefficients, d j,i , contain information on the function behavior on successively finer scales and the coefficients are interpreted as the inner products with the dual basis

f j,i = ⟨f, Φ j,i ⟩ = f (x j,i ) (1.8)
where Φ j,i is the delta distribution at x j,i . Definition of the dual wavelets Ψ j,i such that

d j,i = ⟨f, Ψ j,i ⟩ = f j+1,2i+1 -fj+1,2i+1 (1.9) implies Ψ j,i (x) = Φ j+1,2i+1 - ∑ n β n Φ j,n (x) (1.10) thus f , can summarized in fj+1,2i+1 = st/2 ∑ n=1 β n (f j,i-n+1 + f j,i+n ) (1.11)
where st is the interpolation stencil. For instance, interpolating coefficients for linear interpolation with st = 2 is β 1 = 1/2, and for cubic interpolation with st = 4 are (β 1 , β 2 ) = (9/16, -1/16). The following equation is useful to better understand the relation of indexes between two successive level for a given variable ϕ ϕ j+1,2i+1 ⇐⇒ ϕ j,i+ 1 2

(1.12)

One-dimensional transform

The direct wavelet transform (WT) for a discrete function given on dyadic grid points, is performed simply by replacing the value of the function, with the difference between its value and its interpolated value from the corresponding neighbors at lower levels

d j,i = f j+1,2i+1 -fj+1,2i+1 (1.13)
this must be performed from the finest level J down to j = 1. In contrast the inverse wavelet transform (IWT) starts from the coarsest (lowest) level up to the finest level

f j+1,2i+1 = d j,i + fj+1,2i+1 (1.14)
In this method predictions of the values of the function will be done with interpolating polynomials. The prediction is exact for polynomials of a prescribed degree [START_REF] Domingues | Adaptive multiresolution methods[END_REF]. The most simple example is given by linear (first-order) interpolation

fj+1,2i+1 = f j,i + f j,i+1 2 (1.15)
By using cubic (third-order) interpolation transform of a unit vector (Dirac-delta function δ) in the mid component of a zero vector (of length, e.g., 8192) so the point is belonging to the minimum level. Other unit vectors would give wavelets with the same shapes, but different positions and scales. The smooth scaling functions, e.g., the one obtained with cubic interpolation (see Fig. 1.3) decays rapidly in Fourier space. Close to the left boundary forward interpolation must be used for cubic interpolation, resulting in fj+1,1 = 5f j,0 + 15f j,1 -5f j,2 + f j,3 16 (1.17)

fj+1,2i+1 = -f j,i-1 + 9f j,i + 9f j,i+1 -f j,i+2 16
and on the right boundary a backward stencil must be used as follows fj+1,2 j+1 -1 = f j,2 j -3 -5f j,2 j -2 + 15f j,2 j -1 + 5f j,2 j 16 (1.18) With the use of cubic interpolation, i.e., four point interpolation stencil (st = 4) it is not possible to go lower than j = 3 during direct and inverse wavelet transform. The pseudo-code for direct wavelet transform with linear interpolation of a given function f (x i ) over i = 0, . . . , 2 J points, excluding boundary points is given in the following: 1. DO for j = J to j = 1, step = -1

2.

P ointer := 2 J-j

3.

DO for i = P ointer to i = 2 J -P ointer, step = 2 × P

4.

Interpol := 0.5f (i -P ointer) + 0.5f (i + P ointer)

5.

d(i) := f (i) -Interpol

END DO

END DO

We can note that no intermediate memory is needed in this transformation and the order of operations is O(cN ), where N is the number of points and the constant c depends on the order of interpolation which is usually less than 10. The pseudo-code for inverse transformation is given in the following:

1. DO for j = 1 to j = J, step = 1

2.

P ointer := 2 J-j

3.

DO for i = P ointer to i = 2 J -P ointer, step = 2 × P

4.

Interpol := 0.5f (i -P ointer) + 0.5f (i + P ointer)

5. f (i) := d(i) + Interpol 6. END DO 7. END DO 1.2.

Filtering in wavelet space

Given a threshold parameter for the finest ϵ J or the coarsest ϵ 1 level, data compression will be obtained by thresholding of the detail coefficients, also called nonlinear filtering, applied to wavelet coefficients in wavelet space. After performing the direct transform, wavelet coefficients smaller than a threshold are set to zero and the corresponding point can be eliminated from the set of the points, in other words we can find the value of that point by interpolation and the error remains bounded by the threshold value.

d j,i = { 0 if |d j,i | ≤ ϵ j , d j,i else (1.19)
where [START_REF] Harten | Multiresolution Representation of data: A general framework[END_REF] in which D = 1, 2, 3 is the dimension of the problem, and J denotes the maximum level. A linear filtering also can be used, i.e., all wavelet coefficients above a given scale are set to zero. In the present investigation nonlinear filtering is used. After nonlinear filtering in wavelet space the given function f (x), can be reconstructed f (x), just with the significant wavelet coefficients corresponding to the important points of the function. Those points must be kept to guaranty the boundedness of the error introduced due to filtering and eliminating non necessary points. Following Donoho [START_REF] Donoho | Interpolating wavelet transforms[END_REF], it can be shown that for a sufficiently smooth function f (x)

ϵ j = ϵ J 2 D(j-J) = ϵ 0 2 D(j) (1.
|f (x) -f (x)| ≤ c 1 ϵ 0 (1.21)
This implies that the number of significant wavelet coefficients N s is bounded by ϵ 0 as

N s ≤ c 2 ϵ -D/P W T 0 (1.22)
where P W T is the order of the wavelet, i.e, the number of neighboring points, or the stencil size, used for wavelet construction [START_REF] Sweldens | The lifting scheme: A construction of second generation wavelets[END_REF] and the coefficients c i depend on f (x) (but are of order unity). Note that P W T controls the number of zero moments of the interpolating scaling function. Combining Eq. (1.21) and Eq. (1.22) we have the following bounded error in terms of

N s |f (x) -f (x)| ≤ c 3 N -P W T /D s (1.23)
The estimated error is consistent with numerical experiments for both one-dimensional and two-dimensional cases [START_REF] Vasilyev | An adaptive multilevel wavelet collocation method for elliptic problems[END_REF]. The accuracy of differentiation with wavelet transform applied to a function over a uniform and thresholded grid was examined by Vasilyev and Bowman [START_REF] Vasilyev | Second generation wavelet collocation method for the solution of partial differential equations[END_REF] for the one dimensional case and by Vasilyev [START_REF] Vasilyev | Solving multi-dimensional evolution problems with localized structures using second generation wavelets[END_REF] in multiple dimensions. It was shown that the error bound on the derivative is given by

|D U ni x f (x) -D M R x f (x)| ≤ c 4 N -(P W T -1)/D s ≈ c 4 ϵ (P W T -1)/P W T 0 (1.24)
where D x stands for the derivative operator in the x direction. This relation was verified numerically for both one and two dimensions in [START_REF] Vasilyev | Second generation wavelet collocation method for the solution of partial differential equations[END_REF] and [START_REF] Vasilyev | Solving multi-dimensional evolution problems with localized structures using second generation wavelets[END_REF]. Note that the error bound (1.24) is also correct for the second-order derivative if a symmetric stencil is used [START_REF] Vasilyev | An adaptive multilevel wavelet collocation method for elliptic problems[END_REF]. The discrete L p -norm for error estimation is defined by

L p (f, f ) = [ 1 N N ∑ i=1 |f i -f i | p ] 1 p (1.25)
and L ∞ -norm is defined as

L ∞ (f, f ) = max |f i -f i |, i = 1, . . . , N. (1.26)
The compression rate can be introduced as

% Compression = ( N -N s N ) × 100 (1.27)
where N is the number of total points and N s is the number of significant points after thresholding in wavelet space. The idea of filtering and compression with error analysis will be shown for two examples.

Example: Consider a non-periodic one-dimensional function f (x) over [0, 1]

f (x) =    8.1e 1/4 e -|x-1/2| 0.0 ≤ x < 0.25 9e -|x-1/2| 0.25 ≤ x < 0.75 e -|x-1/2| (16x 2 -24x + 18) 0.75 ≤ x ≤ 1.0 (1.28)
with a discontinuity at x = 0.25, a discontinuity in the first derivative at x = 0.5 and a discontinuity in second derivative at x = 0.75. As a second example we consider a Gaussian function f (x) = exp ( x -0.5 δ

) 2 for x ∈ [0, 1]
their sparse point representations, with the use of cubic (P W T = 4) interpolating wavelet transform, for J = 10, filtered with threshold ϵ = 1 × 10 -3 are illustrated in Fig. 1.4 (a) and (b). Results of the error analysis and compression rates with the linear and cubic interpolating wavelet transforms, by applying different thresholding parameters are reported in Tables 1.1 and 1.2.

WT ϵ % Compression L ∞ -Error L 1 -Error L 2 -Error Linear 1 × 10 -1 97.2 3 × 10 -3 9 × 10 -4 4 × 10 -5 Linear 1 × 10 -2 94.4 6 × 10 -4 2 × 10 -4 8 × 10 -6 Linear 1 × 10 -3 88.4 2 × 10 -4 5 × 10 -5 2 × 10 -6 Linear 1 × 10 -5 52.7 1 × 10 -5 2 × 10 -6 1 × 10 -7 Linear 1 × 10 -6 4.1 9 × 10 -7 2 × 10 -8 4 × 10 -9 Linear 1 × 10 -7 0.2 1 × 10 -9 3 × 10 -12 2 × 10 -12 Linear 1 × 10 -8 0.2 1 × 10 -9 3 × 10 -12 2 × 10 -12 Linear 1 × 10 -10 0.0 1 × 10 -16 6 × 10 -18 7 × 10 -19 Cubic 1 × 10 -1 96.3 7 × 10 -2 6 × 10 -2 1 × 10 -4 Cubic 1 × 10 -2 94.4 6 × 10 -4 2 × 10 -4 8 × 10 -6 Cubic 1 × 10 -3 96.8 3 × 10 -5 7 × 10 -6 3 × 10 -7 Cubic 1 × 10 -4 94.2 6 × 10 -6 8 × 10 -7 5 × 10 -8 Cubic 1 × 10 -6 87.6 1 × 10 -8 3 × 10 -9 1 × 10 -10 Cubic 1 × 10 -8 71.6 5 × 10 -9 1 × 10 -9 5 × 10 -11 Cubic 1 × 10 -10
12.1 1 × 10 -10 7 × 10 -12 7 × 10 -13 Cubic 1 × 10 -12 0.0 7 × 10 -16 9 × 10 -17 5 × 10 -18

Two-dimensional transform

The previously presented one-dimensional wavelet transform can be extended to higher dimensions. For the moment just the two-dimensional transform will be considered, since the three-dimensional transform can be constructed analogously. There are two ways to accomplish this [START_REF] Schneider | Wavelets: Mathematical Theory[END_REF]:

1. Rectangular wavelet transform

Tensor product wavelet transform

Rectangular wavelet transform : Like the Fourier transform for higher dimensions this transform is simply equivalent to applying the one dimensional wavelet transform to the rows and the columns of a matrix or a two-dimensional function. For some applications, this method is more advantageous. Often the notion of a scale has a certain meaning, so one would like to have a unique scale assigned to Tensor product wavelet transform : This transform can be obtained through the tensor product of two one-dimensional multiresolution analysis of L 2 (R), so three different wavelets will be generated

f (x, y) = ∑ jx,
f (x, y) = ∑ j ∑ ix,iy 3 ∑ µ=1 d µ j,ix,iy Ψ µ j,ix,iy (x, y) (1.32)
where

Ψ µ j,ix,iy (x, y) =    Ψ j,ix (x)Φ j,iy (y), µ = 1 Φ j,ix (x)Ψ j,iy (y), µ = 2 Ψ j,ix (x)Ψ j,iy (y), µ = 3 (1.33) and d µ j,ix,iy = ⟨f, Ψ µ j,ix,iy ⟩ = f j/j+1 (i x , i y ) -f µ j/j+1 (i x , i y ) (1.34)
The scale parameter j simultaneously controls the dilatation in x and y directions, and in D dimensions this construction yields to 2 D -1 types of wavelets [START_REF] Schneider | Wavelets: Mathematical Theory[END_REF].

Application to interpolating wavelet transform : For a two-dimensional function f (x, y) with the use of the notion of tensor product, the interpolating wavelet transform will take three different forms for locations with following (i x , i y ) index pair, where i x = 1, . . . , 2 J+1 and i y = 1, . . . , 2 J+1 are the indexes in x and y directions, respectively [START_REF] Bihari | Multiresolution schemes for the numerical solution of 2-D conservation laws[END_REF].

(odd,odd): Will not change (they belong to the lower levels) see Fig. 1.6.

(even,odd): Interpolation along x direction, µ = 1

f 1 j+1 (2i x , 2i y + 1) = st/2 ∑ n=1 β n ( f j (i x + n -1, i y ) + f j (i x -n, i y ) ) (1.35) (odd,even): Interpolation along y direction, µ = 2 f 2 j+1 (2i x + 1, 2i y ) = st/2 ∑ m=1 β m ( f j (i x , i y + m -1) + f j (i x , i y -m) ) (1.36)
(even,even): 2D interpolation at the same level (diagonal direction), µ = 3

f 3 j (2i x , 2i y ) = st/2 ∑ n=1 β n st/2 ∑ m=1 β m ( f j (i x + n -1, i y + m -1) + f j (i x -n, i y + m -1) +f j (i x + n -1, i y -m) + f j (i x -n, i y -m) )
In the above equations the coefficients β i , (i = 1, 2, . . . , st/2) come from standard interpolation methods, where (st -1) is the order of the interpolation, (st being the stencil size). For a few common values of st we have The pseudo-code for the direct wavelet transformation, with linear interpolation, for a matrix or a two-dimensional function f (x ix , y iy ), i x , i y = 1, . . . , N , with N = 2 J + 1 points in each direction including boundary points, is given in the following:

1. DO for j = J to j = 0, step = -1

2. P := 2 J-j 3. DO for i y = 1 + P to i y = N -P , step = 2 × P 4. DO for i x = 1 + P to i x = N -P , step = 2 × P 5. Interpol := 0.25× (f (i x -P, i y )+ f (i x +P, i y )+ f (i x , i y -P )+ f (i x , i y +P )) 6. d(i x , i y ) := f (i x , i y ) -Interpol END DO END DO 7. DO for i y = 1 to i y = N , step = 2 × P 8. DO for i x = 1 + P to i x = N -P , step = 2 × P 9. Interpol := 0.5 × (f (i x -P, i y ) + f (i x + P, i y )) 10. d(i x , i y ) := f (i x , i y ) -Interpol END DO END DO 11. DO for i y = 1 + P to i y = N -P , step = 2 × P 12. DO for i x = 1 to i x = N , step = 2 × P 13. Interpol := 0.5 × (f (i x , i y -P ) + f (i x , i y + P )) 14. d(i x , i y ) := f (i x , i y ) -Interpol END DO END DO

END DO

The inverse transform will be similar by starting from the lowest level j = 0 up to j = J and performing the two-dimensional interpolations after interpolations in x and y directions because in each level two-dimensional interpolations are using values from their own level. In Fig. 1.6 a schematic representation of two-dimensional wavelet transform based on cubic interpolation and tensor product is shown.

Example: Consider a two-dimensional function with a discontinuity In the onedimensional interpolations in x and y directions, (µ = 1, 2) points at the finer level j + 1 (marked •) must interpolated from the points at the coarser level j (marked •), illustrated in (a). In the two-dimensional interpolation, (µ = 3) point at the level j + 1 (marked •) must interpolated from the points at the same level j + 1 (marked •), illustrated in (b), picture taken from [START_REF] Vasilyev | An adaptive multilevel wavelet collocation method for elliptic problems[END_REF]. Table 1.3: Results of the error analysis and compression rates for the (discontinues) function (1.37) with the linear and cubic interpolating wavelet transforms, J = 9.

f (x, y) = { 1 if (x -1/2) 2 + (y -1/2) 2 ≤ 1 0 else (1.
WT ϵ % Compression L ∞ -Error L 1 -Error L 2 -Error Linear 1 × 10 -1
99.

WT ϵ % Compression L ∞ -Error L 1 -Error L 2 -Error Linear 1 × 10 -1 87.6 4 × 10 -2 9 × 10 -3 3 × 10 -5 Linear 1 × 10 -2 58.8 6 × 10 -3 2 × 10 -3 4 × 10 -6 Linear 1 × 10 -3 37.2 5 × 10 -4 1 × 10 -4 3 × 10 -7 Linear 1 × 10 -5 7.7 2 × 10 -6 2 × 10 -7 1 × 10 -9 Linear 1 × 10 -8 0.0 3 × 10 -16 4 × 10 -18 3 × 10 -20 Cubic 1 × 10 -2
99. Chapter 2

An adaptive multiresolution method for evolutionary PDEs

Introduction

In this chapter first the theoretical basis of an adaptive multiresolution method for evolutionary PDEs will be explained in some detail and two time integration methods from the Runge-Kutta family will be presented. Then the space discretization based on finite differences on Cartesian grids coupled with wavelet optimized grid (WOG) adaptation will be described. Finally, the method will be extended to the two-dimensional Navier-Stokes equations in vorticity and stream-function formulation including the penalization term which imposes no-slip and no-penetration boundary conditions, required in the case of complex geometries, e.g., curved walls. Penalization methods allow to impose boundary conditions in the case of complex geometries (even varying in time) without modifying the computational grid.

Evolutionary PDEs

In fluid mechanics we encounter different types of partial differential equations. They come from conservation laws and each of them describes the evolution of a different quantity, i.e., momentum, energy, vorticity and etc. From mathematical point of view these equations can be divided into three main groups, i.e., hyperbolic, parabolic and elliptic. The Euler equations are an example of a hyperbolic system of equations. Two examples for parabolic equations are : the one-dimensional viscous Burgers equation

∂ t u + uu x = νu xx (2.1)
and the two-dimensional vorticity (ω = ∇ × u) transport equation

∂ t ω + (u • ∇)ω = ν∇ 2 ω. (2.2)
As an example for elliptic PDEs, Poisson equation for stream function ψ can be cited. These equations must completed with suitable initial and boundary conditions. In unsteady incompressible flows the governing system of equations, i.e., Navier-Stokes equations, have parabolic character in time (for steady problems they have elliptic character in space). The continuity equation, must be satisfied in all time steps so the evolutionary part of the problem will be the momentum (or vorticity) equations. Hence for instance, we will focus on the time dependent equations in the form of :

-∇ 2 ψ = ω (2.3)
∂u ∂t = RHS (2.4)
where RHS contains all spatial terms either linear or nonlinear. It is also possible to write the equations in conservative form

∂u ∂t = ∇ • F (u) on Ω (2.5)
where F contains all viscous and inviscid fluxes.

Spatial discretization

For discretization of spatial derivatives a second order central finite difference method over Cartesian grids will be used as follows:

First derivative:

∂u ∂x i = u i+1 -u i-1 2h + O(h) 2 (2.6)
Second derivative:

∂ 2 u ∂x 2 i = u i+1 -2u i + u i-1 h 2 + O(h) 2 (2.7)
For diffusion terms with a variable coefficient λ we have:

∂ ∂x ( λ ∂u ∂x ) i = (λ i+1 + λ i )(u i+1 -u i ) -(λ i-1 + λ i )(u i -u i-1 ) 2h 2 + O(h) 2 (2.8)
For maintaining the accuracy in the multiresolution analysis by using symmetric stencils the indices (±1) are replaced by (±d), where d is an integer number, multiplying by h, representing the distance (dh) between the point in i and the nearest active point. For example Eq. (2.7) will be replaced by

∂ 2 u ∂x 2 i = u i+d -2u i + u i-d (dh) 2 + O(h) 2 (2.9)
For a given grid point after finding the nearest active point, all neighbors will be checked with the same distance to verify if they are in the list of active points, if not, we will call them hung points and their values must be interpolated from active points in lower levels, see Fig. 2.1. For two-dimensional problems the same relations will be applied to x and y directions.

Another alternative to the idea of the hung point interpolations is the use of nonsymmetric stencils. In this approach evaluation of the first derivative via a nonsymmetric stencil reads:

∂u ∂x i = a 2 u i+b + (b 2 -a 2 )u i -b 2 u i-a ab(a + b)h + O(h) (2.10)
Second derivative can evaluated by:

∂ 2 u ∂x 2 i = 2au i+b -2(a + b)u i + 2bu i-a ab(a + b)h 2 + O(h) (2.11)
where a and b are integer numbers, multiplied by h, they represent the distances from left and right neighbors respectively [START_REF] Bodeau | Non uniform grids for PDE in finance[END_REF].

Time integration

In this section the time integration for a general evolutionary PDE is described and the results will be assessed in the next chapter. For the moment just the explicit Rung-Kutta family is considered. A classical fourth order Runge-Kutta method with fixed time step and an adaptive time stepping with embedded Runge-Kutta methods will be presented.

Explicit fourth order Runge-Kutta method

Because of high accuracy and straightforward parallelization, the fourth order Runge-Kutta method is one of the best and mostly used methods for integration of ordinary differential equations [START_REF] Press | Numerical Recipes in Fortran 77[END_REF]. By putting all discretized spatial derivatives in the RHS operator one can solve the considered partial differential equation as a system of ordinary differential equations in time. The RHS can be interpreted as the slope of the value of the function in time. This method consists of four stages, including two predictor and two corrector steps, see Fig. 2.2.

First step: From these derivatives the value of the function in the next time step (marked •) can be calculated, picture taken from [START_REF] Press | Numerical Recipes in Fortran 77[END_REF].

k 1 = RHS(u)
u * = u n + ∆t 2 k 1 (2.12)
Second step:

k 2 = RHS(u * ) u * = u n + ∆t 2 k 2 (2.13)
Third step:

k 3 = RHS(u * ) u * = u n + ∆t k 3 (2.14)
Fourth step:

k 4 = RHS(u * ) u n+1 = u n + ∆t 6 (k 1 + 2k 2 + 2k 3 + k 4 ) (2.15)
Usually in addition to one memory location for u, five memory locations are necessary at each grid point for the evaluation of k 1 , k 2 , k 3 , k 4 and u * . However it is possible to release k 2 after the evaluation of u * in the third step, by simply adding the value of

k 2 to k 3 , say k3 ← k 3 + k 2
So it is possible to use again the allocated memory k 2 instead of k 4 during the computation. Hence the fourth step will change into k2 = RHS(u * )

u n+1 = u n + ∆t 6 (k 1 + 2 k3 + k2 ).

Adaptive time stepping with embedded Runge-Kutta methods

For time integration of the considered ODE, instead of using a fixed time step ∆t chosen a priori, it can be advantageous to have a technique that automatically adjusts its size dynamically to advance faster in time when possible and to guarantee the error in solution is bounded by a given value [START_REF] Domingues | An adaptive multiresolution method for parabolic PDEs with time-step control[END_REF]. The time step ∆t must be chosen sufficiently small to satisfy a required precision of the computed results, denoted by δ desired , but it must be sufficiently large to avoid unnecessary computational work. Typically, if u 1 is the approximation of u exact (t + ∆t) estimated by a method of order p, asymptotic developments of the local truncation errors of the form

u exact (t + ∆t) -u 1 = C 1 ∆t p + O(∆t p+1 ) (2.16)
can be used to find the step size required to attain a specified accuracy. However, since the leading constant C is not known a priori, practical error estimates are necessary.

To estimate the local truncation error, the idea is to apply two embedded ODE solvers, one with order p and the other with order p + 1. If u 2 is the approximation of u(t + ∆t) estimated by a method of order p + 1, then, for sufficiently small ∆t we have

u 1 -u 2 ≈ -C 1 ∆t p + C 2 ∆t p+1 ≈ -C 1 ∆t p (2.17)
this yields to estimate

C 1 ≈ u 2 -u 1 ∆t p = δ ∆t p (2.18)
By choosing a RK2 and RK3 method for approximations one gets

u * = u n + ∆t RHS(u n ) u RK2 = 1 2 (u n + u * + ∆t RHS(u * )) (2.19) and u * * = 1 4 (3u n + u * + ∆t RHS(u * )) u n+1 = u RK3 = 1 3 (u n + u * + 2∆t RHS(u * * )) (2.20) where δ = ||u p+1 -u p || ∞ = ||u RK3 -u RK2 || ∞
Hence the step size required to maintain the local truncation error of the first scheme below δ desired is given by ∆t new = ξ∆t old where

ξ = [ δ desired δ ] 1 p
For RK2/3 one chooses p = 2. To prevent the time step from varying too abruptly or to be sure that ∆t new in fact will produce an error less than δ desired , some care is needed. In this method, we cannot go back to the previous time step once the solution at the new time step is computed due to the low storage memory model we are using. Hence it is better to limit the increase of the time step by introducing a so-called safety factor S. The new time step ∆t new is chosen such that

ξ =    1 + S ξ > 1 + S ξ 1 -S > ξ > 1 + S 1 -S ξ < 1 -S (2.21)
where S depends to the problem and must be chosen in a way to not destabilize the time integration, e.g., S = 0.001. Moreover one memory location for u, four memory locations are needed in each grid point [START_REF] Domingues | Adaptive multiresolution methods[END_REF].

Multiresolution analysis

In the following, the multiresolution algorithm for a typical evolutionary PDE is briefly summarized. For more details on the general adaptive MR algorithm see [START_REF] Domingues | Adaptive multiresolution methods[END_REF].

Our multiresolution analysis is based on the cubic interpolating wavelet transform, described in some details in chapter 1. Because of the change in the grid points during each time step, the MR method is not applicable for backward multi step methods in time like Adams-Bashforth, instead, one must choose forward multi step methods like Runge-Kutta family for the time integration of the given PDE. First, depending on the regularity of the solution from previous time step (or the initial condition), we have the values (or the initial values) of an arbitrary dependent variable over a dyadic grid. According to the regularity of the solution, some points have an updated value, we call them active points. To perform a wavelet transform over the active points, we suppose the arrangement of these points over the dyadic grid is in such a manner that by starting from the finest level, down to the coarsest level, all necessary points are in the list of the active points. The wavelet transform is performed to compute the wavelet coefficients of these active points, which are also called details. Afterwards thresholding (nonlinear filtering) of the details (wavelet coefficients) is performed to determine the new active points (i.e., retained grid points are in a one-to-one correspondence with significant wavelet coefficients). Next a so-called safety zone, which will be described in more details further down, must be added to the new active points. As regards the values of some newly added points are unknown, an inverse wavelet transform will be performed just for the points that their values are unknown, that is to say, their values will be interpolated from known points in lower levels by performing an inverse wavelet transform with d j,i = 0. After that spatial derivatives in RHS must be evaluated with the use of the equations described in Section 2.3. Finally, a time evolution is made just for the active points.

The threshold operator (nonlinear filtering) depending on the choice of the threshold ϵ 0 , will select the important points for obtaining compression, while keeping the error introduced by filtering of wavelet coefficients (grid adaptation) below the truncation error of the spatial discretization scheme.

Safety zone

The presence of nonlinear terms in the governing equations of fluid mechanics will imply the creation of smaller scales in the flow, which may lead to energy transfer into the finer scales. Therefore when solving the governing evolutionary PDEs with MR technique, one needs to apply an additional criterion to ensure that the wavelet basis corresponding to a computational mesh is sufficient to approximate the solution throughout the time-integration stage for an evolution problem. To ensure the adequate approximation of the solution during time integration, Liandrat & Tchamitchian (1990) introduced the concept of a safety or an adjacent zone, which adds the points whose wavelet coefficients can possibly become significant during the period of time integration to the current significant points, when the grid remains unchanged. Most current wavelet-adaptation techniques are based on the concept of safety zone. In actual implementations, the safety zone (neighboring wavelets at the same and one above levels of resolution) is added for each significant coefficient. The sufficiently small time step usually chosen by explicit time integration methods, e.g., Runge-Kutta methods, ensures that the regularity of the solution does not propagate outside the safety zone. This typically results in a CFL-like (cell-Reynolds) constraint (Courant et al. 1928), which ensures that no energy at a given resolution scale propagates outside the safety region. The thickness of the safety zone determines the time interval during which the calculations can be carried out without modifying the computational grid. However, for computational efficiency, it was found (e.g., [START_REF] Schneider | Wavelets: Mathematical Theory[END_REF][START_REF] Vasilyev | Solving multi-dimensional evolution problems with localized structures using second generation wavelets[END_REF]) that the safety zone, which includes the immediate neighboring wavelets, is the optimal [START_REF] Schneider | Wavelet methods in computational fluid dynamics[END_REF]. The adaptation strategy for an evolution problem is illustrated in Fig. 2.3. In our simulations the number of points in the safety zone is almost equal to the number of active points.

As explained to allow the basis (grid points) to change, we have to extend the sparse point representation after thresholding. First, in space, we add neighbor points corresponding to neighbor wavelets in the sparse point representation. The number of points to be added depends also on the PDEs wave propagation speed. Then the grid points are able to adjust when the solution is moving in space. Second, in scale, we add neighbor points corresponding to wavelets on the next finer scale in the sparse point representation. The number of points to be added in upper scales depends to order of nonlinear term in the governing PDE, e.g., second order in the case of Navier-Stokes equations and third order in the case of nonlinear Schrödinger equation (NLS, discovered in optics and water waves). This refinement allows the accurate development of the solution in time and space, i.e., four points in one dimension and 16 points in two dimensions must be added but most of them will be common [START_REF] Holmstrom | Solving hyperbolic PDEs using interpolation wavelets[END_REF]. Moreover to the points added in the above procedure another check must be performed to add all necessary points to the current list of the active points for having a consistent direct or inverse WT. This depends on the stencil of the WT, we have presented two methods of interpolation in previous chapter, i.e., linear and cubic interpolation. In one dimension two points for linear interpolation and four points for cubic interpolation are necessary. For two-dimensions depending on the index The second bell-shaped region (---) on the right of the original one marks the locations of the significant wavelet coefficients at the end of the time integration, i.e., at t + ∆t. The safety zone is represented by red bell-shaped region, picture taken from [START_REF] Schneider | Wavelet methods in computational fluid dynamics[END_REF].

composition (odd-even), different situations can occur, see Fig. 1.6. This step will determine all the points that will be used in the current time step to evaluate the solution.

Another point to be considered is the possibility to have, big jumps in the stencil of central finite differences in an adaptive method which may lead to a divergence of the value of the derivatives in some points. To prevent this problem, in each solution point the ratio of the distances of right and left neighbor points must not be more than two. In some literature this is called gradedness of the grid points and maybe it is equivalent to the check for wavelet interpolation stencil of WT based on linear interpolation. Of course the stencil of WT based on cubic interpolation is wider than linear interpolation and so the gradedness is guaranteed. This point will be optional if necessary points for interpolation via WT are already added to the active points.

Multiresolution analysis in summary

Denoting by E(∆t) the discrete time evolution operator, the global algorithm can schematically be summarized by

u n+1 = E(∆t) [ M -1 • S • T (ϵ) • M ] u n (2.22)
where M and M -1 are the direct (WT) and inverse (IWT) wavelet transform operators. T (ϵ) is the thresholding operator and S represents the safety zone operator.

For an Euler explicit time integration we have

E(∆t)u n = u n + ∆t RHS(u n ).
where RHS operator contains all the terms of the considered evolutionary equation to be integrated except time derivative. In summary the algorithm is given in the following:

1. Start from an initial condition over a dyadic grid 2. Apply WT to the active points (from the finest level down to the coarsest level) to compute the wavelet coefficients of the independent variable 3. Perform thresholding T (ϵ) to remove all the points from the list of the active points which their wavelet coefficients are below the corresponding threshold ϵ j

N.B. 3

In the case of the two-dimensional Navier-Stokes equations in velocityvorticity formulation, before calculation of the spatial derivatives it is necessary to solve an elliptic equation, i.e., Eq. (2.27) for updating stream function.

Extension to two-dimensional incompressible Navier-Stokes equations

The governing equations of incompressible flows are the Navier-Stokes equations (2.23 -2.24). In primitive variables momentum equations read

∂u ∂t + (u • ∇)u = - 1 ρ f ∇p + ν∇ 2 u + F (x, t) ∈ Ω × [0, ∞) (2.23)
and the continuity equation corresponds to

∇ • u = 0 , x ∈ Ω (2.24)
where Ω is the spatial domain of interest, given as an open subset of R 3 , which can be bounded or unbounded in general, u(x, t) is the velocity field which assumed to be at least in C 2 with respect to space and in C 1 with respect to time, p(x, t) is the pressure field, ν = µ/ρ f is the kinematic viscosity, ρ f is the fluid density and F(x, t) is a forcing term. Thus, for a complete description of a particular problem, the above equations need to be completed to describe an initial/boundary value problem (IBVP). Hence by specifying an initial condition u(x, 0) = u 0 (x), which we assume to be in C ∞ and divergence free in all of Ω and by giving boundary conditions for velocity and pressure satisfying the global mass conservation constraint one will seek the solution during time evolution [START_REF] Rempfer | On boundary conditions for incompressible Navier-Stokes problems[END_REF]. By choosing respectively U ∞ and L as reference velocity and length for a given problem the Navier-Stokes equations can be written in nondimensional form in which Re = U ∞ L/ν is the Reynolds number and ρ is set to unity:

∂u ∂t + (u • ∇)u = -∇p + 1 Re ∇ 2 u + F (2.25)
However in two-dimensional problems as explained earlier the vorticity and streamfunction formulation has the advantage that it not only eliminates the pressure variable entirely, but also ensures a divergence-free velocity field (continuity equation), if Eq. (2.27) properly satisfied. One will encounter with two scalar quantity, i.e., ψ and ω, instead of velocity vector and the pressure field, thus it makes the computations very efficient. We continue with this formulation, but the concepts are applicable also to primitive variable formulation. By taking the curl of Eq. (2.23) after elimination of terms due to two-dimensional assumption, incompressibility Eq. (2.24) and simplifications due to constant density, one obtains the vorticity transport equation for two-dimensional flows

∂ t ω + (u • ∇)ω = ν∇ 2 ω + ∇ × F (2.26)
where ω = ∇×u = v x -u y denotes the vorticity. The equation is a parabolic equation in time and the velocity components are determined from u = (u, v) = (∂ y ψ, -∂ x ψ) with ψ being the stream function, satisfying

-∇ 2 ψ = ω (2.27)
which is an elliptic

Boundary conditions

The boundary conditions for a curved boundary denoted by (s) moving tangent to its surface with a constant velocity U tan can be written in terms of the stream function ψ, at each boundary section Γ i . The no-penetration boundary condition is equivalent to

∂ψ ∂ τ wall = n • u(s, t) = 0 (Neumann) ↔ ψ| wall = C i (Dirichlet) (2.29)
and the no-slip boundary condition read

∂ψ ∂ n wall = -τ • u(s, t) = ±U tan (Neumann) (2.30)
where C i is a constant for each (Γ i ), τ is the tangent to the wall direction and n is the direction normal to the wall. Free-slip boundary condition is evident and will achieved by ω = 0. For a wall with zero tangential velocity we will have U tan = 0. In a simply-connected domain, C 0 can be taken equal to zero (C 0 = 0). As a result for a fixed horizontal wall the no-penetration and no-slip boundary conditions are (v = u = 0) or (ψ x = ψ y = 0). For simulation of two dimensional incompressible flows, beside the mentioned advantages of vorticity and stream-function formulation one can see some disadvantages. Yet, the main difficulties in the numerical implementation of this formulation come from the boundary conditions [START_REF] Wang | Analysis of finite difference schemes for unsteady Navier-Stokes equations in vorticity formulation[END_REF] the majors among them are as follows:

1. The implementation of the two cited boundary conditions for the stream function ψ simultaneously.

2. For updating the vorticity ω during time integration, there is no definite boundary condition for vorticity.

3. Determining the constants C i at each boundary of "holes" Γ i if the computational domain is multi-connected.

However, several methods have been proposed to update the vorticity boundary condition over a solid boundary moving tangent to its surface with a constant velocity U tan . In the following we cite some of them, which will lead to second and fourth order accuracies :

Thom's formula for the vorticity at wall

ω i,0 = ψ i,0 -ψ i,1 ± hU tan 0.5h 2 + O(h 2 ) (2.31)
Wilkes, Roach, Pearson and Jensen's formula for the vorticity at wall

ω i,0 = 7ψ i,0 -8ψ i,1 + ψ i,2 ± 6hU tan 2h 2 + O(h 2 ) (2.32)
Briley's fourth order formula for the vorticity at wall

ω i,0 = 85ψ i,0 -108ψ i,1 + 27ψ i,2 -4ψ i,3 ± 66hU tan 18h 2 + O(h 4 ) (2.33)
Other relations were proposed by Woods (1954) and Orszag and Israeli (1974). It is very important to know that, the vorticity boundary condition is responsible to enforce no-slip boundary condition. Although non of them cannot fall velocity components to machine zero, but the accuracy in normal to the wall component of the velocity is two order of magnitude more than tangent component. The subject of the vorticity boundary condition has a long history, going back to Thom's formula in 1933 [START_REF] Thom | The flow past circular cylinders at low speeds[END_REF]. In a second order scheme, Thom's formula, Wilkes formula, or some other local formulas can be selected and coupled with a centered difference scheme at the interior points. The advantage of Thom's formula lies in its simplicity as only one interior point of the stream function is involved. Thus its stability is very robust. Yet, it was always very confusing why formulas like Thom's, which seems hopelessly first order by formal Taylor expansion on the boundary, is actually second order accurate. This mystery can be explained by Strang-type high-order expansions [START_REF] Hou | Convergence of a finite difference scheme for the Navier-Stokes equations using vorticity boundary conditions[END_REF]. It was proven in 1964 that for nonlinear hyperbolic or parabolic equations, the L 2 -norm stability for the linearized problem and the smoothness of the exact solution implies that the scheme has full accuracy (e.g., second order) in L ∞ -norm. The main idea in the proof is the construction of high order expansions with respect to the scheme. Meth in his thesis [START_REF] Meth | A vortex and finite difference hybrid method to compute the flow of an incompressible, inviscid fluid past a semi-infinite plate[END_REF], proved the stability for the linearized problem. The theoretical convergence analysis of the second order scheme with Thom's formula on the boundary was given by Hou and Wetton in [START_REF] Hou | Stable fourth-order stream-function methods for incompressible flows with boundaries[END_REF]. It relied upon Strang-type high order expansions, which resulted in much more regularity assumption of the exact solution than needed. We should mention that the stability of Thom's formula cannot be applied to long-stencil formulas automatically, as Wilkes (1963) & Pearson's (1965) formula and other local formulas. It was a doubtful question for a long time: are these long-stencil formulas stable? This question is answered in [START_REF] Wang | Analysis of finite difference schemes for unsteady Navier-Stokes equations in vorticity formulation[END_REF] by performing a simple, clean analysis of the second order scheme using Wilkes formula to determine the vorticity at the boundary. In fact, direct calculations and standard local estimates cannot work it out. The convergence analysis of the compact scheme with Briley's [START_REF] Briley | A numerical study of laminar separation bubbles using the Navier-Stokes equations[END_REF] formula is given in [START_REF] Wang | Analysis of finite difference schemes for unsteady Navier-Stokes equations in vorticity formulation[END_REF], illustrating why the fourth order scheme combined with the third order boundary formula still works. Generally one can say that the order of implementation of the Neumann boundary conditions may be one order less than that of the discretization in the interior points of the solution domain. This choice can keep the global desired accuracy of the discretization method. In the present investigation, Thom's formula at the no-slip boundaries is used for updating the vorticity. For more discussions, on boundary conditions for incompressible Navier-Stokes problems see [START_REF] Rempfer | On boundary conditions for incompressible Navier-Stokes problems[END_REF].

Volume penalization method

For the simulation of curved solid boundaries which do usually not coincide with grid points one can use the volume penalization method proposed by Arquis and Caltagirone [START_REF] Arquis | Sur les conditions hydrodynamiques au voisinage d'une interface milieu fluid -milieux poreux: application à la convection naturelle[END_REF]. It is based on the idea of modeling solid bodies as porous media, thus getting ride of the Dirichlet boundary conditions by considering both the fluid and the solid part as one domain with different permeabilities, so one has a domain in which the solid is embedded. This method will lead to first-order accuracy near the solid boundary. In the Navier-Stokes equations (2.23) in primitive variables or the vorticity transport equation (2.26), penalization term can be added as a force term and thus, it is possible to introduce a solid body in the flow field. The penalization term for unit mass of the fluid reads,

F = - χ η (u -u B ) (2.34)
where u B is the velocity of the immersed body (obstacle) which will be zero for fixed bodies. Penalization parameter η is the porosity (permeability) coefficient of the immersed body. In an explicit time integration ∆t must be smaller than η. Typically we use η = 10 -3 or 10 -4 to ensure the stability of time integration. The mask function χ describes the geometry of the flow domain, see Fig. 2.4

χ(x, t) = { 1 x ∈ Ω s 0 x ∈ Ω f (2.35)
where Ω f represents the domain of the flow and Ω s represents the immersed solid in the domain of the solution. The domain is governed by Navier-Stokes equations in the fluid regions and by Darcy's law in the solid regions, when η → 0. With the use of the penalization method the hydrodynamic forces and the moments acting on the obstacle, which are usually evaluated via surface integrals of the stress tensor σ(u, p) = µ(∇u + (∇u) T )/2 -p I, can computed readily by integrating the penalized velocity over the obstacle's volume, thus we have the forces by

F = ∂Ωs σ • n dl = lim η→0 ρ f η ∫ Ωs χ(u -u B ) ds + ρ f S pen ẍcg (2.36)
and moment in two-dimensions is evaluated by

M xcg = ∂Ωs (x-x cg )×σ •n dl = lim η→0 ρ f η ∫ Ωs χ(x-x cg )×(u-u B ) ds+ ρ f ρ s J cg θcg (2.37) in [N.m],
where J cg is the moment of inertia taken about center of the mass, n is the unit outward vector to ∂Ω s , x cg is the location of the center of gravity of the immersed body, θ is the angle of rotation around the center of gravity (dots denote derivation with respect to time) and S pen is equivalent to the surface of penalized area

S pen = ∫ Ωs χ ds
For a uniform grid in two-dimensions numerically F can be evaluated as

F ≈ ρ f η ∆x∆y Imax ∑ i=1 Jmax ∑ j=1 χ i,j (u -u B ) i,j + ρ f S pen ẍn-1 cg (2.38)
Finally the RHS of the vorticity equation including penalization term reads

RHS = - ∂ψ ∂y ∂ω ∂x + ∂ψ ∂x ∂ω ∂y + ν ( ∂ 2 ω ∂x 2 + ∂ 2 ω ∂y 2 ) + ∇ × F (2.39)
with the penalization term

∇ × F = ∂ ∂y ( χ η ∂ψ ∂y ) + ∂ ∂x ( χ η ∂ψ ∂x ) - ∂ ∂y ( χu B η ) + ∂ ∂x ( χv B η ) (2.40)
For the points for which the complete stencil belongs to the fluid domain (χ = 0) we have ∇ × F = 0, but for the points for which the complete stencil belongs to the solid domain (χ = 1) we have ∇ × F = (ω B -ω)/η, i.e., in the time integration we will have ω n+1 ≈ ω B = 2Ω B , where Ω B is the angular velocity of the embedded body. The penalization term is thus responsible for the vorticity production at the boundaries.

Chapter 3 Applications

Introduction

In this chapter the results obtained by applying the developed adaptive multiresolution method to the one-dimensional viscous Burgers equation and the two-dimensional incompressible Navier-Stokes equations will be presented. A dipole colliding with a straight no-slip wall for a Reynolds of 1000 will be considered as a benchmark, then with the use of the penalization method the collision of dipoles with curved walls, either concave or convex, will be examined for Reynolds= 10000.

Application to the one-dimensional Burgers equation

A well-known example of nonlinear PDEs is the one-dimensional viscous Burgers equation, including both nonlinear advection and diffusion terms.

∂ t u + uu x = νu xx (3.1)
The value of ν determines the ratio of two phenomena in competition. The Eq. (3.1) is completed with suitable initial and boundary conditions. By putting the spatial terms in RHS, using second-order centered finite differences for spatial discretization, the RHS operator reads

RHS(i) = -u i u i+d -u i-d 2d∆x + u i+d -2u i + u i-d d 2 ∆x 2 (3.2)
where for a uniform grid d will be replaced by d = 1. First a solver based on uniform discretization of the Burgers equation via finite differences using Eq. (3.2) is developed as reference solver. An error analysis of the time integration of the Burgers equation (3.1) over a uniform grid is performed to examine the accuracy of the developed code based on the fourth order Runge-Kutta method. The aim is to show the decrease of the error with successive reduction of the time step. A simulation with ∆t max /16 is considered as reference solution to compute the errors, where ∆t max is chosen rather For x ∈ [0, 2] it can be verified that the boundaries are sufficiently far so that their influence are negligible [START_REF] Roussel | Développement d'un algorithme multirésolution adaptatif tridimensionnel pour les équations aux dérivées partielles paraboliques. Application l'étude des instabilités thermo-diffusives de flamme[END_REF] for small time intervals. The initial condition is given by

u(x, 0) = { 1 x ≤ 1 0 x > 1 (3.4)
and Dirichlet boundary conditions at the right and left boundaries are imposed, i.e., u(0, t) = 1 , u(2, t) = 0.

The initial condition and the time evolution of the solution computed with the multiresolution solver for Burgers equation, with threshold ϵ = 10 -3 for ν = 10 -3 and ∆t = 10 -5 , maximum grid level J = 11, are illustrated in Fig. 3.2 (a) and (b), respectively. The exact solution is plotted simultaneously by solid lines. The L ∞ , L 1 and L 2 errors at t = 1.5 are of order 10 -2 , 10 -4 and 10 -5 , respectively. The evolution of the active grid points during the computation, is illustrated in Fig. 3.3 and corresponding active points represented in their level, is illustrated in Fig. 3.4.

With the use of multiresolution code for the Burgers equation (3.1), some qualitatively comparisons of the results with uniform grid solver was performed. The only adjustable parameter which had a considerable impact in the agreement of the results with that of uniform grid solver is the threshold parameter (ϵ). Another error analysis is performed to examine the second-order accuracy of the multiresolution method in space by changing the threshold parameter. To show the conservation of the formal accuracy for the space discretization with adaptive method, the analytical exact solution in Eq. (3.3), is used as a reference solution to compute the errors, different simulations are performed with 2 11 , 2 10 , 2 9 , 2 8 and 2 7 grid points. The computations start from the initial condition showed in Fig. 3.2 (a) at t = 0 and stop at t = 1.5 showed in Fig. 3.3. The time step is chosen smaller than stability limit for 2 11 grid points, to minimize the error due to time integration. Hence the number of iterations for all simulations are the same. Different errors are compared with the theoretical (-2) slope in Fig. 3.5 (a) and (b) for ϵ = 10 -3 and ϵ = 5 × 10 -3 , respectively. The results are in good agreement with the (-2) slope, which confirms that the adaptive multiresolution method preserves the discretization order of the underlying numerical scheme. There is a deviation for ϵ = 5 × 10 -3 , with 2 11 grid points, which is due to big threshold parameter in comparison with truncation error. The truncation error of a second-order method in space x ∈ [0, 2π] is of order 10 -5 for 2 11 grid points. The evolution of the number of active, important (corresponding to the retained points after filtering of wavelet coefficients), safety zone, hung and interpolated points for wavelet transform during the computation with the multiresolution solver for the Burgers equation with ν = 10 -3 and ∆t = 10 -5 with maximum grid level J = 10, for threshold ϵ = 10 -3 and ϵ = 5 × 10 -3 , are illustrated in Fig. 3.6, (a) and (b), respectively. The average number of points after filtering is about 35 and after adding the safety zone, it increases to 63 out of 1025. Hence for this test case we have a compression of points about 95% which is considerable, and the time of computations with the multiresolution solver is at least half in comparison to the uniform grid solver. The gain in CPU-time becomes more especially for fine resolutions, where the overhead due to the additional complexity of the multiresolution algorithm becomes less important. The evolution of the number of active, important (corresponding to the retained points after filtering of wavelet coefficients), safety zone, hung and interpolated points for the wavelet transform during the computation with the multiresolution solver for the Burgers equation, with maximum grid level J = 10 (1025 grid points).

Iteration

Number of Points

Application to two-dimensional incompressible flows

In this section we are considering the two-dimensional incompressible flows by means of the vorticity transport equation (2.26). By putting all of the spatial derivatives in the RHS operator, the discretized form of the vorticity equation including the penalization term (2.34) reads

∂ t ω| i,j = RHS(i, j)
with the use of second-order finite-difference approximations (2.6), (2.7) and (2.8) one gets:

RHS(i, j) = - ψ i,j+d -ψ i,j-d 2d∆y ω i+d,j -ω i-d,j 2d∆x + ψ i+d,j -ψ i-d,j 2d∆x
ω i,j+d -ω i,j-d 2d∆y +ν ( ω i+d,j -2ω i,j + ω i-d,j d 2 ∆x 2 + ω i,j+d -2ω i,j + ω i,j-d d 2 ∆y 2 ) + (χ i+d,j + χ i,j )(ψ i+d,j -ψ i,j ) -(χ i-d,j + χ i,j )(ψ i,j -ψ i-d,j ) 2ηd 2 ∆x 2 + (χ i,j+d + χ i,j )(ψ i,j+d -ψ i,j ) -(χ i,j-d + χ i,j )(ψ i,j -ψ i,j-d ) 2ηd 2 ∆y 2 - (u B i,j+d χ i,j+d -u B i,j-d χ i,j-d ) 2ηd∆y + (v B i+d,j χ i+d,j -v B i-d,j χ i-d,j ) 2ηd∆x
For a uniform grid in the above formula d will be replaced by 1, for the multiresolution solver d will be determined by a search algorithm, to find the nearest active point.

The vorticity equation is advanced in time with a Runge-Kutta technique which was described in Section 2.4.1. Before each evaluation of RHS, one Poisson equation for stream function with vorticity (ω) forcing on the right hand side, i.e., Eq. (2.27) must be solved to update the stream function to the correct value for the interior points. The vorticity at the boundary points must be updated by Thom's formula. The Poisson equation will be solved via an iterative method for elliptic equations, namely the point successive over relaxation (PSOR) method

ψ n+1 i,j = βψ new i,j + (1 -β)ψ old i,j (3.5)
where β is the over relaxation factor, β ∈ [START_REF] Thom | The flow past circular cylinders at low speeds[END_REF][START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF]. So we will have

ψ n+1 i,j = ψ old i,j + β 4 Res (3.6)
by the use of the Gauss-Seidel method

Res = 4 ∑ nb=1 ψ last nb -4ψ old i,j + h 2 ω i,j (3.7)
in which h = ∆x = ∆y and nb denotes the neighbor points with distance d (always their updated value will be used). The norm of the residuum, ||Res|| ∞ must converge below to a prescribed convergence criterion, e.g., 10 -6 . For the finer grids like J = 12 we require smaller criterion, e.g., 10 -12 . For useful discussions in the choice of a proper over relaxation factor we refer to [START_REF] Press | Numerical Recipes in Fortran 77[END_REF]. The hung points are the same for the elliptic and the parabolic part and they must be updated in each iteration. In the present investigation a geometric multi-grid solver for uniform grids will be used to accelerate the convergence of the elliptic part, which is the most time-consuming part of the simulations. But the PSOR solver is used in the multiresolution simulations because for the moment the MG solver is not applicable for adaptive grids. See [START_REF] Vasilyev | An adaptive multilevel wavelet collocation method for elliptic problems[END_REF] for a multilevel adaptive method for elliptic problems and [START_REF] Ding | A stencil adaptive algorithm for finite difference solution of incompressible viscous flows[END_REF] for a MG method used in adaptive solution of incompressible flows via FDM. In terms of memory usage the PSOR solver is more advantageous in comparison to MG solver, thus for the moment the applicability of multigrid solver due to insufficient memory is limited to maximum grid level J = 10 in each direction.

Dipole-wall collision

The famous problem of the dipole collision with a straight wall is considered in this section. We will study the flow in a square domain [0, 2] × [0, 2]. On the four walls of the domain (x, y = 0 & x, y = 2) no-slip and no-penetration boundary conditions are applied. The flow is initialized in the form of two shielded Gaussian mono-polar vortices, their centers placed at a distance 0.2 apart. The vorticity distribution in each monopole is given by

ω(0, x n ) = ω e ( 1 - r 2 r 2 0 ) exp ( - r 2 r 2 0 ) (3.8)
where r 0 is the core radius, r = ||x -x n || with x n being the position of the vortex center. The two isolated monopoles are located at

x 1 = (1, 1.1) and x 2 = (1, 0.9)

Demanding that the root mean square (rms) velocity is initially equal to unity (E = 2) yields the amplitude of each isolated monopole, ω e = ±299.528385375226 [START_REF] Keetels | Fourier spectral and wavelet solvers for the incompressible Navier-Stokes equations with volume-penalization: Convergence of a dipole-wall collision[END_REF]. The core radius of the shielded monopoles is set to r 0 = 0.1. The vorticity amplitude in the radial direction decreases exponentially with r. As a result, the circulation of one isolated monopole calculated over a circular contour around the vortex origin decreases exponentially towards zero for increasing contour radius. Hence, no boundary layers are required at the no-slip walls when constructing the initial flow field. The integral-scale Reynolds number for the initial field is given by

Re = U rms L ν (3.9)
where the characteristic length scale is set to the half-height of the domain, L = 1, and the characteristic velocity to the initial root mean square velocity, U rms = 1. This integral Reynolds number differs slightly from the Reynolds number Re d ≈ 0.8Re based on the dipole translation speed U d and the dipole radius R, see [START_REF] Kramer | Vorticity dynamics of a dipole colliding with a no-slip wall[END_REF]. The time evolution of the dipole is calculated by the developed multiresolution finite difference code with threshold ϵ = 10 -3 and maximum grid level J = 11 for Reynolds 1000. The evolution of the vorticity isolines starting from the initial condition at t = 0 up to t = 1, is shown in Fig. 3.7 and Fig. 3.8. After releasing the initial set of shielded mono-polar vortices, Fig. 3.7 (a), one observes how the cores combine into a dipole, while the surrounding shields are substantially deformed Fig. 3.7 (b). Somewhat later, the dipolar vortex has traveled away in the positive x-direction Fig. 3.7 (c), and the shields left behind combine into another, much weaker dipolar structure that translates in the opposite direction Fig. 3.7 (d). This weak dipole has no considerable impact on the primary dipole and is not discussed further. The first collision with the right wall can be seen in Fig. 3.7 (e) and the formation of secondary vortices is illustrated in Fig. 3.7 (f). Previous studies have revealed that the primary dipole vortex closely resembles the Lamb dipole, i.e., a compact dipolar vorticity structure in a circular area with a linear vorticity stream-function relationship. The rest of the evolution and the collision of the primary and newly generated vortices with right and left walls are illustrated in Fig. 3.8; detachment of vortices from wall is seen in Fig. 3.8 (a), a second collision of newly generated dipoles with the right wall in Fig. and ϵ = 10 -4 , are illustrated in Fig. 3.11, (a) and (b), respectively. The average number of points after filtering is about 3 × 10 5 and after adding the safety zone it is increased to 10 6 . Hence for this test case we have a compression of points about 25% which is not so good in comparison to one-dimensional test case, but the time of the computations with multiresolution solver is at least one over six in comparison to the uniform grid solver. This CPU-time is promising, especially for higher resolutions with large number of grid points, i.e., J larger than nine. The main reason for considerable CPU-reduction is the presence of the elliptic solver in the algorithm in which the order of operations is proportional to O(N 2 ) for a point successive over relaxation method, where N is the number of active points.

Convergence study*

To verify the accuracy of the numerical method the results are compared with the results reported by Clercx et al. [START_REF] Clercx | The normal and oblique collision of a dipole with a no-slip boundary[END_REF] which are computed with a pseudo spectral solver. The method of Clercx et al. is a pseudo spectral method, both the velocity and vorticity are expanded in a truncated series of Fourier polynomials for the periodicdirection and in a truncated series of Chebyshev polynomials for the non periodicdirection. Some invariants of the flow, (i.e. total energy, enstrophy and palinstrophy) which are conserved by the flow dynamics for inviscid fluids (ν = 0), can be assessed in viscous flows, where they will not conserved, but instead varying in time depending on the Reynolds number. Three quantities in the flow field, i.e., total energy E, total enstrophy Z and total palinstrophy*1 P are defined as

E(t) = 1 2 ∫ Ω |u(x, t)| 2 dx (3.10) Z(t) = 1 2 ∫ Ω |ω(x, t)| 2 dx (3.11) * P (t) = 1 2 ∫ Ω |∇ω(x, t)| 2 dx (3.12)
Applying the trapezoidal quadrature formula we get the discrete versions

E(t) = ∆x∆y 2 Imax ∑ i=1 Jmax ∑ j=1 (u 2 i,j + v 2 i,j ) (3.13) Z(t) = ∆x∆y 2 Imax ∑ i=1 Jmax ∑ j=1 (ω i,j ) 2 (3.14) * P (t) = ∆x∆y 2 Imax ∑ i=1 Jmax ∑ j=1 ( ∂ω i,j ∂x
) 2 + ( ∂ω i,j ∂y
) 2 (3.15)
where ∆x = Lx Imax-1 and ∆y = Ly Jmax-1 . For any two-dimensional viscous flow (z > 0) the total energy E(t) decays according to

dE dt = -ν ∫ Ω ω 2 dA = -2νZ. (3.16)
where -2νZ is the energy dissipation [START_REF] Kraichnan | Two-dimensional turbulence[END_REF]. Note that the decay rate is proportional to the total enstrophy, Z(t), which is a measure of the squared vorticity integrated over the domain. Understanding the evolution of the total enstrophy is therefore of crucial importance for explaining the energy decay. For a domain with no-slip boundaries the change in total enstrophy is governed by

dZ dt = -2νP + ν ∂Ω ω(n • ∇ω)ds (3.17)
where n denotes the outer normal vector with respect to ∂Ω. The first term on the right-hand side simply states that the enstrophy decays due to vorticity gradients (palinstrophy) that are present in the flow. The second term represents the vorticity production at the no-slip boundaries involving the vorticity and its gradients, which will give rise to the total palinstrophy. Note that the vorticity influx at the no-slip boundaries is equal to (n • ∇ω). In the case of a square domain with stress-free or periodic boundary conditions the second term on the right-hand side of Eq. (3.17) vanishes. As a result, the total enstrophy cannot increase for a domain with stressfree or periodic boundary conditions and is thus always bounded by its initial value and zero [START_REF] Kramer | Vorticity dynamics of a dipole colliding with a no-slip wall[END_REF]. For a steady flow we have

P = 1 2 ∂Ω ω(n • ∇ω)ds
Different simulations were performed to obtain a grid independent solution, by successive increasing the number of points. To have a stable simulation the time step must be reduced according to the CFL condition, see Table 3.1. A simulation over an uniform grid with a second order (in space) multi-grid solver with maximum grid level J = 10 in each direction, was done in the first step. The evolution of the total kinetic energy, total enstrophy and total palinstrophy for Re = 1000 are compared with the computations of Clercx et al. [START_REF] Clercx | The normal and oblique collision of a dipole with a no-slip boundary[END_REF] in Fig. 3.12 and Fig. 3.13 (a) and (b), respectively. Note that the energy steadily decreases from its normalized initial value of E = 2 towards E ≈ 0.8 at t = 1. At t ≈ 0.35 the first collision of the dipole with the right wall takes place and the kinetic energy decays faster, which is due to the increased enstrophy production (dissipation) on the domain. On the other hand enstrophy starts from its initial value P (0) = 800 and decays until t = 0.2 and then goes up rapidly. The first peak in the enstrophy curve during time takes places at t = 0.35, and thus coincides with the first collision of the dipole with the right wall.

During the first collision the boundary layers create a large amount of vorticity. The enstrophy in the boundary layers is then the main contribution to the total enstrophy. At t ≈ 0.64 another smaller peak is visible in the enstrophy evolution curve, which is due to the second collision of newly generated vortices. The results for total energy are in reasonable agreement with the pseudo-spectral simulations, we can see a systematic deviation in enstrophy and palinstrophy curves especially near the first peak, which is due to insufficient resolution of second-order solver in comparison with the pseudo-spectral solver of Clercx [START_REF] Clercx | Dissipation of kinetic energy in twodimensional bounded flows[END_REF].

Two simulations for Reynolds 1000 with the multiresolution code, using a maximum grid level J = 11, were performed with ϵ = 10 -3 and ϵ = 10 -ulation with pseudo-spectral solver of Clercx is taken as the reference solution. The results of the present computation are illustrated in Fig. 3.16, (a) and (b), respectively. It can be observed that by increasing the number of grid points the curves become closer and closer, we hope the results of J = 12 will match with that of Clercx. Now comparisons of total energy E(t), total enstrophy Z(t) and total palinstrophy P (t) between the uniform grid solver and the multiresolution computation with thresholds, ϵ = 10 -3 and ϵ = 10 -4 , for Reynolds 1000 and with maximum grid level J = 9 will be presented. The comparison of total energy E(t) is plotted in Fig. 3.17, the results for total enstrophy Z(t) and total palinstrophy P (t) are compared in Fig. 3.18, (a) and (b), respectively. The agreement between the uniform grid solver and the multiresolution solver is perfect and the results for ϵ = 10 -3 and ϵ = 10 -4 are almost identical.

The effect of curvature

As mentioned in the previous section the curved walls in a Cartesian grid will take into account with the penalization method. For validation of the method the same benchmark problem, i.e., dipole collision with straight wall, will be considered in a penalized cavity to compare the results with classical boundary conditions (Dirichlet and/or Neumann) as explained in Section 2.6. All of the computations in this section are performed with the uniform multigrid solver because of its efficiency in terms of CPU-time. In the first step, for Re = 1000, the simulations will be performed for quantitative comparisons of total energy E(t), total enstrophy Z(t) and total palinstrophy P (t). The comparisons will be in a same geometry, with the same initial condition. Two category for simulation of a fixed wall will be considered; first the classical no-slip an no-penetration conditions in boundaries of Cartesian grid, second, the velocity in a marginal area of a slightly larger cavity will imposed to zero with penalization method. A schematic representation of the penalized cavity, i.e., the geometry and initial condition are illustrated in Fig. 3.19 (a) and (b), respectively. Comparisons of total energy, total enstrophy and total palinstrophy between the two cited methods for simulation of a rigid wall for Re = 1000 with a maximum grid level of J = 10 in each direction are illustrated in Fig. 3.20, Fig. 3.21 (a) and (b), respectively. The small deviation in the total energy, which is almost invisible, is negligible and below 3%. For both total enstrophy and total palinstrophy a considerable deviation can be observed, which can be explained by the first order nature of implementing a solid wall via the volume penalization method. So one must use a very fine grid for obtaining more accurate results with this method.

After the procedure of validation, for considering the effect of curvature we will continue with Reynolds 10000 because the flow field is more sensitive to the curvature of the hitting wall. Thus this Reynolds number will be used for studying the effect of curvature in the collision of a dipole with concave and convex walls, even if a grid with 1025 2 points is not sufficient to fully resolved all the spatial scales who are present in the flow. But the numerical dissipation of the second-order central method will compensate the effect of eliminated scales to have a stable simulation. For having a reference solution to compare the results of the penalization method for curved walls, first the dipole collision with a straight penalized wall at Reynolds 10000 will be simulated. A schematic representation of the penalized cavity, i.e., the geometry and the initial condition are illustrated in Fig. 3.22 (a) and (b), respectively. To obtain a qualitative idea about the flow and its future destiny the evolution and collision of dipole with a straight penalized wall from t = 0 to t = 2 for Reynolds 10000 is illustrated in Fig. 3.23 and Fig. 3.24. The simulations are performed with a multigrid solver over a uniform grid with maximum level J = 10 in each direction. To study the effect of curvature two configurations will be considered. First the collision of a dipole with a convex wall with radius R = 1.5 where the center of the convex wall is placed at x c = 3.1, y c = 1.0. The initial vortices are placed at distance equal to one from the right wall in all of the simulations. Second the collision with a concave the collision with the walls, this is also independent of the curvature of the hitting wall. The destiny of the small newly generated dipole reflecting from the right wall after the second collision of the dipoles with right wall will be different when arriving to the left wall. One can see the division of this dipole into two vortices in the concave and convex test cases, while it will remain attached in the case of straight wall because it is stronger. At t = 2 in the case of collision with a straight wall the number of vortices are about n = 13, while in the cases of concave or convex walls they are about n = 8. One can see different phenomena in the field such as attachment to and detachment of vortices from the walls, merging of counter-rotating vortices, creation and dissipation of vortex sheets (in two dimensions they will be filaments), creation of new small vortices near the wall with high rotational speed and dissipation of the vortices during the time. Also collision of vortices and dipoles of different sizes with the walls and themselves, interaction of vortices, new dipoles or new isolated vortices moving with different self induced velocities, the effort of vortices to become circular and symmetric (isotropy) especially in the regions far from the wall to minimize the potential. The creation of anisotropy near solid walls can also be seen. Another point to be noted is that the vortices are very persistent and their life time is much longer in comparison to three-dimensional flows, in which the vortex stretching mechanism is responsible for the fast decay of vortices by decreasing their size and increasing the vorticity in the center line of the vortex tubes thus leading to the direct energy cascade. In two-dimensional flows, the dominant mechanism to destroy a vortex is the dissipation which becomes small for high Reynolds numbers. The simulations also exhibit a symmetry of the flow with respect to the y = 1 axis, up to t=4. Then the symmetry is broken due to the instability of the flow. with respect to the computations on a uniform grid, taking into account the overhead necessary for the multiresolution algorithm. For higher resolutions further reduction in CPU-time are expected. The average number of active points over the base grid for two-dimensional problems is higher than that of one-dimensional problems due to safety zone, which is necessary to account for the translation of the solution and the generation of finer scales. In two dimensions 16 points for each active point must be added as safety zone, four at the same and four at one finer scale, while in one space dimension only 4 points are necessary as safety zone, two at the same and two at one finer scale. For the two-dimensional computations the compression rate is typically about 25%, which is thus not so high. However, the CPU time of the multiresolution computations is at least one sixth in comparison with the uniform solver (for maximum grid level J = 10 in each direction). This is due to presence of the elliptic solver, thus an extension to this work will be development of a multigrid solver or Krylov subspace methods for the multiresolution algorithm, to accelerate further the elliptic part of the method. Another suggestion is to develop a wavelet transform based on fifth-order interpolation which will lead to higher compression rates in comparison to the cubic interpolation. This is the case especially for incompressible flows where the flow fields are much smoother than that of compressible flows where shocks may occur in the field. The volume penalization method has been introduced into the finite difference solvers, both on the uniform grid solver and the adaptive one, to be able to compute flows in complex geometries, like curved walls. Comparing the solutions obtained via the penalization method with that of classical no-slip and no-penetration boundary conditions showed that for the dipole-wall collision with a straight wall, for Reynolds 1000 with 1024 2 grid point, we found an average deviation of less than 2% and a maximum deviation of 6% in the total energy evolution. For evolution of the total enstrophy and the total palinstrophy the difference between the two methods of applying the boundary conditions is however larger, a maximum difference (coincide with the first collision) of 30% in total enstrophy and 50% in total palinstrophy evolutions is seen. The explanation is given by the first-order accuracy of the volume penalization method at the wall. Hence small permeability parameters are necessary and fine grids near boundaries have to be used. Furthermore it was shown that the finite difference computations with classical boundary conditions converge towards the pseudo-spectral computations of Clercx et al. For Reynolds 1000 it was found that the grid resolution of J = 11 in each direction yields to reasonable results for energy, however some difference, 12% in enstrophy and in particular 20% for palinstrophy are still present. It is anticipated that for grid resolutions of J = 12 in each direction these differences will become much smaller, but the computational time with out parallel processing is currently too prohibitive. Extension of the volume penalization method to a high order (e.g., second order) immersed boundary method (IBM) will be another forward step. One interesting perspective is the identification and the tracking of the vortex cores in the flow field to study more quantitatively the details of the vortex dynamics. Extending the program to include a dynamic and distributed memory option for parallel computations with message passing interface (MPI) will be another step of the future investigation by using either hash tables or tree data structures. For helpful discussions in this area see [START_REF] Brix | Parallelisation of multiscalebased grid adaptation using space-filling curves[END_REF] and [START_REF] Hejazialhosseini | High order finite volume methods on wavelet-adapted grids with local time-stepping on multicore architectures for the simulation of shock-bubble interactions[END_REF].

(1. 16 )Figure 1 . 2 :

 1612 Figure 1.2: The scaling function Φ 1,i of wavelet transform corresponding to linear interpolation.

Figure 1 . 3 :

 13 Figure 1.3: The scaling function Φ 1,i of wavelet transform corresponding to cubic interpolation.

4 XF 5 Figure 1 . 4 :

 4514 Figure 1.4: Sparse point representation of 1D functions, obtained by WT with cubic interpolation (J = 10), filtered with threshold ϵ = 1 × 10 -3 . The green dots (marked •) show the retained grid points.

Figure 1 . 5 :

 15 Figure 1.5: Schematic representation of two-dimensional wavelet transforms : rectangular WT (left), tensor product WT (right), picture taken from [41].

st = 2 ( 5 st = 4 ( 16 st = 6 (

 254166 linear interpolation) with β 1 = 0.cubic interpolation) with β 1 = 9 16 , β 2 = -1 fifth-order interpolation) with β 1 = 150 256 , β 2 = -25 256 , β 3 = 3 256

37 )

 37 illustrated in Fig. 1.7 (b) and a two-dimensional Gaussian illustrated in Fig. 1.8 (b). Their sparse point representation, with the cubic interpolating wavelet transform, for J = 9, filtered with ϵ = 1×10 -3 is also shown in Fig. 1.7 (a) and Fig. 1.8 (a). Results of error analysis and compression rates with the linear and cubic interpolating wavelet transforms, by applying different thresholding parameters are reported in Tables 1.3 and 1.4.

( a )Figure 1 . 6 :

 a16 Figure 1.6: Schematic representation of the two-dimensional wavelet transform with p = 4, for calculation of the wavelet coefficients d µ j (i x , i y ), µ = 1, 2, 3.In the onedimensional interpolations in x and y directions, (µ = 1, 2) points at the finer level j + 1 (marked •) must interpolated from the points at the coarser level j (marked •), illustrated in (a). In the two-dimensional interpolation, (µ = 3) point at the level j + 1 (marked •) must interpolated from the points at the same level j + 1 (marked •), illustrated in (b), picture taken from[START_REF] Vasilyev | An adaptive multilevel wavelet collocation method for elliptic problems[END_REF].

Figure 1 . 7 :Figure 1 . 8 :

 1718 Figure 1.7: Function (1.37) and its sparse point representation (marked •) with the linear interpolating wavelet transform, J = 9, ϵ = 10 -3 , Compression = 99.4%, L 1 -Error= 10 -15 .

8 6 ×

 6 10 -3 1 × 10 -3 4 × 10 -6 Cubic 1 × 10 -3 98.4 9 × 10 -4 1 × 10 -4 4 × 10 -7 Cubic 1 × 10 -5 58.1 4 × 10 -6 3 × 10 -7 1 × 10 -9 Cubic 1 × 10 -8 6.1 6 × 10 -11 4 × 10 -13 5 × 10 -15 Cubic 1 × 10 -13 0.0 6 × 10 -16 2 × 10 -17 9 × 10 -20

Figure 2 . 1 :

 21 Figure 2.1: Hung point (marked •) and active points (marked •), picture taken from [32].

Figure 2 . 2 :

 22 Figure 2.2: Schematic representation of the fourth-order Runge-Kutta method. In each time step the RHS operator must be evaluated four times: once at the initial point (marked •), twice at trial midpoints (marked •) and once at a trial endpoint (marked •). From these derivatives the value of the function in the next time step (marked •) can be calculated, picture taken from [18].

Figure 2 . 3 :

 23 Figure 2.3: Illustration of the adjacent zone for adaptive wavelet methods in wavelet coefficient space with level (scale) index j and position index i. Solid circles (marked •) indicate the positions of wavelets. The locations of wavelets with significant coefficients kept after thresholding at time t are marked by the orange bell curve (--).The second bell-shaped region (---) on the right of the original one marks the locations of the significant wavelet coefficients at the end of the time integration, i.e., at t + ∆t. The safety zone is represented by red bell-shaped region, picture taken from[START_REF] Schneider | Wavelet methods in computational fluid dynamics[END_REF].

4 .

 4 Add safety zone to the list of the new active points (a) Add neighbor points at the same and one above levels (b) Guarantee the gradedness of the new active points (optional) (c) Add necessary points to the current list of the active points, for having a consistent direct or inverse WT 5. Apply IWT to the new active points to compute the values of the independent variables (or interpolate the values of all newly added points via IWT with d = 0) 6. Perform the time evolution of the independent variable for all the active points (a) Search for the nearest active point to determine dist for all active points (b) Check for the existence of all other neighbors of the active points with distance dist, mark all the missing points as the hung points (c) Interpolate the values of the hung points via IWT with d = 0 (d) compute the spatial derivatives for the given PDE via FDM with symmetric stencils 7. Go to step 2, if T < T end N.B. 1Before interpolation of the values of an independent variables via IWT (from the coarsest level up to the finest level) in some grid points (with wavelet coefficients equal to zero, d = 0), it is necessary to mark all the intermediate necessary points for having a consistent WT, (from the finest level down to the coarsest level) and adding them to the list of the points to be interpolated.N.B. 2For multi-step methods in time integration such as Runge-Kutta family before calculation of spatial derivatives at intermediate steps, the value of u * for the hung points, must be interpolated again from the new values of active points. But (a) and (b) will be done once in each time step.

Figure 2 . 4 :

 24 Figure 2.4: Domain of the solution and the immersed body, Ω = Ω f ∪ Ω s .

5 Figure 3 . 2 :

 532 Figure 3.2: Initial condition (left) and time evolution (right) computed with the multiresolution solver, threshold ϵ = 10 -3 , ν = 10 -3 , ∆t = 10 -5 and a maximum grid level J = 11. The exact solution is represented by (-) solid color lines.

6 (f) t = 1. 5 Figure 3 . 3 : 5 Figure 3 . 4 : 2 (b) ϵ = 5 × 10 - 3 Figure 3 . 5 :

 65335342510335 Figure 3.3: The evolution of the adaptive solution during time (active points represented by •) computed with the multiresolution solver and the exact solution given by Eq. (3.3) for the Burgers equation, ν = 10 -3 , (represented by --), with maximum grid level J = 11 and a threshold parameter ϵ = 10 -3 .

ϵ = 5 × 10 - 3 Figure 3 . 6 :

 510336 Figure 3.6: The evolution of the number of active, important (corresponding to the retained points after filtering of wavelet coefficients), safety zone, hung and interpolated points for the wavelet transform during the computation with the multiresolution solver for the Burgers equation, with maximum grid level J = 10 (1025 grid points).

3 . 8 45 Figure 3 . 7 :(f) t = 1 Figure 3 . 8 :

 384537138 Figure 3.7: The evolution and collision of the vortices (represented by the colored isolines) with walls, maximum grid level J = 11 in each direction, threshold ϵ = 10 -3 , for Reynolds 1000. (a) The initial vorticity distribution; (b) Shedding of the surrounding vorticity shields; (c) Formation of a forward-traveling Lamb-type dipole; (d) Formation of a weaker backward-traveling dipole; (e) First collision with the right wall; (f) Formation of secondary vortices.

Figure 3 . 9 :

 39 Figure 3.9: Evolution of the adaptive grid during the computation of the dipole-wall collision with the multiresolution solver with threshold ϵ = 10 -3 , maximum grid level J = 11 in each direction for Reynolds 1000.

Figure 3 . 10 :- 4 Figure 3 . 11 :

 3104311 Figure 3.10: Evolution of the adaptive grid during the computation of the dipole-wall collision with the multiresolution solver with threshold ϵ = 10 -3 , maximum grid level J = 11 in each direction for Reynolds 1000.

- 4 .

 4 The total enstrophy and palinstrophy during time, with ϵ = 10 -3 , are compared with the computations of Clercx et al. in Fig.3.14, (a) and (b), respectively. The same quantities with ϵ = 10 -4 are compared with that of Clercx et al. in Fig.3.15, (a) and (b), respectively. The accordance of enstrophy and palinstrophy curves during time with Clercx data with J = 11 are better than that of J = 10, which is promising. Another point to be noted is that the results for ϵ = 10 -4 and ϵ = 10 -3 are almost identical and the multiresolution code works properly with ϵ = 10 -3 . Hence we will use ϵ = 10 -3 for the multiresolution computations. The number of grid points and the corresponding time steps and CPU times on an Intel-P4 CPU (1.3 GHz), are given in the

Figure 3 . 12 :Figure 3 . 13 :

 312313 Figure3.12: Comparison of the total energy E(t) between the data from Clercx and the present finite difference computation with a uniform multigrid solver for Reynolds 1000 and maximum grid level J = 10 in each direction.

Figure 3 . 17 :Figure 3 . 18 :

 317318 Figure3.17: Comparison of the total energy E(t) between the uniform grid solver and the multiresolution computation with thresholds, ϵ = 10 -3 and ϵ = 10 -4 for Reynolds 1000 and maximum grid level J = 9 in each direction for all simulations.

Figure 3 . 19 :Figure 3 . 20 :Figure 3 . 22 :

 319320322 Figure 3.19: Schematic representation of penalized cavity (blue represents the fluid domain) where in solid domain (represented by red) η = 10 -3 is imposed. The geometry and the initial vorticity distribution are the same as the benchmark.

(f) t = 1. 0 Figure 3 . 23 :

 0323 Figure 3.23: The evolution and collision of the dipole (vortices represented by colored isolines) with straight penalized walls for Reynolds 10000. Maximum grid level J = 10 in each direction and η = 10 -3 in the penalized margin is imposed. The black lines represent the boundaries of the penalized domain.

Figure 3 . 24 :

 324 Figure 3.24: The evolution and collision of the dipole (vortices represented by colored isolines) with straight penalized walls for Reynolds 10000. Maximum grid level J = 10 in each direction and η = 10 -3 in the penalized margin is imposed. The black lines represent the boundaries of the penalized domain.

Figure 3 . 25 :

 325 Figure 3.25: Mask function (blue represents the fluid domain) with η = 10 -3 inside the convex wall, with R = 1.5, x c = 3.1, y c = 1.0, the initial vortices are placed at distance equal to one from the right wall.

Figure 3 . 26 : 4 (Figure 3 . 27 :

 3264327 Figure 3.26: Mask function (blue represents the fluid domain) with η = 10 -3 inside the concave wall, with R = 1.5, x c = 0.5, y c = 1.0, the initial vortices are placed at distance equal to one from the right wall.

Figure 3 . 28 : 9 (Figure 3 . 29 : 8 ( 8 (

 328932988 Figure 3.28: The evolution and collision of the dipole (the vortices represented by colored isolines) with a convex wall for Reynolds 10000. The center of the mask is placed at (x c = 3.1, y c = 1.0) with R = 1.5, maximum grid level J = 10 in each direction is imposed.

Figure 3 . 30 :Figure 3 . 31 :Figure 3 . 32 :

 330331332 Figure 3.30: The evolution and collision of the dipole (the vortices represented by colored isolines) with a concave wall for Reynolds 10000. The center of the mask is placed at (x c = 0.5, y c = 1.0) with R = 1.5, maximum grid level J = 10 in each direction is imposed.

Table 1 .

 1

	1: Results of the error analysis and compression rates for the function (1.28)
	with the linear and cubic interpolating wavelet transforms, J = 10.
	WT	ϵ	% Compression L ∞ -Error L 1 -Error L 2 -Error
	Linear 1 × 10 -1	96.2	2 × 10 -2 1 × 10 -3 5 × -5
	Linear 1 × 10 -2	93.0	9 × 10 -3 2 × 10 -4 1 × -5
	Linear 1 × 10 -3	86.8	2 × 10 -4 6 × 10 -5 2 × -6
	Linear 1 × 10 -4	72.7	5 × 10 -5 1 × 10 -5 5 × -7
	Linear 1 × 10 -5	42.1	1 × 10 -5 2 × 10 -6 9 × -8
	Linear 1 × 10 -6	0.0	9 × 10 -16 3 × 10 -17 5 × 10 -18
	Cubic 1 × 10 -1	96.3	7 × 10 -2 6 × 10 -2 1 × -4
	Cubic 1 × 10 -2	95.3	9 × 10 -3 5 × 10 -5 1 × -5
	Cubic 1 × 10 -3	94.0	5 × 10 -5 1 × 10 -6 1 × -7
	Cubic 1 × 10 -5	92.4	8 × 10 -7 9 × 10 -8 4 × -9
	Cubic 1 × 10 -7	88.7	2 × 10 -8 5 × 10 -9 2 × 10 -10
	Cubic 1 × 10 -9	67.0	6 × 10 -11 2 × 10 -11 9 × 10 -13
	Cubic 1 × 10 -11	36.7	3 × 10 -12 9 × 10 -13 4 × 10 -14
	Cubic 1 × 10 -13	0.0	5 × 10 -15 9 × 10 -16 4 × 10 -17

Table 1 .

 1 2: Results of the error analysis and compression rates for the Gaussian function, δ 2 = 0.1, with the linear and cubic interpolating wavelet transforms, J = 10.

Table 1 .

 1

		4	10 -15	10 -15	10 -15
	Linear 1 × 10 -7	99.4	10 -15	10 -15	10 -15
	Linear 1 × 10 -13	99.4	10 -15	10 -15	10 -15
	Cubic 1 × 10 -1	98.2	10 -15	10 -15	10 -15
	Cubic 1 × 10 -7	98.2	10 -15	10 -15	10 -15
	Cubic 1 × 10 -13	98.2	10 -15	10 -15	10 -15

4: Results of the error analysis and compression rates for the two-dimensional Gaussian function with the linear and cubic interpolating wavelet transforms, J = 9.

  1 equation in space and will be solved numerically via iterative methods like successive over relaxation or multigrid methods. With the use of auxiliary relations for velocity components it is possible to eliminate velocity vector from Eq. (2.26) and obtain

	∂ω ∂t	+	∂ψ ∂y	∂ω ∂x	-	∂ψ ∂x	∂ω ∂y	= ν	( ∂ 2 ω ∂x 2 +	∂ 2 ω ∂y 2	)	+	( ∂F y ∂x	-	∂F x ∂y	)	(2.28)

Table 3

 3

	.1.

Table 3 .

 3 1: Number of grid points and the corresponding time steps ∆t and CPU time (second) for one RK4 iteration with multigrid and PSOR solver with ||Res|| ∞ ≤ 10 -7 , for Re = 1000 and maximum grid level J in each direction.

	Grid	J	∆t	MG PSOR
	65 × 65	6 1 × 10 -2 0.003 0.01
	129 × 129	7 5 × 10 -3 0.01	0.06
	257 × 257	8 1 × 10 -3 0.07	0.2
	513 × 513	9 5 × 10 -4	0.6	1.2
	1025 × 1025 10 1 × 10 -4	2.5	7
	2049 × 2049 11 5 × 10 -5	-	74
	4097 × 4097 12 1 × 10 -5

  Figure 3.16: Convergence study for the total enstrophy Z(t) and the total palinstrophy P (t) toward the data from Clercx with the present finite difference computation with uniform grid solver for Reynolds 1000 and maximum grid level J = 8, 9, 10, 11 in each direction.

		3500					5E+07
					Clercx et al.					Clercx et al.
					256 2					256 2
		3000			512 2 1024 2		4E+07			512 2 1024 2
					2048 2					2048 2
		2500
	Enstrophy	1500 2000				Palinstrophy	2E+07 3E+07
		1000					1E+07
		500
							0
		0	0.25	0.5	0.75	1	0	0.25	0.5	0.75	1
				Time				Time
			(a) Total enstrophy		(b) Total palinstrophy*
				2
				1.9 1.8			Uniform 512 MR Eps=10 -3 2 MR Eps=10 -4
				1.7
				1.6
			Energy	1.3 1.4 1.5
				1.2
				1.1
				1
				0.9
				0.8
				0	0.25	0.5	0.75	1
						Time

Le cose est développée en FORTRAN et accessible à tous sur demande par mail [55].

The code is developed in FORTRAN and is accessible for all [55].

Perturbations will spread in all directions with the speed of sound which approaches to infinity in the incompressible limit.

*Notice: There is an error in some curves of palinstrophy P (t) noted by *, because of not taking to account the vorticity derivatives at boundaries via forward/backward stencil, which are significant. Otherwise all the curves of palinstrophy P (t) obtained by present study will be above that of Clerx et al. as is in the case of enstrophy. I am apologizing the reader to not be able to modify the figures now but in a probable future publication this will be considered.

DT max / DT

large, e.g., equal to 0.1, to avoid the truncation error from falling in the range of round-off error for ∆t max /16. The other simulations are performed with ∆t max /8, ∆t max /4, ∆t max /2 and ∆t max . The computations start from an initial condition u(x, 0) = sin(x), x ∈ [0, 2π] at t = 0 and stop at t = 1, so the time step and the number of iterations for each simulation is different (boundary conditions are fixed to zero at x = 0, 2π). Different errors as a function of ∆t max /∆t are compared with the theoretical (-4) slope in Fig. 3.1 (a), a good agreement can be seen. Then an error analysis to examine the second-order discretization in space is performed. Again we consider the Burgers equation to show the accuracy of the solver on the uniform grid. A simulation with J = 12 grid points is considered as reference solution to compute the errors, other simulations were performed with 2 11 , 2 10 , 2 9 and 2 8 grid points. The computations start from an initial condition u(x, 0) = sin(x), x ∈ [0, 2π] at t = 0 and stop at t = 1 (with fixed boundary conditions). The time step is chosen small enough, i.e., 4 × 10 -5 , which is stable for time integration with 2 12 grid points, to minimize the error due to time integration, consequently the number of iterations for all simulations are the same. Different errors are compared with the theoretical (-2) slope in Fig. 3.1 (b). Then a one-dimensional adaptive solver based on the multiresolution algorithm presented in Chapter 2 and Eq. (3.2) for the viscous Burgers equation (3.1) is developed. Then we perform multiresolution computations for the viscous Burgers equation (3.1), for which analytical solutions are known. For (x, t) ∈ R × [0, +∞), an analytical solution in an infinite domain is given in the form of wall with radius R = 1.5 will be examined where its center is placed at x c = 0.5, y c = 1.0. The corresponding mask functions and initial conditions are illustrated in Fig. 3.25 and Fig. 3.26, respectively. For having a qualitative idea about the collision of the dipole with a convex wall, the flow evolution from the initial condition t = 0 up to t = 2, is illustrated in Fig. 3.27 and Fig. 3.28 for Reynolds 10000. This will be useful for a qualitative comparison with dipole collision with straight wall. The computations were done with the same solver over the same grid points, i.e., 1025 2 . Next the evolution and the collision of the dipole with a concave wall from the initial condition t = 0 to t = 2, is illustrated in Fig. 3.29 and Fig. 3.30 for Reynolds 10000. Comparisons of total energy, enstrophy and palinstrophy between collision with a straight wall, a convex and a concave wall, with R = 1.5, for Reynolds 10000 and maximum grid level J = 10 in each direction, are shown in Fig. 3.31 and Fig. 3.32 (a) and (b), respectively. In the total energy evolution, the curve corresponding to the straight wall lies between those of the convex and concave walls up to t = 1.1, then it falls below the those of curved walls. But the total energy evolution for the concave wall lies always above the others, which can be explained as follows; For total enstrophy and palinstrophy, the curve for the concave wall usually lies below the others, which can be justified by considering the fact that in that case the surface is minimum, so enstrophy production at walls is reduced with respect to the other cases. Thus by considering Eq. (3.16) the results are qualitatively satisfactory. Some heuristic general features of the flow for different configurations can be summarized as follows; The time of the first collision with the right wall at sufficient high Reynolds numbers, independent of the curvature of the hitting wall is more or less the same, approximately equal to t col = 0.3. There will be a rise in the vorticity of small newly generated vortices after Chapter 4

Conclusions and perspectives

In the current master thesis a self-adaptive code for direct numerical simulation of twodimensional unsteady incompressible flows has been developed. The new algorithm is based on a central second-order finite difference discretization of the governing equations coupled with Harten's point value multiresolution analysis. The multiresolution analysis uses cubic interpolating wavelet transform for grid adaptation. This kind of wavelet transform will impose the necessity of using Cartesian uniform grids. The use of these types of grids will lead to good accuracy, efficiency and simplicity of the solver. Alongside with these advantages, the implementation of curved boundaries in Cartesian grids is a dilemma. The remedy for this problem is to use the idea of immersed boundary methods, e.g., volume penalization method. The multiresolution analysis allows to reduce the number of active grid points significantly, nevertheless the analytical accuracy order of the underlying numerical scheme is preserved. For time integration explicit Runge-Kutta methods of different order are implemented. First the algorithm was applied to the one-dimensional Burgers equation and the errors in time and in space have been assessed in detail. Then the method was extended to the incompressible two-dimensional Navier-Stokes equations expressed in vorticity-stream function formulation. Therewith the near wall behavior of two-dimensional vortical flows has been studied for different geometries, i.e., a dipole-wall collision with either straight or curved walls. For implementation of curved boundaries in a Cartesian grid the volume penalization method has been used. This method treats solid obstacles or walls as porous media, no-slip and no-penetration boundary conditions can thus be enforced by using permeability coefficient of the solid region which tend to zero. The obtained results show that the CPU-time of the adaptive simulations can be significantly reduced with respect to simulations on a regular uniform grid. The only parameter to be adjusted in the adaptive solver is the threshold parameter ϵ, which directly controls the compression rate, hence the introduced error in the nonlinear filtering procedure. The numerical experiments showed that for the considered one and two dimensional problems with maximum grid level less than J = 12 the proper threshold parameter is ϵ = 10 -3 . The average number of active points over the base grid is less than 10% for one-dimensional problems, thus the compression rate of grid points is usually more than 90%, which is promising. The time of the multiresolution computations for one-dimensional problems is reduced by more than a factor two