David Imbert
email: david.imbert.1@free.fr

Sean Mcnamara
email: sean.mcnamara@univ-rennes1.fr

Yves Le

Y Le Gonidec
email: yves.legonidec@univ-rennes1.fr

Fictitious domain method for acoustic waves through a granular suspension of movable rigid spheres

Keywords: granular medium, suspension, fictitious domain method, distributed Lagrange multipliers, movable rigid grains, acoustic wave equation

come

Introduction

In marine Geosciences, acoustic probing is commonly used to detect subsurface reflectors, characterize their physical properties and study geophysical processes such as slope instabilities. In particular, undersea sediment layers constitute granular materials which pore space is filled with waver, i. e. grains are movable bodies. Such layers are acoustically complex materials where sound waves can be multi-scattered and reflected at fluid-grain interfaces. Experiments with model spherical grains under water have been performed to better understand the acoustic behavior of granular media in equilibrium [START_REF] Le Gonidec | Multiscale analysis of waves reflected by granular media: Acoustic experiments on glass beads and effective medium theories[END_REF] or tilted under gravity [START_REF] Zaitesev | Pre-avalanche structural rearrangements in the bulk of granular medium: Experimental evidence[END_REF][START_REF] Thirot | Acoustic emissions in granular structures under gravitational destabilization[END_REF]. These experiments clearly highlight the contribution of the fluid-grain interaction on the acoustic measurements, which can also be affected by inter-particle contacts [START_REF] Tournat | Probing weak forces in granular media through nonlinear dynamic dilatancy: clapping contacts and polarization anisotropy[END_REF].

In the case of pores filled with air, dry granular materials are commonly modeled using discrete element methods such as "Molecular Dynamics" [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF][START_REF] Mcnamara | Méthode Dynamique Moléculaire[END_REF]. Each grain is represented using a rigid and incompressible sphere completely characterized by the position of its center and its radius. The air filling the pores can be neglected because its density is very low, so that the momentum it transmits is negligible. In the case of pores filled with water, the momentum transmitted by the fluid can no longer be neglected. Thus, discrete element methods can not be used to model submerged granular media. Well known numerical models exist to model sound waves in a liquid: for example finite differences [START_REF] Alford | Accuracy of finite-difference modeling of the acoustic wave equation[END_REF][START_REF] Kosloff | Seismic numerical modeling[END_REF] or finite elements [START_REF] Zienkiewicz | The finite element method in engineering science[END_REF][START_REF] Kosloff | Seismic numerical modeling[END_REF]. But to our knowledge, none of them usually consider neither the coupling between the wave and the motion of grains, nor the inter-particles contacts.

The fictitious domain method based on Lagrange multipliers was introduced in [START_REF] Babuška | The finite element method with Lagrange multipliers[END_REF] and started to be popular two decades ago when applied to Navier-Stokes equations [START_REF] Glowinski | A Lagrange multiplier / fictitious domain method for the Dirichlet problemgeneralization to some flow problems[END_REF] and then for scattering problems in acoustics [START_REF] Heikkola | Fictitious domain methods for the numerical solution of two-dimensional scattering problems[END_REF]. The particularity of this method is to work with two independent finite element meshes: the first one regular and time-independant and the second one unstructured and possibly time-dependant. This last attribute was first applied for forced displacement of particles in [START_REF] Glowinski | Fictitious domain methods for incompressible viscous flow around moving rigid bodies[END_REF] and to model vibrations of the membrane of a kettledrum [START_REF] Rhaouti | Time-domain modeling and numerical simulation of a kettledrum[END_REF]. Introduction of distributed Lagrange multipliers [START_REF] Glowinski | On the solution of the Dirichlet problem for linear elliptic operators by a distributed Lagrande multiplier method[END_REF] allowed simulation of freely movable particles through a rigid body motion constraint applied inside particles [START_REF] Glowinski | A distributed Lagrange multiplier / fictitious domain method for particulate flows[END_REF]. The fictitious domain method based on distributed Lagrange multipliers was also developed in [START_REF] Bokil | An operator splitting scheme with a distributed Lagrange multiplier based fictitious domain method for wave propagation problems[END_REF] for the mixed formulation of the acoustic wave equation with Dirichlet conditions.

In this work, we unite a finite element method to a discrete element model through the fictitious domain method, yielding a single method capable of describing acoustic propagation by both the fluid and inter-granular contacts.

To our knowledge, our model introduced in [START_REF] Imbert | Fictitious domain method to model a movable rigid body in a sound wave[END_REF] is the first to couple the acoustic wave equation governing the fluid to the dynamic equations governing the motion of rigid movable grain. However, in that work, we simulated a single rigid sphere in a sound wave, and we compared only its velocity to the analytic solution. In addition to results on the dynamic of the grain, our model allows each grain to generate its own acoustic field when it is under a force due to another incident wave but also due to any other forces acting on it such as a contact force or the gravity. These resulting acoustic fields have never been validated and confronted to analytic solutions.

In the present paper, we present a major development of the model introduced in [START_REF] Imbert | Fictitious domain method to model a movable rigid body in a sound wave[END_REF]. First of all, we take on a new mixed finite elements method based on standard conforming finite elements to fix severe convergence issues of the fictitious domain method. We also change the numerical scheme to improve energy conservation. Finally, the model is now able to work with several grains and we illustrate that with new numerical experiments.

The paper is subdivided into 2 main parts. In section 2 and 3, we are describing the model, i. e. governing equations, the fictitious domain method, and the computational scheme. Section 4 deals with three numerical experiments, involving one, two and many movable grains submerged in water, respectively. With the first experiment, both dynamic and acoustic numerical results are compared with analytic solutions of acoustics and an estimation of the error is given. With the second experiment, we analyse the energy transfers between two oscillating grains separated by water and linked with their own springs. In the last experiment, we highlight the potential of the model to sound a suspension of grains with acoustic pulses.

Ω Γ 1 Γ 2 Γ 3 Γ 4 G 1 ∂G 1 G 2 ∂G 2 G 3 ∂G 3 n 1 t 1 n 2 t 2 n 3 t 3 e x e y n Γ
Figure 1: System with the different domains: domains G k are the domain of each grain k ; and the hatched rectangle is the whole domain Ω and Γ its boundary.

Governing equations

A schematic representation of rigid bodies in a fluid matrix is considered on Fig. 1. Let Ω be the domain which includes a set of grains G and a fluid matrix Ω \ G. The domain G is subdivided into N G circular subdomains G k , one for each grain, such that G = ∪ G k . Let Γ be the boundary of the whole domain Ω, ∂G k the boundary of each domain G k , and therefore let ∂G be the boundary of the domain G. We also define n k the outgoing normal unit vector of the boundary ∂G k , t k the tangential unit vector, and X k the position of the center of the subdomain G k .

We consider that the fluid is perfect, i.e. non viscous, non heat conducting and isentropic. The perfect fluid approximation is commonly used in acoustics because the scatteering and transport of energy are much more important than dissipation due to viscosity. Moreover, the variations of the pressure fields are small compared to their mean value and the density of the fluid ρ 0 and the velocity of sound in the fluid c 0 are constant. These assumptions allow use of the mixed formulation of the acoustic wave equation governing the fluid.

Strong forms

We can write the acoustic wave equation as a system of two equations in terms of pressure field p and velocity field u:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ρ 0 ∂u ∂t + ∇ p = 0 in Ω \ G 1 ρ 0 c 2 0 ∂p ∂t + ∇ • u = 0 in Ω \ G (1a) (1b)
The second system of equations governs the motion of grains:

ρ S k dU k dt = - ∂G k p n k d + F k ∀ k = 1, 2, ..., N G (2)
Here S k and U k are the surface (in a 2D domain) and the velocity of the grain k, respectively. All grains have the same density ρ k = ρ ∀ k = 1, 2, ..., N G . F k can handle any other extra force, for example if the grain is subjected to the gravity, a spring force, or some interaction forces. Contrary to system (1), Eq. ( 2) is discrete and therefore allows use of discrete elements methods such that the Molecular Dynamics algorithm [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF].

Then, to link systems (1) and ( 2), we add the natural jump condition of the acoustic wave equation for movable rigid bodies:

u • n k = U k • n k on ∂G k ∀ k = 1, 2, ..., N G (3)
Eq. ( 3) is essential to couple the grains velocity U k and the acoustic velocity field u.

Eq. ( 3) is just the "slip" boundary condition appropriate for the perfect fluid we are considering, applied to the surface of the grains. The particle rotations do not appear because they couple to the fluid only through tangential viscous stresses exerted by the viscosity of the particle surface. Our fluid is inviscid and therefore this coupling is absent.

Note that the jump condition (3) does not prevent particle rotation. In the special case where F k includes friction forces between grains, then an additional term including angular velocity should be added to the first member of Eq. ( 2). But even in this case, rotations of a sphere in a perfect fluid can not generate any sound.

Another natural consequence of Eq. ( 3) is that only the normal component of the acoustic velocity field is constrained at the interface. Because of the perfect fluid, we have no information on the tangential component of the velocity field. This is consistent with the governing equation Eq. (1) which allows more freedom for the velocity and pressure fields than the Navier-Stokes equation would do.

Two kinds of boundary conditions are set on the outer frontier Γ = ∪ Γ j of Ω. On one hand, Perfectly Matched Layers (PML) [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] are set on the bottom side Γ 1 and the top side Γ 3 of the domain Ω but we will not focus on that point to keep equations simple. See [START_REF] Chew | A 3D perfectly matched medium from modified Maxwell's equations with stretched coordinates[END_REF] for details on the construction of the PML based on a coordinate stretching approach for Maxwell's equations and [START_REF] Bokil | An operator splitting scheme with a distributed Lagrange multiplier based fictitious domain method for wave propagation problems[END_REF] for the acoustic wave equation. These layers end with a Dirichlet condition applied on the pressure field:

p = 0 on Γ 1 ∪ Γ 3 (4)
On the other hand, we set periodic boundary conditions on the right side Γ 2 and the left side Γ 4 of the domain Ω

p| Γ 2 = p| Γ 4 on Γ 2 ∪ Γ 4 (5) u| Γ 2 = u| Γ 4 on Γ 2 ∪ Γ 4 (6)
In that way, we can easily generate plane waves which are useful when comparing with analytic solutions of acoustics.

Weak forms

We write now the weak form of Eqs. (1a) and ( 2) into a single equation with the unknowns

u ∈ H(div, Ω \ G), U k ∈ R 2 ∀ k = 1, ..., N G and p ∈ L 2 (Ω \ G): Ω\G ρ 0 ∂u ∂t • v dx - Ω\G p∇ • v dx + Γ p n Γ • v d - ∂G p n • v d + k ρ S k dU k dt • V k -F k • V k + ∂G k p n k • V k d = 0 ( 7 )
Here v ∈ H(div; Ω \ G) 2 , V k ∈ R 2 , and q ∈ L 2 (Ω \ G) are the respective test functions of u, U k , and p. Eq. ( 7) was called combined equation of motion in [START_REF] Hu | Direct simulation of flows of solid-liquid mixtures[END_REF] and it expresses the conservation of the momentum in the whole domain Ω. We note that the integral on ∂G (last term of the 1st line of ( 7)) and the sum of hydrodynamic forces (last term of the 2nd line of ( 7)) cancel each other out if and only if we extend Eq. ( 3) to the test functions and require v

• n k = V k • n k ∀ k = 1, 2, ..., N G .
To take this condition into account, we define a new functional space WΩ\G that we also called combined variation space in reference to the space introduced in [START_REF] Glowinski | A distributed Lagrange multiplier / fictitious domain method for particulate flows[END_REF] for the Navier-Stokes equation:

WΩ\G = (v, V k ) v ∈ H(div; Ω \ G), V k ∈ R 2 ; v • n k = V k • n k on ∂G k ∀k = 1, 2, ..., N G (8)
Moreover, boundary conditions ( 4) and ( 5) make the integral on Γ vanish, and therefore lead to find (u, U k ) ∈ WΩ\G and p ∈ L 2 (Ω \ G) such as:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Ω\G ρ 0 ∂u ∂t • v dx - Ω\G p ∇ • v dx + k ρ S k dU k dt • V k -F k • V k = 0 Ω\G 1 ρ 0 c 2 0 ∂p ∂t q dx + Ω\G ∇ • u q dx = 0 (9a) (9b) ∀ (v, V k ) ∈ WΩ\G and q ∈ L 2 (Ω \ G).

Fictitious domain method formulation

The fictitious domain method consists in extending Eq. ( 9) defined in the fluid domain Ω \ G to the domain Ω such that the acoustic fields p and u are also defined in the whole regular domain Ω. To do so, we proceed in three steps :

1. first, we extend the combined equation of motion (9a) to the whole domain Ω; 2. then, constraints acting on the grain are enforced in a weak sense using Lagrange multipliers; 3. finally, we extend the linearized continuity equation (9b) to the whole domain using the same constraints

Combined equation of motion for the extended domain

First, we follow the method described in [START_REF] Glowinski | A distributed Lagrange multiplier / fictitious domain method for particulate flows[END_REF] to develop an equation analogous to Eq. (9a) valid in the domain G k . The extension from one grain G k to the whole domain G = G k is straightforward. However, we start with a different hypothesis from [START_REF] Glowinski | A distributed Lagrange multiplier / fictitious domain method for particulate flows[END_REF] because we use a different boundary condition on the grain, namely Eq. (3).

Let (X k , e r , e θ ) be an orthonormal basis with the center of the grain X k as its origin and e r = n k and e θ = t k orthogonal unit vectors. The following development to find an equation for the extended domain G k is split into two steps: one to find an expression involving the projection of u on e r and another one for the projection of u on e θ .

Firstly Eq. ( 3) is extended from ∂G k to G k :

u • e r = U k • e r in G k (10)
Let's call u r the projection of u on e r :

u r = U k • e r ( 11
)
We can derive Eq. ( 11) with respect to the time and multiply each member by ρ 0 :

ρ 0 ∂u r ∂t = ρ 0 , d dt (U k • e r ) in G k (12)
We introduce the combined variation space WG k , counterpart of the space WΩ\G in G k :

WG k = (v, V k ) v ∈ H(div; G k ), V ∈ R 2 ; v • n k = V k • n k on ∂G k (13)
The weak form of Eq. ( 12) can be written:

G k ρ 0 ∂u r ∂t v r dx = G k ρ 0 d dt (U k • e r ) V k r dx (14)
Then, knowing that e r = cos(θ) e x + sin(θ) e y : [START_REF] Glowinski | On the solution of the Dirichlet problem for linear elliptic operators by a distributed Lagrande multiplier method[END_REF] we compute integrals in the right hand side of (15):

G k ρ 0 ∂u r ∂t v r dx = ρ 0 R k 0 dr 2 π 0 d dt U k x cos(θ) + U k y sin(θ) V k x cos(θ) + V k y sin(θ) dθ
G k ρ 0 ∂u r ∂t v r dx = ρ 0 S k dU k dt • V k dx (16)
Now, we do the same for the tangential component of the velocity field u θ = u • t k (see Fig. 1 for the direction of t k ). We begin with another natural jump condition of the wave equation derived from the continuity of the pressure through the interface ∂G k [START_REF] Lombard | Modélisation numérique de la propagation des ondes acoustiques et élastiques en présence d'interfaces[END_REF]:

ρ 0 ∂ ∂t (u • t k ) = ∇ p • t k on ∂G k ( 17
)
As described above, we extend Eq. 17 to G k , multiply by the test function v θ and we integrate using the gradient theorem on a closed surface so that:

G k ρ 0 ∂u θ ∂t v θ dx = 0 (18)
Injecting Eq. ( 18) in Eq. ( 16) leads to:

G k ρ 0 ∂u ∂t • v dx = ρ 0 S k dU k dt • V k dx ∀(v, V k ) ∈ WG k ( 19
)
To go further, we enforce a divergence free constraint of the velocity field inside the grain G k to prevent propagation of waves through the grain:

∇ • u = 0 in G k ( 20
)
The associated constraint ∇ • v = 0 related to Eq. ( 20) has to be joined to the functional space WG k [START_REF] Glowinski | Fictitious domain methods for incompressible viscous flow around moving rigid bodies[END_REF]:

W G k = (v, V k ) v ∈ H(div; G k ), V k ∈ R 2 ; v • n k = V k • n k on ∂G k ; ∇ • v = 0 in G k ( 21
)
In that way, we can add the missing integral to obtain the analog of Eq. (9a) and the problem is now to find (u, U k ) ∈ W G k such as:

G k ρ 0 ∂u ∂t • v dx - G k p ∇ • v dx = ρ 0 S k dU k dt • V k dx ∀(v, V k ) ∈ W G k (22)
Then, performing the symbolic operation "(9a) + k [START_REF] Lombard | Modélisation numérique de la propagation des ondes acoustiques et élastiques en présence d'interfaces[END_REF]" yields to the combined equation of motion for the extended domain with solutions (u,

U k ) ∈ W Ω : Ω ρ 0 ∂u ∂t • v dx - Ω p ∇ • v dx + k (ρ -ρ 0 ) S k dU k dt • V k -F k • V k = 0 ∀(v, V k ) ∈ W Ω (23)
with W Ω the combined variation space extended to Ω:

W Ω = (v, V k ) v ∈ H(div; Ω), V k ∈ R 2 ; v • n k = V k • n k on ∂G k ; ∇ • v = 0 in G k ∀ k = 1, ..., N G (24)

Enforcement of constraints using Lagrange multipliers

In a second phase, the two constraints v • n = V • n and ∇ • v = 0 have to be extracted from the combined variation space [START_REF] Glowinski | Distributed Lagrange multiplier methods for incompressible viscous flow around moving rigid bodies[END_REF] and added explicitly to Eq. ( 23) with Lagrange multipliers. For each unknown Lagrange multiplier added, an equivalent equation have to be joined to Eq. [START_REF] Patankar | A new formulation of the distributed Lagrange multiplier / fictitious domain method for particulate flow[END_REF]. The jump condition of the wave equation ( 3) is enforced using boundary Lagrange multipliers λ ∂G on ∂G and distributed Lagrange multipliers [START_REF] Glowinski | On the solution of the Dirichlet problem for linear elliptic operators by a distributed Lagrande multiplier method[END_REF] λ G are used to enforce the divergence free constraint [START_REF] Chew | A 3D perfectly matched medium from modified Maxwell's equations with stretched coordinates[END_REF]

inside G. The problem is to find u ∈ H(div, Ω), U k ∈ R 2 , λ ∂G ∈ H 1/2 (∂G), λ G ∈ H 1 (G) such as: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Ω ρ 0 ∂u ∂t • v dx - Ω p ∇ • v dx + k (ρ -ρ 0 ) S k dU k dt • V k -F k • V k + k ∂G k λ ∂G (v • n k -V k • n k ) d + G λ G ∇ • v dx = 0 for k = 1, 2, ..., N G ∂G k (u • n k -U k • n k ) μ ∂G d = 0 G ∇ • u μ G dx = 0 ∀v ∈ H(div, Ω), V k ∈ R 2 , μ ∂G ∈ H 1/2 (∂G), μ G ∈ H 1 (G) (25a) (25b) (25c)
with μ G and μ ∂G test functions of λ ∂G et λ G respectively.

Linearized continuity equation for the extended domain

Now, let us include Eq. (9b) into the system [START_REF] Hecht | FreeFem++, Laboratoire Jacques-Louis Lions[END_REF]. We write the weak form of the divergence free constraint [START_REF] Chew | A 3D perfectly matched medium from modified Maxwell's equations with stretched coordinates[END_REF]:

G ∇ • u q dx = 0 ∀q ∈ L 2 (G) (26)
Note that Eq. ( 26) is redundant with Eq. (25c). Consequently, we can directly extend Eq. (9b) to Ω because Eq. ( 25c) is enforced anyway. The final problem consists in finding u ∈ H(div, Ω),

U k ∈ R 2 , p ∈ L 2 (Ω), λ ∂G ∈ H 1/2 (∂G), and λ G ∈ H 1 (G) such as: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Ω ρ 0 ∂u ∂t • v dx - Ω p ∇ • v dx + k (ρ -ρ 0 ) S k dU k dt • V k -F k • V k + k ∂G k λ ∂G (v • n k -V k • n k ) d + G λ G ∇ • v dx = 0 Ω 1 ρ 0 c 2 0 ∂p ∂t q dx + Ω ∇ • u q dx = 0 for k = 1, 2, ..., N G ∂G k (u • n k -U k • n k ) μ ∂G d = 0 G ∇ • u μ G dx = 0 ∀v ∈ H(div, Ω), V k ∈ R 2 , q ∈ L 2 (Ω), μ ∂G ∈ H 1/2 (∂G), μ G ∈ H 1 (G) (27a) (27b) (27c) (27d)
Eq. ( 27) can be solved only if ρ 0 ρ, otherwise the unknown disappears from Eq. (27a). A similar case occurs in [START_REF] Glowinski | A distributed Lagrange multiplier / fictitious domain method for particulate flows[END_REF] for Navier-Stokes equations but another fictitious domain formulation was found to deal with this specific issue in [START_REF] Patankar | A new formulation of the distributed Lagrange multiplier / fictitious domain method for particulate flow[END_REF]. We do not treat the case where the density of the fluid ρ 0 is the same as the density of the grain ρ because it is not relevant in the case of our geophysical applications. Moreover, [START_REF] Patankar | A new formulation of the distributed Lagrange multiplier / fictitious domain method for particulate flow[END_REF] does not refer to an explicit velocity of the grain U k which we need for the discrete algorithm we use.

Computational scheme

Finite elements

With the fictitious domain approach, two independent meshes have to be considered: one for the global domain Ω and one for the domain of grains G. A regular grid T Ω h (see Fig. 2a) is used for the rectangular domain Ω to allow use of optimization techniques for the computation of the pressure and velocity fields p and u. The second grid T G h (Fig. 2b) refers to the domain of grains G: this grid may be unstructured but requires specific quadrature formula to perform the numerical integration on triangles. So we prefer to choose a pseudo-regular grid as in [START_REF] Glowinski | Distributed Lagrange multiplier methods for incompressible viscous flow around moving rigid bodies[END_REF] to keep variational spaces simple. Furthermore, T G h is independent from the grid T Ω h and can be shifted during the simulation. However, the size of the two meshes cannot be chosen independently; there is a compatibility condition that must be satisfied so that the fictitious domain method converges [START_REF] Babuška | The finite element method with Lagrange multipliers[END_REF]. Concretely, we choose in our experiments a mesh ratio of κ ≈ 1.2 which yields results in good agreement with analytic solutions.

We use the free software FreeFEM++ [START_REF] Hecht | FreeFem++, Laboratoire Jacques-Louis Lions[END_REF] to generate meshes and solve the partial differential equations. In fact, FreeFEM++ works in the same way as Matlab does, i.e. with its own highlevel programming scripting language, close to the C language. So it is easy to implement both finite elements and the discrete element method governing the motion of grains. The full access to matrices also allows them to be modified to include constraints of the fictitious domain method.

In [START_REF] Imbert | Fictitious domain method to model a movable rigid body in a sound wave[END_REF], we based our finite element approximation on the lowest order Raviart-Thomas element RT 0 [START_REF] Raviart | A mixed finite element method for 2-nd order elliptic problems[END_REF] which is the natural choice to discretize the functional space H(div; Ω). Despite the validation of the particle velocities in [START_REF] Imbert | Fictitious domain method to model a movable rigid body in a sound wave[END_REF], the examination of the acoustic fields pointed out serious issues with RT 0 element. Indeed, convergence issues were also revealed in [START_REF] Bécache | Convergence results of the fictitious domain method for a mixed formulation of the wave equation with a Neumann boundary condition[END_REF] while using the fictitious domain method with Neumann boundary condition and the RT 0 element. We solve these problems in the present paper by opting for mixed methods using standard conforming finite elements [START_REF] Li | Mixed methods using standard conforming finite elements[END_REF]. Lagrangian finite elements are used to approximate the solutions: the piece-wise constant element P 0 h for the pressure space L 2 (Ω) is defined as:

P 0 h = q h ∈ L 2 (Ω h ) ∀K ∈ T Ω h , q h | K = β K ( 28
)
with K a triangular finite element. The continuous piece-wise linear element P 1 h is used to discretize the velocity space H(div; Ω) in T Ω h :

P 1 h = v h ∈ C 0 (Ω h ) | ∀K ∈ T Ω h , v h | K ∈ P 1 ( 29
)
and also the Lagrange multipliers spaces H 1 (G) in T G h and H 1/2 (∂G) on T ∂G h .

Finite-differences time-domain method 3.2.1. Operator splitting scheme

System ( 27) is too intricate to be solved directly using a standard finite difference scheme. It has to be separated into smaller subproblems using the Marchuk's first order operator splitting scheme [START_REF] Marchuk | Splitting and alternating direction methods[END_REF]. Consider the problem :

∂ϕ ∂t + A (ϕ) = f and ϕ| t=0 = ϕ 0 ( 30
)
where

A = N α=1 A α , A α ≥ 0, N ≥ 2 and A is time-independent.
The absolutely stable algorithm is, for j = 0, 1, ..., N:

ϕ j+ 1 N -ϕ j Δt + A 1 ϕ j+ 1 N = f j+1 1 . . . ϕ j+1 -ϕ j+ N-1 N Δt + A N ϕ j+1 = f j+1 N ( 31
)
with Δt the time step. To apply this algorithm, the system ( 27) is subdivided into two steps discussed in detail below:

1. propagate the acoustic wave in the fluid (step 1a); compute position, velocity and forces acting on the grains (step 1b), 2. enforce constraints applied on the grains (step 2).

We compute the acoustic velocity field and grain velocities at the same time because they are linked through the same velocity space, so they both correspond with step 1 of Marchuk's algorithm. However, resolutions of the acoustic wave equation and the equation of motion are totally independent inside step 1.

Step 1a: wave equation with semi-implicit Euler

Firstly, we solve the wave equation system by finding u h ∈ P 2 1 h and p h ∈ P 0 h such as:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Ω 1 ρ 0 c 2 0 ∂p h ∂t q h dx + Ω ∇ • u h q h dx = 0 Ω ρ 0 ∂u h ∂t • v h dx - Ω p h ∇ • v h dx = 0 (32a) (32b)
∀ v h ∈ P 2 1 h and q h ∈ P 0 h . We can solve successively equations (32a) and (32b) using the following first order semi-implicit Euler scheme :

1. find p n+ 1 2 h ∈ P 0 h such as: Ω 1 ρ 0 c 2 0 p n+ 1 2 h -p n h Δt q h dx + Ω ∇ • u n h q h dx = 0 ∀ q h ∈ P 0 h (33) 2. find u n+ 1 2 h
such as:

Ω ρ 0 u n+ 1 2 h -u n h Δt • v h dx - Ω p n+ 1 2 h ∇ • v h dx = 0 ∀ v h ∈ P 2 1 h ( 34
)
with t = Δt min . In the semi-implicit Euler scheme, the minimal time step Δt min is related to the minimal spatial step h min by:

Δt min ≤ α c 0 h min ( 35
)
where α is a positive number.

Step 1b: equation of motion with Verlet algorithm

Independently to the computation of the wave equation ( 33) and ( 34), we can solve the part of the system corresponding to the equation of motion. Find U k ∈ R 2 ∀k = 1, 2, ..., N G such as:

(ρ -ρ 0 ) S k dU k dt • V k -F k • V k = 0 ∀ V k ∈ R 2 (36)
We use the Molecular Dynamics algorithm [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF][START_REF] Mcnamara | Méthode Dynamique Moléculaire[END_REF] to compute forces, integrate the equation of motion for each grain and govern interactions between rigid grains. We choose to solve the equation of motion in three parts using the Verlet algorithm [START_REF] Verlet | Computer "experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules[END_REF][START_REF] Swope | A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters[END_REF].

1. First, position of the grain X

n+ 1 2 k
is computed:

X n+ 1 2 k = X n k + Δt U n k + Δt 2 A n k (37)
2. Secondly, all forces F depending on the positions are computed:

A n+ 1 2 k = F X n+ 1 2 k (ρ -ρ 0 ) S k (38)
3. Finally, the velocity is deduced from the acceleration values:

U n+ 1 2 k = U n k + Δt 2 A n+ 1 2 k + A n k (39)

Step 2: enforcement of constraints with conjugate gradient algorithm

In the second time step of Marchuk's operator splitting scheme, the jump condition of the acoustic wave equation is applied on the boundary ∂G and the rigid body motion constraint is applied inside the grain G. This leads to find u n+1 h

∈ P 2 1 h , U n+1 k ∈ R 2 , λ n+1 ∂G h ∈ P 1 h and λ n+1 G h ∈ P 1 h such as: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Ω ρ 0 u n+1 h -u n+ 1 2 h Δt • v h dx + k ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ (ρ -ρ 0 ) S k U n+1 k -U n+ 1 2 k Δt • V k ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ + k ∂G k λ n+1 ∂G h (v h • n k -V k • n k ) d + G λ n+1 G h (∇ • v h ) dx = 0 N G times ∂G k u n+1 h • n k -U n+1 k • n k μ ∂G h d = 0 G ∇ • u n+1 h μ G h dx = 0 (40a) (40b) (40c) ∀ v h ∈ P 2 1 h , V k ∈ R 2 , μ ∂G h ∈ P 1 h and μ G h ∈ P 1 h .
The discretized system (40) can also be represented as matrices:

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ A 1 0 C T 1 C T 3 0 A 2 C T 2 0 C 1 C 2 0 0 C 3 0 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ u n+1 U n+1 λ ∂G λ G ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ u n+ 1 2 U n+ 1 2 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (41)
To simplify the notation, we introduce block matrices A, C; and vectors b, d, and λ:

A = A 1 0 0 A 2 C = C 1 C 2 C 3 0 b = u n+1 U n+1 d = u n+ 1 2 U n+ 1 2 λ = λ ∂G λ G ( 42
)
and we rewrite system (41) as:

A C T C 0 b λ = d 0 (43)
Such system can be solved using Uzawa's method for the conjugate gradient algorithm [START_REF]Studies in linear and non-linear programming[END_REF]. We initialize the algorithm with an initial approximation λ 0 and we construct an approximation of the solution b for i = 0 until convergence by successively solving:

A b i+1 = d -C T λ i (44)
and computing:

λ i+1 = λ i + β C b i+1 (45)
with β a scalar. The optimal rate of convergence of the algorithm is given by [START_REF] Elman | Inexact and preconditioned Uzawa algorithms for saddle point problems[END_REF]:

β = 2 μ 1 + μ 2 (46)
with μ 1 and μ 2 the minimum and maximum eigenvalues of the Schur complement of sytem (43) C A -1 C T .

Velocity integration

Independently of Marchuk's algorithm and after solving the system (40b), we integrate the final velocity U n+1 k of the grains to get their next position. We simply integrate the velocity for the whole time step to retrieve updated position of grains:

X n+1 k = X n k + Δt 2 U n+1 k + U n k (47)
Then the mesh T G h is updated with the new positions X n+1 k without modifying the regular mesh T Ω h (Fig. 2).

Numerical experiments

In this part, we apply the method developed in sections 2 and 3 on three examples. On one hand, one movable grain is subject to an incident plane wave, and numerical results on dynamics and the acoustics are compared to analytic solutions. On the other hand two movable grains fixed by springs but without any incident plane wave illustrates energy transfers between fluid and grains. Finally we test a more realistic configuration of a random granular suspension with hundreds of grains to illustrate the possibilities of the model.

Movable grain in sound

To study the error compared to an analytic solution, we send a plane wave on a movable grain (Fig. 3). We first define a harmonic plane wave coming from the top (with θ i = -π/2 the incident angle) and defined as:

p i (r, θ, t) = A e i(k r cos(θ+θ i ) -ω t) (48)
with r the distance from the center of the grain and θ the observation angle. The pulsation ω and the wavenumber k are related to the frequency by f = ω 2 π = k c 0 2 π . The diameter of the grain is d = 2 R = 1 mm. At the beginning of the experiment, the center of the grain is located at coordinates (5; 4.5) mm. The size of the domain Ω is 10 × 10 mm and the size of PML layers on the top and bottom is 5 h Ω . The source line is located at y src = 8.3 mm and there are PML layers of thickness 0.4 mm on the top and on the bottom of the domain.

The wave velocity in the fluid medium is c 0 = 1,500 m/s. The density of the fluid medium is ρ 0 = 1,000 kg/m 3 and the density of the material of grains is: ρ = 2,500 kg/m 3 . The frequency of the incident wave is f = 1.5 MHz which corresponds to a of λ/d = 1. The amplitude of the incident wave is A = 1.5 Pa in terms of pressure or 1 μm/s in terms of velocity.

We discretize the domain Ω with N 2 = 120 × 120 nodes and so h Ω = 1/12 mm and the maximum mesh size of the fictitious domains satisfy h G ≈ 1.2 √ 2 h Ω . We set parameter α of Eq. ( 35) related to the CFL condition of Euler semi-implicit scheme for the time discretization to α = 0.65/ √ 2. The duration of the simulation is t f = H/c 0 with H the height of the domain Ω.

Pressure and velocity fields

The reference analytic solution for the pressure field scattered by a movable rigid disc can be written in cylindrical coordinates as [START_REF] Potel | Acoustique Générale -Équations différentielles et intégrales, solutions en milieux fluide et solide, applications[END_REF]:

p s (r, θ) = ∞ n=0
B n H (1) n (k r) cos (n (θ

-θ i )) ∀r ≥ R, 0 ≤ θ ≤ 2 π ( 49
)
with n an integer number and H (1) n is the Hankel function of the first kind. The incident pressure field can also be rewritten as a sum of Bessel functions of the first kind J n :

p i (r, θ) = A ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ J 0 (k r) + 2 ∞ n=1 i n J n (k r) cos (n (θ -θ i )) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ∀r ≥ R, 0 ≤ θ ≤ 2 π ( 50
)
B n is a constant determined using Eqs. ( 49) and (50), and the boundary conditions applied on the grains defined in Eq. ( 3):

B n = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ -A J 0 (k R) H (1) 0 (k R) when n = 0 -2 A i k R J 1 (k R) -ρ 0 ρ J 1 (k R) k R H (1) 1 (k R) -ρ 0 ρ H (1) 1 (k R) when n = 1 -2 A i n J n (k R) H (1) n (k R) when n ≥ 2 (51)
Then the total pressure field is the sum of the incident field p i and the scattered field p r :

p(r, θ) = p i (r, θ) + p r (r, θ) (52)
Figure 4: Numerical solution of the total pressure field computed using the fictitious domain method (4a) versus the analytic solution (4b) for h Ω = 1/12 mm. Fig. 4a shows the total pressure field computed with the fictitious domain method and Fig. 4b shows the analytic solution computed using Eqs. ( 49) and (50). When comparing figures 4a and 4b, we observe that the numerical solution presents good agreements with the analytic solution.

Our model also provides the acoustic velocity field. The analytic solution of the total velocity field can be deduced from the solution of the pressure field Eq. (52):

u(r, θ) = -i ω c 0 ∇ p(r, θ) (53)
Figure 5 represents numerical and analytic solutions of the total velocity field. We can see that the two velocity fields are quite similar. The difference at the top of the domain is a numerical artifact due to the source line that generates a plane wave in two opposite directions. We define the relative error e Φ for any field Φ during a period T as:

e Φ = sup t ref ≤t≤t ref +T ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ Φ num -Φ an L 2 Φ an L 2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (54)
with t ref = H/c 0 the time for the wave to propagate through the entire medium. The relative error e Φ is computed inside a square of size [1.5; 8.5] × [1; 8] mm to avoid PML and the effect of lateral periodic boundary conditions. Relative errors computed for fields p, u x and u y are reported in Table 1. The growth of the error on the acoustic fields is not optimal, revealing the weakness of Lagrangian finite elements when used in mixed methods [START_REF] Li | Mixed methods using standard conforming finite elements[END_REF]. The agreement of the numerical solution with the analytic solution (Figs. 5 and4) shows nevertheless that the fictitious domain method converges with the element P 1 -P 0 on the contrary of the Raviart-Thomas element used in [START_REF] Imbert | Fictitious domain method to model a movable rigid body in a sound wave[END_REF].

U y [nm/s] N 2 =120 2 N 2 =240 2 N 2 =480 2 (a) y-component 0 1 2 3 4 5 t [μs] -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 U x [nm/s] N 2 =120 2 N 2 =240 2 N 2 =480 2 (b) x-component
Figure 6: Velocity of the grain function of the time computed using the fictitious domain method for several discretizations.

Velocity of the grain

The main feature of the model is to take into account the dynamic of the grain submitted to an incident wave. We plot on Figure 6 the velocity of the grain as a function of the time for three discretizations at t = 2.2 μs. We observe on Fig. 6a that grains start to oscillate at the same frequency as the incident wave, when the wavefront reaches its location.

The x-component of the velocity of the grain can also be assimilated to an absolute error because it should be U x = 0 according to the analytic solution. So we note on Fig. 6b that the decrease of maximal error is better for the velocity of the grain than for the norm of the velocity field. Indeed, the scale of the grain is larger than the one of the finite elements and the grain is therefore less sensitive to small perturbations of the pressure field.

Two movable grains fixed by springs

We perform a second experiment to show that the model can describe acoustic emissions directly coming from one grain subject to an external force F k . Most numerical methods in acoustics can not work without an incident wave. In this experiment, we link each grain with its own spring of stiffness k G = 100 GN/m (to keep approximately the same frequency as in the previous example f = √ (k/m)/(2 π)), so that the external force F k from Eq. ( 2) p. 3 becomes:

F k = -k G Y k -Y k eq y (55)
with y the unitary vector of the y-axis and Y k eq the equilibrium position of the spring of the grain k. The bottom grain has an initial displacement from equilibrium ΔY 0 (t = 0) = +1 pm and the top grain is at its equilibrium position (see scheme Fig. 7a). Both are put in the middle of a domain of size 15 × 17 mm and at a distance of 1 mm of each other. We set PML on each side of the box but we stop the simulation before the wavefront reaches the layer. Figure 7b shows the numerical pressure field for a discretization step h Ω = 1/12 mm at the final time t = 4.3 μs.

We examine now to the energy inside the system which is composed, according to system [START_REF] Bécache | Convergence results of the fictitious domain method for a mixed formulation of the wave equation with a Neumann boundary condition[END_REF], on one hand of acoustic energy in Ω :

E a = Ω 1 2 ρ 0 c 2 0 p 2 dx + Ω 1 2 ρ 0 u 2 dx (56)
and on the other hand of the mechanical energy of each grain k:

E m k = 1 2 k G Y k -Y k eq 2 + 1 2 (ρ -ρ 0 ) S k U 2 k (57)
Energies of the wave and the two grains are plotted in Figure 8. First we observe that the total energy in the system is well conserved considering the first order numerical scheme used for the time discretization. We observe that at t = 0, all the energy is in the grain 0 (in blue Fig. 8) in the form of mechanical potential energy due to the initial stretching of the spring. Then the grain 0 starts to oscillate and generates an acoustic field so its mechanical energy is partially converted to acoustic energy (in red Fig. 8). At t = 1 μs, the wavefront of the wave coming from the grain 0 reaches the other grain 1 and a part of the acoustic energy is converted to mechanical energy of grain 1 (in green Fig. 8). Then because the grain 1 is starting to oscillate, it is also radiating and gives back some mechanical energy in the form of acoustic energy.

As a comparison, we do the same experiment replacing water by air. Results reported in table 2 confirms that a numerical model able to link the wave equation and the equation of motion is essential to take in account transfers of energy through both water and grains. On the contrary, it is insignificant in the air.

Acoustic pulse propagation through a granular suspension

The last experiment involves a random suspension of 400 movable grains subject to a unique pulse. We show that the model is able to work with many grains, contrary to the one developed in [START_REF] Imbert | Fictitious domain method to model a movable rigid body in a sound wave[END_REF]. The source is on a horizontal line located at y src = 4.4 cm and the pulse is a the 4th derivative of a Gaussian with the central frequency f c = 150 kHz. We choose a pulse to follow easily the propagation of the wave through the medium. We consider a case where the characteristic wavelength λ/d = 10 because displacement of grains is maximum at λ/d 1. The domain Ω is a rectangle 2 cm width and 4.8 cm height discretized by N x × N y = 108 × 260 and so h Ω = 5/27 mm. The boundary conditions are the same as in section 4.1: PML on the top and on the bottom; and periodic boundary conditions on the left side and right sides. Fig. 9 show snapshots of the pressure field (a), the velocity field (b) and positions, velocities and acceleration of each grain (c). The thin figure between the velocity field snapshot and the pressure field snapshot shows the mean value of the y-component of the velocity field and helps the reader to follow the propagation of the pulse. The variations of pressure (Fig. 9a) around the grains determine the value of the hydrodynamic force acting on them. Since the hydrodynamic force is the unique force, its value is proportional to the acceleration A (in gray scale Fig. 9c).

Moreover, by observing red arrows in Fig. 9c, we can see that the wave sets in motion grains under its acoustic field. The result of the ratio between the amplitude of the velocity of one grain and the amplitude of the velocity field is provided by the following 2D analytic solution for low wavelength (inspired from the 3D version in [START_REF] Lamb | Hydrodynamics, sizième Edition[END_REF]) :

U u = 2 ρ 0 ρ + ρ 0 ≈ 0.57 (58)
Comparison of Figs. 9b and9c shows that numerical results are similar to the one predicted by the analytic solution.

Then, we can see Fig. 10 that the pulse propagates through the granular suspension from the top to the bottom. The mean value of the amplitude the velocity field Fig. 10b shows that the amplitude of the wave is decreasing as the energy of the wave is transmitted to the grains Fig. 10c. We also observe on Fig. 10a several other reflections after the wave pass. Thus, we can use solution (58) only at the beginning of the propagation when grains only subject to the incident wave not yet reflected waves from their neighbors.

Conclusions and discussion

Our results show that this fictitious domain method can couple the acoustic wave equation and the equation of motion for rigid grains. We show this using three numerical experiments. First, both acoustic fields and dynamics of grains computed using this model satisfy the analytical solution of a plane wave diffracted by a movable grain. In the second case, energy is transmitted between two grains and the fluid without loss of total energy. The last experiment shows that this method can also be expanded to hundreds of grains and used in the framework of acoustic experiments on granular suspensions.

The strength of this model is the ability of each grain to generate its own acoustic field when subject to an external force or due to an incident wave. As far as we know, there are no other acoustic model capable of coupling dynamics of rigid grains and waves. The rigidity of grains implies an infinite velocity inside the grains but still a finite velocity through the contacts network of a granular packing. In our knowledge, our model is the only one able to record acoustic emissions related to any forces acting on the grains of a granular medium.

Nevertheless interesting comparative studies could be done by observing changes in the reflected wave signal while changing parameters related to contact forces or materials of the grains. In the particular framework of compact granular media, interaction forces are defined according to linear or elastic models to simulate repulsion and friction between grains. So waves generated by grains are then directly related to the parameters of these models. Therefore we plan to use our model to study acoustic waves emitted by a granular medium during a gravitational destabilization.

Figure 2 :

 2 Figure 2: The regular grid T Ω h (on the left) is the grid where the physical variables u and p are computed and the unstructured grid T G h (on the right) is the domain of the Lagrange multipliers.

Figure 3 :

 3 Figure 3: Scheme of the experiment with 1 grain

Figure 5 :

 5 Figure 5: Numerical solution of the norm of the total velocity field computed using the fictitious domain method (5a) versus the analytic solution (5b) for h Ω = 1/12 mm.

  Scheme of the experiment (b) Pressure field at t = 4.3 μs

Figure 7 :

 7 Figure 7: Numerical experiment of two movable grains fixed by springs. Fig. 7a shows a simplified representation of the experiment and Fig. 7b the pressure field at the final time t = 4.3 μs.

Figure 8 :

 8 Figure 8: Evolution of the distribution of energy in the system containing two glass grains fixed by springs and surrounded by water

  Velocity field (u y and u )

Figure 9 :

 9 Figure 9: Propagation of a pulse through a granular suspension at t = 20 μs. Fig. (a) shows the acoustic pressure field, Fig. (b) the mean of the y-component and the norm of the velocity field and Fig. (c) positions, velocities (represented by red arrows) and acceleration (in grayscale) of each grain.

Figure 10 :

 10 Figure 10: Propagation of a pulse through a granular suspension at t = {25; 30; 35} μs.

Table 1 :

 1 Relative errors computed for pressure and velocity fields for different discretizations /12 mm 1.44 × 10 -1 2.38 × 10 -1 1.47 × 10 -1 240 2 1/24 mm 9.19 × 10 -2 1.28 × 10 -1 9.01 × 10 -2 480 2 1/48 mm 8.00 × 10 -2 1.09 × 10 -1 7.93 × 10 -2

	N 2	h Ω	e p	e u x	e u y
	120 2 1

Table 2 :

 2 Comparison of average energy distribution in water and in the air during the last period

	of oscillation
		water	air
	c 0	1,500 m/s	340 m/s
	ρ 0 1,000 kg/m 3 1.2 g/m 3
	E a	95.9 %	0.6 %
	E m 0	1.7 %	99.4 %
	E m 1	2.4 %	≈ 0 %

(a) Numerical solution (b) Analytic solution

✩ This work was supported by the ANR project STABINGRAM No. 2010-BLAN-0927-01. 1 D. Imbert's work was supported by the grant No. 211-14-231/ARED from the Région Bretagne.


Ivan Romanytsia

Jean-Paul Viricelle
email: viricelle@emse.fr

Philippe Vernoux

Christophe Pijolat

Application of advanced morphology Au-X (X = YSZ, ZrO 2 ) composites as sensing electrode for solid state mixedpotential exhaust NOx sensor

Keywords: Normal, Interligne : simple (Par défaut) +Corps (Calibri), Anglais (États Unis), Vérifier l

Among various NOx sensors developments, mixed potential sensor based on Yttria Stabilized Zirconia (YSZ) with a simple planar architecture Au/YSZ/Pt is of practical interest. Au composites electrodes were investigated to improve sensing performances . Such potentiometric solidstate gas sensors -were fabricated by screen-printing and tested for NO 2 detection. Electrochemical impedance spectroscopy has shown that the addition of YSZ in the Au electrode decreased the polarization resistance in air. In addition, the (Au+10 wt% YSZ)/YSZ/Pt sensor has a shorter response time and higher sensitivity to NO 2 ( range 20-100 ppm) at 450-550 °C in comparison with a reference Au/YSZ/Pt sensor. On the other hand, the addition of non-conductive ZrO 2 doesn't significantly modify the electrochemical property but strongly damages the sensor responses toward NO 2 .

Introduction

Continuous stringent legislations in automotive exhausts oblige the vehicles manufacturers to optimize combustion operations and to monitor pollutants such as nitrogen oxides [START_REF]Regulation No 715/2007 of the European Parliament and of the Council of 20[END_REF][START_REF] Docquier | Combustion control and sensors: a review Progress in energy and[END_REF]. From 2014, Euro 6 standards require the utilization of a NOx (NO and NO 2 ) catalytic after-treatment device for Diesel engines. Selective catalytic reduction (SCR) of nitrogen oxides by ammonia or urea [START_REF] Docquier | Combustion control and sensors: a review Progress in energy and[END_REF][START_REF] Pârvulescu | Catalytic removal of NO[END_REF] is the most advanced technology used in heavy diesel engine trucks but also in Diesel passenger cars. Conversion of NOx depends on oxygen partial pressure as well as on NH 3 /NO and NO/NO 2 ratios, and could reach above 99% at 573 K [START_REF] Pârvulescu | Catalytic removal of NO[END_REF].In this context, fast and accurate on-board diagnostics of NOx concentration in the exhaust are strongly required.

Due to this strong demand, there are numerous studies dealing with the development of NOx sensors for on-board detection in light and heavy duty vehicles. In this paper, we will focus only on potentiometric gas sensor designed with two electrodes and working in the same atmosphere, according to the mixed potential principle. A recent development of such a sensor is reported in [START_REF] Wang | Sensing Exhaust NO2 Emissions Using the Mixed Potential Principle[END_REF] . Most of these sensors are composed of a solid electrolyte based on yttriastabilized zirconia (YSZ) [START_REF] Wang | Sensing Exhaust NO2 Emissions Using the Mixed Potential Principle[END_REF][START_REF] Gao | Improvement of the NOx selectivity for a planar YSZ sensor[END_REF][START_REF] Somov | Gas analysis with arrays of solid state electrochemical sensors: implications to monitor HCs and NOx on exhaust[END_REF][START_REF] Zhang | Electrochemical studies on cell M/YSZ/Pt (M=Pt, Pt-Ga2O3) in NO, O2, N2 mixtures[END_REF][START_REF] Yin | Highly sensitive mixed-potential-type NO2 sensor using porous double layer YSZ substrate[END_REF][START_REF] Zosel | Au-oxide composites as HC-sensitive electrode material for mixed potential gas sensors[END_REF][START_REF] Miura | Mixed potential type sensor using stabilized zirconia and ZnFeO4 sensing electrode for NOx detection at high temperature[END_REF], but other electrolyte may be used : β-alumina (2Na 2 O-11Al 2 O 3 ) [START_REF] Guillet | Development of a gas sensor by thick film technology for automotive applications: choice of materials-realization of a prototype[END_REF] or samaria-doped-ceria (SDC) [START_REF] Pijolat | Detection of CO in H2-rich gases with a samarium doped ceria (SDC) sensor for fuel cell applications[END_REF]. Two electrodes, a reference one and a sensing one are associated with the electrolyte. These two electrodes must present different electrocatalytic activity to detect NO x such as the couples Pt-Au [START_REF] Gao | Improvement of the NOx selectivity for a planar YSZ sensor[END_REF][START_REF] Guillet | Development of a gas sensor by thick film technology for automotive applications: choice of materials-realization of a prototype[END_REF], Pt-MnCr 2 O 3 [START_REF] Yin | Highly sensitive mixed-potential-type NO2 sensor using porous double layer YSZ substrate[END_REF], Pt -(Au-oxide) [START_REF] Zosel | Au-oxide composites as HC-sensitive electrode material for mixed potential gas sensors[END_REF][START_REF] Zosel | Mixed potential gas sensor with short response time[END_REF], or Pt with other oxide [START_REF] Wang | Sensing Exhaust NO2 Emissions Using the Mixed Potential Principle[END_REF]. The signal of a two electrode potentiometric sensor, in a complex gas environment, shows a non-Nernstian behavior due to the establishment of a mixed potential [START_REF] Wang | Sensing Exhaust NO2 Emissions Using the Mixed Potential Principle[END_REF][START_REF] Gao | Improvement of the NOx selectivity for a planar YSZ sensor[END_REF][START_REF] Somov | Gas analysis with arrays of solid state electrochemical sensors: implications to monitor HCs and NOx on exhaust[END_REF][START_REF] Yin | Highly sensitive mixed-potential-type NO2 sensor using porous double layer YSZ substrate[END_REF] on the sensing electrode. This latter is linked to the superposition of parallel electrode reactions occurring at the triple phase boundaries (TPB) [START_REF] Yin | Highly sensitive mixed-potential-type NO2 sensor using porous double layer YSZ substrate[END_REF][START_REF] Zosel | Au-oxide composites as HC-sensitive electrode material for mixed potential gas sensors[END_REF][START_REF] Miura | Mixed potential type sensor using stabilized zirconia and ZnFeO4 sensing electrode for NOx detection at high temperature[END_REF][START_REF] Yin | Highly sensitive mixed-potential-type NO2 sensor using porous double-layer YSZ substrate[END_REF], such as the oxygen electrode reaction and the electrochemical reduction (oxidation) of NO 2 (NO). The sensitivity of NO 2 is always higher than that of NO. In addition, by using the Au/Pt couple of electrodes interfaced on YSZ, NO 2 and NO exhibit opposite responses, positive for NO and negative for NO 2 [START_REF] Gao | Improvement of the NOx selectivity for a planar YSZ sensor[END_REF].

The morphology of the two electrodes plays a significant role in the number of TPB [START_REF] Yin | Highly sensitive mixed-potential-type NO2 sensor using porous double layer YSZ substrate[END_REF][START_REF] Yin | Highly sensitive mixed-potential-type NO2 sensor using porous double-layer YSZ substrate[END_REF], then on the kinetic of the electrode reactions and finally on sensing performances. For instance, a sensor with two Pt based electrodes showing asymmetric morphologies can be effective for NO 2 detection [START_REF] Yang | High temperature potentiometric NO2 sensor with asymmetric sensing and reference Pt electrodes[END_REF][START_REF] Guillet | The influence of the electrode size on the electrical response of a potentiometric gas sensor to the action of oxygen[END_REF]. Composite electrodes made of Pt or Au mixed with an oxide [START_REF] Zhang | Electrochemical studies on cell M/YSZ/Pt (M=Pt, Pt-Ga2O3) in NO, O2, N2 mixtures[END_REF][START_REF] Zosel | Mixed potential gas sensor with short response time[END_REF][START_REF] Chaoyang | Preparation of nanostructured Pt-YSZ composite and its application in oxygen potentiometric sensor[END_REF][START_REF] Song | Electrochemical decomposition of NO over composite electrodes on YSZ electrolyte[END_REF][START_REF] Striker | Effect of nanocomposite Au-YSZ electrodes on potentiometric sensor response to NOx and CO[END_REF] are also strongly investigated to improve the sensing performances. Recent studies have reported the sensing performances of Au+YSZ [START_REF] Striker | Effect of nanocomposite Au-YSZ electrodes on potentiometric sensor response to NOx and CO[END_REF] and Au-NiO [START_REF] Zosel | Mixed potential gas sensor with short response time[END_REF][START_REF] Plashnitsa | Zirconia-based electrochemical gas sensors using nano-structured sensing materials aiming at detection of automotive exhausts[END_REF] composite electrodes.

Cross sensitivity to exhaust gases as CO and hydrocarbons can be eliminated by using a catalytic filter [START_REF] Gao | Improvement of the NOx selectivity for a planar YSZ sensor[END_REF][START_REF] Yang | High temperature amperometric total NOx sensors with platinum-loaded zeolite Y electrodes[END_REF] deposited on the top of the electrodes. This catalytic layer can fully oxidize these molecules into CO 2 which is not detected by the sensor. It can also catalyze the NO oxidation into NO 2 to reach the thermodynamic equilibrium. Therefore, the NO/NO 2 ratio is only dependent on the temperature for a fixed partial pressure of oxygen. Our aim is to develop a potentiometric sensor highly sensitive to NO 2 equipped with a catalytic filter. Then, at a constant tem-perature and oxygen partial temperature, this system could provide concentrations of NO 2 , NO and overall NO x . This paper describes the sensing performances for NO 2 detection of YSZ-based solid state planar sensors using Pt as a reference electrode and Au-composites as sensing electrodes. The objective was to increase the porosity and then the number of TPBs of the sensing electrodes. These latter were composed of Au mixed either with YSZ or ZrO 2 . The addition of YSZ into the Au electrode could provide a mixed electronic and ionic conductivity into the sensing electrode in order to delocalize the TPBs in the overall volume, while the incorporation of ZrO 2 should only act on the porosity. Electrochemical performances for the oxygen electrode reaction were investigated by impedance spectroscopy and compared with those of pure Au sensing electrode. Finally, the sensing performances for NO 2 were evaluated between 450 and 550°C.

Experimental

Solid state planar sensors (Figure 1) were fabricated by screen-printing (Aurel C890) on alumina substrates (CoorsTek). A functional and porous layer 5x5 mm of solid electrolyte (YSZ: (ZrO 2 ) 0.92 (Y 2 O 3 ) 0.08 provided by Superconductive) was deposited by screen-printing (200 mesh) from a home-made ink containing 4 g of YSZ powder, 1.7 g of organic binder (ESL V400A) and 1 g of solvent (ESL 404) mixed 10 minutes in a three-roll mill (Exakt 80E). Then, the samples were dried for 15 min in air at 120°C and sintered at 1380 °C in air for 2 hours with a ramp of 5°C/min. Platinum and gold electrodes, of 1 mm width and 4 mm length, were screen-printed on the surface of the YSZ layer using a commercial paste of Pt (ESL 5545) and Au (ESL 8880-H).

The distance between the electrodes was 2 mm. The deposited electrodes were further dried for 15 minutes at 120°C and then sintered at 850°C for 2 hours in air with a ramp of 5°C/min. The gold and platinum electrodes were denoted as the sensing (SE) and the reference (RE) electrode, respectively. Gold based composite electrodes were prepared by mixing the Au commercial paste (ESL 8880-H) with 10 wt% YSZ (denoted as AuY) or with ZrO 2 (denoted as AuZ) homemade inks.

YSZ ink is the same as for electrolyte. ZrO 2 ink was composed of 4 g of ZrO 2 powder (Merck), 1.7 g of organic binder (ESL V400A) and 1 g of solvent ESL (59 wt % of ZrO 2 ). This ZrO 2 ink (10 wt%) was mixed in a three roll mill with the Au commercial paste. The depositions were further dried for 15 minutes at 120°C and then sintered at 850°C for 2 hours in air with a ramp of 5°C/min. Final amount of ceramic fraction in composite electrodes is 9.5 %wt.

The powders (YSZ and ZrO 2 ) were characterized by using laser granulometry (MASTER-SIZER 2000 MALVERN) and BET measurements with nitrogen (ASAP 2020). Specific area of YSZ powder after calcination at 850 °C in air was found to be 15 m²/g while that of ZrO 2 was 3.8 m²/g. Mean particle size of agglomerates was 1.6 µm and 0.288 µm, for YSZ and ZrO 2 , respectively.

Number and nature of impurities in the commercial ink of Pt (ESL 5545) and Au (ESL 8880-H) were controlled by ICP (Jobin-Yvon JY138 Ultrace). The presence of 1500 ppm of Cu and around 30 ppm of Al, Na, and Zn have been detected in the Au paste whereas significant loadings of Ti ( 29000 ppm), Zn (9000 ppm), Fe (2500-3500 ppm), Cr (1500 ppm) and Al (1500 ppm) have been analyzed in the Pt ink sample.

Scanning electron microscopy (JEOL JSM 6500F) was used to investigate the morphology of the electrodes and of the YSZ layer. Samples were placed in epoxy resin and polished to obtain cross section. Surface of samples was covered with a conductive carbon layer (Quorum

Q150R ES).

The control of the sensor heating was performed with a Pt-resistance element (ESL 5545) placed on opposite side of the alumina support (Figure 1A). Before testing, all sensors were aged at 800°C for 120 hours in air. Symmetric sensors were also made to investigate the electrochemical performances of the Au based electrodes (Figure 1B). Three series of symmetric sensors were prepared by screen-printing with the same procedure than used for the electrolyte layer and the electrodes: Au/Au as a reference sample, AuY/AuY and AuZ/AuZ.

Resistivity of the YSZ porous layers and electrochemical properties of the Au composite electrodes were measured by using electrochemical impedance spectroscopy (EIS, VersaSTAT 3 -Princeton Applied Research). Each symmetric sensor was placed in a tubular furnace (CAR-BOLITE) in ambient air. AC signal amplitudes of 150 mV were applied between the two electrodes in the frequency range 10 5 Hz -10 -3 Hz. EIS measurements were performed between 550°C and 700°C with a step of 25°C both during heating and cooling down. Impedance diagrams were recorded during a temperature plateau of 3 h. EIS measurements were not possible below 550°C, due to the extremely high resistivity (higher than 10 MΩ) of the YSZ porous layer.

Fitting of the Nyquist plots was realized in ZView interface.

The sensing performances for NO 2 of the solid state planar sensors (Figure 1A) were measured in a synthetic gas bench (Figure 2) with a "base gas" containing 12 vol% O 2 , 1.5 vol% H 2 O, balanced with N 2 . Amounts of 20, 50 and 100 ppm of NO 2 for 15 min were added in the base gas. The total flow rate was 15 L h -1 per cell, simultaneously four sensors were tested. For heating control and signal processing, an electronic card coupled with a LabView interface was used. Response and recovery times of the sensor responses were compared taking into account that the dead volume in the sensor cell is around 30 cm 3 . Therefore, considering the flow rate in the sensor cell (15 l h -1 ), around 35 seconds are necessary to deeply change the NO 2 concentration. Response (and recovery) times were measured and correspond to the time to reach (or to recover) 90% of the maximum sensor response.

Sensors have been tested between 450 and 550 °C. The open-circuit voltage (OCV) was measured as follow:

OCV = E RE -E SE (1)
where E RE was the potential of the reference electrode (Pt) and E SE was the potential of the sensing electrode (Au, AuY or AuZ). The responses of the sensors toward NO 2 , ∆V, were calculated according to the following relationship:

∆V=OCV NO2 -OCV (2)
when OCV NO2 was the stabilized sensor signal measured during NO 2 injection.

3. Results and discussion.

SEM characterizations

Figures 3 and4 display that the YSZ layer is extremely porous, probably due to the presence of evaporated binder (30 wt%) and solvent in the initial paste. However, the surface of this layer evidences percolation between micrometric YSZ grains ensuring bulk ionic conductivity.

Electrochemical characterizations

Electrochemical Impedance Spectroscopy was used to investigate the electrochemical properties of the symmetric sensors Au/Au, AuY/AuY and AuZ/AuZ. Experiments below 550°C

were not possible to carry out due to the extremely high resistivity of the sensors, of the order of 10 MΩ. This is due to the high porosity of the electrolyte layer. Indeed, literature data reports high increase of the YSZ electrolyte resistance with the porosity [START_REF] Kleitz | Microstructure blocking effects versus effective medium theories in YSZ[END_REF][START_REF] Bultel | Electrochemical and catalytic properties of porous Pt-YSZ composites[END_REF].

Figure 5 shows a typical Nyquist plot of the Au-Au sample at 551°C in air. Impedance spectra exhibit two well-defined semicircles that could be separated and resolved. The first semicircle, in the high frequency (HF) domain, represents the resistivity of the YSZ solid electrolyte layer (R SE ) which is mainly due to the blocking of the conduction by the pores [START_REF] Kleitz | Microstructure blocking effects versus effective medium theories in YSZ[END_REF].

The second contribution, in the low frequency (LF) range (from 10 2 to 10 -3 Hz), corresponds to the overall polarization resistance (Rpol) at the two symmetrical electrodes linked with the kinetic of the electrode reactions [START_REF] Zhang | Electrochemical studies on cell M/YSZ/Pt (M=Pt, Pt-Ga2O3) in NO, O2, N2 mixtures[END_REF][START_REF] Chaoyang | Preparation of nanostructured Pt-YSZ composite and its application in oxygen potentiometric sensor[END_REF][START_REF] Bultel | Electrochemical and catalytic properties of porous Pt-YSZ composites[END_REF][START_REF] Wang | A study of factors that influence zirconia/platinum interfacial impedance using equivalent circuit analysis[END_REF].

The SE semi-circle can be reasonably fitted with a parallel R SE C SE element [START_REF] Zhang | Electrochemical studies on cell M/YSZ/Pt (M=Pt, Pt-Ga2O3) in NO, O2, N2 mixtures[END_REF][START_REF] Wang | A study of factors that influence zirconia/platinum interfacial impedance using equivalent circuit analysis[END_REF] with R CE and C SE , the resistance and the capacitance of the YSZ layer (Figure 5). The resistance of the YSZ layer is extremely high, around 2*10 6 Ω at 550°C and is almost similar for the three sensors (Figure 6). This underlines the good repeatability of the preparation method. This value is higher than the ones reported in previous studies [START_REF] Kleitz | Microstructure blocking effects versus effective medium theories in YSZ[END_REF][START_REF] Bultel | Electrochemical and catalytic properties of porous Pt-YSZ composites[END_REF]. This could be related to the high porosity of the screen printed YSZ layer. The activation energy is similar for the three symmetrical sensors, at around 1.13 eV. This value is in good agreement with that reported for porous YSZ membranes impregnated with 1 wt% Pt [START_REF] Bultel | Electrochemical and catalytic properties of porous Pt-YSZ composites[END_REF]. On the other hand, Kleitz et al [START_REF] Kleitz | Microstructure blocking effects versus effective medium theories in YSZ[END_REF] have found similar activation energies, at around 1.04 eV for both bulk and pores conductivity.

Modeling of LF part of impedance spectra can be performed with a RC element [START_REF] Zhang | Electrochemical studies on cell M/YSZ/Pt (M=Pt, Pt-Ga2O3) in NO, O2, N2 mixtures[END_REF][START_REF] Song | Electrochemical decomposition of NO over composite electrodes on YSZ electrolyte[END_REF], but better fits were obtained by using a parallel combination of a resistance (Rpol) with a constant phase element (CPEpol) [START_REF] Wang | A study of factors that influence zirconia/platinum interfacial impedance using equivalent circuit analysis[END_REF][START_REF] Woo | Impedance characterization of a model Au/Yttria-Stabilized Zirconia/Au electrochemical cell in varying oxygen and NOx concentration[END_REF][START_REF] Woo | Effect of electrode composition and microstructure of impedancemetric nitric oxide sensors based on YSZ electrolyte[END_REF][START_REF] Woo | Effect of electrode material and design on sensitivity and selectivity for high temperature impedancemetric NOx sensors[END_REF]. Figure 7 depicts the Arrhenius diagram of the polarization reresistance in air between 550 and 700 °C. From a general point of view, polarization resistances of the two symmetrical electrodes are one order of magnitude lower than resistances of the YSZ layer. The two samples Au/Au and AuZ/AuZ show similar values of polarization resistances as well as of the apparent activation energy (1.54 ± 0.03 eV). Therefore, the enhancement of the electrode porosity with the addition of an insulating material such as zirconia seems to not significantly modify the oxygen electrode kinetic. Opposite results were reported in the literature data for Pt composites electrodes containing Ga 2 O 3 [START_REF] Zhang | Electrochemical studies on cell M/YSZ/Pt (M=Pt, Pt-Ga2O3) in NO, O2, N2 mixtures[END_REF] deposited on a dense YSZ sheet. The presence of the oxide phase was found to inhibit the cathodic electrode reaction.

The electrode composite of Au with YSZ shows the highest electrochemical performances (Figure 7) and the lowest apparent activation energy (1.43 ± 0.3 eV). For instance, at 550°C, the value of R POL decreases from 5 10 5 Ω for Au and AuZ down to 3 10 5 Ω for AuY. The apparent activation energy values at the Au/YSZ interface in oxygen, reported in the literature data, are in the order of 0.9 -1 eV [START_REF] Woo | Impedance characterization of a model Au/Yttria-Stabilized Zirconia/Au electrochemical cell in varying oxygen and NOx concentration[END_REF][START_REF] Pil Yoon | Characteristic of cathodic polarization at Pt/YSZ interface without the effect of electrode microstructure[END_REF]. For instance, Woo et al. [START_REF] Woo | Impedance characterization of a model Au/Yttria-Stabilized Zirconia/Au electrochemical cell in varying oxygen and NOx concentration[END_REF] have recently performed experiments on a symmetrical Au/YSZ/Au sensor using thin Au plates contacted on a porous spray-coated YSZ layer. In the temperature range 600-700°C, they found an activation energy value of 1.02 eV under oxygen.

Our experimental data evidence high Ea values. Blocking process of the charge transfer may occur at the Au/YSZ interface which could be linked with the extremely high YSZ layer porosity as well as with the presence of impurities (Cu ~1500 ppm, Al, Na, Zn ~30ppm) in the Au paste.

Sensor performances

OCV measurements in the base gas.

The evolution of OCV in the base gas (12% O 2 , 1.5% H 2 O) was monitored as a function of time (Figure 8) at 450°C. These responses correspond to potential differences between the Pt electrode and the Au-based one (Figure 1a). Initial signal abrupt variations are due to the temperature regulation. Regarding the Au sensor, OCV values significantly increase with time during the first two hours and then gradually raise but never reach a steady-state. After 18 hours on stream, the OCV value of the Pt/Au sensor is + 57 mV. The OCV, in oxygen atmosphere, between Au and Pt electrodes is generally attributed to a difference in the O 2 electrode kinetic rate linked with different oxygen activity [START_REF] Mori | Detection of offensive odorant in air with a planar-type potentiometric gas sensor based on YSZ with Au and Pt electrodes[END_REF]. A positive value indicates that the oxygen activity is higher on Pt than on Au, in good agreement with the morphology of the two electrodes (Figure 4). The two Au composites present similar behaviors with an initial increase followed with a gradual decay of the OCV with time. Responses of AuZ and AuY sensors reach a steady-state value of +10 mV and +2 mV after 6.5 and 13 hours on stream, respectively. These results underline that the electrode porosity increases the surface specific area and the number of TPB, active for the oxygen adsorption. Therefore, the electrochemical properties of AuY for the oxygen electrode reaction become closest to those of the Pt electrode, in good agreement with EIS data (Figure 7).

NO 2 sensing performances

NO 2 sensing performances such as the response time and the signal stability were monitored as a function of the morphology and the composition of the sensing electrodes. In oxygen-rich atmosphere, in presence of NO 2 , a mixed potential is established between the oxygen electrode reaction and the electrochemical reduction of NO 2 [START_REF] Pârvulescu | Catalytic removal of NO[END_REF][START_REF] Chevallier | Non-Nernstian planar sensors based on YSZ with a Nb2O5 electrode[END_REF].

As already reported in the literature [START_REF] Gao | Improvement of the NOx selectivity for a planar YSZ sensor[END_REF][START_REF] Striker | Effect of nanocomposite Au-YSZ electrodes on potentiometric sensor response to NOx and CO[END_REF], whatever the temperature and the nature of the electrodes, the responses to NO 2 and NO are opposite. As expected, the sensors responses to NO 2 are negative (Figure 9). A 450°C, the AuZ sensor presents a strange and non usable behavior, suggesting extremely slow and complex interactions of NO 2 with the porous electrode.

On the opposite, the AuY sensor shows stable and reversible responses. In addition, this electrode exhibits the highest sensing performances with sensitivity to NO 2 of around -52 mV/decade. The recovery time of pure Au electrode is extremely high, most probably because a slow diffusion of NO 2 into the dense Au layer.

By increasing temperature to 500°C and 550 °C, the sensitivity decreases for all sensors. For instance on AuY, the response drops from 42 mV/decade at 500°C down to 28 at 550°C. In parallel, response and recovery times are strongly decreased (Figure 10). Responses of both Au and AuZ electrodes become reversible. However, the response time of AuZ remains much longer. Whatever the temperature, AuY exhibits the highest sensitivity and the lowest response time. Sensitivity of AuY sensor is two times higher than that of Au at 500°C and 550°C.

Nevertheless, the AuY response tends to decrease with time during the NO 2 injections.

Conclusions

The objective of this study was to study the influence of gold composite electrodes on the performances of a mixed potential Au/YSZ/Pt sensor. The effect of YSZ and ZrO 2 additions to the gold electrode on the response toward NO 2 was investigated between 450°C and 550°C.

The initial gold electrode was quite dense. SEM characterizations evidenced that the porosity and then the number of triple phase boundaries were increased when ZrO 2 and YSZ were incorporated into the Au electrode. However, with ZrO 2 addition, polarization resistance was not modified compared to reference Au/YSZ/Pt sensor. In this case, the presence of porosity in the electrode damages the sensor performances both in terms of NO 2 sensitivity and of signal stability. On the contrary, as shown by electrochemical impedance measurements, the polarization resistance was significantly reduced in the case of Au-YSZ electrode. Such a result indicates that the electrochemical properties for the oxygen electrode reaction were only improved for the Au-YSZ composite electrode because the ionic conductivity of YSZ can delocalize the electrode reaction in the overall electrode volume. Furthermore, the addition of YSZ improved the sensitivity toward NO 2 as well as the response time of Au-based electrodes.

As a conclusion, porous Au-YSZ composite electrodes in a YSZ-based sensor have a beneficial effect compared to pure gold electrode.

Figure captions

Figure

  Figure4Ashows that the Au sensing electrode is dense and totally covers the surface of the electrolyte. The interactions between the gold film and the porous YSZ layer are poor. In this case, TPB electrode/electrolyte/gas is limited by the perimeter of the Au dense film. Addition to the gold paste of 10 wt% of ZrO 2 or YSZ ink, produces strong modifications of the electrodes morphology (Figures4B, 4C). The Au-Z composite electrode presents open pores in which the gas can diffuse from the electrode surface to the electrolyte layer. The diameter of these pores varies from 0.5 to 5 µm. Metallic phase contain a wide size distribution of Au grains from 1 to 20 µm. Agglomerates of ZrO 2 larger than 2 µm can be observed inside the composite electrode Au-Z. Au particles in the Au-Y composite layer present a similar size from 1 to 20 µm. YSZ particles,

For example, Figure 11

 11 allows to compare responses and recovery times of the three sensors during addition of 100 ppm NO 2 at 500°C. Response time for AuY is the shortest one, at approximately 50 s. The Au sensor shows a longer response time around 120 s. The AuZ sensor is extremely slow: even after 15 min, the signal is not stable. Whatever the electrode, several minutes are necessary to recover the initial response after the end of the NO 2 exposure at 500°C. The utilization of Au-Y composite can strongly reduce this recovery time which is 8 min instead of 45 min for both Au and AuZ. These long recovery times are probably linked with a strong adsorption of NO 2 onto the electrodes. At higher temperature (550°C) and for a lower concentration of NO 2 (20 ppm), the recovery times become reasonable (considering the dead volume of the cell) on AuY (130 s) and Au (70 s).

Figure 1 :Figure 2 :Figure 3 :Figure 4 :

 1234 Figure 1: A) Schematic drawing of solid state gas sensors and B) symmetric samples for EIS measurements Figure 2: Test bench1 -Humidifier, 2 -Mixing chamber, 3 -Control card/potentiostat Figure 3: Surface SEM images of the screen-printing YSZ layer

Figure 5 :Figure 6 :Figure 7 :

 567 Figure 5: Experimental and fitted Nyquist plot of Au-Au sample at 550°C in air (numbers with arrow indicate the log 10 of the measuring frequency) Figure 6: Arrhenius diagrams of the solid electrolyte resistance (R SE )

Figure 8 .Figure 9 :Figure 10 :Figure 11 :Figure 2 :

 8910112 Figure 8. Sensor Pt/YSZ/SE baseline at 450°C in 12% O2, 1.5% H2O

Figure 3 :Figure 4 :

 34 Figure 3: Surface SEM images of the screen-printing YSZ layer

Figure 5 :Figure 6 :

 56 Figure 5: Experimental and fitted Nyquist plot of Au-Au sample at 550°C in air (numbers with arrow indicate the log 10 of the measuring frequency)

Figure 7 :

 7 Figure 7: Arrhenius diagrams of the polarization resistance (R POL ) of the three symmetric sensors Au/Au, AuZ/AuZ and AuY/AuY.

Figure 8 .Figure 9 :Figure 10 :Figure 11 :

 891011 Figure 8. Sensor Pt/YSZ/SE baseline at 450°C in 12% O 2 , 1.5% H 2 O

Mis en forme : Police :(Par défaut) +Corps (Calibri)