Philippe Destuynder
email: philippe.destuynder@cnam.fr

Caroline Fabre
email: caroline.fabre@u-psud.fr

Olivier Wilk
email: olivier.wilk@cnam.fr

Few numerical aspects of the Neumann-Kelvin's model with capillary in hydrodynamics

Keywords: AMS Subject Classification: Primary: 74F10, 35L85, 35R35; Secondary: 76B15, 35L70, 35Q35 water waves, instabilities, wave equations, spectral theory, non-linear models

A discussion and an improvement on the Neumann-Kelvin's model are suggested in this paper. This model is used in the simulation of progressive wave phenomenon. But as mentioned by several authors [CHEN 01], [YAN 00], [RIC 96], [LIH 78], this model is ill posed unless a capillary energy is introduced. The mathematical explanation is that a compactness inversion occurs if the capillary forces are omitted. Theoretical arguments and numerical simulations are used in the following which aim at giving an explanation of what happens from a mechanical point of view.

Introduction

Let us consider a bounded connected and non empty open set in R 3 denoted by Ω the boundary of which is ∂Ω. This set is a water pool with a bottom Γ 0 , a free surface Γ f and two sides Γ 1i and Γ 10 through which a water flow is respectively entering and leaving the pool (see figure 1). The unit normal along the boundary outwards Ω is ν. A submarine can be included in the pool and its boundary is Γ b . Another possibility is to have a pebble stone bottom. A steady flow is entering Ω through Γ 1i at the normal velocity U and getting out 1 through Γ 10 . There are two kinds of boundary conditions on Γ 1 depending if one considers the steady flow or the transient one. Concerning the steady flow the most realistic condition consists in fixing the normal velocity. But in the transient analysis, one can consider that the acoustic pressure in the wake is not perturbed. Nevertheless, the mathematical analysis is similar in each case. A comparaison in the numerical tests shows that the changes of boundary condition on the transient analysis are not meaningfull as far as the time delay of the simulation is small enough in order to avoid any return of surface waves. Another strategy is to use non-reflecting boundary conditions on Γ 1 . But this last possibility is a little bit complex because there are several kinds of waves on the surface and the water itself is assumed to be incompressible. Therefore, this point is not discussed in this paper. Because the fluid is assumed to be unviscid, it is assumed that the flow velocity can be modelled by a potential function Φ. From the incompressibility, one can ensure that -∆Φ = 0 in Ω. Let us denote by η the normal displacement of a geometrical point on the free surface of the water Γ f . We assume that the boundary conditions outside Γ 1 are

Entry

             ∂Φ ∂ν = 0 on Γ b ∪ Γ 0 , ∂Φ ∂ν = ∂η ∂t + ∇ s Φ.∇ s η, on Γ f , (1)
The boundary conditions on Γ 1 are defined separately in the following for the steady component and the transient one.

The stationary potential -say Φ 0 -is solution of the following set of equations:

                           -∆Φ 0 = 0, ∂Φ 0 ∂ν = U (e 1 , ν) on Γ 1 , ∂Φ 0 ∂ν = 0 on Γ b ∪ Γ 0 , ∂Φ 0 ∂ν = 0 on Γ f . (2)
Because Φ 0 is harmonic in Ω it is necessary that the Fredholm 's condition should be satisfied:

∂Ω ∂Φ 0 ∂ν = 0. (3)
In the present analysis, this relation is in fact a condition on the geometry of the boundary Γ 1 which traduces the incompressibility of the fluid. It can be explicited as follows and it is easy to ensure that it is satisfied:

Γ1 (e 1 , ν) = 0. (4)
In most classical analysis of progressive waves, the fonction Φ 0 is estimated by x 1 U which is the right solution for a parallelepipede. Let us notice that this is not correct for instance as soon as a submarine or a pebble stone bottom is considered inside Ω. These cases are those studied in the following. The solution method used in this paper for computing Φ 0 is a finite element one. But this step could also solved by a boundary integral method The Fredholm's condition (3) is not necessary for the transient state, as far as a Dirichlet boundary condition is prescribed on Γ 1 . One practical advantage is that it avoids to add an additional constraint on η. But this could be mathematically possible without new difficulties excepted in the notations. In order to split the steady state and the transient state let us set:

Φ = Φ 0 + ϕ, (5)
where Φ 0 is the solution of (2). Because Φ 0 is proportional to the boat velocity, one sets:

Φ 0 = U Φ0 , ( 6
)
where Φ0 is the solution of (2) with U = 1. The transient term ϕ should be solution of:

                         -∆ϕ = 0 in Ω, ϕ = 0 on Γ 1 , ∂ϕ ∂ν = 0 on Γ 0 ∪ Γ b , ∂ϕ ∂ν = ∂η ∂t + U ∇ s Φ0 .∇ s η + ∇ s ϕ.∇ s η on Γ f . (7)
The last boundary condition is a non-linear one because of the term ∇ s ϕ.∇ s η and it is omitted in the linear approximation. It will be discussed later on how to take into account non-linear terms when an instability occurs in the linear model and therefore the linearization wouldn't be justified. Several possibilities occur for the boundary condition satisfied by ϕ on Γ 1 . The one used in (7) traduces a vanishing acoustic pressure. Another possibility would be to prescribe that the total pressure (linearized expression) is zero (assuming also that the steady velocity is normal to Γ 1 ):

∂ϕ ∂t + U ∂ϕ ∂ν = 0. ( 8
)
Clearly there are many other possibilities in the definition of the boundary condition that should be satisfied by ϕ on this part of the boundary ∂Ω.

Furthermore it is certainly more realistic to use different boundary conditions at the flow-intake and the exit even if the mathematical analysis is similar. But their influence is restricted as far as the time delay for the computation doesn't allow the waves to reach the boundary of the open set Γ f .

The equilibrium of the free surface

Let us denote by p the pressure in the fluid. From Bernoulli theorem (see for instance G. Duvaut [DUV 90] or L. Landau-E. Lipschitz [LAL 71]), one obtains on the free surface of the water -say Γ f (x 3 = 0)-by considering only the linearized expression in η and ϕ of the pressure (p 0 is the pressure in the air over the free surface):

p = p 0 -̺ ∂ϕ ∂t - ̺ 2 |∇ϕ| 2 -̺U ∇ s Φ0 .∇ s ϕ -̺ U 2 2 |∇ s Φ0 | 2 -̺gη + σ∆ s η. (9)
In this expression σ is the capillary constant, g is the gravity and ̺ is the mass density of the fluid. The notation ∆ s is the Laplace operator restricted to the free surface.

Remark 1. One could consider a more complete expression for the capillary forces (the energy of which is proportional to the variation of surface) which is given by:

div s ( ∇ s η 1 + |∇ s η| 2 ).
But on the one hand the first additional non-linear term is a third order one, and on the other hand, this would lead us to use the space BV which allows discontinuities and which is much more complicated in the approximation method (see G. Aubert and P. Kornprobst [Aub 02]).

One can linearize the formula (9) by cancelling the term |∇ϕ| 2 or by linearizing it around a value ϕ 0 for instance corresponding to a stationary (time independent) flow different from the one induced by Φ 0 . A digging effect of the free surface due to the term -̺ U 2 2 |∇ s Φ0 | 2 can also appear in this stationary flow and it is denoted by η 0 . This so-called digging component is solution of:

     η 0 ∈ H 1 0 (Γ f ), -σ∆ s η 0 + ̺gη 0 + ̺U 2 ∇ s Φ0 .∇ s G(∇ s Φ0 .∇ s η 0 ) = - ̺U 2 2 |∇ s Φ0 | 2 , (10)
where G(f ) = ϕ| Γ f is the restriction on Γ f of the solution ϕ of the model ( 7)

where the boundary condition on Γ f is replaced by ∂ϕ ∂ν = f . This operator G plays an important role in the following. For example, if Φ0 = x 1 and if one restricts the preceding model to a one dimensional case, one obtains the following model for the digging effect:

   -σ∆ s η 0 + ̺gη 0 + ̺U 2 ∂G( ∂η 0 ∂x 1 ) ∂x 1 = - ̺U 2 2 . ( 11
)
In general, the expression of the digging effect depends on the steady flow which is described by Φ 0 . For instance it could be much more meaningful for shallow water and with different boundary condition on Γ 1 . For sake of simplicity in the writing, let us consider here that both η 0 and ϕ 0 can be neglected. Therefore the linearized pressure is given by:

p = p 0 -̺ ∂ϕ ∂t -̺U ∇ s Φ0 .∇ s ϕ -̺ U 2 2 |∇ s Φ0 | 2 -̺gη + σ∆ s η. (12)
Let us adopt for instance homogeneous Dirichlet boundary conditions for η on ∂Γ f . The equilibrium of the free surface is therefore explicited by (p 0 which acts on both side of Γ f is eliminated):

-σ∆ s η + ̺gη + ̺ ∂ϕ ∂t + ̺U ∇ s Φ0 .∇ s ϕ + ̺ U 2 2 |∇ s Φ0 | 2 = 0. ( 13
)
2 Formulation of the model in η

Let us introduce the functional space in which η is looked for:

V = H 1 0 (Γ f ). ( 14
)
By multiplying equation ( 13) by an arbitrary element of the space V , one obtains:

∀v ∈ V, σ Γ f ∇ s η.∇ s v + ̺g Γ f ηv + ̺ Γ f ∂ϕ ∂t v + ̺U Γ f ∇ s Φ0 .∇ s ϕv = - ̺U 2 2 Γ f |∇ s Φ0 | 2 v. (15)
Initial conditions should be added to this model. Let us set:

η(0) = η 0 and η(0) = η 1 on Γ f . (16)
One convenient possibility in the analysis of this model is to eliminate the function ϕ from equation (7). In fact, as far as the non-linear term ∇ s ϕ.∇ s η is omitted in (7), ϕ depends linearly on η and for eliminating ϕ, one can introduce the following operator denoted by G which is the so-called added mass operator. Let g be an arbitrary function defined on the boundary Γ f . The element G(g) is the restriction to Γ f of the solution to the following model:

                         -∆Φ = 0 in Ω, ∂Φ ∂ν = 0 on Γ 0 ∪ Γ b , Φ = 0 on Γ 1 , ∂Φ ∂ν = g on Γ f . (17)
It is proved in [DCF 09] that G -1 is an isomorphisme between the spaces H1/2 (Γ f ) and its dual. Definitions and details on these fractional spaces can be found in R. Adams

|||g||| 1/2 = inf v = g on Γ f v = 0 on Γ f ||v|| 1,Ω .
From the linearity of G one can write (let us recall that the non-linear term ∇ s ϕ.∇ s η is omitted):

ϕ |Γ f = G( ∂η ∂t ) + U G(∇ s η.∇ s Φ0 ). ( 18
)
By introducing this expression into the equation (13) which traduces the equilibrium of the free surface, one obtains formally (an integration by parts has been used and boundary terms have disappeared because v = 0 on ∂Γ f ):

                 ∀t ∈ [0, T ], η(t) ∈ H 1 0 (Γ f ), ∀v ∈ H 1 0 (Γ f ), m s ( ∂ 2 η ∂t 2 , v) + 2U c( ∂η ∂t , v) + a(η, v) = 0, η(0) = η 0 , ∂η ∂t (0) = η 1 on Γ f . (19)
with the notations:

                                           ∀η, v smooth enough : m s (η, v) = ̺ Γ f G(η)v, c(η, v) = ̺ 2 Γ f G(∇ s Φ0 .∇ s η)v -∇ s Φ0 .∇ s v G(η) -∆ s Φ0 G(η)v , a(η, v) = Γ f σ∇ s η.∇ s v + ̺gηv -̺U 2 Γ f G(∇ s Φ0 .∇ s η)∇ s Φ0 .∇ s v + ∆ s Φ0 G(∇ s Φ0 .∇ s η)v . (20)
In many examples the effect of the term ∆ s Φ0 is numerically neglectible as far as they are no surface ship. This point is observed in the numerical tests, even if it is not zero. First of all let us check few useful properties of the bilinear forms m s , c and a. Let us notice that they are well defined if η and v are functions lying in the space H 1 0 (Γ f ). It is also worth noting that the classical symmetry property holds for G:

∀η, v ∈ L 2 (Γ f ), Ω ∇Φ(η).∇Φ(v) = Γ f G(η)Φ(v) = Γ f G(v)Φ(η), (21)
which implies the symmetry of the bilinear form m s . This is also true for the three first terms in the definition of a. One can also claim that the bilinear form c restricted to its two first terms is anti-symmetrical. This is due to fact that it represents a gyroscopic coupling. In fact, ∇ s η is a rotation of the unit normal ν to the free surface Γ f . But the main properties are the coerciveness properties which are summarized here-after. In order to make sense to the expressions which appear in the definitions of the bilinear forms one assumes that Φ0 ∈ C 3 (Γ f ), (see P. Grisvard [GRI 86] or M. Borsuk and V. Kondratiev [KON 06] for the justification) because one has in this case:

                           | Γ f ∆ s Φ0 G(η)v| ≤ c|||η||| -1/2,Γ f |||v||| -1/2,Γ f ,
and:

| Γ f -∆ s Φ0 G(∇ s Φ0 .∇ s η)v| ≤ c|||v||| 1/2,Γ f ||∇ s Φ0 .∇ s η|| -1/2,Γ f ≤ c ′ ||η|| 1/2,Γ f |||v||| 1/2,Γ f (22)
Let us mention a result [DCF 09] which enables to justify the model used.

Theorem 1. Let us assume that min(̺g, σ) > 0. Then, there exists constants c i > 0 and a scalar λ 1 > 0 such that:

∀η, v ∈ H 1 0 (Γ f ) : a(v, v) ≥ c 0 ||v|| 2 1,Γ f -c 5 U 2 ||v|| 2 1/2,Γ f ≥ c 6 (λ 1 -U 2 )||v|| 2 1,Γ f , m s (v, v) ≥ c 1 |||v||| 2 -1/2,Γ f , |m s (η, v)| ≤ c 2 |||η||| -1/2,Γ f |||v||| -1/2,Γ f , |a(η, v)| ≤ c 3 ||η|| 1,Γ f ||v|| 1,Γ f , |c(η, v)| ≤ c 4 ||η|| 1/2,Γ f ||v|| 1/2,Γ f .
The main application concerns the existence and uniqueness of a solution to the coupled model (19) for any velocity U and for ad hoc initial conditions η 0 and η 1 .

Remark 2. A natural question is to obtain the best constant λ 1 which appears in the previous theorem. This is discussed in the following. It leads to the first critical velocity which would be equal to √ λ 1 .

Remark 3. One can consider several kind of boundary conditions on ∂Γ f concerning η. For instance, one can set:

∂η ∂ν + αη = 0 sur ∂Γ f , (23)
which is a Robin's boundary condition. If α = 0 there is a perfect gliding of the water, if α = ∞ there is an adhesion of the water which is the case considered here. But the influence on the results is not important as far as the time is small enough in order to avoid rebounds of waves on the boundary ∂Γ f . One can also consider a condition as:

∂η ∂t + γ ∂η ∂ν + δη = 0, ( 24
)
and so on. The most suitable boundary condition for η is certainly a transparency one in order to avoid reflection. But it would also be useful to use such a condition for ϕ on Γ 1 .

About the static instabilities

The existence (and uniqueness) theorem of a solution in η is true for any U . But the stability requires the hypothesis: U < √ λ 1 where λ 1 is a constant which appears in the coerciveness on the space H 1 0 (Γ f ) of a(., .) (see theorem 1). Let us discuss the best value for the constant λ 1 in the particular case where the term ∆ s Φ0 is zero (neglected). In this case, the bilinear form a is symmetrical and c is antisymmetrical. Two cases have to be considered depending if there is a capillary term or not (σ > 0 and σ = 0).

Let us introduce the eigenvalue model:

                           find (λ, w) ∈ R + × H 1 0 (Γ f ) such that: ∀v ∈ H 1 0 (Γ f ) : λ̺ Γ f G(∇ s Φ0 .∇ s w)∇ s Φ0 .∇ s v = σ Γ f ∇ s w.∇ s v + ̺g Γ f wv, ̺ Γ f G(∇ s Φ0 .∇ s w)∇ s Φ0 .∇ s w = ̺ Ω |∇Φ(∇ s Φ0 .∇ s w)| 2 = 1. ( 25
)
The non-trivial solutions of:

v ∈ H 1 0 (Γ f ), ∇ s Φ0 .∇ s η = 0,
play a particular role in the stability analysis. For instance in one dimension, one can prove (see [DCF 09] that the only solution is zero. But this result is true in many cases excepted may be for very particular geometries of the surface Γ f . From general spectral theory of linear operators one can deduce the following theorem:

Theorem 2. Let us assume for sake of simplicity that the set:

K 0 = {v ∈ H 1 0 (Γ f ), ∇ s Φ0 .∇ s v = 0},
is reduced to {0}. Two cases are considered in which the results are inverted. 1 • If σ > 0, there exists a countable set of elements denoted by

(λ n , w n ) in R + * × H 1 0 (Γ f ) solution of (

25). Each term of the sequence (ordered by increasing values) λ n has a finite multiplicity (∞ is the only accumulation point). The family {

w n √ λ n } is an Hilbert basis of the space H 1 0 (Γ f ). The smallest eigenvalue denoted by λ 1 will be the best constant in theorem 1.

2 • If σ = 0.
Let us restrict our analysis to the one dimensional case. The result is opposite to the one of the first case. There exists a countable set of solutions (λ n , w n ) ∈ R + * ×L 2 (Γ f ) and the terms of the sequence λ n are ordered by decreasing values (the largest value is denoted by λ 1 ). The multiplicity of each term is finite and 0 is the only accumulation point. The family {w n } is an Hilbert basis of the space L 2 (Γ f ).

Proof. Let us begin with σ > 0. The result is a direct consequence of the spectral theory for linear operator as it is presented P.A. Raviart and J.M. Thomas [RATO 83]. In fact the bilinear form:

(η, v) → σ Γ f ∇ s η.∇ s v + ̺g Γ f ηv,
is symmetrical, continuous and coercive on the space H 1 0 (Γ f ). Furthermore the bilinear and symmetrical form:

(η, v) → Γ f G(∇ s Φ 0 .∇ s η)∇ s Φ 0 .∇ s η, is continuous on the space H 1/2 (Γ f ).
From the properties of G -1 which is an isomorphism [DCF 09] between the spaces H1/2 (Γ f ) and its dual, one obtains:

       ∃c > 0, such that ∀η ∈ H1/2 (Γ f ), Γ f G(∇ s Φ 0 .∇ s η)∇ s Φ 0 .∇ s η ≥ c|||∇ s Φ 0 .∇ s η||| 2 -1/2,Γ f .
The continuous mapping:

j : v ∈ H 1 0 (Γ f ), → j(v) = |||∇ s Φ 0 .∇ s v||| -1/2,Γ f , (26)
is a norm on the space H 1 0 (Γ f ) (hypothesis). But even if it is a norm, it is not always equivalent to the one of the space H1/2 (Γ f ). The completed space of H 1 0 (Γ f ) equiped with the norm (26), is at least included into [ H1/2 (Γ f )] ′ and contains H1/2 (Γ f ). The natural embedding from H 1 0 (Γ f ) into H1/2 (Γ f ) being compact, it is also compact from H 1 0 (Γ f ) into this completed space. The assumptions of the spectral theory formulated in [RATO 83] are satisfied and the result is proved if σ > 0.

Let us consider the second case: σ = 0 which is much more tricky. Now the bilinear and symmetrical form:

(η, v) → ̺g Γj ηv,
is only equivalent to the square of the norm in the space L 2 (Γ f ). It has been assumed that j defined at (26) is a norm on the space H 1 0 (Γ f ). The embedding from H1/2 (Γ f ) into L 2 (Γ f ) being compact, the variational eigenvalue problem is inverted. This completes the proof of theorem 2 by applying the spectral theory adapted to the variational formulation [RATO 83]. ✷ Remark 4. In the first case where σ > 0 one can characterize the constant λ 1 which appears in theorem 1 as the smallest eigenvalue of the model (25). Therefore there are no confusion in the notation. Let us underline that λ 1 which is the smallest eigenvalue solution of (25) satisfies:

∀v ∈ H 1 0 (Γ f ), λ 1 ̺ Γ f G(∇ s Φ 0 .∇ s v)∇ s Φ 0 .∇ s v ≤ σ Γ f |∇ s v| 2 + ̺g Γ f v 2 . ( 27
)
Hence:

       ∀δ ∈ [0, 1], ∀v ∈ H 1 0 (Γ f ), a s (v, v) ≥ inf(σ, ̺g)(1 -δ)||v|| 2 1,Γ f +̺(δλ 1 -U 2 ) Γ f G(∇ s Φ 0 .∇ s v)∇ s Φ 0 .∇ s v. ( 28
)
Thus λ 1 is the best constant in theorem 1 and it enables to characterize the critical velocity U c corresponding to the instability of the linear wave model in which the term ∆ s Φ0 is neglected, by:

U c = λ 1 . (29)
If U ≥ U c the coerciveness of the stiffness bilinear form on the spaceH 1 0 (Γ f ) is lost. But it is still true for the complementary of the finite dimensional space spanned by the eigenvectors corresponding to the eigenvalues λ n such that U ≥ √ λ n .

Let us now assume that σ = 0. The coerciveness of a s is no more true on the space

H 1 0 (Γ f ) but only on L 2 (Γ f ) which is not contained in H 1/2 00 (Γ f ).
In fact, the contrary can be true and this would be a so-called mathematical inversion of the compactness. For any velocity U > 0, there exists an infinite number of instable eigenmodes. They are more and more local (the wave length is smaller and smaller). Therefore, in a numerical approximation, the more the mesh is refined, the larger is the number of instabilities. This remark which has already been formulated in a different numerical framework by Xiao-bo Chen [CHEN 02], condemns the Neuman-Kelvin's model without capillary. In fact, the variational model ( 19) is fully instable for U > 0 excepted for a finite dimensional space spanned by a finite number of eigenmodes (see theorem 2) the eigenvalues of which -say λ n solution of (25) with σ = 0-would satisfy:

λ n ≥ U 2 . ( 30
)
One can draw a strange conclusion: if U is small enough and if the mesh size in a numerical approximation is large enough, the Neumann-Kelvin is stable. This is really disturbing because the true model is not. One has a numerical filtering of the instabilities due to mesh size which would be too large.

The conclusion of this section is that the Neumann-Kelvin model is non physical (instable) as far as the capillary is not taken into account (excepted if U = 0). This result is obviously in contradiction with the usual argument that the capillary is very small and can be neglected compared to the gravity effect.

Remark 5. The dimension of the instable space corresponds exactly to the number of eigenvalues λ n solution of (25) and such that:

λ n ≤ U 2 .

Numerical discussion in one dimension

Our goal is mainly to give a numerical illustration of the theoretical results obtained for the progressive wave model. The accuracy of the tests must be sufficient in order to be able to claim that our conclusions are valid. This justifies that all the tests are restricted to one dimensional models. There are four steps. The first one is to compute the steady flow in order to check the influence of the depth of a submarine with respect to the free surface of the water or the influence of a random distribution of pebble stones on the bottom of the pool. In a second step the eigenvalue problem which characterizes the static instabilities is discussed numerically. The influence of the mesh size, of the position of the submarine and of the capillary coefficient are examined. In a third step few theoretical results concerning the evolution of the spectral density are compared with the numerical results. Finally the fourth step concerns a simulation of the progressive waves on the free surface of the water.

Computation of the steady flow Φ 0

The solution method used is a finite element method. Conformal first degree polynomials have been used with a triangular mesh for Ω. The results are quite standard but the goal is to show that the velocity on the free surface is different from U because of the effect of the submarine or the pebble stone bottom. Something similar occurs for a shallow water model with irregular bottom. The coarser meshes used are represented on figure 2. The point is important because the fact that Φ 0 is different from x 1 + c destroys the possibility to use an integral formulation as it is done in most of the published papers (see for instance

Convergence of the critical eigenvalues

The stability of the progressive wave model depends on the positiveness of the eigenvalues of the model defined at (25). Hence a numerical study is carried out in order to check the accuracy of the method which enables to compute the critical velocities U c defined as the square root of these eigenvalues. The lowest defines the first critical velocity of the steady flow. One important feature of this mechanical phenomenon is that the spectral density is quite large and thus 1/2 00 (Γ f ). Secondly, the operator is not local and couples all the point of the free surface Γ f through the operator G. For instance, if a lack of convergence appears near a corner inside of Ω it slows down the convergence everywhere and thus the eigenvalue model ( 25) is also degraded. The smallest eigenvalue is the square of the critical flow velocity at which a dynamic instability can occur.

Influence of the depth of the submarine

The computation shown in the following are performed with a refined mesh in the previous analysis. The evolution of the ten smaller eigenvalues (square of the critical velocities) have been plotted for several depths of the submarine from the free surface on figure 7. The laplacian ∆ s Φ0 is kept. The smallest eigenvalue is smaller than in the analytical case used for the analysis of the convergence in the previous subsection when the submarine is quite close to the free surface and the inverse if it is quite deep. This is due to the Venturi effect between the submarine and the free surface for small depth and between the submarine and the bottom for large depth. The results obtained for σ = 0 are also drawn on figure 7. One can see that the eigenvalues are much smaller. In a theoretical analysis of the eigen value model (25), it has been proved that for σ = 0 the sequence of the eigenvalues is upper bounded and the lower bound is zero (for the continuous model and are close to zero for a given finite element mesh). The lower graph on figure 7 shows that if the submarine is deep enough, the eigenvalues are very close from one to the other. And according to the theoretical analysis, this cluster of values tends to zero when the mesh size tends also to zero (for any depth of the submarine). At the opposite for σ > 0 the sequence of eigenvalues tends to the infinity and is lower bounded by a strictly positive number. Thus one can see that even if the capillary force are small their global energy is not. Hence it is not possible to cancel them in a realistic model for progressive waves.

Influence of the depth of the pool with pebble stone bottom

The critical velocities depend on the shape of the pebble stone bottom but also on the depth of the pool. Similarly to the case of the submarine, the evolution of the spectrum has been computed and compared. The results are plotted on figure 8 depending if the capillary phenomenon is taken into account or not. Furthermore, the computation with and without the laplacian term have been compared but no meaningful influence of this term has been observed.

About the numerical scheme used for progressive wave

Two time step schemes are discussed in this subsection but one has to be careful about the stability property. The test examples correspond to a submarine which is quite close to the free surface of the water. It is represented on the figures of section 3.4 (see figures 9 and 10).

A first scheme for a flow velocity such that U < U c .

Let us use the progressive wave model as follows where V h is the approximation space for the displacement η (the laplacian ∆ s Φ0 is neglected for the theoretical discussion):

             find η n+1 ∈ V h such that ∀v ∈ V h , m s (η n+1 -2η n + η n-1 , v) + 2U ∆tc s (η n+1 -η n , v) + ∆t 2 2 a s (η n+1 + η n , v) = F (v). ( 31
)
Setting v = η n+1 -η n in this relation, one obtains (the gyroscopic term disappears because C s is antisymmetrical):

m s ( η n+1 -η n ∆t , η n+1 -η n ∆t ) + 1 2 a s (η n+1 , η n+1 ) = m s ( η n -η n-1 ∆t , η n+1 -η n ∆t ) + 1 2 a s (η n , η n ). (32)
Finally, from the Schwarz's inequality (which is optimal) one deduces that:

1 2 m s ( η n+1 -η n ∆t , η n+1 -η n ∆t ) + 1 2 a s (η n+1 , η n+1 ) ≤ 1 2 m s ( η n -η n-1 ∆t , η n -η n-1 ∆t ) + 1 2 a s (η n , η n ). (33)
This inequality ensures the stability of the scheme as soon as U < U c because the terms above represent the energy of the system. But if U ≥ U c the bilinear form a s is no more positive. Let use an Hilbert basis of the space H 1 0 (Γ f ) solution of the eigenvalue problem:

           find w, λ ∈ H 1 0 (Γ f ) × R such that: ∀v ∈ H 1 0 (Γ f ), λm s (w, v) = a s (w, v), m s (w, w) = 1. (34)
The family of solutions {w n } is an Hilbert basis of the space [ H1/2 (Γ f )] ′ which is ordered such that:

λ 1 ≤ λ 2 ≤ . . . ≤ λ n ≤ λ n+1 ≤ . . .
Furthermore the basis is also conjugate versus the bilinear form a s . Let us denote by P the smallest integer such that λ P > 0. Let us set:

η n = n≥1 α n j w j . (35)
Hence from (32) one deduces that:

∀j ≥ 1, ∀n ≥ 1, | α n+1 j -α n j ∆t | 2 + λ j 2 |α n+1 j | 2 = ( α n+1 j -α n j ∆t ) ( α n j -α n-1 j ∆t ) + 1 2 λ j |α n j | 2 . ( 36
)
This implies the stability for any component α n j such that j ≥ P . But for j < P one only has from an induction argument (where the constant c depends on initial conditions):

| α n+1 j -α n j ∆t | 2 -|λ j ||α n+1 j | 2 ≤ c. (37)
Because this estimate is optimal (Schwarz's inequality) it proves that the solution η n j is instable. One can check it by considering the particular sequences which satisfies (37) and which tends to the infinity as soon as ∆t |λ j | < 2:

α n+1 j = α n j 1 -∆t |λ j | . ( 38
)
But the numerical test performed on the progressive wave model are still more convincing. The simulation have performed with and without the surface laplacian but the result are quite identical (the results presented are obtained with this term). Furthermore the influence of the capillary is analyzed and one can see that the influence is fundamental, even for flow velocities U smaller than the first critical value U c . The results have been plotted on figure 9 for σ = 7, 5 10 -2 N m -2 (with capillary) and σ = 0 (without capillary), ̺ = 10 3 Kg. One can observe that as soon as the critical velocity (here U c = 0.161 m/s) is overtaken, the scheme is fully instable and much more than the physical solution. This transfer of energy is due to the gyroscopic coupling (because it is the only mean to do that) which extracts energy from the stable eigenmodes and delivers it onto the instable eigenmodes. For this reason one can claim that the scheme (33) should be avoided for the development of the instability. This result which is true for the linear model is also true for the non linear progressive wave model.

A more appropriate scheme of the instable wave model

U ≥ U c .
Let us now consider another time step scheme for the wave equation by setting (let us again consider that the laplacian δ s Φ0 is neglected for the theoretical discussion) :

             find η n+1 ∈ V h such that ∀v ∈ V h , m s (η n+1 -2η n + η n-1 , v) + U ∆tc s (η n+1 -η n-1 , v) + ∆t 2 2 a s (η n+1 + η n-1 , v) = F (v). ( 39
)
Figure 9: Instable time-evolution of the energy for the linear model for the scheme (31) and for several velocities (U > U c and U < U c ). Top for σ = 0 and bottom for σ = 7, 5 10 -2 .

Setting v = η n+1 -η n-1 in this relation, one obtains (the gyroscopic term disappears because C s is antisymmetrical):

m s ( η n+1 -η n ∆t , η n+1 -η n ∆t ) + 1 2 a s (η n+1 , η n+1 ) = m s ( η n -η n-1 ∆t , η n -η n-1 ∆t ) + 1 2 a s (η n-1 , η n-1 ). (40)
In fact this quantity is an approximation of the energy E(n) of the mechanical system (progressive waves) at time t = n∆t is obtained as far as η n and η n±1 are close, because:

E(n) = 1 2 m s ( η n+1 -η n ∆t , η n+1 -η n ∆t ) + 1 2 a s (η n , η n ). ( 41
)
The important point is that this expression is constant for the scheme (39). Hence it keeps the energy or at least a quantity which is equivalent to it. Therefore even for instable movement this scheme is better adapted. The numerical tests are plotted on figure 10. The values of the physical constants are the same as in the scheme (31). One can observe that the stabilities observed are close to the true physical solution as far as a capillary energy is kept. The justification is that only the instable eigenmodes generate instable approximate solutions.

Visualization of the surface wake

Let us consider in this section two simulations of the progressive water waves occuring at the free surface. The first case corresponds to a submarine moving forwards from the right to the left at a given velocity and quite close to the free surface (see figure 11). The second one concerns the influence of a rocky bottom (see figure 12). The geometry of the bottom is the one shown on figure 2. It is a superposition of several sinusoïdes. It is an artificial approximation of real pebble stone bottom. But it should be point out that a slight random perturbation of the pebbles stones on the bottom can induce a large variation on the shape of the instabilities as far as the wave lenghts are close to the one of the instabilities represented by the eigenvectors solution of (25). But this would be another study. One important point is that the instabilities on the surface of the water are localized (at the beginning) above the submarine and everywhere in the second case concerning the pebble stone bottom. Furthermore the pebble stone distribution has been chosen such that the wave lenght is different from the one of the eigenmodes of the pool in order to avoid artificial resonances. In both cases the movement is initiated by the digging effect η 0 due to the term -̺ 2 |∇ s Φ 0 | 2 at the right hand side of the wave equation (10).

The case of the submarine

One can observe that the surface becomes instable at the vertical of the submarine. Thus the instability is transfered everywhere on the surface. This is due to the fact that a Venturi effect appears between the surface of the water and the submarine. Therefore the velocity is locally larger and the surface is sucked downwards in a first step. Then a wave appears in the wake. This phenomenon is always true but it is much more amplified over the critical velocity. Because the spectral density is high around the first eigenvalue which is the square of the first critical velocity, the mechanism becomes more and more complex for larger velocities which imply several eigenvectors.

Visualization of the surface wave for the pebble stone bottom

Several velocities have been considered for the steady flow. The depth is fixed and the four first graphs correspond to a velocity smaller than the critical value obtained from the smallest eigenvalue of (25). The simulation is performed with the linear model and the shape of the bottom is represented on the last graph of these graphics. The two following correspond to instable velocities. One can observe on figure 12 that the progressive waves depend strongly on the value of this velocity compared to the first critical value. But because of the high spectral density several instable eigenvectors are concerned for instance if U ≥ 0.2m/s. This is justified from figure 6. Hence the instable movement which is generalized on all the free surface is a very complex one.

The wake with a non-linear term

The two cases (submarine and pebble stone bottom) are discussed in this subsection. For U ≥ √ λ 1 the linearized model of the surface is instable. Simulations where the non-linear term:

̺ 2 Γ f [|∇ s ϕ| 2 + | ∂η ∂t + U ∇ s Φ0 .∇ s η| 2 ],
has been added to the equation (19). In the numerical computation this term is computed using backward time step difference in order to be explicit on the non linear term. One could object that there are other non-linear terms which have been omitted. Therefore the simulation given hereafter are only an indication of what can happen.

The submarine (non-linear) The visualization of the wake observed is plotted on figure 13 for the submarine. One can observe an amplification of the solution compared to the linear one, even when the the linearized model is stable. The pebble stone bottom (non-linear) One point in the non linear case is that the instability seems to be more structured that in the linear case. There can be a reason. The non linear term that has been introduced (which is not complete) selects among the instabilities those which are non-linearly stable and those which are instable.

Concerning the non-linear analysis, the results published up to now mainly concern a simplified model in which there are no boundary condition and furthermore there is no submarine or pebble stone bottom. The first mathematical contribution is the paper by T. Brooke Benjamin and J. E. Feir [BEN 67] who detected the possibility of an instability on a water surface flow under the effect of a periodic paddle excitation at a given frequency. Most of the references (both theoretical and applied ones) [IOS 96], [DIKH 99], [CHA 07] concern a finite dimensional approximation and mainly without capillary. Furthermore the numerical tests are mostly restricted to case where one can use the integral representation for the velocity potential of the solution, which is not adapted to the geometry discussed here.

Conclusion

In this paper a modification of the Neumann-Kelvin's wave model for the instable progressive wave is discussed from the numerical point of view. Let us underline the main points.

1 The role of the capillary is explained in mathematical framework for the physical equations for progressive waves. It is proved that the classical Neumann-Kelvin's model is fully instable and should be modified in order to be mechanically reliable.

2 It is shown that the best way to take the capillary into account is to use a formulation where the normal displacement η of the free surface appears as a main unknown variable. In particular cases one can eliminate the velocity potential ϕ making use of the added mass operator. This operator plays a basic role in the balance energy between the gravity and the capillary. On simple examples the first instability due to the velocity U occurs when the energy of the capillary reach the same level as the one of the gravity (see the example treated J. Lighthill [LIH 78]). But in more general cases, this critical velocity is characterized by an eigenvalue problem. As far as one has to consider a coupling of the sea with a ship, the formulation in η is more appropriate than the one in ϕ. For instance the slamming between the shell of the ship and the surface of the water can be formulated with respect to η. Furthermore the mathematical model of progressive waves formulated in η -even if it is equivalent to the one in ϕ-has better properties for the numerical implementation.

There are other physical mechanisms which should be added to the progressive wave model. For instance the vortex have been eliminated in the potential model, but the viscosity should play a role on the damping of wave and the characterization of the critical velocity U c . The compressibility of the water is certainly acceptable inside the water, but much less for the surface waves. The coupling with the atmosphere is a fundamental problem and it can change a lot the behaviour of the progressive waves, mainly if temperature gradient is meaningful.

Our goal was only to focus on the validity (in fact the non validity) of the Neumann-Kelvin's model) for which a modification is suggested in order to overcome classical difficulties mentioned in many scientific contributions.

Figure 1 :

 1 Figure 1: Main geometrical notations

Figure 2 :

 2 Figure 2: The coarser meshes used for the submarine (left graph) and for the pebble stone bottom (right graph). isovalues of Φ 0 have been plotted for two distinct situations (a submarine and a pebble stone bottom) on figure 3 . It appears that ∇ s Φ 0 is not constant. The variations of ∇ s Φ 0 and ∆ s Φ 0 have been plotted on figure 4.

Figure 3 :

 3 Figure 3: The isovalues of Φ 0 are plotted for the submarine on the left and the one for the pebble stone bottom on the right.

Figure 4 :

 4 Figure 4: On the left are the results for the submarine and on the right for the pebble stone bottom.

Figure 5 :

 5 Figure 5: Convergence of the spectrum of (25) versus the mesh size for the rectangular pool without submarine (the analytical solution is known). The figure 6 corresponds to a pebble stone bottom with three different ratio between the sides of the domain Ω. In both cases the convergence is quite slow. It is much slower than in a classical membrane model for instance and even in a Steklov's one [NEC 67]. The reason is double: first of all the compactness is only based on the compact embedding of the Sobolev's spaces H 1 (Γ f ) into H

Figure 6 :

 6 Figure 6: Convergence of the spectrum of (25) versus the mesh size for the pebble stone bottom (the analytical solution is unknown).

Figure 7 :

 7 Figure 7: Influence of the depth of a submarine on the critical velocities (square) σ = 7, 5 10 -2 for top curve and σ = 0 for bottom curve.

Figure 8 :

 8 Figure 8: Influence of the depth of a pebble stone bottom pool on the critical velocities (square) σ = 7.510 -2 for top curve, σ = 0 for bottom curve.

Figure 10 :

 10 Figure 10: Instable time-evolution of the energy for the linear model for the scheme (39) for several velocities (U > U c and U < U c ): top without capillary, bottom with σ = 7.510 -2

Figure 11 :

 11 Figure 11: Wake above a submarine for different velocities of the steady flow: 0.1 m/s ≤ U ≤ 0.2 m/s 24

Figure 12 :Figure 13 :

 1213 Figure 12: Progressive waves for a small depth: 0.14 m/s ≤ U ≤ 0.24 m/s

Figure 14 :

 14 Figure 14: Non linear solution of the progressive wave model for several velocities U = 0.16m/s, 0.18m/s, 0.2m/s. The simulation concerns the pebble stone model.

η Wave H L d (depth) U e 1


Tilo Mentler
email: t.mentler@inf.hochschule-trier.de

Philippe Palanque
email: palanque@irit.fr

Susanne Boll
email: susanne.boll@uni-oldenburg.de

Chris Johnson
email: c.w.johnson@qub.ac.uk

Kristof Van Laerhoven

Control Rooms in Safety Critical Contexts: Design, Engineering and Evaluation

Keywords: Control room, Pervasive computing, Safety, Usability, UX, Security, Resilience, Dependability

Human-Computer Interaction (HCI) research has been focussing on the design of new interaction techniques and the understanding of people and the way they interact with computing devices and new technologies. The ways in which the work is performed with these interactive technologies has arguably been less of a focus. This workshop aims at addressing this specific aspect of Human-Computer Interaction in the control rooms domain. Control rooms are crucial elements of safety-critical infrastructures (e.g., crisis management, emergency medical services, fire services, power supply, or traffic management). They have been studied in terms of Human-Computer Interaction with respect to routine and emergency operations, human-machine task allocation, interaction design and evaluation approaches for more than 30 years. However, they are dynamic and evolving environments with, for instance, the gradual introduction of higher levels of automation/autonomy. While state of the art control rooms are still characterized by stationary workstations with several smaller screens and large wall-mounted displays, introducing mobile and wearable devices as well as IoT solutions could enable more flexible and cooperative ways of working. The workshop aims at understanding how recent technologies in HCI could change the way control rooms are designed, engineered and operated. This workshop is organized by the IFIP WG 13.5 on Human Error, Resilience, Reliability and Safety in System Development.

1

Overview and goals

While Human Factors approaches [START_REF] Carvalho | Human factors approach for evaluation and redesign of human--system interfaces of a nuclear power plant simulator[END_REF] and User-Centered Design [START_REF] Ulrich | Example user centered design process for a digital control system in a nuclear power plant[END_REF] have been used to design and evaluate control rooms for a long time, the interaction technologies they offer are usually several steps behind what can be found in other areas such as home entertainment or gaming. Control rooms are deployed in very different domains such as crisis management, emergency medical services, intensive care units, fire services, power supply, maritime navigation, or traffic management.

More recent work has been focussing on the use of new technologies in the context of control rooms exploiting new interaction techniques and new interaction technologies introduced in research contributions from the Human-Computer Interaction area. Speech was considered early as an input modality [START_REF] Baber | Feedback requirements for automatic speech recognition in the process control room[END_REF] while the issues raised by its deployment was questioned in [START_REF] Huber | Does speech technology have a place in the control room[END_REF]. Auditory information was used as an alerting system in addition to the traditional display of information in [START_REF] Patterson | Auditory Warning Sounds in the Work Environment [and Discussion[END_REF] even though studies have shown [START_REF] Dehais | Failure to detect critical auditory alerts in the cockpit: evidence for inattentional deafness[END_REF] that human brain filters out that information in case of high workload or stress. Tactile feedback was introduced in cockpits more recently [START_REF] Nojima | Cockpit display using tactile sensation[END_REF] (but this feedback directly used in combination with information displays) while tangible interactions were introduced as a mean to bring interactions closer to the ones on physical knobs [START_REF] Müller | Back to tangibility: a post-WIMP perspective on control room design[END_REF]. Multi-touch interactions have also been recently studied both in control rooms [START_REF] Selim | EGrid: supporting the control room operation of a utility company with multi-touch tables[END_REF] and interactive cockpits [START_REF] Cockburn | Turbulent Touch: Touchscreen Input for Cockpit Flight Displays[END_REF] raising multiple issues. Ambient displays [START_REF] Cobus | Multimodal head-mounted display for multimodal alarms in intensive care units[END_REF] [ref] and head mounted displays [START_REF] Claudius Stratmann | LAME: light-controlled attention guidance for multi-monitor environments[END_REF] [ref] have been used to add a digital layer of information on top of control room elements to support situation awareness and attention in safety critical environments. This introduction of nonstandard (at the time) modalities, led to multimodal interactions (both in terms of input and output) but they are still considered a challenge in 2016 [START_REF] Heimonen | Multimodal interaction in process control rooms: are we there yet?[END_REF] while the early work from Bolt in that domain was introduced more than 25 years earlier [START_REF] Bolt | Put-that-there": Voice and gesture at the graphics interface[END_REF]. Multimodality as input was however introduced successfully in industrial control rooms [START_REF] Fagerlönn | Designing a Multimodal Warning Display for an Industrial Control Room. 12th International Audio Mostly Conference on Augmented and Participatory Sound and Music Experiences (AM '17)[END_REF] or military aircraft cockpits [START_REF] Bastide | A model-based approach for real-time embedded multimodal systems in military aircrafts. 6th int[END_REF] while multimodal output was largely used for alerting flying crews [START_REF] Selcon | Multi-Modal Cockpit Warnings: Pictures. Words. or Both[END_REF].

Beyond interaction technologies, the centralised nature of control rooms was questioned [START_REF] Juhlin | Decentralizing the control room: mobile work and institutional order[END_REF] and mobile solutions were proposed to support control and monitoring activities on the move. However, the introduction of these technologies in critical systems raised immediately the issue of their dependability and their security as soon as the environment loses its closeness nature.

While state of the art in deployed control rooms is still characterized by stationary workstations with several smaller screens and large wall-mounted displays, introducing mobile and wearable devices as well as IoT solutions could enable more flexible and cooperative ways of working. However, turning control rooms into pervasive computing environments raises user-related challenges such as usability and user experience, system-related ones such as reliability and dependability and more global ones such as safety and security [START_REF] Palanque | Beyond usability for safety critical systems: how to be sure (safe, usable, reliable, and evolvable)[END_REF]. However, it is clear that control rooms must evolve and integrate advanced interaction technologies [START_REF] Roth | Integrating digital and conventional human system interface technology: Lessons learned from a control room modernization program[END_REF] so that operations can be improved.

This workshop promotes sharing of experiences in designing, implementing, and evaluating interactive systems in control rooms. We are especially interested by contributions presenting theories, methods, and approaches for considering them as pervasive computing environments and interdependencies between all the properties listed above. In this workshop we consider control rooms in their broader sense including (ship) bridges, cockpits, and operating rooms. This workshop is organized by the IFIP Working Group 13.5 on Human Error, Resilience, Reliability, Safety and System Development.

2

The specificities of control rooms

These control rooms are characterized by their operations that must be resilient to adverse events and unforeseen situations including human errors, system failures, environmental adversary conditions and malicious attacks. One of the specific challenges in their development is the demonstration that they meet requirements imposed by regulatory authorities.

As discussed in [START_REF] Boring | SIG: the role of human-computer interaction in next-generation control rooms. CHI '05 Extended Abst[END_REF] control rooms present very specific aspects that make them very different from other interactive systems:  They allow monitoring and control of complex cyber-physical systems;  Their operations are usually defined and regularly assessed by regulatory authorities;  They are operated by trained professional users that are usually qualified by dedicated organizations;  They usually involve multiple operators who have to coordinate their actions to reach their goals;  Operators have different levels of authority which grant different access to functions and data;  Operator errors might have catastrophic consequences on people and goods;  Failures in the display and control system might have catastrophic consequences on people and goods;  Their deployment requires the acquisition of a certification from external organizations;  They are operated in close worlds reducing security risks. This list is far from being exhaustive but it gives a set of characteristics to be taken into account when considering the design of new types of control rooms and the constraints that have to be addressed.

Target audience and expected outcomes

This workshop is open to everyone who is interested in the aspects related to the design, the engineering, the evaluation, the deployment, the maintenance and the certification of control rooms. We expect a high participation of IFIP working group 13.5 members.

We invite participants to present position papers describing real-life case studies that illustrate, if possible, how a new technology would enhance control room operations. They could also highlight the trade-offs between two or more properties of interactive systems such as user experience and dependability or usability and security. The way the new technology will be addressing some known or envisioned problem in a control should be presented in the contribution. We are also interested in methods, theories and tools for control rooms development as long as they address some user interface properties. Accepted position papers are published in INTERACT 2021 adjunct conference proceedings. We also expect to discuss at the workshop how to disseminate individual contributions to the community in a special issue in a journal or in an edited volume.

Structure of the workshop

This proposal encompasses a full-day workshop organized around presentation of position papers and working activities in small groups. From the set of contributions, a subset of selected case studies is invited to be presented at the beginning of the workshop and is used to support the discussion that follows. The morning sessions are dedicated to welcoming participants and presenting case studies. Participants are invited to comment the case studies and to report similar experiences. The afternoon sessions are devoted to interactive sessions, where participants are engaged to work in small groups to propose solutions to the issues raised by the case studies presented in the morning. Proposed solutions will be compiled and compared. Based on the lessons learned, participants will draft an agenda of future work that can be accomplished.

5 Workshop organizers Tilo Mentler is

  a professor of Human Computer Interaction and User Experience at Trier University of Applied Sciences. His research is focused on human-centered design in safety-critical contexts (e.g. mobile devices and mixed reality in healthcare, novel approaches to critical infrastructure). Currently, he works on control rooms as pervasive computing environments and examines the role of user experience in safetycritical settings. Prof. Mentler is chair of the special interest group "Usable Safety & Security" within the German Informatics Society (GI), member of the IFIP Advisory Board of GI and has been the GI representative in the IFIP Domain Committee on IT in Disaster Risk Reduction. Philippe Palanque is Professor in Computer Science at University of Toulouse 3, where he leads the Interactive Critical Systems research group. Since the late 1980s he has been working on the development and application of formal description techniques for interactive systems. He has worked on research projects at the Centre National d'Études Spatiales (CNES) for more than 10 years and on software architectures and user interface modelling for interactive cockpits in large civil aircraft (funded by Airbus). The main driver of Philippe's research over the last 20 years has been to address in an even way usability, safety and dependability in order to build trustable safetycritical interactive systems. As for conferences he is TPC co-chair of EICS 2021 and is a member of the ACM CHI steering committee. He is a member of CHI academy, the chair of the IFIP TC 13 committee on Human-Computer Interaction, and secretary of IFIP WG 13.5. Susanne Boll is a full professor for Media Informatics and Multimedia Systems at the University of Oldenburg and a member of the board of the OFFIS-Institute for Information Technology. Her research focus lies in the field of human-computer interaction (HCI). In her ongoing research, she is designing novel interaction technology for a respectful and beneficial cooperation of human and technology in a future automated world. She works on novel interaction methods for a safe cooperation and humans in safety-critical automated environments such as automated driving, health care and production. Her scientific results have been published in top venues in her field such as ACM CHI, MobileHCI, AutomotiveUI, as well as internationally recognized journals. Prof. Boll was named a Fellow of the German Informatics Society in October 2020. She was named a Distinguished Member of the ACM in 2019 and is an elected member of acatech, The German National Academy of Science and Engineering. Chris Johnson is Professor and Pro Vice Chancellor of Engineering and Physical Sciences at Queen's University Belfast. Over the last 20+ years, he has authored more than 200 peer reviewed publications, including one of the first textbooks on Accident and Incident Reporting. He has held fellowships from NASA, the US Air Force and Navy. He is a Scientific Advisor to the EC SESAR JU on the future of air traffic management and helped author guidelines on incident reporting for EUROCONTROL and for the European Railway Agency. He also helped the UK Department for Transport to develop the national cyber security strategy for aviation. In 2021, he is one of some fifteen expert witnesses retained to support the public inquiry into the Grenfell Tower fire. Kristof Van Laerhoven is Professor in Ubiquitous Computing at the University of Siegen, Germany. His research interests span the areas of wearable and distributed sensing systems that focus on machine learning challenges, such as recognising what human users are doing, what they are focusing on, and how stressed they are. He is wearable department co-editor at IEEE Pervasive, is co-editor of Springer Adaptive Environments, associate editor for ACM IMWUT, and was general co-chair for ACM UbiComp/ISWC in 2020. More information can be found on http://ubicomp.eti.unisiegen.de .


Lisa Royer

Iryna Makarchuk

Simon Hettler

Raul Arenal

Tristan Asset

Benjamin Rotonnelli

Antoine Bonnefont

Elena Savinova

Benoit P Pichon
email: benoit.pichon@unistra.fr

Prof P Benoit

Pichon

Core-shell Fe3O4@CoFe2O4 nanoparticles as high-performance anode catalysts for enhanced oxygen evolution reaction

Water electrolysis is a promising and environmentally friendly means for renewable energy storage. Recent progress in the development of anion exchange membrane (AEM) opened new perspectives for high-performance anode catalysts based on transition metal oxides (TMOs) for the sluggish anodic oxygen evolution reaction (OER). Here, we report on core-shell nanoparticles (Fe3O4@CoFe2O4) which allow combining an electrocatalytic shell (CoFe2O4) with a conductive core (Fe3O4). Such an original approach significantly minimizes critical Co content in the catalyst and avoid unstable conductive carbon black. The catalyst shows an exceptional OER activity per Co unit mass (2800 A/gcobalt at 1.65 V vs. RHE). Along with the core-shell structure, the size of the Fe3O4 core is a critical parameter. A large conductive Fe3O4 core is beneficial for the enhancement of OER.

Water electrolysis provides an environmentally friendly way to store renewable energies through the H2 production. While proton exchange membrane water electrolysis is well suited to the renewable energy storage, it relies on the utilization of scarce and expensive Ir to accelerate the sluggish kinetics of the anodic oxygen evolution reaction (OER). Thanks to the recent progress in the development of anion-exchange membranes (AEM), transition metal oxides (TMOs) appear as promising substitutes for noble metals. [1][2][3][4] Various approaches were proposed to enhance their OER activity, notably varying the composition, particle size, defect concentration, or, recently, forming core-shell structures. [5][6][7][8][9] Owing to their composition-dependent and widely tunable properties, TMOs with spinel structures became very attractive. While Co3O4 and CoOOH thin films deposited on an Au(111) substrate were reported to significantly enhance the OER activity, 10 core-shell nanoparticles (NPs) combining a gold core and a TMO shell reached 10 000 A/gTM at 1.65 V vs. RHE. Such a high OER activity of TMOs may result from an electronic effect of gold that favors the oxidation of TM up to the 4+ state. 11,12 The use of core-shell nanoparticles might be an efficient way to reduce the amount of expensive and rare electrocatalysts by circumscribing them to a thin shell, 13,14 the exclusive use of cost-effective and abundant materials is the next challenge for the worldwide dissemination of water electrolysis.

TMO NPs are of significant interest since the catalytic activity can be greatly enhanced by adapting their size and shape, i.e. the active surface area. The valence of the metallic cations and their chemical environment are also critical to enhance the OER activity 3,15,16 . Beyond the fine control of these parameters, recent advances in nanoparticle synthesis in liquid media avoid aggregation which also favors high surface area, i.e. enhanced mass-weighted catalytic activity. Consequently, the challenge lies in controlling the formation and the chemical structure of NPs in order to systematically study their electrochemical properties.

Although most of TMOs are insulating, magnetite (Fe3O4) is conductive due to electron hopping between Fe (II) and Fe (III) 17 . Therefore, it may be considered as a promising support material to avoid carbon black, which is often added to circumvent the lack of conductivity of TMOs catalysts, but is unstable under the anode operation conditions. Furthermore, iron cations in TMO phases may result in synergetic effects with the surrounding atoms, thus enhancing the OER 18,19 . Hence, nanostructures combining Fe3O4 as a conductive core and a TMO as catalyst shell would be promising high performance noble metal-free OER catalysts.

Here we report on an original approach to design a new type of carbon-free, noble-metal-free nanostructured material for the OER. The core-shell structure of Fe3O4@CoFe2O4 NPs allows combining a conductive magnetite core and a catalytically active cobalt ferrite shell. Such a cost-effective material resulted in an excellent OER activity per unit mass of Co. The fine control of the core size and the shell thickness led us to investigate the influence of the core-shell structure on the OER activity which is markedly enhanced by a larger Fe3O4 core.

Core-shell NPs were synthesized by the thermal decomposition of metal complexes in hightemperature-boiling solvents (around 300 °C) as we recently reported 20 (see SI for details). Two different approaches were applied: (i) The seed-mediated growth of a CoFe2O4 shell onto the surface of pristine Fe3-O4 nanoparticles which led to samples CS-1, CS-3, and CS-4. (ii) The diffusion of cobalt cations into vacancies at the surface of partially oxidized pristine Fe3-O4 nanoparticles which led to sample CS-2. 21 The core size were varied by adjusting the experimental condition upon synthesis of Fe3-O4 nanoparticles as we reported earlier. 29 The shell thickness was modified by adjusting the amount of Co and Fe metal complexes with a Co/Fe molar ratio of 0.5. An increase of the mount of Co and Fe complexes by 2.5 from CS3 and CS4 resulted in a volume increase of the shell by two times which agree with the larger size of pristine Fe3-O4 nanoparticles used for CS4. For instance, an increase of the Conventional transmission electron microscopy (TEM) images showed faceted-spherical shapes and narrow size distributions both for pristine and for core-shell NPs (Figure 1). Granulometry measurements showed stable colloidal suspensions of non-aggregated NPs thanks to the efficient coating of oleic acid which was added in the reaction medium (Figure S2). Energy-dispersive X-ray (EDX) spectroscopy confirmed the presence of Co between 12 and 19 at. % with respect to Fe in core-shell NPs (Table S1). For CS-1,2,3, these values are much higher than those corresponding to the size variation between pristine and core-shell NPs (Table S1) which is indicative of the diffusion of cobalt in the Fe3-O4 core as we have reported earlier. 22 X-ray diffraction (XRD) patterns recorded for core-shell NPs display peaks that are all indexed to the spinel crystal structure (Figure S3). Although Fe3O4 and CoFe2O4 phases cannot be distinguished because of similar cell parameters (8.396 Å and 8.3919 Å, respectively), 20 cell parameters of coreshell NPs are relatively close to that of Fe3O4 (Table S2). This is ascribed to the low extent of the surface oxidation of pristine nanoparticles which are protected by the CoFe2O4 shell when exposed to air. The high crystallinity was confirmed by lattice fringes all-across core-shell NPs as observed in scanning transmission electron (STEM) micrographs recorded in the high angle annular dark field (HAADF) mode (Figure S4). The local composition of these NPs was further investigated by electron energy-loss spectroscopy (EELS) analysis in STEM (Figure 2 and Figure S6). Spatially-resolved EELS maps show the homogeneous distribution of O in NPs. In contrast, the Co content is the highest at the edge while the Fe content is at its lowest. Cross-section profiles show that the Co fraction reaches a maximum value of 15 -20 % with respect to Fe and O on the edge while it quickly goes down to 5 % (CS-2 and CS-4) and 2 % (CS-3) in the center of the NP, in agreement with a much thinner shell than the core size. Considering stoichiometric CoFe2O4 (14 % at. Co) and the resolution of the measurement (5 Å), Co is certainly distributed as a graded concentration (Co1+xFe2-xO4) which decreases from the NP surface to the center. This result is supported by the variation of Co Content on a longer distance than the size variation calculated from TEM micrographs. A slight increase of O content can also be observed in the edge of nanoparticles which can be correlated to some cationic vacancies and surface defects as usually observed for such nanoparticles. Complementary information on the average near-surface composition of NPs were obtained by performing X-ray photoelectron spectroscopy (XPS) measurements. The variation of the photon energy from 850 to 1350 eV allowed changing the inelastic mean free path of the photoelectron emission (λ) from 0.6 to 1.3 nm depth. Considering that 68 % of the emitted photoelectrons arise from a thickness λ and 99 % emerge from 3 λ, 23,24 the Fe/Co ratio was calculated from the area of Co2p and Fe2p peaks which were recorded for CS-1 and CS-3 (Figure 2 G). While the Fe/Co molar ratio is close to stoichiometric CoFe2O4 at short λ, the Fe content gradually increases with the photon energy. These results confirm that the Co content at the surface of core-shell NPs is much higher than expected for CoFe2O4. It may be ascribed to the decomposition of the Co stearate which happens at a higher temperature than the Fe stearate. 25 According to the La Mer theory, the iron oxide starts growing before cobalt oxide at the surface of pristine NPs. It results in a gradient of the Co concentration in the shell. Nevertheless, both decomposition temperatures (of the iron and cobalt stearate) are close enough to avoid phase segregation since no CoO phase was observed in the XRD patterns 22 (Figure S3). The electrochemical properties of core-shell NPs were studied in order to determine their OER activity. Cyclic voltammograms were recorded after deposition of a thin layer of NPs on a glassy carbon electrode. In order to accurately compare our results with the literature, current-potential curves were normalized in three different ways : to the mass of metal oxide, to the mass of cobalt, and to the NPs surface area (calculated from the TEM mean diameter), see Figure 3 A and S8. In contrast to pristine Fe3-δO4 NPs, core-shell NPs display an important OER activity at 1.65 V vs. the reversible hydrogen electrode (RHE). Hence, the OER activity of NPs unambiguously arises from the cobalt ferrite shell. This was confirmed by Tafel slopes which were extracted from the anodic scans of cyclic voltammograms (Figure 3 B). All the core-shell NPs all display Tafel slopes of about 60 mV/decade which are much lower than those of non-active pristine NPs. 15,26,27 Considering the low NP loading on the electrode, the Co mass-and surface-weighted OER activities were measured at a relatively high potential of 1.65 V vs. RHE (Figure 3C,D). The activity per metal oxide unit mass can be found in Figure S8. While all core-shell NPs exhibit an excellent OER activity, the largest CS-3 and CS-4 NPs display Co-mass weighted-activities (2300 and 2800 𝐴. 𝑔 𝑐𝑜𝑏𝑎𝑙𝑡 -1

, respectively) which are twice those of CS-1 and CS-2 (about 1300𝐴. 𝑔 𝑐𝑜𝑏𝑎𝑙𝑡 -1

). This trend is even clearer when considering the surface-weighted activity which is increased 3-4-fold. It is worth noting that the Coweighted activity of the largest Fe3O4@CoFe2O4 core-shell NPs (CS-4) greatly outperforms that of TMO catalysts, either with uniform or core-shell structures, which were reported in the literature (Table S3). Indeed only CoOx(OH)y NPs 28 and core-shell NPs with a gold core 6,8 (Au@CoFeOx and Au@NiCo2S4) display higher efficiency although they are mixed with unstable carbon black. The difference in activity between CS-1 and CS-2 vs. CS-3 and CS-4 likely originates from the core size (7.5 -9.1 vs. 12.7 -13.6 nm). Indeed, the surface of Fe3-O4 NPs is usually oxidized thus generating an insulating layer assimilated to maghemite 29 . Therefore, a larger core is likely to be essential for efficient OER since it results in a higher fraction of conductive magnetite. 30,31 This hypothesis is supported by the low activity of Fe3O4@CoO, consisting of a 4.5 nm core as reported in Ref. 32 Furthermore, as mentioned in the literature on gold catalysts, we also presume some synergistic electronic effects between Fe and Co cations which may contribute to the OER activity. In this purpose, the OER will be nearly studied in operando by means of soft X-ray spectroscopy.

To sum up, Fe3O4@CoFe2O4 core-shell NPs with a narrow size distribution and uniform shape were successfully synthesized through two alternative approaches (diffusion and crystal growth). They present excellent OER activities per unit mass of Co (from 1300 to 2800 A/gcobalt.) which greatly outperform the OER activities of transition metal oxide nanoparticles reported in the literature. Along with the core-shell structure, the size of the Fe3O4 conductive core and the Co loading at the NP surface are critical parameters for efficient OER activity. Thus, TMO based core-shell nanostructures with fine control of the chemical structure is a viable approach for high performance alkaline OER electrocatalysts.

Figure 1 .

 1 Figure 1. TEM micrographs of core-shell NPs: CS-1 (A), CS-2 (C), CS-3 (E), CS-4 (G) and the corresponding size distributions of the core-shell NPs and their pristine core NPs.

Figure 2 .

 2 Figure 2. Scanning transmission electron microscopy (STEM) high angle annular dark field (HAADF) images of CS-2 (A, B, C) and CS-3 (D, E, F). Spatially-resolved electron energy loss spectroscopy (EELS) analyses (B, E) showing the spatial distributions of oxygen O (red), iron Fe (green) and cobalt Co (blue). Quantitative cross-section profiles (C, F) recorded along the arrows depicted in B, E, respectively. The left y-scale refers to percentage s of Fe, Co and O atoms. Depth profiling analysis performed for different incident energy (G) on CS-1 and CS-3, the photon energy of the X-ray source was varied from 850 to 1350 eV in order to analyse Co2p and Fe2p peaks. Schematic illustration of the nanoparticle depth (H) probed as a function of the inelastic mean free path ().

Figure 3 .

 3 Figure 3. A) Cyclic voltammetry measurements performed on CS-1, CS-2, CS-3 and CS-4 deposited on a glassy carbon electrode (oxide loading below 1µg cm -2 ), in 0.1M NaOH, scan rate 10mV/s. The current has been divided by the oxide mass (B) Tafel slopes of these NPs (in mV/decade). Panels C and D show activity of these NPs at 1.65V vs. RHE per cobalt unit mass (C) and per cm 2 (D).

Acknowledgements

Project financially supported by the Foundation for Frontier Research in Chemistry. We thank the Helmholtz-Zentrum Berlin für Materialien und Energie for the allocation of synchrotron radiation beamtime. The authors acknowledge Juan Velasco Velez for his help on the synchrotron beamline as well as Vinavadini Ramnarain, Barbara Freis and Théo Lucante for the TEM images. STEM and EELS measurements have been performed at the Laboratorio de Microscopías Avanzadas (LMA) at the Universidad de Zaragoza (Spain). R.A. and S.H. acknowledge the funding from the Spanish MICINN (project grant PID2019-104739GB-100/AEI/10.13039/501100011033), from the Government of Aragon (project DGA E13-20R) and from the European Union H2020 program "ESTEEM3" (823717).

Keywords

Nanoparticles, Oxygen Evolution Reaction, Cobalt, Core-shell, Catalyst

Supplementary information

Nanoparticle synthesis, TEM images, XRD patterns, cyclic voltammetry, impedance measurements, EELS and STEM HAADF micrographs

Notes

The authors declare no competing financial interests.


The impact of

published or not. The documents may come

Topic of research

The impact of overexploitation of groundwater resources on the resilience of agricultural farms in semi-arid zones

Nsiri N ¹,²,³ Zaatra R ¹ , Kleftodimos G ¹ , Belhouchette H ¹ , Drogué S ²,³ ¹CIHEAM-IAMM Mediterranean Agronomic Institute of Montpellier, 34093 Montpellier, France

² Montpellier Institut Agro , ³ UMR MoISA

Objective

The main objective of the research is to evaluate the impact of groundwater overexploitation on the resilience of agricultural households in Morocco. We identified 3 farm-types in the area; intensive production system based mainly on vegetables, semi-intensive cereal monoculture households and one perennial crops.

• Water resources in Morocco are rather well known, but limited, irregular, and fragile.

  In order to study the resilence of agricultural farms in the Souss Massa region (Morroco), we used DAHBSIM bio-economic model (Komarek et al. 2017). It is based on mathematical programming methods and maximizes the expected utility of household income.

ResultsReferences:

  Bouchaou et al., 2011/ Hssaisoune et al., 2020/ Komarek et al., 2017 / El Ansari et al., 2020/ Malki et al., 2017

Table 1 :

 1 Income variation with precipitation after simulationOur approach, to analyze the resilience of farmers, focuses on the household level where the main decision-making is taking place.

 South of the Atlas mountains  Semi-arid to arid climate  Average rainfall of 200 to 250 mm/y  Quasi-absence of surface water  Importance of groundwater resources  Water consuming activities

.73 112,433,980Table 1 :

 1 Farm income and pumping costs

		Mm3/y
	• The expansion of irrigated	100
	agricultural increased the groundwater land has use, resulting in the overexploitation of local aquifers. • Water scarcity is expected to have a negative impact on food production and threaten the resilience of the local agricultural system	50 0 -50 7372-100 Indicator Scenario Business As Average 1972 1996 1998 2003 2007 2015 Figure 1 : Evolution of groundwater withdrawal and water balance of reference (Sc_REF) 2020 Usual (BAU) (2035) annual cost of degradation Water cost(dh/ m3) 26.57 61.93 40.85 43.7 Pumping costs (dh/m) 1827.53 2268.65 388.01 441.12 Farm 34243.8 26871.14 3686.36 Bilan Irrigation Pumping Drinking water and irrigation pumping	Cost of degradation Sc_REF -BAU 6,536,000,000 -	)
		incom	7
		(dh/farm)