Kevin Bioche

Amanda Pieyre

Guillaume Ribert
email: ribert@coria.fr

Franck Richecoeur

Luc

Luc Vervisch

The role of gravity in the asymmetry of flames in narrow combustion chambers

Keywords: Narrow channel, Flame symmetry, Experiments, Numerical simulation, Gravity effect, Flame wall interaction

come

Introduction

In small-scale combustion, a vast amount of studies consider the problem as symmetric, either with respect to the middle plane in a flat channel or to the centerline in a tube [START_REF] Raimondeau | Modeling of high-temperature microburners[END_REF][START_REF] Kim | A numerical study on propagation of premixed flames in small tubes[END_REF][START_REF] Li | A numerical study on premixed micro-combustion of CH 4-air mixture: effects of combustor size, geometry and boundary conditions on flame temperature[END_REF][START_REF] Xie | Effects of non-catalytic surface reactions on the CH 4-air premixed flame within micro-channels[END_REF][START_REF] Ju | Studies of the effects of radical quenching and flame stretch on mesoscale combustion[END_REF][START_REF] Short | Asymptotic and numerical study of variable-density premixed flame propagation in a narrow channel[END_REF][START_REF] Ronney | Analysis of non-adiabatic heat-recirculating combustors[END_REF][START_REF] Gauthier | An evaluation of numerical models for temperature-stabilized CH4/air flames in a small channel[END_REF][START_REF] Gauthier | Effect of external heat loss on the propagation and quenching of flames in small heat-recirculating tubes[END_REF][START_REF] Norton | A CFD study of propane/air microflame stability[END_REF][START_REF] Norton | Combustion characteristics and flame stability at the microscale: A CFD study of premixed methane/air mixtures[END_REF][START_REF] Kizaki | Effect of radical quenching on CH 4/air flames in a micro flow reactor with a controlled temperature profile[END_REF]. However, recent numerical [START_REF] Tsai | The asymmetric behavior of steady laminar flame propagation in ducts[END_REF][START_REF] Pizza | Threedimensional simulations of premixed hydrogen/air flames in microtubes[END_REF][START_REF] Kurdyumov | Lewis number effect on the propagation of premixed flames in narrow adiabatic channels: Symmetric and non-symmetric flames and their linear stability analysis[END_REF][START_REF] Jiménez | DNS study of the propagation and flashback conditions of lean hydrogen-air flames in narrow channels: symmetric and nonsymmetric solutions[END_REF] and experimental [START_REF] Zamashchikov | Some features of gas-flame propagation in narrow tubes[END_REF][START_REF] Ju | Effects of channel width and Lewis number on the multiple flame regimes and propagation limits in mesoscale[END_REF] studies reveal the apparition of non-symmetric flame shapes for various conditions. Depending on the configuration, the physical phenomena susceptible to break the symmetry are numerous and the objective of this work is to focus on the role of gravity.

Aside from gravity, at least three generic types of instabilities have been discussed in the literature, which could potentially promote a non-symmetrical laminar flame behaviour, namely Darrieus-Landau [START_REF] Darrieus | Propagation d'un front de flamme[END_REF][START_REF] Landau | On the theory of slow combustion[END_REF], viscosity-induced Saffman-Taylor [START_REF] Saffman | The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid[END_REF] and diffusive-thermal [START_REF] Turing | The chemical basis of morphogenesis[END_REF][START_REF] Sivashinsky | Diffusional-thermal theory of cellular flames[END_REF] instabilities. The Darrieus-Landau instability is an amplification of the flame surface by hydrodynamic disturbances coming from gas expansion [START_REF] Matalon | The Darrieus-Landau instability of premixed flames[END_REF]. It was introduced to explain the predominance of slanted against symmetrical tulip-shaped flames, for a tube diameter exceeding a hundred times the flame thickness δ F [START_REF] Tsai | The asymmetric behavior of steady laminar flame propagation in ducts[END_REF]. However, because this instability requires a minimal characteristic length to develop, such an asymmetry was not observed for lower diameters as those considered in this study (about 15δ F ). On another hand, the viscosity-induced Saffman-Taylor (S-T) instabilities [START_REF] Saffman | The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid[END_REF][START_REF] Kang | Effects of heat and momentum losses on the stability of premixed flames in a narrow channel[END_REF] are typical of combustion in narrow channels, as a result of viscous-induced pressure gradients, which are favoured by small channel dimensions [START_REF] Joulin | Influence of momentum and heat losses on the large-scale stability of quasi-2d premixed flames[END_REF]. In the development of the diffusive-thermal instabilities [START_REF] Turing | The chemical basis of morphogenesis[END_REF][START_REF] Sivashinsky | Diffusional-thermal theory of cellular flames[END_REF], the Lewis number (Le) was found as a key parameter. For Le < 1, the existence of non-symmetric flames is discussed in various studies [14-17, 27, 28]. Besides, in the context of a thermo-diffusive model, Le > 1 flames with very low mass flow rate are also found asymmetric [START_REF] Kurdyumov | Lewis number effect on the propagation of premixed flames in narrow adiabatic channels: Symmetric and non-symmetric flames and their linear stability analysis[END_REF] and vibratory instability of the planar flame in the tube can develop [START_REF] Clavin | One-dimensional vibratory instability of planar flames propagating in tubes[END_REF]. Unity Lewis flames have been found symmetric when approaching the micro-combustion scale [START_REF] Gauthier | Effect of external heat loss on the propagation and quenching of flames in small heat-recirculating tubes[END_REF] and a stability analysis including heat loss also reports that unity Lewis number symmetric flames are stable [START_REF] Kurdyumov | Propagation of symmetric and non-symmetric premixed flames in narrow channels: Influence of conductive heat-losses[END_REF]. All of these studies do not include the gravity in the analysis.

In small-scale combustion, gravity is often neglected because the Froude number (Fr), inversely proportional to the characteristic length of the system, is expected larger than unity. For instance for methane burning in air under a fuel lean condition,

Fr = S o L 2 g × i = 1.59 , (1)
where g = 9.81 m • s -2 denotes the acceleration of gravity, S o L = 28 cm • s -1 is the methane/air adiabatic flame speed at equivalence ratio 0.8 considered in this study [START_REF] Vagelopoulas | Further considerations on the determination of laminar flame speed with the counterflow twin-flame technique[END_REF] and i = 5 mm is the characteristic length of the tube diameter or channel width. Nonetheless, having S o L 2 > g × i indicates that gravitational effects do not fully drive the flame dynamics, but not necessarily that they can be completely neglected.

Experiments of flame propagation in horizontal tubes show the appearance of slanted flames with a preferential direction [START_REF] Zamashchikov | Some features of gas-flame propagation in narrow tubes[END_REF][START_REF] Ju | Effects of channel width and Lewis number on the multiple flame regimes and propagation limits in mesoscale[END_REF]. So far this observation was explained by a flame stabilisation mechanism influenced by the free convection in the air surrounding the tube. A decrease of the Nusselt number toward the top of the cylinder, as the outside boundary layer thickens, has been reported in thermal and flow plume developing around heated tubes [START_REF] Saitoh | Bench mark solutions to natural convection heat transfer problem around a horizontal circular cylinder[END_REF][START_REF] Acharya | Natural Convection Heat Transfer From a Short or Long, Solid or Hollow Horizontal Cylinder Suspended in Air or Placed on Ground[END_REF][START_REF] Wang | Transient laminar natural convection from horizontal cylinders[END_REF]. The plume insulates the cylinder from the surrounding air and the reaction zone in the upper part of the tube is then submitted to a lower heat transfer coefficient. Then gravity is driving the flame shape indirectly by affecting the flow surrounding the tube. A direct influence of gravity inside wide tubes was also discussed in the literature using a model equation for flame propagation [START_REF] Kazakov | Analytical study in the mechanism of flame movement in horizontal tubes[END_REF], with buoyant induced effects affecting the flame shape inside the tube.

The objective of this paper is to progress on the understanding of these 'internal' gravity effects occurring inside the narrow combustion chamber by combining experiments and numerical simulations. The experimental and numerical set-up are first described. Then, the overall flames properties are reported, flames which are asymmetrical in all the conducted experiments. Two cases are studied, a flame freely propagating above an iso-thermal cold wall and a flame fixed in the laboratory frame, i.e. stabilised by the incoming flow, the latter subjected to strong coupling with wall heat transfer. The experimentally observed flame shape is reproduced by the simulation only in the presence of gravity and the breaking of symmetry is explained by the change of behaviour of the baroclinic torque across the flame front after introducing gravity.

Experimental and numerical set-up

Experimental configuration

The experimental set-up shown in Fig. 1 consists in a i = 5 mm inner diameter and e w = 1 mm thick quartz tube placed in horizontal position. The methane/air fresh mixture at a temperature T o = 300 K flows from the left to the right at an equivalence ratio of 0.8, with a Reynolds number of about 50. The length of the tube is 200 mm. The mixture is first ignited at the free end of the tube x e = 150 mm, then the mass flow rate is decreased for the flame to enter the tube, which may be stabilised at a given position for a bulk velocity of the incoming flow balancing its burning velocity.

The quartz was chosen for its thermal and optical properties, allowing for a direct visualization of the flame and its propagation. The diameter was fixed right above the quenching diameter of the mixture, according to the analysis reported in [START_REF] Daou | Influence of conductive heat-losses on the propagation of premixed flames in channels[END_REF]. Two mass flow controllers (EL-FLOW Bronkhorst) with 0-700 Nml/min range for air and 0-70 Nml/min for methane are enslaved to command the velocity and equivalence ratio of the flow.

A mixture chamber prior to the quartz tube entrance (not shown in Fig. 1) secures the premixing of the gas. With a precision of 0.35 mL/min, the mass flow controllers allow to adjust the bulk velocity to stabilize the propagating flame. The position and the motion of the flame in the tube is tracked by a camera (Nikon 7000D -Lens micro Nikkor 60 mm f/2.8G) recording videos, which are then decomposed in series of individual pictures. To observe accurately the flame front, and especially the CH * chemiluminescence, a band-pass interference filter centered at 430.0 nm with a fwhm of 10 nm (±2 nm) and 40% minimum peak transmission is used.

Case Wall seen by the flame The flame is examined, (i) stabilised by an adjusted incoming flow matching the overall burning rate and (ii) during an upstream propagation after decreasing the mass flow rate. Because the heat diffusivity is much larger in the gas than in the solid, in case (ii) the flame propagates along an almost iso-thermal wall. In case (i), the stabilisation mechanism of the steady flame depends on complex heat exchanges with the wall, including heat transfer between the gas and the wall, conduction of heat in the solid wall and convective heat transfer between the outside wall and the environment [START_REF] Jiménez | DNS study of the propagation and flashback conditions of lean hydrogen-air flames in narrow channels: symmetric and nonsymmetric solutions[END_REF][START_REF] Bioche | Premixed flame-wall interaction in a narrow channel: Impact of wall thermal conductivity and heat losses[END_REF]. The slowest thermal exchanges occurring in this configuration is heat conduction through the quartz wall, featuring a characteristic time t cd = e 2 w /α s , where α s is the quartz thermal diffusivity [START_REF] Momentive | Thermal Properties of Fused Quartz URL[END_REF]. A flame residence time t res = δ F /V p is defined from

V bulk [cm • s -1 ] V p [cm • s -1 ] t cd /t res S L [cm • s -1 ] (i)
δ F = α g /S o L = 94.
2 µm, a characteristic flame thickness and -V p x, the absolute velocity of the flame front in the laboratory frame, where α g denotes the fresh gases diffusivity. For t cd /t res < 1, the budget of heat fluxes reaches a steady state with a non-isothermal wall. In practice this is observed experimentally when V p → 0, representative of a stabilised flame (case (i)). For the ratio t cd /t res >> 1, the wall temperature at the flame position stays of the order of 300 K and the flame propagates over an iso-thermal wall (case (ii)). These two operating points are summarised in Table 1. The burning velocity, S L = V bulk + V p , where V bulk is the bulk velocity of the incoming fresh gases, increases by 9.17% between the flame moving upstream and the stabilised one. This enhancement of the burning velocity results from heat retrocession by the wall to the flow slightly upstream of the stabilised flame front, a mechanism which has been revisited recently in [START_REF] Bioche | Premixed flame-wall interaction in a narrow channel: Impact of wall thermal conductivity and heat losses[END_REF]. The increase in amplitude of the flame speed depends on the amount of heat release and therefore on the equivalence ratio of the mixture.

Numerical set-up

To verify that the studied gravity mechanisms do not depend significantly on the very detail of the flow configuration, and also because experiments in the flame tube show no sign of azimuthal effects, the most simple case of a planar channel is numerically studied in two-dimension, with a channel height i = 5 mm and quartz walls of thickness e w = 1 mm, thereby with characteristic lengths similar to the experimental ones.

The fully compressible form of the unsteady conservation equations of mass, momentum and total sensible energy are integrated with a finite volume method on a Cartesian grid with the SiTCom-B flow solver [START_REF] Domingo | Large-Eddy Simulation of a lifted methane-air jet flame in a vitiated coflow[END_REF][START_REF] Domingo | DNS and approximate deconvolution as a tool to analyse onedimensional filtered flame sub-grid scale modeling[END_REF][START_REF] Bouheraoua | Large-eddy simulation of a supersonic lifted jet flame: Analysis of the turbulent flame base[END_REF][START_REF] Duboc | Evaluation of chemistry models on methane/air edge flame simulation[END_REF]. Fourth-order skew-symmetric-like scheme for the convective fluxes [START_REF] Ducros | Large-eddy simulation of the shock/turbulence interaction[END_REF] and fourth-order centered scheme for the viscous and diffusive fluxes, are used for spatial integration. Time advancement is performed with a Runge-Kutta scheme of order four.

Due to CFL stability restrictions and spatial resolution, the time step is limited to 16 ns approaching the steady state. One-dimensional NSCBC [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF] are employed for inlet and outlet boundary conditions. The molecular transport properties of the gaseous mixture are computed following the Curtiss and Hirschfelder approximations [START_REF] Curtiss | Transport properties of multicomponent gas mixtures[END_REF]. A two-way flow/solid coupling is organised to solve heat transfers between the gas and the wall. The coupling procedure proposed in [START_REF] Duchaine | Development and assessment of a coupled strategy for conjugate heat transfer with large eddy simulation: application to a cooled turbine blade[END_REF] is employed for parallel computation in the two solvers. The alternate direction implicit Douglas-Gunn method [START_REF] Douglas | On the Numerical Integration of ∂ˆ2u∂xˆ2+∂ˆ2u∂yˆ2=∂u∂t by Implicit Methods[END_REF] is adopted to solve for the temperature in the solid. The exterior wall surface exchanges energy with the surrounding air at 300 K, with a heat transfer convective coefficient 30 W • m -2 • K -1 and by radiation, through a gray body hypothesis with an emission coefficient decaying linearly from 0.95 at 290 K to 0.75 at 1800 K. The wall thermal conductivity and capacity are the ones of fused quartz tabulated versus temperature from [START_REF] Momentive | Thermal Properties of Fused Quartz URL[END_REF]. The density of the solid wall is fixed at 2200 kg • m -3 . Due to the large difference in the characteristic time scales of internal energy evolution in the solid and the flow, a de-synchronization method is employed when converging toward the steady state solutions [START_REF] Koren | Modélisation des transferts de chaleur couplés pour la simulation multi-physique des chambres de combustion[END_REF]. More details concerning the set of equations solved along with the boundary conditions of the present case studied may be found in [START_REF] Bioche | Premixed flame-wall interaction in a narrow channel: Impact of wall thermal conductivity and heat losses[END_REF].

The origin of the axial coordinate (x = 0) is set at the flame position taken as the peak heat release rate on the axis of symmetry. The mesh extends from -51 mm to the left in the fresh gases to 15 mm to the right in the burnt gases. This domain is long enough to capture the upstream heat diffusion through the wall while ensuring a zero velocity gradient in the streamwise direction at the inlet. The mesh is composed of regular squares of resolution δ x = 25 µm from the inlet down to x = 9 mm. It is then progressively coarsened in the longitudinal direction down to the outlet, with a geometric coefficient of 1.0025. The spanwise mesh resolution does not vary in this Cartesian grid and is fixed to δ y = δ x , resulting in a mesh composed of 512k cells.

A specific procedure is applied to rapidly determine the inlet bulk velocity, so that the incoming flow balances the two-dimensional flame burning velocity. As in previous works [START_REF] Bioche | Premixed flame-wall interaction in a narrow channel: Impact of wall thermal conductivity and heat losses[END_REF][START_REF] Ruetsch | Effects of heat release on triple flame[END_REF], this is achieved by measuring the progression velocity of the methane iso-surface relative to the flow, at the location of its peak burning rate, to adjust the incoming velocity accordingly.

A reduced chemical mechanism composed of 17 species and 53 reactions has been specifically developed based on the GRI-1.2 mechanism [START_REF] Frenklach | GRI-Mech-An optimized detailed chemical reaction mechanism for methane combustion[END_REF] (Table 1 of supplemental material) using the ORCh (Optimised and Reduced Chemistry) approach [START_REF] Jaouen | Automatic reduction and optimisation of chemistry for turbulent combustion modelling: Impact of the canonical problem[END_REF], targeting the flame speeds and species profiles of freely propagating premixed flames for various levels of heat loss up to quenching and auto-ignition in homogeneous reactors at various initial temperatures. The reference adiabatic flame speed with the detailed scheme for the equivalence ratio φ = 0.8 is S o L = 28.40 cm • s -1 and the reduced scheme leads to S o L = 28.56 cm • s -1 . Representative species and temperature profiles compared in the one-dimensional flames at equivalence ratio φ = 0.8, between the detailed and reduced schemes are given in the supplemental material. Wall boundary conditions similar to the two experimental cases of Table 1 are considered, i.e. (i) heat-conductive with all heat transfers active and (ii) iso-thermal (T = 300 K). In the twodimensional iso-thermal channel with gravity, the flame burning velocity is S L = 24.95 cm • s -1 , to become 27.09 cm • s -1 with wall heat-transfer. As it should, the absolute values of these flame burning velocities differ between the axi-symmetric tube in the experiment (Table 1) and the two-dimensional channel in the simulation. However, the relative increase in burning velocity between iso-thermal and non-isothermal wall cases are quite close, 9.17% in the experiment and 8.57 % in the simulation, which brings some confidence in the retained strategy.

Analysis of gravity effects

Case (i): Flame stabilised with heat-conductive wall

The flame stabilised by an incoming flow exactly balancing its burning rate is considered at first (case (i) of Table 1). The CH * chemiluminescence is collected in the experiment over an exposure time of 2 s and the mean flame emission is represented in Fig. 2 left. The reaction zone makes an angle with the vertical of 24 • , whereas the numerical simulation without gravity reports a fully symmetrical flame shape (Fig. 2 top-right). Adding the gravity force in the Navier-Stokes equation solved, the flame takes a tilt as in the experiment (Fig. 2 bottom-right), even though the simulation is planar and the experiments axi-symmetrical. The angle of the reaction zone with the vertical is 30 • in the simulation.

In the case with gravity, the top and bottom close-to-wall edge-flame shapes differ (Fig. 2). The topology and the relative progression velocity of these edge-flames benefit from the preheating of the gases upstream of the flame after diffusion of heat inside the wall [START_REF] Bioche | Premixed flame-wall interaction in a narrow channel: Impact of wall thermal conductivity and heat losses[END_REF]. Streamwise profiles of velocity and temperature taken at a distance of 0.7 mm from the top and bottom channel walls are now analysed. The streamwise component of the velocity is larger at the bottom due to the confinement of the flow by the concave flame shape (Fig. 3(a)). This higher velocity goes with smaller residence times and thus a less efficient preheating by the wall thermal boundary layer, leading to smaller temperature levels ahead of the bottom edge-flame (Fig. 3(b)). In both experiments and simulations, the heat release rate is found larger at the bottom than at the top wall edge-flame (Fig. 2). As pointed out in [START_REF] Choi | Response of flame speed to positively and negatively curved premixed flames[END_REF], the burning rate is intensified by the curvature of concave flame fronts. This enhancement also explains the reduced quenching distance to the wall at the bottom edge-flame.

Because an eventual modification by gravity of the free convection surrounding the channel is not included in the simulation, another mechanism driven by gravity is at play inside the channel.

To isolate this mechanism from heat transfer inside and outside the wall, the flame propagating over an iso-thermal wall is further examined (case (ii) of Table 1).

Case (ii): Flame propagating over a quasi-isothermal wall

In the experiment, starting from a stabilised flame (case (i)), the mass flow rate is lowered to reach the operating point of case (ii), in which the flame proceeds inside the tube over quasiisothermal walls. The diagnostics reported above are applied, the recording starts ten seconds before the mass flow rate modification and is pursued up to a steadily propagating flame. The time evolution of the flame angle with the vertical determined from the fitted ellipse is given in Fig. 5. Starting at a steady state with an angle of 24 • , the propagation state is reached at which the angle is 3 • , yielding a relative decrease in inclination of 87.5%. This decrease of angle to the vertical (reduction of inclination) is also observed in the simulation. Snapshots of the iso-thermal wall simulations with and without gravity are shown in Fig. 6. The angle taken by the flame with gravity and iso-thermal wall is 6 • , also much less pronounced than in case (i) ( 30 • ).

The relative decrease in inclination of flame simulated is of 80%. The breaking of the symmetry is thus reported in both experiments and simulations with iso-thermal and cold (T =300 K) wall. Figure 6: Simulation of iso-thermal wall (case (ii)) of Table 1), with and without gravity. Dashed-line: lines used for probing the vorticity budget.

Analysis of baroclinic-torque response to gravity in 2D narrow channel

Kazakov [START_REF] Kazakov | Analytical study in the mechanism of flame movement in horizontal tubes[END_REF] discussed in a detailed theoretical analysis the effect of gravity on confined flames. Body forces generate pressure and density gradients and among the numerous coupling between density and pressure gradient present in flames, the baroclinic torque is a well-established source of vorticity in curved reaction zones [START_REF] Clavin | Combustion Waves and Fronts in Flows[END_REF]

B = 1 ρ 2 ∇ρ × ∇P = 1 ρ 2 ∂ρ ∂x ∂P ∂y - ∂ρ ∂y ∂P dx z , (2)
where ρ is the density and P is the pressure. x is the streamwise coordinate (flow direction), y is the transverse coordinate (gravity acceleration is -gy) and z the coordinate normal to the channel plane. The balance equation for the vorticity ω = ∇ × u reads

∂ω ∂t = -(u • ∇) ω i + (ω • ∇) u ii -ω (∇ • u) iii + B iv + ∇ × 1 ρ ∇ • τ v , (3)
where u is the velocity vector and τ is the viscous tensor. In this two-dimensional case, there is no change in vorticity due to vortex stretching and the term (ii) in Eq. ( 3) is zero. (i) is the transport of vorticity by convection, (iii) is vorticity stretching by density change and (v) is the transport of vorticity by viscous effects.

In the isothermal case without gravity, approaching the edge-flame close to the wall, the pressure and density gradients in the y-direction feature opposite sign on both sides of the channel centreline, whereas the density and pressure gradient in the x-direction stay the same. The baroclinic torque given by relation (2) therefore changes its sign on both sides of the channel centreline. This is verified in Fig. 7 displaying the amplitude of the baroclinic torque along the dashed-lines seen in Fig. 6(a). These plots are versus a reaction progress variable defined from the CO 2 mass fraction normalised by its value in the fully burnt gases (Y b CO 2 = 0.122576), also collected along the dashed lines of Fig. 6(a). The baroclinic torque in Fig. 7 is positive on upper part of the channel (Fig. 7(a)) and negative on the bottom part (Fig. 7(b)), confirming the symmetrical effects pictured by the white arrows in Fig. 6(a). All the terms contributing to the vorticity budget in Eq. (3) are shown in Fig. 7, along with the total budget, which sums up to zero as expected for both gravity (lines) and no-gravity (symbols) steady cases. In the case without gravity, all terms change their sign between the upper and bottom parts of the narrow channel, leading to a fully symmetric stable reaction zone. Because the density difference between fresh and burnt gases and flow acceleration are very close with or without body forces, the amplitudes of the various terms of Eq. ( 3) are only slightly affected by gravity. However, accounting for the effect of gravity, the terms are then not fully symmetric with respect to the channel axis (compare symbols with lines in Fig. 7). Actually, to elucidate the mechanism that makes the flame to rotate when gravity is introduced, it is necessary to examine the transient when the body forced is added. The introduction of gravity leads to an additional contribution to the baroclinic torque (2), which may be approximated as

(i) (iii) (iv) (v)
1 ρ 2       ∂ρ ∂x ∂P ∂y g - ∂ρ ∂y g ∂P ∂x       . (4)
Isolating the effect of gravity, its premier impact is to promote flow stratification with negative density and pressure gradients in the vertical direction, (∂P/∂y) g < 0 and (∂ρ/∂y) g < 0. Considering an isentropic flow at rest subjected to gravity, the relative variation of pressure and density defines the speed of sound c 2 = (∂P/∂ρ) > 1. Therefore, the amplitudes of the density and pressure gradients in the vertical direction and due to gravity may be ranked as (∂P/∂y) g < (∂ρ/∂y) g < 0. Across the premixed reaction zone, (∂P/∂x) < 0, (∂ρ/∂x) < 0 and (∂T/∂x) > 0. The ranking in pressure and density gradient evolves across the flame front (Fig. 8). In the upstream part of the flame front (∂ρ/∂x) < (∂P/∂x) < 0, while further downstream (∂P/∂x) < (∂ρ/∂x) < 0. These gradients in the streamwise direction are orders of magnitudes larger than those in the vertical direction, therefore they may be assumed weakly affected by the addition of gravity. Combining these observations, the baroclinic torque induced by gravity (Eq. ( 4)) should be largely positive in the upstream part of the flame front and decrease after, yielding globally a positive enhancement of the baroclinic torque across the reaction zone. To verify this simple scaling, at the time t = 18 µs after adding gravity, the source of vorticity is analysed by computing the variation of all the terms of Eq. ( 3). This variation is measured in the simulation between their steady state value without gravity (Fig. 7) and their value at t = 18 µs after gravity addition. In both the upper and the bottom edge-flame close to wall, a positive source of baroclinic torque is indeed observed, corresponding to a positive source of vorticity (Fig. 9). This addition of vorticity represents a few percent of the overall baroclinic torque and is located across the reaction zones, where the longitudinal pressure and density gradients occur. In the top part of the channel, the streamlines are then less deviated toward the centreline with gravity when crossing the flame (Fig. 10), leading to their spreading upstream of the flame with a local flow deceleration, followed by an upstream flame movement (see Fig. 11, dashed line). In the bottom part of the channel, the opposite mechanism is found, with a highest concentration of the streamlines by the added vorticity (Fig. 10

Conclusion

The effect of introducing gravity in flames propagating in a narrow-channel is studied from numerical simulations with comparison against experiments conducted in a narrow-tube (channel -0.7 -0.5 -0. of internal width or tube diameter of i = 5 mm). Experimental measurements and simulation reports similar trends for flames freely propagating above iso-thermal walls or stabilised by the incoming flow with strongly coupled heat exchanged with and within the wall. Both flames are asymmetrical, with an inclination that is less pronounced in the iso-thermal wall case. In the absence of modification by gravity of the convection on the outside channel wall, the flame is still tilted, with an anchoring at the upper wall. The heat-retrocession via conduction in the wall and the thermal boundary layer in the fresh gases upstream of the flame, increase the inclination of the slanted flames. Following previous works [START_REF] Kazakov | Analytical study in the mechanism of flame movement in horizontal tubes[END_REF], the response of the baroclinic torque to gravity is explored and results confirm its driving role in such narrow combustion-systems.

Figure 1 :

 1 Figure 1: Schematic view of the experimental set-up composed of a quartz tube of i = 5 mm inner diameter and e w = 1 mm thickness, inside which a lean methane/air mixture is injected. The inlet (x i = -50 mm) is connected to the mass flow controllers while the outlet (x e = 150 mm) is in open-air, where the flame is first ignited to then propagate inside the tube.

Figure 2 :

 2 Figure 2: Case (ii) of Table 1. Experiment (flame in tube): Mean chemiluminescence CH * . Simulation (flame in channel): Wall temperature and heat release rate.

Figure 3 :

 3 Figure 3: Streamwise velocity and temperature distribution at a distance of 0.7 mm from the top and bottom walls. x flame : maximum of heat release on the probed line (top or bottom). No-gravity: symmetrical flame.

  Figure 4 shows flame images taken initially and then subsequently at 10 and 15 seconds. The flame inclination is significantly reduced.

Figure 4 :

 4 Figure 4: Experimental snapshots at t = 0 s, t = 10 s and t = 15 s. White dotted line: flame detected contour. Red dashed line: fitted ellipse.

Figure 5 :

 5 Figure 5: Time evolution of the flame angle to the vertical. t = 0 denotes the decrease in mass flow rate. Dotted: Stabilised flame with heat-conductive wall (case (i)). Triangle: Propagating flame with iso-thermal wall (case (ii)).

Figure 7 :

 7 Figure 7: Vorticity budget versus Y CO 2 /Y b CO 2 along the dashed-lines of 6. Lines: With gravity. Symbols: Without gravity. Circles: (i) of Eq. 3 . Triangles: (iii). Squares: (iv) Baroclinic torque. Crosses: (v). Stars: budget.

Figure 8 :

 8 Figure 8: Pressure and density streamwise gradients versus Y CO 2 /Y b CO 2 along the dashed-lines of figure 6(a).

Figure 9 :

 9 Figure 9: Deviation of vorticity balance versus Y CO 2 /Y b CO 2 along the dashed-lines of Fig. 6. Dashed-dot: (i) of Eq. 3. Dotted: (iii). Solid: (iv) Baroclinic torque. Dashed: (v).

  (b)), leading to a local flow acceleration and a downstream flame movement (see Fig.11, dotted line). The net result is a reaction zone that is pushed downstream at the bottom and pulled upstream at the top to evolve towards a slanted shape.

Figure 10 :

 10 Figure 10: Case (ii). Streamlines and flame contours of 10, 30, 50 and 70% of max heat release rate.

Figure 11 :

 11 Figure 11: Velocity distribution at a distance of 0.7 mm to the top and bottom wall. x flame : maximum of heat release on the probed line. No-gravity: symmetrical flame.

Table 1 :

 1 Experimental conditions. Lean premixed methane/air flame (φ = 0.8). V bulk x is the bulk flow velocity. -V

		Heat-conductive	23.10	≈ 0	≈ 0	23.10
	(ii)	Iso-thermal	20.70	0.46	53	21.16

p x is the flame front velocity in the laboratory frame (see Fig.

1

for axis). t cd is a characteristic heat conduction time in the solid. t res is a flame residence time. S L is the burning velocity.

Acknowledgments

This work was granted access to the HPC resources of IDRIS, CCRT and CINES under the allocation A0032B07358 made by GENCI (Grand Equipement National de Calcul Intensif). The computing resources of CRIANN were also used. This research is founded by the Agence Nationale de la Recherche in the framework of the project ANR-14-CE05-0030 -MAPEE: Microcombustion Assistée par Plasma et Excès d'Enthalpie.


L V B Diop

O Isnard

B Li

X H Luo

H Wang

W J Ren

S Yano

C W Wang

J S Gardner

K D Liss

P Miao

S H Lee

T Kamiyama

R Q Wu

Y Kawakita

Phys. Rev. B Z D Zhang

Magnetically driven giant negative thermal expansion covering room temperature in Hf0.875Ta0.125Fe2

Keywords: A. Intermetallic compound, C. Neutron powder diffraction, D. Negative thermal expansion, D. Magnetoelastic transition

Neutron diffraction is used to reveal the origin of the unusual thermal expansion properties of the itinerant-electron system Hf0.875Ta0.125Fe2. Hf0.875Ta0.125Fe2 shows temperature-induced magnetic transitions from ferromagnetic (FM) to antiferromagnetic (AFM) order and then to the paramagnetic (PM) state upon heating. The FM-AFM transition proceeds in a stepwise fashion, as a first-order phase transition, and is accompanied by an extremely anisotropic lattice collapse. The unit cell volume shrinks abruptly, but only in the basal plane; the dimension along the c-axis varies almost continuously. Hf0.875Ta0.125Fe2 exhibits a 0.37 % spontaneous volume contraction across the FM-AFM transition, where a giant negative thermal expansion (NTE), with a volumetric thermal expansion coefficient -74×10 -6 K -1 , is observed over a temperature interval of ΔT~50 K around room temperature.

Introduction

It is general knowledge that most materials expand when they are heated and contract upon cooling under constant pressure (positive thermal expansion, PTE). However, only a few materials behave inversely, i.e., shrink upon heating while expand upon cooling, which is called negative thermal expansion (NTE). Materials featuring the NTE properties are of great importance for both scientific research and modern technological applications [1][2][3][4] . NTE systems have enormous industrial merit because they can be used to tune and control the overall thermal expansion of materials. NTE materials are generally combined with usual PTE materials to fabricate the zero thermal expansion (ZTE) composites [5][6][7][8][9][10][11] , which are desirable in many fields such as high-precision optics and aerospace [12][13][14][15] . NTE materials are mainly used as the thermalexpansion inhibitor for making the composites with precisely tailored coefficient of thermal expansion. Therefore, there are a number of important potential technical applications of NTE materials used as, for example, high-precision optical mirrors, optical fiber reflective grating devices, athermalizer for Bragg grating, printed circuit boards and machinery parts.

Up to the present time materials of various kinds have been identified and investigated for their potential as NTE materials. In addition to the well-studied ZrW2O8 family of materials 4 , NTE effects have also been observed in silicates β-LiAlSiO4 [START_REF] Sleight | [END_REF] , cyanides Cd(CN)2 17 , ReO3 18 , CuO nanoparticles 3 , ScF3 [19][20][21] , PbTiO3-based ferroelectrics 22 , antiperovskite manganese nitrides [23][24] , and as well as reduced layered ruthenate Ca2RuO3.74 25 . Nevertheless, only a very limited number of NTE materials serve as high-performance thermal-expansion compensators, due to low coefficient of thermal expansion, the relatively narrow NTE temperature window, thermal expansion anisotropy, as well as low mechanical and/or electrical insulating properties. Traditional NTE systems with flexible framework in the crystal structure (e.g., ZrW2O8 and ScF3) usually display NTE over a wide temperature range up to 1000 K, but the related coefficient of linear thermal expansion, , can hardly go beyond -10×10 -6 K -1 .

Unfortunately, their NTE coefficient, is so small that a large amount of addition is needed in order to obtain high-performance ZTE composites. In practical applications, a large coefficient of thermal expansion, which can effectively compensate the thermal expansion of PTE matrix but with little influence on its original physical properties, is thus desirable for NTE material.

Colossal NTE has been discovered during the last decade in various classes of materials, particularly systems associated with a metal-insulator 25,26 , intermetallic charge-transfer 27,28 or ferroelectric transition 29 . Most recently, the utilization of a magnetic phase transition accompanied by large volume reduction upon heating is deemed as a promising avenue towards discovery of giant NTE. The effectiveness of this approach was recognized widely by the giant NTE reported for the Mn3AN-based antiperovskite manganese nitrides (A: the main group elements or 3d elements) which use volume change due to the magnetic transformation, i.e., magnetovolume effects. The giant NTE of the antiperovskite manganese nitrides strongly influenced subsequent research, leading to the discovery of many phase-transition-type NTE intermetallic systems such as MnCo0.98Cr0.02Ge 30 and La(Fe,Si,Co)13 31 , which exhibit large magnetovolume effects. For instance, a large coefficient of linear thermal expansion  = -26×10 -6 K -1 was achieved in the itinerant-electron compound LaFe10.5CoSi1.5 by proper chemical doping in La(Fe,Si)13-based alloys. In the magnetic systems, the magnetovolume effect may significantly influence the overall thermal expansion. A prominent example is the Invar alloys, which present very low PTE or NTE below the magnetic ordering Néel or Curie temperatures. The magnetic moment-volume instability was argued to cause NTE, which suppresses the usual PTE and consequently leads to the Invar effect. It is well established that the unit-cell volume of a ferromagnetic (FM) spin configuration is usually larger than that of an antiferromagnetic (AFM) one; hence giant NTE could be discovered at the FM-AFM magnetic transition.

To this end, one interesting example is the itinerant-electron system Hf1-xTaxFe2, when 0.1 < x < 0.3. The alloys in this compositional range crystallize in the hexagonal MgZn2-type structure (C-14 Laves phase, space group P63/mmc) in which the Fe atoms are located on two inequivalent Wyckoff positions (2a and 6h) with the Hf/Ta atoms occupying the 4f site. The compounds undergo a first-order ferromagnetic (FM) to antiferromagnetic (AFM) phase transition upon heating [32][33][34][35] . This FM-AFM first-order phase transition is critical for technological applications, as it is accompanied by strong responses to relatively minor stimuli, such as large magnetocaloric, magnetoresistive, magnetostrictive or external pressure effects.

In the present work, we report giant room-temperature NTE in polycrystalline Hf0.875Ta0.125Fe2, chosen in the most interesting part of the (Hf,Ta)Fe2 magnetic phase diagram with coexistence of FM and AFM behavior. We investigate the magnetovolume effects, such as the spontaneous volume magnetostriction and the pressure effects, in correlation with magnetic properties due to the Fe sublattice. We report how first-order phase transition correlates to instabilities of magnetic structures in this itinerant-electron compound. To determine the spin configurations and understand their coupling with the lattice, we conducted a detailed high-intensity neutron powder diffraction study as a function of both temperature and applied pressure.

Methods

The alloy of composition Hf0.875Ta0.125Fe2 was prepared by melting high-purity constituting metals (at least 99.95%) in a high frequency induction furnace under a purified atmosphere of argon. The sample was melted three times with the button flipped over after each melting, which is usually sufficient to achieve compositional homogeneity. The so-obtained ingot was wrapped in tantalum foil, sealed in an evacuated quartz tube and annealed in a resistive furnace as follows. It was heated up to 1000 °C at a rate of 5 °C/minute and kept at this temperature for one week. Then, the power was switched off and the furnace quickly cooled down to room temperature. Neutron powder diffraction (NPD) measurements were carried out at the Institut Laue Langevin (ILL) in Grenoble, France. About 3 g of powder were introduced into a cylindrical vanadium sample holder (D = 6 mm, H = 5 cm) and mounted on the stick of a He cryostat, whose contribution to the powder diffraction patterns was eliminated using a radial oscillating collimator. Several patterns were collected at selected temperatures ranging between 2 and 400 K on the high-intensity diffractometer D1B (λ = 2.52 Å and λ/2 = 1.26 Å) with a detector angular range coverage 5° ≤ 2θ ≤ 128° which is especially suited for magnetic structure determination. All measurements were performed upon heating after temperature stabilization with typical acquisition times of 20 minutes per isotherm. A second set of powder diffraction patterns was collected while ramping temperature at 1 K/minute. The Rietveld analysis was performed using the FullProf program 36 , which allows the simultaneous refinement of structural and magnetic profiles. NPD spectra were also recorded under hydrostatic pressure up to 1 GPa on D1B instrument. The pressure cell was made of Ti-Zr alloys. In order to avoid diffraction from the cell itself a composition of Ti0.525Zr0.475 mass content was used enabling the compensation of the coherent scattering length bc of Ti and Zr (bc Ti = -3.438 fm, bc Zr = 7.16 fm). The aluminium container was filled with fluorinert F77 as a pressure transmitting medium.

Results and discussion

For a deeper insight into the peculiar behavior of this system, a precise knowledge of both crystallographic and magnetic structure is required. With that aim, we have carried out NPD studies to establish the magnetic structures of this compound and its temperature dependence. A detailed Rietveld analysis has been performed on diffraction patterns collected at selected representative temperatures in the PM, AFM and FM states, in order to refine the crystal structure and magnetic ordering -Figs. 1a-1c. The pattern recorded at 375 K (PM state) is characteristic of the nuclear Bragg scattering only and all the Bragg peaks could be identified to arise from the original C-14 Laves phase structure (space group P63/mmc) found at room temperature.

The analysis of the 2 K pattern reveals that the MgZn2-type crystal structure is preserved down to 2 K; no additional reflections have been observed as compared to the pattern collected in the PM phase. Hf0.875Ta0.125Fe2 shows a ferromagnetic ordering of Fe moments below 310 K. The contribution of the magnetic order to the diffraction pattern is most clearly reflected by the presence of magnetic scattering on several Bragg peaks such as (100), (101) and The Fe magnetic moments have been refined independently for the two non-equivalent sites and lie within the ab-plane of the hexagonal lattice. At 2 K, the refined moments of Fe at the 2a and 6h positions amount to 1.23 and 1.36 B respectively. They are converted to 2.66 B for one chemical formula, a value very consistent with the spontaneous magnetization 2.82 B at 4 K for Hf0.875Ta0.125Fe2 33 . Hf0.875Ta0.125Fe2 is antiferromagnetic only in a limited temperature interval between TFM-AFM = 310 K and TN = 334 K. Over the same temperature range, the patterns display the same Bragg reflections as that in the PM regime, no additional magnetic peaks are detected. The analysis of the position of the magnetic reflections reveals that the corresponding AFM structure is described by the magnetic propagation vector k = (0, 0, 0) meaning that the magnetic unit cell is the same as the crystallographic one. The best refinement at 320 K is achieved for a spin configuration in which Fe-6h magnetic moments align ferromagnetically within the same layer while the coupling between successive Fe-6h layers is antiferromagnetic. In this AFM structure, the determined moment for Fe-6h atoms at 320 K is around 0.62 , the Fe-2a atoms carry no ordered magnetic moment. The magnetic structures of Hf0.

K, FM

The thermal dependencies of the lattice parameters a and c, and the unit-cell volume V are displayed in Fig. 2 and Fig. 3, respectively. Below TFM-AFM, a cell parameter stays nearly constant, while c gradually decreases upon cooling. At TFM-AFM, the a lattice constant exhibits a large drop a/a ≈ 0.27 % and the system gets transformed into the AFM state. On the contrary, c shows only a tiny change at TFM-AFM. The first-order nature of the magnetic transition is reflected in the steep decrease of a parameter rather than the small anomaly in c lattice constant.

It can be clearly seen that there exists a temperature range in which the cell dimensions in the basal plane expand with decreasing temperature. This signals a strong negative thermal expansion, leading to the linear thermal expansion coefficient  = (1/a)(Δa/ΔT) = -54×10 -6 K -1 over a temperature window of 50 K. This value is one order of magnitude larger than that for commercial NTE materials. [Typically,  is of the order of 10 -6 K -1 , although for ZrW2O8,  = -9×10 -6 K -1 (Ref. 16).] The first-order FM-AFM magnetic phase transition is accompanied with a significant spontaneous volume magnetostriction V/V ≈ 0.37 %. This large magnetovolume effect is due to the abrupt change in the local magnetic moment, i.e., the amplitude difference of the Fe magnetic moment between the FM and AFM state. K -1 . In the here studied compound, the volumetric thermal expansion coefficient is at least an order of magnitude larger than the value found in isostructural system (Hf,Nb)Fe2 Because of the observed intimate relation between the crystal structure and the magnetic behavior in this itinerant-electron system, the application of hydrostatic pressure, thus changing the interatomic distances, should modify the transition temperatures, stabilizing the low-volume AFM phase. To probe the pressure dependence of the hexagonal lattice, NPD spectra have been collected up to 1 GPa at room temperature. NPD patterns recorded under pressures reveal that the hexagonal symmetry of the atomic arrangement remains unaltered within the investigated compression range; no indication of a structural phase transition was found up to the highest applied pressure. The lattice constants decrease linearly with pressure (Fig. 4), and values of da/dP = -10.1×10 -3 Å GPa -1 and dc/dP = -17.0×10 -3 Å GPa -1 are obtained. The contraction of the hexagonal lattice is anisotropic and more pronounced along the c direction whereas this lattice parameter was much less sensitive to the magnetic ordering temperature. The pressure dependence of the unit-cell volume is illustrated in Fig. 5. The pressure induced decrease of the volume is linear with a slope of dV/dP = -10.6×10 -1 Å 3 GPa - 1 , leading to a compressibility value of  = 0.62×10 -2 GPa -1 . The obtained value compares very well with the compressibility for neighbouring Hf0.825Ta0.175Fe2 composition ( = 0.57×10 -2

GPa -1 Ref. 37).

Conclusion

In summary, the results presented here elucidate how NTE for the itinerant-electron system Hf0.875Ta0.125Fe2 is intimately correlated with the magnetic structure. NPD studies This FM-AFM transition is of first order and is accompanied by a stepwise drop of the magnetic moment and by an isomorphic lattice collapse of a peculiar anisotropic kind. The results unambiguously demonstrate that the disappearance of the magnetic moment of the Fe atoms at the 2a site plays a crucial role in the NTE of this system and sheds light on searching for NTE materials by tuning the competing magnetic orders.

FIG. 1 .

 1 FIG.1. Refinements of the neutron diffraction patterns recorded for Hf0.875Ta0.125Fe2 at 375 K (a), 320 K (b) and 2 K (c). In the 375 K pattern, the top row of Bragg markers is for the Hf0.875Ta0.125Fe2 nuclear contributions, the second row corresponds to the position of the vanadium sample environment contribution. At 2 and 320 K the top and second rows of Bragg markers are referring to the nuclear and magnetic contributions of the Hf0.875Ta0.125Fe2 phase. The third row is corresponding to the nuclear contribution from the vanadium sample environment.

  One can observe a peculiarly anisotropic character of the volume collapse at the magnetic transformation: the lattice shrinks suddenly, but only in the basal plane. The dimension along c-axis varies almost continuously across the transition. The 0.37% volume drop occurs over a temperature interval of ΔT~50 K; the unit-cell volume shows a prominent NTE effect over this temperature span. It clearly indicates that the lattice expansion induced by the magnetic transition is larger than the thermal contraction due to the temperature decrease. Based on the refined values, we calculated the volume (bulk) thermal expansion coefficient  = (1/V)(ΔV/ΔT) and the corresponding value is  = -74×10 -6

  phenomena responsible for the giant NTE observed over a 50 K temperature intervalvolumetric thermal expansion coefficient -74×10 -6 K -1 . It is shown that this unusually large NTE is induced by the magnetoelastic transition from FM to AFM states.

Acknowledgements

The authors gratefully acknowledge the Institut Laue Langevin for access to the neutron facility.


Matthieu Virot

Laurent Venault

Philippe Moisy

Sergey I Nikitenko

Sonochemical Behaviour Of Pu

Sonochemical Behaviour of Pu(III), Pu(IV), and Pu(VI) in Aqueous Nitric Solutions

Sonochemical Behaviour of Pu(III), Pu(IV), and Pu(VI) in Aqueous Nitric Solutions

Matthieu Virot*1) , Laurent Venault 2) , Philippe Moisy 2) and Sergueï I. Nikintenko 1)

1) ICSM-UMR5257 CNRS/CEA/UM2/ENCSM, Site de Marcoule, 30207 Bagnols/Cèze, France 2) CEA/DEN/DRCP, Radiochemistry and Process Department, 30207 Bagnols/Cèze, France

Reprocessing of spent nuclear fuel is currently carried out via the well-established PUREX process which is a hydrometallurgical treatment performed in nitric acid and based upon redox operations followed by selective solvent extraction of U and Pu. 1 In order to adapt this reprocessing technology to the future nuclear reactors but also to reply to the society concerns, new alternatives are nowadays under investigation to minimise the quantity of generated radioactive wastes, replace obsolete chemicals, increase safety and propose new routes of treatments also supposed resistant against proliferation. Sonochemistry uses the propagation of ultrasonic waves in liquid media to generate physico-chemical transformations in homogeneous solutions and heterogenous systems. 2 The effects generally attributed to ultrasound do not result from direct interactions of acoustic waves with molecules but rather from acoustic cavitation, which is the nucleation, growth and rapid implosive collapse of vapour filled micro-bubbles. At collapse, these bubbles lead to extreme transient conditions of several thousands of degrees and hundreds of atmospheres while keeping the solution at ambient temperature and pressure.

Sonication may therefore induce chemical transformations of the liquid medium involving redox reactions in solution. 2 For instance, sonolyse of pure water allows the accumulation of H 2 O 2 in solution while the sonication of aqueous nitric solutions produces NO x and HNO 2 . In the presence of nitrous scavengers, H 2 O 2 can also accumulate in nitric solutions with an increased formation yield when compared to pure water due to a different radical mechanism. Sonolysis of various Pu solutions was studied in diluted aqueous nitric solutions (20 kHz, Ar, V= 50 mL, 0.26-0.52 W.mL -1 ) while varying the experimental conditions. The sonolysis of Pu(III) in 1 M HNO 3 allows to observe its oxidation to Pu(IV) resulting from HNO 2 generation by HNO 3 sonolysis. In these conditions, Pu(IV) cannot be reduced or oxidized to another specie. However, in the presence of anti-nitrous reagents, Pu(IV) is found to be reduced to Pu(III) with a rate increasing with the acoustic intensity.

  The concentration of Pu(VI) in aqueous nitric solutions in the presence of sulfamic acid was also found to decrease as a function of sonication time. The reduction allows to observe Pu(V) accumulation which disproportionate to Pu(IV) and Pu(VI). The reduction of Pu(VI) is however found to accelerate with sonication time which is explained by H 2 O 2 -reduction of Pu(IV) to Pu(III). Pu(III) may thereafter react with Pu(VI) in agreement with the well-known equilibrium Eq. 3 which explains the acceleration observed for the different species and allows the simultaneous observation of four different oxidation states of Pu in ~1 M HNO 3 .

		The extended sonication
	leads to an equilibrium state where Pu(III) and Pu(IV) coexist in various concentrations. This
	behaviour is attributed to the in-situ generation of H 2 O 2 in agreement with what observed without ultrasound (Eq.1-2). 1
	2 Pu 4+ + H 2 O 2 → 2 Pu 3+ + O 2 + 2 H + 2 Pu 3+ + H 2 O 2 + 2 H+ → 2 Pu 4+ + 2 H 2 O	(Eq. 1) (Eq. 2)
	PuO 2 2+ + Pu 3+ ↔ PuO 2 + + Pu 4+ (Eq. 3)

L. R. Morss, N. Edelstein and J. Fuger. The Chemistry of the Actinide and Transactinide Elements, Springer, Dordrecht,

2 K. S. Suslick and D. J. Flannigan, Ann.Rev. Phys. Chem. 2008, 59, 659.


Marc Dambrine

Djalil Kateb

ON THE WEAK MATERIAL APPROXIMATION IN LEVEL-SET METHODS

Keywords: 2000 Mathematics Subject Classification. Primary 49Q10; Secondary 34A55, 49Q12 shape optimization, stability, second order shape derivative, level-set method, weak material approximation

come

Introduction

Shape optimization consists in finding the shape of a domain which minimizes an objective function or a criterion. Classical objectives in mechanical engineering are for example, to maximize the torsional rigidity of a body or to minimize its compliance. Such criteria depend on the computation of a state function. In such cases, it usually solves a second-order elliptic partial differential equation that contains the physic under consideration.

Since the pioneering work of Hadamard, a shape calculus has been developed leading to shape derivatives and to subsequent optimization methods. Among the numerous difficulties one has to face during the design of such algorithms, the problem of the eventual change of the number of connected components of the current domain found a satisfactory answer thanks to the level set method introduced by S. Osher and J.A. Sethian ( [START_REF] Osher | Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF]). The idea is to define a domain of R d as the set of points where a real valued function defined on R d takes its non-possitive values. The geometry can then be discretized as a function on a regular grid and not thanks to control points. The flow of the shape gradient then formally leads to an Hamilton-Jacobi equation for the evolution of the domain coupled with the state function. Recently, P. Cardaliaguet and O. Ley have performed a first step in the theoretical justification of this approach on an example ( [START_REF] Cardaliaguet | Some flows in shape optimization[END_REF][START_REF] Cardaliaguet | On the energy of a flow arising in shape optimization[END_REF]. However, the level set grid has no particular reason to be useful for the computation of the state function.

In the context of mechanical structure optimization, the boundary value conditions on the part of the structure subject to optimization are usually traction free conditions. In [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF], G. Allaire, F. Jouve and A.-M. Toader suggest to take advantage of this particular boundary condition to use the level set grid for the numerical computation of the state.

The weak material approximation consists in filling the outside of the domain by a material of weaker conductivity and replace the state function by the solution of the bimaterial equation stated in the whole level set grid. The physical motivation is that the weak material behaves almost like an insulating layer leading to a reasonable approximation of the state that should be computed. This approximation gives very good results for all the usual criteria used in structure optimization and is now a reference method.

However, shape optimization problems are usually unstable. Homogenization is a well known example of unstable shape problem connected to a lack of existence of minimizers. However, even in the case where a smooth minimizer exists, the question of stability of such a shape is difficult. In the recent years, the analysis of order two conditions has been studied: A positive second order shape derivative at a critical shape does not imply stability. To sum up, two distinct behaviours of the numerical procedures of optimization have been encountered.

The first example under analysis was the Dirichlet energy. In [START_REF] Descloux | Stability of the solutions of the bidimensional magnetic shaping problem in abscence of surface tension[END_REF], J. Descloux showed that the two norms discrepancy phenomenon appears in a electromagnetic shaping problem: the shape hessian at a critical shape is coercive but only in strictly weaker norm than the norm of differentiability. In [START_REF] Dambrine | About stability of equilibrium shapes[END_REF], M. Dambrine and M. Pierre have proved that weak information was sufficient to insure stability at the continuous level. This first result was extended in [START_REF] Eppler | On convergence in elliptic shape optimization[END_REF] for a numerical scheme.

The second example comes from an inverse problem in electrical impedance tomography reformulated as a question of shape optimization. In [START_REF] Eppler | A regularized newton method in electrical impedance tomography using hessian information[END_REF][START_REF] Afraites | Detecting perfectly insulated obstacles by shape optimization techniques of order two[END_REF][START_REF] Afraites | On second order shape optimization methods for electrical impedance tomography[END_REF], more situations have been studied. The shape hessian at the global minimum is compact and the optimization procedures are severely ill-posed. Hence, appropriate strategies of regularization are required. Since the sequence of eigenvalues of the hessian tend towards zero, there exists a strong link between the number of degrees of freedom used to parametrize the shape and the numerical precision required for the computation of the state. Indeed, the number of shape parameters impose a threshold. Then, if the precision of the computed solution for the state is less than the threshold, the errors committed on the approximated criterion can destroy the well-posedness of the optimization problem.

In this paper, we present a model problem where the parameter of control is the shape of a structure which presents the worst behaviour. We are concerned by a shape optimization problem suggested by the design of micromechanisms introduced in [START_REF] Bernardoni | Outils et méthode de conception de structures mécaniques à déformations et actionnements répartis[END_REF] and studied in [START_REF] Allaire | A level-set method for vibration and multiple loads in structural optimization[END_REF][START_REF] De Gournay | Velocity extension for the leve-set method and multiple eigenvalues in shape optimization[END_REF]. The criterion under consideration has been observed to lead to numerical instabilities (see the comments on the gripping mechanism in [START_REF] De Gournay | Velocity extension for the leve-set method and multiple eigenvalues in shape optimization[END_REF]). However, we remain for simplicity in the scalar case. The first goal of the paper is to explain rigorously these instabilities by studying the second shape derivative of the criteria in a continuous setting. In particular, in our toy problem, the hessian at the absolute minimizer vanishes. As a consequence of this flatness of the objective to minimize, the weak material approximation is prohibited. Indeed, the systematic errors committed on the computation of the state generate instabilities. The second objective of this manuscript is to convince the reader of the interest of second order shape derivatives to understand the behaviour of numerical algorithms of minimization.

The paper is organized as follows. In Section 2, we present the model problem and state the results of this work. Section 3 is devoted to the computation of shape derivatives while the instabilities are discussed in section 4. The technical results of shape calculus that we need are postponed into an appendix.

2.

Setting of the problem and main results.

The shape optimization problem. Let D be a smooth domain of R d , d ≥ 2. We suppose that the boundary ∂D of D has three particular disjoint parts Γ c , Γ d and Γ o as illustrated in Figure 1. Concerning the admissible shapes Ω ∈ U ad , smooth C 2,α subdomain

Input flux Γ c Target potential u 0 Γ o Γ d Γ m D Ω Fixed Dirichlet boundary 1 Figure 1. Geometrical settings. of D, we will suppose that Γ c ∪ Γ d ∪ Γ o ⊂ ∂Ω.
Since those parts of ∂Ω are fixed, the part subject to optimization is the remainder

Γ m = ∂Ω \ (Γ c ∪ Γ d ∪ Γ o ).
In the sequel, we use the convention that a bold character denotes a vector. For example, n denotes the normal vector pointing outside of Ω. If h denotes a deformation field, it can be written as h = h τ + h n n on Γ m . Note that h τ is a vector while h n is a scalar quantity. The admissible deformation fields have to preserve the complement of Γ m in ∂Ω and therefore the space of admissible fields is

H = {h ∈ C 2,α (D, D), h n = 0 on Γ c ∪ Γ d ∪ Γ o }.
This problem recovers the following physical meaning: how should one design the structure to find the closest potential distribution to a reference one u 0 on Γ o while imposing the flux on another part of the boundary ? More precisely, the state function u solves the following boundary value problem:

(1)

         -∆ u = f in Ω u = 0 on Γ d , ∂ n u = 0 on Γ m ∪ Γ o , ∂ n u = g on Γ c .
where f and g are C ∞ . Both f and g do not live in the natural space but to a stronger one since additional regularity is required in order to define the shape derivatives. Once Γ c = ∅, this problem admits a unique solution u ∈ H 1 (Ω).

Let us describe our optimization problem.We consider the criterion J LS defined by

(2)

J LS (Ω) = 1 2 Γo |u Ω -u 0 | 2 dσ
where u 0 ∈ H 1 2 (Γ 0 ). The function u Ω solves (1) and hence depends on Ω. For the lightness of notations, we will omit this dependency and denote this function by u. We are interested in the minimization of J LS (Ω) with respect to the shape, where the minimum is taken on the set of admissible domains U ad . Consequently, the only optimized part of the shape boundary is Γ m . The Dirichlet boundary conditions are only here to ensure uniqueness of solutions.

For an arbitrary data u 0 , existence of a minimizer to J LS is not clear. Nevertheless, this question does not enter the scope of this work and we will focus on a well posed problem. We fix an admissible domain Ω * and choose for data u 0 exactly the trace on Γ o of the solution to the boundary value problem (1) set on Ω * . This assumption ensures that the optimization problem has a global minimizer in the admissible domains that is Ω * .

Differentiability results for the state u. To compute the shape derivatives, we have to recall some notions of shape optimization. The following result concerns the first order derivative of the state functions u.

         ∆u = 0 in Ω, u = 0 on Γ d , ∂ n u = 0 on Γ c ∪ Γ o , ∂ n u = div τ (h n ∇ τ u) + f h n on Γ m . (3)
(ii) If k ≥ 3, the state u has a second order shape derivative u = D 2 u(Ω; h 1 , h 2 ) in H 1 (Ω) solution of the boundary value problem (4)

                       ∆u = 0 in Ω, u = 0 on Γ d , ∂ n u = 0 on Γ c ∪ Γ m , ∂ n u = div τ (h 2,n ∇ τ (u) 1 + h 1,n ∇ τ (u) 2 + h 1τ .(Dn h 2τ )∇ τ u) -div τ ((h 1τ .∇ τ h 2,n + ∇ τ h 1,n .h 2τ )∇ τ u) +div τ (h 2,n h 1,n (2Dn -HI)∇ τ u) +h 1,n h 2,n ∂ n f + f (h 1τ Dnh 2τ -h 1τ .∇ τ h 2,n -h 2τ ∇ τ h 1,n ) , on Γ m .
In (4), (u) i denotes the first order derivative of u in the direction of h i as given in (3), Dn stands for the second fundamental form of the manifold Γ m and H stands for its mean curvature.

Differentiability of the objective. Once the differentiability of the state function has been established, the chain rule provides the differentiability with respect to the shape of the criterion. As usual for Least Squares objective, this derivative can be simplified thanks to an adjoint state that will be denoted by p.

Proposition 1. Let Ω be an open smooth subset of R d (d ≥ 2) with a C 2,α boundary.
Then, for all admissible directions h ∈ H, the shape derivative of the criteria DJ(∂Ω; h) is given by

(5) DJ LS (Ω; V ) = Γm (f p -∇ τ u.∇ τ p) h n
where the adjoint state p is solution of the adjoint problem

(6)          ∆ p = 0 in Ω, p = 0 on Γ d , ∂ n p = u -u 0 on Γ o , ∂ n p = 0 on Γ m ∪ Γ c .
One can also be interested by the computation of higher order derivatives and specially by the second derivative or shape hessian. To that end, we will need the shape derivative of the adjoint state p obtained as a consequence of Theorem 1.

Corollary 1. Let Ω be an open smooth subset of R d (d ≥ 2).
Then, the adjoint state function p is shape differentiable. Its shape derivative p = Dp(Ω; h) ∈ H 1 (Ω) solves the boundary value problem

(7)          ∆p = 0 in Ω p = 0 on Γ d , ∂ n p = 0 on Γ c ∪ Γ o , ∂ n p = div τ (h n ∇ τ p) on Γ m .
To investigate the properties of stability of this cost function, we now consider the quadratic form associated to the shape hessian at a critical shape Ω c . Theorem 2. Let Ω c be an open smooth subset of R d (d ≥ 2) with a C 3,α boundary. Let h 1 and h 2 be two deformation fields in H. The objective J LS is twice differentiable with respect to the shape and its second derivative in the direction h is given by:

(8) D 2 J LS (Ω c ; h, h) = Γm (f p -∇u .∇ τ p -∇u.∇ τ p ) h n + ∂ n (f p -∇ τ u.∇ τ p)h 2 n
Then, we specify Ω c : We consider the global minimizer Ω * . Claim 1. If Ω * realizes the absolute minimum of the criterion J LS , then

(9) D 2 J LS (Ω * , h, h) = 0.
Equation ( 9) means that the objective is very flat around the minimizer. This fact has two main consequences.

Claim 2. On the continuous criterion, the shape hessian at the global minimizer is not coercive. This means that this minimizer may not be a local strict minimum of the criterion.

Claim 3. The second consequence concerns any numerical scheme used to obtain this optimal domain Ω * . In order to capture the right behaviour of the shape gradient or hessian, one should compute with an extreme precision the state function and its derivative. In the case of a convex objective, the level of the threshold is naturally given by the size of the lowest eigenvalues of the discretized hessian at the critical shape. If the approximation of u inserted in the expression of the shape derivative is computed with a larger error, the computed hessian may have eigenvalues with a non positive real part. In such a case, the numerical descent schemes cannot converge since the discrete objective has no minimum but a saddle point.

The approximation by a weak material. One of the main advantage of the level set method in structure optimization is to avoid remeshing of the moving geometry at each step of the evolution. In [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF], this is achieved thanks to the weak material approximation. In mathematical terms, one approximates the solution u to the boundary value problem (1) by the solution u ε of the new boundary value problem [START_REF] Delfour | Shapes and Geometries : Analysis, Differential Calculus, and Optimization[END_REF]

         -div (σ ε ∇) u ε = f in D u ε = 0 on Γ d , ∂ n u ε = 0 on ∂D \ Γ c , ∂ n u ε = g on Γ c .
where D is a domain such that

(11) Ω ⊂ D and Γ d ∪ Γ o ∪ Γ c ⊂ ∂D.
The right hand side f is extended by zero outside of Ω and the distribution of conductivity σ ε is then such that (12)

σ ε = χ Ω + εχ D\Ω .
A very important question is then to know how good is the approximation of u by u ε ? By variational techniques, F. Murat and F. Jouve have shown ( [16]) that

u -u ε H 1 (Ω) ≤ C(Ω, D, f, g) √ ε.
In fact, numerical simulations show a better order of approximation. We tackle this question with asymptotic methods to recover the order ε.

Theorem 3. Let D and Ω be as in [START_REF] Descloux | Stability of the solutions of the bidimensional magnetic shaping problem in abscence of surface tension[END_REF]. Let u (resp. u ε ) be the solution of the boundary value problem (1) (resp. [START_REF] Delfour | Shapes and Geometries : Analysis, Differential Calculus, and Optimization[END_REF]). There is a constant C(Ω, D, f, g) and a value ε 0 > 0 such that

(13) u -u ε H 1 (Ω) ≤ C(Ω, D, f, g)ε.
rolls for all ε ∈ (0, ε 0 ).

The limit of this approach is then clear: in order to keep reasonable matrix to inverse, the value of the parameter ε cannot be to small: the problem becomes more and more badly conditioned when ε tends to 0. This means that the numerical errors committed on the computation of the state function are small but not small enough. With the view to our shape optimization problem with a flat objective, this lack of precision is a serious obstacle to the use of the method. The bad situation depicted in Claim 3 appears in the worst case: the hessian provides no possible error.

A better but still highly unstable case is the inverse problem presented in [START_REF] Afraites | Detecting perfectly insulated obstacles by shape optimization techniques of order two[END_REF] where the hessian with respect of the shape at the global minimizer is non-negative but compact. This phenomenon explains why a parametrized model of the inclusion is used in [START_REF] Afraites | Detecting perfectly insulated obstacles by shape optimization techniques of order two[END_REF]. It provides the threshold and the mesh used to compute the state function is then appropriately defined.

Justification of the shape derivatives.

The section is devoted to the proof of Theorems 1 and 2. We follow the usual strategy to prove differentiability in shape optimization. Computations made in this section require some classical facts in the context of shape calculus. We present them in Appendix A where we also introduce the notations.

3.1. Derivatives of the state function. Proof of Theorem 1. Proof of Theorem 1: First order derivatives. Existence of derivatives is well known in that case. Hence we only explain how to derive [START_REF] Allaire | A level-set method for vibration and multiple loads in structural optimization[END_REF]. Let H 1 Γ d (Ω) the subspace of H 1 (Ω) made of functions that vanishes on Γ d . We write the weak formulation of ( 1)

∀φ ∈ H 1 Γ d (Ω), Ω ∇u.∇φ + Γc gφ = Ω f φ.
We then differentiate these integrals with respect to the shape in the direction of h to

obtain ∀φ ∈ H 1 Γ d (Ω), Ω ∇u .∇φ = ∂Ω (f φ -∇u.∇φ) h n .
To be able to interpret the boundary integral as a boundary condition, we need to remove the derivative of φ in the right hand side. To that end, we use the boundary conditions and the definition of admissible deformation fields

∂Ω ∇u.∇φ h n = ∂Ω (∇ τ u.∇ τ φ + ∂ n u∂ n φ) h n = ∂Ω ∇ τ u.∇ τ φ h n .
Thanks to the integration by part formula (37), we get

∂Ω ∇ τ u.∇ τ φ h n = - Γm div τ (h n ∇ τ u) φ.
We conclude by density in L 2 (Γ m ) of the traces of the test functions φ. Second order derivatives. We compute the second derivative by considering two admissible deformations h 1 , h 2 ∈ H that will describe the small variations of ∂ω. J. Simon shows that the second derivative F (∂Ω; h 1 , h 2 ) of F (∂Ω) is defined as a bounded bilinear operator satisfying ( 14)

F (∂Ω; h 1 , h 2 ) = (F (∂Ω; h 1 )) h 2 -F (∂Ω; Dh 1 h 2 ).
For more details, the reader can consult the lecture of Murat and Simon or the book of Pierre and Henrot [START_REF] Henrot | Variation et optimisation de formes[END_REF]. We split the proof into two distinct parts: in a first time, we prove existence of the order two derivative, then, in a second time, we show that this derivative solves the boundary value problem (4).

First step: Existence of the second derivative.

Let us begin the proof. Let h 1 , h 2 ∈ H be two vector fields. The direction h 1 being fixed, we consider u1,h 2 the variation of u1 with respect to the direction h 2 . We recall that the material derivative u1 of u in the direction h 1 satisfies [START_REF] Hettlich | A second degree method for nonlinear inverse problems[END_REF] ∀v ∈ H 1 0 (Ω),

Ω ∇ u1 .∇v = Ω ∇u.A h 1 ∇v + Ω div (f h 1 ) v.
Let φ 2 : Ω → Ω be the diffeomorphism defined by φ 2 (x) = x + h 2 (x) and we set

ψ 2 = φ -1 2 . Setting ω h 2 = {x + h 2 (x), x ∈ ω}, Ω h 2 = {x + h 2 (x), x ∈ Ω} = Ω we get (16) Ω h 2 ∇ u1,h 2 .∇v = Ω h 2 ∇u h 2 .A h 1 ∇v + Ω h 2 div (f h 1 ) v
where u h 2 is the solution of the original problem with ω h 2 instead of ω. Making the change of variables x = φ 2 (X), we write the integral identity (16) on the fixed domain Ω

Ω σ∇ u1,h 2 . Dψ 2 (Dψ 2 ) T det(Dφ 2 ) ∇v = Ω ∇ũ h 2 . Dψ 2 A h 1 (Dψ 2 ) T det(Dφ 2 ) ∇v + Ω det(Dφ 2 ) div (f h 1 )v. ( 17
)
with the notations ũ = u • φ 2 and A h 1 = A h 1 • φ 2 . Since the material derivative u1 of u with respect to the direction h 1 satisfies

Ω ∇ u1 .∇v = Ω ∇u.A h 1 ∇v + Ω div (f h 1 ) v.
The difference of ( 15) and [START_REF] Kisch | The domain derivative and two applications in inverse scattering theory[END_REF] gives

Ω ∇ u1,h 2 -u1 .∇v = Ω ∇ u1,h 2 . I -Dψ 2 (Dψ 2 ) T det(Dφ 2 ) ∇v + Ω σ∇ũ h 2 . Dψ 2 A h 1 (Dψ 2 ) T det(Dφ 2 ) -A h 1 ∇v + Ω (∇ũ h 2 -∇u).A h 1 ∇v + Ω det(Dφ 2 ) div (f h 1 ) -div (f h 1 ) v.
We quote from [START_REF] Kisch | The domain derivative and two applications in inverse scattering theory[END_REF] and [START_REF] Hettlich | A second degree method for nonlinear inverse problems[END_REF] the following general asymptotic formulae that contains information on the deformations fields.

det(Dφ

i ) -1 -div (h i ) ∞ = O( h i 2 C 2 ), Dψ i (Dψ i ) T det(Dφ i ) -I + A h i ∞ = O( h i 2 C 2 ), div (h 1 ) det(Dφ 2 ) -div (h 1 ) -div (h 1 ) div (h 2 ) -∇div (h 1 ).h 2 ∞ = O( h 2 2 C 2 ), Dψ 2 A h 1 (Dψ 2 ) T det(Dφ 2 ) -A h 1 + Dh 2 A h 1 +A h 1 (Dh 2 ) T -div (h 2 ) A h 1 -(A h 1 ) (h 2 ) ∞ = O( h 2 2 C 2 )
. Making the adequate substitutions, we easily check that the material derivative of u1 with respect to h 2 exists. This derivative, denoted by ü1 , satisfies ( 18)

Ω ∇ü 1 .∇v dx = Ω ∇ u1 .A h 2 ∇v + ∇ u2 .A h 1 ∇v -∇u.A∇v + Ω div (h 2 div (f h 1 )) v.
where A is defined in (42). Second step: Derivation of (4) by formal differentiation of the boundary conditions.

The aim of this section is to retrieve the expression of the flux ∂ n u by computing the normal derivatives of each of the expressions ∇u .n and ḋiv τ (h 1,n ∇ τ u). Since

∇u .n = div τ (h 1,n ∇ τ u) + f h 1,n = h 1,n ∆ τ u + ∇ τ h 1,n .∇ τ u + f h 1,n
then, after taking the material derivative of the two sides of the above identity we obtain

∇u .n = ḋiv τ (h 1,n ∇ τ u) + ḟ h 1,n = ḣ 1,n ∆ τ u + h 1,n ∆ τ u + ∇τ h 1,n .∇ τ u + ∇ τ h 1,n . ∇ τ u + ḟ h 1,n + f ḣ 1,n . (19)
In order to avoid lengthy computations, we shall concentrate on each normal derivative appearing in the above formula. Some results are straightforward and their proof will be then left to the reader. Thanks to Proposition 4 , we conclude that

∇τ h 1,n = -∇ τ (h 1 .∇ τ h 2,n ) + (D 2 h 1,n .h 2 ) τ -∇h 1,n . ṅ n -∇h 1,n .n ṅ.
In the same manner, we also get

∇ τ u = ∇ τ u 2 + (σD 2 u.h 2 ) τ -∇ τ u. ṅ n -∇ τ un ṅ.
Hence, we can write

ḣ1 .n = h 2 .∇h ,n -∇ τ h 2,n .h 1τ ; ḋiv τ (h 1,n ∇ τ u) = ḣ 1,n ∆ τ u + h 1,n ∆ τ u + ∇τ h 1,n .∇ τ u + ∇ τ h 1,n . ∇ τ u, = -∇ τ (h 1 .∇ τ h 2,n ) + (D 2 h 1,n .h 2 ) τ -∇h 1,n .n ṅ ∇ τ u + h 1,n ∆ τ u +∇ τ h 1,n ∇ τ u 2 + (D 2 u.h 2 ) τ + (h 2 .∇h 1,n -∇ τ h 2,n .h 1τ ) .∆ τ u.
It remains to simplify the terms A = (D 2 u.h 2 ) τ .∇ τ h 1,n and B = ∇ τ u.(D 2 h 1,n .h 2 ) τ . We obtain:

A = -[σ∇ τ u].(Dn ∇ τ h 1,n )h 2,n + [σD 2 u]h 2τ .∇ τ h 1,n B = ∇ τ (∂ n h 1,n ).[σ∇ τ u]h 2,n -[σ∇ τ u].(Dn ∇ τ h 1,n )h 2,n + (D 2 h 1,n .h 2τ ).[σ∇ τ u].
We tackle the computation of (∂ n u ) . We first expand ḋiv

τ (h 1,n ∇ τ u). ḋiv τ (h 1,n ∇ τ u) = ḣ1,n ∆ τ u + ∇τ h 1,n .∇ τ u = ḣ 1,n ∆ τ u + h 1,n ∆ τ u + ∇τ h 1,n .∇ τ u + ∇ τ h 1,n . ∇ τ u.
After substitution, one gets

ḋiv τ (h 1,n ∇ τ u) =div τ (h 1,n ∇ τ u 2 ) + div ((h 2,n ∂ n h 1,n -∇ τ h 2,n .h 1τ )∇ τ u) + ∇ τ u.(D 2 h 1,n .h 2τ ) -∂ n h 1,n ∇ τ u.(Dn h 2τ ) -2h 2,n ∇ τ u.(Dn ∇ τ h 1,n ) + D 2 uh 2τ .∇ τ h 1,n + ∆ τ u∇ τ h 1,n .h 2τ + h 1,n ∆ τ u -∆ τ u 2 . ( 20
)
We then compute ∂n u 1 . From the expression of ṅ, we get after some straightforward computations:

(21) ∂n u 1 = ∂ n (u 1 ) 2 + (D 2 u 1 h 2 ).n + ∇ τ u 1 .(Dn h 2τ -∇ τ h 2,n ).
Finally, we compute ∂ n (u 1 ) 2 . We deduce from (20), (21) and formula

∇τ h 1,n = ḋiv τ (h 1,n ∇ τ u) + ḟ h 1,n + f ḣ 1,n
a first expression of the normal derivative of the second order shape derivative

∂ n (u 1 ) 2 = ḋiv τ (h 1,n ∇ τ u) -(D 2 u 1 h 2 ).n -∇ τ u 1 .(Dn h 2τ -∇ τ h 2,n ) = div τ (h 1,n ∇ τ u 2 + (h 2,n ∂ n h 1,n -∇ τ h 2,n .h 1τ ) ∇ τ u) -∂ n h 1,n ∇ τ u.(D n h 2τ ) + ∇ τ u 1 . (∇ τ h 2,n -Dn h 2τ ) + (D 2 h 1,n h 2τ ).∇ τ u + D 2 uh 2τ .∇ τ h 1,n + ḟ h 1,n + f ḣ 1,n -2h 2,n (Dn ∇ τ u).∇ τ h 1,n + ∇ τ h 1,n .h 2τ ∆ τ u + h 1,n ∆ τ u -∆ τ u 2 -(D 2 u 1 h 2 ).n. We have ḟ h 1,n + f ḣ 1,n = ∇f.h 2 + f (h 2 .∇h 1,n -∇ τ h 2,n .h 1τ ) . Taking account of the following calculation, -(D 2 u 1 h 2 ).n + ∇ τ u 1 .∇ τ h 2 , n = -h 2,n D 2 u 1 n + D 2 u 1 h 2τ .n + ∇ τ u 1 .∇ τ h 2 , n, = h 2,n (∆ τ u 1 + H∂ n u 1 ) + ∇ τ u 1 .∇ τ h 2,n -(D 2 u 1 h 2τ ).n, = div τ (h 2,n ∇ τ u 1 ) + Hh 2,n ∂ n u 1 -(D 2 u 1 h 2τ ).n; we rewrite ∂ n (u 1 ) 2 as ∂ n (u 1 ) 2 = div τ (h 1,n ∇ τ u 2 + h 2,n ∇ τ u 1 + (h 2,n ∂ n h 1,n -∇ τ h 2,n .h 1τ ) ∇ τ u) -(∇ τ u 1 + ∂ n h 1,n ∇ τ u) .(Dn h 2τ ) + Hh 2,n ∂ n u 1 -(D 2 u 1 h 2τ ).n + (D 2 h 1,n h 2τ ).∇ τ u + D 2 uh 2τ .∇ τ h 1,n -2h 2,n ∇ τ h 1,n .(Dn ∇ τ u) + ∇f.h 2 + ∆ τ u∇ τ h 1,n .h 2τ + h 1,n ∆ τ u -∆ τ u 2 + f (h 2 .∇h 1,n -∇ τ h 2,n .h 1τ ) . ( 22
)
This formula remains hard to handle. To get a more convenient one, we have to work a bit more. First, we derive tangentially to the direction h 2 the boundary identity

∂ n u 1 = h 1,n ∆ τ u + ∇ τ h 1,n .σ∇ τ u + f h 1,n
This leads to:

(D 2 u 1 h 2τ ).n+(Dn h 2τ ).∇ τ u 1 = ∇ τ h 1,n .h 2τ ∆ τ u + h 1,n ∇ τ ∆ τ u.h 2τ + (D 2 h 1,n h 2τ ).∇ τ u -∂ n h 1,n ∇ τ u.(Dn h 2τ ) + D 2 uh 2τ .∇ τ h 1,n + ∇ τ f.h 2τ h 1,n + f ∇ τ h 1,n .h 2τ . (23)
From ( 41) and subtracting ( 23) from ( 22), we can write

∂ n (u 1 ) 2 = div τ (h 1,n ∇ τ u 2 + h 2,n ∇ τ u 1 + (h 2,n ∂ n h 1,n -∇ τ h 2,n .h 1τ ) ∇ τ u) + div τ (h 1,n h 2,n (HI -2Dn).∇ τ u) -h 1,n (∇ τ ∆ τ u.h 2τ + ∆ τ ∇ τ u.h 2τ ) + h 1,n ∇ τ div τ (h 2τ ) .∇ τ u -div τ Dh 2 + (Dh 2 ) T ∇ τ u τ + ∂ n f h 2,n h 1,n + f (h 2,n ∂ n h 1,n -∇ τ h 2,n .h 1τ ) . (24)
From (41), we obtain

(25) ∆ τ u = ∆ τ u+∇ τ div τ (h 2τ ) .∇ τ u+∇ τ (Hh 2,n ).∇ τ u-div τ Dh 2 + (Dh 2 ) T ∇ τ u τ ,
and using the relation between the material and shape derivative, we get

∆ τ u = ∆ τ u + ∇ (∆ τ u) .h 2 and ∆ τ u = ∆ τ u + ∆ τ (∇u.h 2 ) .
Injecting these relations in (25) and applying them for h 2τ , we get

∆ τ (∇ τ u.h 2τ ) + ∇ τ div τ (h 2τ ) .∇ τ u = ∇ τ ∆ τ u.h 2τ + div τ Dh 2 + (Dh 2 ) T ∇ τ u τ .
Thanks to this last fact, expression (24) gives:

∂ n (u 1 ) 2 = div τ (h 1,n ∇ τ u 2 + h 2,n ∇ τ u 1 + (h 2,n ∂ n h 1,n -∇ τ h 2,n .h 1τ ) ∇ τ u) +div τ (h 1,n h 2,n (HI -2Dn).∇ τ u) + ∂ n f h 2,n h 1,n + f (h 2,n ∂ n h 1,n -∇ τ h 2,n .h 1τ ) .
To conclude, we use the following formula

∂ n u 1,2 = ∂ n (u 1 ) 2 -∂ n u Dh 1 h 2 .
where

∂ n u Dh 1 h 2 =div τ ((h 2,n n.∇h 1,n + ∇ τ h 1,n .h 2τ -h 1τ .Dn h 2τ )∇ τ u) + f (h 2,n ∂ n h 1,n + ∇ τ h 1,n .h 2τ -h 1τ .Dnh 2τ ) .
Finally, we obtain:

∂ n u 1,2 =div τ (h 2,n ∇ τ u 1 + h 1,n ∇ τ u 2 ) -div τ ((h 1τ .∇ τ h 2,n + ∇ τ h 1,n .h 2τ )∇ τ u) -div τ (h 2,n h 1,n (2Dn -HI)∇ τ u) + div τ (h 1τ .Dn h 2τ )∇ τ u) + ∂ n f h 2,n h 1,n -f (∇ τ h 2,n .h 1τ + ∇ τ h 1,n .h 2τ -h 2τ .Dnh 1τ )
Proof of Corollary 1:

It is completely similar to the proof of Theorem 1 without the volume right hand side f . Since Γ c remains fixed, no derivative of u appears in that boundary condition.

Derivatives of the criterion. Proof of Theorem 2. Proof of Proposition 1:

We differentiate the definition of the criterion on the fixed domain Γ c and get

DJ LS (Ω; h) = Γo u (u -u 0 ).
We introduce the adjoint state p defined in (6) to write

DJ LS (Ω; h) = Γo u ∂ n p = ∂Ω u ∂ n p.
Since both u and p are harmonic, a double integration by parts makes possible to permute the normal derivation and leads to

DJ LS (Ω; h) = ∂Ω ∂ n u p.
We insert the boundary conditions of (3) and use integration by part formula (36) to compute the gradient given in [START_REF] Bernardoni | Outils et méthode de conception de structures mécaniques à déformations et actionnements répartis[END_REF].

Proof of Theorem 2:

To compute the second derivative of DJ LS , we follow the classical scheme for computing the derivative for a boundary integral of a function ψ. The integral on Γ t m can be brought back from the moving boundary Γ t m to Γ m thanks to the change of variables

1 t I(Γ t m , h) -I(Γ m ) = 1 t Γm ψ • T t ω(t) -ψ,
where

ω(t) = |det(DT t )| (DT -1 t ) * n . The useful formula ω (0) = div τ (h) is stated in [10, 14].
In the case of J LS , we apply this approach to ψ = f p h n -∇ τ u.∇ τ p h n to get

d dt DJ LS (Ω; h) |t=0 = d dt Γm B(t)∇u t .∇p t + f p h • T t .n t |t=0 .
where we set

u t = u(Γ t m )•T t , p t = p(Γ t m )•T t , n t = n(Γ t m )•T t and B(t) = ω(t)DT -1 t (DT -1 t ) * . Taking account ṅ = -(Dh * ) τ and B (0) = div τ (h) I -Dh + (Dh) * ,
we split the derivative of the shape derivative into three terms

d dt DJ LS (Ω; h) |t=0 = A 1 (Ω; h) + A 2 (Ω; h) + A 3 (Ω; h)
where

A 1 (Ω; h) = - Γm [∇ τ u.∇ τ p + ∇ τ ṗ.∇ τ u + B (0)∇p.∇u] h n , A 2 (Ω; h) = Γm (f p -∇u.∇p) (Dhh.n -h.(Dh * n) τ ) , A 3 (Ω; h) = Γm (f p div τ (h) + f ṗ + p h.∇f ) h n = Γm (f p div τ (h) + f p + h.∇(pf )) h n .
At the optimum, we know that A 2 (Ω * ; h) = 0. Hence

D 2 J LS (Ω * ; h, h) = d dt DJ LS (Ω; h) |t=0 = A 1 (Ω * , h) + A 3 (Ω * , h), = - Γm (∇ τ u.∇ τ p + ∇ τ ṗ.∇ τ u + B (0)∇p.∇u) h n + (f p div τ (h) + f ṗ + ph.∇f ) h n .
Let us show how to simplify the expression. First of all, we write

∇ τ u.∇ τ p = ∇u + D 2 u h + (Dh) * ∇u .∇ τ p,
and we then gather all the terms so that

∇ τ u.∇ τ p + ∇ τ ṗ.∇ τ u = ∇u .∇ τ p + ∇u .∇ τ p + h.∇ (∇ τ u.∇ τ p) + (Dh + (Dh) * ) ∇p.∇u,
and then that

A 1 (Ω; h) = -∇u .∇ τ p -∇u .∇ τ p -h, ∇ (∇ τ u.∇ τ p) -div τ (h) ∇ τ u.∇ τ p.
We use Euler's Equation to remark that

h.∇ (f p -∇ τ u.∇ τ p) = h n ∂ n (f p -∇ τ u.∇ τ p) .
4. Analysis of stability.

4.1. Shape hessian at a critical shape. When Ω c is a critical shape, the shape hessian takes the expression

D 2 J LS (Ω c ; h, h) = Γm (∇u .∇ τ p + ∇u.∇ τ p -f p ) h n + ∂ n (f p -∇ τ u.∇ τ p)h 2 n
This expression is difficult to deal with since the adjoint state is not well known at an arbitrary critical shape. Let us consider the particular case Ω c = Ω * where the adjoint is simple.

Proof of Claim 1:

Since at Ω * , u = u 0 by definition, the adjoint state p cancels identically in Ω * and from Corollary 1 so does its derivative p . Inserting these facts in [START_REF] Dambrine | About stability of equilibrium shapes[END_REF], we get the claim. 4.2. On the weak material approximation. Our proof of Theorem 3 is inspired by Schwarz's method in domain decomposition. To that end, we introduce the operators

P : H 1/2 (Γ m ) → H 1 (D \ Ω), φ → v; and Q : H -1/2 (Γ m ) → H 1 (Ω), ψ → w;
where v and w solves the boundary value problems (26)

         ∆ v = 0 in D \ Ω, v = 0 on ∂Γ d , ∂ n v = 0 on Γ c ∪ ∂D, v = φ on Γ m .
and

         ∆ w = 0 in Ω, w = 0 on ∂Γ d , ∂ n w = 0 on Γ c , ∂ n w = ψ on Γ m .
The classic theory of elliptic PDE learns us that both P and Q are bounded. Let T Γm is the usual trace operator from H 1 (Ω) into H 1/2 (Γ m ). We then define the operators

P : H 1 (Ω) → H 1 (D \ Ω), u → P (T Γm (u)) ; and Q : H 1 (D \ Ω) → H 1 (Ω), v → Q(∂ n v).
Note that Q is not defined on the whole H 1 (D \ Ω) since the normal derivative ∂ n v should have a sense. We will only apply the operator Q to harmonic functions. In that case, the normal derivative is defined in the weak sense by duality. Moreover, if v is harmonic in

D \ Ω, then there is a constant C(D \ Ω) such that ∂ n v H -1/2 (∂Ω) ≤ C(D \ Ω) v H 1 (D\Ω)
and then

(27) Q(v) H 1 (Ω) ≤ C(D \ Ω) Q v H 1 (D\Ω) .
Proof of Theorem 3:

In fact, we will prove that u ε writes as an power series in the parameter ε. The idea is to split u ε into its values U ε on Ω and its values V ε on D \ Ω and then to use an double ansatz as follows

(28) U ε = ∞ i=0 ε i U ε,i and V ε = ∞ i=0 ε i V ε,i .
The error estimate ( 13) is then easily deduced from the convergence of the series defined in (28). Let us prove that these series really converge. We first determine the coefficients. Plugging these ansätze in the jump relations through the interface Γ m for the solution of the boundary value problem [START_REF] Delfour | Shapes and Geometries : Analysis, Differential Calculus, and Optimization[END_REF] [u ε ] = 0 and [σ ε ∂ n u ε ] = 0, we get

U ε = V ε ⇒ ∞ i=0 ε i U ε,i = ∞ i=0 ε i V ε,i , ∂ n U ε = ε∂ n V ε ⇒ ∞ i=0 ε i ∂ n U ε,i = ∞ i=1 ε i ∂ n V ε,i-1 .
By a term by term identification, we obtain ∂ n U ε,0 = 0 and the relationships

U ε,i = V ε,i , (29) ∂ n V ε,i = ∂ n U ε,i+1 . (30)
hold on Γ m for all i ≥ 1.

As fist consequence, U ε,0 solves the boundary value problem (1): it then holds U ε,0 = u. The term of order zero of V ε is then completely characterized thanks to (29) as solution of the boundary value problem (26) with φ = T Γm U ε,0 where T Γm is the trace operator from H 1 (Ω) into H 1/2 (Γ m ). We conclude that V ε,0 = P (U ε,0 ).

In fact, we obtain an iterative construction of the coefficients V ε,i and U ε,i in a Schwarz method's way. For all i ≥ 0, (29) implies that V ε,i = P (U ε,i ) and (30

) that U ε,i+1 = Q(V ε,i ).
The coefficients in (28) are then completely defined by

U ε,0 = u, U ε,i+1 = Q P (U ε,i ) for i ≥ 0; and V ε,0 = P (u), V ε,i+1 = P Q(V ε,i ) for i ≥ 0.
Since the operator P is bounded and thanks to (27), the expansions converge in the spaces H 1 (Ω) and H 1 (D\Ω) for ε small enough, for example such that ε < 1/ C(D \ Ω) P Q .

Appendix A. Elements of shape calculus

For the reader's convenience, we recall some basic facts from shape optimization without proof. We refer to [START_REF] Delfour | Shapes and Geometries : Analysis, Differential Calculus, and Optimization[END_REF][START_REF] Henrot | Variation et optimisation de formes[END_REF] for references and full proofs. Let h be a deformation field in C 2 (Ω, R d ) with h C 2 < 1. We set T t (h, .) = Id + th and denote by Ω t the transported domain Ω t = T t (Ω). It is well known that, for sufficiently small values of t, T t (h, .) is a diffeomorphism in R d . To avoid heavy notations, we will misuse the notation T t instead of T t (h, .).

The key concepts we need are material and shape derivatives. For any vector field h ∈ H, we define the material derivative of the domain functional y = y(Ω) at Ω in an admissible direction h as the limit

(31) ẏ(Ω; h) = lim t→0 y(Ω t ) • T t -y(Ω) t ,
Similarly, one can define the material derivative ẏ(∂Ω; h) for any domain functional y = y(∂Ω) which depends on ∂Ω. Another kind of derivative, called the shape derivative of y(Ω) in the direction h, occurs. It is viewed as a first local variation. Its definition is given by the following Definition 1. The shape derivative y = y (Ω; h) of a functional y(Ω) at Ω in the direction of a vector field h is given by (32) y = ẏ -h.∇y.

We will need in the sequel to manipulate the tangential differential operators on a manifold. For the reader's convenience, we recall from [START_REF] Delfour | Shapes and Geometries : Analysis, Differential Calculus, and Optimization[END_REF][START_REF] Henrot | Variation et optimisation de formes[END_REF] some definitions and also some useful rules of calculus.

Definition 2. The tangential divergence of a vector field

h ∈ C 1 (R d , R d ) is given by (33) div τ (h) = div (h) -Dh.n.n,
where the notation Dh denotes the Jacobian matrix of h. When the vector h ∈ C 1 (∂Ω, R d ) is defined on ∂Ω, then the following notation is used to define the tangential divergence

(34) div τ (h) = div h -(D h.n).n,
where h stands for an arbitrary C 1 extension of h on an open neighbourhood of ∂Ω.

We introduce the notion of tangential gradient ∇ τ of a smooth scalar function f in C 1 (∂Ω, R d ).

Definition 3. Let an element f ∈ C 1 (∂Ω, R d ) be given and let f be an extension of f in the sense that f ∈ C 1 (U ) and f | ∂Ω = f and where U is an open neighbourhood of ∂Ω. Then the following notation is used to defined the tangential gradient

(35) ∇ τ f = ∇ f | ∂Ω -∇ f .n n on ∂Ω.
The details for the existence of such an extension can be found in [START_REF] Delfour | Shapes and Geometries : Analysis, Differential Calculus, and Optimization[END_REF]. Let us remark that these definitions do not depend on the choice of the extension. Furthermore, one can show the important relation (36

) ∂Ω ∇ τ f.F = - ∂Ω f div τ (F) ,
for all elements f ∈ C 1 (∂Ω) and all vector fields

F ∈ C 1 (∂Ω, R d ) satisfying F n = F, n = 0.
In general, the condition above F n = 0 is not always satisfied. The extension of this integration by parts formula to fields with a normal vector component involves curvature. First, we point out that the curvature is connected to the normal vector via the tangential divergence operator. Recall that the mean curvature of ∂Ω is defined as H = div τ (n). Making use of the form of div τ (n) on the boundary, one shows straightforwardly the following statement. We assume now that the domain Ω has a C 3 boundary. The simplest second-order derivative is the Laplace Beltrami operator; it is defined as follows thanks to the following usual chain rule.

Definition 4. Let f ∈ H 2 (∂Ω). The Laplace-Beltrami ∆ τ of f is defined as follows

(38) ∆ τ f = div τ (∇ τ f ) .
There is a relation connecting the Laplace operator and the Laplace-Beltrami operator. Let us denote by ∂ 2 nn f = (D 2 f.n).n where D 2 f stands for the hessian of f . Proposition 3. Let Ω be a domain with a boundary ∂Ω of class C 3 . For all functions f ∈ H 3 (Ω), it holds (39) ∆f = ∆ τ f + H∂ n f + ∂ 2 nn f, on ∂Ω. We need to compute shape and material derivative of special vector fields: the outer unit normal vector n, the tangential gradient and the Laplace-Beltrami operator applied to a function. While the derivative of the normal vector is obtained by a straightforward calculus, we have to transport from ∂Ω t to ∂Ω the Laplace-Beltrami operator and the tangential gradient in order to compute the other derivatives. We recall here facts proved in [START_REF] Afraites | On second order shape optimization methods for electrical impedance tomography[END_REF] We denote by n the gradient of the signed distance to ∂Ω. This is an unitary extension of the unitary normal vector n at ∂Ω which is smooth in the vicinity of ∂Ω. This extension furnishes a symmetric Jacobian Dn that satisfies Dn n = 0 on ∂Ω. The direction h will be supposed to be in C 2 (R d , R d ) or in C 2 (∂Ω, R d ). We recall Proposition 4.

• The material derivative ṅ of the normal vector n at Ω in the direction of a vector field h ∈ C 1 (R d , R d ) is given by ṅ = -∇ τ (h.n) + Dn h τ , where h τ = hh.n n.

• The shape boundary n in the direction of h is given by n = -∇ τ (h.n).

• For all functions f ∈ C 2 (R 3 ) and directions h ∈ C 2 (∂Ω, R 3 ), one has

(40) ∇ τ f = ∇ ḟ + (D 2 f h) τ -∇f.n ṅ -∇f. ṅ n • Let f ∈ D(R d ).
The material derivative of ∆ τ f in the direction h is given by

∆ τ f = ∆ τ ḟ +∇ τ f.∇ τ [div τ (h τ )] + ∇ τ (Hh n ).∇ τ f -div τ Dh + (Dh) T ∇ τ f τ (41)
In the sequel, we will use some technical formulae given in [START_REF] Afraites | On second order shape optimization methods for electrical impedance tomography[END_REF]. Given a smooth vector field h, we denote Here, the notation (A h 1 ) (h 2 ) stands for the matrix defined by its elements ((A h 1 ) (h 2 )) k,l = ∇(((A h 1 ) ) k,l ).h 2

A h = Dh + Dh T -

Theorem 1 .

 1 Let Ω be an open smooth subset of R d (d ≥ 2) with a C k,α boundary. Let h and h 1 , h 2 be deformation fields in H.

( i )

 i If k ≥ 2, the state function u is shape differentiable. Its shape derivative u = Du(Ω; h) ∈ H 1 (Ω) solves the boundary value problem

Proposition 2 .

 2 Let Ω be an open subset of R 3 with a C 2 boundary. For any unitary extension N of n on a neighbourhood of ∂Ω, one has div (N ) = H on ∂Ω.Assume that the manifold ∂Ω has no borders. If F ∈ H 2 (∂Ω) 3 and f ∈ H 2 (∂Ω), then we have(37) ∂Ω ∇f.F + f div τ (F) = ∂Ω (∇f.n + Hf ) F.n.

  div (h) I Given two smooth vector fields h 1 and h 2 , we set(42) A = Dh 2 A h 1 + A h 1 Dh 2 T -A h 1 div (h 2 ) -(A h 1 ) (h 2 ), and b = (h 2 .∇u)A h 1 ∇v + (h 2 .∇v)A h 1 ∇u -((A h 1 ∇u).∇v)h 2 .

(

For any function u solution of (1) in Ω and for every test function φ ∈ D(Ω), we have (45) div (A h 1 ∇u) = ∆(h 1 .∇u) + div (f h 1 ) .