Jérôme Crevatin
email: jerome.crevatin@eleve.enac.fr

Garrett Smith
email: garrett.smith@cosmicaspacelines.com

INTEGRATION OF SUBORBITAL SPACEFLIGHT AT A COMMERCIAL AIRPORT IN FRANCE

This paper investigates a preliminary evaluation of the suitability of airports in France for reusable suborbital spaceflight operations. The authors propose a methodology and tool to easily perform preliminary spaceport site evaluation, which considers runway length & surface, local weather conditions, existing air & ground traffic, and safety for the uninvolved public. All of these have consequences on site selection for a spaceport.

The authors briefly present the potential legal framework for suborbital flights in France with the objective of finding an acceptable regulatory framework. Due to the absence of explicit requirements, the U.S. regulations established by the FAA-AST have been used as a guideline to frame the scope of the investigations.

A methodology has been developed and a weighting established for various criteria to rate any selected airport. Matlab and C scripts were authored to automate data intensive aspects of the analysis. A simplified aerodynamic model for two representative space planes were developed to provide trajectory data for nominal and nonnominal flight. Data tables were created through manual entry and automated databases extraction to provide the input for the overall analysis.

This paper shows a case study of the airport in Montpellier, France, with a flight path directly over the Mediterranean sea. The overall study shows that several airports in France are acceptable and justify further investigations regarding their feasibility.

I. INTRODUCTION

Commercial suborbital spaceflight is rapidly developing in the United States, with initial operations foreseen within the next few years. After a successful introduction in the U.S. operators will be looking to expand abroad. Unfortunately, commercial operations in Europe, a promising market, might be hindered due to the lack of technical and regulatory preparations in this field in Europe.

France has always been a major actor in space history. The French national space agency, CNES, is recognized around the world for launching the ATV to the ISS on Ariane 5 for example, and the French military (DGA) for launching missiles. Before the selection of Kourou, the CNES studied the south-easterly facing Mediterranean coastline in the Languedoc region to build their launch site. These examples underline the fact that France has always been present for space activities for operations and innovation.

Suborbital spaceflight has many special characteristics: very high vertical speeds during the climb, potentially explosive propellants, glide approach and unpowered landing. All of these make suborbital flights very different from traditional aviation operations and often lead to conflicting needs for airport utilisation.

The authors believe that these conflicts can be appropriately managed and thus are investigating existing airports in France to identify preferred locations. Use of an existing airport also has the benefit of reducing costs associated with runway and terminal development for a green field project.

This study was originally conducted during an engineering internship for the main author, under the technical guidance of the co-author, and hosted by the Destination Etoiles foundation with financing by DLA Piper. The overall ranking of the French airports remains confidential and proprietary to the internship hosts although the aggregate analysis in this paper gives some insight into the promising nature of several airports.

II. SCOPE OF WORK

A very pragmatic approach was taken to limit the scope of research to the main areas of interest for integrating suborbital spaceflight into airport operations in order to ensure the work could be completed within the timeframe of the internship. Due to the very broad domains of potential investigation, the authors chose to concentrate on the essential technical characteristics as illustrated in Figure 1. The regulatory environment served as a boundary for the investigation by providing key parameters to satisfy.

Environmental impacts such as pollution, noise or wildlife were ignored for this study even though the authors acknowledge that they are an essential element in the choice of an appropriate spaceport.

The key objective was to develop a tool to support the application of the method. After choosing a space plane, this software tool is able to rank the airports for suborbital flight operations. The tool could be adapted to rank airports in other countries with additional data.

III. REGULATORY FRAMEWORK

In Europe, there is no harmonized regulation similar to FAA-AST (Federal Aviation Administration, Office of Commercial Space Transportation) regulations for suborbital spaceflight. Smith argues that harmonization of activities is not foreseen in Europe for spaceflight operations due to the Treaty on the Functioning of the European Union and that even if suborbital operations are considered aviation activities, then this remains a national competence. [1] Thus the authors briefly reviewed the French legal situation. The key question to ask is whether this is a space or an aviation activity to identify which French authority to approach. Unfortunately, this question has not yet been answered so both potential regulatory frameworks are briefly presented including their drawbacks.

DGAC Sensational Flights Law

In France, a decree has been adopted concerning the "Vols à Sensations" (thrill seeking flights). The ministry decree dated on the 8 th February of 2012, has been adopted to set a legal framework on the use of the "Airbus ZeroG", operating parabolic flights to reproduce microgravity environment during few seconds and for high performance jet trainer aircraft. This legal framework could be applied to suborbital flights. Unfortunately, the law states that a "Titre de Navigabilité" (title of navigability) must be provided to the DGAC for the flight operations application. Unfortunately, it appears that an FAA-AST launch license does not qualify as a title of navigability. [START_REF] Desvallees | [END_REF] The law text states that the flight must be operated from a point A to A, and applies to any aircraft operations which are publicly advertised. A safety analysis must be conducted, concerning people on the ground and other aircraft in case of non-nominal events.

French Space Operations Law

This "Loi des Operations Spatial" (LOS) n° 2008-518 dated 3 June 2008 is used for launch vehicle and satellite operations, notably for Ariane 5, Soyuz and Vega launches from Kourou in French Guiana. The competent authority in France is the CNES.

This regulation provides a clear legal framework but there are no explicit requirements which operators need to meet. Every operation must prepare an application including safety management and CNES will review it on a case-bycase basis. This law provides requirements for insurance guarantees and the split of responsibility between the operator and the state. It also ensures the safety of the uninvolved public and compliance with international obligations.

Space planes can also be used as launchers and must respect in this case the Space Operations Law. Nevertheless, human spaceflights are specifically excluded from this regulation, which currently

Developments in International Cooperation

The FAA and CNES have signed a Memorandum of Cooperation (MOC) on 16 June 2015 to cooperate on research and development related to the safety of private sector orbital space launches and re-entry activities. Unfortunately suborbital spaceflight is explicitly excluded from the scope of the memorandum. Fortunately for suborbital operators interested in France, the joint press release issued by the FAA and CNES stated that the MOC "...could pave the way for development of future arrangements between the FAA and France related to the regulation of suborbital reusable vehicles...". [3]

IV. DEVELOPMENT OF THE METHOD

The study focuses on technical aspects only, meaning that economical purposes such as market research is not evaluated in the method. Three major factors are analyzed for airport evaluation as shown in Figure 2. First of all, the technical characteristics of the airport must respect spaceplane operational constraints: the runway length or its type of surface. Moreover, an analysis has to be conducted to ensure the safety of people on the ground in case of incident. That part takes into account the population density around the airport. The last criterion but not the least, is the integration in the current air and ground traffic. That part ensures the best cohabitation between air and space operations.

Additional factors have to be studied. The airport safety is a relevant element of the method. It analyzes and ensures the safety of travellers in commercial terminals of the airport but also airport employees and people living in the airport's neighbourhood. The weather is also a fundamental criterion that must not be under-estimated. If the sky is cloudy all year, suborbital activities could not be operated as frequently as operators would like.

All these criteria and how they are analyzed are detailed in following parts.

How to rank airports ?

In addition to evaluating the possibility to use an airport or not for suborbital activity, the method has to rank the best-selected airports. Each criterion will be graded and this grade has a weight.

Then, the weighted average is computed to make the overall ranking. To be included in the ranking, airports should respect minimum requirements. For example, if any trajectory is found to ensure the safety of habitants on the ground, the airport is eliminated. In the opposite, a too-short runway for operational use could be ranked, but receives a low grade due to the construction needed to extend the runway.

There are 4 sub-grades, one for each category. The coefficients and weights for each grade have been defined to represent the importance of each factor and the consequences. But this is a subjective evaluation with the aim to make a ranking. If the analysis is made with the same coefficients and weights, the result is credible due to the tendency.

Method's goals

The method is a "pre-selection". For this reason, the approach is global and may be studied deeper in each parts for an airport-focused study. The method has been reviewed by experts in their domains who approve the approach while acknowledging the shortcomings requiring further study. The final objective is to show that suborbital flights have the potential to be operated in France with all the safety constraints without building a dedicated spaceport.

V. TRAJECTORY DETERMINATION

The first step is to define the nominal path for the space plane. Trajectory is calculated from flight equations. The goal is to compute all the information that will be needed for the following parts of the study. Along the trajectory, flight data

Computing the trajectory

The flight trajectory is firstly defined in two dimensions (Figure 3), computed in the plane: altitude versus downrange. The flight could be divided in five phases, detailed below. For each phases of the flight, aerodynamic equations are written. The lift equation [1] and the thrust equation [2] deal with aerodynamic purposes. Mass equation is added [3] and both position equations [4 and 5]. All equations issued are extracted from "Flight Mechanics" course given by G. Hattenberger (ENAC) and "Spaceflight Dynamics" course book written by William E. Wiesel [START_REF] Hattenberger | Flight Mechanics[END_REF] [START_REF] Wiesel | [END_REF] . Some assumptions are made for this computation. The specific impulse is considered constant in the simplified approach. Other assumptions are made: the angle of attack is assumed constant and defined for each phase of the flight, a friction force is added for the take-off to model the interaction between the runway and wheels.

Equations :

[1]
where with : -m is the mass of the spacecraft, -V its speed, -g is the gravitational constant (g=9,81 m.s -2 ), -γ is the climb angle, -S is the lifting area, -ρ is the air density, function of altitude, -C L is the lift coefficient, [2] where :

where : -Isp is the specific impulse, -C D is the drag coefficient,

[3]
where β is the constant mass flow.

[4 and 5]

To compute this differential system, external data are needed like air density, C L and C D , which is function of z, the altitude. To calculate air density, the standard atmosphere table is used. The software used to compute the trajectory is Matlab.

The different phases along the flight are detailed below: 0 -Take-off (t 0 -> t 1 ) During take off, the contact of the wheels in the ground adds an extra force on the spacecraft. This is modelled by a friction force, which is proportionally inverted to the speed.

1 -1 st segment of climb (t 1 -> t 2 ) During the first phase, the climb angle is constant and equals to γ 1 . The value of this climb angle depends on the aircraft. It is an input data. In this phase, the space plane accelerates. The climb angle is still constant and equals this time to γ 2 (with γ 2 higher than γ 1 ). The value of this climb angle depends also of the aircraft. The plane continues to accelerate until the end of the combustion.

3 -Ballistic shot (t b -> t g ) At t b , the engines power off. At this altitude, there is effectively no atmospheric drag. The trajectory is a ballistic path. The space plane continues the climb until reaching the apogee. Then, the descent starts and the spacecraft accelerates freely. With the reentry in atmosphere, the space plane slows down to reach the descent angle.

4 -Glided phase (t g -> t end ) At t g , the plane continues its descent with a fixed constant descent angle γ d . The plane decelerates.

The trajectory is computed to a point near the airport at a certain altitude. Then, the plane circles to the airport. Some space planes are air launched from an airplane. In this case, the take-off phase is not computed and the flight is analyzed from the dropped altitude with the first segment of climb.

Data computed with the previous explanations create a straight line in the downrange direction, varying only in altitude. Then, a 180-degrees turn is added to make the trajectory go back to the airport. In the method, one turn is allowed after the take-off to avoid populated area. This turn is automatically computed in the following part (non nominal study) to find the best trajectory as a function of population density. The trajectory is now in three dimensions.

Converting of the trajectory

The trajectory has been defined by solving differential equations. The output is a list of points, in a local frame. The points defining the trajectory must be converted in a global coordinate system. The LLA (Latitude, Longitude, Altitude) system is chosen with the reference WGS84, used for GPS systems.

For a dropped spacecraft, the dropping airplane can reach a precise and desired position but the operator cannot choose a point by chance. In an operational point of view and for safety reasons, the operator must take into account two important things: in case of ignition failure just after the drop, the space plane must be able to reach an airport, gliding, and the second one is that the space plane must be able to come back at the launching airport after apogee. The first condition is not analyzed in the tool but it is part of the method and must be checked manually. The second one is checked by the tool. The trajectory is oriented to the airport when re-entry has been made. Then the altitude over the airport is checked if it is high enough or not to accept the trajectory.

Results are stored in excel files. At this time, data can be used to visualize the trajectory and can be used for the following studies: non-nominal flight and traffic integration.

Visualizing of the trajectory

The idea is to visualize clearly the trajectory in GIS software such as Google Earth. A file readable by this kind of software (extension .kml) is written with the Matlab code. The trajectory can be seen in the software, starting at the airport to the apogee, before coming back to the airport (Figure 6).

Re-entry

The re-entry is a critical part of the flight path. Spaceplane behaviours during this phase are not really easy to predict. First of all, the theory from the re-entry aerodynamics is needed to compute the flight path. Then, the trajectory is drawn and defined in function of the performances parameters of the space plane (speed, angle of attack) and also the position around its gravity centre.

The re-entry must be taken into account for the non-nominal part because it is a critical part of the flight, but for the trajectory, the impacted elliptical area defined in the next part already takes into account a non nominal event during this phase.

VI. NON-NOMINAL STUDY

The purpose of this part is to ensure the safety of people on the ground. First, failure probabilities are defined using existing data on space and air operations to evaluate the acceptable risk level. Then, the computation of debris impact is done to be compared with population density.

Analysis description

First, the probability of having a power loss is defined. The value is around 1%, uniformly distributed on combustion time.

In reality, this probability is little higher at the beginning of the flight, after ignition but to simplify calculations, it is ignored. The plane glides back to the airport or proceeds to an emergency landing in a non-populated area minimizing injuries for people on board and on the ground if the plane remains usable.

Two different cases appear when the plane is unusable: dislocation or explosion. The worst case is the explosion due to the energy released during the explosion giving relative speeds to debris causing a high dispersion of debris on the ground. Both cases are considered significant for the safety ground study and plane dislocation is included in the explosion case. This is conservative. Unusable plane case appears every 1.000 launches, regarding Figure 7. P flight is equals to 10 -3 . [6] Moreover, the probability P 2 is the one to kill at least one person knowing that a non-nominal event happened at instant t. It has to be computed all along the propelled phase. It is characterized by the Casualty Expectation (E C ) [7] :

[6]
where : -P(i,j) is the probability that the debris i hits the cell j, -N j is the number of habitants in the cell j, -A c,i is the impacted surface of debris i (Casualty Area), -S j is the surface of cell j.

The difficulty of this method is to define P(i,j). This probability is usually found running a Monte-Carlo simulation. For a first approach, it is reasonable to simplify the method aiming to reduce computing time during the analysis. For a second run, the Monte-Carlo has to be chosen.

Simplified method

The simplification deletes the Monte-Carlo simulation. At an instant t, the explosion is simulated. The energy delivered by the explosion is function of the remaining fuel mass. In fact, the energy "E" is proportional to the TNT equivalent mass (given in kilograms) of the fuel and propellants remaining at the moment of the explosion defined in Joules in the Equation [7]. To compute the equivalent TNT mass, the chart given by the FAA-AST [7] is used. Many assumptions are made. The energy received by all debris is the same. The explosion energy is not fully converted to kinetic energy. Only 50% of the energy is transmitted as a blast, remaining energy for thermal and radiation [8]. But the debris doesn't go as fast as the blast. In the end, it is assumed that 10% of the initial energy is transmitted to the debris. Then the relative speed for a debris can be found with equation [8].

[8]

where : -E k : kinetic energy for one debris which is equals to the total kinetic energy divided by the number of debris -m i : mass of debris i

To compute debris trajectories, only few extreme debris are analyzed to increase computing performances. Debris are modelled around the space plane gravity centre. The computation uses the same equations described in part IV for aerodynamics and impacts on the ground can be found. Two assumptions are made: debris do not bounce on the ground and people protection due to buildings is not taken into account. An elliptic zone is modelled around the impact points and converted in the LLA system. This elliptic area is named S t . The sum of all the elliptic areas S t gives the S area, which can be visualized in Google Earth (Figure 8). The Casualty Area (A C ) is estimated for each space plane. This is the total surface covered by all the debris increased by 36 centimetres. It is really difficult to estimate this value which depends on the number of pieces constituting the plane. The fuselage surface is estimated with available data and wings surfaces are known. Inside pieces has to be added, such as seats, computers, tanks. But the Casualty Area depends also with the numbers of debris.

Knowing the LLA coordinates of S, the number of people living in the area can be found. The national database is used. The mesh is composed of 200-meters-side-squares. The coordinate system used is the LAEA Europe -ESPG3035 based on ETRS89. It is a (x,y) system, using meters. The coordinates of surface S have to be converted in this system and the number N of habitants in the area can be found.

The number to evaluate is equal to (N*A C )/S (equation [9]). This is the number of person reached (N R (t)) by debris, for one simulation (at a given instant t).

[9]

A limit must not be crossed. A test is realized at each instant t along the propelled phase. This limit must respect the safety probability for space, characterised by E C,MAX [7] : E C,MAX = 30.10 -6 .

Previously, it has been shown that the probability that a catastrophic event happens (P flight ) is equals to 10 -3 . E C can be written, in equation [10]: [10] This equation gives us a maximum of = 30.10 -3 for the average of N R. To know if a trajectory is acceptable or not, the average of N R (

) must be computed for the whole combustion phase. The value can be so high that it is a waste of time to simulate during all the combustion phase. A test is made to gain time. The mean is computed for each instant t (= ), assuming that during the whole remaining combustion time, nobody is reached by debris.

If

< , the test is valid and the study can begin for the next time step. If not, this trajectory cannot be accepted. This process is explained in the following part.

Trajectory Optimization

If the risk exceeds the limit, the initial trajectory must be defined again including turns to avoid the highly-populated area. If all the solutions had been tested and the population criteria never reached, the airport is not usable for suborbital flights. The process tests at each instant t if a non-nominal event The program computes in this case a new trajectory including a turn before the instant t where the probability was too high. The default value is 20 degrees. The program computes alternatively a new trajectory to the right and then to the left of the original trajectory. If a complete U-turn has been realized and the test is still accepted, the airport is eliminated. Analyzed data are stored in a file to show the number of habitants, the surface impacted and the probability, for each time step.

Operational point of view

The non-nominal study analyzes both take-off directions for departure. It means that two nominal trajectories have to be found for horizontal T/O planes. If a trajectory is accepted for one take-off direction only, then the airport should not be penalized due to possible downwind T/O.

Re-entry

The re-entry is one of the most critical points during a suborbital flight. The consequences are high temperatures and pressures on the plane. Suborbital planes re-entry is quite different due to the relatively low speed of the vehicle around Mach 4 ; it was Mach 25 for the Space Shuttle.

Anyway, the re-entry should be taken into account in the P flight . For safety reason, the re-entry has consequences on plane flying above the space plane if an explosion or dislocation occurs. The process is detailed in the next part.

VII. TRAFFIC INTEGRATION

This part deals with the integration of suborbital flights in current air traffic. Vertical and high velocity manoeuvres of those spacecraft make the integration difficult. The actual way to manage the ATC is more horizontal and the air traffic controller is not used with vertical path around Mach 3. A process must be established to ensure safety during the spaceplane's ascent and descent.

VII.I Operational issues

Most of space planes have emptied their tanks before reaching the apogee and will not be able to perform a go around in case of missed approach or, in case of extraordinary closure of the airfield. Some spacecrafts have a special engine system which can be re-ignited and let the plane perform one go-around to integrate again the circuit for a last landing. This imposes one statement: the runway must be closed from the moment that the space plane takes-off until its landing has been completed. The commercial traffic is impacted during the runway closure.

Alternatives can be considered:

-Two (or more) parallel runways in the airport. One runway can remain opened. But, the probability that the two runways may be closed at the same time has to be analyzed.

-An alternate airport could be chosen. It has to be clearly defined for each airfield. From the take-off until a decision point, the alternate airport is closed, or at least one runway is closed. At the decision point, the pilot of the space plane chooses the regular or the alternate airport. The one chosen is closed, or at least one runway, until the landing completed. The other one can be opened.

-For dropped space plane, the closing time of the "rescue airport" and "final airport" runways is the dropping time.

The method has to analyze the impact of a runway closure. Airports with a low number of movements per day are the most appropriate to operate suborbital spaceflights.

VII.II Avoiding en-route airplanes during climbing and gliding phases

Cruising airplanes follow routes or airways, materialized by virtual points. The major part of commercial planes movements can be known using the airways positions. But some planes don't really follow those airways due to ATC needs to avoid collision or direct-route trajectories. The selection of an airport must take into account the airways above and around the runways. The idea is to find a vertical corridor, wide enough to ensure the safety of operations but impacting as little as possible the air traffic. This corridor must include the fact that a non nominal flight can occur and the debris impact zone has to be taken into account.

A no-cross zone is defined relative to the trajectory. The number of airplanes and the frequency crossing this zone is analyzed. A database composed of all flights routes for one typical day has been created with the help of different CRNA (French Air Traffic Control Centres).

The zone envelope takes into account the trajectory of the spacecraft but also the trajectories of debris in case of incident. The zone is activated during all the propelled time and remains closed until the re-entry has been completely done, at least. The number of airplanes crossing the virtual zone and their crossing time can be found and analyzed. If the zone is activated, those planes should modify their route to avoid it under the orders of ATC. The less aircrafts diverted, the best is for suborbital activity for airport selection.

If the ascent phase is complicated due to the high speed of the vehicle, the way back is not so simple. The plane goes slower, but without any propulsion system, meaning that it is restricted regarding altitude changes. When it approaches again the altitude around 15 kilometres, the air traffic control must be in contact with the pilot to manage the situation with other aircrafts.

VII.III Ground traffic impact

Commercial planes and general aviation use also taxiways and runways. As mentioned in part A, the eventuality of the runway closing has to be taken into account. To quantify the impact on commercial operations, the number of aircraft movements can be a relevant indicator [START_REF] Dgac | [END_REF] .

The very restrictive consequences of ground traffic in comparison of air traffic should be shown in the grading. The air traffic can be diverted and the time spent to do this little change on their route is around the minute. But on the ground, when suborbital operations start, planes wait in the air or on the ground but there is not other solutions, except if there is an other runway.

VIII. ADDITIONAL CRITERIA VIII.I Airport characteristics : runway length and strength

This criterion links space planes performances and airport characteristics, especially runway length and strength. The take-off run is evaluated during the trajectory modelling. But one thing has to be considered : an engine failure shortly after take-off. Even if the non nominal study takes into account the population density, this part of the method is more for the safety of people inside the space plane.

Runway length

If the space plane has a total engine failure, it can glide. But it has to be at a certain altitude to come back and align on the runway to perform an unpowered landing. The distance computed takes into account the critical point at the limit where the space plane cannot perform the manoeuvre described just before. From this critical point, the space plane glides to the runway but keeping the take-off axis and can land safely. For a horizontal take-off spaceplane, the needed distance is evaluated at 2600 meters. Under this value, the runway must be extended and the grade is lower for the ranking.

For an air launched space plane, the needed distance is the one necessary for the carrying plane to take-off with the space plane attached on or under it. This distance is estimated around 3000 meters.

Runway strength -ACN/PCN

Runway strength is never a limiting parameter. The relative small weights of the space planes do not impose restrictions on this criterion.

VIII.II Ground safety study

Commercial airplanes burn kerosene when spaceplanes use propellant mixes such as liquid oxygen (LOX) and RP-1 (fuel), or liquid hydrogen for example. This implies that safety distances between public and spacecrafts must be defined and respected. Precise processes during refuelling or in case of LOX leak have to be established.

Safety distances

Three things are considered if an explosion happens: the pressure wave, the heat wave and the debris dispersion. The distance needed is directly taken from the regulation, proportional to the NEW (Net Equivalent Weight). The American legislation defined by the FAA-AST was taken into account for the study. To analyze if an airport is ready to welcome such machines, distances between any point where the space plane is susceptible to be and commercial terminals, roads or maintenance and servicing buildings must be analyzed. This study requires a huge database for airport buildings. The choice that has been made is to draw the exclusion zone. Then, it is possible to visualize it on GoogleEarth or GeoPortail with the building database.

Special processes

Due to the fact that space planes could evolve near commercial aircraft, some special processes must be defined. In the exclusion zone, it is necessary that nobody except the tow pilot is present. Lighters, cigarettes, telephones, electricity systems, or any objects which could provoke spark are excluded from the exclusion zone. During refuelling, the operation must be done over concrete even if the airport pavement is made of asphalt. A special concrete area must be built. During taxiing, the doors are unlocked and "passengers" unattached until the power on.

One of the most critical event and not usual for commercial aviation is a liquid oxygen (LOX) spill. According to the studies [10][11] , one conducted by the NASA, the cohabitation of LOX and asphalt could be explosive in the opposite of the concrete and LOX contact where nothing happen. A surveillance system for LOX spill must be built. For this reason, non-concrete runways are penalized in the grading for LOX-use planes.

VIII.III Local weather conditions

Space planes operate under VFR conditions, meaning that the pilot uses his eyes to control the plane and is not allowed to go through clouds. In fact, the local weather conditions must be favourable to ensure that operations can be done as many days as possible. The unique constraint is the impossibility to go through the clouds. Meaning that a partially cloudy sky can be reasonable to fly. The rain is not a constraint directly to operate or not the flight, but will increase the needed runway length for landing.

The method analyzes three major weather factors. The first one is the rain, the second is the period of sunshine and then, the wind. For each factor, a mean value is taken into account but also the number of days regarding special events. For example, the first value for sunshine level is the average sunshine in hours per year. The second factor is the number of days where the sunshine is lower than 20%. It means that this day, a suborbital activity is probably not possible. Due to gliding approach and unpowered landing, it is better to have also lower winds and the most regular as possible, meaning that gusts are not appreciated.

VIII.III USING THE TOOL AND ANALYZING RESULTS

Airport selection and data

Due to the very selective criterion of nonnominal flight density limit, it was not necessary to study all the airports in France. The position of each thresholds and runway orientation had been stored in the data file using Google Earth tools for selected airports. The runway length and width have been taken from the SIA (Aeronautical Information Service). Airports too far from the wanted criterion have been directly eliminated. Airports with long runways, near a low-populated area and a lowtraffic density have been chosen where the weather conditions fitted well with suborbital operations. The following airports have been evaluated : Some airports could appear far from criteria chosen (Nice for traffic for example), but other criteria could compensate the one which is not good for the ranking. This is for this reason that Nice or Marseille have been chosen for example.

If the user wants to study another airport, he just has to enter airport characteristics in the data file and select the airport when he uses the tool.

Traffic data

Traffic data were essential for the project. The 5 different French en-route control centres provide the data from their databases. Planes are identified by radars of the control centre. Many information is given in the database : position, time, horizontal and vertical speeds, altitude, Flight Level requested, heading, departure and arrival airports, call sign and aircraft type. Only the time information is very useful for the study.

Population density data

For the non-nominal study, the population density was used. To be precise, a mesh built with squares of 200 meters-side-length is used by the tool. The data comes from the INSEE (National Institute for Statistics and Economic Studies).

Same process can be done with other databases of world's countries.

IX. TOOL RUNNING

The tool developed under Matlab code has two big different uses. The first one is an execution of the method, the user just have to select a space plane and the tool will run to rank best airports for this space plane. For the second aspect, the tool is used to check a specific criterion. The user can choose an airport and a space plane, then he has the possibility to draw the trajectory, to find a trajectory respecting the non nominal criteria or checking the probabilities with a defined trajectory, to analyze traffic flow in the airport area, technical characteristics or weather conditions.

Output files are created when the method is used. Two different types of files are made. Excel files are created for stored data : flight parameters, non nominal detailed information and traffic data representing crossing airplanes.

The other kind of files is the .kml one. These files can be opened with a GIS software and specially Google Earth to visualize the results of the study. Three different things can be shown in Google Earth. The first one is the trajectory for both QFUs, then the S area can be visualized in green and the last thing is for the ground safety study, the no-cross zone is drawn around the runway.

X. CASE STUDY -Montpellier's Airport

The aim of the case study is to show how the method runs through the tool developed. It can show the limits of the method.

X.I Airport global situation

Montpellier is located in the Mediterranean coast. This is a touristic area with warm and sunny weather. For space, it is a dynamic region where space school and previous studies on possible space launch site building took place. The choice of Montpellier has been made because it is a quite good place to operate suborbital flights. Furthermore, Montpellier was already studied for suborbital operations by the co-author and other members of the 3AF (Association Aéronautique & Astronautique de France).

X.II Checking airport characteristics

Each airplane has its own performances for take-off and landing. This preliminary study just checks if the needed distance is available or not at the airport.

-For the horizontal take-off planes, the needed distance is 2600 m (blue distance on Figure 24). -For air-dropped space planes (in fact their carryplanes), the needed distance is 3000 m. (red distance on Figure 10). At Montpellier, the available distance is 2600m. That gives us a grade of 100% for the horizontal take-off spaceplane and 73,4% for air-launched spaceplane. In addition, the type of runway is significant. For Montpellier, it is an asphalt runway. The coefficient for the grade will be 0,75 for LOXuse on an asphalt runway versus 1 for a concrete runway. This coefficient gives the grade of 75% for horizontal-take-off space planes and remains 73,4% for air-dropped ones.

X.III Finding a trajectory

The first task to perform is to find a nominal trajectory that will be accepted in case of non nominal event for each kind of space plane. Regarding horizontal take-off space plane and before running the simulation, the QFU 12L is logically the best one to use in term of ground safety, due to the presence of the sea. The tool helps to determine if the maximum average N R for safety ground will be reached or not. In the other direction, there is the city of Montpellier. It seems to be difficult to find a trajectory acceptable. For dropped-space plane, the coast is close enough to the runway to drop the plane above the sea, reaching anybody in case of non-nominal event. The dropping point is close enough if a failure ignition happens and then the space plane can reach Montpellier, gliding.

The QFU 12L is first analyzed. To start with, a trajectory is computed right in the direction of the runway, no turn added. Rapidly, the tool performs iterative turns to avoid the dense area of la Graudu-Roi. This trajectory is acceptable for QFU 12L (Figure 11). For the QFU 30R, the idea was to come back as early as possible to fly over the sea. A 230° turn has been created to the right, to avoid Montpellier. The city of Mauguio was well avoided but then some cities remained in the S area. Nevertheless, with the same process as before, the average N R was computed and accepted at the limit for this trajectory.

X.IV Analyzing air & ground traffic

Now, trajectories are defined. Traffic analysis can start. A volume is created around the space plane to include debris dispersion in case of explosion. This volume is 19 km wide, 100 km long and has a height of 100 km. 112 planes had been found in this area for one day. No slot of 25 minutes has been found too. The Air traffic airport grade is 60%. For the ground traffic, there are 15 557 movements per year (for the year 2014). The ground traffic grade is 90%. The global traffic grade is 75%.

X.V Weather conditions

For the wind study, there are 3 major directions for the wind : both in the runway axis and the last one is in the North-North-East direction, which is in crosswind direction. This direction is the major one from November to February, during 32% of the time up to 17 knots. The grade of 49,8 % is obtained. For the sunshine level, there are 2668 hours of sunshine per year (national average of 1664 hours/year), equivalent of 222 sunny days. The sunshine grade is 87,4 %. The last criterion is rain, and the average mean is about 578 mm/year against 652 mm/year for the French mean. It must be underlined that huge thunderstorms increase this figures. (cf. flooding in Hérault region). The rain grade is equals to 83,4 %. The global weather grade is 80,6 %, which ranks the airport in 6 th position for the weather conditions.

X.VI Ground safety

The tool draws the exclusion zone in Google Earth. Then, the analysis must be done manually. The exclusion zone is shown in Figure 12. On the threshold 12L (extreme left circle) the distance is not respected with the road and the house.

The road must be closed when the space plane is on the extremity of the runway and habitants have to be considered protected inside the house.Trails located in the North-East of the runway should be closed. Except that, everything is fine for Montpellier ground safety study. Planes must not be on the taxiway when suborbital operations start.

X.VII Conclusion

The Montpellier's airport is a quite good place for suborbital operations. The traffic is not so dense and particularly on the ground. The weather conditions are excellent with a warm and sunny weather, but sometimes windy. The runway is long enough for horizontal take-off space planes but should be increased by 400 metres for dropped ones. The ground safety is not so bad, but can be improved by doing small works. Montpellier is quite ready to welcome suborbital operations.

XI. LIMITATIONS AND IMPROVEMENTS

As mentioned in the beginning of this paper, the method and the tool developed for it have been made to do a pre-selection of airports. It is a global approach taking into account many criteria. The first limitation, if it is one, is the fact that it is not a study which can be used for operational approvals. Nevertheless, some suborbital consequences have too be studied also : propellants storage, noise study, pollution forecast for example.

Method

Lacking criteria should be integrated in the next version of the study. Environmental aspects are very relevant nowadays and they are burning issues for the society.

The Non-Nominal part has to be more precise. This is the part that should be improved the most. Statistical approach must be used for debris dispersion. The computation method should be reviewed too: the fact that people are inside a building or not is not taken into account.

The second big work is to improve aerodynamic model. The idea is to create directly the trajectory in 3D. It allows multi-turn trajectory and the turn consequences (decrease of apogee altitude) can be quantified. The re-entry part of the trajectory should be studied deeper and integrated in the tool.

Weather part can be improved with more precise data. The idea is to determine number of slots usable to operate per day.

Tool

The tool is operational and works well. The functions defined at the beginning of the internship have been created and compute well things established as objectives. The main problem of the tool is its speed performances, it is very slow. To solve the problem, many parts of the tool (maybe all) should be written in an other language (C, Fortran). Nevertheless, the huge databases used for computations and the need to have an interface with the user lead to the fact that using an other language could be the best idea.

XII. DISCUSSION

The goal was to prove that suborbital operations can be suited in France. Of course, it is easier to build a dedicated platform in the middle of the desert to welcome suborbital activities and if something happens, there is no risk for reaching people on the ground. Does it mean that a space plane will never take-off from France ? From Europe ? No. Europeans are late but the things are moving and quite rapidly, regarding the UK for example. French administration is known to be slow, but "Destination Etoiles" and "Cosmica Spacelines" are stepping on the gas. Some airports are well-designed to fit with suborbital activities. The runways are long enough and the traffic not so dense. Suborbital flights can boost the airport development and attract people in the area.

To respect the confidentiality of results, airports are designated by a random letter. Following rankings have been made. But suborbital flight cannot just be a gift from the sky. Investments must be done to secure the operations, the users and workers at the airport, the population under the flight path and all the passengers inside commercial aircrafts. The perfect airport does not exist at this time, no one had reached the 100% grade.

Airport

The runway length and type is the easier thing to adapt, most of the time it is just a problem of money. The traffic integration can be solve with the idea proposed just above : shared-activity between airports.

The non-nominal consequences and ground safety are more difficult to make them evolve. This part deals with the regulation evolution : possibility to reduce exclusion zone and new way to compute the risks with a more acceptable criterion. The other solution is to make people move from critical area, which is, I think, quite impossible and more difficult in a country as France.

Other criteria not studied in this report are very relevant and must be analyzed. First of all, there is pollution. People are very sensible to environmental issues and the use of LOX or other propellants should be quantified deeply. What are the consequences if a ground LOX (or anything else toxic) tanks breaks into parts? What is the quantity of pollution during initial climb? The noise is the second big aspect. Will local population accept it? Some places are ready near the Mediterranean Sea to welcome suborbital flights. Technically, it is feasible. Soon, the authors are sure that an individual will reach the edge of space from Europe, making him think that space is so close ; 100 kilometres, it is nothing at the Earth scale. And the proximity is underlined if the launch is made near cities and life, not in a desert where you have already left the every-day's world.

XIII. CONCLUSION

Regulation is a big part of the study. An evolving regulation could change the ranking: lower minimum safety distances or higher threshold for non-nominal part. Technology could also bring some changes regarding grades. The space planes may take-off faster and could use smaller runways or climb better and the no-cross zone becomes smaller, impacting less commercial aircraft. The implementation of basic-airplane engine to let the runway open during suborbital operations could be also a solution.

Suborbital activities are just at their beginning and they will rapidly develop themselves and it has been shown that it is possible in France. The South part of France is inevitably the best one with the Mediterranean Sea, a warm and sunny weather and lots of airports in a touristic region.

The future selected airports to operate suborbital flights will be the most valuable to launch orbital flights. Will the South of France will be connected to New-York in less than one hour ?

Figure 1 -

 1 Figure 1 -Scope of study

Figure 2 -

 2 Figure 2 -Parts of the method

Figure 4 -

 4 Figure 4 -Aerodynamic forces (www.astop.com)

Figure 5 -Figure 6 -

 56 Figure 5 -Different frames

Figure 7 -

 7 Figure 7 -Flight failures (CNES, C. Bonnal)

Figure 8 -

 8 Figure 8 -S area visualization and debris impact modeling

Figure 9 -Figure 10 -

 910 Figure 9 -Visualization of population density through the .kml file at Nice

Figure 11 -

 11 Figure 11 -S area at Montpellier

Figure 13 -

 13 Figure 13 -Results for autonomous spaceplanes


Loïc Lacouture

Local error analysis for the Stokes equations with a singular source term

Keywords: Stokes problem, Dirac measure, Stokeslet, finite element method, local error estimates

The solution of the Stokes problem with a punctual force in source term is not H 1 ˆL2 and therefore the approximation by a finite element method is suboptimal. In the case of Poisson problem with a Dirac mass in the right-hand side, an optimal convergence for the Lagrange finite elements has been shown on a subdomain which excludes the singularity in L 2 -norm by Köppl and Wohlmuth, and, independently, Bertoluzza and co-authors have proved a quasi-optimal convergence in H s -norm, for s ě 1. Here we show a quasi-optimal local convergence in H 1 ˆL2norm for a P k {P k´1 -finite element method, k ě 2, and for the P 1 b{P 1 . The error is still analysed on a subdomain which does not contain the singularity. The proof is based on local Arnold and Liu error estimates, a weak version of Aubin-Nitsche duality lemma applied to the Stokes problem and discrete inf-sup conditions. These theoretical results are generalized to a wide class of finite element methods, before being illustrated by numerical simulations.

Introduction.

This paper is about the accuracy of the finite element method to solve the Stokes problem with a punctual force in source term. Let us consider this following problem

$ & % ´△u `∇p " δ x 0 F in Ω, divpuq " 0 in Ω, u " 0 on BΩ, (1)
where Ω Ă R 2 is a square, and δ x 0 F denotes the punctual force F located at x 0 P Ω such that distpx 0 , Ωq ą 0. Our interest in Problem (1) is motivated by the modeling of the movement of thin structures in a viscous fluid, such as flagella connected to bacteria or cilia involved in the muco-ciliary transport in the lung [START_REF] Fulford | Muco-ciliary transport in the lung[END_REF]. Indeed, for instance in the lung, the cilium is very thin and a direct simulation with a graded mesh would be too expensive. In the asymptotics of a zero diameter cilium and an infinite velocity, the cilium is thus replaced by a line Dirac of forces in source term. In order to ease the computations, the line Dirac of forces is approached by a sum of punctual forces distributed along the cilium [START_REF] Lacouture | A numerical method to solve the Stokes problem with a punctual force in source term[END_REF]. Finally, by linearity of the Stokes problem, the analysis of the subsequent problem reduces to Problem [START_REF] Araya | A posteriori error estimates for elliptic problems with Dirac delta source terms[END_REF].

In dimension 2, Problem [START_REF] Araya | A posteriori error estimates for elliptic problems with Dirac delta source terms[END_REF] has no H 1 pΩq 2 ˆL2 pΩq-solution. Although the numerical solution can be defined in this case, the H 1 pΩq-error (respectively L 2 pΩq-error) for the velocity (respectively the pressure) has no sense, and the L 2 -estimates of the velocity cannot be derived like in the regular case without suitable modifications.

Let us notice that the scalar version of this problem, the Poisson Problem with a Dirac mass in right-hand side, has already been widely studied. It occurs in many applications from different areas like in optimal control of elliptic problems with state contraints [START_REF] Casas | Control of an elliptic problem with pointwise state constraints[END_REF] or in the mathematical modeling of electromagnetic fields [START_REF] Jackson | Classical electrodynamics[END_REF]. Problems of this type are found in controllability for elliptic parabolic equations [START_REF] Casas | Parabolic control problems in measure spaces with sparse solutions[END_REF][START_REF] Casas | Spike controls for elliptic and parabolic PDEs[END_REF][START_REF] Leykekhman | Optimal error estimates for finite element discretization of elliptic optimal control problems with finitely many pointwise state constraints[END_REF] and in parameter identification problems with pointwise measurements [START_REF] Rannacher | A priori error estimates for the finite element discretization of elliptic parameter identification problems with pointwise measurements[END_REF]. In this case, Babuska proved in [START_REF] Babuška | Error-bounds for finite element method[END_REF] for two-dimension smooth domain an L 2 pΩq-convergence of order h 1´ε , ε ą 0, where h is the mesh size, and Scott has shown in [START_REF] Scott | Finite element convergence for singular data[END_REF] an a priori error estimates of order 2 ´d 2 , where the dimension d is 2 or 3. Casas has got the same result in [START_REF] Casas | L 2 estimates for the finite element method for the Dirichlet problem with singular data[END_REF] for general Borel measures on the right-hand side.

To the best of our knowledge, there is no finite element method convergence result for the Stokes problem with a punctual force in source term. Moreover, in applications, the punctual force (or the Dirac measure) at the point x 0 is often a model reduction approach and the finite element solution does not need to approximate precisely the exact solution at the point x 0 . Thus, it is interesting to estimate the error on a fixed subdomain which does not contain the singularity. In the case of the Poisson problem, Köppl and Wohlmuth recently showed in [START_REF] Köppl | Optimal a priori error estimates for an elliptic problem with Dirac right-hand side[END_REF] a quasi-optimal local convergence for low order in L 2 -norm for Lagrange finite elements and optimal local convergence for higher orders. A quasi-optimal local convergence in H s -norm, s ě 1 and an optimal local convergence in the case of low order have been proved in dimension 2 in [START_REF] Bertoluzza | Local error estimates of the finite element method for an elliptic problem with a dirac source term[END_REF]. In this paper, we establish still in dimension 2 local error estimates for the Stokes problem with a punctual force in source term, Problem (1), and show a quasi-optimal convergence in H 1 ˆL2 -norm. The proof is based on the Arnold and Liu Theorem [START_REF] Arnold | Local error estimates for finite element discretizations of the Stokes equations[END_REF] that establishes local error estimates for finite element discretizations of the Stokes equations with regular source term. It is written in the case of the P k {P k´1 elements for k ě 2, and the MINI finite element method P 1 b{P 1 if k " 1. No graded meshes are required for these results and they imply that there is no pollution effects far from the singularity.

The paper is organized as follows. Our main result is formulated in Section 2 followed by the Arnold and Liu Theorem [START_REF] Arnold | Local error estimates for finite element discretizations of the Stokes equations[END_REF], an important tool for the proof presented in Section 3. Our theoretical results are generalized in Section 4, before being illustrated in Section 5 by some numerical simulations.

Main results.

In this section, we first set all the notations used in this paper. Then, we formulate our main result and give an important tool for the proof: the Arnold and Liu Theorem. For the sake of clarity, this result is first set and proved in the particular case of the P k {P k´1 finite elements (k ě 2) and the P 1 b{P 1 elements. It will be generalized in Section 4.

Notations.

For a domain D, we will denote by } ¨}s,q,D (respectively | ¨|s,q,D ) the norm (respectively semi-norm) of the Sobolev space W s,q pDq, and by } ¨}s,D (respectively | ¨|s,D ) the norm (respectively semi-norm) of the Sobolev space H s pΩq. Letters in bold refer to a vector of R 2 or a function with values in R 2 .

For the numerical solution, let us introduce a family of quasi-uniform simplical triangulations T h of Ω, where h is the meshsize. For the approximation spaces V k h and W k h , we will use the

P k {P k´1 x 0 Ω Ω 0 Ω 1 BΩ 0 BΩ 1 mesh Figure 1: Domains Ω 0 and Ω 1 .
finite elements, for k ě 2, defined as

V k h " v h P C p Ωq 2 ˇˇv h| T P P k rT s, @T P T h ( , W k h " ! p h P C p Ωq ˇˇp h| T P P k´1 rT s, @T P T h ) ,
and if k " 1, we will use the MINI finite element method P 1 b{P 1 , where P 1 b denotes the continous piecewise linear and bubble functions. For a finite element T , the bubble function b is defined by bpxq "

# λ T 1 pxqλ T 2 pxqλ T 3 pxq if x P T, 0 else,
where λ T 1 , λ T 2 and λ T 3 are the barycentric coordinates of x in relation to the triangle T . We fix two subdomains of Ω, called Ω 0 and Ω 1 , such that Ω 0 ĂĂ Ω 1 ĂĂ Ω and x 0 R Ω 1 (see Figure 1). We consider a mesh which satisfies the following condition: Assumption 1. For some h 0 , we have for all 0 ă h ď h 0 (see Figure 1),

Ω m 0 Ş Ω c 1 " H, where Ω m 0 " ď T PT h T Ş Ω 0 ‰H T.

Statement of our main results.

Our main result is given by Theorem 1. The rest of the paper is mostly concerned by the proof, the generalization and the illustration of this theorem.

Theorem 1. Consider Ω 0 ĂĂ Ω 1 ĂĂ Ω satisfying Assumption 1, k ě 1, 1 ď q ă 2, let pu, pq P W 1,q 0 pΩq ˆLq 0 pΩq be the solution of Problem (1) and pu h , p h q its Galerkin projection onto V k h ˆW k h satisfying ş Ω p h " 0 and ż Ω ∇pu ´uh q :: ∇η ´żΩ pp ´ph qdivpηq " 0 for all η P V k h , ż Ω divpu ´uh qξ " 0 for all ξ P W k h .

(

) 2
Under the assumption that pu, pq P H k`1 pΩ 1 q 2 ˆHk pΩ 1 q, there exists h 1 such that if 0 ă h ď h 1 , we have,

}u ´uh } 1,Ω 0 `}p ´ph } 0,Ω 0 ď CpΩ 0 , Ω 1 , Ωqh k a | ln h|.
2.3 Regularity of the solution pu, pq.

In this subsection, we focus on the singularity of the solution, which is the main difficulty in the study of this kind of problems. In dimension 2, Problem (1) has a unique weak solution pu, pq P W 1,q pΩq 2 ˆLq 0 pΩq for all q P r1, 2r. Indeed, the 2d Stokeslet denoting pu δ , p δ q is defined as (see for instance [START_REF] Pozrikidis | Boundary integral and singularity methods for linearized viscous flow[END_REF])

u δ pxq " 1 4π ˆ´ln }x}I 2 `x ¨tx }x} 2 ˙F, p δ pxq " x ¨F 2π}x} 2 .
(3)

The Stokeslet pu δ , p δ q satisfies in D 1 pR 2 q " ´△u δ `∇p δ " δ 0 F, divpu δ q " 0, so that the Stokeslet pu δ p¨´x 0 q, p δ p¨´x 0 qq contains the singular part of pu, pq, the solution of Problem [START_REF] Araya | A posteriori error estimates for elliptic problems with Dirac delta source terms[END_REF]. As it is done in [START_REF] Araya | A posteriori error estimates for elliptic problems with Dirac delta source terms[END_REF] in the case of the Poisson problem, the solution pu, pq can be built by using a suitable lift procedure which consists in adding to pu δ , p δ q a corrector term pw, πq P H 1 pΩq 2 ˆL2 0 pΩq, which satisfies the following problem:

$ & % ´△w `∇π " 0 in Ω, divpwq " 0 in Ω, w " ´uδ p¨´x 0 q on BΩ.
Then, the solution is given by: upxq " u δ pxq `wpxq "

1 4π ˆ´ln }x}I 2 `x ¨tx }x} 2 ˙F `wpxq, ppxq " p δ pxq `πpxq " x ¨F 2π}x} 2 `πpxq.
Moreover, it is easy to show that u δ R H 1 0 pΩq 2 and p δ R L 2 pΩq. Actually, we can specify how the quantity |u δ | 1,q,Ω goes to infinity when q goes to 2, with q ă 2 (which will be noted q Ñ ă 2).

According to (3), estimating |u

δ | 1,q,Ω when q Ñ ă 2 is reduced to estimate |u δ | 1,q,B
, where B denotes the ball Bpx 0 , 1q: we can easily show that there exists C ą 0 depending only on F such that @1 ď q ă `8, u δ P L q pΩq and |∇u δ | ď C }x} , and so, using polar coordinates, we get for q ă 2,

|u δ | q 1,q,Ω " ż B |∇u δ pxq| q dx ď ż B C q }x} q dx " C q ż 1 0 ż 2π 0 1 r q´1 dθdr " 2πC q 1 2 ´q .
Finally, there exists C ą 0 independent of q such that, for 1 ď q ă 2,

|u δ | 1,q,Ω ď C ? 2 ´q . (4)
In the same way, we can easily show that there exists C ą 0 independent of q such that, for 1 ď q ă 2,

|p δ | 0,q,Ω ď C ? 2 ´q . (5)
2. [START_REF] Bertoluzza | The discrete commutator property of approximation spaces[END_REF] Arnold and Liu Theorem.

Before stating Arnold and Liu Theorem, let us enumerate the assumptions that the finite element spaces V k h and W k h have to satisfy so that the theorem is true. Assumption 2. Given two fixed concentric spheres B 0 and B with B 0 ĂĂ B ĂĂ Ω, there exists an h 0 such that for all 0 ă h ď h 0 , we have for some integers k 1 and k 2 :

B1 For any 1 ď ℓ, for each v P H ℓ pBq 2 , there exists η P V k h such that }v ´η} 1,B ď Ch r 1 ´1}v} ℓ,B , r 1 " minpk 1 `1, ℓq.
For any 0 ď s, for each π P H s pBq, there exists ξ P W k h such that }π ´ξ} 0,B ď Ch r 2 }π} s,B , r 2 " minpk 2 `1, sq.

Moreover, if v P H 1 0 pB 0 q 2 (respectively π vanishes on BzB 0 ) then η (respectively ξ) can be chosen to satisfy η P H 1 0 pBq 2 (respectively ξ vanishes on ΩzB).

B2 Let ϕ P C 8 0 pB 0 q, v h P V k h and π h P W k h , then there exist η P V k h Ş H 1 0 pBq and ξ P W k h with supp ξ Ă B such that }ϕv h ´η} 1,B ď Cpϕ, B, B 0 qh}v h } 1,B , }ϕπ h ´ξ} 0,B ď Cpϕ, B, B 0 qh}π h } 0,B .
B3 For each 0 ă h ď h 0 there exists a domain B h with B 0 ĂĂ B h ĂĂ B such that for any 0 ď ℓ, for all v h P V k h and π h P W k h , we have

}v h } 1,B h ď Ch ´1´ℓ }v h } ´ℓ,B h , }π h } 0,B h ď Ch ´ℓ}π h } ´ℓ,B h .
B4 There exists β ą 0 such that for all 0 ă h ď h 0 , there is a domain B h , with B 0 ĂĂ B h ĂĂ B for which

inf π h PW k h supp π h ĂB h sup v h PV k h supp v h ĂB h ş B h divpv h qπ h |π h | 0,B h |v h | 1,B h ě β ą 0.
We now state the following theorem by Arnold and Liu [START_REF] Arnold | Local error estimates for finite element discretizations of the Stokes equations[END_REF], a key tool in the forthcoming proof of Theorem 1.

Theorem (Arnold and Liu [START_REF] Arnold | Local error estimates for finite element discretizations of the Stokes equations[END_REF]). Consider Ω 0 ĂĂ Ω 1 ĂĂ Ω, V k h and W k h satisfy Assumption 2. Suppose that pv, πq P H 1 pΩq 2 ˆL2 pΩq satisfies pv, πq | Ω 1 P H ℓ pΩ 1 q 2 ˆHℓ´1 pΩ 1 q for some ℓ ą 0. Suppose that pv h , π h q P V k h ˆW k h satisfies ş Ω π ´πh " 0 and ż Ω ∇pv ´vh q :: ∇η ´żΩ pπ ´πh qdivpηq " 0 for all η P V k h , ż Ω divpv ´vh qξ " 0 for all ξ P W k h .

Let t be a nonnegative integer. Then there exist a constant C ą 0 and a real h 1 ą 0 depending only on Ω 1 , Ω 0 , and t, such that if 0 ă h ď h 1 we have

}v ´vh } 1,Ω 0 `}π ´πh } 0,Ω 0 ďCph r 1 ´1}v} ℓ,Ω 1 `hr 2 ´1}π} ℓ´1,Ω 1 `}v ´vh } ´t,Ω 1 `}π ´πh } ´t´1,Ω 1 q,
where r 1 " minpk 1 `1, ℓq, r 2 " minpk 2 `2, ℓq, and k 1 , k 2 as in Assumption B1.

Assumptions B1 and B3 are quite standard and satisfied by a wide class of finite element spaces, including all finite element spaces defined on quasi-uniform meshes [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. The parameters k 1 and k 2 play respectively the role of the order of approximation of the spaces V k h and W k h . In our paper, for k ě 2, we will have k 1 " k and k 2 " k ´1, and for k " 1, we will have k 1 " k 2 " k " 1. Assumption B2 is less common but also satisfied by a wide variety of approximation spaces, including the P 1 bfinite elements [START_REF] Arnold | Local error estimates for finite element discretizations of the Stokes equations[END_REF]. Actually, for Lagrange finite elements, a stronger property than assumption B2 is shown in [START_REF] Bertoluzza | The discrete commutator property of approximation spaces[END_REF]: let 0 ď s ď ℓ ď k, ϕ P C 8 0 pBq and v h P V k h , then there exists η P V k h such that

}ϕv h ´η} s,B ď Cpϕqh ℓ´s`1 }v h } ℓ,B . (6)
Applied for s " ℓ " 1, inequality [START_REF] Brezis | Analyse fonctionnelle[END_REF] gives assumption B2. When B h " Ω, Assumption B4 is the standard stability condition or discrete inf-sup condition of the Stokes elements. It usually holds as long as B h is a union of elements.

Remark 1. The assumption pv, πq P H 1 pΩq 2 ˆL2 pΩq is not necessary, but it ensures that the finite element solution pv h , π h q is well-defined. In the Dirac right-hand side case, as V k h Ă C pΩq, the discrete solution pu h , p h q is well-defined and Arnold and Liu Theorem holds.

Proof of Theorem 1.

This section is devoted to the proof of our main result: Theorem 1. First, we show a weak version of Aubin-Nitsche duality lemma (Lemma 1), then we establish two discrete inf-sup conditions (Lemmas 2 and 3), and finally we use these results to prove Theorem 1.

3.1 Aubin-Nitsche duality lemma with a singular source term.

The proof of Theorem 1 is based on Arnold and Liu Theorem. In order to estimate the quantities }u ´uh } ´t,Ω 1 and }p ´ph } ´t´1,Ω 1 , we will first show a weak version of Aubin-Nitsche Lemma in the case of the Stokes Problem with a singular source term. Lemma 1. Consider f P W ´1,q pΩq 2 " pW 1,q 1 0 pΩq 2 q 1 , 1 ă q ă 2, and let pw, πq P W 1,q 0 pΩq ˆLq pΩq be the unique solution of

$ & % ´△w `∇π " f in Ω, divpwq " 0 in Ω, w " 0 on BΩ.
Let pw h , π h q be the Galerkin projection of pw, πq in

V k h ˆW k h . For any integer 0 ď t ď k ´1, }w´w h } ´t,Ω `}π ´πh } ´t´1,Ω ď Ch 2p1{q 1 ´1{2q h t`1 p|w ´wh | 1,q,Ω `|π ´πh | 0,q,Ω q . ( 7
)
Proof. We aim at estimating, for t ě 0, the H ´tpΩq-norm and the H ´t´1 pΩq-norm respectively of the errors w ´wh and π ´πh :

}w ´wh } ´t,Ω " sup ϕPC 8 0 pΩq 2 1 }ϕ} t,Ω ˇˇˇż Ω pw ´wh q ¨ϕˇˇˇˇ( 8)
}π ´πh } ´t´1,Ω " sup

ψPC 8 0 pΩq 1 }ψ} t`1,Ω ˇˇˇż Ω pπ ´πh qψ ˇˇˇ( 9)
The Galerkin projection pw h , π h q satisfies ż Ω π ´πh " 0 and ż Ω ∇pw ´wh q :: ∇η ´żΩ pπ ´πh qdivpηq " 0 for all η P V k h , ż Ω divpw ´wh qξ " 0 for all ξ P W k h .

(

) 10
Consider ϕ P C 8 0 pΩq 2 and let pw ϕ , π ϕ q P H t`2 pΩq ˆHt`1 pΩq be the solution of $ & % ´△w ϕ `∇π ϕ " ϕ in Ω, divpw ϕ q " 0 in Ω, w ϕ " 0 on BΩ.

Existence and uniqueness of the solution to this problem are given in [START_REF] Temam | Theory and numerical analysis[END_REF] (see Chapter I, §2), and we have the estimate

}w ϕ } t`2,Ω `}π ϕ } t`1,Ω ď C}ϕ} t,Ω . (11)
In dimension 2, by the Sobolev injections established for instance in [START_REF] Brezis | Analyse fonctionnelle[END_REF], we have

H t`2 pΩq Ă W 1,q 1 pΩq, H t`1 pΩq Ă L q 1 pΩq, (12)
for all q 1 in r2, `8r. Thus ż Ω pw ´wh q ¨ϕ " ´żΩ pw ´wh q ¨△w ϕ `żΩ pw ´wh q ¨∇π ϕ " ż Ω ∇pw ´wh q :: ∇w ϕ ´żΩ divpw ´wh qπ ϕ .

By adding [START_REF] Casas | Spike controls for elliptic and parabolic PDEs[END_REF] in the last equality, we get for any η P V k h and any ξ P W k h , ż Ω pw ´wh q ¨ϕ " ż Ω ∇pw ´wh q :: ∇pw ϕ ´ηq ´żΩ divpw ´wh qpπ ϕ ´ξq `żΩ divpηqpπ ´πh q.

By definition of w ϕ , divpw ϕ q " 0 on Ω, so ż Ω pw ´wh q ¨ϕ " ż Ω ∇pw ´wh q :: ∇pw ϕ ´ηq ´żΩ divpw ´wh qpπ ϕ ´ξq `żΩ divpη ´wϕ qpπ ´πh q ď|w ´wh | 1,q,Ω `|w ϕ ´η| 1,q 1 ,Ω `|π ϕ ´ξ| 0,q 1 ,Ω |π ´πh | 0,q,Ω |w ϕ ´η| 1,q 1 ,Ω .

Now let us deal with the pressure estimate. For any ψ P C 8 0 pΩq, we denote by r ψ the function

r ψ " ψ ´1 |Ω| ż Ω ψ.
By definition, it is easy to see that r ψ satisfies ż Ω r ψ " 0 and @t ě 0, } r ψ} t`1,Ω ď CpΩq}ψ} t`1,Ω .

We can now establish the result for the pressure: consider ψ P C 8 0 pΩq and let pw ψ , π ψ q P H t`2 pΩq Ĥt`1 pΩq be the solution of $ & % ´△w ψ `∇π ψ " 0 in Ω, divpw ψ q " r ψ in Ω, w ψ " 0 on BΩ, See [START_REF] Temam | Theory and numerical analysis[END_REF] (Chapter I, §2) for the existence and the uniqueness, and the following estimate

}w ψ } t`2,Ω `}π ψ } t`1,Ω ď C} r ψ} t`1,Ω ď C}ψ} t`1,Ω . (13)
Moreover,

ż Ω π ´πh " 0, so that ż Ω pπ ´πh qψ " ż Ω pπ ´πh q r ψ `1 |Ω| ż Ω ψ ż Ω π ´πh " ż Ω pπ ´πh q r ψ.
By the Sobolev injections recalled in [START_REF] Fulford | Muco-ciliary transport in the lung[END_REF], and the Galerkin projection property [START_REF] Casas | Spike controls for elliptic and parabolic PDEs[END_REF], we can write for all η P V k h , ż Ω pπ ´πh qψ " ż Ω pπ ´πh q r ψ " ż Ω pπ ´πh qdivpw ψ q " ż Ω pπ ´πh qdivpw ψ ´ηq `żΩ ∇pw ´wh q :: ∇η.

Then, for all v P W 1,q 0 pΩq ż Ω ∇w ψ :: ∇v ´żΩ π ψ divpvq " 0, so, with v " w ´wh , and for any ξ P W k h , ż Ω pπ ´πh qψ " ż Ω pπ ´πh qdivpw ψ ´ηq `żΩ ∇pw ´wh q :: ∇pη ´wψ q `żΩ π ψ divpw ´wh q " ż Ω pπ ´πh qdivpw ψ ´ηq `żΩ ∇pw ´wh q :: ∇pη ´wψ q `żΩ pπ ψ ´ξqdivpw ´wh q ď|π ´πh | 0,q,Ω |w ψ ´η| 1,q 1 ,Ω `|w ´wh | 1,q,Ω ´|w ψ ´η| 1,q 1 ,Ω `|π ψ ´ξ| 0,q 1 ,Ω ¯.

Finally, for any pη 1 , ξ

1 q P V k h ˆW k h , ż Ω pw ´wh q ¨ϕ ď|w ´wh | 1,q,Ω `|w ϕ ´η1 | 1,q 1 ,Ω `|π ϕ ´ξ1 | 0,q 1 ,Ω |π ´πh | 0,q,Ω |w ϕ ´η1 | 1,q 1 ,Ω , (14)
and for any pη 2 , ξ

2 q P V k h ˆW k h , ż Ω pπ ´πh qψ ď|π ´πh | 0,q,Ω |w ψ ´η2 | 1,q 1 ,Ω `|w ´wh | 1,q,Ω ´|w ψ ´η2 | 1,q 1 ,Ω `|π ψ ´ξ2 | 0,q 1 ,Ω ¯. ( 15
)
In order to estimate

|w ϕ ´η1 | 1,q 1 ,Ω , |w ψ ´η2 | 1,q 1 ,Ω , |π ϕ ´ξ1 | 0,q 1
,Ω and |π ψ ´ξ2 | 0,q 1 ,Ω , we will need the following result: Proposition 1 (Girault, Raviart, Corollary A.2, page 97 [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF]). Let T h be a family of quasi-uniform simplicial triangulations of Ω Ă R 2 , where h is the meshsize. For any 0 ď m ď t `1 ď k, for any mesh element T in the family, for any v P W k`1,q 1 pΩq, q 1 ě 2 real,

|v ´Πh v| m,q 1 ,T ď Ch 2p1{q 1 ´1{2q h t`2´m |v| t`2,2,T , (16)
where Π h v is the P k -interpolant of the function v.

Up to now and until the end of this proof, we will take

η 1 " Π h w ϕ and η 2 " Π h w ψ P V k h , ξ 1 " r Π h π ϕ and ξ 2 " r Π h π ψ P W k h ,
where Π h v is the P k -interpolant of the function v and r Π h v is the P k´1 -interpolant of the function v. By [START_REF] Lacouture | A numerical method to solve the Stokes problem with a punctual force in source term[END_REF], with m " 1, 0 ď t ď k ´1, for all T finite element in the family,

|w ϕ ´η1 | 1,q 1 ,T ď Ch 2p1{q 1 ´1{2q h t`1 |w ϕ | t`2,2,T , (17)
|w ψ ´η2 | 1,q 1 ,T ď Ch 2p1{q 1 ´1{2q h t`1 |w ψ | t`2,2,T ,
and with m " 0,

|π ϕ ´ξ1 | 0,q 1 ,T ď Ch 2p1{q 1 ´1{2q h t`1 |π ϕ | t`1,2,T , |π ψ ´ξ2 | 0,q 1 ,T ď Ch 2p1{q 1 ´1{2q h t`1 |π ψ | t`1,2,T .
We denote the triangles of the mesh by tT i u i"1,¨¨¨,N , and we set a " pa i q i and b " pb i q i , where a i " |w ϕ ´η1 | 1,q 1 ,T i and b i " |w ϕ | t`2,2,T i .

By [START_REF] Leykekhman | Optimal error estimates for finite element discretization of elliptic optimal control problems with finitely many pointwise state constraints[END_REF], we have, for all i in rr1, N ss,

a i ď Ch 2p1{q 1 ´1{2q h t`1 b i .
We recall the norm equivalence in R N for 0 ă r ă s,

}x} ℓ s ď }x} ℓ r ď N 1{r´1{s }x} ℓ s ,
with here N " Ch ´2. As 2 ă q 1 , we have }b} ℓ q 1 ď }b} ℓ 2 . Then, we can write

|w ϕ ´η1 | 1,q 1 ,Ω " }a} ℓ q 1 ď Ch t`1 h 2p1{q 1 ´1{2q }b} ℓ q 1 ď Ch t`1 h 2p1{q 1 ´1{2q }b} ℓ 2 ď Ch t`1 h 2p1{q 1 ´1{2q |w ϕ | t`2,2,Ω .
Similarly, we get

|w ψ ´η2 | 1,q 1 ,Ω ď Ch t`1 h 2p1{q 1 ´1{2q |w ψ | t`2,2,Ω , |π ϕ ´ξ1 | 0,q 1 ,Ω ď Ch t`1 h 2p1{q 1 ´1{2q |π ϕ | t`1,2,Ω , |π ψ ´ξ2 | 0,q 1 ,Ω ď Ch t`1 h 2p1{q 1 ´1{2q |π ψ | t`1,2,Ω ,
and by ( 11) and ( 13), we get

|w ϕ ´η1 | 1,q 1 ,Ω ď Ch t`1 h 2p1{q 1 ´1{2q }ϕ} t,2,Ω , |w ψ ´η2 | 1,q 1 ,Ω ď Ch t`1 h 2p1{q 1 ´1{2q }ψ} t`1,2,Ω , |π ϕ ´ξ1 | 0,q 1 ,Ω ď Ch t`1 h 2p1{q 1 ´1{2q }ϕ} t,2,Ω , |π ψ ´ξ2 | 0,q 1 ,Ω ď Ch t`1 h 2p1{q 1 ´1{2q }ψ} t`1,2,Ω ,
Finally, the proof is ended by combining ( 8), ( 9), ( 14), [START_REF] Köppl | Optimal a priori error estimates for an elliptic problem with Dirac right-hand side[END_REF], and the last inequalities.

Corollary 1. Let pu h , p h q P V k h ˆW k h be the Galerkin projection of the solution pu, pq of Problem (1), for any 0 ă ε ă 1,

}u´u h } ´k`1,Ω `}p ´ph } ´k,Ω ď Ch ´εh k p|u ´uh | 1,qε,Ω `|p ´ph | 0,qε,Ω q ,
where q ε P r1, 2r is defined by

q ε " 2 1 `ε ˆand so q 1 ε " 2 1 ´ε ˙. (18)
Proof. We will apply Lemma 1 with F " δ x 0 F, w " u, π " p and t " k ´1. We can explicit inequality (7):

2 ˆ1 q 1 ε ´1 2 ˙" 2 ˆ1 ´ε 2 ´1 2 ˙" ´ε, (19)
and so, it follows

}u´u h } ´k`1,Ω `}p ´ph } ´k,Ω ď Ch ´εh k p|u ´uh | 1,qε,Ω `|p ´ph | 0,qε,Ω q .
3.2 Discrete inf-sup conditions in L qε -norm.

Section 3.3 is devoted to estimate of |u ´uh | 1,qε,Ω and |p ´ph | 0,qε,Ω . In that prospect, we need to establish two discrete inf-sup conditions.

Lemma 2. With q ε and q 1 ε defined in [START_REF] Pozrikidis | Boundary integral and singularity methods for linearized viscous flow[END_REF], the approximation space V k h defined by

V k h " " v h P V k h ˇˇˇż Ω divpv h qp h " 0, @p h P W k h * ,
satisfies the following discrete inf-sup condition:

inf u h P V k h sup v h P V k h ş Ω ∇u h :: ∇v h |u h | 1,qε,Ω |v h | 1,q 1 ε ,Ω ě Ch ε .
Proof. The bilinear form apu, vq " ż Ω ∇u :: ∇v is continuous and coercive on H 1 0 pΩq, so for V k h vector subspace of H 1 0 pΩq, we have the inf-sup condition: inf

u h P V k h sup v h P V k h ş Ω ∇u h :: ∇v h |u h | 1,Ω |v h | 1,Ω ě α ą 0,
where α only depends on Ω. We recall the following inverse inequality:

Proposition 2 (Ciarlet, Theorem 3.2.6, page 140 [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]). Let T h a family of quasi-uniform simplicial triangulations of Ω Ă R d , where h is the meshsize. For

v h P V k h , 1 ď r, s ă `8, 0 ď ℓ ď m, ˜ÿ T PT h |v h | r m,r,T ¸1{r ď Ch ´drmaxt0,1{s´1{rus h ´pm´ℓq ˜ÿ T PT h |v h | s ℓ,s,T ¸1{s
.

We will apply this to any v h P V k h Ă C pΩq, with d " 2, m " l " 1, s " 2 and r " q 1 ε to get:

|v h | 1,q 1 ε ,Ω ď Ch ´2p1{2´1{q 1 ε q |v h | 1,2,Ω " Ch ´ε|v h | 1,2,Ω ,
using equality [START_REF] Rannacher | A priori error estimates for the finite element discretization of elliptic parameter identification problems with pointwise measurements[END_REF]. Moreover, for any

u h P V k h , |u h | 1,qε,Ω ď C|u h | 1,2,Ω ď C sup v h P V k h ş Ω ∇u h :: ∇v h |v h | 1,2,Ω ď Ch ´ε sup v h P V k h ş Ω ∇u h :: ∇v h |v h | 1,q 1 ε ,Ω Finally, inf u h P V k h sup v h P V k h ş Ω ∇u h :: ∇v h |u h | 1,qε,Ω |v h | 1,q 1 ε ,Ω ě Ch ε .
The second discrete inf-sup condition we need is given by the following lemma: Lemma 3. With q ε and q 1 ε defined in [START_REF] Pozrikidis | Boundary integral and singularity methods for linearized viscous flow[END_REF], the approximations spaces V k h and W k h satisfy the following discrete inf-sup condition:

inf p h PW k h sup v h PV k h ş Ω divpv h qp h |p h | 0,qε,Ω |v h | 1,q 1 ε ,Ω ě Ch ε .
Proof. The proof is similar to the proof of Lemma 2. According to Assumption B4,

inf p h PW k h sup v h PV k h ş Ω divpv h qp h |p h | 0,Ω |v h | 1,Ω ě β ą 0.
According to Proposition 2, for any

v h P V k h , |v h | 1,q 1 ε ,Ω ď Ch ´ε|v h | 1,2
,Ω . So, we have, for any p h P W k h and q ε ă 2,

|p h | 0,qε,Ω ď C|p h | 0,Ω ď C sup v h PV k h ş Ω divpv h qp h |v h | 1,2,Ω ď Ch ´ε sup v h PV k h ş Ω divpv h qp h |v h | 1,q 1 ε ,Ω
Finally, we get inf Proof. First, we will estimate |u h | 1,qε,Ω in terms of |u| 1,qε,Ω . As we have divpuq " 0 on Ω, by [START_REF] Arnold | Local error estimates for finite element discretizations of the Stokes equations[END_REF] we have ż Ω divpu h qq h " 0, @q h P W k h , and so, u h P V k h . According to Lemma 2, there exists v h P V k h such as |v h | 1,q 1 ε ,Ω " 1, and We can now estimate |p ´ph | 0,qε,Ω .

p h PW k h sup v h PV k h ş Ω divpv h qp h |p h | 0,qε,Ω |v h | 1,q 1 ε ,Ω ě Ch ε .
|u h | 1,
Lemma 5. Let pu h , p h q P V k h ˆW k h be the Galerkin projection of the solution pu, pq of Problem (1), for any small enough real ε ą 0, |p ´ph | 1,qε,Ω ď Ch ´2ε p|u| 1,qε,Ω `|p| 0,qε,Ω q .

Proof. The proof is similar to the velocity case: according to Lemma 3, there exists

v h P V k h such as |v h | 1,q 1 ε ,Ω " 1 and |p h | 0,qε,Ω ď Ch ´ε ż Ω divpv h qp h .
By (2), we have ż Ω divpv h qp h " ´żΩ ∇pu ´uh q :: ∇v h `żΩ divpv h qp.

By applying Lemma 4, as |v

h | 1,q 1 ε ,Ω " 1, we get |p h | 0,qε,Ω ď Ch ´ε ˆ´ż Ω ∇pu ´uh q :: ∇v h `żΩ divpv h qp ˙, ď Ch ´ε p|u ´uh | 1,qε,Ω `|p| 0,qε,Ω q ,
ď Ch ´2ε p|u| 1,qε,Ω `|p| 0,qε,Ω q .

3.4 Proof of Theorem 1.

We can now prove Theorem 1.

Proof. The functions u and p are analytic on Ω 1 , so the quantities }u} k`1,Ω 1 and }p} k,Ω 1 are bounded.

Let us note that in this case pu, pq R H 1 0 pΩq 2 ˆL2 0 pΩq, but Remark 1 allows us to apply Arnold and Liu Theorem. For k 1 " k and

k 2 " " 1 if k " 1, k ´1 if k ě 2,
and l " k `1 and t " k ´1, we have

}u ´uh } 1,Ω 0 `}p ´ph } 0,Ω 0 ď Cph k `}u ´uh } ´k`1,Ω 1 `}p ´ph } ´k,Ω 1 q.
By combining Corollary 1, Lemmas 4 and 5, and inequalities (4) and ( 5), we get

}u ´uh } ´k`1,Ω 1 `}p ´ph } ´k,Ω 1 ď Ch k h ´3ε ? 2 ´qε .
According to [START_REF] Pozrikidis | Boundary integral and singularity methods for linearized viscous flow[END_REF], with ε ď 1,

1 ? 2 ´qε " ? 1 `ε ? 2ε ď 1 ? ε , therefore, taking ε " | ln h| ´1, }u ´uh } 1,Ω 0 `}p ´ph } 0,Ω 0 ď Ch k a | ln h|,
which ends the proof of Theorem 1.

4 General case.

Theorem 1 and its proof have been written in the particular case of the P k {P k´1 finite element method, k ě 2, and the P 1 b{P 1 elements (which corresponds to the case k " 1). But we can state a more general result.

First, let us focus on the assumptions: let T h be a family of quasi-uniform simplicial triangulations of Ω, let V k 1 h and W k 2 h be two approximation spaces satisfying Assumption 2. We will also assume that V k 1 h P C pΩq: this assumption ensures that the finite element solution is well-defined. Moreover, we will need two more assumptions, they will play the role of Propositions 1 and 2: Assumption 3. Given B Ă Ω, consider q 1 ě 2, there exists an h 0 such that for all 0 ă h ď h 0 , we have for some positive integers k 1 and k 2 : r B1 For any 1 ď ℓ, for each v P H ℓ pBq 2 , there exists η P V k 1 h such that, for any mesh element T Ă B, |v ´η| 1,q 1 ,T ď Ch dp1{q 1 ´1{2q h r 1 ´1|v| ℓ,2,T , r 1 " minpk 1 `1, ℓq.

For any 0 ď s, for each π P H s pBq, there exists ξ P W k 2 h such that, for any mesh element T Ă B, |π ´ξ| 0,q 1 ,T ď Ch dp1{q 1 ´1{2q h r 2 |π| s,2,T , r 2 " minpk 2 `1, sq.

r B3 For all v h P V k 1 h , for any mesh element T P T h , we have }v h } 1,q 1 ,T ď Ch 2p1{q 1 ´1{2q }v h } 1,2,T .

Assumptions r B1 and r B3 are also satisfied by a wide class of finite element spaces, including all finite element spaces defined on quasi-uniform meshes [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. They are actually common generalisations of Assumptions B1 and B3.

We can now state the following result:

Theorem 2. Consider Ω 0 ĂĂ Ω 1 ĂĂ Ω satisfying Assumption 1, 1 ď q ă 2, let pu, pq P W 1,q 0 pΩq Lq 0 pΩq be the solution of Problem (1) and pu h , p h q its Galerkin projection onto V k 1 h ˆW k 2 h satisfying ş Ω p h " 0 and ż Ω ∇pu ´uh q :: ∇η ´żΩ pp ´ph qdivpηq " 0 for all η P V k h , ż Ω divpu ´uh qξ " 0 for all ξ P W k h .

Under the assumption that pu, pq P H k 0 `1pΩ 1 q 2 ˆHk 0 pΩ 1 q, there exists h 1 such that if 0 ă h ď h 1 , we have, }u ´uh } 1,Ω 0 `}p ´ph } 0,Ω 0 ď CpΩ 0 , Ω 1 , Ωqh k 0 a | ln h|,

where k 0 " minpk 1 , k 2 `1q.

About the error in L 2 pΩ 0 q-norm for the velocity, Figure 12 suggests that the P k {P k´1 finite element method (or P 1 b{P 1 if k " 1) converges at the order k `1 on Ω 0 . This result has only been observed numerically but it is still an open question.

3. 3 Lemma 4 .

 34 Estimates of |u ´uh | 1,qε,Ω and |p ´ph | 0,qε,Ω .Following Corollary 1, the quantities |u ´uh | 1,qε,Ω and |p ´ph | 0,qε,Ω have to be estimated to prove Theorem 1. We will apply the last two results to bound them in terms of |u| 1,qε,Ω and |p| 0,qε,Ω . Let pu h , p h q P V k h ˆW k h be the Galerkin projection of the solution pu, pq of Problem (1), for any small enough real ε ą 0, |u ´uh | 1,qε,Ω ď Ch ´ε p|u| 1,qε,Ω `|p| 0,qε,Ω q .

Elements P 1 b{P 1 Elements P 2 {P 1 Elements P 3 {P 2 OrderFigure 10 :

 1213210 Figure10: Estimated order of convergence for the H 1 pΩ 0 q-norm of the velocity.

  qε,Ω ď Ch ´ε ż We conclude with the triangulary inequality, |u ´uh | 1,qε,Ω ď |u| 1,qε,Ω `|u h | 1,qε,Ω ď Ch ´ε p|u| 1,qε,Ω `|p| 0,qε,Ω q .

	Moreover, equality (2) gives
	ż	ż
		Ω	∇u :: ∇v h	´żΩ	divpv h qpp ´ph q.
	Now, v h P V k h , so	ż

Ω ∇u h :: ∇v h . Ω ∇u h :: ∇v h " Ω divpv h qp h " 0. Finally, as |v h | 1,q 1 ε ,Ω " 1,

we get |u h | 1,qε,Ω ď Ch ´ε ˆżΩ ∇u :: ∇v h ´żΩ divpv h qp ˙, ď Ch ´ε p|u| 1,qε,Ω `|p| 0,qε,Ω q .

Proof. We will not develop the complete proof here because it is essentially the same as the proof of Theorem 1 (see Section 3). But we will explain two differences between the both proofs:

' the result of Lemma 1 holds in this case, but for 0 ď t ď minpk 1 ´1, k 2 q.

' the result of Corollary 1 becomes

where k 0 " minpk 1 , k 2 `1q.

The end of the proof is the same.

Numerical illustrations.

In this section, we present some computations which illustrate the theoretical results proved in this paper.

Concentration of the error around the singularity. First, we define Ω as the unit square, Ω " r0, 1s 2 . and solve the following Stokes problem with F " t r1, 1s and x 0 " p0.5, 0.5q,

where u δ is the 2d Stokeslet defined in (3). [START_REF] Araya | A posteriori error estimates for elliptic problems with Dirac delta source terms[END_REF], Problem (20) has non homogeneous Dirichlet boundary conditions, but in this case, the exact solution is known: u δ . Thus, it is easier to get some information on the error. Estimated orders of convergence. For this second example, the domain Ω is still the unit square, and Ω 0 is defined as the following portion of Ω, Ω 0 " tx P Ω : }x ´x0 } 2 ą 0.4u , where x 0 " p0.5, 0.5q. We fix F " t r1, 1s and solve Problem 1 for different mesh sizes h with the P 1 b{P 1 , P 2 {P 1 and P 3 {P 2 finite element methods.

Remark 2. Unlike Problem

Figure 10 (respectively Figure 11) presents the estimated orders of convergence for the H 1 pΩ 0 qnorm of the velocity (respectively the L 2 pΩ 0 q-norm of the pressure) for these three methods. The convergence far from the singularity (i.e. on Ω 0 ) is the same as in the regular case: the P k {P k´1 method (or the P 1 b{P 1 method if k " 1) converges at the order k on Ω 0 in H 1 -norm for the velocity and in L 2 -norm for the pressure, as proved in this paper. Let us just note that there is an over-convergence in pressure for the P 1 b{P 1 elements: the estimated order of convergence is approximately 2, greater than the convergence expected by Theorem 1.

Elements P 1 b{P 1 Elements P 2 {P 1 Elements P 3 {P 2 Order = 2.17