Julie Stal-Le Cardinal

Franck Marle
email: franck.marle@ecp.fr

Project: the Just Necessary Structure to Reach your Goals

Keywords: managing projects, processes, design, competence assignment, dysfunctions

This paper proposes a definition process of the project structure which should be constructed in order to reach the objectives and to deliver the final results.

The first part of the paper specifically describes this process, with its inputs, tools and methods and outputs. The inputs of this process, initial situation, objectives and environment, are described in part 2, with research proposals on environment management. Other research proposals about tools and methods are presented in part 3, focused on the topics of scope and activity definition and resource assignment. Concrete propositions stemming from research works and from their application at VALLOUREC and PSA PEUGEOT CITROEN will be presented all along the paper.

Introduction

A project is a temporary endeavor undertaken to reach some objectives and to deliver some results. The project is then a change vector in companies, markets and society. It consists of a start, an initial situation, and of a finish, with results that change the company's situation, in terms of internal performance, product offer, communication tools, and so on.

Project management consists of the whole concepts, methods and tools, in order to bring the project from start to finish, which means in order to reach the objectives and to deliver the results of the project.

The project lifecycle can be described in several ways, for instance with the PMI® processes: initiating, planning, executing and controlling, closing. The planning process consists of identifying, estimating and preparing the whole activities that have to be executed in order to reach the objectives and to deliver the results. Its output is the project plan, which means the project structure, deliverables, activities and resources, estimated on time, cost and quality dimensions.

According to a survey executed by the Standish Group International in 2000 [26], 80% of project successes or failures may be linked to planning, like bad scope definition, bad stakeholder analysis, bad activity decomposition and bad resource assignment.

The planning process is then very important and not under control. Very important, because the decisions made upstream in planning phase have much bigger consequences in the downstream execution phase: failing to plan is planning to fail. Not under control because it is impossible to predict the future, we can just make forecasts and estimations. But as the available information is not sure, and as the context is very complex, uncertain and changing, the plan is always not reliable.

The scope of the paper is about the process from the idea, the objectives to the project structure (planning), and not from the structure to the results (the execution). We call it the project structure definition process.

Project Structure Definition

A project is a transformation process, from an initial to an expected final situation, evolving in an often complex and changing environment.

The initial situation is composed of:

-The company's corporate strategy, -Historic and standards: global or specific to the company, -Initial resources available in the company: human resources, skills, knowledge, material resources, money,

-Constraints and assumptions

The description of the final situation gives the project objectives. Project objectives may result from or be consistent with the corporate strategy, and may be a response to a constraint (legal, social, ecological) or an opportunity / risk (technological or financial for instance).

They can be expressed with customer requirements analysis like functional analysis for product development. Therefore, to develop a new product, for instance, the beginning of the project consists in defining the expected functions, the expected delivery time and the budget.

In figure 1, the process is displayed with a loop, in order to show the existence of the project objectives at the beginning, and their comparison with project results after execution.

The project environment is composed with the whole stakeholders involved in the project lifecycle.

The transformation from initial to final situation can be made by planning and executing some project activities, using project resources and organizing in a logical way. It may include product-oriented activities, human management activities, and quality or procurement Everything is uncertain and may change during the project lifecycle: objectives, environment, resources, constraints and assumptions, and so should the structure change. It is not in the scope of this paper to see when and on which conditions the structure must change.

Inputs of the process

We consider here any project that takes place in an industrial environment. The definition of the project structure has to take into account parameters presented in the three following parts: the initial situation, the objectives and the environment. In each part, we are going to describe each of these parameters, and analyze their dysfunctions and their consequences. For the last one, the project environment, a proposal is made from a research work to help to manage environment during the project.

Initial situation: to know where we are today

Description

Information, here, depends on the situation of the company, and depends on the initial need of the project. The project should be aligned with the corporate strategy, so the project manager must be able to show how his project contributes to the global strategy achievement of the performing organization. The project manager should take into account additional information, like historic and standards, initial available resources, constraints and assumptions.

Dysfunctions and consequences

The initial need is seldom clearly formulated. It is the responsibility of the project manager, or the management or sponsor, if the project manager is not yet assigned, to clarify and refine this need in order to describe the project stakes, the justification of the project, and elaborate the first detailed specifications.

As the project manager is part of the description of the initial situation, an assignment dysfunction may involve dysfunctions for the project. For example, a project manager assigned too late or unofficially, or assigned to a project he does not have the skills for, may involves delays, under-quality and non satisfaction of the project objectives and of the customer.

Moreover, if the company does not have any high level indicator to know (in terms of performance, competencies, procedures and quality) where it is, it will be difficult then to know if means are well-planned with regard to the objectives, to the initial situation and to the environment.

Proposals

There is no proposal to better control the definition of any initial situation in this paper. The assumption is that the company has made the decision to launch a project, taking into account its own situation and the need. The company is then ready to give resources and expects results of the project.

Objectives: to know where we want to go

Description

Objectives are defined here as the expected end of the project, the goals that need to be reached. To define objectives involves delimiting the scope of the project, including its frontiers. What is not in the scope is not in the project, and reciprocally.

Objectives may be clearly defined with a project requirements definition, or a project specification. It is a formal document where may be detailed: the context and stakes of the project, the goals, the major deliverables, the customer(s), the frontiers, the global cost and deadline and the organization type (project manager, matrix or functional or project team, sponsor, …). It may be more or less detailed, and it may include functional or technical specifications. The most important thing is that it should correspond to the customer needs and expectations.

Dysfunctions and consequences

Dysfunctions and consequences are relative to the definition of perimeter, the existence of specifications of the project and the determination of available resources. Approaching the limits of the project, knowing what we do and what we do not do in the project is often an important moment of negotiation. It allows to enrich or to decrease the expected objectives.

The earlier on a project we approach this clarification of the objectives, the less danger we have of skidding and not achieving them [START_REF] Goh | Benchmarking the Learning Capability of Organizations[END_REF].

A dysfunction can be characterized by a difference between the expected action and the realized action and thus a non-achievement of the objectives [START_REF] Cardinal | A Study of Dysfunctions within Decision-Making Process[END_REF]. The objectives can be not SMART (Specific, Measurable, Achievable, Realistic, Time-related), as proposed by Drucker [START_REF] Drucker | The Practice of Management[END_REF], which reveals a dysfunction.

According to the Standish Group survey [26], the influence of objectives definition is important on project success or failure. People have declared this as a cause for 21% of the successes and 32% of the failures, decomposed into complete (or not), realistic (or not) and stable (or not) specifications.

The main dysfunction here is a difference between the project specifications and the customer needs and expectations, for instance if one forgets one need, or puts a specification not needed by the customer.

Proposals

The list of dysfunctions and their consequences is in itself one proposal, but no other solution is introduced. Proposals about objectives definition and control are not in the scope of the paper.

Environment: to know what is around us

According to a Standish Group survey [26], the influence of environment is important on project success or failure. People have declared environment as a cause for 30% of the success and 28% of the failures. This part details some facts about complexity of interactions inside a project and outside the project, with its environment. Some proposals are made about interactions description, using a project interactions model [START_REF] Marle | Information Model and Methods in order to Help Decision-Making in Project Management[END_REF].

Description

In a project, there are a lot of exchange flows, of information, decisions and material or financial items, which make the project complex, as described by Baccarini [4] or Williams [START_REF] Williams | The Need for New Paradigms for Complex Projects[END_REF]. For example, with N persons, there are N*(N-1)/2 possibilities of communication flows. These persons, or organizational structures, have different interests and different influence capacities, positive or negative, and these interests and influences may change during the project. The influence between two persons may be in only one way, or reciprocal. The action of identifying the environment of a project, in order to anticipate what influences, positive or negative, may affect the project or may be affected by the project, is called stakeholder management.

Dysfunctions and consequences

If only one stakeholder is not taken into account, forgotten or neglected, the risks of the project are increased. A project can be initiated, modified or stopped by influence of only one stakeholder.

For some stakeholders, below is an example of possible dysfunction and its potential consequence:

-The customer/user: if not taken into account, the customer may not validate the project, or may not pay, or may not want to work with us any more.

-Executive management, project sponsors, program director, steering committee and executive committee: they have the power to modify or to stop the project, and they may do it if they do not have proper and periodic reporting about the project and its stakes.

-Other projects: a project shares priorities and resources with other projects. Not to be aware of the other projects, and not anticipating our own resource needs in regard to other projects needs may conduct to a global constraint impossible to solve.

-Other items, like departments, suppliers, sub contractors, partners, …

Proposals

The project interactions model

The most known tool is the stakeholder analysis. We have elaborated, in collaboration with PSA Peugeot Citroën, a model describing the interactions existing inside a project, and between the project and its environment. This model is detailed in [START_REF] Marle | Information Model and Methods in order to Help Decision-Making in Project Management[END_REF], and is standard wherever we are in the project hierarchy. It has been built in order to complete stakeholder analysis, and to be applicable at each level, sub-project or work package, and not only at the project level. It is composed of seven types of interactions, and of some objects, like project, objective, deliverable, activity, process and actor:

-The hierarchical link: exists in WBS and PBS for example.

-The sequential link: exists in every project schedule, and shows time dependency.

-The contribution link: a result of object O1 (a project) may contribute to the result of object O2 (an objective, or another project).

-The influence link: O1 (sponsor) may modify O2 (project), or the result of O1 (project) may impact O2 (users).

-The resource link: two objects share the same critical resource.

-The proximity link: O1 looks like O2, which may give an opportunity to reuse experience and best practices.

-The exchange link: is only an exchange of information and data.

The integration of external constraints

In concurrent engineering, the goal is to integrate downstream constraints in the upstream activities, in order to avoid rework and waste. The goal is exactly the same in identifying the six other types of interactions: we shall integrate the whole constraints of the environment before to endure them, because the consequence of a mistake in planning is multiplied by ten or hundred in execution.

A last key point of this part is that environment influences the project structure. If one stakeholder is forgotten, the structure will be different. If the importance given to one stakeholder is different, the structure will be different. The proposal is then to integrate environment during the project structure definition process, in order to perform right stakeholder management during the rest of the project. Our model helps to identify, characterize and visualize interactions, as detailed in [START_REF] Marle | Information Model and Methods in order to Help Decision-Making in Project Management[END_REF], and a software prototype was performed and tested in PSA Peugeot Citroën, on some projects of their progress plan. It was a 250-projects plan consisting of the whole actions of internal improvement, transformation, what we can call support. The test showed the complexity of this project network, and demonstrated the difficulty to manage this properly. It had an impact on project structure definition, but also on multi-project management. Now, there are much fewer and much bigger internal projects in this company to limit complexity in the management of such a project network.

Tools and Methods

We have first defined the project structure concept, then the inputs of the process that leads to this structure. We are now going to propose tools and methods to help building the structure of a project. The definition process of any project structure is a decision process made of two main sub-processes: the decomposition of the project and the resource assignment. In this chapter, each sub-process is described, then examples of dysfunctions are analyzed and our propositions of improvement are presented and justified. A last point not developed in this paper is about risk management, which consists of identifying, assessing and responding to the whole potential events, positive or negative, that may affect the project, because the project structure will be modified, enriched and refined.

Decomposition process

Description

Decomposition is a cognitive and complex problem-solving process. A non-human tool can not make decomposition, except by repeating previous projects or by using The input is often called « father » and the outputs « sons ». The outputs are not unique; it may be several decomposition possibilities. So, decomposition is also a decision-making process, a choice between decomposition alternatives.

Dysfunctions and consequences

-Frontiers definition problems between two sub-objects: authority balance, responsibility sharing, time waste due to bad scope limits, rework or work made two times … -Project decomposition may be not consistent with the existing organization decomposition: waste of time due to organizational and personal conflicts, -The execution of the whole sub-objects may not give the object: under-quality, delays, scope variance, non satisfaction of objectives and customer, -The reticence to use formalized and rigorous methods for a day-to-day action: but decomposing a project is far more complex than decomposing the tasks for planning a holiday trip.

Consequences of a bad decomposition on project performance are big, because the rest of the planning process is based on it, like cost, resource and time estimating, quality, communications and procurement activities. So, it involves cautiousness and use of risk mitigation actions, like the methods proposed below.

Proposals

Proposals are on potential decomposition generation and evaluation. The project manager's role is then to apply the decomposition decision, by communicating to concerned people, and finally to keep experience and lessons learned for next projects and next decompositions.

Generation step

Based on other approaches, like [START_REF] Pahl | Engineering Design[END_REF], [START_REF] Thomson | Requirements Engineering -Laying the Foundations for Successful Design[END_REF] and [START_REF] Zeng | An Environmental Decomposition-based Approach to Design Concept Generation[END_REF], we have defined a four-step solution generation process for project decomposition, consisting of:

-The identification of specific parameters of the object that describe its context. It allows to know what can be reused from past for this project, and what will be reusable from this project to future projects, -The information research into historic and standards, -The decomposition itself of the object, whether by using one or more decomposition criterion, described below, or by using creativity methods, -The identification of additional interactions with other objects, which are directly involved by the decomposition process, like described below.

Decomposition criteria

The proposal for decomposition criterions is expressed into a list of eight elements:

-by phase, -by profession, -by department, -by technical domain, -by objective or function, -by mechanical or structural sub-system, -by geographical location, -by addressee, in terms of population or process.

Other criteria lists are available [START_REF] Hubka | Functions Revisited. 13 th International Conference on Engineering Design[END_REF].

There are two types of decomposition: mono-and multi-criteria. Every project decomposition is mono or multi-criteria, with the criteria of the list described above.

For example, a new bicycle development project may be decomposed like this: project management, handlebar, fork, transmission, wheels, safety equipment, mechanics, prototyping, manufacturing, distribution, marketing/sales. This decomposition is a mix of organic (handlebar, fork, transmission, and wheels), technical (project management, safety, mechanics) and sequential decomposition (prototyping, manufacturing, distribution, marketing/sales). It involves more risk, on the frontiers definition. The other possibility would be to decompose first by organs, then by technical domain, and finally by phase, but the project structure would have three levels instead of only one, which is heavier. The best decomposition is the best compromise between scope, time and organizational parameters, which means the less risky.

Additional interactions

Decomposition involves other interactions than the only hierarchical link between father and sons. These other links are between the sons or between one son and another object elsewhere in the project. But there is no reciprocity: if there is a link between two objects, it does not imply that it is an additional link created by decomposition.

Mono-criteria and homogenous decomposition creates between the sons a link corresponding to these criteria: for instance, decomposing a project among its phases involves a sequential link between some of these phases. Criteria "by phase" => additional "sequential" link.

A multi-criteria decomposition creates multiples links corresponding to each criterion. A multi-criteria heterogeneous decomposition creates automatically between the sons several types of links, corresponding to each criterion. The decomposition of a project into electronics, mechanics, data processing and project management creates contribution links (from project management to the other technical areas), sequential links (from electronics to data processing) and exchange links (between mechanics and electronics).

These links should not be forgotten, as they are created in an invisible and indirect way.

Evaluation step

Every decomposition choice should be done only from an « a priori » evaluation, before the project is executed:

-Relative evaluation by comparing several alternatives, -Absolute evaluation if only one decomposition is proposed.

The aim of this evaluation is to assess if the decomposed structure has a chance or not to reach the project objectives. In order to realize this assessment, a list of questions has been developed, split in 7 main questions. The complete list is detailed in [START_REF] Marle | Information Model and Methods in order to Help Decision-Making in Project Management[END_REF]:

-Is the decomposition complete? -Is the decomposition homogeneous? -Does the number of sons lie between three and seven? -Is the decomposition innovating or was it already tested in preceding projects? -Is the decomposition consistent among the complete project? -What are the additional links generated with the hierarchical link? -What is the flexibility of the decomposition? Its robustness?

The problem with standard decomposition is that it does not take into account additional interactions, it is not evaluated, and it does not allow the team members to be creative in the construction of their own decomposition. It should be used only as a basis for inspiring decomposition, not as a rigid tool [START_REF] Mulenburg | Report of Research Examining the Characteristics of Managers of Complex Contemporary Projects in the NASA[END_REF]. Standardization should not be at the scale of the company but more at a project scale. Payne [START_REF] Payne | Company-wide Project Management: the Planning and Control of Programmes of projects of different type[END_REF] proposes, for instance, the planning and control of programs of projects of different type.

Synthesis: tips for a right decomposition

This is developed in [START_REF] Marle | Plan more by Planning less: Contributions to a Recursive Planning Process[END_REF] and [START_REF] Marle | Information Model and Methods in order to Help Decision-Making in Project Management[END_REF], and was implemented in a simplified manner in PSA Peugeot Citroën as a procedure included in project planning:

-Try to do homogeneous decompositions, -Always evaluate the chance of the decomposed structure to succeed. By comparison if there are several alternatives, but it is not necessary, -Identify the additional links generated by the decomposition, -Take into account the parameters which may influence the decomposition result of the decomposition process, like the context of the project, the pressure on the project manager, the liberty and autonomy, the innovation degree, and the existing standards.

Assignment process

Our work and proposals about the assignment process is based upon a study we have made within the Vallourec Group. During this collaboration, we have worked with two units and more than 20 projects during 2 years. The purpose was to help in project management and more specifically to help in the choice of people as project manager or member of a project team. This is a key stake, and Turner [START_REF] Turner | On the Nature of the Project as a Temporary Organization[END_REF] describes a "project as a temporary organization, (…) as an agency for assigning resources to the management of change within the functional organization".

Description

In a company, projects are the translation of the strategic axes into actions. At the strategic level of a company, decision-makers have to choose people/actor to be project managers and then, a project manager has to build up his project team, he can choose actors from the company or external people. An example is the choice of contractors during the procurement phase [START_REF]Project Management Body Of Knowledge[END_REF]. We define an actor as a human being among company means. Material resources, software, hardware, people are part of the means in a company. What distinguishes an actor from the other means is that an actor has competencies quantified by a level, can make decisions and is able to characterize the impact of his action in advance.

So as to take part in a project, people/actors have to be chosen as responsible for some actions. The success of the project depends among others on the choice of the actors in the project [START_REF] Belout | Factors Influencing Project Success: the Impact of Human Resource Management[END_REF]. The "choice of actor" is a decision made by a decision-maker that consists in selecting, evaluating and choosing a person to accomplish an action. Mezher [START_REF] Mezher | Embedding Critics in Decision-Making Environments to Reduce Human Errors[END_REF] considers the decision as "a process which generates and evaluates alternatives and which makes choices among them."

Dysfunctions & consequences

The Standish Group survey already introduced before gives some statistics about the importance of human resources in projects: sufficient and competent resources represent more than 10% of successes and failures of the studied projects [26]. We propose here an analogy with the maintenance area [START_REF] Artikis | A Stochastic Model for Proactive Risk Management Decisions[END_REF] where "a failure is the stochastic cessation of an entity aptitude to accomplish a required function." According to Villemeur [START_REF] Villemeur | Sûreté de fonctionnement des systèmes industriels, fiabilité, facteurs humains, informatisation[END_REF], "after a failure identification, the entity is considered out of order. A breakdown is always due to a failure." By analogy, we consider an actor in a company as an entity whose function is to decide. [START_REF] Reginato | Project Management as a Core Competency[END_REF] consider a project actor as someone who can make decisions including the whole parameters of the project, managerial and technical. If the decision-maker is not able to make his decision, then there is a dysfunction, because he is not able to realize his function.

Reginato and Ibbs

We, therefore, propose the following definition of a dysfunction in a decision process: a dysfunction is a stochastic cessation of an actor's aptitude to make a decision. If the action expected to be accomplished has not been realized, then there is a dysfunction. The gap between the effective result and the negotiated objective to be reached is called the dysfunction value.

Dysfunctions are decomposed in elementary dysfunctions and represented along the DTL (Decision Time Line, [START_REF] Cardinal | A Study of Dysfunctions within Decision-Making Process[END_REF]). Each elementary dysfunction is related to a DTL step, independent from each others and cannot be studied as a combination of dysfunctions. The gap due to a dysfunction is then evaluated. Using our model, it is possible to give a temporal characterization of dysfunctions as well as a functional characterization that allows classification of various dysfunction types. Concerning actors' choice, we make the hypothesis that the reasons of the dysfunctions result from problems of competence.

Competences concern the decision-maker or the chosen actor. We suggest classifying these types of dysfunctions on a model of competence. Other models exist, like [START_REF] Crawford | Profiling the Competent Project Manager[END_REF], [START_REF] Gareis | Competences in the Project-Oriented Organization[END_REF] and [START_REF] Nakatani | Coordinated Interfaces for Real-Time Decision Making in Hierarchical Structures[END_REF], or selection criteria definition, like [START_REF] Midler | L'auto qui n'existait pas, management des projets et transformation de l'entreprise[END_REF]:

-Knowledge: know-who, know-why, know-what, -Know-how: technology, techniques, knacks, -Attitude: behavior, will, identity.

Proposals

The target process

The target process is a process to follow for a specific type of decision, choice of actor, and which helps to avoid dysfunctions, as illustrated in figure 2. But even if the different steps are properly followed, risks of dysfunctions still exist, nothing can completely erase them.

How to avoid dysfunctions? The decision card

The previous target process becomes a reality in an index card performed by the project managers which serves for discovering the risks a priori (context and quantification of tasks to perform), for estimating the decisions a posteriori (steps 2 to 6 of the target process) and for preserving the information with the aim of statistics, for example (result analysis, recommendations and capitalization).

Once such a card is completed, project managers, in the Vallourec Group, are able to analyze the reasons of the dysfunctions. With a consequent number of cards, it is possible then to study the general trend for a given company to cope with dysfunctions.

The capitalization database provides a means for dysfunction prediction, by capturing the firm's experience from prior decision-making tasks. Software accessing the capitalization database and the decision-making model can give an early warning to a decision-maker when an action or decision-making strategy is likely to lead to a known dysfunction and consequently a problem in quality, cost, or delay.

Conclusion

The aim of this paper is to help in the definition of a project structure that will correspond to the expected deliverables for a given project.

The main message is that building a correct project structure is achievable and gives more guarantee for success. On contrary, a bad structure is a guarantee for failure: failing to plan is planning to fail. We, then, show the feasibility of formalizing and managing the project structure definition process, by using innovating concepts and tools, in the area of interactions management, decomposition process and resource assignment process. This work is relevant for two reasons:

-Some recent statistics show the importance of projects and project management in companies and societies, and the importance of the planning phase on the rest of the project.

-The implication of two big companies as VALLOUREC and PSA PEUGEOT CITROEN in this research work, with their concern about the efficiency of their actions, shows that this topic is in their current preoccupations.

Practically, the concepts presented here have been transformed into industrial tools, like procedures and software modules. The reader could then gain advantage by using concrete templates, communication modules and decision-making tools.

The ideas of future developments are about -Evaluation of the impact of the structure on the project success (evaluation "a priori"), -Periodic re-estimation of the probability of success with the current structure (continuous evaluation), -A procedure for correcting a structure with a low success probability and implementation of the corrections in an ongoing project (corrective actions).

3 "

 3 activities. This transformation is what is called project structure in this paper. It is possible to use tools and methods to make this structure definition. Examples of tools and methods are Work Breakdown Structure, Gantt charts, resource assignment matrix and project baselines (performance, cost). The focus of the paper is about project structure definition. The figure 1 gives illustration of the project structure definition process, and its location in the global project lifecycle. The focus of the paper is only to study the local structure definition process. It may be implemented in any of the existing standards about project management and project lifecycle, like PMI®, IPMA, six-sigma for example.The initial situation, objectives and environment are described as the inputs of the project structure definition process, in part 2 "Inputs of the process", with a research focus on environment management. The research proposals about tools and methods are in part Tools and methods", focused on the topics of scope and activity definition, and resource assignment.

  standard templates. It is the operation of expressing an object into several smaller ones. Some examples of decomposition are the Work Breakdown Structure or Product Breakdown Structure. The major deliverables of the project or the main components of the product are decomposed into smaller, more manageable items. It allows better managing and controlling of each lowest-level element, called work package for WBS. It allows assigning each work package to one organizational and financial unit.

  The capitalization database can also be used to conduct a statistical study of dysfunctions within the firm: what are the risky steps in the decision-making process, what are the main sources of dysfunctions, what are the main effects… As far as Vallourec is concerned, statistics have been made to detect the risky steps in the decision-making process concerning a particular type of project. The risky steps are indeed different from a product development project to a human resources project, for instance.

Figure 1 .

 1 Figure 1. The location of project structure definition in the global project lifecycle

Figure 2 .

 2 Figure 2. The target process as a core process for a decision-maker

Corresponding author. Tel: 33-1-41131568. e-mail address: franck.marle@ecp.fr


Souleymane Kadri Harouna
email: souleymane.kadri_harouna@univ-lr.fr

Val Érie Perrier
email: valerie.perrier@imag.fr

DIVERGENCE-FREE WAVELET PROJECTION METHOD FOR INCOMPRESSIBLE VISCOUS FLOW ON THE SQUARE

Keywords: Divergence-free wavelets, Navier-Stokes simulation, projection method, Dirichlet boundary condition

. Here the temporal discretization is borrowed from the projection method. The projection operator is defined through a discrete Helmholtz-Hodge decomposition using divergence-free wavelet bases: this prevents the use of a Poisson solver as in usual methods, whereas improving the accuracy of the boundary condition. The stability and precision order of the new method are stated in the linear case of Stokes equations, confirmed by numerical experiments. Finally the effectiveness, stability and accuracy of the method are validated by simulations conducted on the benchmark problem of lid-driven cavity flow at Reynolds number Re = 1000 and Re = 10000.

1. Introduction. The characterization of turbulent flows is a continuing challenge encountered in several scientific areas. Physically, turbulent flows are characterized by the presence of many phenomena at different scales in interaction and with rapid variations in time and space. The mathematical equations that model turbulent flows are the incompressible Navier-Stokes equations, which are derived from newtonian laws in the context of hydrodynamics [START_REF] Temam | Navier Stokes Equations[END_REF]:

   ∂v ∂t -ν∆v + (v • ∇)v + ∇p = 0, ∇ • v = 0, (1.1)
on Ω ⊂ R d , an open domain with smooth or piecewise smooth boundary Γ = ∂Ω. Here v ∈ R d denotes the velocity vector field, p ∈ R is the pressure and ν > 0 is the kinematic viscosity. We focus in this paper on the two-dimensional equations (d = 2), the extension of our method to dimension three being straightforward.

To take into account the physics of the problem, we suppose that the fluid is confined in Ω, and does not cross the boundary Γ. In such case, the velocity field v must be tangential to the boundary:

v • n = 0 on Γ, (1.2)
where n denotes the outward normal to Γ. For solid boundary, the viscous flow satisfies a no-slip condition:

v = 0 on Γ, (1.3)
whereas a non homogeneous Dirichlet boundary condition: v = g on Γ, (1.4) holds for the lid driven cavity test case. More general boundary conditions may be handled, but will not be considered in this article.

The construction of performing numerical schemes is very important for effective models of prediction. The main difficulty in the numerical resolution of Navier-Stokes equations comes from the nature of equations which are nonlinear. Otherwise, physical boundary condition on v (1.2) and (1.3 or 1.4) can be simply imposed considering the velocity-pressure formulation. In such formulation, the projection method has for advantage to decouple the computation of the velocity v and this of the pressure p [START_REF] Chorin | Numerical simulation of the Navier-Stokes equation[END_REF][START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires II[END_REF].

In this article, we introduce a variant of the projection method based on the Helmholtz-Hodge decomposition [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF] of the main equation, explicitly computed in the wavelet domain. This will allow to avoid the principal drawbacks of the classical projection method.

We first introduce the divergence-free function space, with free-slip boundary condition, in order to deal with the incompressibility constraint:

H div (Ω) = {u ∈ (L 2 (Ω)) 2 : ∇ • u = 0, u • n| Γ = 0}.
(1.5)

By Stokes theorem, the space H div (Ω) is orthogonal to any gradient function in (L 2 (Ω)) 2 [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF]. Then, projecting the Navier-Stokes equations (1.1) onto H div (Ω) yields:

   ∂v ∂t + P[-ν∆v + (v • ∇)v] = 0, ∇ • v = 0, (1.6)
where P denotes the orthogonal Leray projector from (L 2 (Ω)) 2 to H div (Ω). According to the Helmholtz-Hodge decomposition of the term -ν∆v + (v • ∇)v, the pressure should satisfy the following equation:

∇p = -ν∆v + (v • ∇)v -P[-ν∆v + (v • ∇)v]. (1.7)
Now the difficulty relies in the integration in time of (1.6). The conventional projection method consists in a splitting of the operator ∂ ∂t -νP∆: a first step computes a -non divergence-free-velocity v * with the operator ∂ ∂t -ν∆, followed by a correction step, which consists in the projection of v * onto the divergence-free function space, using a Poisson solver with suitable boundary conditions.

We follow another point of view: in the simplest case of periodic boundary conditions, the first equation of (1.6) becomes:

∂v ∂t -ν∆v + P[(v • ∇)v] = 0, (1.8)
and the pressure p is recovered via:

∇p = (v • ∇)v -P[(v • ∇)v].
(1.9)

Remark that this formulation is very close to the projection method, since the numerical resolution of (1.8) reduces to a heat kernel integration, with as source term the projection of the nonlinear term P[(v • ∇)v]. This approach is common in spectral method [START_REF] Orszag | Numerical simulation of turbulence[END_REF], where P is explicit in Fourier domain. It was also used in [START_REF] Deriaz | Direct Numerical Simulation of Turbulence using divergence-free wavelets[END_REF][START_REF] Deriaz | Divergence-free and curl-free wavelets in 2D and 3D, application to turbulent flows[END_REF], where a divergence-free wavelet based resolution method was introduced. Using a finite difference scheme in time (for example backward Euler), the method is summarized as follows: starting with v n , compute v n+1 by

(1 -νδt∆)v n+1 = v n -δt P[(v n • ∇)v n ], (1.10)
where the term (v n • ∇)v n is explicitly computed on mesh grid points. Each time step requires the computation of the projection P[(v n • ∇)v n ], which is done through an iterative algorithm using alternatively the div-free and curl-free oblique wavelet projectors [START_REF] Deriaz | Orthogonal Helmholtz decomposition in arbitrary dimension using divergence-free and curl-free wavelets[END_REF]. This method gives rise to sparse divergence-free wavelet representation of the velocity, and adaptive discretizations can be derived. Furthermore, one can cite other wavelet-based methods in the velocity-vorticity formulation [START_REF] Farge | Wavelet transforms and their applications to turbulence[END_REF][START_REF] Muller | Adaptive Multiscale Schemes for Conservation Laws[END_REF][START_REF] Schneider | Wavelet methods in computational fluid dynamics[END_REF].

Our objective in the next coming sections is to provide an effective numerical method that extends the works of Deriaz-Perrier [START_REF] Deriaz | Direct Numerical Simulation of Turbulence using divergence-free wavelets[END_REF], more flexible for desired boundary conditions. First remark that in the case of physical boundary conditions (1.2) and (1.3 or 1.4), the situation becomes more complicated, since the projector P does not more commute with the Laplacian operator:

P(∆v) = ∆P(v).
Usual approaches consist in computing the divergence of (1.7), which leads to the well known Poisson equation for the pressure p:

∆p = ∇ • [(v • ∇)v], (1.11)
with Neumann boundary conditions.The resulting equations are then solved by standard methods for heat and Poisson equations, for which a large number of works exist [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF][START_REF] Temam | Navier Stokes Equations[END_REF].

On the other hand, we directly consider equations (1.6-1.7): the spatial approximation will handle divergence-free wavelet bases on square/cubic domains satisfying homogeneous boundary conditions, recently constructed (see [START_REF] Kadri-Harouna | Ondelettes pour la prise en compte de conditions aux limites en turbulence incompressible[END_REF][START_REF] Kadri-Harouna | Effective construction of divergence-free wavelets on the square[END_REF][START_REF] Stevenson | Divergence-free wavelet bases on the hypercube: Free-slip boundary conditions, and applications for solving the instationary Stokes equations[END_REF]). In such bases the Helmholtz-Hodge decomposition and the projector P become explicit and can be computed [START_REF] Kadri-Harouna | Helmholtz-Hodge Decomposition on [0, 1] d by Divergence-free and Curl-free Wavelets, Curves and Surfaces[END_REF]. We then deduce another formulation of the projection method for Navier-Stokes equations: contrarily to usual approaches [START_REF] Bell | A second-order projection method for the incompressible Navier-Stokes equations[END_REF][START_REF] Chorin | Numerical simulation of the Navier-Stokes equation[END_REF][START_REF] Kim | Application of a fractional-step method to incompressible Navier-Stokes equations[END_REF][START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires II[END_REF], the method we develop in next coming sections will not use a Poisson solver.

The paper is organized as follows. Section 2 recalls the setting of divergence-free wavelet bases on the square satisfying homogeneous boundary conditions, and the computation of Leray-Hopf projectors. Section 3 presents the wavelet-based projection method for the unsteady Stokes problem, with a convergence result, validated by numerical experiments. This method is extended for Navier-Stokes equations in Section 4, and compared to benchmark results on the well known lid driven cavity problem.

2. Divergence-free wavelets and Leray-Hopf projector. This section introduces the construction principles and main properties of divergence-free wavelet bases on the square [0, 1] 2 , and associated Leray-Hopf projector. Practical details of the construction recalled below can be found in [START_REF] Kadri-Harouna | Effective construction of divergence-free wavelets on the square[END_REF].

2.1. Divergence-free wavelets. Beginning with the seminal works of Lemarié-Rieusset and collaborators [START_REF] Jouini | Analyse multirésolution biorthogonale sur l'intervalle et applications[END_REF][START_REF] Lemarié-Rieusset | Analyses multi-résolutions non orthogonales, commutation entre projecteurs et dérivation et ondelettes vecteurs à divergence nulle[END_REF], several constructions of divergence-free (and curlfree) wavelets on the square have been provided [START_REF] Kadri-Harouna | Effective construction of divergence-free wavelets on the square[END_REF][START_REF] Stevenson | Divergence-free wavelet bases on the hypercube: Free-slip boundary conditions, and applications for solving the instationary Stokes equations[END_REF][START_REF] Stevenson | Divergence-free wavelets on the hypercube: General boundary conditions[END_REF][START_REF] Urban | Using divergence-free wavelets for the numerical solution of the Stokes problem[END_REF]. All are based on one-dimensional biorthogonal multiresolution analyses linked by differentiation and integration. The construction [START_REF] Kadri-Harouna | Effective construction of divergence-free wavelets on the square[END_REF] follows the two following main steps:

(i) Construction of two biorthogonal multiresolution analyses of L 2 (0, 1), denoted by (V 1 j , Ṽ 1 j ) and (V 0 j , Ṽ 0 j ) satisfying:

d dx V 1 j = V 0 j and Ṽ 0 j = { x 0 f (t)dt : f ∈ Ṽ 1 j } ∩ H 1 0 (0, 1). (2.1)
Each space is spanned by scaling function biorthogonal bases:

V 1 j = span{ϕ 1 j,k ; 0 ≤ k ≤ N j -1} and Ṽ 1 j = span{ φ1 j,k ; 0 ≤ k ≤ N j -1},
and

V 0 j = span{ϕ 0 j,k ; 0 ≤ k ≤ N j -2} and Ṽ 0 j = span{ φ0 j,k ; 0 ≤ k ≤ N j -2},
with dimension parameter N j = 2 j + c (for some small c ∈ N). For ε = 0, 1, the scaling functions ϕ ε j,k can be written as ϕ ε j,k = 2 j/2 ϕ ε (2 j x -k) inside the interval [0, 1], where ϕ ε is a compactly scaling function on R, but this is no more true near the boundaries 0 and 1 (idem for φε j,k ). In practice, the scale index j must be greater than some index j min , to avoid boundary effects [START_REF] Monasse | Orthogonal Wavelet Bases Adapted For Partial Differential Equations With Boundary Conditions[END_REF]. The biorthogonality between bases writes:

ϕ ε j,k / φε j,k ′ = δ k,k ′ .
Biorthogonal wavelet spaces (W 1 j , W 1 j ) are defined as usual in the biorthogonal wavelet formalism [START_REF] Cohen | Wavelets on the Interval and Fast Wavelet Transforms[END_REF]:

W 1 j = V 1 j+1 ∩ ( Ṽ 1 j ) ⊥ , W 1 j = Ṽ 1 j+1 ∩ (V 1 j ) ⊥ .
These spaces are generated by finite dimensional biorthogonal wavelet bases on the interval:

W 1 j = span{ψ 1 j,k ; 0 ≤ k ≤ 2 j -1} and W 1 j = span{ ψ1 j,k ; 0 ≤ k ≤ 2 j -1}. Biorthogonal wavelet bases {ψ 0 j,k } j≥jmin of W 0 j = V 0 j+1 ∩ ( Ṽ 0 j ) ⊥ and { ψ0 j,k } j≥jmin of W 0 j = Ṽ 0 j+1 ∩ (V 0 j )
⊥ are then simply defined by respectively differentiating and integrating the wavelets bases of (W 1 j , W 1 j ) j≥jmin [START_REF] Jouini | Analyse multirésolution biorthogonale sur l'intervalle et applications[END_REF][START_REF] Kadri-Harouna | Effective construction of divergence-free wavelets on the square[END_REF]:

ψ 0 j,k = 2 -j (ψ 1 j,k ) ′ and ψ0 j,k = -2 j x 0 ψ1 j,k . (2.2)
Homogeneous Dirichlet boundary conditions can be easily imposed on (V 1 j , Ṽ 1 j ) by removing scaling functions that reproduce constant functions at edges 0 and 1, prior biorthogonalization [START_REF] Masson | Biorthogonal spline wavelets on the interval for the resolution of boundary problems[END_REF][START_REF] Monasse | Orthogonal Wavelet Bases Adapted For Partial Differential Equations With Boundary Conditions[END_REF]. Then, the spaces

V d j = V 1 j ∩ H 1 0 (0, 1) = span{ϕ 1 j,k ; 1 ≤ k ≤ N j -2}, (2.3)
and

Ṽ d j = Ṽ 1 j ∩ H 1 0 (0, 1) = span{ φ1 j,k ; 1 ≤ k ≤ N j -2}, (2.4)
provide a biorthogonal multiresolution analysis for H 1 0 (0, 1), see [START_REF] Masson | Biorthogonal spline wavelets on the interval for the resolution of boundary problems[END_REF].

For the construction of divergence-free wavelets in (H 1 0 (Ω)) 2 , we will also need to impose homogeneous Neumann boundary conditions to both spaces V 1 j and Ṽ 1 j :

V dd j = {f ∈ V 1 j ; f (0) = f (1) = 0 and f ′ (0) = f ′ (1) = 0}, Ṽ dd j = {f ∈ Ṽ 1 j ; f (0) = f (1) = 0 and f ′ (0) = f ′ (1) = 0}.
Following [START_REF] Masson | Biorthogonal spline wavelets on the interval for the resolution of boundary problems[END_REF][START_REF] Monasse | Orthogonal Wavelet Bases Adapted For Partial Differential Equations With Boundary Conditions[END_REF], such space can be easily constructed by removing two scaling functions at boundaries, that allows to reconstruct the constant and polynomial of degree 1, then:

V dd j = span{ϕ 1 j,k ; 2 ≤ k ≤ N j -3} and Ṽ dd j = span{ φ1 j,k ; 2 ≤ k ≤ N j -3}. (2.5)
It can be proven that:

d dx V dd j ⊂ V 0d j , where V 0d j = H 1 0 ∩ V 0 j ,
Ṽ 0d j ⊂ Ṽ dd j .
The wavelet bases associated to the multiresolution analyses (V d j , Ṽ d j ), (V 0d j , Ṽ 0d j ) and (V dd j , Ṽ dd j ) are constructed using classical approaches [START_REF] Masson | Biorthogonal spline wavelets on the interval for the resolution of boundary problems[END_REF]. We denote these bases by {ψ d j,k , ψd j,k }, {ψ 0d j,k , ψ0d j,k } and {ψ dd j,k , ψdd j,k } respectively and for each basis, a fast wavelet transform algorithm exists with a linear complexity [START_REF] Cohen | Wavelets on the Interval and Fast Wavelet Transforms[END_REF][START_REF] Masson | Biorthogonal spline wavelets on the interval for the resolution of boundary problems[END_REF][START_REF] Monasse | Orthogonal Wavelet Bases Adapted For Partial Differential Equations With Boundary Conditions[END_REF][START_REF] Kadri-Harouna | Effective construction of divergence-free wavelets on the square[END_REF]. The approximation order of each space is linked to the number of vanishing moments of its biorthogonal wavelets. Since the space V 1 j contains polynomials up to degrees r -1, it corresponds to an approximation order r, and r vanishing moments for the biorthogonal wavelets ψd j,k . Due to the differentiation relation (2.1), the space V 0 j has for approximation order r -1, and r -1 vanishing moments for the biorthogonal wavelets ψ0d j,k . Similarly, the wavelets ψ d j,k and ψ 0d j,k have r and r + 1 vanishing moments respectively.

(ii) Divergence-free wavelet construction

Let Ω = [0, 1] 2 . To construct divergence-free wavelet bases, remark first that the divergence-free space H div (Ω) with free-slip boundary condition (1.5) is the curl of H 1 0 (Ω) stream functions [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF]:

H div (Ω) = {u = curl(Ψ) : Ψ ∈ H 1 0 (Ω)}.
By (2.3-2.4) the tensor spaces (V d j ⊗V d j ) j≥jmin provide a biorthogonal MRA of H 1 0 (Ω), and divergence-free scaling functions on Ω = [0, 1] 2 are constructed by taking the curl of scaling functions of

V d j ⊗ V d j : Φ div j,k := curl[ϕ d j,k1 ⊗ ϕ d j,k2 ] = ϕ d j,k1 ⊗ (ϕ d j,k2 ) ′ -(ϕ d j,k1 ) ′ ⊗ ϕ d j,k2 , 1 ≤ k 1 , k 2 ≤ N j -2.
Accordingly, anisotropic divergence-free wavelets on [0, 1] 2 are constructed by taking the curl of the three types of scalar anisotropic wavelets associated to V d j ⊗ V d j (j ≥ j min ):

Ψ div,1 j,k := curl[ϕ d jmin,k ⊗ ψ d j2,k2 ], Ψ div,2 j,k := curl[ψ d j1,k1 ⊗ ϕ d jmin,k ], Ψ div,3 j,k := curl[ψ d j1,k1 ⊗ ψ d j2,k2 ], for j = (j 1 , j 2 ), j 1 , j 2 ≥ j min , k = (k 1 , k 2 ) ∈ I j ,
with:

I j = {0, 1, . . . , 2 j1 -1} × {0, 1, . . . , 2 j2 -1}.
Then, one can prove [START_REF] Kadri-Harouna | Effective construction of divergence-free wavelets on the square[END_REF] that:

H div (Ω) = span{Φ div jmin,k ; k = (k 1 , k 2 ) with 1 ≤ k 1 , k 2 ≤ N jmin -2} ⊕ span{Ψ div,ε j,k ; ε = 1, 2, 3, j = (j 1 , j 2 ) with j 1 , j 2 ≥ j min , k ∈ I j },
which means that each vector function u of H div (Ω) has a unique divergence-free wavelet decomposition:

u = k c div jmin,k Φ div jmin,k + j,k ε=1,2,3 d div,ε j,k Ψ div,ε j,k , (2.6)
with the norm-equivalence:

u 2 L 2 ∼ k |c div jmin,k | 2 + j,k ε=1,2,3 |d div,ε j,k | 2
, and a linear complexity for the computation of the coefficients [START_REF] Kadri-Harouna | Effective construction of divergence-free wavelets on the square[END_REF]. For j ≥ j min , the approximation of u on the finite dimensional divergence-free space

V div j = span{Φ div j,k } 1≤k1,k2≤Nj -2 = (V d j ⊗ V 0 j ) × (V 0 j ⊗ V d j ) ∩ H div (Ω), (2.7)
writes using the orthogonal projector:

P div j (u) = k c div jmin,k Φ div jmin,k + |j|<j k ε=1,2,3 d div,ε j,k Ψ div,ε j,k . (2.8)
The approximation error is linked to the approximation order of spaces V 1 j . If spaces V 1 j contain polynomials up the degree r -1, then V 0 j contain polynomials up the degree r -2 and for all u ∈ (H s (Ω)) 2 with 0 ≤ s ≤ r -1, the following Jackson type estimation holds:

u -P div j (u) L 2 ≤ C2 -js u H s . (2.9)
For homogeneous boundary conditions, i.e. in (H 1 0 (Ω)) 2 , the divergence-free function space is slightly different, and will be denoted by:

H div,0 (Ω) = {u ∈ (H 1 0 (Ω)) 2 : ∇ • u = 0} = (H 1 0 (Ω)) 2 ∩ H div (Ω)
The space H div,0 (Ω) is a closed subspace of (H 1 0 (Ω)) 2 , then we have the following decomposition:

(H 1 0 (Ω)) 2 = H div,0 (Ω) ⊕ H div,0 (Ω) ⊥ ,
which is orthogonal for the scalar product of (

H 1 0 (Ω)) 2 : (u, v) (H 1 0 (Ω)) 2 = (∇u, ∇v) (L 2 (Ω)) 2 .
It is easy to prove [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF] that:

H div,0 (Ω) ⊥ = {(-∆) -1 ∇q : q ∈ L 2 (Ω)}, (2.10)
where v = (-∆) -1 f denotes the solution of -∆v = f , with Dirichlet homogeneous boundary conditions.

Since H div,0 (Ω) ⊂ H div (Ω), a multiresolution analysis of H div,0 (Ω) is then provided by the spaces:

V div,0 j = V div j ∩ (H 1 0 (Ω)) 2 = (V d j ⊗ V 0 j ) × (V 0 j ⊗ V d j ) ∩ (H 1 0 (Ω)) 2 ∩ H div (Ω).
By definition (2.7) of the space V div j , we have

V div,0 j = span{Φ div,0 j,k } 2≤k1,k2≤Nj -3 , with: Φ div,0 j,k := curl[ϕ dd j,k1 ⊗ ϕ dd j,k2 ] = curl[ϕ 1 j,k1 ⊗ ϕ 1 j,k2 ], 2 ≤ k 1 , k 2 ≤ N j -3.
corresponding to scaling functions on [0, 1] satisfying both Dirichlet and Neumann boundary conditions at 0 and 1 [START_REF] Monasse | Orthogonal Wavelet Bases Adapted For Partial Differential Equations With Boundary Conditions[END_REF]. Similarly, replacing in the previous construction wavelets of W d j by those of W dd j allows to construct wavelets of H div,0 (Ω). Let now P div,0 j be the L 2 -orthogonal projector from (H 1 0 (Ω)) 2 onto V div,0 j . Again, for all u ∈ (H s (Ω)) 2 with 1 ≤ s ≤ r -1, the following Jackson type estimation holds, for some C > 0:

u -P div,0 j (u) H 1 0 ≤ C 2 -j(s-1) u H s . (2.11)
2.2. Leray-Hopf projector computation. We present below the basics to compute in practice the Leray-Hopf projector P, using divergence-free wavelet bases, first in the free-slip boundary conditions case (see [START_REF] Kadri-Harouna | Helmholtz-Hodge Decomposition on [0, 1] d by Divergence-free and Curl-free Wavelets, Curves and Surfaces[END_REF] for more details), then in the case of Dirichlet boundary conditions.

For any vector field u ∈ (L 2 (Ω)) 2 , the Helmholtz-Hodge decomposition theorem states that, there exist unique u div ∈ H div (Ω) and q ∈ H 1 (Ω) with Ω q = 0, such that:

u = u div + ∇q with u div = P(u). (2.12)
This decomposition is orthogonal, and P : (L 2 (Ω)) 2 → H div (Ω) denotes the Leray-Hopf orthogonal projector. The explicit computation of u div can be obtained using the divergence-free wavelet basis constructed in Section 2 (with the unified notation Ψ div j,k ):

H div (Ω) = span{Ψ div j,k }, ∀ j, k, ∇ • Ψ div j,k = 0 and Ψ div j,k • n = 0.
Searching u div in terms of its divergence-free wavelet series:

u div = j,k d div j,k Ψ div j,k , (2.13)
and by the orthogonality Ψ div j,k ⊥ ∇q in (L 2 (Ω)) 2 , we obtain:

u, Ψ div j,k = u div , Ψ div j,k .
Accordingly the computation of coefficients (d div j,k ) is reduced to the resolution of a linear system:

M div (d div j,k ) = ( u, Ψ div j,k ), (2.14)
where M div denotes the Gram matrix of the basis {Ψ div j,k }. The computation of the matrix elements and the right hand side ( u, Ψ div j,k ) in (2.14) is described in [START_REF] Kadri-Harouna | Helmholtz-Hodge Decomposition on [0, 1] d by Divergence-free and Curl-free Wavelets, Curves and Surfaces[END_REF]. Since the Ψ div j,k are "curl" functions, the matrix M div is no more than the stiffness matrix of the 2D Laplacian operator, on the scalar wavelet basis of (V d j ⊗ V d j ) in H 1 0 (Ω). This matrix, in wavelet basis, admits an optimal (diagonal) preconditioning, and we showed in [START_REF] Kadri-Harouna | Helmholtz-Hodge Decomposition on [0, 1] d by Divergence-free and Curl-free Wavelets, Curves and Surfaces[END_REF] that the linear system (2.14) is solved with a quasi-linear complexity (with respect to the size of the wavelet coefficient matrix) using a preconditioned conjugate gradient method.

For homogeneous Dirichlet boundary condition problems, we need in this case to compute the L 2 -orthogonal Leray-Hopf projector P div,0 : (H 1 0 (Ω)) 2 → H div,0 (Ω). Since H div,0 (Ω) is also spanned by a divergence-free wavelets basis (satisfying homogeneous Dirichlet boundary conditions), the computation of P div,0 will be obtained using a similar method than the computation of P.

Let {Ψ div,0 j,k } be the divergence-free wavelet basis of H div,0 (Ω). The (L 2 (Ω)) 2orthogonal projector P div,0 is then defined by:

P div,0 (u) = j,k d div,0 j,k Ψ div,0 j,k , u ∈ (H 1 0 (Ω)) 2 , (2.15)
with

(d div,0 j,k ) = M -1 div,0 ( u, Ψ div,0 j,k L 2 ), (2.16)
and M div,0 denotes the Gram matrix of the basis {Ψ div,0 j,k }, which is symmetric positive definite.

A L 2 -orthogonal decomposition writes now in (H 1 0 (Ω)) 2 as:

u = u div,0 + u ⊥ div,0 with u div,0 = P div,0 (u) (2.17)
The orthogonal part u ⊥ div,0 = u -P div,0 (u) satisfies, for all v ∈ H div,0 (Ω):

u -P div,0 (u), v L 2 = u, v L 2 -P div,0 (u), v L 2 = u, v L 2 -u, P div,0 (v) L 2 = 0,
Using De Rhams theorem [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF], there exists p ∈ H 1 (Ω) such that u -P div,0 (u) = ∇p. Every function u ∈ (H 1 0 (Ω)) 2 can be decomposed as:

u = P div,0 (u) + (u -P div,0 (u)) = P div,0 (u) + ∇p. (2.18)
where ∇p satisfies homogeneous boundary conditions.

3. Modified projection method based on divergence-free wavelets. The purpose in this section is to introduce a wavelet based numerical scheme for incompressible viscous flows calculations. The method can be seen as a variant of the projection method [START_REF] Bell | A second-order projection method for the incompressible Navier-Stokes equations[END_REF][START_REF] Chorin | Numerical simulation of the Navier-Stokes equation[END_REF][START_REF] Kim | Application of a fractional-step method to incompressible Navier-Stokes equations[END_REF][START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires II[END_REF], where we replace the operator splitting by the exact Helmholtz-Hodge decomposition of the intermediate velocity field, computed using divergence-free wavelets. Then, we prevent some numerical difficulties and drawbacks related to the computation of the pressure at each time step with artificial boundary conditions done in the classical approaches [START_REF] Bell | A second-order projection method for the incompressible Navier-Stokes equations[END_REF][START_REF] Kim | Application of a fractional-step method to incompressible Navier-Stokes equations[END_REF].

3.1. General principles of divergence-free wavelet schemes for the Stokes equations. The use of divergence-free wavelet bases in incompressible flow calculations began with the works of Urban et al [START_REF] Dahmen | Adaptive Wavelet Methods-Basic Concepts and Applications to the Stokes Problem, Wavelet Analysis-Twenty Years Developments[END_REF][START_REF] Urban | Using divergence-free wavelets for the numerical solution of the Stokes problem[END_REF], for the resolution of the stationary Stokes problem:

   -ν∆v + ∇p = f , ∇ • v = 0, (3.1)
in Ω = [0, 1] 2 , with periodic or homogeneous Dirichlet boundary conditions.

The main advantage of using divergence-free wavelet basis in the resolution of Stokes equations is the direct representation of the incompressibility constraint of the flow. Following Urban's works [START_REF] Urban | Using divergence-free wavelets for the numerical solution of the Stokes problem[END_REF][START_REF] Urban | Wavelets in Numerical Simulation[END_REF] problem (3.1) is solved using a variational approach, and a Galerkin approximation using divergence-free wavelets as trial functions. The velocity field v is searched in terms of its divergence-free wavelet coefficients:

v(x) = j,k d j,k Ψ div j,k (x). (3.2)
Replacing (3.2) into (3.1), the computation of coefficients d j,k is done by solving the linear system:

ν j,k d j,k ∇Ψ div j,k , ∇Ψ div j ′ ,k ′ = f , Ψ div j ′ ,k ′ , ∀ j ′ , k ′ . (3.3)
The divergence-free wavelet stiffness matrix [ ∇Ψ div j,k , ∇Ψ div j ′ ,k ′ ] is symmetric and the associated bilinear form is coercive [START_REF] Urban | Using divergence-free wavelets for the numerical solution of the Stokes problem[END_REF]. The problem is thus reduced to an elliptic problem on some finite-dimensional divergence-free function space V div j : therefore standard error estimations hold, derived from the projection error (2.9, 2.11). In addition, the formulation (3.3) has the advantage to eliminate directly the pressure p which is computed by a post processing procedure [START_REF] Urban | Using divergence-free wavelets for the numerical solution of the Stokes problem[END_REF].

In comparison with classical approaches based on finite differences, finite elements or (non divergence-free) wavelet methods [START_REF] Canuto | Stabilized wavelet approximations of the Stokes problem[END_REF][START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF], equation (3.3) presents the advantage of reducing the number of degree of freedom: only-scalar-coefficients {d j,k } are computed instead of one type of coefficients per components of the velocity v. Moreover, adaptive strategies can be applied to only compute significant wavelet coefficients, and optimal preconditioning for the stiffness matrix can be provided explicitly [START_REF] Cohen | Numerical Analysis of Wavelet Methods[END_REF][START_REF] Kadri-Harouna | Helmholtz-Hodge Decomposition on [0, 1] d by Divergence-free and Curl-free Wavelets, Curves and Surfaces[END_REF].

For unsteady problems, recently Stevenson [START_REF] Stevenson | Divergence-free wavelet bases on the hypercube: Free-slip boundary conditions, and applications for solving the instationary Stokes equations[END_REF] proposed a new theoretical variational formulation of the Stokes equations. The method of Stevenson is an extension of Urban's method to the unsteady problem using divergence-free wavelets satisfying a free-slip boundary condition.

We adopt in the present article another point of view, and the method we present in next coming sections will use a Galerkin formulation with standard wavelet bases for the approximation of the velocity, and the Leray-Hopf projector P computed with divergence-free wavelets. This leads to a new projection method formulation, using divergence-free wavelets satisfying Dirichlet boundary condition. The advantage is that classical wavelet methods can be used to solve the diffusion problem and the incompressibility constraint is incorporated via the projector P. The approach will include one phase devoted to the temporal discretization and a second one to spatial discretization.

3.2. Modified projection method for unsteady Stokes equations. We consider in this section the unsteady Stokes problem, with no-slip boundary conditions:

           ∂v ∂t -ν∆v + ∇p = f , v = 0 on ∂Ω, ∇ • v = 0. (3.4)
The time discretization of system (3.4) is obtained through a finite difference method. Without loss of generality, given a time step δt and considering the approximation v n (x) ≈ v(x, nδt), the backward Euler scheme leads to:

v n+1 -v n δt -ν∆v n+1 + ∇p n+1 = f n , ∇ • v n+1 = 0. (3.5)
However, scheme (3.5) is inefficient since it requires, at each time step, the solution of coupled equations for (v n+1 , p n+1 ). Now, let us introduce new variables ṽn ∈ (H 1 0 (Ω)) 2 and Φ n by setting ṽn+1 = v n+1 + ∇Φ n+1 . One can prove that ṽn+1 verifies the following system:

ṽn+1 -v n δt -ν∆ṽ n+1 + ∇[p n+1 - 1 δt Φ n+1 + ν∆Φ n+1 ] = f n , v n+1 = P div,0 (ṽ n+1 ).
(3.6) This equation can be split by defining the pressure p n+1 such that:

p n+1 - 1 δt Φ n+1 + ν∆Φ n+1 = 0, (3.7)
and (3.6) reduces to

(1 -νδt∆)ṽ n+1 = v n + δtf n , v n+1 = P div,0 (ṽ n+1 ). (3.8)
Remark that the change of variable of this formulation is similar to a Gauge formulation applied to the discrete equation (3.5). The pressure is replaced by a Gauge variable Φ, linked to p by equation (3.7), which replaces a continuous heat equation appearing in the continuous formulation [START_REF] Wang | Convergence of Gauge method for incompressible flow[END_REF]. However, the difference between our method and the Gauge method lies in the choice of boundary conditions. Instead of imposing ∇Φ n+1 • τ = 0 or ∇Φ n+1 • n = 0 as in the Gauge method, we choose boundary conditions on the auxiliary field ṽn+1 equal to those of v n , which leads to: ṽn+1 = 0 and ∇Φ n+1 = 0, on ∂Ω, then we get the desired boundary conditions v n = 0 on ∂Ω (note that in our method Φ has no physical meaning). The resolution of (3.8) requires the resolution of a heat equation followed by a projection step, described in Section 2.2. The spacial discretization is obtained through a variational Galerkin method on suitable tensorial wavelet bases.

Spatial discretization of Stokes equations.

The spatial approximation of the velocity will be performed using the multiresolution analysis of (H 1 0 (Ω)) 2 that contains the divergence-free wavelets. This multiresolution analysis is constituted by the vector spaces

V d j = (V d j ⊗ V 0d j ) × (V 0d j ⊗ V d j ).
We now consider a Galerkin formulation of equations (3.8): at a given resolution j, the components of the approximate solution v n j = (v n 1,j , v n 2,j ) are searched under the form of finite dimensional wavelet series in V d j :

v n 1,j = |j|<j,k d 1,n j,k ψ d j1,k1 ⊗ ψ 0d j2,k2 and v n 2,j = |j|<j,k d 2,n j,k ψ 0d j1,k1 ⊗ ψ d j2,k2 , (3.9)
and similarly for ṽn

j with coefficients [ d1,n j,k ] and [ d2,n j,k ].
Denoting by a δt the bilinear form:

a δt (v, w) = Ω vw + νδt Ω ∇v • ∇w , (3.10)
and by P div,0 j the restriction of the Leray projector P div,0 to V d j , the Galerkin formulation of equations (3.8) writes at each time-step n: knowing the approximate solution v n j ∈ V div,0 j and the r.h.s.

f n j ∈ V d j , find ṽn+1 j ∈ V d j , v n+1 j ∈ V div,0 j s.t. ∀w j ∈ V d j , a δt (ṽ n+1 j , w j ) = Ω v n j w j + δt Ω f n j w j , v n+1 j = P div,0 j (ṽ n+1 j ). (3.11)
Now the tensorial structure of wavelet bases of V d j allows to factorize at each time step the stiffness matrix of the discrete heat operator (1 -δtν∆), like in [START_REF] Charton | A Pseudo-Wavelet Scheme for the Two-Dimensional Navier-Stokes Equations[END_REF]. Such factorization remains in the context of alternated direction implicit (ADI) methods, and consists, for small α, to use the approximation:

(1 -α∆) ≈ (1 -α ∂ 2 ∂x 2 )(1 -α ∂ 2 ∂y 2 ). (3.12)
Remark that this ADI factorization (3.12) neglects the fourth order derivatives and this has some implication for the unsteady boundary conditions. In such case, the spatial discretization accuracy has to be improved by increasing the number of vanishing moments of the wavelet family. In this work we will only consider steady Dirichlet boundary condition. Thus, in (3.11) we only have to invert the Galerkin matrix of the one-dimensional heat operator (1 -νδt ∂ 2 ∂x 2 ), which is done once before starting the time integration procedure. Finally, the computation of matrices of coefficients [ d1,n+1 j,k ] and [ d2,n+1 j,k ] from those of v n reduces to solve the following systems:

A d δt [ d1,n+1 j,k ] A 0d δt = M d [d 1,n j,k ] M 0d + δt M d f n j1 M 0d (3.13)
and

A 0d δt [ d2,n+1 j,k ] A d δt = M 0d [d 2,n j,k ] M d + δt M 0d f n j2 M d , (3.14)
where A ε δt and M ε correspond respectively to the stiffness matrix of operator (1νδt ∂ 2 ∂x 2 ) and the Gram matrix of the one-dimensional wavelet basis of {V ε j } ε=d,0d :

A d δt = [ ψ d j,k , ψ d j ′ ,k ′ + νδt (ψ d j,k ) ′ , (ψ d j ′ ,k ′ ) ′ ] and M d = [ ψ d j,k , ψ d j ′ ,k ′ ],
and

A 0d δt = [ ψ 0d j,k , ψ 0d j ′ ,k ′ + νδt (ψ 0d j,k ) ′ , (ψ 0d j ′ ,k ′ ) ′ ] and M 0d = [ ψ 0d j,k , ψ 0d j ′ ,k ′ ].
Elements of the form

ψ d j,k , ψ d j ′ ,k ′ or (ψ d j,k ) ′ , (ψ d j ′ ,k ′ )
′ are analytically computed, by solving eigenvalue problems [START_REF] Beylkin | On the representation of operator in bases of compactly supported wavelets[END_REF][START_REF] Monasse | Orthogonal Wavelet Bases Adapted For Partial Differential Equations With Boundary Conditions[END_REF].

We summarize below the resolution algorithm. Starting with an initial velocity ṽ0 (x) = v 0 (x) = v(0, x), compute its wavelet coefficients [ d1,0 j,k ] and [ d2,0 j,k ] in V d j [START_REF] Kadri-Harouna | Ondelettes pour la prise en compte de conditions aux limites en turbulence incompressible[END_REF]. For 1 ≤ n ≤ N , repeat:

Step 1: Find [ d1,n+1 j,k ] and [ d2,n+1 j,k ] solution of (3.13, 3.14).

Step 2: Find [d div,n+1 j,k ] solution of

M d [d div,n+1 j,k ] R d + R d [d div,n+1 j,k ] M d = M d [ d1,n+1 j,k ] A 0 d -(A 0 d ) T [ d2,n+1 j,k ] M d ,
where

R d = [ (ψ d j,k ) ′ , (ψ d j ′ ,k ′ ) ′ ] and A 0 d = [ ψ 0d j,k , (ψ d j ′ ,k ′ ) ′ ] [22].
Step

3: Compute [d 1,n+1 j,k ] and [d 2,n+1 j,k ] from [d div,n+1 j,k
] using the change of basis between {(ψ d j,k ) ′ } and {ψ 0d j,k } [START_REF] Kadri-Harouna | Effective construction of divergence-free wavelets on the square[END_REF].

As the matrices A 0d δt and A 0 δt are inverted once for all before starting the algorithm, Step 1 is thus only a matrix-matrix multiplication. If j denotes the maximal space resolution, i.e, N = (2 j + 1) 2 grid points, the theoretical complexity of this step is at most O(2 3j ). Step 2 corresponds to v n+1 j = P div,0 j (ṽ n+1 j ), this is solved with a preconditioned conjugate gradient method. Since M div,0 is part of M div , which is the stiffness matrix of the 2D Laplacian operator on the scalar wavelet basis of (V d j ⊗ V d j ), the complexity of Step 2 can not exceed the complexity of a Poisson solver on the same wavelet basis (quasi-linear complexity [START_REF] Cohen | Numerical Analysis of Wavelet Methods[END_REF][START_REF] Cohen | Wavelet methods for second order elliptic problems -preconditioning and adaptivity[END_REF][START_REF] Kadri-Harouna | Helmholtz-Hodge Decomposition on [0, 1] d by Divergence-free and Curl-free Wavelets, Curves and Surfaces[END_REF]). The last step is a change of basis, with a linear complexity. Finally the theoretical complexity of the whole algorithm is at most O(2 3j ), if the whole set of wavelet coefficients is considered. In practice the effective complexity is much lower, due to the sparse structure of both the wavelet coefficients and wavelet operators.

Stability and consistency analysis.

This section analyzes the stability and consistency of the modified projection method for the unsteady Stokes problem. For sake of simplicity we take f n = 0 in (3.8) and suppose that v n , ṽn ∈ (H 1 0 (Ω)) 2 are regular enough. Proof. To prove the stability of scheme (3.8), a standard energy estimate will be used with ṽn+1 as test function, thanks to the Dirichlet boundary conditions on ṽn+1 . Taking the inner product of equation (3.8) with 2ṽ n+1 , we obtain

ṽn+1 2 L 2 + ṽn+1 -v n 2 L 2 -v n 2 L 2 + 2νδt ∇ṽ n+1 2 L 2 = 0. (3.15) Since ṽn+1 = v n+1 + ∇Φ n+1 , which is an L 2 -orthogonal decomposition in H 1 0 (Ω) 2 , equation (3.15) simplifies: v n+1 2 L 2 + v n+1 -v n 2 L 2 + 2 ∇Φ n+1 2 2 -v n 2 L 2 + 2νδt ∇ṽ n+1 2 L 2 = 0. which leads to: v n+1 2 L 2 -v n 2 L 2 + 2 ∇Φ n+1 2 L 2 + 2νδt ∇ṽ n+1 2 L 2 ≤ 0 (3.16)
and completes the proof.

The convergence of the method is a consequence of the following theorem:

Theorem 3.2. Let v be a smooth solution of Stokes equations with smooth initial data v 0 (x) and let v δt be the numerical solution of the semi-discrete modified projection method (3.7) and (3.8), then:

v -v δt L ∞ ([0,T ];L 2 ) ≤ C 1 δt, ∇v -∇v δt L ∞ ([0,T ];L 2 ) ≤ C 2 δt 1/2 .
Remark that the modified projection method has the same convergence order than the backward Euler scheme, which means that the projection step preserves the time discretization order. A similar result holds for more accurate time discretization schemes like the Crank-Nicolson one.

Proof. Let v n+1 be the solution of (3.7) and (3.8) computed from v n (x) = v(x, nδt). Let ǫ n+1 = v(x, (n + 1)δt)v n+1 be the consistency error. Thus, ǫ n+1 is linked to ṽn+1 by: ṽn+1 = v(•, (n + 1)δt) + ∇Φ n+1ǫ n+1 .

(3.17)

Replacing (3.17) in (3.6) we obtain:

-ǫ n+1 + νδt∆ǫ n+1 + δt∇p n+1 + v(x, (n + 1)δt) -v n -νδt∆v(x, (n + 1)δt) = 0. (3.18)
Using now the Taylor series expansions:

v(x, (n + 1)δt) = v n (x) + δt ∂v ∂t (x, nδt) + O(δt 2 ), ∆v(x, (n + 1)δt) = ∆v(x, nδt) + O(δt) then: v(x, (n + 1)δt) -v n (x) -νδt∆v(x, (n + 1)δt) = δt ∂v ∂t (x, nδt) -ν∆v(x, nδt) + O(δt 2 ).
Since v is the exact solution of Stokes equation, we have ∂v ∂t -ν∆v = -∇p, and (3.18) rewrites:

-ǫ n+1 + νδt∆ǫ n+1 = δt∇[p(x, nδt) -p n+1 ] + O(δt 2 ). (3.19)
By definition, ǫ n+1 ∈ (H 1 0 (Ω)) 2 and is divergence-free: ∇ • ǫ n+1 = 0. Taking -ǫ n+1 as a test function in (3.19) yields:

ǫ n+1 2 L 2 + νδt ∇ǫ n+1 2 L 2 = O(δt 2 )( Ω ǫ n+1 ). Since Ω = [0, 1] 2 is bounded, then Ω ǫ n+1 ≤ ǫ n+1 L 2 and
ǫ n+1 2 L 2 + νδt ∇ǫ n+1 2 L 2 ≤ Cδt 2 ǫ n+1 L 2 ,
which implies

ǫ n+1 L 2 ≤ Cδt 2 , ∇ǫ n+1 L 2 ≤ C √ ν δt 3/2 .
The method is then consistent and using the stability result, it converges with a first order convergence in time. The second estimate follows from the fact that one can prove that ∇v n L 2 is bounded using (3.5).

The spacial consistency error depends on the regularity s of the solution v and the approximation order r provided by spaces V div,0 j (2.11), sucht that 0 ≤ s ≤ r -1.

Theorem 3.3. Let v be a s-smooth solution of Stokes equations with smooth initial data v 0 (x) and let v n j be the wavelet numerical solution in V div,0 j of the semidiscrete modified projection method (3.11) computed from v 0 j = P div,0 j (v 0 ), then:

v -v n j L ∞ ([0,T ];L 2 ) ≤ C 1 (δt + 2 -js ), (3.20)
√ δt ∇v -∇v n j L ∞ ([0,T ];L 2 ) ≤ C 2 (δt + 2 -j(s-1) ). (3.21)
where the constant C 1 , C 2 are independant of δt, n and j.

Proof. The proof can be viewed as a modified Cea's lemma. Let ǫ n j and ǫn j be the errors defined by ǫ n j = v nv n j and ǫn j = ṽnṽn j , where ṽn , v n are the solutions of (3.8) and ṽn j , v n j the solutions of (3.11), computed at resolution j and time nδt. Remark that ∇ • ǫ n j = 0. At step n, we suppose that v n j = P div,0 j (v n ), and we study the consistency error ǫ n+1 j provided by (3.11). We have:

a δt (ṽ n+1 , w) = Ω v n w, v n+1 = P div,0 (ṽ n+1 ), ∀ w ∈ (H 1 0 (Ω)) 2 , a δt (ṽ n+1 j , w j ) = Ω v n j w j , v n+1 j = P div,0 j (ṽ n+1 j ), ∀ w j ∈ V d j .
Taking w = w j ∈ V d j , we obtain:

a δt (ǫ n+1 j , w j ) = Ω
ǫ n j w j and ǫ n+1 j = P div,0 (ṽ n+1 ) -P div,0 j (ṽ n+1 j

).

Now, considering w j ∈ (V div,0 j ) ⊥ ∩ V d j , then w j = ∇φ for some φ and one has Ω ǫ n j w j = 0, which implies: a δt (ǫ n+1 j , w j ) = 0.

Moreover since v n j = P div,0 j (v n ), then Ω v n w j = Ω v n j w j for all w j ∈ V div,0 j , which leads to a δt (ǫ n+1 j , w j ) = 0. Finally, we have

a δt (ǫ n+1 j , w j ) = 0, ∀w j ∈ V d j . (3.22)
Taking w j = u j -ṽn+1 j in (3.22), which satisfies w j ∈ V d j for any u j ∈ V d j , we get:

0 = a δt (ǫ n+1 j , u j -ṽn+1 j ) = a δt (ǫ n+1 j , u j -ṽn+1 ) + a δt (ǫ n+1 j , ṽn+1 -ṽn+1 j ),
which can be rewritten as:

a δt (ǫ n+1 j , ǫn+1 j ) = a δt (ǫ n+1 j , ṽn+1 -u j ). (3.23)
From the continuity of the bilinear for a δt (3.10), we deduce that:

a δt (ǫ n+1 j , ṽn+1 -u j ) ≤ (1 + νδt) ǫn+1 j H 1 ṽn+1 -u j H 1 , ∀ u j ∈ V d j .
and

a δt (ǫ n+1 j , ṽn+1 -u j ) ≤ ǫn+1 j L 2 ṽn+1 -u j L 2 + νδt ǫn+1 j L 2 ṽn+1 -u j H 2 , ∀ u j ∈ V d j ,
On the other hand, using the Poincaré-Friedrichs [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF] inequality leads to:

νδt ǫn+1 j 2 H 1 a δt (ǫ n+1 j , ǫn+1 j ),
and we also have:

ǫn+1 j 2 L 2 ≤ a δt (ǫ n+1 j , ǫn+1 j ).
Then, (3.23) gives:

√ νδt ǫn+1 j H 1 ṽn+1 -u j H 1 , ∀ u j ∈ V d j .
For u j = P d j (ṽ n+1 ) (projection onto V d j ), the usual Jackson's inequality yields:

√ νδt ǫn+1 j H 1 2 -(j-1)s ṽn+1 H s 2 -(j-1)s ṽ0 H s , (3.24)
where we have used the stability properties (3.15, 3.16).

Similarly, one has:

ǫn+1 j L 2 2 -js + νδt 2 -j(s-2) ṽ0 H s . (3.25) Since V div,0 j = V d j ∩ H div,0 (
Ω), see [START_REF] Kadri-Harouna | Effective construction of divergence-free wavelets on the square[END_REF], we have P div,0 j (ṽ n+1 j ) = P div,0 (ṽ n+1 j

) and the two errors are linked by the relation:

ǫ n+1 j = P div,0 (ǫ n+1 j ) ⇒ ǫ n+1 j L 2 ≤ ǫn+1 j L 2 and ǫ n+1 j H 1 ≤ ǫn+1 j H 1 .
Finally, to prove (3.20) it suffices to write, for all n:

v -v n j L 2 = v -v n + v n -v n j L 2 ≤ v -v n L 2 + v n -v n j L 2 (δt + 2 -js + νδt 2 -j(s-2) )
(δt + 2 -js ), using theorem 3.2 and (3.25), as ν2 -j(s-2) ≤ 1. Similar arguments allow to prove (3.21) using the H 1 norm, which completes the proof.

Example.

To investigate the convergence rates of the wavelet modified projection method for the unsteady Stokes problem (3.1), two numerical tests are conducted: the first one to evaluate the time discretization error and the second one to evaluate the spatial discretization error. As exact solution, we used:

   v 1 (x, y, t) = 1 8π 2 ν (1 -e -8π 2 νt )[cos(2πx) sin(2πy) -sin(2πy)], v 2 (x, y, t) = -1 8π 2 ν (1 -e -8π 2 νt )[sin(2πx) cos(2πy) - p(x, y, t) = 0.5e -t [cos(2πx) -cos(2πy)]. (3.26)
This solution satisfies Dirichlet homogeneous boundary conditions v |∂Ω = 0, where Ω = [0, 1] 2 . The right-hand side term f is computed appropriately to ensure that (3.26) is the exact solution of system (3.1). The 1D wavelet generators of (V 1 j , Ṽ 1 j ) are biorthogonal splines with r = r = 3, which means that the approximation order of V 0 j is 2; this corresponds to a bidimensional approximation order of at most 2.

For the time discretization we chose two implicit methods that are easy to implement: backward-Euler and Crank-Nicholson schemes. We recall that the order of these schemes are 1 and 2 respectively. However, any accurate numerical scheme would have been be used without restriction on the order. Tab. 3.1 shows the L 2 -error between the exact solution projected onto V j (with a space resolution fixed at j = 10 and ν = 2 -j ), and the numerical solution of (3.11) to the theoretical order (1 or 2) of the chosen time scheme. For the evaluation of the spatial discretization error, we considered a regular grid. Following theorem 3.3, we have to chose a time step very small compared to the spatial resolution: we took δt = 5.10 -4 and ν = 2 -10 . Tab. 3.2 shows the spatial error at grid points for the simulation final time T = 2. Remark that, as the solution is C ∞ , the convergence rate given by Tab. of our wavelet family (equal to 2 in our spline approximation for the L 2 -error), and we lose one order for the H 1 -error.

4. Modified projection method for unsteady Navier-Stokes equations.

4.1. Divergence-free wavelet schemes for Navier-Stokes equations. Divergence-free wavelet schemes for the numerical resolution of Navier-Stokes equations were first introduced by Deriaz and Perrier [START_REF] Deriaz | Direct Numerical Simulation of Turbulence using divergence-free wavelets[END_REF][START_REF] Deriaz | Divergence-free and curl-free wavelets in 2D and 3D, application to turbulent flows[END_REF], in the case of periodic boundary conditions. In this section, we will consider Dirichlet boundary conditions, and we will present a new scheme for Navier-Stokes equations, in the projection method setting [START_REF] Chorin | Numerical simulation of the Navier-Stokes equation[END_REF][START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires II[END_REF], but replacing the original projection step by the divergence-free wavelet Leray projector, as in Section 3.1.

Temporal discretization of Navier-Stokes equations.

Since the original works of Chorin [START_REF] Chorin | Numerical simulation of the Navier-Stokes equation[END_REF] and Temam [START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires II[END_REF], the projection method has become highly popular for the numerical resolution of Navier-Stokes equations in velocity pressure formulation, notably with physical, such as no-slip, boundary conditions. Several versions of the projection method exist, according to the chosen pressure boundary condition [START_REF] Bell | A second-order projection method for the incompressible Navier-Stokes equations[END_REF]. Without loss of generality, we focus here on the second order boundary approximation in time, called projection method with accurate pressure boundary condition [START_REF] Guo-Liu | Projection Method I: Convergence and Numerical Boundary Layers[END_REF]. The principle (time) steps of this method is summarized bellow [START_REF] Guo-Liu | Projection Method I: Convergence and Numerical Boundary Layers[END_REF]:

• Prediction step: knowing v n ∈ H div,0 (Ω), compute an intermediate velocity field v * such that:    v * -v n δt + (v n+1/2 • ∇)v n+1/2 = ν∆ v * +v n 2 , v * = 0, on ∂Ω, (4.1)
with

(v n+1/2 • ∇)v n+1/2 = 3 2 (v n • ∇)v n - 1 2 (v n-1 • ∇)v n-1 . (4.2)
• Correction step: project v * onto the divergence-free functions space to get

v n+1 ∈ H div,0 (Ω):            v n+1 = v * -δt∇p n+1/2 , ∇ • v n+1 = 0, ∇p n+1/2 • n = -n • [∇ × (∇ × v * )], on ∂Ω, (4.3)
with ∇p n+1/2 = 3 2 ∇p n -1 2 ∇p n-1 .

In classical approaches [START_REF] Bell | A second-order projection method for the incompressible Navier-Stokes equations[END_REF][START_REF] Chorin | Numerical simulation of the Navier-Stokes equation[END_REF][START_REF] Kim | Application of a fractional-step method to incompressible Navier-Stokes equations[END_REF][START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires II[END_REF] the computation of the velocity v n+1 needs first to solve a Poisson equation:

δt∆p n+1/2 = ∇ • v * , (4.4)
with suitable boundary condition. Since v * and v n+1 vanish on the boundary ∂Ω, we should have ∇p n+1/2 = 0 on ∂Ω, which means:

∇p n+1/2 • n = 0 and ∇p n+1/2 • τ = 0 on ∂Ω. (4.5)
Any Poisson solver can not incorporate both boundary conditions (4.5) at the same time, the problem becoming overdetermined. In practice, to solve (4.4), only the Neumann boundary condition for the pressure ∇p n+1/2 • n = 0 is considered [START_REF] Guo-Liu | Projection Method I: Convergence and Numerical Boundary Layers[END_REF][START_REF] Kim | Application of a fractional-step method to incompressible Navier-Stokes equations[END_REF], which leads to a loss of accuracy at the boundary [START_REF] Guo-Liu | Projection Method I: Convergence and Numerical Boundary Layers[END_REF]. Moreover, preferring one of these two boundary conditions (4.5) distorts the Dirichlet homogeneous boundary condition satisfied by v:

∇p n+1/2 • n = 0 ⇒ v n+1 • τ = ∇p n+1/2 • τ = 0, and
∇p n+1/2 • τ = 0 ⇒ v n+1 • n = ∇p n+1/2 • n = 0.
As done for the Stokes equations, the objective of the following paragraph will be to define a divergence-free wavelet based modified projection method for Navier-Stokes equation that moreover satisfies (4.5), this will be a major difference with the methods based on a Poisson solver.

Analyzing differently the problem, one can take advantage of the boundary conditions like (4.5), the Helmholtz-Hodge decomposition in (H 1 0 (Ω)) 2 and the new construction of divergence-free wavelets with boundary conditions, to derive a new correction step for (4.1). Indeed, let Φ n+1/2 be a scalar potential in L 2 (Ω) satisfying:

ṽn+1 = v n+1 + ∇Φ n+1/2 , with ṽn+1 ∈ (H 1 0 (Ω)) 2 . ( 4
   ṽn+1 -v n δt + (v n+1/2 • ∇)v n+1/2 = ν∆ ṽn+1 +v n 2
, ṽn+1 = 0, on ∂Ω.

(4.7)

• Correction step:

       v n+1 = P div,0 (ṽ n+1 ), p n+1/2 = 1 δt Φ n+1/2 -ν 2 ∆Φ n+1/2 . (4.8)
where P div,0 denotes the L 2 -orthogonal projector from (H 1 0 (Ω)) 2 onto H div,0 (Ω) introduced in Section 2.2. Remark that the Navier-Stokes formulation (4.7) is no more than a change of variables, whereas the classical projection method is an operator splitting, which implies a loss of precision.

The spatial discretization uses a Wavelet Galerkin formulation as for the Stokes equations in Section 3.1: the components of v n and v * are searched in the form of a finite wavelet series defined by (3.9). Following [START_REF] Bell | A second-order projection method for the incompressible Navier-Stokes equations[END_REF][START_REF] Kim | Application of a fractional-step method to incompressible Navier-Stokes equations[END_REF], the nonlinear term (v n+1/2 • ∇)v n+1/2 (4.2) is computed explicitly with finite differences on the mesh grid points, and this explicit treatment imposes a CFL condition on the time step, see [START_REF] Deriaz | Direct Numerical Simulation of Turbulence using divergence-free wavelets[END_REF]. The order of this finite difference scheme has to be at least equal to the approximation order of the wavelet approximation (namely the parameter r), to not reduce the global spatial accuracy of the method.

Exemple.

To investigate the accuracy and spatial convergence rate of the modified projection method (4.7) and (4.8), we performed a convergence study on regular grids. The exact solution is that given in (3.26) with ν = 1 (using a suitable forcing term). As for the Stokes problem, the 1D wavelet generators of (V 1 j , Ṽ 1 j ) are biorthogonal splines with r = r = 3, which corresponds to a bidimensional spatial approximation order of at most 2. We used a time scheme of order 2 (Crank Nicolson), with a time step δt = 5.10 -4 , chosen small enough such that the error induced by the time discretization is negligible with respect to the spatial discretization error. In Fig 4 .1 we plotted the L 2 , L ∞ and H 1 norms of the velocity error in terms of the grid step in log-log scale, at the simulation final time T = 1. As it was stated in Th.3.3 for the Stokes problem, since the solution is C ∞ , Fig 4 .1 gives a spatial convergence rate which saturates to the number of vanishing moments of our wavelet family (equal to 2 in our spline approximation for the L 2 -error), and we lose one order for the H 1 -error.

Lid driven cavity flow.

To validate the divergence-free wavelet projection method on Navier-Stokes equations, we focus on the classical lid-driven cavity problem on regular grids. This problem has been investigated by many authors since the pioneer work of [START_REF] Botella | Benchmark spectral results on the lid-driven cavity flow[END_REF][START_REF] Ghia | High-Re solutions for incompressible flows using Navier-Stokes equations and a multigrid method[END_REF]. Recently, Bruneau and Saad [START_REF] Bruneau | The 2D lid-driven cavity problem revised[END_REF] provide new simulations at high Reynolds numbers, obtaining highly accurate benchmark results, using a multigrid solver with a special emphasis on the discretization of the convection term, for which a high space resolution is used: j = 10 or j = 11.

The objective in this section is to compare the results obtained with the scheme (4.7)-(4.8) to those of [START_REF] Botella | Benchmark spectral results on the lid-driven cavity flow[END_REF][START_REF] Bruneau | The 2D lid-driven cavity problem revised[END_REF][START_REF] Ghia | High-Re solutions for incompressible flows using Navier-Stokes equations and a multigrid method[END_REF], in order to evaluate the accuracy and performance of this new method.

The governing equations are the Navier-Stokes equations the (free-slip) non homogenous boundary condition on Γ = ∂Ω:

v |Γ = g(x, y, t) = g(x, 1, t) = (-1, 0) ∀t ∈]0, T [, ∀x ∈]0, 1[, g(x, y, t) = (0, 0) ∀t ∈]0, T [, ∀(x, y) ∈ Γ, y = 1. (4.9)
The initial velocity being chosen as v(x, y, 0) = 0 in Ω and the Reynolds number corresponds to Re = 1 ν in equation (1.1).

Boundary condition (4.9) means that for all t > 0, the edge scaling functions and wavelet coefficients of v can be computed once from those of g. In practice, to incorporate (4.9) on v n , it suffices to impose ṽ|Γ = g in (4.7) and replace P div,0 by P div in (4.8). Since the new horizontal velocity ṽ1 does not satisfy homogeneous Dirichlet boundary condition, an homogenization technique is used for this component to solve (4.7), see [START_REF] Kadri-Harouna | Ondelettes pour la prise en compte de conditions aux limites en turbulence incompressible[END_REF]. The wavelet basis generators of (V 1 j , Ṽ 1 j ) are biorthogonal splines with three vanishing moments for both space family: r = r = 3. The advection term (v n+1/2 • ∇)v n+1/2 is computed with a finite difference method of order 3 on a regular grid, the same approximation order as that of the wavelet spaces. For the steady convergence state, Fig. 4.2 shows the plot of the middle horizontal and vertical profiles of the velocity obtained with the present method for j = 7 spatial resolution and Re = 1000, and compared to the results of [START_REF] Bruneau | The 2D lid-driven cavity problem revised[END_REF] obtained with j = 10 spatial resolution. The vorticity contour and the divergence-free scaling functions coefficient for this Reynolds number Re = 1000 are plotted on Fig. 4.3 and their values again confirm the convergence of the method. Tab. 4.1 and Tab. 4.2 show the values of these profiles computed with j = 7 and j = 8 spatial resolution for present method, j = 10 spatial resolution for Bruneau and Saad [START_REF] Bruneau | The 2D lid-driven cavity problem revised[END_REF], compared to pioneers results of [START_REF] Ghia | High-Re solutions for incompressible flows using Navier-Stokes equations and a multigrid method[END_REF] and [START_REF] Botella | Benchmark spectral results on the lid-driven cavity flow[END_REF]. In [START_REF] Botella | Benchmark spectral results on the lid-driven cavity flow[END_REF], the spatial discretization is done using a spectral method with N = 160 Chebyshev polynomials.

For the simulation at moderate Reynolds number Re = 10000, to avoid corner problem, a regularized velocity v(x, 1, t) = (-16x 2 (1 -x) 2 , 0) is used. The value Re = 10000 is a critical test case, many studies were performed to see for instance if a steady solution is achieved or not, or if it has periodic behavior [START_REF] Bruneau | The 2D lid-driven cavity problem revised[END_REF]. The results obtained with the divergence-free wavelet based projection method prove that it remains stable and give rise to a numerical solution. Fig 4.5 shows the divergence-free scaling function coefficients contours, for the simulation time T = 80, which correspond to the contours of the stream function coefficients on the scalar multiresolution analysis

V 1 j ⊗ V 1 j (j = 8
). The associated divergence-free wavelet coefficients isovalues are also plotted on Fig. 4.5, which, in this case, gives evidence of the sparsity of the solution, in the wavelet representation. This sparse repartition of the wavelet coefficients is confirmed by Fig. 4.6 which produces the repartition of the divergence-free wavelet coefficients greater, in absolute value, than a thresholding parameter ǫ. (Bruneau and Saad [4]): Re = 1000 and j = 7. Since the wavelet based numerical schemes provide sparse representations of the Navier-Stokes solutions [START_REF] Schneider | Wavelet methods in computational fluid dynamics[END_REF], an important issue to improve the numerical complexity should be to implement dynamical adaptive techniques. At present, the Navier-Stokes simulation for the lid driven cavity needs about 2s of CPU time to run one iteration, using 512 2 grid points on the square [0, 1] 2 , with a Matlab code on a conventional computer. Adaptive procedure should be considered with the present schemes (4.7) and (4.8), which is confirmed by Fig. 4.4 where only almost 1% of the divergence-free wavelet coefficients are retained for small values of ǫ. To implement adaptive method on schemes (4.7) and (4.8), first notice that (4.7) is a classical elliptic equation for the velocity ṽn+1 , for which optimal adaptive wavelet methods have been designed [START_REF] Cohen | Numerical Analysis of Wavelet Methods[END_REF], linear with respect to the number of significant wavelet coefficients, i.e. great that some small ǫ. The key-point at that time is the computation of the (explicit) nonlinear term, which for the moment is computed on the regular grid points. For (4.8), the Leray-Hopf projection also reduces to Laplace-Dirichlet problem, in wavelet basis: as previously, this part can also benefit from adaptive schemes developed by [START_REF] Cohen | Numerical Analysis of Wavelet Methods[END_REF]. The only points to be considered are first the operator between the divergence-free wavelet coefficients [d div,n j,k ] and the "standard" wavelet one [d 1,n j,k ] and [d 2,n j,k ], which is diagonal (same indices), see [START_REF] Kadri-Harouna | Effective construction of divergence-free wavelets on the square[END_REF]. The second point concerns the boundary: to avoid a loss of accuracy, even in adaptive scheme, the same number of edge wavelet functions has to be kept at each level (or scale index). This point is well documented in reference [START_REF] Cohen | Wavelet methods for second order elliptic problems -preconditioning and adaptivity[END_REF].

y v 1 , j = 7 v 1 , j = 8 v 1 , j = 9 v 1 , Ref. [3] v 1 ,Ref. [4] v 1 ,
x v 2 , j = 7 v 2 , j = 8 v 2 , j = 9 v 2 , Ref[3] v 2 , Ref[4] v 2 ,

Conclusion.

In this paper we have constructed a divergence-free wavelet based projection method for the numerical resolution of Stokes and Navier-Stokes equations with Dirichlet boundary conditions. Stability and consistency of the method were proved for the unsteady Stokes equations and verified numerically for the Navier-Stokes equa- tions, with application to the benchmark problem of the lid-driven cavity flow, using regular grids.

The main interest of the method lies in the correction step where the Dirichlet boundary condition on the divergence-free solution is exactly satisfied, without imposing boundary conditions on the pressure field, which is not the case in classical approaches. An important issue to even more reduce the computational complexity, is to develop an adaptive method, taking advantage of the compression property provided by wavelet bases, which, for elliptic problems, leads to optimal numerical solvers [START_REF] Cohen | Numerical Analysis of Wavelet Methods[END_REF]. An extension to dimension three of the present method is on-going: a 3D divergence-free and curl-free wavelet based Helmholtz-Hodge decomposition on the hypercube already exists [START_REF] Kadri-Harouna | Helmholtz-Hodge Decomposition on [0, 1] d by Divergence-free and Curl-free Wavelets, Curves and Surfaces[END_REF]. Its incorporation in the correction step of a modified

Theorem 3 . 1 .

 31 The modified projection method (3.7) and (3.8) for the Stokes equations is unconditionally stable.

. 6 )

 6 Substituting this change of variable in a Crank-Nicolson time scheme, (4.1) and (4.3) are replaced by the following new prediction and correction steps for the intermediate velocity ṽn+1 and pressure computation:• Prediction step:

Figure 4 . 1 .

 41 Figure 4.1. Spatial errors at grid points according to δx (log-log scale), scheme (4.7) and (4.8) for the simulation final time T = 1.

Figure 4 . 2 .

 42 Figure 4.2. Steady state horizontal velocity v 1 (left) and vertical velocity v 2 (right) profiles in the middle of the cavity. Solid line (present work) and circle(Bruneau and Saad [4]): Re = 1000 and j = 7.

Figure 4 . 3 .

 43 Figure 4.3. Vorticity contour (left) and divergence-free scaling function coefficients contour (right). Steady state for Re = 1000 and j = 7.

Figure 4 . 4 .

 44 Figure 4.4. Evolution in time of the ratio of divergence-free wavelet coefficients up to a fixed ǫ. The maximal space resolution is j = 8, the coarse decomposition level is j min=3 and Re = 1000.

Figure 4 . 5 .

 45 Figure 4.5. Contour of divergence-free scaling function coefficients at j = 8 (left) and divergence-free wavelet coefficients isovalues (right) for a coarse decomposition level j min = 3. The simulation final time is T = 80 and Re = 10000.

  ǫ = 10 -7 .

Figure 4 . 6 .

 46 Figure 4.6. Divergence-free wavelet coefficients with absolute values greater than epsilon for a maximal space resolution j = 8 and a coarse decomposition level j min = 3. The simulation final time is T = 80 and Re = 10000.

Table 3 . 1

 31 , in terms of the discretization time step δt. As expected (theorem 3.2), the time discretization convergence rate is close -error 30.09E -6 71.72E -7 14.48E -7 16.91E -8 2.0031 Time discretization ℓ 2 -error according to the time step δt, for final time T = 2.

			Backward-Euler
	δt	0.1	0.05	0.025	0.0125	Order
	L 2 -error 25.24E -3 12.70E -3 63.71E -4 31.90E -4 0.9922
			Crank-Nicholson
	δt	0.1	0.05	0.025	0.0125	Order
	L 2

  3.2 saturates to the number of vanishing moments -error 32.11E -7 80.23E -8 20.03E -8 49.79E -9 2.026 L 2 -error 19.84E -7 49.18E -8 12.26E -8 30.49E -9 2.023 H -error 32.11E -7 80.27E -8 20.06E -8 50.17E -9 2.000 L 2 -error 19.84E -7 49.21E -8 12.29E -8 30.72E -9 2.020 H 1 -error 13.38E -5 58.00E -6 26.54E -6 14.05E -6 1.088

				Backward-Euler
	j	6	7	8	9	Order
	L ∞ 1 -error	13.38E -5 57.99E -6 26.49E -6 15.51E -6 1.045
				Crank-Nicholson
	j	6	7	8	9	Order
	L ∞

Table 3 . 2

 32 Spatial discretization errors at grid points according to the resolution j, for final time T = 2.

Table 4 . 1

 41 Velocity v 1 values in the middle of the cavity at the steady state for Re = 1000.

	Ref. [18]

Table 4 . 2

 42 Velocity v 2 values in the middle of the cavity at the steady state for Re = 1000.

	Ref[18]

3D projection method will be the subject of a forthcoming paper.