Rabah Rabah
email: rabah.rabah@mines-nantes.fr

Grigory M Sklyar
email: sklar@univ.szczecin.pl

Exact observability and controllability for linear neutral type systems

Keywords: Observability, Controllability, Duality, Neutral type systems

Introduction

Approximate and spectral controllability and the corresponding dual notions of observability for delay systems of neutral type were widely investigated at the end of the last century (see books by [START_REF] Bensoussan | Representation and control of infinite-dimensional systems[END_REF] and [START_REF] Salamon | Control and observation of neutral systems[END_REF] and references therein). The duality between these notions for systems of neutral type is not so trivial. The main reason is that the dual or adjoint system is not obtained directly by simple transposition. It is necessary to consider the duality using some hereditary product proposed first for retarded systems and later for neutral type systems (see [START_REF] Delfour | The structural operator F and its role in the theory of retarded systems. II[END_REF][START_REF] Delfour | The structural operator F and its role in the theory of retarded systems. I[END_REF][START_REF] Salamon | On controllability and observability of time delay systems[END_REF] and [START_REF] Salamon | Control and observation of neutral systems[END_REF] for example). In this context, the important technique of the so called structural operator was used. It enables some explicit formulations for duality between approximate controllability, spectral controllability and the same notions for observability and the characterizations of that concepts. We shall consider some of them in the context of our framework.

The infinite dimensional setting has been developed essentially for exact controllability and often for neutral type systems without distributed delays. The exact observability problem has been less studied. In [START_REF] Metel ′ Skiȋ | Criteria for the constructive identifiability and complete controllability of linear time-independent systems of neutral type[END_REF][START_REF] Khartovskiȋ | The complete controllability problem and its generalization for linear autonomous systems of neutral type[END_REF] and [START_REF] Khartovskiȋ | Complete controllability and controllability for linear autonomous systems of neutral type[END_REF] an approach is described based on the reconstruction of a part of the state for the case of a neutral type system with discrete delays. A duality condition with null controllability is given. The time of controllability (and of possible reconstruction of a part of the observed state) is estimated sufficiently large, without more precision.

The present paper is concerned with exact observability which is related to the notion of exact controllability developed in the paper of the authors [START_REF] Rabah | The analysis of exact controllability of neutral-type systems by the moment problem approach[END_REF] as an extension of other results obtained essentially for neutral type systems with discrete delay [START_REF] O'connor | On the function space controllability of linear neutral systems[END_REF][START_REF] Bartosiewicz | Approximate controllability of neutral systems with delays in control[END_REF]. The semigroup approach used by the authors in [START_REF] Rabah | The analysis of exact controllability of neutral-type systems by the moment problem approach[END_REF] is based on the model introduced in [START_REF] Burns | Linear functional-differential equations as semigroups on product spaces[END_REF] in the product space M 2 (see the definition below, in this section). In the infinite dimensional setting described in [START_REF] Rabah | The analysis of exact controllability of neutral-type systems by the moment problem approach[END_REF], exact controllability means reachability of the operator domain because reachability of all the state space is not possible by finite dimensional control. Hence, as it may be expected, the dual notion of observability is also adapted. Here, the approach using the structural operator is not used. Considering the adjoint system in the operator form, that is in an infinite dimensional framework, we construct a transposed neutral type system corresponding to the adjoint system in Hilbert space. This relation between the adjoint semigroup and the obtained neutral type system is different from that of the model given in [START_REF] Burns | Linear functional-differential equations as semigroups on product spaces[END_REF].

The notions of exact controllability and observability are important because they imply exponential stabilizability or exponential convergences for possible estimators.

The results obtained in [START_REF] Rabah | The analysis of exact controllability of neutral-type systems by the moment problem approach[END_REF] use the approach of moment problems and allows the minimal time of exact controllability to be determinated. The main contribution of our study is to specify how duality may be used in a nontrivial context and to deduce the characterization of exact observability and also the minimal time of observation.

We consider the neutral type system given by the equation

ż(t ) = A -1 ż(t -1) + 0 -1 [A 2 (θ) ż(t + θ) + A 3 (θ)z(t + θ)] dθ (1)
where A -1 is a constant n × n-matrix, and A 2 , A 3 are n × n-matrices whose elements belong to L 2 (-1, 0). If we introduce the linear operator L :

H 1 ([-1, 0]; R n ) -→ R n , defined by L f = 0 -1 A 2 (θ) f ′ (θ) + A 3 (θ) f (θ)dθ (2)
then the system may be written concisely as

ż(t ) = A -1 ż(t -1) + Lz t , z t (θ) = z(t + θ).
This system may be represented, following the approach developed in [START_REF] Burns | Linear functional-differential equations as semigroups on product spaces[END_REF], by an operator model in Hilbert space given by the equation

ẋ = A x, x(t ) = v(t ) z t (•) , ( 3
)
where A is the infinitesimal generator of a C 0 -semigroup S(t ) = e A t given in the product space

M 2 = M 2 (-1, 0; R n ) def = R n × L 2 (-1, 0; R n )
and defined by

A x(t ) = A v(t ) z t (•) = Lz t (•) dz t (θ)/dθ , (4)
with the domain

D(A ) ⊂ M 2 given by v ϕ(•) : ϕ(•) ∈ H 1 , v = ϕ(0) -A -1 ϕ(-1) , (5)
where

H 1 = H 1 ([-1, 0]; R n ).
We consider the finite dimensional observation

y(t ) = C x(t ) ( 6
)
where C is a linear operator and y(t ) ∈ R p is a finite dimensional output. There are several ways to design the output operator C [START_REF] Salamon | Neutral functional differential equations and semigroups of operators[END_REF][START_REF] Salamon | Control and observation of neutral systems[END_REF][START_REF] Metel ′ Skiȋ | Criteria for the constructive identifiability and complete controllability of linear time-independent systems of neutral type[END_REF]. One of our goals in this paper is to investigate how to design a minimal output operator like

C x(t ) = C z(t ) or C x(t ) = C z(t -1), (7)
where C is a p × n matrix. More general outputs, for example with several and/or distributed delays are not considered in this paper. We want to use some results on exact controllability in order to analyze, by duality, the exact observability property in the infinite dimensional setting like, for example, in [START_REF] Tucsnak | Observation and control for operator semigroups[END_REF].

The operator C defined in ( 7) is linear but not bounded in M 2 . However, in both cases it is admissible in the following sense:

T 0 C S(t )x 0 2 R n dt ≤ κ 2 x 0 2 M 2 , ∀x 0 ∈ D(A ),
because it is bounded on D(A ). We recall that if x 0 ∈ D(A ) then S(t )x 0 ∈ D(A ), t ≥ 0 (see for example [START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF]). In fact, C is admissible in the resolvent norm:

x 0 -1 = R(λ, A )x 0 = (λI -A ) -1 x 0 M 2 , λ ∈ ρ(A ).
This is a consequence of the fact that C is a closed operator and takes value in a finite dimensional space (see [START_REF] Tucsnak | Observation and control for operator semigroups[END_REF]Def. 4.3.1] and comments on this Definition).

Definition 1.1. Let K be the output operator

K : M 2 -→ L 2 (0, T ; R p ), x 0 -→ K x 0 = C S(t )x 0 .
The system ( 1) is said to be approximately observable (or observable) if ker K = {0} and exactly observable (or continuously observable [START_REF] Salamon | Control and observation of neutral systems[END_REF]) if

K x 0 2 L 2 = T 0 C S(t )x 0 2 R p dt ≥ δ 2 x 0 2 M 2 , ( 8
)
for some constant δ.

3 This is the classic definition. In the case of a neutral type system with a finite dimensional output [START_REF] Khartovskiȋ | The complete controllability problem and its generalization for linear autonomous systems of neutral type[END_REF] the exact observability in this sense is not possible. It may be possible if we consider another topology for the initial states x 0 . Unlike approximate observability, which does not depend on the topology, exact observability depends essentially on the topology in the space. We can expect that, the given neutral type system is not exactly observable if we consider x 0 ∈ D(A ), with the norm of the graph and no longer in the topology of M 2 . Taking in account the result on exact controllability, it seems that (8) must be changed by taking a weaker norm for x 0 , namely the resolvent norm (λI -A ) -1 x 0 and considering the extension of the operator K to the completion of the space with this norm. In fact, we obtain the observability in the initial norm but we need some delay in the observation in the general case.

Exact observability can be investigated directly, but another way is to use the duality between exact observability and exact controllability. In [START_REF] Rabah | The analysis of exact controllability of neutral-type systems by the moment problem approach[END_REF] the conditions of exact controllability were given for the controlled system

ż(t ) = A -1 ż(t -1) + Lz t + Bu(t ).
In order to use the duality between observability and controllability, we need to compute the adjoint operator K * in the duality with respect to the pivot space M 2 in the embedding

X 1 ⊂ X = M 2 ⊂ X -1 , (9)
where X 1 = D(A ) with the graph norm noted x 1 and X -1 is the completion of the space M 2 with respect to the resolvent norm

x -1 = (λI -A ) -1 x M 2 . The duality relation is 〈K x 0 , u(•)〉 L 2 (0,T ;R p ) = x 0 , K * u(•) X 1 ,X d -1 , ( 10
)
where X d -1 is constructed as X -1 with A * instead of A (see [START_REF] Tucsnak | Observation and control for operator semigroups[END_REF] for example). Our purpose is to compute the adjoint operator

K * : L 2 (0, T ; R p ) → X d -1 .
The abstract formulation is well known. Exact controllability is dual with exact observability in the corresponding spaces with the corresponding topologies. It is expected that the operator K * corresponds to a control operator for some adjoint system.

We then need the expression of the adjoint state operator A * and the corresponding adjoint system in the same class: the class of neutral type systems. As it will be shown, the situation is not so simple. This is the object of Section 2. In Section 3 we return to the duality relation with the explicit expression of the adjoint system after formulations of exact controllability results. As the adjoint neutral type system has a slightly different structure, we give an explicit relation between the new neutral type system and the original one. After that we can give the expression of duality between exact controllability and exact observability. This enables to formulate the characterization of exact observability and to give the minimal time of observability. Some illustrative examples are given.

For the sake of completeness, we recall some results on approximate controllability (from [START_REF] Salamon | Control and observation of neutral systems[END_REF] and [START_REF] Bartosiewicz | Approximate controllability of neutral systems with delays in control[END_REF]) and formulate the duality with the corresponding notion of observability in our framework.

The adjoint system

In this section we give the expression of the adjoint system corresponding to the adjoint operator A * as the operator A corresponds to the system (1). Let us recall first the expression of the adjoint operator A * and its spectrum σ(A * ).

Proposition 2.1. ([16]) The adjoint operator A * is given by

A * w ψ(•) = (A * 2 (0)w + ψ(0) - d[ψ(θ)+A * 2 (θ)w ] dθ + A * 3 (θ)w , ( 11
)
with the domain D(A * ):

w, ψ(•) : ψ(θ) + A * 2 (θ)w ∈ H 1 , A * -1 A * 2 (0)w + ψ(0) = ψ(-1) + A * 2 (-1)w . ( 12
)
σ(A * ) consists of eigenvalues, roots of the equation det ∆ * (λ) = 0, where

∆ * (λ) = λI -λe -λ A * -1 - 0 -1 e λs A * 3 (s) + λA * 2 (s) ds. (13)
The adjoint operator A * in [START_REF] Bartosiewicz | Approximate controllability of neutral systems with delays in control[END_REF] seems to be different from a state operator generated by a neutral type system. However, we can construct a system of neutral type corresponding, in some sense, to the given adjoint operator.

Theorem 2.2. Let x be a solution of the abstract equation

ẋ = A * x, x(t ) = w(t ) ψ t (θ) . ( 14
)
Then the function w(t ) is the solution of the neutral type equation

ẇ (t + 1) = A * -1 ẇ(t ) + 0 -1 A * 2 (τ) ẇ(t + 1 + τ) + A * 3 (τ)w(t + 1 + τ) dτ. ( 15
)
Proof. Our purpose is to find the corresponding neutral type equation in R n . Equation ( 14) may be written as

∂ ∂t w(t ) ψ t (θ) = A * 2 (0)w(t ) + ψ t (0) - ∂[ψ t (θ)+A * 2 (θ)w (t )] ∂θ + A * 3 (θ)w(t ) Let us put r (θ) = A * 2 (θ)w + ψ(θ) and r (t , θ) = A * 2 (θ)w(t ) + ψ t (θ) = A * 2 (θ)w(t ) + ψ t (θ). ( 16
)
Then the operator A * may be rewritten as

A * w r (θ) -A * 2 (θ)w = r (0) -dr (θ) dθ + A * 3 (θ)w , ( 17
)
and the differential equation ẋ = A * x as a system of two equations:

∂ ∂t w(t ) r (t , θ) -A * 2 (θ)w(t ) = r (t , 0) -∂r (t ,θ) ∂θ + A * 3 (θ)w(t ) . ( 18
)
The second equation of this system may be written as a partial differential equation:

∂ ∂t r (t , θ) + ∂ ∂θ r (t , θ) = A * 2 (θ) ẇ (t ) + A * 3 (θ)w(t ). ( 19
)
The general solution of this equation is

r (t , θ) = f (t -θ) + θ 0 A * 2 (τ) ẇ(t -θ + τ) + A * 3 (τ)w(t -θ + τ) dτ, (20)
where f (tθ) is the solution of the homogeneous equation obtained from [START_REF] Wonham | Linear multivariable control: a geometric approach, 3rd Edition[END_REF]:

∂ ∂t r (t , θ) + ∂ ∂θ r (t , θ) = 0.
and the second term is a particular solution of [START_REF] Wonham | Linear multivariable control: a geometric approach, 3rd Edition[END_REF].

The first equation of the system [START_REF] Rudin | Functional analysis, 2nd Edition[END_REF] gives

ẇ(t ) = r (t , 0). (21)
From (20) (obtained from the second equation), putting θ = 0, we get with (21)

ẇ(t ) = r (t , 0) = f (t ). ( 22
)
Then (20) and (22) allow r (t , θ) to be written as follows:

r (t , θ) = ẇ(t -θ) + θ 0 A * 2 (τ) ẇ(t -θ + τ) + A * 3 (τ)w(t -θ + τ) dτ. ( 23
)
From the definition of the domain D(A * ) we obtain A * -1 r (0) = r (-1). For the function r (t , θ), this condition reads

r (t , -1) = A * -1 r (t , 0) = A * -1 ẇ(t ) (24)
and by (23) we have

r (t , -1) = ẇ(t + 1) - 0 -1 A * 2 (τ) ẇ(t + 1 + τ) + A * 3 (τ)w(t + 1 + τ) dτ. ( 25
)
Finally, from (24) and (25), we obtain the dual equation

ẇ (t + 1) = A * -1 ẇ(t ) + 0 -1 A * 2 (τ) ẇ(t + 1 + τ) + A * 3 (τ)w(t + 1 + τ) dτ. ( 26
)
On the other hand the solution of equation ( 18) is

e A * t x 0 = w(t ) ψ t (θ) = w(t ) r (t , θ) -A * 2 (θ)w(t ) , ( 27
)
where w(t ) is the solution of equation ( 26). If x 0 ∈ X then it is a mild solution.

This result may also be formulated, by simple duality (transposition), in the following way.

Theorem 2.3. Let x be a solution of the abstract equation

ẋ = Ã x, x(t ) = w(t ) ψ t (θ) ,
where the operator à is defined by

à w ψ(•) = A 2 (0)w + ψ(0) - d[ψ(θ)+A 2 (θ)w ] dθ + A 3 (θ)w ,
with the domain

D( Ã ) = w, ψ(•) : ψ(θ) + A 2 (θ)w ∈ H 1 , (A -1 A 2 (0) -A 2 (-1)) w = ψ(-1) -A -1 ψ(0) .
Then the function w(t ) is the solution of the neutral type equation

ẇ(t + 1) = A -1 ẇ(t ) + 0 -1 [A 2 (τ) ẇ(t + 1 + τ) + A 3 (τ)w(t + 1 + τ)] dτ. (28)
Let us now specify the relation between the solutions of neutral type equations ( 28) and [START_REF] Bensoussan | Representation and control of infinite-dimensional systems[END_REF]. Let us put w(t ) ψ t (θ) = e à t x0 = e à t w(0) ψ 0 (θ) , and

v(t ) z t (θ) = w(t + 1) -A -1 w(t ) w(t + 1 + θ) = e à t v(0) z 0 (θ) = e à t ξ 0 ,
where z 0 (θ) = w(t + 1) and v(0) = z 0 (0) -A -1 z 0 (-1). Our purpose is to give the explicit relation between the initial conditions x0 and ξ 0 :

x0 = w(0) ψ 0 (θ) , ξ 0 = v(0) z 0 (θ) .
The formal relation between these vectors is

x0 = w(0) ψ 0 (θ) = F ξ 0 = F w(1) -A -1 w(0) w(θ + 1) .
Theorem 2.4. The operator F representing the relation between initial conditions x0 and ξ 0 corresponding to the neutral type systems ( 1) and ( 28) is linear bounded and bounded invertible from X 1 to M 2 .

Proof. Let us calculate the explicit expression for the linear operator F . From ( 23) and ( 16), taking in account that we consider here the operator à instead of A * , we obtain r (0, θ) = ψ 0 (θ)

+ A 2 (θ)w(0) = ẇ(-θ) + θ 0 [A 2 (τ) ẇ(τ -θ) + A 3 (τ)w(τ -θ)] dτ, (29)
which can be written as

r (0, θ) = ẇ(-θ) + θ 0 [A 2 (θ -s) ẇ(-s) + A 3 (θ -s)w(-s)] ds.
Putting w(-s) = s 0 ẇ (-σ)dσ + w(0), we get r (0, θ) -

θ 0 A 3 (θ -s)ds • w(0) = ẇ(-θ) + θ 0 A 2 (θ -s) ẇ(-s) + A 3 (θ -s) s 0 ẇ(-σ)dσ ds. ( 30
)
This may be represented by the expression r (0, θ) -

θ 0 A 3 (θ -s)ds • w(0) = (I + V ) ẇ(-s), ( 31
)
where V is the Volterra operator defined by

V µ(•) = θ 0 A 2 (θ -s)µ(s) + A 3 (θ -s) s 0 µ(σ)dσ ds.
The operator V is a compact linear operator from

L 2 (-1, 0; R n ) to L 2 (-1, 0; R n ) with a spectrum σ(V ) = {0}. This implies that the operator I + V is bounded invertible on L 2 (-1, 0; R n ).
Let us now represent the operator F as a composition of operators according to the following commutative diagram v(0) z 0 (θ)

F P / / w(0) ẇ(-θ) Q w(0) ψ 0 (θ) R o o w(0) (I + V ) ẇ(-s)
where, as explained above, v(0) z 0 (θ) = w(1) -A -1 w(0) w(θ + 1) and

(I + V ) ẇ(-s) = r (0, θ) + -θ 0 A 3 (θ + s)ds • w(0),
where r (0, θ) is given as in (29). The operators P : X 1 -→ M 2 and R : M 2 -→ M 2 are bounded invertible. Moreover, as

I + V : L 2 (-1, 0; R n ) -→ L 2 (-1, 0; R n ) is bounded invertible, then Q : M 2 -→ M 2 is also bounded invertible.
We also need the following property of the bounded operator F -1 .

Proposition 2.5. For λ = σ( Ã ), the operator

F -a = F -1 (λI -Ã )
can be extended to a bounded (and bounded invertible) operator from M 2 to M 2 .

Proof. We need to prove that

F -1 (λI -Ã ) x0 ≤ C x0 , x0 ∈ D( Ã ),C > 0, ( 32
)
where • is the initial norm in M 2 . Let L 0 and D 0 be the subspaces

L 0 = {(0, ψ(•)) : ψ(•) ∈ L 2 (-1, 0; R n )}, D 0 = L 0 ∩ D( Ã ).
It is clear that D 0 is of finite co-dimension n, and this implies that it is enough to prove the relation (32) for x0 ∈ D 0 . Let x0 = (0, ψ 0 (•)) ∈ D 0 . The action of the operator F -a = F -1 (λI -Ã ) may be decomposed according to the following diagram

0 ψ 0 (•) (λI -Ã ) -→ -ψ 0 (0) λψ 0 (•) + ψ0 (•) R -1 -→ -ψ 0 (0) Ψ 0 (•) ↓ F -a ↓ Q -1 w(1) -A -1 w(0) w(θ + 1) P -1 o o ψ 0 (0) ẇ(-θ)
where

Ψ 0 (θ) = λψ 0 (•) + ψ0 (•) + -θ 0 A 3 (θ + s)ds • ψ 0 (0),
and the function w(•) is determinated from the equation (obtained from (29)):

λψ 0 (θ) + ψ0 (θ) -A 2 (θ)ψ 0 (0) = ẇ(-θ) + θ 0 [A 2 (τ) ẇ(τ -θ) + A 3 (τ)w(τ -θ)] dτ,
with the initial condition w(0) = -ψ 0 (0). Integrating the last equation from 0 to -1 -θ and taking in account the initial condition, we obtain

ψ 0 (-1 -θ) + λ -1-θ 0 ψ 0 (τ)dτ - -1-θ 0 A 2 (τ)dτ • ψ 0 (0) = -w(1 + θ) + -1-θ 0 τ 0 A 2 (s) ẇ(s -τ) + A 3 (s)w(s -τ)ds dτ.
Using a transformation in the double integration and the initial condition, we get

ψ 0 (-1 -θ) + λ -1-θ 0 ψ 0 (τ)dτ = -w(1 + θ) + -1-θ 0 A 2 (τ)w(1 + θ + τ) + τ 0 A 3 (s)w(s -τ)ds dτ,
and this can be written as

(I + V 1 )ψ 0 (-1 -θ) = (-I + V 2 )w(1 + θ), θ ∈ [-1, 0], where V 1 and V 2 are Volterra operators from L 2 (-1, 0; R n ) to L 2 (-1, 0; R n ).
Both operators have spectra concentrated at {0}. Then

w(1 + θ) = (-I + V 2 ) -1 (I + V 1 )ψ 0 (-1 -θ) = W ψ 0 (-1 -θ),
where W is a bounded invertible operator on L 2 (-1, 0; R n ). This enables the final expression for the operator F -a = F -1 (λI -Ã ) on the set D 0 to be obtained:

F -a 0 ψ 0 (θ) = W ψ 0 (-1 -θ) θ=0 -A -1 W ψ 0 (-1 -θ) θ=-1 W ψ 0 (-1 -θ) = w(1) -A -1 w(0) w(1 + θ) .
Taking in account A -1 w(0) = -A -1 ψ 0 (0) = -ψ 0 (-1), we can rewrite

W ψ 0 (-1 -θ) θ=0 -A -1 W ψ 0 (-1 -θ) θ=-1 = W ψ 0 (-1 -θ) θ=0 + ψ 0 (-1).
As a Volterra operator is quasinilpotent, we can write:

(-I + V 2 ) -1 = -∞ k=0 V k 2 and W = -∞ k=0 V k 2 (I + V 1 )
. This gives

W ψ 0 (-1 -θ)| θ=0 + ψ 0 (-1) = - ∞ k=1 V k 2 (I + V 1 ) + V 1 ψ 0 (-1 -θ) θ=0 = Gψ 0 (-1 -θ),
and G is a linear operator from a dense set of

L 2 (-1, 0; R n ) to R n . Since equalities V 1 ϕ(-1 -θ) θ=0 = - 0 -1 ϕ(τ) dτ, V 2 ϕ(-1 -θ) θ=0 = - 0 -1 A 2 (τ)ϕ(-1 -τ) τ 0 A 2 (s)ϕ(τ -s) ds dτ
define bounded operators from L 2 (-1, 0; R n ) to R n , we conclude that the operator G can be extended by continuity to a linear bounded operator from L 2 (-1, 0; R n ) to R n . As a consequence, the operator F -a is extended to a bounded (in the norm of M 2 ) operator, defined on L 0 by the formula

F -a 0 ψ(θ) = Gψ(-1 -θ) W ψ(-1 -θ) . ( 33
)
Let us observe that the subspace L 0 as well as its image F -a L 0 have codimension n in the space M 2 (a complement subspace for them is R n × {0} ⊂ M 2 ). Moreover, the mapping L 0

F -a
-→F -a L 0 is bijective. On the other hand, the subspace D 0 has codimension equals n in D( Ã ), F -a is defined on D( Ã ) by the formula

F -a = F -1 (λI -Ã ), ( 34
)
and the mapping D( Ã )

F -a -→D( Ã ) is also bijective. Besides, D( Ã ) + L 0 = D( Ã ) + F -a L 0 = M 2 .
Comparing all these facts, we conclude that the operator F -a , given by formulas (33), (34), can be extended to a bounded bijective operator on M 2 .

A direct consequence of this proposition is the following corollary.

Corollary 2.6. For all x ∈ D(A ), for λ = σ( Ã ), we have

c x ≤ (λI -Ã ) -1 F x ≤ C x ,
where • is the norm of the space M 2 .

The control system and duality

Consider the controlled neutral type system

ż(t ) = A -1 ż(t -1) + Lz t + Bu(t ), ( 35
)
where u(t ) ∈ L 2 (0, T ; R m ) is a m-dimensional control vector-function. This system may be represented by an operator model in Hilbert space given by the equation

ẋ = A x + Bu(t ), x(t ) = v(t ) z t (•) , ( 36
)
where Bu = (Bu, 0) is linear and bounded from R m to M 2 . We can note that B is not bounded from M 2 to X 1 because Bu = (Bu, 0) ∉ D(A ) if Bu = 0.

Exact controllability

Let us denote by R T ⊂ M 2 the reachable subspace of the system (36):

R T = R T u(•) = T 0 e A t Bu(t )dt : u(t ) ∈ L 2 (0, T ; R m ) ,
where R T : L 2 -→ M 2 is a linear bounded operator. As was pointed out in [START_REF] Ito | A linear quadratic optimal control for neutral systems[END_REF] and [START_REF] Rabah | The analysis of exact controllability of neutral-type systems by the moment problem approach[END_REF], R T ⊂ D(A ) for all T > 0. This implies that exact controllability may be defined as follows.

Definition 3.1. The system (36) is exactly controllable if R T = D(A ).

The abstract condition of exact controllability is (see [START_REF] Rudin | Functional analysis, 2nd Edition[END_REF] for example)

T 0 B * e A * t x 2 R m dt ≥ δ 2 x 2 X d -1 , ∀x 0 ∈ D(A * ), (37)
which means that the operator R T : L 2 -→ X 1 is onto. Here the space X d -1 is the completion of the space X = M 2 with respect to the norm

x X d -1 = (λI -A * ) -1 x M 2 , λ ∉ σ(A * ).
For the system (35) the condition of exact controllability is given by the following theorem (see [START_REF] Rabah | The analysis of exact controllability of neutral-type systems by the moment problem approach[END_REF]).

Theorem 3.2. The system (35) is exactly controllable at time T if and only if, for all

λ ∈ C, the following two conditions are verified i) rank λI -λe -λ A -1 - 0 -1 e λs [λA 2 (s) + A 3 (s)] ds B = n, ii) rank (λI -A -1 B) = n. The time of controllability is T > n 1 (A -1 , B).
The integer n 1 (A -1 , B) is the controllability index of the pair (A -1 , B) (see [START_REF] Wonham | Linear multivariable control: a geometric approach, 3rd Edition[END_REF]). If the delay is h, then the critical time is n 1 h.

Let us now consider the dual notion of observability for the adjoint system. The condition (37) is equivalent to the exact observability of the observed system

ẋ = A * x, y = B * x (38)
and the corresponding neutral type system is the system [START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF]. Then the conditions (i)-(ii) of Theorem 3.2 are necessary and sufficient for the exact controllability of the adjoint system (38). But what is the corresponding property for the associate neutral type system (26)? This question will be investigated in the following paragraph.

Duality

Consider the transposed controlled neutral type system

ż(t ) = A * -1 ż(t -1) + L * z t + C * u(t ), (39)
where

L * f = 0 -1 A * 2 (θ) f ′ (θ) + A * 3 (θ) f (θ)dθ.
Let A † be the generator of the semigroup e A † t generated by this equation (39). We cannot consider A * for this system because this operator does not correspond directly to this system as infinitesimal generator of the semigroup of solutions. The domain D(A † ) of the operator A † is given by v z(•)

: z ∈ H 1 ([-1, 0]; C n ), v = z(0) -A * -1 z(-1) . The spectrum of A † is σ(A † ) = {λ : ∆ * (λ) = 0} = σ(A * ). Let X † 1 be D(A † ) with the norm x X † 1 = (λI -A † )x M 2 , λ ∉ σ(A † ),
which is equivalent to the graph norm. Consider now the reachability operator for this system

R † T u(t ) = T 0 e A † t C * 0 u(t )dt .
From the properties of the operator R T , we can deduce that R † T is linear, bounded from L 2 (0, T ; R ) to X † 1 . The exact controllability for the system (39) can be formulated as

R † T = ImR † T = X † 1 .
The conditions of exact controllability for this system (39) may be obtained directly from Theorem 3.2.

Let us now consider the corresponding space X † -1 of linear functionals on X † 1 as the completion of the space X = M 2 with respect to the norm

x X † -1 = (λI -A † ) -1 x M 2 , λ ∉ σ(A † ).
We then have the embedding

X † 1 ⊂ X = M 2 ⊂ X † -1 . (40)
Then, for x ∈ X † 1 and y ∈ X † -1 , the functional acts as

x, y X † 1 ,X †d -1 = (λI -A † ) -1 x, (λI -A † * ) -1 y X , (41)
where A † * is the adjoint of the operator A † in M 2 , and the space X †d -1 is constructed as X †d -1 with A † * instead of A † (see [START_REF] Tucsnak | Observation and control for operator semigroups[END_REF] for example).

Let us note that the operator A † * is in fact the operator à defined in Section 2 (see Theorem 2.3 and later). We shall use the properties obtained for this operator.

Let us now consider the adjoint R † * T of R † T with respect to the duality induced by X † 1 and X †d -1 with the pivot space X = M 2 . Let x 0 ∈ X , then

R † T u(•), x 0 X 1 ,X d -1 = R † T u(•), x 0 X = T 0 e A † t C * 0 u(t )dt , x 0 X = T 0 C * 0 u(t ), e A † * t x 0 X dt .
Suppose now that x 0 = w(0)

ψ 0 (θ) ∈ D(A † * ).
Then, as a consequence of the results in Section 2, namely from (27) but for the operator A † * = Ã , we obtain

e A † * t x 0 = w(t ) ψ t (θ) = w(t ) r (t , θ -A 2 (θ)w(t ) , t ≥ 0. Hence, R † T u(•), x 0 X = T 0 〈u(t ),C w(t )〉 R p dt = u(•), R † * T x 0 L 2 .
On the other hand, we can write x 0 = F ξ 0 (cf. Proposition 2.5), where

ξ 0 = v(0) z 0 (θ) = z 0 (0) -A -1 z 0 (-1) z 0 (θ) = w(1) -A -1 w(0) w(θ + 1) ,
and then e A t ξ 0 = w(t + 1) -A -1 w(t ) w(t + 1 + θ) .

Suppose now that the relation ( 44) is verified for all ξ 0 ∈ D(A ), then from Corollary 2.6 we obtain

K ξ 0 L 2 ≥ δ (λI -A † * ) -1 F ξ 0 M 2 ≥ δc =δ 1 ξ 0 M 2 , ξ 0 ∈ D(A ).
This inequality can be extended by continuity to ξ 0 ∈ M 2 :

K ξ 0 L 2 ≥ δ 1 ξ 0 M 2 , ∀ξ 0 ∈ M 2
Conversely, suppose that the preceding relation is verified. For ξ 0 ∈ D(A ), and from Corollary 2.6, we get

ξ 0 ≥ 1 C λI -A † * -1 F ξ 0 M 2 ,
and then

K ξ 0 L 2 ≥ δ 1 C λI -A † * -1 F ξ 0 M 2 .
This is the relation (44) with δ = δ 1 /C . As the relations (43) and ( 44) are equivalent, the first assertion of the theorem is proved.

To prove item 2 of the theorem, it is sufficient to remark that the condition det A -1 = 0 is equivalent to the fact that the operator e A is bounded invertible (e A t is a group), and then the relations (43) and (44) are equivalent.

From this result and from Theorem 3.2 we can formulate the condition of exact observability.

Theorem 3.4. 1. The system [START_REF] Bensoussan | Representation and control of infinite-dimensional systems[END_REF] with the output y = C z(t -1) is exactly observable over [0, T ] if and only if (i) For all λ ∈ C, rank ∆ * (λ) C * = n, where ∆ * (λ) is defined in [START_REF] Salamon | Neutral functional differential equations and semigroups of operators[END_REF].

(ii) For all λ ∈ C, rank λI -A * -1 C * = n, (iii) T > n 1 (A * -1 ,C * ), where n 1 is the index of controllability for the pair (A * -1 ,C * ). 2. If det A -1 = 0, then Assertion 1 is verified for the output y(t ) = C z(t ).

Approximate controllability and observability

Let us now formulate the result on approximate controllability and observability. For more general duality relations between approximate controllability and observability of neutral type, we refer to the book [START_REF] Salamon | Control and observation of neutral systems[END_REF]. We give here a precise formulation in the light of our results on adjoint systems. Let us first recall the definition of approximate controllability.

Definition 4.1. The system (36) is approximately controllable at time T if

cl R T = M 2 ,
where cl R T is the closure of the attainable set R T at time T .

Sometimes approximate controllability is defined as cl

T >0 R T = M 2 ,
however for neutral type systems this notion of approximate controllability (as exact controllability) means that there is an universal time of controllability T 0 > 0 (see [START_REF] Salamon | Control and observation of neutral systems[END_REF]), i.e. such that:

cl

R T 0 = cl T >0 R T .
According to the relation (42) and the definition of observability we obtain the following result on duality between approximate controllability and observability.

If det

A -1 = 0, Assertion 1 of the theorem is verified for the output y(t ) = C x(t ) = C z(t ), and for the same time T .

Proof. The proof is a direct consequence of the definitions and (42).

The conditions of approximate observability may be obtained from the conditions of approximate controllability by duality. For our system, such conditions were obtained in [START_REF] Bartosiewicz | Approximate controllability of neutral systems with delays in control[END_REF] in the space W 2 1 ([-1, 0]; R n ). In our notations, this means that the reachability set R T 0 for the system (36) is dense in D(A ) with the norm of the graph. It is equivalent to the density of R T 0 in the space M 2 . In [START_REF] Salamon | Control and observation of neutral systems[END_REF], it is shown that approximate observability and approximate controllability for such neutral type systems are dual and this does not depend on the state space, and then duality holds in the space M 2 [2, Corollary 4.2.10].

From the necessary and sufficient conditions in [11, Th. 2] and Theorem 4.2, it is easy to see that approximate observability holds if the following two conditions are verified:

1. ∀λ ∈ C, rank ∆ * (λ) C * = n, 2. rank A * -1 C * = n. Note that the second condition is not necessary. The two conditions are verified when exact observability holds (the second condition is the condition (ii) of Theorem 3.4 for the particular case λ = 0). This emphasizes the difference between the concepts of exact and approximate observability.

Examples

Let us give some simple examples to illustrate our results. Example 1. Consider the system ż(t ) = ż(t -1), where z(t ) ∈ R n , n > 1, with two possible outputs y 0 (t ) = C 0 x(t ) = z(t ), y 1 (t ) = C 1 x(t ) = z(t -1).

The conditions of observability are verified, and the system is exactly observable for the output y 0 or y 1 .

Example 2. Consider the system ż1 (t ) = 0 ż2 (t ) = ż2 (t -1) + z 1 (t -1),

where z(t ) = (z 1 (t ), z 2 (t )) ∈ R 2 , with two possible output y 0 (t ) = C 0 x(t ) = z(t ), y 0 (t ) = C 1 x(t ) = z(t -1).

The system with the output y 1 is exactly observable for the time T > 1 and not observable for T = 1. The system with the output y 0 is not observable for any time T > 0.

Conclusion

For a large class of linear neutral type systems which include distributed delays we give the duality relation between exact controllability and exact observability. The characterization of exact observability is deduced.

grant N 514 238 438.

Let K be the output operator introduced in Definition 1.1. Then

, which implies for all x 0 ∈ X :

We can now formulate our main result on duality between exact controllability and exact observability. Theorem 3.3. 1. The system [START_REF] Bensoussan | Representation and control of infinite-dimensional systems[END_REF] with the output

if and only if the adjoint system (39) is exactly controllable at time T , i.e.

If det

A -1 = 0, Assertion 1 of the theorem is verified for the output

and for the same time T .

Proof. Let us recall that the exact controllability of the system (39) may be formulated by the equality ImR † T = X † 1 . Then, taking in account the embedding (40) and the duality product (41), we can write the condition of exact controllability as (see [START_REF] Rudin | Functional analysis, 2nd Edition[END_REF] for example):

Let ξ 0 = F -1 x 0 ∈ D(A ), where F : D(A ) → M 2 is the bounded invertible operator defined in Section 2. Let us remember that we have, from (42):

Then the inequality (43) is equivalent to


Philippe Angot

Jean-Paul Caltagirone

Pierre Fabrie

A spectacular solver of low-Mach multiphase Navier-Stokes problems under strong stresses

. This procedure, used as an approximate divergence-free velocity projection step, yields a velocity divergence vanishing as O(ε δt), δt being the time step. It only requires a few iterations of preconditioned conjugate gradients whatever the spatial mesh step h, if the penalty parameter ε is chosen sufficiently small up to machine precision, e.g. ε = 10 -14 . These methods prove to be efficient, fast and robust to accurately compute incompressible or low Mach multiphase flows under strong stresses: large mass density, viscosity or anisotropic permeability jumps, strong surface tension inducing large interface deformations, or with open boundary conditions, whereas other methods either cannot reach the suitable mesh convergence and run slower or simply crash.

Introduction: time-splitting methods for Navier-Stokes

We consider the numerical solution of unsteady incompressible viscous flows with variable density in a bounded domain Ω ⊂ R 3 , the frontier ∂ Ω := Γ = Γ D ∪ Γ N (Γ D ∩ Γ N = / 0) being subjected to a Dirichlet boundary condition v = v D on Γ D and a given traction stress vector σ (v, p)• n := -p n + µ(∇v + (∇v) t )• n = g on Γ N . Many splitting projection methods [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] of Chorin-Temam's type require, at each time step t n = n δt, the solution of the pressure Poisson equation to get the divergence-free velocity field v n+1 from a standard first-order velocity prediction step giving v n+1 : However, our experience with many Navier-Stokes solvers leads us to formulate the following conjecture. We cannot prove it rigorously, but several arguments explained further are likely to confirm it; see also [START_REF] Ph | Analysis for the fast vector penaltyprojection solver of incompressible multiphase Navier-Stokes/Brinkman problems[END_REF] for a more precise analysis.

     div δt ρ n+1 ∇φ n+1 = div v n+1 , ∇φ n+1 • n |Γ D = 0, φ n+1 |Γ N = 0. ⇒    ρ n+1 v n+1 -v n+1 δt + ∇φ n+1 = 0, φ n+1 = p n+1 -p n , div v n+1 = 0. (1)
Conjecture 1 (Numerical solution of Navier-Stokes equations). The main drawbacks of splitting projection methods arise from the introduction of the scalar pressure Poisson equation (1) involving a spatial derivative of mass density which inherently limits the consistency by degrading the original vector formulation of the Navier-Stokes problems.

This yields the first crucial point considered in the development of the new family of vector penalty-projection (VPP) methods.

Key idea 1 (Vector-penalty projection methods for Navier-Stokes equations.) To keep a fully vector formulation of the numerical methods for the solution of Navier-Stokes problems, it is essential to compute at each time step an accurate and curl-free approximation of the pressure gradient, and then reconstruct the pressure field if necessary. Thus, the primary unknows should be the velocity and pressure gradient vectors (v, ∇p). This is the objective of the VPP methods proposed in [START_REF] Ph | A spectacular vector penalty-projection method for Darcy and Navier-Stokes problems[END_REF][START_REF] Ph | A new fast method to compute saddlepoints in constrained optimization and applications[END_REF][START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF] to directly calculate the curl-free component

v n+1 := v n+1 -v n+1 of v n+1 .
It is not so surprising remembering that Euler 1 introduces the so-called Euler equations from the works of Bernouilli 2 and D'Alembert 3 by defining F p := -∇p as the pressure force which induces the fluid motion.

Fast Helmholtz-Hodge decompositions in bounded domains

For a given vector field v ∈ L 2 (Ω ) := L 2 (Ω ) 3 , later chosen as the predicted velocity (-v n+1 ), let us now consider the orthogonal Helmholtz-Hodge decomposition :

v = v φ + v ψ , with v φ = ∇φ , v ψ = rot ψ, div ψ = 0, v ψ • n |Γ = 0. (2)
The last normal boundary condition is necessary to ensure the L 2 -orthoganality and thus the uniqueness of the components v φ and v ψ , the scalar potential φ and vector potential ψ being determined up to an additive constant. Then, for v ∈ H div (Ω ) we have in Ω , supposed to be at least connected and possibly simply connected : Key idea 2 (Fast Helmholtz-Hodge decompositions of L 2 -vector fields.)

div v φ = div v, rot v φ = 0, v φ • n |Γ = v• n on Γ . and rot v ψ = rot v, div v ψ = 0, v ψ • n |Γ = 0 on Γ . (3)
The curl-free component v φ and the divergence-free (solenoidal) component v ψ of v satisfying (3) can be accurately and efficiently calculated with a penalization method by respectively solving the vector penalty-projection problem (V PP n ) and the rotational penalty-projection (RPP τ ) problem below, as proposed in [START_REF] Ph | Fast discrete Helmholtz-Hodge decompositions in bounded domains[END_REF]. Moreover, the computation can be extremely fast whatever the spatial mesh step h if the penalty parameter ε > 0 is chosen sufficiently small up to machine precision. Indeed, these problems take advantage of the splitting penalty method for saddlepoint [START_REF] Ph | A new fast method to compute saddlepoints in constrained optimization and applications[END_REF] in order to get adapted right-hand sides when ε → 0, and the effective conditioning becomes independent of ε and h, as confirmed numerically in [START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF][START_REF] Ph | Fast discrete Helmholtz-Hodge decompositions in bounded domains[END_REF][START_REF] Cheaytou | Vector penalty-projection methods for outflow boundary conditions with optimal second-order accuracy[END_REF].

(V PP n )    ε v ε φ -∇ div v ε φ = -∇ (div v) in Ω v ε φ • n |Γ = v• n on Γ ⇒      v ε φ = 1 ε ∇ div (v ε φ -v) , φ ε = 1 ε div (v ε φ -v), v ε φ = ∇φ ε (RPP τ ) ε v ε ψ + rot rot v ε ψ = rot (rot v) in Ω (rot v ε ψ ) ∧ n |Γ = (rot v) ∧ n on Γ ⇒      v ε ψ = 1 ε rot rot (v -v ε ψ ) , ψ ε = 1 ε rot (v -v ε ψ ), v ε ψ = rot ψ ε Theorem 3 (Optimal accuracy for (V PP n ) and (RPP τ ) problems).
For any v ∈ H div (Ω ) ∩ H rot (Ω ) and all ε > 0, the unique solution v ε φ or v ε ψ to the (V PP n ) or (RPP τ ) problems respectively, satisfies the optimal error estimates below: [START_REF] Ph | Analysis for the fast vector penaltyprojection solver of incompressible multiphase Navier-Stokes/Brinkman problems[END_REF]Sect. 3], Ω being connected or simply connected for (RPP τ ). The convergence towards 3-D Navier-Stokes weak solutions for constant density is also proved in [START_REF] Ph | Convergence results for the vector penalty-projection and twostep artificial compressibility methods, Discrete and Continuous Dynamical Systems[END_REF] for a continuous version of the (VPP ε ) method when ε → 0.

(i) v φ -v ε φ 1 + div (v -v ε φ ) 1 + φ -φ ε 2 ≤ c 1 (Ω ) v 0 ε, (ii) v ψ -v ε ψ 1 + rot (v -v ε ψ ) 1 + ψ -ψ ε 2 ≤ c 2 (Ω ) v 0 ε. Proof. See [4],
The condition number κ of the discrete left-hand side operators for these problems varies as κ = O(1/(ε h 2 )). Thus, the associated matrices are ill-conditioned in the usual sense for ε 1 since the number of iterations of a preconditioned conjugate gradient solver verifies N iter ≤ O( √ κ), which is the bound in the worst case for an arbitrary right-hand side (Arrhs). However, the amazingly fast convergence history observed in [3, Fig. 3] and [4, Fig. 2] for these P.D.E. problems with adapted right-hand sides (Adrhs) can be explained by defining the notion of effective condition number of the whole linear system; see [2, Corollary 1.2 & 1.3]. Indeed, the linear system itself is extremely well-conditioned in this latter case.

Let us look at a simple example where f ∈ H -1 (Ω ) or u ∈ H 1 0 (Ω ) are given:

(Arrhs) ε u ε -∆ u ε = f in Ω , u ε = 0 on Γ . or (Adrhs) ε u ε -∆ u ε = -∆ u in Ω , u ε = 0 on Γ . (4)
In this simple case, since the operator -∆ : H 1 0 (Ω ) -→ H -1 (Ω ) is a (bi-continuous) isomorphism between Hilbert spaces, we can even take here ε = 0 and it is clear that u 0 = u and u εu 1 ≤ c(Ω ) u 0 ε for (Adrhs). Now, let A h := -∆ h , a symmetric positive definite matrix, be the discrete Laplacian operator with homogeneous Dirichlet boundary condition, e.g. by standard finite differences. Then, the matrix system A ε,h :

= ε I + A h verifies κ := cond 2 (A ε,h ) -→ ε→0 cond 2 (A h ) = O(1/h 2 ).

Nevertheless using Neumann series, since A ε,h = A h (I + ε A -1

h ), we get the asymptotic expansion of the discrete solution u ε,h in both cases when ε < 1/ A -1 h :

(Arrhs) (ε I + A h ) u ε,h = f h , ⇒ u ε,h = A -1 h f h + O(ε). or (Adrhs) (ε I + A h ) u ε,h = A h u h , ⇒ u ε,h = u h + O(ε). (5)
Hence for ε 1, the zero-order term in these expansions is a good approximation of u ε,h , but its computation for (Arrhs) requires the solution of the linear system A h u 0,h = f h with a condition number as O(1/h 2 ), whereas u 0,h = u h for (Adrhs) and therefore, the effective condition number is bounded independently on both ε and h in this case.

In [2, Theorem 1.1 & Corollary 1.3], a generalization of such a result is proved for (V PP n ) problems, although the ∇(div ) operator is not invertible. A similar result also holds for (RPP τ ) problems using the identity -∆ u = rot (rot u) -∇(div u) for all vector field u in 3-D.

Vector penalty-projection methods for Navier-Stokes equations

Let us now consider the incompressible multiphase Navier-Stokes model during the time interval (0, T ) and supplemented with suitable initial and boundary conditions:

     ρ (∂ t v + (v• ∇)v) -2div (µ(ρ) d(v)) + ∇p = f in Ω × (0, T ), div v = 0 in Ω × (0, T ), ∂ t ρ + div (ρ v) = 0 in Ω × (0, T ). (6)
Here, d(v) := (∇v + (∇v) t ) /2 denotes the strain rate tensor and f = ρ g+σ st κ n Σ δ Σ includes the gravity force and surface tension on interfaces Σ , κ being the curvature of Σ and δ Σ the Dirac measure supported on Σ , which yields the immersed stress jump interface condition on Σ :

[[σ (v, p)• n]] Σ = σ st κ n Σ .
Then as proposed in [START_REF] Ph | A spectacular vector penalty-projection method for Darcy and Navier-Stokes problems[END_REF][START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF], using a phase function ϕ to track the motion and deformation of Σ by VOF-PLIC or level-set methods, e.g. [START_REF] Sarthou | Eulerian-Lagrangian grid coupling and penalty methods for the simulation of multiphase flows interacting with complex objects[END_REF] and the references therein, the fast (VPP ε ) method reads with a first-order linearly implicit scheme:

                               ρ(ϕ n ) v n+1 -v n δt + (v n • ∇) v n+1 -2div µ(ϕ n ) d( v n+1 ) + ∇p n = f n , ρ(ϕ n ) v n+1 -v n δt - 1 ε ∇ div v n+1 = 1 ε ∇ div v n+1 , with v n+1 • n |Γ = 0, v n+1 = v n+1 + v n+1 ,
velocity and pressure gradient corrections

φ n+1 = p n+1 -p n from: ∇φ n+1 = ∇ p n+1 -p n = -ρ(ϕ n ) v n+1 -v n δt ϕ n+1 -ϕ n δt + v n+1 • ∇ϕ n = 0. (7)
Key idea 4 (Vector penalty-projection methods for variable density flows.)

In the fast (VPP ε ) prediction-correction method [START_REF] Ph | Convergence results for the vector penalty-projection and twostep artificial compressibility methods, Discrete and Continuous Dynamical Systems[END_REF] for non-homogeneous Navier-Stokes equations, the mass density ρ is only included in the diagonal of the discrete velocity penalty-projection step and a cheap diagonal preconditioning is therefore recommended. This ensures not only the space-time consistency with the pressure gradient correction by keeping a fully vector formulation of the scheme, but also does not introduce any spatial derivative of ρ. This important feature makes the robustness of the method insensitive to large variations of density [START_REF] Ph | A spectacular vector penalty-projection method for Darcy and Navier-Stokes problems[END_REF][START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF]. Besides, the fully vector formulation of (VPP ε ) versions with a formally secondorder time scheme yields a natural expression of outflow boundary conditions with a given stress vector (traction or pseudo-traction). It ensures the optimal secondorder accuracy in time and space, both for the velocity and pressure fields [START_REF] Cheaytou | Vector penalty-projection methods for outflow boundary conditions with optimal second-order accuracy[END_REF].

For low-Mach number flows, typically when M < 0.2, the parameter ε has a physical meaning since it can be related as below [START_REF] Ph | Analysis for the fast vector penaltyprojection solver of incompressible multiphase Navier-Stokes/Brinkman problems[END_REF] to the Mach number M := V /c (V being a given reference velocity and c the speed of acoustic waves), the velocity divergence vanishing as O(χ T ) where χ T , χ S are the isothermal or isentropic compressibility coefficients of the fluid and γ := c p /c v ≥ 1:

ε δt = χ T = γ χ S = γ M 2 ρ V 2 , or also γ M 2 = ρ V 2 (ε δt) 1. ( 8
)
4 Sharp test cases with large density jumps or surface tension

The (VPP ε ) method, spatially discretized with the finite volume method on the cartesian staggered MAC grid or generalized MAC unstructured meshes, was tested and validated in [START_REF] Ph | A spectacular vector penalty-projection method for Darcy and Navier-Stokes problems[END_REF][START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF] against usual standard benchmark problems. In particular, the comparison with the scalar incremental projection (SIP) or Uzawa augmented Lagrangian (UAL) methods was excellent for the first benchmark in [START_REF] Hysing | Quantitative benchmark computations of two-dimensional bubble dynamics[END_REF]; see [START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF]Fig. 4].

As shown in the video during the talk, we have computed with the (VPP ε ) method combined with VOF-PLIC [START_REF] Sarthou | Eulerian-Lagrangian grid coupling and penalty methods for the simulation of multiphase flows interacting with complex objects[END_REF], the dispersed two-phase dynamics of air bubbles in a liquid melted steel at a temperature of about 850 • C. Here we have a large mass density ratio ρ l /ρ g = 8000, a viscosity ratio of µ l /µ g = 6 and a large surface tension constant σ st = 1.5 N/m which induces very large shape deformations of the bubble. A few results are shown in Fig. 1, either with an initial cylindrical bubble centered inside the steel column which is then rising along the median, or with a non centered initial bubble which rises by oscillating from the left to the right side.

Nevertheless for dispersed bubbly flows, it is difficult to compare our numerical method with others since most of them have difficulties to compute results with a suitable mesh convergence when the mass density ratio exceeds several hundreds; see [9, second benchmark]. Thus, to evaluate and validate the robustness of the (VPP ε ) method with respect to large density or viscosity ratios, we computed in [START_REF] Ph | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF]Fig. 5] the motion of a heavy solid ball which freely falls vertically in air with the gravity force f = ρ s g. Here we have : ρ s /ρ g = 10 6 and µ s /µ g = 10 17 .

Indeed, the SIP method crashes after a few time iterations. The UAL method is still able to compute the flow but with a larger velocity divergence and the computation is far more expensive than with (VPP ε ).

Fig. 1

 1 Fig. 1 2-D air bubble dynamics in melted steel (at about 850 • C) with (VPP ε ) method, ε = 10 -8 : cylindrical initial diameter d = 1 cm in a steel column L = 3 cm × H = 10 cm, ρ l /ρ g = 8000, µ l /µ g = 6, σ st = 1.5 N/m, g = 9.81 m/s -Isolines of the VOF-PLIC function at different times: LEFT: bubble initially centered in the column -RIGHT: bubble initially non-centered.