Abdelmajid Amine
email: amine@univ-paris12.fr

Shérazade Gatfaoui
email: gatfaoui_sherazade@yahoo.fr

A DYNAMIC ANALYSIS OF TRUST IN CUSTOMER-BANK RELATIONSHIP: A RETROSPECTIVE CASE STUDIES APPROACH

Keywords:

This research analyzes trust building within the customer/bank relationship. It introduces the time dimension to study the mechanisms and the development processes of trust, both institutional and interpersonal, during the customer/bank relationship. A qualitative and longitudinal study was undertaken based on retrospective case studies conducted within a French bank. This research highlights the general development process of interpersonal trust and the existence of two modes of trust construction -linear and by stagesdepending on the occurrence of incidents, their nature and their resolution by the bank. Our results show also a shift from institutional to interpersonal trust that occurs during the customer/bank relationship.

Introduction

Trust is a notion that has provoked considerable interest from both researchers and practitioners in the fields of psychology, sociology, economics and more recently marketing and management. The concept is at the core of numerous studies in particular with regard to its nature, its status, and how it is built up in business relationships.

If we try to draw up a synthesis of theoretical work on trust, it is possible to identify two perspectives in studying this concept. The first approach, a 'static' one, has tried to study and measure the concept of trustits dimensions, antecedents and consequencesat a particular moment T. The second approach is more dynamic, but much less developed in management, and aims at ensuring a better understanding of the processes by which trust is built over the course of time (Ring and [START_REF] Peter | Structuring cooperative relationships between organizations[END_REF][START_REF] Zajac | From Transaction Cost to Transactional Value Analysis: Implications for study of Interorganizational Strategies[END_REF].

In this study we address the following main question : "How trust is built up during the customer/bank relationship ? What are the mechanisms that contribute to the emergence and the reinforcement of trust over time ?"

Drawing on the literature, the first part of the paper shows the difficulty of conceptualizing trust through the examination of its various definitions, types and development processes. We then present the research methodology based on a qualitative and longitudinal study conducted in a French bank. Finally the results are discussed and the main implications, including the time dimension in the analysis of trust, are drawn.

WHY IS IT RELEVANT TO STUDY TRUST THROUGH A DYNAMIC APPROACH?

In the marketing literature, interpersonal trust is variously defined as a belief, an expectation, a feeling, a behavioral intention or a willingness to rely on an exchange partner within perceived risky situations [START_REF] Moorman | Relationship between providers and users of market research : The dynamics of trust within and between organizations[END_REF]Chaudhari and Holbrook, 2001)). Some authors consider that trust is made up either of two componentscredibility and benevolence [START_REF] Doney | An examination of the nature of trust in buyer-seller relationships[END_REF][START_REF] Ganesan | Determinants of long-term orientation in buyer-seller relationship[END_REF]Morgan and Hunt, 1994)or of threecredibility, honesty and benevolence [START_REF] Hess | Construction and assessment of a scale to measure consumer trust[END_REF]. The definitions provided do not explicitly include the question of the dynamic of trust, its development on the basis of the history or shared experience of the partners, or to any incidents encountered during the customer/firm relationship. This account seems consistent in so far as the studies have often aimed at defining and measuring trust at a given moment (T) by means of measurement scales, while ignoring the consumer's perception in a qualitative and longitudinal perspective [START_REF] Dwyer | Developing buyer-seller relationships[END_REF][START_REF] Peter | Development processes of cooperative intreorganizational relationships[END_REF][START_REF] Zajac | From Transaction Cost to Transactional Value Analysis: Implications for study of Interorganizational Strategies[END_REF]. Faced with the growing number of conceptions and definitions of trust, certain authors have also wondered whether trust should be viewed as interpersonal, based on individuals knowing one another, rather than as institutional, based on the reputation of the firm or institution [START_REF] Williamson | Calculativeness, trust and economic organization[END_REF][START_REF] Young | The role of trust and cooperation in marketing channels: A preliminary study[END_REF].

The literature review also allows us to identify various determinants and consequences of trust according to its contexts of use, namely industry, services, retailing, and bank-to-bank or bank-to-customer relationships [START_REF] Andaleeb | Trust and dependence in channel relationships: Implications for satisfaction and perceived stability[END_REF][START_REF] Geyskens | An investigation into the joint effects of trust and interdependance on relationship commitment[END_REF][START_REF] Moorman | Factors affecting trust in market research relationships[END_REF]. However, some variables such as satisfaction, cooperation and communication appear either as antecedents or consequences of trust. This problem occurs in particular because the authors don't indicate the stage of development of the relation being studied.

Thus, analyzed from a static view of trust, these determinants and consequences hide/mask the reality and complexity of trust and precisely the mechanisms of its construction during a relation. Nevertheless, because trust is by nature a dynamic and contextualized concept, it is relevant to take into account the time dimension of the development of the relation when we analyze its underlying processes. One can hypothesize that this dynamic vision of trust may reconcile the former conflicting findings about the double status of some factors (as both antecedents and consequences) that emerge in the static approach and thus contribute to a better understanding of the mechanisms of trust building. This shift in the conception of trust must be accompanied by a change in the method used to analyze this construct.

This line of inquiry is all the more interesting since it concerns the banking services sector characterized by a high level of perceived risk ,depending on the intangibility, simultaneity of production and consumption, and variability of services.

RESEARCH METHODOLOGY

A qualitative and longitudinal study was adopted, based on 18 retrospective case studies [START_REF] Yin | Case study research, design and methods[END_REF] conducted within a French bank. Each case focuses on the relationship between a customer and his financial adviser, at the time when the study was carried out. The critical incident method enables us to trace back the historical relationship between a customer and his adviser [START_REF] John | The critical incident technique[END_REF]Bitner, Booms and Treteault, 1990 ;[START_REF] Keaveney | Customer switching behavior in service industries: An exploratory study[END_REF]. Further, identifying the incidents that come up between a bank and its customers as well as the way they are handled allow us to understand more clearly the development processes of trust and the mechanisms by which it is built. Finally, the dyadic approach allows us to trace back the historical relationship as perceived by the customer and his financial adviser [START_REF] Brown | A Gap Analysis of Professional Service Quality[END_REF][START_REF] Deshpandé | Corporate Culture, Customer Orientation, and Innovativeness in Japanese Firms : A Quadrad Analysis[END_REF][START_REF] John | The reliability and validity of key informant data from dyadic relationships in marketing channels[END_REF].

The retrospective in-depth interviews were set up one to one and allow us to trace back the customer's bank history. The interviews were set up in seven different branches of the bank in three departments in Paris (Paris-Ouest, Paris-Est and Paris Hauts-de-Seine). Customers were selected on the basis of the four main criteria expounded in the literature: product and services portfolio, the duration of the relationship with the bank, the frequency of contacts between the customer and the adviser, and the occurrence of incidents during the customer/bank relationship. Clarifying conceptions of trust in the private customer/bank relationship necessarily involved choosing those customers having a reasonably substantial bank product and services portfolio. Customers needed to have sufficient experience of bank services to be able to discuss trust, and to have been in the adviser's customer portfolio for more than a year.

In addition, the customers selected were required to have maintained regular contact with their existing adviser. Financial advisers were also required to have dealt at least once with both negative events (either incidents in the customer's private life or critical incidents provoked by the bank) and positive events in the course of their relation with the customer.

Incidents in the customer's private life may have given rise to problems at the level of his financial situation and of the relationship between the customer and the bank (for example : a period of unemployment, a divorce, a death). Incidents provoked by the bank may involve an error involving the transfer of funds to another customer's account, a "forced" sale of a financial product or service, or a refusal to a request for a loan.

In the course of the interviews, we questioned nine customers aged from 28 to 45 and nine customers aged from 51 to 69, and their financial advisers aged from 23 to 54. Customers aged from 28 to 45 are part of the active customers segment and those aged from 51 to 69 are part of the senior customers segment, as defined by the bank's classification scheme. The 18 retrospective cases studies include 81 interviews (31 financial adviser interviews and 45 customer interviews and 5 branch manager interviews). Finally, the interviews were recorded and transcribed, then analyzed.

FINDINGS & RESULTS DISCUSSION

We carried out intra-case and inter-case content analyses, highlighting three principal results.

First, we demonstrate the recursive functioning/nature of the interpersonal trust-building process (i). This kind of trust appears as the major dimension of reliance/confidence in the customer-bank relationship. Then we propose two modes of interpersonal trust construction, depending on the history of the banking relation (ii). Third, we reveal some transfer mechanisms from institutional trust to interpersonal trust in the course of the relation between the customer and the bank/adviser (iii).

i-Recursive process of interpersonal trust building

Our research shows/explains how some variables are identified in the literature either as antecedents or consequences of (interpersonal) trust. By introducing the time factor in our analysis, we succeed in explaining the double status (antecedent and consequence of trust) of certain variables such as satisfaction, communication or coordination, depending on the moment when the investigation is made and the development stage of the relation. The construction process of interpersonal trust takes a recursive form.

For instance, if a customer is satisfied with an initial experience (at moment T), his satisfaction may generate trust in the adviser and in the bank, which in turn may reinforce his satisfaction (becoming a consequence) in the next encounter (at moment T+1) (antecedent=>consequence=>antecedent ...). The following quotations illustrate this recursive process of interpersonal trust building: "Because I am satisfied with my past experience with my adviser, it is easier to have confidence in his proposals. And when things go well, I am satisfied and pleased, and I rely on him now and for the future."

Communication between customers and bank advisers and coordination of their behavior seem also to contribute to building interpersonal trust, and are in turn improved by it.

These findings show that when we focus on the behaviors and experiences of customers and advisers, we better understand the mechanisms and modes of interpersonal trust building.

They also corroborate the results of past studies, which stress that good communication between partners strengthens the tendency to mutual trust.

ii-Modes of interpersonal trust building

Our content analyses highlight two modes/types of construction of interpersonal trust, depending on the relational history between customers and bank advisers. The first mode of building interpersonal trust takes a linear form (figure 2). This mode appears when customers, in their relation to the bank, are faced with "minor" incidents or incidents that have been favourably resolved by the adviser (10 cases out of 18). In these cases, both customers and advisers agreed to say that "they enjoy and maintain a reciprocal relation based on trust"

and "that maintaining/managing the relation is quite easy". However, these partners of exchange are conscious that "institutional trust (customer/bank) has not yet been really tested". Thus customers expect that "the bank will answer positively to future needs (e.g. a credit request)" and the advisers state that they "pay attention and remain aware of the evolution of our customers' needs and expectations".

Figure 1-Linear mode of construction of interpersonal trust

E1 E2 E3 E4 E5 E6 E7 E8 E9 E10 Experience 1 (period T)………………………………………………………… … ... .. ..… Experience x (period T+x) (satisfaction) (satisfaction)
The second mode of interpersonal trust building "by stages" relates to those situations where the incidents were mainly caused by customers, and were considered to be serious and/or were not satisfactory resolved by the bank (8 cases up to 18). These cases were characterized either by an entry in relation with the present bank following incidents with a previous bank or non-satisfying resolution of major incidents encountered during the relation with the bank being studied. Thus, customers report "becoming mistrustful" toward the credibility of the adviser or the bank.

When senior management (the branch manager, for example) intervenes to solve an incident caused by the adviser, the customer calls into question interpersonal trust but not institutional confidence. These customers end up doubting and reconsider their confidence progressively in the course of their experience with the new advisers whom they will subsequently meet.

Interpersonal trust is then put to the test in relation to the behavior of the advisers during these various experiences. In such cases, the advisers maintain that "management of the relation with these customers is a daily battle". They need to be "very attentive to this type of customer, who is much more demanding than those who have not encountered a serious incident", in particular for "major customers" in terms of profitability to the branch.

Figure 2-Transition from "linear" to "by stages" mode of interpersonal trust building

"Linear" mode of interpersonal trust building "By stages" mode of interpersonal trust building before the occurrence of mismanaged incidents following mismanaged incidents

E1 E2 E3 E4 E5 E6 E7 E8 Experience n°1, 2, 3, 4, (T+1,2,3,4) (T+5) Experience n°6 (T+6) Experience n°7 (T+7) Experience 8 (T+8) …etc. (Satisfaction) Incidents (Satisfaction) (Satisfaction) (Satisfaction) not properly resolved (dissatisfaction)
Expresses the linear mode of trust building Expresses the by stages mode of trust construction

iii-Transfers of institutional confidence towards interpersonal trust during the banking relation

Finally, content analysis shows some transfer mechanisms of institutional confidence towards interpersonal trust during the banking relation, according to the age and the banking relation history of the customers.

Segment of "active" customers (28-45 years old)

Customers belonging to the "active segment" attached more importance to institutional confidence at the moment they enter in relation with the bank (7 cases out of 12). Thus they did not consider the relation with the adviser to be important, either because they did not start their relation with the bank following incidents unresolved by a previous provider or because they consider that they have little experience of banks (they were "students" or "recently began their working life").

Nevertheless, during their relation with the present bank, these customers became aware of the importance of interpersonal trust to such an extent that some of them followed their adviser when he moved to another branch. This transfer of trust was observed for customers having met with incidents which were satisfactorily resolved either at the time of entering or in the course of the relation with the current bank. So even if they say they retain institutional confidence in the bank itself, their reliance on the adviser has become more significant.

Customers belonging to the "senior" segment (over 50)

These customers also explain their entry into relation with the present bank in terms of the perceived credibility of the institution. Like the "active" customers, they were not aware of the importance of the relation with the adviser (4 cases out of 6) ). However, following their long experience of the banking relation with various advisers, these customers tend to put past incidents into perspective. They say "they have learned about the internal functioning of the bank" and incident management.. Thus, when a problem occurs which they might previously have considered to be serious, they say they are "calmer" because they know "that their adviser will solve it successfully". Even if this does not happen, they think "that the bank manager will be available if a problem arises with the adviser". They talk about their "learning about the relation with advisers and the bank". Remember that in these cases the customers have not encountered incidents deemed to be "very serious" or that "incidents have on the whole been satisfactorily handled" by their advisors and/or their superiors in the course of their relation with the bank (5 cases out of 6).

With the passage of time and through experience considered largely to be positive, these customers reveal their trust both in their adviser and their bank. Because they "trust the bank which recruits good advisers", these customers do not consider it appropriate to follow these advisers if they move to another bank.

CONCLUSION

To analyze trust in the customer/bank relationship, we adopted a time perspective which contrasts with earlier studies works that mainly use a static approach. This approach enables us to improve understanding of the mechanisms of trust building by identifying two modes of construction. The first mode shows that a relation based on trust can be built through a "linear" process when this relation is not marred by incidents or when incidents are correctly managed if they do occur. A second method of trust construction "by stages" characterizes relations that suffer from significant incidents which were not satisfactory resolved by the bank. In this case, confidence is put to the test and may be progressively renegotiated by the partners in their course of their ongoing relationship. The inclusion of the time dimension in our analysis also makes it possible to reconcile some contradictory results arising in the literature in relation to satisfaction now considered as an antecedent or consequence of trust.

Lastly, the study highlights how institutional trust is transferred into interpersonal trust in the banking relation, depending on the age and the banking history of the customers.

The comprehension of the modes of trust building and the mechanisms of transfer of institutional to interpersonal trust are key levers for improving the retention and loyalty of the various customer segments. This result may initiate reflection on a qualitative segmentation of customers in the banking sector, including the nature and the history of the relation. In regard to the primacy of interpersonal trust, another implication of our research relates to the management of mobility of personnel in contact with customers. The current strategy of French banks in relation to staff turnover does not enable the role of advisers in the construction and reinforcement of the trust relationship to be fully taken advantage of.

Our research evidently has its limitations. One limitation relates to the bias of hindsight which may effect the respondents' accounts. Ideally, data collection would be based on observing customer-adviser behavior in real time over a succession of interactions. Further, our customer sampling criteria undoubtedly drew upon individuals with profiles that were able to cast greater light on our understanding of the trust relation.

  , durant la relation client/banque. Une étude qualitative et longitudinale a été entreprise, basée sur des études de cas rétrospectives, conduite dans une banque française. Cette recherche met en évidence le processus général de développement de la confiance interpersonnelle et l'existence de deux modes de construction de la confiance -linéaire et par paliers -selon l'occurrence d'incidents, de leur nature et de leur résolution par la banque. Nos résultats montrent aussi l'avènement d'un basculement de la confiance institutionnelle vers la confiance interpersonnelle au cours de la relation client/banque. Mots clefs : relation client/banque, processus de construction de la confiance, confiance interpersonnelle et institutionnelle, études de cas rétrospectives.


Frédéric Golay

Philippe Helluy
email: helluy@isitv.univ-tln.fr1998

Numerical simulation of a viscous compressible fluid based on a splitting method

In this paper, a new approximation to solve the equations of a stationary, viscous and slightly compressible fluid is proposed. This method is based on a splitting of the equation into an elliptic and a hyperbolic part.

Introduction

In several recent papers, it appears that great progress has been made on the understanding of the compressible Navier-Stokes equations. For example, for the system -µ∆v -(λ + µ)∇div v + ∇(Kρ γ ) = -div(ρv ⊗ v) + ρg [START_REF] Padula | l p approach to steady flows of viscous compressible fluids in exterior domains[END_REF] div(ρv) = 0 inΩ (2)

v = 0 on∂Ω (3)
which corresponds to a 2D isentropic flow with speed v and density ρ, A.

Novotný in [START_REF] Novotný | Some remarks to the compactness of steady compressible isentropic navier-stokes equations via the decomposition method[END_REF] was able to prove the existence of a solution in W 1,2 (Ω) × L γ (Ω) when the density force g is in L ∞ and γ > 1. One of the main tool in his proof is a decomposition of the operator. The same idea is used in an another framework by Novotný and M. Padula in [START_REF] Padula | l p approach to steady flows of viscous compressible fluids in exterior domains[END_REF]. In this approach, the first step is to express the speed with the help of the Helmholtz decomposition, which reads:

v = u + ∇ϕ div u = 0
The system (3) is then decomposed into three parts:

• A Stokes like system for the divergence-free part of the speed u and the effective pressure P = ρ γ -(λ + 2µ)divv

• A Neumann problem for the potential part of the speed ϕ

• A non-linear transport equation for the density ρ A crucial point in the proof is that the effective pressure is more regular than ρ γ and divv separately.

Recently, M. Padula has discovered a more natural decomposition of the equations. This new splitting permitted her to sketch a simplified existence and uniqueness result for the isothermal flow (which corresponds to γ = 1 or to a linearization of the isentropic pressure law around ρ = ρ 0 ), under the assumption of the smallness of the data. Her work is presented in [START_REF] Padula | Recent contributions and open questions in the mathematical theory of viscous compressible fluids[END_REF]. Roughly speaking, the system (3) is now split into two parts:

• A Stokes like system for the speed v and the effective pressure P = ρ -(λ + 2µ)divv

• A linear transport equation for the density ρ

An important difficulty in studying the compressible Navier-Stokes system (3) is that it is neither elliptic nor hyperbolic. With this splitting, elliptic and hyperbolic features are now clearly separated, and the theoretical study is simplified.

In this paper, we intend to show that this splitting can also be very fruitful numerically speaking. Indeed, the Stokes and the incompressible Navier-Stokes systems have been extensively studied. Building numerical approximations of these systems is not easy because of the possibility of spurious numerical solutions, but robust methods now exist (for a description of the theory, see Brezzi-Fortin [START_REF] Fortin | Mixed and hybrid finite element methods[END_REF], Bernardi-Maday [START_REF] Maday | Approximation spectrale des problèmes aux limites elliptiques[END_REF]...). On the other hand, many schemes have been developed for the numerical treatment of the transport equations.

When numerically solving the compressible Navier-Stokes equations two independent dimensionless numbers play a crucial role: the Mach number M , which measures the compressibility effects, and the Reynolds number R, which measures the viscosity effects. Even with a small Reynolds number (big viscosity), in some physical configurations, the Mach number may be big. In that case the transport phenomena can not be damped by the viscosity, because there is no second derivative of ρ in the system (3)). In the numerical treatment, oscillations may appear without a carefully chosen approximation. The reverse situation is analog: if the Mach number is small, we are almost in the case of incompressible fluid. Therefore, the approximate density may present the same pathologies than the pressure in the Stokes problem.

In that context, the splitting permitted us to apply a very natural numerical scheme to solve the compressible Navier-Stokes equations, considered as a coupled problem of two systems of partial differential equations. We shall also present an iterative algorithm to numerically solve this coupled problem.

For the approximation we decided to use two very classical and robust finite element methods:

• The quadratic Crouzeix-Raviart element (also known as the P2-bubble + P1 element) for the approximation of the Stokes-like problem.

• The Streamline Upwind Petrov Galerkin (SUPG) method of Hughes for the transport problem.

Of course many other approximation techniques may have been used (like spectral elements for the Stokes problem or high order finite volume for the transport problem) with the same splitting of the equations.

The paper is organized as follow:

• We first rewrite the equations under a dimensionless form. Thus, the assumption of smallness of g is equivalent to an assumption of smallness of some dimensionless quantities. We are then able to show how the splitting can be used to derive an existence and uniqueness result when g is small enough and γ = 1 (isothermal fluid). This result is sketched in the paper of Padula [START_REF] Padula | Recent contributions and open questions in the mathematical theory of viscous compressible fluids[END_REF], but we found it interesting to provide a detailed exposition. Indeed, some ideas of the proof are new, and we are convinced that they will also help the numerical analysis of compressible flows. It must be pointed out here that the global existence of a solution is still an open question in the case of isothermal flow.

• The technique of the previous proof is then adapted to construct an iterative algorithm to solve [START_REF] Brézis | Analyse fonctionnelle[END_REF]. It is shown to be convergent when g is small enough.

• Then, an approximation scheme, adapted to compressible flows, and numerical results are presented as an application of the previous theory. We did not perform error estimate or convergence proof for our approximation. Anyway, we present convincing numerical experiments which prove that our approximation scheme has good precision properties.

2 Dimensionless form of the equations

Equations

We present a model of viscous compressible fluid, based on an isothermal assumption. This allows us to present our method in a simple way. But it can be extended to more realistic fluids and boundary conditions. We consider a stationary, viscous, compressible and isothermal fluid, in a bounded open set Ω of R 2 or R 3 . The speed v and the density ρ satisfy the Navier-Stokes equations (see [START_REF] Lifchitz | Mécanique des fluides[END_REF]):

-µ∆v -(λ + µ)∇div v + ∇(Kρ γ ) = -div(ρv ⊗ v) + ρg (4) div(ρv) = 0 inΩ (5) v = 0 on∂Ω (6)
K, λ and µ are constant, > 0, and g is a vector field. The pressure law is of the form p = Kρ γ and here, γ = 1 because the fluid is isothermal.

Dimensionless equations

Our goal is to prove an existence and uniqueness result on [START_REF] Pitkäranta | Finite element methods for linear hyperbolic problems[END_REF] based on an assumption of smallness of g. But in the proof, it is equivalent and more elegant to use an assumption of smallness of some dimensionless constants such as the Reynolds number or the Mach number. Therefore we shall first write (6) in dimensionless form. We consider the following physical dimensions

• [L] : length unit • [T ] : time unit • [M ] : mass unit
The dimensions of the physical parameters units (6) are:

x spatial position [L] µ dynamic viscosity [M ][L] -1 [T ] -1 λ second viscosity [M ][L] -1 [T ] -1 K pressure law coeficient [L] 2 [T ] -2 v speed [L][T ] -1 ρ mass density [M ][L] -3 g force density [L][T ] -2
Let L be a characteristic length of our problem (i.e., the diameter of Ω), G a characteristic force density (i.e. the mean value of |g| on Ω, or any other norm of g), and ρ 0 , the mean value of the density on Ω. We then set x = Lx :

∂ ∂x = 1 L ∂ ∂x
And by taking V = √ GL, which is a characteristic speed of the problem, with v = v/V , we have

-µ ∆ v -(λ + µ )∇ div v + ∇ (K σ ) = -div (ρ v ⊗ v ) + (1 + σ )g (7) div (ρ v ) = 0 in Ω (8) v = 0 on ∂Ω (9)
where

• µ = µ ρ0V L • ρ = ρ ρ0 • λ = λ ρ0V L • K = K V 2 • g = g/G
• Ω is the image of Ω in the transformation x = x/L.

The Reynolds number is here

R = ρ 0 V L µ
and the Mach number is

M = V √ K
When all the parameters are fixed, it is clear that the smallness of g is equivalent to the smallness of R and M . In order to simplify the notations we will omit the primes in the sequel. Setting α = λ/µ, we obtain

- 1 R (∆v + (1 + α)∇div v) + 1 M 2 ∇ρ = -div(ρv ⊗ v) + ρg (10) div(ρv) = 0 inΩ (11) v = 0 on∂Ω ( 12
)
3 Theoretical study 3.1 Preliminary

Notations

Ω is an open set of R d (d = 2 or 3). We denote by ∂Ω its boundary, which is supposed to be of regularity

C 3 . For 1 ≤ p ≤ ∞, L p (Ω) is the usual Banach space of functions u defined on Ω such that |u| p is summable. ||u|| 0,p = Ω |u| p 1 p
For k > 0, W k,p is the usual Sobolev space of functions in L p (Ω) such that all the derivatives up to order k are in L p (Ω). If k is an integer, the norm in W k,p is:

||u|| k,p =   k i=0 |α|=i Ω |D α u| p   1 p Here, α = (α 1 • • • α d ) ∈ N d , |α| = α 1 + • • • + α d , and the notation D α u means ∂ |α| u ∂x α1 1 • • • ∂x α d d
For k ≥ 1 and p ≥ 2, the trace on ∂Ω of elements of W k,p is well defined and we can set

W k,p 0 = u ∈ W k,p , u = 0 on∂Ω
we also define

W k,p = u ∈ W k,p , Ω u = 0
If u is a vector valued function and if each component of u is in a space E, we shall write u ∈ E, rather than u ∈ E d . The Sobolev injection theorems will be frequently applied (see, for example, Brézis [START_REF] Brézis | Analyse fonctionnelle[END_REF], Adams [START_REF] Adams | Sobolev Spaces[END_REF]...).

C will be a generic constant, > 0, which depends only on Ω.

Here, our study will be restricted to rather strong solutions of ( 6). More precisely, we shall assume that g is in the Sobolev space W 1,2 . We look for a solution such that

• the speed v is in W 3,2 0
(thus, it is in C 1 (Ω) thanks to the Sobolev injection theorem and because the space dimension d ≤ 3)

• the density ρ is in W 2,2 (thus σ ∈ C(Ω))
Our approach can be extended to less regular solutions, i.e. with v in W 2,p 0 , p > 2 (see [START_REF] Padula | l p approach to steady flows of viscous compressible fluids in exterior domains[END_REF]).

Some well known results

Let's recall first some regularity properties that will be used below.

Stokes problem

The Stokes problem reads

-∆u + ∇P = F (13) div u = g (14)
u = 0 on ∂Ω (15)
it is well set for F and g in L 2 (Ω), and (u, P ) in W 1,2 × L 2 . Moreover, thanks to the ellipticity of ( 15) and the regularity of Ω, the following estimates hold (they are proved for example in Galdi [START_REF] Galdi | An introduction to the mathematical theory of the Navier-Stokes equations[END_REF]):

u 3,2 + P 2,2 ≤ C F 1,2 + g 2,2 u 2,2 + P 1,2 ≤ C F 0,2 + g 1,2
Using the relation

∇ ∧ ∇∧ = -∆ + ∇∇•
this estimate is still true for the boundary value problem

∇ ∧ ∇ ∧ u + ∇P = F (16) div u = g ( 17
) Ω P = Ω g = 0 (18) u = 0 on ∂Ω (19)
Transport equation Let v be a regular vector field on Ω (at least in W 3,2 ). We suppose that v • n = 0 on the boundary and that v is small in W 3,2 norm. The transport problem σ + div (σv) = P (20)

is well set for P in L 2 (Ω) and σ in L 2 (Ω) (see Lax-Philips [START_REF] Philips | Local boundary conditions for dissipative symmetric linear differential operators[END_REF], Rauch [START_REF] Rauch | Symmetric positive systems with boundary characteristic of constant multiplicity[END_REF]). The following regularity estimates are proved by Novotný in [START_REF] Novotný | About the steady transport equation: l p -approach in domains with smooth boundaries[END_REF] (for Ω regular)

σ k,2 ≤ C P k,2 div (σv) k,2 ≤ C P k,2 0 ≤ k ≤ 2
Taking the derivatives of (20), these inequalities are easy to prove for regular solutions (but the construction of regular sequences converging to the weak solution of (20) is quite technical).

The splitting method

splitting of the density

Because it is easier to work with functions of zero mean value, we first set ρ = 1 + σ then, Ω σ = 0. σ will be called in the sequel the perturbation of density.

splitting into a linear and a non linear problems

Let F (σ, v) = -div((1 + σ)v ⊗ v) + (1 + σ)g (21)
we then have

- 1 R (∆v + (1 + α)∇div v) + 1 M 2 ∇σ = F (σ, v) div(v) + div(σv) = 0 inΩ v = 0 on∂Ω
It is then natural to study, for fixed w and G, the linear problem

- 1 R (∆v + (1 + α)∇div v) + 1 M 2 ∇σ = G (22) div(v) + div(σw) = 0 inΩ (23) v = 0 on∂Ω (24)
The resolution of this problem defines a linear mapping Λ w : G -→ (σ, v). The solutions of ( 12) are then the fixed points of the non linear mapping (τ, w)

-→ G = F (τ, w)) -→ Λ w (G) = (σ, v).

Study of the linear problem splitting

We use here the approach initiated by Padula in [START_REF] Padula | Recent contributions and open questions in the mathematical theory of viscous compressible fluids[END_REF]. The difficulty is that the linear problem ( 24) is neither elliptic nor hyperbolic. We shall split it into two easier problems: a standard Stokes problem on v and a transport problem on σ. The coupling of these two problems will be achieved through a new unknown

P = R M 2 σ -(2 + α)divv.
We want to study

- 1 R (∆v + (1 + α)∇div v) + 1 M 2 ∇σ = G div(v) + div(σw) = 0 inΩ v = 0 on∂Ω
w is a vector field, such that w = 0 on the boundary of Ω. Moreover, if v is a solution of ( 6), it is easy to check that div v = 0 on ∂Ω. Thus, we suppose that it is also true for w, because at the end of the fixed point procedure, we must have w = v.

Using the relation

∇ ∧ ∇∧ = -∆ + ∇∇•, we can write 1 R (∇ ∧ ∇ ∧ v + ∇P ) = G div v = -div(σw) inΩ v = 0 on∂Ω Ω P = 0
Our second equation is then

σ + div σ M 2 R (2 + α)w = M 2 R P (25)
We have seen that it is natural to set w = 0 and div w = 0 on ∂Ω. Therefore, we define the following Banach space

W 3,2 0,0 = w ∈ W 3,2 0 , w = 0 et div w = 0 on ∂Ω For w in W 3,2 0,0 , we also introduce W 2,2 w = σ ∈ W 2,2 , div (σw) ∈ W 2,2
It is a Banach space, with the norm

σ 2,2,w = σ 2,2 + div σw 2,2
or with the equivalent norm (because div σw = 0 on ∂Ω)

σ 2,2,w = σ 2,2 + ∆div σw 0,2 (26)
The following result can then be stated

Theorem 1 there exists a constant β > 0, which depends only on Ω, such that if M 2 /R < β then problem (24) admits one and only one solution for G in W 1,2 (Ω) and w in W 3,2 0,0 , w 3,2 < 1. Moreover the following estimates hold

v 3,2 ≤ CR G 1,2 σ 2,2 ≤ CR G 1,2
Proof: Define the map L : σ → σ , by

• v and P are solutions of 1 R (∇ ∧ ∇ ∧ v + ∇P ) = G (27) div v = -div(σw) inΩ (28) v = 0 on∂Ω (29)
Ω P = 0 ( 30
)
• σ is solution of

σ + div σ M 2 R (2 + α)w = M 2 R P (31)
Obviously, if (v, σ) is a solution of (24), then it also verifies (30) and (31). Conversely, if σ is a fixed point of L, and if v is given by (30), then (v, σ) is solution of (24).

Then we have to find the fixed points of L. More precisely, we are going to prove that if M 2 R w is small enough, L is actually a contracting mapping from W 2,2 w into itself. This will prove the existence and uniqueness of σ in that space.

• L is from W 2,2 w into W 2,2 w : Indeed, if σ is in W 2,2 w , then div(σw) is in W 2,2
. Moreover, G belonging to W 1,2 , after resolution of the Stokes problem, v is in W 3,2 and P in W 2,2 . Then,

σ is in W 2,2
w .

• L is a contraction L(σ 1 ) -L(σ 2 ) 2,2,w should be evaluated as a function of σ 1 -σ 2 2,2,w , where σ 1 and σ 2 are two elements of W 2,2 w . v i , P i , and σ i , are the notations for the quantities associated to σ i (i = 1, 2) after resolution of (30) and (31).

From (30), we get

v 1 -v 2 3,2 + P 1 -P 2 2,2 ≤ C σ 1 -σ 2 2,2,w and ∆(P 1 -P 2 ) = 0
The resolution of the transport problem (31) gives

σ 1 -σ 2 2,2 ≤ C M 2 R P 1 -P 2 2,2 div (σ 1 -σ 2 )w 2,2 ≤ C P 1 -P 2 2,2
Those two inequalities hold under the condition that M 2 R w is small enough. However, they are not sufficient to prove the contraction. Actually they lead to

σ 1 -σ 2 2,2,w ≤ C(1 + M 2 R ) σ 1 -σ 2 2,2,w
and we cannot assume that the constant C, which depends on Ω, is small enough.

To overcome this difficulty, let us take the laplacian of the transport equation

∆(σ 1 -σ 2 ) + ∆div (σ 1 -σ 2 ) M 2 R (2 + α)w = M 2 R ∆(P 1 -P 2 ) = 0
According to the Leibniz formula, we have

∆(σ 1 -σ 2 ) + div ∆(σ 1 -σ 2 ) M 2 R (2 + α)w = M 2 R r
where r is a sum of products of derivatives of (σ 1 -σ 2 ) up to order 2 and derivatives of w up to order 3. It reads as a transport equation for ∆(σ 1 -σ 2 ). Using the regularity estimates, and the Sobolev inequalities, we get

div ∆(σ 1 -σ 2 )w 0,2 ≤ C r 0,2 ≤ C σ 1 -σ 2 2,2 w 3,2
On the other hand, div ∆(σ 1 -σ 2 )w = ∆div (σ 1 -σ 2 )w + r , where, as before, r is a sum of products of derivatives of (σ 1 -σ 2 ) up to order 2 and derivatives of w up to order 3. Hence, we obtain

σ 1 -σ 2 2,2 + ∆div (σ 1 -σ 2 )w 0,2 ≤ C M 2 R 1 + w 3,2 σ 1 -σ 2 2,2,w
and, by the norm equivalence (26),

L(σ 1 ) -L(σ 2 ) 2,2,w ≤ C M 2 R (1 + w 3,2 ) σ 1 -σ 2 2,2,w L is then a contraction in W 2,2 w , for M 2 R (1 + w 3,2 ) small enough, for instance if w 3,2 ≤ 1 and M 2 R ≤ β = 1 2C • A posteriori estimate
Let σ be the unique fixed point of L. We have

σ 2,2,w -L(0) 2,2,w ≤ σ -L(0) 2,2,w = L(σ) -L(0) 2,2,w ≤ C M 2 R σ 2,2,w moreover, σ 0 = L(0) is solution of σ 0 + div σ 0 M 2 R (2 + α)w = M 2 R P 0 then σ 0 2,2,w ≤ C P 0 2,2
and P 0 is solution of (30) with σ = 0, thus

P 0 2,2 ≤ CR G 1,2
and for M 2 R small enough, we get

σ 2,2,w ≤ CR G 1,2
Finally, we deduce, thanks to (30)

v 3,2 ≤ CR G 1,2
and the proof is completed.2

Study of the global problem

We are now in a position to give our main result for the problem (6).

Theorem 2 Let B be the ball

B = (τ, w) ∈ W 2,2 × W 3,2 0,0 , τ 2,2 + w 3,2 ≤ 1
there exists a constant β depending only on Ω such that if R + M 2 R ≤ β, then problem (6) admits one and only one solution in B.

Proof :

We consider the mapping H : (τ, w) → (σ, v), defined as follow :

• Following definition (21), we set G = F (τ, w)

• then (σ, v) is the solution of the previous linear problem (24).

Each solution (σ, v) of ( 6) is a fixed point of H. And each fixed point of H is a solution of [START_REF] Pitkäranta | Finite element methods for linear hyperbolic problems[END_REF].

• Let us show that H(B) ⊂ B If M 2
R is small enough, we know that H is well defined, and that

σ 2,2 + v 3,2 ≤ CR F (τ, w) 1,2
Sobolev inequalities give

F (τ, w) 1,2 ≤ C(1 + τ 2,2 )( w 2 3,2 + g 1,2 ) thus, for R small enough, τ 2,2 + w 3,2 ≤ 1 =⇒ σ 2,2 + v 3,2 ≤ 1 • Let us show that H is a contraction from B to B
With obvious notations, we have

σ 1 -σ 2 2,2 + v 1 -v 2 3,2 ≤ CR F (τ 1 , w 1 ) -F (τ 2 , w 2 ) 1,2
moreover, if (τ 1 , w 1 ) and (τ 2 , w 2 ) are in B, then the Sobolev injections give

F (τ 1 , w 1 ) -F (τ 2 , w 2 ) 1,2 ≤ C (1 + g 1,2 ) τ 1 -τ 2 2,2 + w 1 -w 2 3,2
for R small enough, H is thus a contraction from B to B.2

Iterative algorithm

With the aid of the previous theory, an algorithm of resolution of (6) will now be given. The construction we have described to prove theorem (2) is actually made up of two interacting iterative procedures. For practical reasons, we have to modify it in order to build a simpler algorithm.

Presentation

Let E be the set

E = (v, σ) ∈ W 3,2 0,0 × W 2,2 , div (σv) ∈ W 2,2 (32)
(E is not a linear space) We consider the following mapping

Γ : E -→ E (v, σ) -→ (v , σ ) (33)
defined as follow

• (v , P ) ∈ W 3,2 0,0 × W 2,2 is the solution of 1 R (∇ ∧ ∇ ∧ v + ∇P ) = F (σ, v) (34) div v = -div (σv) in Ω (35) v = 0 on∂Ω ( 36
) Ω P = 0 (37) • σ ∈ W 2,2 is the solution of σ + div σ M 2 R (2 + α)v = M 2 R P (38)
A fixed point of Γ is a solution of ( 6) and every solution of ( 6) is a fixed point of Γ.

Remarks

The advantages of the algorithm (37), (38) are:

• (37) and ( 38) are linear problems

• the theory of each problem is well known and a great amount of numerical methods have been developed to solve it.

• because it is decomposed into several problems, ( 6) can be solved, even if it is too huge to stand alone in memory.

• only the transport linear system has to be factorized at each iteration

The drawbacks are:

• the convergence may be slow and we cannot expect more than linear convergence.

• the algorithm is divergent if the Mach number or the Reynolds number are too big.

This algorithm should thus be considered as a first step to a faster and more efficient resolution algorithm. Actually, we decided to adopt it because it permitted us very easily to reuse already written finite element software.

Stability

We show here the stability of our iterative scheme Theorem 3 Let η > 0, be small enough. Let E be defined as in (32). Suppose that

g 1,2 ≤ η Let (v, σ) ∈ E, such that v 3,2 + σ 2,2 ≤ η div(σv) 2,2 ≤ η If we set (v , σ ) = Γ(v, σ) (c.f. definition (33)), then there exists β > 0, such that for R + M 2 R < β, we also have v 3,2 + σ 2,2 ≤ η div(σ v ) 2,2 ≤ η
and the constants β and η depend only on Ω.

proof this result is the consequence of the following estimates (already proved before)

• F (σ, v) 1,2 ≤ C(1 + σ 2,2 )( v 2 3,2 + g 1,2 ) • v 3,2 + P 2,2 ≤ C R F (σ, v) 1,2 + div (σv) 2,2 • ∆P 0,2 ≤ CR F (σ, v) 1,2 • σ 2,2 ≤ C M 2 R P 2,2 • ∆div (σ v ) 0,2 ≤ C ∆P 0,2 + σ 2,2 v 3,2 .2

Convergence

Theorem 4 Let η > 0, be small enough. Let E be defined as in (32). Suppose that

1,2 ≤ η Let (v 0 , σ 0 ) ∈ E, such that v 0 3,2 + σ 0 2,2 ≤ η div(σ 0 v 0 ) 2,2 ≤ η If we set (v n+1 , σ n+1 ) = Γ(v n , σ n ) (c.
f. definition (33), then there exists β > 0,

such that for R + M 2 R < β, the sequence (v n , σ n ) converges in W 2,2 0,0 × W 1,2
towards the unique solution in E (noted (v, σ)) of ( 6).

proof Denote by (v, σ) the unique solution of ( 6) in E. Existence and uniqueness of that solution is given by theorem 2. P = R M 2 σ+div (σ(2 + α)v) is the effective pressure related to that solution.

We then have

1 R (∇ ∧ ∇ ∧ (v n+1 -v) + ∇(P n+1 -P )) = F (σ n , v n ) -F (σ, v) div (v n+1 -v) = -div (σ n v n ) + div (σv) inΩ v n+1 -v = 0 on∂Ω Ω P n+1 -P = 0
therefore, using the previous stability result and the regularity of the Stokes problem, we deduce

v n+1 -v 2,2 + P n+1 -P 1,2 ≤ CR σ n -σ 1,2 + v n -v 2,2 +Cη σ n -σ 1,2 + v n -v 2,2
It is worthy to notice here that because of the term -div (σ n v n )+div (σv) we lost an order of derivative in the estimate. A stability result in a stronger norm gives us a convergence result in a coarser norm.

We then consider the transport equation

σ -σ n+1 + div (σ -σ n+1 ) M 2 R (2 + α)v = M 2 R (P n+1 -P ) + div σ n+1 M 2 R (2 + α)(v n+1 -v) thus σ n+1 -σ 1,2 ≤ C M 2 R P n+1 -P 1,2 + σ n+1 2,2 v n+1 -v 2,2
and for R + M 2 R + η small enough, the algorithm is convergent.2

4 Numerical study

Choice of the approximations

We shall now describe the finite element methods used to solve (38) and (37).

The speed and the pressure of the Stokes problem and the density of the transport problem will be computed on the same mesh (with the same numbering of the elements but not necessarily of the vertices). In this way, the quantities of a problem are easily accessed by the other problem.

Stokes problem

The classical quadratic finite element of Crouzeix-Raviart (described in Cuvelier [START_REF] Steenhoven | Finite element methods and Navier-Stokes equations[END_REF], Crouzeix-Raviart [START_REF] Raviart | Conforming and nonconforming finite element methods for solving the stationnary stokes equations[END_REF], Pironneau [START_REF] Pironneau | Méthode des éléments finis pour les fluides[END_REF]...), and the following weak formulation are used: To solve 1 (for u ∈ W 1,2 , P ∈ L 2 ) It must be remembered that in the case of a triangular quadrangular element, the unknowns are the two components of the speed at the vertices and the centroid of the element, and the pressure and its gradient at the centroid. The approximation is thus quadratic and continuous for the speed, and linear and discontinuous for the pressure. The precision is then quadratic for the speed (in the W 1,2 norm) and for the pressure (in the L 2 norm), when the solution is regular enough. Stability and convergence of this approximation are studied for example in Cuvelier [START_REF] Steenhoven | Finite element methods and Navier-Stokes equations[END_REF], Pironneau [START_REF] Pironneau | Méthode des éléments finis pour les fluides[END_REF], Crouzeix-Raviart [START_REF] Raviart | Conforming and nonconforming finite element methods for solving the stationnary stokes equations[END_REF], Brezzi-Fortin [START_REF] Fortin | Mixed and hybrid finite element methods[END_REF], Bernardi-Maday [START_REF] Maday | Approximation spectrale des problèmes aux limites elliptiques[END_REF]... The size of the associated linear system can be reduced because the unknowns at the centroids of the elements can be eliminated (2 speed components and 3 pressure unknowns) without loss of precision. Those unknowns can be recovered a posteriori.

-∆u + ∇P = F div u = h inΩ u =
Moreover, this finite element method can obviously be adapted to the nonlinear incompressible Navier-Stokes system. A non-linear system has to be solved (for example with the Newton method).

Transport equation

A weak formulation of Hughes [START_REF] Brooks | Streamline upwind petrov-galerkin formulation for convection dominated flow with particular emphasis on the incompressible navier-stokes equations[END_REF], the Streamline Upwind Petrov Galerkin (or SUPG in short), is used to solve the transport problem, with quadratic quadrangular finite elements. For σ ∈ W 1,2 , P ∈ W 1,2 , the general transport equation reads:

σ + div (σw) = P σ = σ 0 on Γ - Γ -= {x ∈ ∂Ω, n(x) • w(x) < 0} Let B be the bilinear form on W 1,2 × W 1,2 B(σ, τ ) = Ω (σ + div (σw)) (τ + ε∇τ • w) - Γ - στ (w • n) the formulation is then B(σ, τ ) = Ω P (τ + ε∇τ • w) - Γ - σ 0 τ (w • n) ∀τ ∈ W 1,2
ε is a small parameter which stabilizes the approximation. It should be of order h, where h is the maximal diameter of the finite elements (for a theoretical study, see Johnson & Al. [START_REF] Pitkäranta | Finite element methods for linear hyperbolic problems[END_REF]). For a quadratic finite element basis, the error can be proved to be of order h 5/2 in a mesh-dependant norm, intermediate between the L 2 and the H 1 norm. Numerically, the error behaves like h 3 in the L 2 norm, and like h 2 in the H 1 norm.

Numerical results

Numerical convergence study

Our first test is devoted to a numerical precision study. We shall try to guess the convergence orders for the speed and the density.

Consider the two disks

D 1 = (x, y) ∈ R 2 , x 2 + y 2 ≤ 1 D 2 = (x, y) ∈ R 2 , x 2 + y 2 ≤ 4 and set Ω = D 2 \ D 1 .
The following speed vector field, and density are then considered

v = 2 3 (r - 1 r ) e θ (39)
ρ = 1 + σ = Ar -2a 2 /K exp a 2 2K (r 2 - 1 r 2 ) ( 40
)
where

r 2 = x 2 + y 2 e r = 1 r x y e θ = 1 r -y x a = 2/3 K = 10 = 1/M 2
The constant A is given by the relation

Ω ρ = 1
We found A 0.992426... it is then easy to check that (σ, v) is a solution to

∆v + ∇div v + 1 M 2 ∇σ = -div(ρv ⊗ v) (41) div((1 + σ)v) = 0 inΩ (42) Ω σ = 0 (43) v = e θ on ∂D 2 (44) v = 0 on ∂D 1 (45)
A numerical resolution of (45) as been performed for 1/M 2 = 10 and for several mesh stepsizes h. The following table (see also figure 4.2.3) gives the H 1 and the L 2 norm of the error for the speed and the perturbation of density σ = ρ -1.

h H 1 speed error L 2 density error H 1 density error 1 0,00789573 0,000397334 0,012118 1/2 0,0020248 0,0000979495 0,0060295 1/4 0,000490715 0,0000243944 0,00306567

The error curves in a logarithmic scale are almost lines, and their slopes gives the approximation order H 1 speed order L 2 density order H 1 density order 2,04482213558421 2,00548812857888 0,975835962359618 It thus appears that the precision of our approximation is the same as for the incompressible flow. There is no loss of precision due to the transport phenomena. It is also an indication that we could have use linear finite element to solve the transport problem, with the same precision.

Numerical stability study

The stability of the iterative algorithm is now studied on the previous problem (45).

The algorithm (37), (38) has been tested for different values of M . In order to obtain faster convergence, the described algorithm must be slightly modified. We implicit, at each iteration, the right hand side of (37), and actually solve This modification allows us to achieve convergence with fewer iterations. Of course, each iteration requires more CPU time because the Jacobian matrix in the Newton method has to be factorized. This approach is used in our test, but keeping a factorization for several iterations is certainly more efficient.

The results are presented on figure 4.2.3 (where, as above, K = 1/M 2 ). We have plotted the discrete quadratic norm of σ n+1 -σ n divided by the discrete quadratic norm of σ 1 -σ 0 during the iterations.

We notice that the algorithm is divergent for M 2 ≥ 1/20 The convergence (or the divergence) is obviously linear and the logarithm of the contracting constants is given by the slope of the linear part of the curves. Actually, if we suppose that

σ n+1 -σ n = (C + ε n ) σ n -σ n-1
where ε n is a sequence which tends rapidly to zero, it is easy to prove that, asymptotically log σ n+1 -σ n σ 1 -σ 0 = n log C + C 0 thus, we are able to plot log C against the constant K. We obtain the plot of figure 4.2.3. The points are very close to a line of slope -1. It is a numerical proof that the contracting constant is of the form C = AM 2 (when the Reynolds number is fixed). This result is in accordance with the previous theory where we have seen the importance of the ratio M 2 R . It is also an indication that theorem 4 could probably be improved. In that theorem we were only able to prove that the contracting constant is of the form C(R + M 2 R + η). On the other hand, a very simple relaxation method permits us to make the iterative algorithm converge up to M 2 = 2: instead of Γ, we consider the following mapping Γ : • σ is solution of

E -→ E (v,
σ + div σ M 2 R (2 + α)v = M 2 R P
• Finally, we set σ = (1 -ω)σ + ωσ, with 0 < ω < 1 and ω close to 1 when the Mach number increases.

We obtain the results of figure 4.2.3 (with ω = 0.9).

1ΩF w ∀w ∈ W 1 L 2 u

 12 Of course, the following compatibility relation is supposed to hold (ν is the outward normal to Ω on ∂Ω): ∂Ω u b • ν = Ω h and the following weak formulation A(u, w) + B(w, P ) = = u b on∂Ω which is restricted to the finite dimension space spanned by the finite element basis of Crouzeix-Raviart.

1 R

 1 (∇ ∧ ∇ ∧ v + ∇P ) = F (σ, v ) (rather than F (σ, v)) (46) div v = -div (σv) in Ω (47) v = e θ on ∂D 2

  σ) -→ (v , σ ) defined as follow• v and P are solutions of 1 R (∇ ∧ ∇ ∧ v + ∇P ) = F (σ, v ) div v = -div (σv) in Ωand subject to the usual boundary conditions.

Figure 1 :

 1 Figure 1: Approximations errors

Figure 2 :Figure 3 :

 23 Figure 2: Convergence of the density

An another example

In our last example, we consider

and the problem

which corresponds to a rotating elliptic cylinder in a 2D compressible fluid limited by a circular cylinder. The equations are written in the referential linked to the ellipse. It explains the presence of the term g = 1 4 r which is the centrifugal acceleration.

We fix M 2 = 1/20 and we obtain the results of figures 4.2.3, 4.2.3, 4.2.3 and 4.2.3 for different Reynolds numbers. For R = 1, we find sup σ 0.3, and inf σ -0.3. To plot, we chose a scale between -0.06 and 0.06. Indeed, the compressibility effects are localized near the extremities of the ellipse. With a larger scale, only the variations close to these extremities would have been visible. For the other values of R, our scales respect the variations of σ (which are of order 0.05).We explain our results by the fact that the shear forces due to the viscosity diminish when the Reynolds number increases. When R → ∞, the ellipse behaves as a circle of radius 1.


J.-C Le Gac
email: legac@shom.fr

Y Stephan

M Asch

P Helluy

A VARIATIONAL APPROACH FOR GEOACOUSTIC INVERSION USING ADJOINT MODELING OF A PE APPROXIMATION MODEL WITH NON LOCAL IMPEDANCE BOUNDARY CONDITIONS

The adjoint model method of control theory is known to give accurate and efficient data assimilation processes in oceanography and meteorology. However, it has rarely been applied in underwater acoustics for inversion purposes. Based on the back-propagation of the mismatch between observations and their predictions, the adjoint model can produce the corrections to the associated direct forward model input parameters. In this paper, the adjoint of a parabolic equation propagation model with non local impedance boundary conditions at the water sediment interface is used in order to determine an acoustically equivalent representation of the seabed. The bottom is represented by these boundary conditions that play the role of the control parameter.

Introduction

Bottom properties are essential in shallow water acoustics for prediction of the transmission losses encountered in sonar applications. In the last decade inversion methods based on signal processing approaches (matched field processing, inversion of broadband signals) have been applied extensively to estimate geoacoustic features of the bottom (see [START_REF] Chapman | Benchmarking geoacoustic inversion methods -Special Issue[END_REF][START_REF] Hermand | Broad-band geoacoustic inversion in shallow water from waveguide impulse response measurements on a single hydrophone: Theory and experimental results[END_REF] for a survey of the state of the art). The optimization techniques that have so far been used within these methods are mostly perturbative approaches based on the iterative adjustment of the input parameters of a forward model. The forward model is run for a great number of input parameter values until the mismatch between the observations (like the spatial structure of the complex pressure field or the temporal structure of the propagated signal) and their predictions is low. This often results in highly computationally intensive methods.

In this paper, we present an alternative technique based on a variational approach coming from the control theory. This technique aims to solve optimization problems based on partial differential equations. The theory is detailed in [START_REF] Lions | Contrôle optimal des systèmes gouvernés par des équations aux dérivées partielles[END_REF][START_REF] Lions | Contrôlabilité exacte, perturbations et stabilisation des systèmes distribués[END_REF]. Variational methods by adjoint modeling are tools that have so far been developed in the framework of inverse modeling of physical systems, in particular in the domains of geophysics and molecular physics. Their use has been generalized especially to meteorology and oceanography for the purposes of data assimilation [START_REF] Talagrand | Variational assimilation of meteorological observations with the adjoint vorticity equation -Part I[END_REF][START_REF] Courtier | Variational assimilation of meteorological observations with the adjoint vorticity equation -Part II[END_REF][START_REF] Leredde | On initial boundary conditions and viscosity coefficient control for Burgers' equation[END_REF], model adjustment [START_REF] Schröter | Variational assimilation of GEOSAT data into an eddy-resolving model of the Gulf-Stream extension area[END_REF], data inversion [START_REF] Plessix | Waveform inversion of reflection seismic data for kinematic parameters by local optimization[END_REF][START_REF] Jameson | Aerodynamic design via control theory[END_REF][START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observation: Theoretical aspects[END_REF][START_REF] Belmiloudi | A control method for assimilation of surface data in a linearized Navier-Stokes-type problem related to oceanography[END_REF] and sensitivity studies [START_REF] Cacucci | Sensitivity theory for nonlinear systems. Part I: Non linear functional analysis approach[END_REF]. Variational methods have been recently used for acoustic inversion [START_REF] Asch | An adjoint method for geoacoustic inversion[END_REF][START_REF] Hursky | Adjoint-assisted inversion for shallow water environment parameters[END_REF][START_REF] Asch | Adjoint modeling for geoacoustic inversion -I: PE approximation with local impedance boundary conditions[END_REF][START_REF] Gac | Adjoint modeling for geoacoustic inversion -II: PE approximation with non local impedance boundary conditions[END_REF]. In [START_REF] Asch | An adjoint method for geoacoustic inversion[END_REF], we proposed an optimal control approach based on the use of a parabolic equation (PE) propagation model in which a local impedance condition at the water-sediment interface plays the role of the control parameter. The geoacoustic concept that was used was that of the acoustically equivalent medium. Rather than seeking the geoacoustic parameters in the physical parameter space, we computed the sea bottom conditions in a way that was not directly comparable with a ground truth model but that was sufficient to predict the field of acoustic transmission losses. These first developments and a feasibility study allowed us to demonstrate the controllability of the PE approximation suggesting the variational approach as a rather promising technique for geoacoustic inversion. However, the use of local impedance conditions is rather limited (though realistic in the framework of aerial acoustics [START_REF] Robertson | Low frequency sound propagation modeling over a locally reacting boundary with the parabolic approximation[END_REF]), since they can only be used for local reacting media. For underwater acoustics, a downgoing radiation condition must be imposed on the transmitted component of the field since the seabed is penetrable to sound waves, especially at low frequencies. In PE models, this can be done by appending an absorbing layer [START_REF] Jensen | Computational Ocean Acoustics (AIP Series in Modern Acoustics and Signal Processing[END_REF] or alternatively by applying a non local boundary condition (NLBC). This last method that was first formulated and implemented in a finite difference code by Papadakis [START_REF] Papadakis | Impedance formulation of the bottom boundary condition for the parabolic equation model in underwater acoustics[END_REF] is rather attractive. Since then, numerous different forms of the NLBC concept have been proposed (see [START_REF] Yevick | Non local boundary conditions for finite difference parabolic equation solvers[END_REF] for some references).

In this paper, we apply the NLBC formulation described in [START_REF] Yevick | Non local boundary conditions for finite difference parabolic equation solvers[END_REF] to generate realistic boundary conditions and then adopt a generalized form of these NLBC conditions for the rest of the study. In section 2, the inverse continuous physical problem is presented. We define the associated optimal control problem and compute the analytical form of the gradient of the cost function by using the adjoint model of a PE propagation model. The use of the adjoint model and the gradient of a cost function within a gradient based optimization loop is then presented. In section 3, we present some numerical simulations showing the data assimilation of acoustical data for geoacoustic inversion purposes.

The continuous problem

The variational approach being new in the framework of geoacoustic inversion, we adopt a rather simple monochromatic matched field processing inverse technique. The acoustical data we thus use (or the observations) are the spatial structure of the complex pressure field measured along a vertical line array.

The direct forward model

The direct forward model is the standard PE approximation [START_REF] Tappert | The parabolic approximation method[END_REF]. The model recognized as the Standard PE is given for a harmonic point source of time dependence For the boundary conditions at the water-sediment interface H z = , we adopt the NLBC formulation proposed by Yevick and Thomson [START_REF] Yevick | Non local boundary conditions for finite difference parabolic equation solvers[END_REF]. In that paper, the authors proposed a procedure for obtaining NLBC's directly from the z-space Crank-Nicolson solvers for both the narrow-angle standard and the wide-angle Claerbout PE's. The calculation of the wave field

) exp( t iω - in cylindrical coordinates ) , ( z r independent of the azimuthal angle θ : . 0 ) 1 ( 2 2 2 0 2 2 0 = - + ∂ ∂ + ∂ ∂ u n k z u r u ik
( ) , 0 , 0 = ∆ + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ Γ - ∂ ∂ H r r u i z b w ρ ρ (3)
where ρ is the density, the subscript b indicates the bottom, and the subscript w indicates the sea water. The numerical implementation of this NLBC requires that the operator 0 Γ is Taylor-expanded in T , which can be finally expressed as: The above model was shown to give accurate simulations of the acoustic pressure field for homogeneous semi-infinite fluid or solid half-spaces. In the following study, this reference model is used to generate realistic NLBC's related to some given geoacoustic parameters. However, due to the difficulty to directly apply the variational approach to this reference model, we chose to generalize Eq.( 4) in the following form:

[ ] [ ], , ) 1 ( , ) 1 ( 1 1 , 0 ∑ + = ∆ - + = ∆ + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ - ∂ ∂ J j j H r j J u g i H r J u i z β β
[ ]. , 0 , ) ( ) ( R r H z r F u r i z ∈ ∀ = = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ - ∂ ∂ β ( 5 )
For the numerical results presented in section 3, the initial and the true NLBC parameters β and F are first calculated with the reference model to obtain simulations that represent realistic environments of underwater acoustics. As a synthesis, the continuous direct model can be written as:

[ ] [ ] [ ] [ ] [ ] ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ ∈ ∀ = = - ∂ ∂ ∈ ∀ = = ∈ ∀ = = ∈ ∀ ∀ = - + ∂ ∂ + ∂ ∂ R r H z r F u r i z u R r z u H z r z z S u H R z r u n k z u r u ik s , 0 , ) ( ) ( , 0 , 0 0 , 0 , 0 ) , ( , 0 * , 0 , 0 ) 1 ( 2 2 2 0 2 2 0 β (6)
where R is the distance between the source and the vertical line array on which the acoustic pressure field is measured, and ) , ( s z z S is the acoustic source at s z z = .

The inverse problem formulation

The control procedure consists of finding the optimal controls ] , [ , obtained from the direct model Eq.( 6):

( ) . ) , ( , ; , 2 1 ) , ( 0 2 ∫
= = - = H z z d dz z R u z R F u F J β β (7)
The problem is to find:

[ ] { } ad G F F J ∈ , ), , ( min
β β (8)
where ad G is the set of admissible controls. We note that Eq.( 7) is the simplest possible cost function and that we can add regularization terms to it. , the concatenation of some perturbations of the control parameters. The necessary condition for the existence of a minimum is given by the following theorem.

Let us define

Theorem 1: If J attains a (local) minimum at a point ad G ∈ * ψ , then 0 ) ; ( * = ∈ ∀ φ ψ δ φ J G ad , where t J t J J t ) ( ) ( lim ) ; ( 0 ψ φ ψ φ ψ δ - + = →
is the directional Gâteaux derivative of J at the point ψ along the direction φ .

Following the same developments as in [START_REF] Asch | An adjoint method for geoacoustic inversion[END_REF], we start by calculating the variation of J :

, ) 0 ( ) , ( ) , ( ) , ; ( ) , ( ) , ; ( lim 2 1 ) ( ) ( lim ' ' 2 2 0 0 g J t z R u z R u z R u z R t u t J t J d d t t = = - - - + = - + → → φ ψ ψ φ ψ ψ φ ψ (9)
where g is the function defined by: ). ( ) (

φ ψ t J t g + = (10)
We introduce the real-valued scalar product:

) ( , v u e v u ℜ = , ( 11
)
where v is the complex conjugate of v .

Assuming that ) , ( ' φ ψ J exists for all φ and is continuous and linear inφ , the gradient ) (φ J ∇ can be defined as the linear form: . ),

J J ∈ ∀ ∇ = φ φ φ φ ψ ( 1 2 ) ( ) , ( ' ad G
The optimization problem Eq.( 8) cannot be resolved directly by using 12). The variation can be written in the following form:

[ ] [ ] [ ] [ ] [ ] ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ ∈ ∀ = = - ∂ ∂ ∈ ∀ = = ∈ ∀ = = ∈ ∀ ∀ = - + ∂ ∂ + ∂ ∂ R r H z F u i z u R r z u H z r z z S u H R z r u n k z u r u ik s , 0 , ) ( ) ( , 0 , 0 0 ) ( , 0 , 0 ) , ( ) ( , 0 * , 0 , 0 ) ( ) 1 ( ) ( ) ( 2 2 2 0 2 2 0 ψ β ψ ψ ψ ψ ψ ψ (14) and [ ] [ ] [ ] [ ] [ ] ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ ∈ ∀ = + = + + - ∂ + ∂ ∈ ∀ = = + ∈ ∀ = = + ∈ ∀ ∀ = + - + ∂ + ∂ + ∂ + ∂ R r H z tf F t u tb i z t u R r z t u H z r z z S t u H R z r t u n k z t u r t u ik s , 0 , ) ( ) ( ) ( , 0 , 0 0 ) ( , 0 , 0 ) , ( ) ( , 0 * , 0 , 0 ) ( ) 1 ( ) ( ) ( 2 2 2 0 2 2 0 φ ψ β φ ψ φ ψ φ ψ φ ψ φ ψ φ ψ ( 1 5 )
By subtracting Eq.( 15) from Eq.( 14), dividing the result by t and letting t tend to zero, the system satisfied by w , called the tangent linear model (TLM) is given by:

[ ] [ ] [ ] [ ] [ ] ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ ∈ ∀ = - - ∂ ∂ ∈ ∀ = ∈ ∀ = ∈ ∀ ∀ = - + ∂ ∂ + ∂ ∂ R r f H r u ib H r w i z H r w R r r w H z z w H R z r w n k z w r w ik , 0 ) , ( ) , ( ) , ( , 0 0 ) 0 , ( , 0 0 ) , 0 ( , 0 * , 0 , 0 ) 1 ( 2 2 2 0 2 2 0 β (16)
The aim is to use the TLM to evaluate ) , ( ' φ ψ

J

. By introducing a function ) , ( z r p , and by using the TLM, ) (φ J ∇ can be identified directly from Eq.( 12). The new variable ) , ( z r p , called the adjoint variable, satisfies the following system (that is, of the same degree of complexity as the TLM satisfied by w :

[ ] [ ] [ ] [ ] [ ] ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ ∈ ∀ = + ∂ ∂ ∈ ∀ = ∈ ∀ - = ∈ ∀ ∀ = - + ∂ ∂ + ∂ ∂ R r H r p i z H r p R r r p H z u u z R p H R z r p n k z p r p ik d , 0 0 ) , ( ) , ( , 0 0 ) 0 , ( , 0 ) , ( , 0 * , 0 , 0 ) 1 ( 2 2 2 0 2 2 0 β ( 1 7 )
This model, called the adjoint model is running backwards in r from R r = to 0 = r . It allows us to calculate ' J . Taking the scalar product of the TLM Eq.( 16) with the adjoint variable p and integrating by parts:

, 0 , ) 1 ( 2 2 2 0 2 2 0 ∫ ∫ = - - ∂ ∂ - ∂ ∂ z r drdz p w n ik z w k i r w ( 1 8 )
we finally find, after some analytical developments, that the derivative of J with respect to φ is given by: . ,

∫ ∫ = = - - = H z r H z r ip f k p u b k J φ ψ 2 1 , 2 1 ) , ( , 0 , 0 '
This latest equation can be reformulated in terms of the gradient of the cost function:

, 2 2 , 0 0 H z r k ip k p u J = ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ - - = ∇ ( 20
)
where p is the solution of the adjoint equation obtained from the application of the theorem.

The optimization procedure

By using the successive resolution of both the forward direct model and the backward adjoint model, we can therefore calculate the exact gradient of the cost function relative to the control parameter [ ]

) ( ), ( ) ( r F r r β ψ =
. Once the gradient of the cost function J ∇ is known, we can seek a (local) minimum of ) (ψ J

. The global optimization procedure is summarized in Fig. 1. In order to accelerate the convergence, we use a classical conjugate gradient method of Polack-Ribière type rather than the simplest method of steepest descent. The algorithm that is used is described in [START_REF] Frandsen | Unconstrained optimization[END_REF] and proceeds in two steps. First, the direction of minimization is computed following the Polack-Ribière algorithm that exploits the gradient of the cost function Eq. [START_REF] Tappert | The parabolic approximation method[END_REF]. Then the step length along the minimization direction is computed thanks to a soft linesearch algorithm that applies the Wolfe criteria. The interested reader is referred to [START_REF] Frandsen | Unconstrained optimization[END_REF] for a precise description of the implemented algorithm. The numerical results that are presented in this section rely on the finite-difference implementation of the direct and the adjoint models presented above, and of the associated discrete version of the gradient of the cost function. Before applying the inversion procedure, a gradient check (not shown here) was performed.

For a frequency of 500 Hz, we present the results obtained with an isovelocity sound speed profile in the water column of 1500 m s -1 . The true non local impedance boundary conditions were initially calculated for a sandy reflecting bottom a compression speed of 1600 m s -1 , an attenuation of 0.5 dB/λ and a density of 1.8 g cm -3 . The initial ones were obtained with a soft sediment (mud) characterized by a compression speed of 1505 m s -1 , the other parameters remaining unchanged. No more than 5 or 6 iterations of the optimization algorithm were needed for the algorithm to converge. However, due to the drastic convergence criteria we adopted for our calculations, we reached a maximum of 10 iterations. The corresponding computations, performed on a high-end PC (Pentium IV, 2 GHz), required around 20 minutes for a complete inversion. By optimizing the convergence criteria, this could be reduced by a factor of two.

Figure 2 gives an overall representation of the reconstruction of the field in the entire range-depth domain. Although the initial field (top) was quite different from the true field (middle), the field calculated after inversion (bottom) agrees very well with the true one. The improvement of the acoustic pressure filed is drastic, the initial error for ranges greater than 400 m being reduced from 30-40% to only 2% approximately. 3 is devoted to the representation of the parameters. We can observe that the accuracy of the F parameter is excellent after the inversion process. On the contrary, the β parameter is badly determined, which is certainly due to the ill-posed nature of the inverse problem. The use of a wider angle PE approximation is susceptible to improve our results since the angular limitation of the standard PE approximation imposes that the chosen Gaussian source has no interaction with the bottom near the source. This assumption is under investigation.

Conclusions

We have generalized the variational approach that was presented in [START_REF] Asch | An adjoint method for geoacoustic inversion[END_REF]. Here we treat the case of a PE propagation model constrained by NLBC's to simulate realistic environments in underwater acoustics. The obtained results are promising. They show that for a reasonable computational cost punctual acoustical data can be assimilated (here, vertical array data) which results in a global and drastic improvement of the prediction of complex pressure fields. The other central concept underlying our variational approach is that of the acoustically equivalent medium. This results in a rather abstract but concise representation of the geoacoustic properties of the seabed via its non local boundary conditions. These NLBC's are not easily related to the classical geoacoustic parameters of the bottom but are able to simulate the same acoustical behavior, at least for the simple but realistic environments considered in our simulations. In sonar applications, one of the main objectives is to be able to predict the acoustical losses in the water column. For this purpose the concept of acoustically equivalent medium is relevant and NLBC's are therefore good potential candidates as equivalent models. The presented method is rather new and should be studied in order to improve the generality of the approach. Further works should be involved in the use of sparse linear arrays and wider angle PE approximations or other propagation models. The application of the method will also be studied in the case of more complex seabed properties (multi-layered sub-bottoms). More importantly, and following the last tendencies in geoacoustic inversion, multi-tone and broadband signals will be considered in order to physically constrain the inverse problem and to aim at a unique inverted solution.

are

  be slowly varying with range. In order to obtain a well-posed initial boundary value problem, we add an initial pressure Dirichlet conditions (perfectly reflecting interface).

  along the water-sediment interface is obtained in terms of the known field at the previously calculated range values from 0 to r by expanding the approximate vertical wave number operator for the downgoing field in powers of the translation operator is the index of refraction of the bottom. The authors obtain the following impedance boundary condition:

  on the standard matched field processing between a measured field, the vertical line array at a given range R and the corresponding

  then have to write down the conditions satisfied by dt du w / =. In order to do this, we write the two systems satisfied by )

Figure 1 .

 1 Figure 1. Optimization algorithm -J is minimized until some convergence criteria are attained.

Figure 2 .

 2 Figure 2. Visualization of the amplitude of the initial (top), true (middle) and inverted (bottom) u fields for f=500 Hz. After the assimilation process, the inverted and the true fields are nearly identical.

Figure 3 .

 3 Figure 3. Comparison of the initial, the true and the inverted control parameters -(a)(b) Inversion of the real and imaginary parts of the β parameter -(c)(d) Inversion of the real and imaginary parts of the F parameter.

Figure

  Figure 3 is devoted to the representation of the parameters. We can observe