S Bordère

J.-P Caltagirone

A multi-physics and multi-time scale approach for modeling fluid-solid interaction and heat transfer

Keywords: multi-physics, multi-time scale, fluid-solid interaction, acoustic wave propagation, heat transport, monolithic approach

A novel non-conservative formulation for equations governing thermo-mechanical phenomena is developed to address multi-material and multi-physics issues. The first key point is that this formulation achieves a unifying equation for compressible viscous fluid flow and elastic solid deformation.

The second is that the thermo-mechanical equations are both written with velocity and thermal flux variables to solve them simultaneously. With that formulation, interaction conditions at the fluid-structure interface become implicit and state equation is no longer necessary. Multi-time scale problems are solved, from the time scale of acoustic and thermoacoustic wave propagation, to longer time scale of fluid flow and thermal diffusion.

I. INTRODUCTION

Many physical phenomena in materials science, biomedecine or process engineering are dependent on compressible or incompressible fluid mechanics, acoustics or non-linear wave propagation [START_REF] Lamb | Hydrodynamics[END_REF], heat transfer, elastic solid mechanics. These physical phenomena are very often coupled leading to complex problems such as solid particle behavior in a fluid flow [START_REF] Lorenz | Corrected momentum exchange method for lattice boltzmann simulations of suspension flow[END_REF][START_REF] Kanarska | Mesoscale simulations of particulate flows with parallel distributed lagrange multiplier technique[END_REF], plasma flow [START_REF] Caruyer | Modeling the first instant of the interaction between a liquid and a plasma jet with a compressible approach[END_REF], fluid-structure interaction for aerospace [START_REF] Farhat | Provably second-order time-accurate looselycoupled solution algorithms for transient nonlinear computational aeroelasticity[END_REF][START_REF] Cavagna | Application of navier-stokes simulations for aeroelastic stability assessment in transonic regime[END_REF], or biomedical applications [START_REF] Sankaranarayanan | Flow induced by acoustic streaming on surface-acoustic-wave devices and its application in biofouling removal: a computational study and comparisons to experiment[END_REF][START_REF] Barker | Scalable parallel methods for monolithic coupling in fluid-structure interaction with application to blood flow modeling[END_REF][START_REF] Vukicevic | Role of inertia in the interaction between oscillatory flow and a wall-mounted leaflet[END_REF][START_REF] Quaini | Validation of a 3d computational fluid-structure interaction model simulating flow through an elastic aperture[END_REF], thermo-elastic fracture [START_REF] Pant | Numerical simulation of thermo-elastic fracture problems using element free galerkin method[END_REF], transmission and reflection of elastic waves through an interface between dissimilar materials [START_REF] Singh | Reflection and transmission of couple longitudinal waves at a plane interface between two dissimilar half-spaces of thermo-elastic materials with voids[END_REF], thermoelastic shock wave diffraction with a multimaterial interface where fluid compressibility cannot be ignored [START_REF] Cavagna | Application of navier-stokes simulations for aeroelastic stability assessment in transonic regime[END_REF][START_REF] Sambasivan | Ghost fluid method for strong shock interactions part 1: Fluid-fluid interfaces[END_REF][START_REF] Sambasivan | Sharp interface simulations with local mesh refinement for multi-material dynamics in strongly shocked flows[END_REF][START_REF] Banks | Deforming composite grids for solving fluid structure problems[END_REF][START_REF] Breuer | Fluid-structure interaction using a partitioned semi-implicit predictor-corrector coupling scheme for the application of large-eddy simulation[END_REF].

The usual approaches for modeling such multi-physics problems are based on the standard equations that govern materials that often differ with the time scale considered. For instance, for short time scale characteristics in two-phase flows, thermal energy and material compressibility are considered through the Euler equations, whereas isothermal incompressible viscous flow is considered through Navier-Stokes equation for longer time scale steady-state flows.

Similarly, for temperature step variation, Fourier's equation and a compressible form of the Navier-Stokes equations are used to model thermoacoustic waves propagation. When a large discrepancy between the time scale constants of the coupled governing equations is present, for instance that of sound velocity and thermal diffusion, solving process needs to consider the time increment consistent with the lower time constant to avoid numerical instabilities, thus leading to time-consuming computation. Thus, for slow heating involving large time scale modeling, incompressible Navier-Stokes and heat transport equations are used since fluid flow convection and thermal diffusion time constants have closer order of magnitude.

Then, when the mechanical governing equations of the materials involved are not closely related which is the case for a fluid and an elastic solid, interface coupling also implies additional boundary equations at the interface. Simulation of such fluid-structure interaction problems called as fluid-structure interaction can be addressed through numerical methodologies which sequentially solve each governing equation with data transfer through the boundary conditions which can be temperature, pressure, shear stress or interface velocity depending on the problem [START_REF] Cavagna | Application of navier-stokes simulations for aeroelastic stability assessment in transonic regime[END_REF][START_REF] Breuer | Fluid-structure interaction using a partitioned semi-implicit predictor-corrector coupling scheme for the application of large-eddy simulation[END_REF][START_REF] Yigit | Grid movement techniques and their influence on laminar fluidstructure interaction computations[END_REF][START_REF] Degroote | Partitioned simulation of the interaction between an elastic structure and free surface flow[END_REF][START_REF] Degroote | Performance of partitioned procedures in fluid-structure interaction[END_REF][START_REF] Kassiotis | Partitioned solution to fluid-structure interaction problem in application to free-surface flows[END_REF][START_REF] Kassiotis | Nonlinear fluid-structure interaction problem. part i: implicit partitioned algorithm, nonlinear stability proof and validation examples[END_REF][START_REF] Jaiman | Combined interface boundary condition method for unsteady fluid-structure interaction[END_REF]. Some other approaches introduce a pseudo-simultaneous solving process for velocity and displacement fields in the fluid and solid domains, respectiveley, by implementing implicit iterative sequential coupling in order to converge at any time step and thus to avoid non-physical interface oscillations [START_REF] Barker | Scalable parallel methods for monolithic coupling in fluid-structure interaction with application to blood flow modeling[END_REF][START_REF] Quaini | Validation of a 3d computational fluid-structure interaction model simulating flow through an elastic aperture[END_REF]. This latter approach is defined as monolithic in contrast to the previous partitioned ones, but does not strictly solve system variables simultaneously, which will insure accuracy and stability of the interface coupling. The Arbitrary Lagrangian Eulerian (ALE) [START_REF] Farhat | Provably second-order time-accurate looselycoupled solution algorithms for transient nonlinear computational aeroelasticity[END_REF][START_REF] Sankaranarayanan | Flow induced by acoustic streaming on surface-acoustic-wave devices and its application in biofouling removal: a computational study and comparisons to experiment[END_REF][START_REF] Sambasivan | Sharp interface simulations with local mesh refinement for multi-material dynamics in strongly shocked flows[END_REF][START_REF] Huerta | Viscous flow with large free surface motion[END_REF][START_REF] Bathe | Finite element developments for general fluid flows with structural interactions[END_REF][START_REF] Turek | Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow[END_REF][START_REF] Bathe | A mesh adaptivity procedure for cfd and fluid-structure interactions[END_REF][START_REF] Murea | Arbitrary lagrangian eulerian approximation with remeshing for navier-stokes equations[END_REF][START_REF] Richter | Finite elements for fluid-structure interaction in ale and fully eulerian coordinates[END_REF] method which can be coupled to both monolithic and partitioned approaches has the advantage of discretizing strictly separated material. It can be used for viscous flow governed by Navier-Stokes equation as well as elastic solid equation. This method based on a moving mesh lead nevertheless to distorted elements giving rise to re-meshing developments for efficient calculations [START_REF] Bathe | A mesh adaptivity procedure for cfd and fluid-structure interactions[END_REF].

Finite elements methods have been largely improved in order to increase the convergence [START_REF] Kim | The finite element method enriched by interpolation covers[END_REF][START_REF] Kim | Towards a procedure to automatically improve finite element solutions by interpolation covers[END_REF] or to enhance the stress prediction [START_REF] Jeon | The mitc3 shell finite element enriched by interpolation covers[END_REF]. On the opposite, eulerian method based on a cartesian grid leads to a non-conforming interface, which has given rise to many numerical developments either within partitioned [START_REF] Banks | Deforming composite grids for solving fluid structure problems[END_REF][START_REF] Jaiman | Combined interface boundary condition method for unsteady fluid-structure interaction[END_REF] or monolithic [START_REF] Kloppel | Fluid-structure interaction for non-conforming interfaces based on a dual mortar formulation[END_REF][START_REF] Robinson-Mosher | A symmetric positive definite formulation for monolithic fluid structure interaction[END_REF] procedures in order to increase the accuracy of coupling the two materials at the interface. To overcome such inherent limitations of the numerical methods based on moving meshes and non-conforming interfaces, the method of finite spheres was recently developed [START_REF] De | The method of finite spheres[END_REF][START_REF] Hong | Coupling and enrichment schemes for finite element and finite sphere discretizations[END_REF] and applied to elastic wave propagation problems [START_REF] Ham | The method of finite spheres for wave propagation problems[END_REF].

Although numerical method developments or improvements are essential to deal with complex multi-material behaviors involving coupled governing equations, the development of new mathematical formulations is essential as well in order to avoid mismatch between governing equations for short and long time scale modeling, and to obtain strong fluid-structure coupling at the interface by solving strictly simultaneously the governing equations. An advance with such mathematical formulations was made in the particular case of isothermal multi-fluid systems involving either compressible or incompressible flows [START_REF] Caltagirone | A multiphase compressible model for the simulation of multiphase flows[END_REF] and for the modeling of trans-critical path from supercritical to subcritical states [START_REF] Amiroudine | A lagrangian-eulerian compressible model for the trans-critical path of near-critical fluids[END_REF]. In this paper we aim to develop the complete mathematical formulation for compressible and incompressible viscous flow, elastic solid deformation and heat transfer allowing to solve simultaneously for the velocity, displacement and temperature fields of the multi-material system. More precisely, this formulation has to combine the six main points discussed below.

• Compressible and "incompressible" fluids and solids are characterized by physical constants which are not strictly related, for instance Lamé's coefficients of an elastic solid and that of a fluid. The first point is thus to unify the material intrinsic properties which are involved in the conservative equations in order to reduce the different mechanical equations to a unique equation for strong coupling across the interface. This implies that every fluid will have to be considered as compressible or almost incompressible, in relation to its experimental thermodynamic coefficient.

• The second point is to be able to formulate the unified mechanical equation and the thermal equation with two vectorial variables in order to be able to solve for them simultaneously through a single linear system.

• The third point concerns time scaling. Obviously, it is the shorter time constants of the physical phenomena involved (acoustic, dynamic or thermal) that have to be considered. But, the formulation will have to give the possibility to obtain the solution using a higher time constant, allowing continuous evolution from a compressible fluid at the time scale of wave propagation to an incompressible fluid at the time scale of steady state evolution.

• As pointed out previously, the model has to deal with compressible fluids. In the literature highly compressible phenomena are resolved through a conservative formulation where the momentum ρV, may result in different types of discontinuities for density ρ and velocity V. While it is consistent to choose ρV as a variable when it varies slightly as for one-phase flows, this is no longer the case for two-phase flows where one phase is highly compressible and the other not. Indeed density can vary greatly across the interface whereas velocity varies continuously and slightly. So, the third point is to formulate a model with the velocity vector as variable, density being related to mass flux. Studies dealing with strongly two-phase shocked flows in the framework of a conservative formulation [START_REF] Sambasivan | Ghost fluid method for strong shock interactions part 1: Fluid-fluid interfaces[END_REF] prove to be insufficiently accurate to deal with two-phase continuous flows.

• If density is related to mass flux, we can discuss the suitability of using the state equation to determine density knowing temperature and pressure. Numerous numerical studies use state equations formulated to be valid over a broad interval of variable variation, particularly the pressure. Even in the case for the stiffened equation of state [START_REF] Banks | Deforming composite grids for solving fluid structure problems[END_REF][START_REF] Abgrall | On essentially non-oscillatory schemes on unstructured meshes: Analysis and implementation[END_REF][START_REF] Saurel | Numerical study of cavitation in the wake of a hypervelocity underwater projectile[END_REF] non-physical variable smoothing at the interface is obtained. The fourth point is thus to take the counterpart to the use of state equations by considering the thermodynamic coefficients. They can be more easily obtained from experimental data tables and semi-empirical laws and will not modify the formulation of the mathematical model as using the state equations does.

• Finally, the mathematical model proposed has to be independent of the numerical methodology allowing the use of either moving mesh or ALE methods where the different materials are strictly separated or Cartesian grids where interface tracking methods can be Lagrangian [START_REF] Tryggvason | A front-tracking method for the computations of multiphase flow[END_REF][START_REF] Raad | The three-dimensional eulerian-lagrangian marker and micro cell method for the simulation of free surface flows[END_REF][START_REF] Vincent | Eulerian-lagrangian multiscale meth-ods for solving scalar equations -application to incompressible two-phase flows[END_REF] or Eulerian [START_REF] Lopez | An improved plic-vof method for tracking thin fluid structures in incompressible two-phase flows[END_REF][START_REF] Osher | Level set methods: An overview and some recent results[END_REF].

II. UNIFYING MODEL FOR FLUID FLOW AND ELASTIC SOLID DEFORMA-TION COUPLED WITH HEAT TRANSFER

Let us consider a two-phase domain Ω delimited by a surface Γ. The interface between the two phases is noted Σ (Fig. 1). These two phases can be either isotropic elastic solids or Newtonian viscous fluids. The case where one phase is an elastic solid and the other a Newtonian viscous fluid is considered in the model as well. No mass exchange through the interface is involved, leading to a divariant system. The thermodynamic state can thus be described at any time t and any material point M with two intensive variables, for instance the absolute temperature T (M ) and the pressure p(M ). Density ρ(M ) and every thermodynamic coefficient such as specific heat c v (M ), isothermal compressibility coefficient χ T (M ), constant volume compressibility coefficient α(M ) and coefficient of thermal expansion β(M ) are thus entirely defined from these two variables. These thermodynamic coefficients are defined below together with the relationship between them Eq. 1. All these coefficients can be determined experimentally whatever the method used.

                           α = 1 p ∂p ∂T ρ β = - 1 ρ ∂ρ ∂T p χ T = 1 ρ ∂ρ ∂p T c p = ∂h ∂T p c v = ∂e ∂T ρ - ∂ρ ∂T p ∂T ∂p ρ ∂p ∂ρ T = β α p χ T = 1 (1)
The modeling of elastic materials requires a mechanical variable, the residual shear stress tensor τ (M ) (see paragraph 2.3.2).

For multi-material flows, either Lagrangian [START_REF] Maire | A second-order cell-centered lagrangian scheme for two-dimensional compressible flow problems[END_REF] or Euler-Lagrange [START_REF] Maire | A cell-centred arbitrary lagrangian-eulerian (ale) method[END_REF] approaches are used. Note, that the framework of our approach is to first develop the thermo-mechanical model within the Lagrangian representation. Thus, the material point temperature T (M ), pressure p(M ) and density ρ(M ) at time t will be determined from previously known values T 0 (M ), p 0 (M ) and ρ 0 (M ) at time t 0 (with t = t 0 +dt) and from the knowledge of the material derivative of these variables dX(M ) dt (with X = T, p, ρ), which depends on the applied velocity field V and heat flux Φ. It is in the final step that the Eulerian variables will be determined from the obtained Lagrangian variables.

A. Lagrangian form of pressure, temperature and density

As a first step, we focused on the Lagrangian formulation of pressure, temperature and density which can be established from their material derivatives. For simplicity, the indication of the material M -point for the variables is omitted.

The material derivative of density is obtained directly from the mass conservation equation:

dρ dt = -ρ ∇ • V (2)
The material derivative of temperature can be first expressed as follows.

dT dt = 1 c v de dt + 1 ρ 2 c v   T   ∂p ∂T   ρ -p   dρ dt (3)
Aiming to rewrite this equation as a function of velocity and flux vectors, we have: i) substituted for the material derivative of energy, the term defined from the conservation equation of the internal energy Eq. ( 4) ,

ρ de dt = -∇ • Φ -p∇ • V + q + φ (4)
where q and φ are the internal heat production and dissipation function, respectively;

ii) substituted for the material derivative of density the term of Eq. (2);

iii) substituted for the partial derivative of pressure with respect to temperature, the thermodynamic coefficients defined in the equation system (1) which leads to Eq. ( 5).

  ∂p ∂T   ρ = α p = β χ T (5)
We thus obtain the following equation.

dT dt = - β T ρ c v χ T ∇ • V - 1 ρ c v (∇ • Φ -q -φ) (6)
The material derivative of pressure can be written as a function of its partial derivatives as follows.

dp dt = ∂p ∂ρ T dρ dt + ∂p ∂T ρ dT dt (7)
To rewrite this equation as a function of velocity and flux vectors, we have: i) substituted for material derivatives of density and temperature the terms in equations

(2) and ( 6) respectively;

ii) substituted for partial derivative of pressure with respect to density, the thermodynamic coefficients defined in equation system (1) leading to Eq. ( 8);

∂p ∂ρ T = 1 ρ χ T (8)
iii) substituted for partial derivatives of pressure with respect to temperature, the thermodynamic coefficients in reference to Eq. ( 5).

We thus obtain the following equation.

dp dt = - 1 χ T + β 2 T ρ c v χ 2 T ∇ • V - β ρ c v χ T (∇ • Φ -q -φ) (9)
Finally, we can write the Lagrangian formulations of density, temperature and pressure at time t as a function of the divergences of velocity vector ∇ • V and heat flux vector ∇ • Φ,

                         ρ = ρ 0 e -dt ∇•V T = T 0 -dt β T 0 ρ 0 c v χ T ∇ • V -dt 1 ρ 0 c v (∇ • Φ -q -φ) p = p 0 -dt 1 χ T + β 2 T 0 ρ 0 c v χ 2 T ∇ • V -dt β ρ 0 c v χ T (∇ • Φ -q -φ) (10)
where the thermodynamic coefficients are defined at the M -material point with temperature T 0 and pressure p 0 at time t 0 .

Note that no state equation f (ρ, p, T ) = 0 is needed to calculate density, since it is defined exactly from the velocity vector divergence ∇ • V.

B. Vector form of the governing equations

The governing equations of the evolving system are composed of two equations [START_REF] Pant | Numerical simulation of thermo-elastic fracture problems using element free galerkin method[END_REF] which are the conservation of momentum ρV defined through Cauchy's equation and the conservation of heat flux Φ. These equations are written within the Lagrangian representation.

           ρ dV dt = ∇ • σ + ρ g + f τ p dΦ dt + Φ = -k ∇T (11)
Here, k is the thermal conductivity, g the gravity vector, f a volume force (for instance those induced by surface tensions) and σ is the stress tensor which depends on the rheological behavior of the materials involved. The parameter τ p in Eq. ( 11) corresponds to the time constant of phonon dynamics and has some importance in the numerical solution of hyperbolic systems in which time constants are very small (∼ 10 -9 s) [START_REF] Cattaneo | A form of heat equation which eliminates the paradox of instantaneous propagation[END_REF][START_REF] Vernotte | Some possible complications in the phenomena of thermal conduction[END_REF]. When the time constant is larger, the heat flux conservation equation converges towards the classical Fourier equation.

C. Elastic solid and fluid behaviors

The stress tensor σ can be split into spherical and non-spherical parts by introducing the pressure p and the shear stress tensor τ :

σ = -p I + τ (12)
The expression of the shear stress tensor depends on the material rheology involved.

Newtonian viscous behavior

For Newtonian viscous behavior, τ is expressed as a function of the shear strain rate tensor εs depending on the first-order derivative of the velocity vector V,

τ = 2µ εs = 2µ 1 2 (∇V + ∇ t V) - 1 3 ∇ • V I ( 13
)
where µ is the shear viscosity coefficient.

Isotropic elastic behavior

For an isotropic elastic solid, τ is expressed as a function of the shear strain tensor ε s depending on the first-order derivative of the displacement vector U in the case of small deformations,

τ = 2µ E ε s = 2µ E 1 2 (∇U + ∇ t U) - 1 3 ∇ • U I (14)
where µ E is the first Lamé's coefficient.

For an evolving system, the displacement U at time t can be rewritten as

U = U 0 + V dt
where U 0 is the displacement vector at time t 0 . Then, defining the shear stress τ 0 at time t 0 as a function of the U 0 -vector displacement using Eq. ( 14), the shear stress tensor τ can be written as a function of the velocity vector.

τ = τ 0 + 2µ E dt 1 2 ∇(V) + ∇ t (V) - 1 3 ∇ • (V) I (15)
D. Complete formulation of the governing equations

The final system (16) of the complete governing equations is constituted of the two conservative equations [START_REF] Pant | Numerical simulation of thermo-elastic fracture problems using element free galerkin method[END_REF] within which: i) the shear stress tensor at time t = t 0 + dt is introduced using the constitutive laws for viscous fluid Eq. ( 13) and elastic solid Eq. ( 15) which have been unified Eq. ( 17) using the coefficient µ * (µ * = µ for the Newtonian viscous fluid and µ * = µ E dt for the isotropic elastic solid), and the residual shear stress tensor τ 0 (that equals zero for Newtonian viscous fluids). The only equations which are solved are the two vectorial equations of system [START_REF] Breuer | Fluid-structure interaction using a partitioned semi-implicit predictor-corrector coupling scheme for the application of large-eddy simulation[END_REF];

ii) the temperature T and pressure p at time t are introduced using the Lagrangian formulation previously defined in system [START_REF] Quaini | Validation of a 3d computational fluid-structure interaction model simulating flow through an elastic aperture[END_REF]. The solution of the two equations corresponding to the unknown velocity field (or incremental displacement field) and heat flux field, also includes the pressure and temperature variations between the equilibrium state at time t 0 and that at time t 0 + dt, in thermodynamic consistency. The state equation for pressure determination is no more needed;

iii) only the density ρ 0 remains defined at time t 0 .

                             ρ 0 dV dt = -∇ p 0 -dt 1 χ T + β 2 T 0 ρ 0 c v χ 2 T ∇ • V -dt β ρ 0 c v χ T (∇ • Φ -q -φ) +ρ 0 g + f +∇ •   τ 0 + 2µ *   1 2 (∇V + t ∇V) - 1 3 ∇ • V I     τ p dΦ dt + Φ = -k ∇ T 0 -dt β T 0 ρ 0 c v χ T ∇ • V -dt 1 ρ 0 c v (∇ • Φ -q -φ) (16)
τ = τ 0 + 2µ *   1 2 ∇V + ∇ t V - 1 3 ∇ • V I   (17)
Since pressure p and temperature T are defined at time t within the governing equations using formulations of system [START_REF] Quaini | Validation of a 3d computational fluid-structure interaction model simulating flow through an elastic aperture[END_REF], these variables are thus simply updated from the resulting velocity and flux divergences using the as introduced formulations. Similarly, but for the elastic solid phase only, the shear stress is updated using the equation ( 17) (with µ * = µ E dt), the one that was introduced within the momentum conservation equation. Density ρ is updated, for its part, from a classical incremental process using the corresponding equation of system [START_REF] Quaini | Validation of a 3d computational fluid-structure interaction model simulating flow through an elastic aperture[END_REF].

Through the introduction of thermodynamics within the governing equations, original contribution of this model appears since temperature T 0 and pressure p 0 at time t 0 are usually considered in these equations written at time t. Since compressibility is clearly made explicit, modeling of "incompressible" fluids is done without using the additional equation ∇ • V = 0, but simply introducing the real isotherm compressibility coefficient. Low values of the compressibility coefficient (for instance χ T = 4.5 10 -10 P a -1 for water) will constrain the velocity divergence to approach zero.

Finally, as regards the dissipation function φ, for elastic solids φ = 0; for viscous fluids, where τ = τ 0 = 0 and µ * = µ it gives,

φ = 2µ (D - 1 3 tr(D) I) 2 (18)
where D = 1 2 (∇V + t ∇V).

E. Additional step for a Eulerian representation and system closure

Up to now, the model equations have been defined within a Lagrangian representation.

An Eulerian representation of the model, only requires one additional step which consists in the advection of the scalar Lagrangian variables X from their material derivative through the relation :

∂X ∂t = dX dt -V • ∇X (19)
where

X = ρ, T, p, φ, τ ij , V i (i = 1, 3; j = 1, 3).
If we differentiate the Lagrangian variables, again using the notation X(M ) previously defined from the Eulerian variables X(x) where x is the position of the geometric point, Eq. ( 19) gives,

X(x) -X 0 (x) dt = X(M) -X 0 (M) dt -V(M ) • ∇X(x) with X 0 (x) = X 0 (M) at time t 0 .
After simplification, the Eulerian variables are defined from the Lagrangian variables through the following equation:

X(x) = X(M ) -dt V(M ) • ∇X(x) (20)
The thermodynamic coefficients (χ T , β, c v ) and the physical constants (µ E , µ, k) are updated from the new Eulerian temperature T (x) and pressure p(x), as well as the heat production q(x).

III. FOCUS ON THE UNIFIED MECHANICAL MODEL FOR FLUID FLOW AND ELASTIC SOLID DEFORMATION

In the framework of this complete thermo-mechanical compressible model, this paper first focuses on the consistency of the unified modeling of compressible or "incompressible" fluid flow and elastic solid deformation. Within this purpose, isothermal conditions with no consideration of gravity and volume force are studied, which reduced the Lagrangian equation system ( 16) to the simplified unified equation defined below. Pressure update Eq. ( 10) is also simplified, density and shear stress updates, Eqs. ( 10) and ( 17), remaining unchanged.

                                                   ρ 0 dV dt = -∇ p 0 - dt χ T ∇ • V + ∇ •   τ 0 + 2µ *   1 2 (∇V + ∇ t V) - 1 3 ∇ • V I     p = p 0 - dt χ T ∇ • V ρ = ρ 0 e -dt ∇•V τ = τ 0 + 2µ E dt   1 2 (∇V + ∇ t V) - 1 3 ∇ • V I   for elastic solid only (21)
The consistency of the unified model will be demonstrated through three steps. As a first step, only mono-phase systems will be considered in order to validate the equation either for compressible fluids or elastic solids. In order to exhibit the physical contribution of the time characteristic statement dt in the equation, two conditions of time characteristics will be considered, dt > L /c and dt L /c, in relation to the characteristic length of the system L and the sound velocity c in the material involved.

As a second step, a two-phase system, involving two elastic solids, will be considered in order to evaluate how accurately the model deals with force equilibrium at the interface.

Then as the third step, a fluid-elastic solid system will be considered to show that this single equation can deal with fluid-structure interaction problems involving incompressible (almost incompressible) fluids. It is important here to explain why this formulation can be substantially used to model incompressible fluids since compressibility is introduced in the formulation. In fluid mechanics, two main types of methods are used to model incompressible flows. First the projection methods [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] of predictor-corrector type, where the Navier-Stokes equation is resolved first without any penalty, and then the solution projected on a free divergence field. The Vector Penalty Projection [START_REF] Angot | Fast discrete helmholtz-hodge decompositions in bounded domains[END_REF], corresponding to a vectorial scheme, is particularly efficient for the resolution of variable density flows. In the second type of resolution methods, the incompressibility penalty is introduced within the discretized momentum conservation equation. One of the more powerful methods is the Augmented Lagrangian method of Fortin and Glowinski [START_REF] Fortin | Méthodes de lagrangien augmenté. Application à la résolution numérique de problèmes aux limites[END_REF][START_REF] Fortin | Chapter 1 augmented lagrangian methods in quadratic programming[END_REF][START_REF] Fortin | Chapter 2 application to the stokes and navier-stokes equations[END_REF], where the equation and penalty, introduced through the term ∇(-r∇ • V ) , are solved simultaneously within the same vectorial system. Subsequent improvements have made this method competitive with the projections methods for solving the velocity-pressure coupling [START_REF] Caltagirone | Discrete Mechanics[END_REF]. For large values of the r-parameter the divergence is constrained to be close to the computer precision. In 2011, in the framework of multi-fluids flows only and isothermal conditions, it is thermodynamics which is introduced for the first time within the momentum conservation equation through the formulation of the pressure p (p = p 0 -dt /χ T ∇ • V ) giving rise to the physical meaning of the r-parameter (r = dt /χ T ) [START_REF] Caltagirone | A multiphase compressible model for the simulation of multiphase flows[END_REF]. This term also includes implicit mass conservation (see the development of the general lagrangian formulation of pressure in section II paragraph A where the mass conservation equation is used). The low experimental values of χ T for "incompressible" media corresponding to high dt /χ T -values thus insures the compressibility condition at the computer precision.

A. Numerical method and discretization of the Lagrangian-Eulerian formulation

The numerical methods that can be used to solve coupled equations in the framework of multi-material domains are numerous. The phase differentiation can be described by numerical schemes which strictly separate the phases, such as the ALE methods, or by Eulerian methods through a phase field function, such as the volume of fluid methods [START_REF] Breuer | Fluid-structure interaction using a partitioned semi-implicit predictor-corrector coupling scheme for the application of large-eddy simulation[END_REF] and other methods of interface tracking.

The present formulation could be solved using various numerical schemes since the physical model is written in a two-step formulation: 1) a lagrangian formulation step described in section II.D allowing the development of moving mesh methods (for instance ALE method); and 2) an additional formulation step described in section II-E to determine the eulerian variables from the advection of the lagrangian variables (Eq. ( 20)) allowing the development of fixed grid methods.

Discretization of the lagrangian momentum conservation equation

In this paper, the unified equation Eq. ( 21) was discretized in time and space by an implicit volume method on a staggered grid. This choice was first motivated by the numerical expertise we had in the discretization of the Augmented Lagrangian method [START_REF] Caltagirone | Discrete Mechanics[END_REF] that was extended to our formulation. This numerical scheme is certainly not the more efficient to deal with sharp interface coupling, but is nevertheless powerful to deal with compressible and incompressible media. Indeed, for finite volume schemes, we have to distinguish collocated grid from staggered one. When the unknown variables are located at the same point of the grid (collocated), some numerical instabilities can be observed since the calculation of pressure gradient results from non-adjacent points. In that case, some pressure smoothing and other numerical techniques are necessary such as the well-known one of Rhie and Chow [START_REF] Rhie | Numerical study of the turbulent flow past an air foil with trailing edge separation[END_REF]. Harlow and Welch [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF] were the first to locate the velocity components between two close neighboring pressure points. The staggered grid is powerful to deal with high pressure gradients and multi-phase flows [START_REF] Meier | Comparison between staggered and collocated grids in the finitevolume method performance for single and multi-phase flows[END_REF]. Staggered grid is also useful in this article in order to calculate velocity divergence using adjacent points.

In 2D, a direct solver MUMPS was chosen to solve the linear system. In 3D, an iterative biconjugate gradient stabilized BICG-StabII [START_REF] Vincent | An adaptative augmented lagrangian method for three-dimensional multimaterial flows[END_REF] was implemented. The spatial discretization employed centered schemes. This multiphase model is independent of the way the interface is located, which enables the use of any Lagrangian or Eulerian interface tracking method.

Here, the Eulerian Volume of Phase method (VOF) [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF][START_REF] Scardovelli | Interface reconstruction with least-square fit and split eulerianlagrangian advection[END_REF][START_REF] Aniszewski | Volume of fluid (vof) type advection methods in twophase flow: A comparative study[END_REF] is used, where the two-phase system is numerically built using a phase function C varying within the interval [0, 1].

The value of the phase C -function is 0 in one phase and 1 in the other. The interface between the two phases is defined as C = 0.5. Every physical constant A of the system is defined in relation with the C-function in two ways depending on the problem: i) the discontinuous way when the interface displacement is small, where A = A (0) for C ≤ 0.5 and A = A (1) for C > 0.5 (the exponents (0) and ( 1) are assigned to the first and second phases, respectively); and ii) the continuous way when the interface displacement is significant, where (1) .

A = (1 -C) A (0) + C A
As a first step, the intermediate unknown Lagrangian velocity field V at time step n + 1 is determined from the Eulerian variables at the time step n, after a time increment ∆t, using the following time discretization of Eq. ( 21).

ρ n V -V n t = -∇   p n -∆t 1 χ n T ∇ • V   + ∇ • τ n + 2µ * 1 2 ∇V + ∇ t V - 1 3 ∇ • V I (22)
The determination of Lagrangian velocities V from Eq. ( 22) only requires the knowledge of Eulerian density ρ n , pressure p n , residual shear stress τ n and isotherm compressibility coefficient χ n T at time step n. The intermediate unknown Lagrangian variables p, ρ and τ can then be determined from the knowledge of the velocity divergence ∇ • V through the discretization of the pressure update equation defined in the equation system [START_REF] Kassiotis | Nonlinear fluid-structure interaction problem. part i: implicit partitioned algorithm, nonlinear stability proof and validation examples[END_REF].

                                   p = p n -∆t 1 χ n T ∇ • V ρ = ρ n e -∆t ∇•V τ = τ n + 2µ *   1 2 (∇V + ∇ t V) - 1 3 ∇ • V I   (23)
This first order updating of pressure is exact only when the isotherm compressibility term χ T is independent of pressure. Otherwise, second order accuracy can be obtained in determining the Lagrangian pressure at time step n+1 by substituting for p n its development at time step n but using χ n T instead of χ n-1 T . It gives :

p = p n-1 -t 1 χ n T (∇ • V + ∇ • V n ) (24)

Discretization of the advection equations

Finally the intermediate scalar Lagrangian variables X (X = ρ, p, τ ij , V i with i = 1, 3; j = 1, 3) are advected using time discretization of Eq. ( 20) to determine the Eulerian

variables X n+1 (X n+1 = ρ n+1 , p n+1 , τ n+1 ij , V n+1 i with i = 1, 3; j = 1, 3) at time step n + 1,
which leads to the following equation.

X n+1 = X -t V • ∇X n+1 (25)
A Lax-Wendroff TVD scheme with a Superbee limiter for spatial discretization is used to accurately deal with sharp variations of scalar variables at the interface.

The phase function C n+1 at time step n + 1 is determined from the advection of the phase function C n at time step n in relation to the Eulerian velocity field V n+1 (Eq. ( 26))

using the same time and space discretization schemes as previously defined for scalar variable advection.

C n+1 = C n -t V n+1 • ∇C n (26)
The physical constants (χ n+1 T , µ * n+1 ) are then updated at time step n + 1 from the C n+1 function within the discontinuous or the continuous criteria at the interface described above, and from th values of the intensive variables T n+1 and p n+1 .

B. Isothermal one-phase compression

The relevance of the unified equation for fluid flow and elastic solid deformation is first demonstrated through basic problems of fluids and elastic solid compressions and simple shear deformations of elastic solids. They are simulated using a cubic cavity of length side L=1 m containing either one or two phases at initial pressure p i and density ρ i .

Compressions are simulated by injecting at a velocity V 0 , normal to the lower cavity surface, the material involved at this surface (Fig. 2a). Then, the system volume remains constant, it is the density ρ, pressure p and shear stress tensor τ which vary. The velocity field is kept at zero value at the upper surface of the cavity, and symmetric conditions are set at the four lateral surfaces (Fig. 2a). For simulation of simple shear deformation, the lower surface is displaced with a velocity V 0 collinear to the X-axis. The velocity field is maintained at zero value at the upper surface and periodic conditions are set at the four lateral surfaces (Fig. 2b). The physical characteristics of fluid and solids are shown in table I. For solids, Lamé's 1st and 2nd coefficients µ E and λ E are calculated using Young's modulus E and Poisson's coefficient ν (Eq. [START_REF] Murea | Arbitrary lagrangian eulerian approximation with remeshing for navier-stokes equations[END_REF].

             µ E = E 2(1 + ν) λ E = Eν (1 + ν)(1 -2ν) (27)
The coefficient of isothermal compressibility χ T is deduced from µ E and λ E coefficients using Eq. [START_REF] Richter | Finite elements for fluid-structure interaction in ale and fully eulerian coordinates[END_REF].

1

χ T = 2µ E + 3λ E 3 ( 28
)
The physical constants of Table I show that cork is a perfect compressible solid (ν = 0), and rubber a quasi-incompressible solid (ν ≈ 0.5) with a coefficient of isothermal compressibility close to that of water. Isothermal compressibility coefficient,χ T (Pa -1 )

9.87 × 10 -6 6.00 × 10 -7 6.00 × 10 -10 TABLE I: Physical constants of the three materials considered for the basic compression and shear deformation problems.

Analytical solution

The one-phase compressions within the cubic cavity (Fig. 2a) proceed through constant injection velocity V 0 . Analytical solutions can be obtained when steady compression is assumed which leads to uniform velocity divergence ∇ • V = -V 0 /L within the overall system with time. In these conditions, the integration of density conservative equation (Eq. ( 2))

gives the evolution of mass density with injection time at any point of the system:

ρ = ρ i e V 0 L (t-t i ) (29)
where ρ i is density at initial time t i .

The evolution of pressure depends on the phase involved through the coefficient of isothermal compressibility χ T .

For elastic solids, χ T is assumed to be constant with compression pressure. This approximation leads, when integrating Eq. ( 9) simplified for isothermal conditions, to a linear dependence of pressure with time through the relation :

p -p i = V 0 χ T L (t -t i ) (30)
where p i is the pressure of the elastic solid at initial time t i .

For fluids, considered to be perfect gasses, the coefficient of isothermal compressibility χ T can be substituted by the reverse pressure term 1/p, which leads, after integration of Eq. ( 9) simplified for isothermal conditions, to the following evolution of pressure with time.

p = p i e V 0 L (t-t i ) (31)
For elastic solids only, the evolution of the shear stress tensor τ with compression time gives,

τ = τ i + 2 µ E      1 3 V 0 L 0 0 0 -2 3 V 0 L 0 0 0 1 3 V 0 L      (t -t i ) ( 32
)
where τ i is the shear stress tensor at time t i .

A normalized shear stress tensor is defined as τ = 3 2 τ : τ , leading to scalar evolution of the shear stress with time :

τ -τ i = 2µ E V 0 L (t -t i ) ( 33
)
where τ i is the normalized shear stress at initial time t i . For viscous flow, τ is reduced to zero at any time.

Numerical simulation

Compressions of air, cork and rubber are compared until injection time t f =100 s by considering the same injection speed V 0 = 0.001 m/s. The 3D simulations were carried out on a 21 x 21 x 21 mesh grid for different time increments t (100 s, 10 s, 1 s, 0.1 s, 0.01 s).

During compression, density, pressure and normalized shear stress (for elastic solids only) increase over time, and are verified to be uniform throughout the entire volume even for the first time step resolution. The comparison between the numerical evolution and the analytical evolution of relative ρ/ρ i density (Eq. ( 29)) is presented in Fig. 3 for the three materials studied. The evolution is shown to be independent of the material involved as predicted by Eq. ( 29). Moreover, the simulated evolution strictly fits the exponential analytical one,

showing that the model accurately accounts for density variation by solving velocity divergence. A state equation is thus no longer needed with this compressible model whatever the materials involved. Similarly, the evolutions of pressure and normalized shear stress during compression for the elastic materials perfectly match the analytical ones in reference to Eq. ( 31) and Eq.

(33), respectively, as shown in Fig. 4 for air and Fig. 5 for cork. For air, since the coefficient of isothermal compressibility is pressure dependent, the evolution of pressure with compression time deviates from linearity (Fig. 4). The update of pressure using Eq. ( 24) gives, as expected, a second order convergence with decreasing time increment (fig. 6). A relative error of around 10 -8 was obtained for the lowest time increment t = 0.01 s.

FIG. 6: Second order convergence rate of numerical air pressure towards the analytical solution with time increment.

In contrast, for elastic material, the coefficient of isothermal compressibility is assumed to be pressure independent, giving rise to linear evolutions of pressure and normalized shear stress with time (Fig. 5). The relative errors obtained with these linear evolutions are presented in table II for time increment t = 1 s. These relative errors are close to computer precision showing that, as expected, exact solutions for pressure and normalized shear stress are obtained (see section 4). Note that the steady compressions obtained were simulated with time increments greater than the characteristic time of the system, L /c L , where c L is the velocity of longitudinal elastic wave propagation. The velocity values for cork and rubber are indicated in table IV, the value for air being c L = 1 / √ ρχ T = 293 m/s . So, within this time criterion t > L /c L , compression modeling is independent of time discretization, apart from the numerical precision when the coefficient of isothermal compressibility is pressure dependent.

C. Isothermal simple shear deformation of elastic solids

Analytical solution

The shear deformation of elastic solids proceeds through the constant velocity displacement V 0 of the lower surface in the x-axis direction (Fig. 2b). Density and pressure do not vary during shear deformation. Assuming steady evolution leading to a uniform velocity gradient within the overall system, the time integration of the shear stress tensor gives,

τ = µ E      0 -V 0 L 0 -V 0 L 0 0 0 0 0      (t -t i ) (34)
from which we can determine the evolution of the normalized shear stress tensor,

τ -τ i = √ 3µ E V 0 L (t -t i ) ( 35
)
where τ i is the normalized shear stress at initial time t i . For viscous flows, τ is reduced to zero at any time.

Numerical simulation

The 3D simulations of shear stress deformations were carried out on a 21 x 21 x 21 mesh grid, for a constant velocity V 0 = 0.001 m/s and for different time increments t (100 s, 10 s, 1 s). The physical constants of elastic solids and initial conditions used are indicated in table I. During deformation, the normalized shear stresses of cork and rubber are verified to be uniform throughout the material. The simulated evolutions of shear stress within cork and rubber systems are presented in Fig. 7 until the deformation time t f =100 s. They are shown to perfectly match the exact evolutions given by Eq. ( 35). At the final time t f =100 s, exact solutions are obtained for shear deformations in reference to the low relative error (table III). Note that, as for compressions, the steady shear deformations obtained are modeled with time increments greater than the time characteristics of the system, i.e. t > L /c T , where c T is the velocity of the transverse elastic wave propagation (see Table IV TABLE III: Normalized shear stress calculated within the material for the final deformation time t f =100 s ; corresponding relative errors obtained for time increment t = 1 s.

D. Elastic wave propagation

In order to illustrate the physical meaning of the time dependent term dt/χ T in Eq. ( 21), compression and shear deformations are simulated by considering time increments much lower than the system time characteristics t L /C. We can thus expect to get some insight into the time scale of wave propagation.

Analytical solution

Elastic wave propagation results in the oscillating exchange of kinetic energy into elastic energy. For isotropic elasticity, elastic waves are the addition of two types of waves: longitudinal waves where the material displacement field is collinear with the direction of wave propagation; and the transverse waves where the displacement field is perpendicular to the direction of wave propagation. Wave propagation velocity depends on the elastic constants of the material, and the type of waves. For longitudinal waves, the theoretical propagation velocity c th L is,

c th L = 2µ E + λ E ρ ( 36
)
and for transverse waves, the theoretical propagation velocity C th T is defined as follows.

c th T = µ E ρ ( 37
)

Numerical simulation

The 3D simulations of elastic wave propagation were carried out on a 4 x 64 x 4 mesh grid, considering cork and rubber as elastic materials.

As a first step, the propagation of a compression pulse was studied by applying an injection velocity V 0 = 0.0001 m/s during the first time increment only, the injection velocity being equal to zero during the remaining simulation time. For time increment t = 3 × 10 -5 s much lower than the characteristic time t = L /c th L = 1.45 10 -3 s a non-steady evolution of the pressure was obtained as shown in Fig. 8 for rubber. Indeed, we can observe that the pressure at the center point of the system evolves through periodic oscillations which become attenuated with time until the equilibrium state. This is the result of the initial wave compression which propagates in the direction of the displacement field (Y-direction), reflects off the upper surface, then propagates backwards to reflect off the lower surface, and so on until the equilibrium state. The oscillation period T L corresponds to the time lapse required for the wave to pass through the (Y = 0.5)-plane of the material and return to it after being reflected, that is the time for the wave to cover the distance L. From the measurement of the period (Fig. 8), the longitudinal propagation velocity can be calculated through the relation c cal L = L T L . The calculated value for rubber is shown in table IV together with that obtained with cork as a perfectly compressible material.

FIG. 8: Pressure variation at the center point of the rubber system during longitudinal elastic wave propagation; time increment, t = 3 × 10 -5 s.

In the second step, the propagation of a shear deformation pulse was studied by displacing the lower material surface in the X-direction at the velocity V 0 = 0.0001 m/s during the first time increment only. Here again, the time increment t = 2×10 -3 s is much smaller than the characteristic time t = L /c th T = 0, 19 s which allowed wave propagation to be followed. So, the shear wave propagates in the Y-direction perpendicularly to the displacement field, which characterizes transverse wave propagation. Fig. 9 presents the evolution of the τ xy -element of the shear stress tensor at the center point of the rubber system. The periodicity of the curve shows, here again, that the transverse wave reflects successively off the upper and lower surfaces of the cavity. From the measurement of the oscillation period T T , the transverse propagation velocity can be calculated through the relation c cal T = L T T . The calculated value is indicated in table IV together with that obtained on cork.

FIG. 9: Variation, during transverse wave propagation, of the τ xy -element of the shear stress tensor at the center point of the rubber system; time increment, t = 2 × 10 -3 s.

The good correspondence between the simulated wave propagation velocities compared to the calculated ones (table IV) clearly shows that the model proposed perfectly solves elasticity behavior at the time scale of elastic wave propagation. As expected, the longitudinal wave velocities are greater than the transverse ones, especially for quasi incompressible materials where the value of the 2nd Lamé coefficient λ E is high (or the Poisson coefficient ν tends to 0.5) which is the case for rubber. Two elastic solids are stacked within the cavity (fig. 2a) in such a way that the interface is parallel to the XZ-plane. The thickness of the lower and upper materials are thus defined by xL and (1 -x)L, respectively with 0 < x < 1. The injection velocity V 0 is constant during compression and we assume a steady state leading to uniform displacement divergences ∇ • U (1) and ∇ • U (2) within phases ( 1) and ( 2) at each time t. Indexes 1 and 2 are assigned to the lower and upper phases, respectively. The upper phase is assumed to be incompressible to avoid interface displacement during compression. In these conditions, Eq. ( 38) can be written as follows.

x∇ • U (1) + (1 -x)∇ • U (2) = - V 0 L (t -t i ) ( 38
)
From the general equation of the interface force equilibrium :

k τ (2) ik -τ (1) ik = (p (2) -p (1) ) n i (39)
we can write Eq. ( 40) by substituting for p (i) and τ

(i)
yy the terms in Eq. ( 41) for i = y (see Eq. ( 28) for the χ T dependence with Lamé's coefficients).

(2µ

(2) E + λ (2) E )∇ • U (2) = (2µ (1) E + λ (1) E )∇ • U (1) (40)
Combining Eqs. [START_REF] Amiroudine | A lagrangian-eulerian compressible model for the trans-critical path of near-critical fluids[END_REF] and Eq. ( 40), ∇ • U (1) and ∇ • U (2) can be expressed at time t as a function of the constants in question,

∇ • U (2) = - (2µ (1)
E + λ (1) E ) V 0 L (t -t i ) x(2µ (2)
E + λ (2) E ) + (1 -x)(2µ (1)
E + λ (1) E ) ∇ • U (1) = - (2µ (2)
E + λ (2) E ) V 0 L (t -t i )
x(2µ

(2) E + λ (2) E ) + (1 -x)(2µ (1)
E + λ (1)
E ) giving the analytical solution for pressure, normalized shear stress and the τ yy -element of the shear stress tensor through the following relations with i = 1, 2.

p (i) = - 1 χ (i) T ∇ • U (i) τ (i) = 2 µ (i) E ∇ • U (i) τ (i) yy = 4 3 µ (i) E ∇ • U (i) (41)

Numerical simulation

The 3D simulations of the two-material compression were carried out on a 4 x 65 x 4 mesh grid considering cork as the lower material and rubber as the upper one (Fig. 10).

FIG. 10: System geometry and boundary conditions for the two-elastic solid compression.

Equal thicknesses for the two phases is considered (x = 0.5), as well as a constant injection speed V 0 = 0.0001 m/s. Similar boundary conditions as the one-phase compression studies are considered. The values of the cork constants µ

E , λ

E , χ (1)
T , and rubber constants µ

E , λ

E , χ (2)
T are indicated in table I. These constants are assigned to the discretisation nodes using the phase function C, assuming discontinuity across the interface. The time increment, t = 1 s is more than three orders of magnitude greater than the characteristic time L /2 c T of the cork phase (that of rubber being lower) yielding a linear steady compression within the two phases. The pressure and normalized shear stress profiles are reported in Figs. These profiles clearly show the uniformity of pressure and normalized shear stress throughout each material, the exact values being obtained in reference to the low relative error (table V). As expected, pressure and shear stress discontinuities are obtained at the interface in accordance with the interface equilibrium condition (Eq. ( 39) with i = y). Note that force equilibrium at the interface is not enforced as a boundary condition, but is an implicit so-lution of the unified equation. This force equilibrium criterion is shown to be obtained with good accuracy (see last column of table V).

Material p -p i p -p th p th -p i τ -τi τ -τth τth -τi τ yx -τ yxi τ (2) yy -τ (1)
yy -(p (2) -p (1) )n y

(Pa) (Pa) (Pa) (Pa)
cork (1) 33234 3 × 10 -13 99701 8 × 10 -14 -66467 1.8 × 10 -7 rubber (2) 99674 1 × 10 -12 39.87 1 × 10 -10 -27

TABLE V: Calculated pressure and normalized shear stress within the two materials at the final compression time t f = 100 s; corresponding relative errors and interface equilibrium criteria obtained for time increment t = 1 s.

F. Complex problem of fluid-elastic solid interaction

Problem description

Here, the unifying model for fluid flow and elastic solid deformation is put to test through the study of a complex fluid-solid interaction. This problem is based on previous 2D studies [START_REF] Breuer | Fluid-structure interaction using a partitioned semi-implicit predictor-corrector coupling scheme for the application of large-eddy simulation[END_REF][START_REF] Degroote | Partitioned simulation of the interaction between an elastic structure and free surface flow[END_REF][START_REF] Turek | Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow[END_REF] of laminar incompressible flow around an elastic object resulting, as a function of parameter settings, in either a steady state solution or self-induced oscillations of the object. The 2D configuration of the system is addressed here in 3D by considering a system thickness of 0.01 m. The system setup (Fig. 13) is a long channel of length L = 2. A parabolic profile is set at the left channel inflow,

V 0 (0, Y, Z) = 1.5 V 0 Y (H -Y )   H 2   2 = 6.0 0.1681 Y (0.41 -Y ) (42)
where V 0 is the mean inflow velocity. Smooth increase of velocity profile in time is used,

V 0 (t, 0, Y, Z) =                  V 0 (0, Y, Z) 1 -cos( π 2 t) 2 if t < 2.0 s V 0 (0, Y, Z) otherwise (43)
where V 0 (0, Y, Z) is the velocity profile given in Eq. ( 42).

The Neumann condition, dV /dX = 0, is prescribed for the outflow, and the no-slip condition V = 0 is prescribed for the fluid on the top and bottom walls. A penalty of the velocity field, V = 0, is set within the entire obstacle [START_REF] Khadra | Fictitious domain approach for numerical modelling of navier-stokes equations[END_REF]. The physical constants of the fluid and elastic object and the two studied inflow mean velocities are indicated in table VI. To sum up, similar settings are considered compared to the 2D benchmarking study [START_REF] Turek | Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow[END_REF],

except for i) the compressibility characteristics of the fluid which is perfectly incompressible in the benchmarking study whereas the physical coefficient of compressibility of water is considered in this study; ii) the no-slip boundary condition at the fluid/elastic solid interface, since the unified model formulation does not require any explicit boundary condition at the fluid/solid interface.

Numerical results

The movement of the elastic flag during channel flow is analyzed through the Xdisplacement U X and the Y-displacement U Y of the A (0.6 m; .0.2 m)-point at its very right end. Fig. 14 shows these U X, U Y -displacements for the parameter setting of case 1 (Table VI). After early time oscillations, the steady state evolution is obtained. It is characterized by A-point coordinates very close to those calculated in the benchmark study [START_REF] Turek | Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow[END_REF] as shown in Table VII. [START_REF] Turek | Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow[END_REF] .

For the parameter setting of case 2 (Table VI), a non-steady state evolution was obtained with periodic oscillations of the flag as shown in Fig. [START_REF] Banks | Deforming composite grids for solving fluid structure problems[END_REF] with the Y-displacement curve of the A-point. The frequency, the mean value (average of the minimum and maximum values) and the amplitude (difference of the minimum or maximum values from the mean) of the periodic evolution were calculated (Table VIII). These period characteristics are close to those computed in the benchmarking study (also see Table VIII). In this part, we study thermal transport within a monophase fluid system, and its coupling with mechanics through the thermal expansion of the fluid. So, the complete governing equation system [START_REF] Breuer | Fluid-structure interaction using a partitioned semi-implicit predictor-corrector coupling scheme for the application of large-eddy simulation[END_REF] in which the residual shear stress τ 0 equals zero has to be solved.

This apparent simple problem nevertheless involves different time scale phenomena such as: i) phonon dynamics with a very short time constant τ p ≈ 10 -12 s, ii) thermoacoustic propagation time constant τ c = L /c L , where L is the characteristic length of the system and c L the sound celerity, ii) the fluid flow time constant τ n = L/V 0 where V 0 is the characteristic velocity of the fluid flow, and iv) the thermal diffusion time constant τ d = ρ c p L 2 /k. We studied two time scale thermal transport characteristics, those of thermoacoustics and of thermal diffusion, so that the heat equation of system ( 16) can be reduced to Fourier's equation.

Φ + k ∇ T 0 -dt β T 0 ρ 0 c v χ T ∇ • V -dt 1 ρ 0 c v (∇ • Φ -q -φ) = 0 (44)
A. Numerical scheme and flux equation discretization

Lagrangian heat fluxes Φ are determined from the vectorial solving process of Eq. [START_REF] Osher | Level set methods: An overview and some recent results[END_REF] which results in the time dicretization of Eq. ( 44) using the known temperature and velocities at time step n.

Φ + k ∇ T n -dt β n T n ρ n c n v χ n T ∇ • V n -dt 1 ρ n c n v (∇ • Φ -φ n ) = 0 ( 45
)
Here φ n is the dissipation function at time step n defined from velocity vector V n using Eq. ( 18)

Then Lagrangian velocities V at time step n + 1 are determined from V n and Φ using Eq. ( 22), where the predictive Lagrangian pressure term,

  p n -∆t 1 χ n T ∇ • V  
is extended to the following ones.

p n -dt 1 χ n T + β n 2 T n ρ n c n v χ n 2 T ∇ • V -dt β ρ n c n v χ n T (∇ • Φ -φ n )
The Lagrangian pressure p, temperature T and density ρ at time step n + 1 are then re-actualized from the Lagrangian velocity and flux vectors using the following equations.

                         ρ = ρ n e -dt ∇•V T = T n -dt β n T n ρ n c n v χ n T ∇ • V -dt 1 ρ n c n v (∇ • Φ -φ n ) p = p n -dt 1 χ n T + β n 2 T n ρ n c n v χ n 2 T ∇ • V -dt β n ρ n c n v χ n T (∇ • Φ -φ n ) (46)
The Lagrangian dissipation function [START_REF] Degroote | Partitioned simulation of the interaction between an elastic structure and free surface flow[END_REF] φ was calculated from the Lagrangian velocity vector V and the thermodynamic coefficients at time step n.

The Eulerian variables p n+1 , T n+1 , ρ n+1 , φ n+1 , V n+1 and Φ n+1 at time step n + 1 were determined from the calculated Lagrangian variables p, T , ρ, φ, V and Φ following the same advection scheme described in section III-A. Similarly, the thermodynamic coefficients The simulation was carried out on a one-dimensional 1024 mesh grid. Time increment ∆t = 10 -5 s was chosen to be two orders of magnitude lower than the thermoacoustic characteristic time τ c = L/c L ≈ 3. 10 -3 s where c L = 1/ √ ρ i χ T i = 430 m.s -1 is the velocity of elastic wave propagation at initial time. The sudden heating of the air in contact with the wall produces thermoacoustic wavelets in agreement with previous studies [START_REF] Farouk | Numerical study of thermoacoustic waves in an enclosure[END_REF][START_REF] Aktas | A numerical study of the generation and propagation of thermoacoustic waves in water[END_REF]. These wavelets propagate toward each other in the direction of the cavity center, as the temperature profile at time t = 10 -3 s shows (Fig. 17). At that time, the maximum temperature amplitude of the wave is very low, not exceeding 0.04 K around the initial temperature T i = 400 K, whereas the air temperature near the wall reaches 400 K. As regards the corresponding pressure profile in Fig. 18, wavelets are also observed with a strict similar wavelength and amplitude decay. No phase difference is observed. Strong coupling between temperature and pressure thus occurs for thermoacoustic wavelets. Temperature depends strongly on air expansion or compression since low temperature domains correspond strictly to tensile pressure domains and conversely high temperature domains correspond to compression domains (comparison between Fig. 17 and Fig. 18). This is not the case far from the wavelets and especially near the walls where pressure is close to its initial value whereas temperature is very high. In this latter case, it is the density decrease which compensates for the temperature increase.

Finally, from the position of the first peak base, x p = 0.34 m, at time t = 10 -3 s (Figs. 17

Thermal diffusion time scale solution

The simulation was carried out on a one dimensional 1024 mesh grid. Time increment ∆t = 0.1 s was chosen to be thirty times greater than the thermoacoustic characteristic time Theoretical solution for final steady state is given by the thermodynamic equilibrium state.

τ c = L/c L ≈ 3. × 10 -
It gives uniform air density, temperature and pressure within the cavity, that is : ρ th f = ρ th i , T th f = 400 K and p th f = ρ th f r T th f = 135100 Pa (r = 287 J.kg -1 K -1 for air). Mass conservation is maintained to within computer accuracy. The calculated pressure p f = 135076 Pa which is not deduced from the state equation, is obtained with a relative error of 1. In the framework of this complete thermo-mechanical formulation, we first focus in this paper on the consistency of the unified governing equation for compressible fluid flow and elastic solid deformation by considering isothermal problems. First of all, the same equation is shown to deal with multi-time scale problems with the study of elastic solid deformation. Indeed, for time scale analysis lower than the characteristic time of acoustic wave propagation, we succeed in following the longitudinal and transverse elastic waves with an accurate wave velocity, whereas for longer analysis time, the compression process was obtained with steady evolution. For two-phase systems, the model is shown to implicitly solve the interface equilibrium criterion accurately, even in the case of the fluid-structure interaction. There is no longer a need for interface boundary conditions inherent to governing equations differing in type across the interface. Note that "incompressible" fluid flow can be modeled with the formulation by simply considering the experimental compressibility coefficient of the phase involved along with time scale analysis longer than the characteristic time of acoustic wave propagation.

Similarly, through the study of air heating in an enclosure, the complete thermomechanical formulation was shown to succeed in catching thermoacoustic waves for analysis time scales shorter than the characteristic time of acoustic wave propagation. For longer time scale analysis, steady evolution of density, pressure and temperature is obtained until the thermodynamic equilibrium state is reached. This new mathematical formulation of the thermo-mechanical equations including thermodynamics brings major assets for the modeling of multi-phase and multi-physics problems. Currently limited to isotropic linear elastic solids, we aim to extent this formulation to anisotropic elasticity and inelasticity such as the elasto-plastic behavior. Finally, if the purpose of this paper was not to propose the more efficient numerical scheme, this is of course of great importance and we will have to think about further developments in order to deal with more complex fluid-structure problems that we can found in the literature [START_REF] Kassiotis | Nonlinear fluid-structure interaction problem. part i: implicit partitioned algorithm, nonlinear stability proof and validation examples[END_REF][START_REF] Bathe | Finite element developments for general fluid flows with structural interactions[END_REF][START_REF] Bathe | A mesh adaptivity procedure for cfd and fluid-structure interactions[END_REF][START_REF] Kassiotis | Nonlinear fluid-structure interaction problem. part ii: space discretization, implementation aspects, nested parallelization and application examples[END_REF].

Figure captions

FIG. 1 :

 1 FIG.1: Evolution of the Lagrangian thermo-mechanical variables for a material point M between time t 0 to time t = t 0 + dt

FIG. 2 :

 2 FIG. 2: System geometry and boundary conditions for the compression of fluid and solids (a) and the shear deformation of solids (b).

FIG. 3 :

 3 FIG. 3: Comparison between simulated and analytical exponential evolutions of relative density of air, cork and rubber during compression (time increment t = 1s); the dashed line corresponds to a linear evolution.

FIG. 4 :FIG. 5 :

 45 FIG. 4: Comparison between simulated and analytical exponential evolutions of pressure during air compression (time increment t = 1s); the dashed line corresponds to a linear evolution.

12 rubber 1 . 67 ×

 12167 Material p -p i (Pa) p -p th p th -p i τ -τi (Pa) τ -τth τth -τi cork 1.67×10 5 4 × 10 -12 5.00 × 10 5 4 × 10 -10 8 2 × 10 -13 6.67 × 10 4 2 × 10 -13

FIG. 7 :

 7 FIG.7: Comparison between the simulated and the analytical linear evolutions of the normalized shear stress during simple shear deformation of cork (circles) and rubber (squares) for time increment t = 1 s.

  for c T -values of cork and rubber).

11 and 12 29 FIG. 11 :FIG. 12 :

 1112291112 FIG. 11: Pressure profile along the (X=0.5 m;Y; Z=0.5 m)-line at the final compression time t f = 100 s.

2 FIG. 13 :

 213 FIG.13: System geometry for the study of fluid-elastic solid interaction.

  Compressibility coefficient, χ T (Pa -1 ) 4.5 × 10 -10 4.5 × 10 -10 Mean inflow velocity, V 0 (m.s -1 ) 0.2 1 TABLE VI: Physical constants of the fluid and elastic flag involved in the fluid-structure interaction problems. The parameter settings lead to Reynolds number values (Re = 2 ρ V 0 r µ ) low enough to induce laminar flow (Re = 20, 100, for the first and the second parameter settings, respectively). Note that the time increment t = 0.001 s is not sufficiently small to follow the longitudinal elastic wave propagation within the fluid since the characteristic times L /c L = 1.7 10 -3 s and H /2 c L = 1.4 10 -1 s (with c L = 1 /ρχ T ) are of the same order of magnitude or greater. So a very small effect of fluid compressibility will have to be expected. The calculation was performed on a 415 x 110 x 2 mesh grid, with refinement in the vicinity of the elastic flag and the obstacle. The material constants are assigned to the system discretization nodes using the phase function C and assuming the continuity of the constants across the interface (see paragraph 4).

FIG. 14 :

 14 FIG. 14: U X, U Y -displacements of the A-point of the elastic flag during channel flow for parameter setting of case 1 (TableVI).

FIG. 15 :

 15 FIG. 15: Y-displacement of the A-point of the elastic flag during the channel flow for parameter setting of case 2 (TableVI).

Fig. ( 16 )

 16 Fig.[START_REF] Breuer | Fluid-structure interaction using a partitioned semi-implicit predictor-corrector coupling scheme for the application of large-eddy simulation[END_REF] shows a snapshot of the flag morphology during the self-induced oscillations. As expected, compressive values are calculated under the convex interface while tensile values are calculated under the concave interface. In the fluid, as in previous study[START_REF] Degroote | Partitioned simulation of the interaction between an elastic structure and free surface flow[END_REF], compressive areas with comparable intensities are localized in the left part of the obstacle and above the flag when going down, and tensile areas are localized on the right part of the obstacle and below the flag. Note that the end part of the flag becomes slightly rounded. This non physical deformation is the consequence of the numerical VOF representation of the interface which tends to smooth interfaces.

  (χ n+1 T , β n+1 , c n+1 v , µ n+1 , k n+1 ) were updated from the phase function C and the values of the intensive variables T n+1 and p n+1 . For a perfect gas we have χ n+1 T = 1/p n+1 and β n+1 = 1/T n+1 . B. Air cavity heating 1. System description A 2D square cavity of width L = 1 m, containing air at pressure p i = 101325 Pa, is heated from initial temperature T i = 300 K to T f = 400 K. The horizontal walls of the cavity are assumed to be adiabatic and the vertical walls are kept isothermal at temperature T f . The initial air density ρ i , initial isothermal compressibility coefficients χ T i and viscosity coefficients µ are indicated in Table I. Such thermo-elastic coupling also requires knowledge of the initial coefficient of thermal expansion β i = 3.33 × 10 -3 K -1 , the specific heat c v = 719 J.kg -1 .K -1 and the thermal conductivity k = 2.63 × 10 -2 W.m -1 .K -1 . The µ, c v , kcoefficients are assumed to be constant within the temperature and pressure intervals studied, unlike the χ T and β coefficients which are pressure and temperature dependent. Assuming air to be a perfect gas, we have: β = 1/T and χ T = 1/p. The value of thermoacoustic characteristic time τ c ≈ 3. × 10 -3 s and that of the thermal diffusion τ d ≈ 4.5 × 10 4 s 2. Thermoacoustic time scale solution

and 18 )FIG. 17 :FIG. 18 :

 181718 FIG. 17: Temperature profile T -T i at t = 10 -3 s along the line perpendicular to the vertical heating wall of the cavity set at 400 K; time increment ∆t = 10 -5 s.

  3 s in order to obtain pressure convergence at each time step. As expected, uniform pressure is calculated throughout the cavity which increases with time until the equilibrium state is reached, Fig.19. The different features for the evolution of relative temperature and pressure at the center point of the cavity, i.e. sharp increase for pressure and delayed increase for temperature clearly shows that, early on, temperature variations depend on local temperature gradient whereas pressure is consistent with the global mean temperature value.

  FIG. 19: Evolutions of relative pressure p * = (p -p i )/(p f -p i ) and temperature T * = (T -T i )/(T f -T i ) at the cavity center for large time scale increment ∆t = 0.1 s

Fig. 1 :

 1 Fig.1: Evolution of the Lagrangian thermo-mechanical variables for a material point M between time t 0 to time t = t 0 + dt.

Fig. 2 :

 2 Fig. 2: System geometry and boundary conditions for the compression of fluid and solids (a) and the shear deformation of solids (b).

Fig. 3 :

 3 Fig. 3: Comparison between simulated and analytical exponential evolutions of relative density of air, cork and rubber during compression (time increment t = 1s); the dashed line corresponds to a linear evolution.

Fig. 4 :

 4 Fig. 4: Comparison between simulated and analytical exponential evolutions of pressure during air compression (time increment t = 1s); the dashed line corresponds to a linear evolution.

Fig. 5 :

 5 Fig. 5: Comparison between simulated and analytical linear evolutions of pressure and normalized shear stress during cork compression (time increment t = 1s).

Fig. 6 :

 6 Fig.6: Second order convergence rate of numerical air pressure towards the analytical solution with time increment.

Fig. 7 :

 7 Fig. 7: Comparison between the simulated and the analytical linear evolutions of the normalized shear stress during simple shear deformation of cork (circles) and rubber (squares) for time increment t = 1 s.

Fig. 8 :

 8 Fig.8: Pressure variation at the center point of the rubber system during longitudinal elastic wave propagation; time increment, t = 3 10 -5 s.

Fig. 9 :

 9 Fig.9: Variation, during transverse wave propagation, of the τ xy -element of the shear stress tensor at the center point of the rubber system; time increment, t = 2 10 -3 s.

Fig. 10 :

 10 Fig. 10: System geometry and boundary conditions for the two-elastic solid compression.

Fig. 11 :

 11 Fig. 11: Pressure profile along the (X=0.5 m;Y; Z=0.5 m)-line at the final compression time t f = 100 s.

Fig. 12 :

 12 Fig. 12: Normalized shear stress profile along the (X=0.5 m;Y; Z=0.5 m)-line at the final compression time t f = 100 s.

Fig. 13 :

 13 Fig. 13: System geometry for the study of fluid-elastic solid interaction.

Fig. 14 :

 14 Fig. 14: U X, U Y -displacements of the A-point of the elastic flag during channel flow for parameter setting of case 1 (TableVI).

Fig. 15 :

 15 Fig. 15: Y-displacement of the A-point of the elastic flag during the channel flow for parameter setting of case 2 (TableVI).

Fig. 16 :

 16 Fig. 16: Pressure map in the surroundings of the flag at the simulation time t = 17.1 s of

Fig. 17 :

 17 Fig.17: Temperature profile T -T i at t = 10 -3 s along the line perpendicular to the vertical heating wall of the cavity set at 400 K; time increment ∆t = 10 -5 s.

Fig. 18 :

 18 Fig.18: Pressure profile p -p i at t = 10 -3 s along the lines perpendicular to the vertical heating wall of the cavity; time increment ∆t = 10 -5 s.

TABLE II :

 II Calculated pressure and normalized shear stress within the elastic solids at the final compression time t

f = 100 s; corresponding relative errors obtained for time increment t = 1 s.

TABLE IV :

 IV Comparison between theoretical and calculated wave propagation velocities for the two basic cases of longitudinal and transverse elastic waves.

	E. Two-phase compression
	1. Analytical solution

  Table VI).

	Case 1	U X (m)	U Y (m)
	benchmark	2.27 × 10 -5 8.21 × 10 -4
	present result 2.26 × 10 -5 8.29 × 10 -4

TABLE VII :

 VII Comparison of the U X and U Y displacements of the A-point during the steadystate evolution with that of the benchmarking study

  Table VI).

	Case 2	frequency mean value amplitude
	benchmark	2.0	1.23 × 10 -3 8.06 × 10 -2
	present result	1.9	1.1 × 10 -3 7.6 × 10 -2

TABLE VIII :

 VIII Characteristics of the Y-periodic displacement of the A-point during the non-

	steady state evolution; comparison with the same characteristics computed in the bench-
	marking study [25]

Acknowledgments

The authors thank the Aquitaine Regional Council for the financial support dedicated to a 256-processor cluster investment, located in the I2M-TREFLE laboratory.


Lucile Cornu

Manuel Gaudon

Evgeniy Ilin

Cyril Aymonier

Philippe Veber

Alain Garcia

Myrtil Kahn

Yohan Champouret

Veronique Jubera
email: jubera@icmcb-bordeaux.cnrs.fr

Luminescence of sensitive materials: towards new optical sensing

Keywords: Sensitive materials, luminescence, X-chrome, photoactive materials, sensors, data storage

In the last decades, considerable efforts have been carried out to develop new tools and knowledge in the domain of functionalized materials, for application ranging information, lighting, communication, energy, optical sources or detection. To couple an optical answer with another property is possible through the follow-up of the luminescence. The control of the structural symmetry, oxidation state or surface chemistry enables the chemists to precisely tune the emission. Illustrations will be provided on inorganic powder and crystal materials.

1-INTRODUCTION

Functionalized materials have been drawing significant interest in the domain of application ranging information, lighting, communication, energy, optical sources or detection [START_REF] Ronda | Luminescence, from theory to applications[END_REF]. The control of the structural symmetry, oxidation state or surface chemistry enables the chemists to precisely tune the emission. For instance: (i) A red to green luminescence switch occurres for manganese doped spinel treated between 1200°C and 1350°C. By modifying the stoichiometry, it is possible to change the temperature range of this luminescence spectrum switch [START_REF] Cornu | [END_REF]3]. (ii) UV irradiation of cerium doped indium elpasolite matrix leads to the decrease of the doping element emission in favor of a bright orange emission which results from the reduction of indium. The redox process is reversible under adequate irradiation or heat treatment temperature. The sensitized materials are stable over two decades, showing hence perfectly bi-stable state [4,5]. (iii) Emission of ZnO is strongly related to the synthesis route. Low operating temperature can result in the appearing of surface and bulk defects. Their existence impacts not only the emission wavelength but also its stability under specific atmosphere [6][7][8][9]. Tunable emission properties of these materials make them potential candidates for highly sensitive thermal, gas or photosensors.

2-EXPERIMENTAL RESULTS AND DISCUSSION

Color switch due to luminescent cation migration [2,3]

ZnAl 2 O 4 spinel compound is an interesting material with a specific cubic crystallographic structure which possesses 64 tetrahedral sites and 32 octahedral sites. A doping with luminescent metallic cations such as transition metals or rare earth elements can be used to gain satisfying optical properties for flat panel displays [10] or as pigments in glazes or stonewares [11]. Indeed, one can modulate the location of the doping element by changing the experimental condition. The proportion of divalent cations in octahedral site is called inversion rate (and can vary between 0 and 1. A direct spinel structure is constituted by only 8 of the tetrahedral sites occupied by the Zn and 16 of the octahedral sites occupied by Al cations and the oxygen anions fully occupy the 32e positions [12]. In ZnAl 2 O 4 , a part of Zn 2+ cations is located in the 16d sites and a part of Al 3+ cations is located on the 8a sites. In literature the  parameter ranges between 0 <  < 0.05 [13][14][15] for the ZnAl 2 O 4 matrix.

Because it has been demonstrated that the inversion rate is strongly linked to the experimental conditions [START_REF] Lenestour | Corrélation structure-propriétés d'absorption UV-Vis-IR associée aux états de valence du cuivre dans des oxydes à base de zinc de type spinelle et würtzite[END_REF][START_REF] Matsui | [END_REF][18][19][20], we propose here to use a Pechini synthesis to have access to a large crystallite size distributio n to impact the  ratio [3]. To control the inversion factor of doped ZnAl 2 O 4 materials may be of great interest especially when the luminescent properties are affected by the location of the emitting center, the crystal field strength . By consequence, divalent manganese was selected because of its appropriate optical transition (3d-3d) and its appropriate oxidation state and radius size.

The solubility limit of divalent manganese was controlled by XRD pattern and TEM -EDX analyses. Even if the 30% doped material clearly shows a secondary aluminum-rich phase (-Al 2 O 3 phase n°00-046-1212), XRD pattern and TEM -EDX show a pure spinel phase with experimental Zn:Mn:Al concentrations close to the theoretical ones for the 15% doped-samples. Considering these results, manganese concentration was fixed at 0.5% to avoid precipitation of the alumina phase and potential concentration effect which may quench the manganese luminescence. Since it has not been possible to discriminate Zn from Mn ions because of the low doping rate and their vinicity in the periodic None significant variation of the a cell parameter was observed as function of the temperature but a significant increase of the crystallite size is reported. The crystallites grow from about 20 nm for the lowest temperature up to 120 nm for the highest one with a decrease of the global inversion rate from 16% to 2.6%. One can deduce that the annealing temperature impacts the cationic distribution within the spinel structure. This global inversion rate taking into account the manganese distribution, the influence of the doping element should be visible on the optical properties.

The emission properties of Mn-doped ZnAl 2 O 4 depend on the symmetry of the coordination polyhedron within the spinel network in which the manganese ion is introduced. A location in a tetrahedral environment is (b)

generally associated to a green luminescence of this element. A higher crystal field effect as in octahedral site shifts the emission to lower energy. Moreover, a blue shift of the metallic zinc-oxygen charge transfer band is detected between 1200°C and 1350°C on the diffuse reflectance spectra, this last thermal treatment leading to a total migration of the divalent cation into the fourfold coordination site [3]. The optical properties are illustrated for these two temperatures in Figure 2. The 1350°C treatment makes possible the obtaining of a total direct Mn-doped spinel i.e.,with Mn ions fully located in the tetrahedral sites. Many papers have reported the green emission of this transition metal in a tetrahedral environment [21][22][23]. The main 4 T 1 → 6 A 1 transition is observed at 508 nm. In addition vibronic vibrations are detected at 514 nm, 520 nm and 529 nm, respectively. For an excitation at 426 nm, red emission peaks ranging between 640 nm and 730 nm were attributed to trivalent chromium impurities coming from the raw reactant. The 1200°C sample possesses an inversion rate equal to about 5%. By consequence, a significant amount of manganese is expected to be in the octahedral cationic sites. The spectrum obtained under a 426 n m excitation shows a shift of the emission to higher wavelengths. A large band centered at 738 nm is clearly detected at room temperature. The comparison with our previous study on the non-doped matrix confirms the fact that this red emission is related to the location of manganese in the partially indirect spinel host lattice [3].

For this annealing temperature, the luminescence could be also linked to the occurrence of punctual defects such as the presence of oxygen vacancies in the vicinity of Zn Al antisite. But defect emission is not observed anymore under a 426 nm excitation wavelength. A part of manganese ions could be localized either in distorted five -fold coordinated polyhedra or in regular AlO 6 octahedra, substituting for Al 3+ . One can conclude that the Pechini synthesis has successfully stabilized divalent manganese ions in the Al 3+ site with coordination number higher than 4 which makes possible the observation of a red emission for the Mn-doped ZnAl 2 O 4 spinels. In addition one can say that the large component of the excitation spectrum indicates that the polyhedron of the doping element gets a well-defined and regular geo metry. Due to thermal treatments, divalent cations including the doping element, progressively migrate from the sixfold-coordinated sites to the fourfold-coordinated ones. It could be concluded that increasing the heat treatment temperature, the modification of the emitting properties is clearly visible as the emission switch es from the red to the green range.

Color switch due to redox phenomenon [5]

In the A 2 BMX 6 elpasolite phases, A and B belong to the alkali group, M is a transition metal or a rare earth element, and X is an halogen [24][25][START_REF] Tanner | Transition Metal and Rare Earth Compounds III)[END_REF]. Doped with active cation, these materials are studied for various applications as light-emitting diodes, laser devices, [START_REF] Pedro | [END_REF][28][29] scintillators and detectors, [4,[30][31] and for data storage devices. The elpasolite material crystallizes in the cubic system (Fm-3m space group) and can be described as a double perovskite with cationic ordering. The A cation possesses a twelvefold-coordination polyhedron (central cuboctahedral site, 8c Wyckoff position) whereas the B alkali cation and the M transition metal/rare earth element regularly alternate along the unit cell axis in a sixfoldcoordinated site (octahedral site).

Concerning the Ce-doped Rb 2 KInF 6 elpasolite compound, the spectral distributions reflected the substitution of In 3+ by Ce 3+ in the regular octahedra but also the fact that a minor part of Ce 3+ ions is located in the cuboctahedron. [5] Under UV irradiation, a redox phenomenon between the Ce/In ions occurs leading to a reversible luminescence switch based on the following equilibrium: 2Ce 3+ +In 3+ ⇄2Ce 4+ +In + . Such a phenomenon is rare and its origin is due to the peculiar crystal structure and composition of the studied material.

The understanding of the luminescent properties requires the doping element cation location. For such a purpose, DFT calculations were performed (VASP code). Different structural models were considered for the Ce 3+ to In 3+ substitution (calculation, within a 2a2bc 80-atoms supercell) in order to evaluate what is the more appropriate substitution site for trivalent cerium. The converged data confirm that cerium cations are preferentially located on the transition metal site (indium site) but this creates an oxidation stress on these latter ions with an equilibrium valence equal to 3.36. The rubidium site represents a second location of cerium element. The corresponding coordination polyhedron relaxes to form an eightfold-coordinated site but the resulting valence is equal to 2.63, i.e. in this site, the Ce 3+ ions are submitted to reductive stress. Luminescent spectroscopy was performed on a crushed single crystal grown under controlled atmosphere by the Bridgman-Stockbarger method [5]. As expected the main emission observed under a 315 nm source is associated with the location of Ce 3+ in In 3+ octahedral site as illustrated on Figure 3a. The most intense band peaks at 480 nm. The complete excitation spectrum is composed of different bands due to the splitting of the 5d orbitals under the crystal field effect. The calculated Stokes shift is equal to 10,714 cm -1 . In this matrix, the trivalent cerium emission is significantly red-shifted compared with previous literature studies [32,33] and the Stokes shift is more pronounced. A second weaker signal located at 340 nm is also detected in the Ce-doped Rb 2 KInF 6 materials (Figure 3b) under a 273 nm excitation wavelength. The corresponding Stokes shift is equal to 7630 cm -1 . We can discard the existence of a secondary phase contribution because of the XRD control performed on a crushed asgrown crystal which showed a pure phase. The final hypothesis is the location of trivalent cerium in a second environment. This effect has already been observed by B. F. Aull et al. in rare earth-doped Rb 2 NaYF 6 materials; the authors announced a larger site with a higher coordination [32]. We can conclude that this secondary emission corresponds to the contribution of trivalent cerium in the rubidium cuboctahedral site. Under a 315 nm, the main blue -green emission decreases in favor of a new orange luminescence at 650 nm. This new emission band was attributed by Chaminade et al. and M.A. Buñuel et al. to the 5s 1 5p 1 5s² radiative transition of monovalent indium, the UV light making possible the displacement of the following equilibrium: 2Ce 3+ +In 3+ ⇄2Ce 4+ +In + [4,34]. This phenomenon is reversible and an excitation in the monovalent indium absorption band at 255 nm results in an increase of the cerium contribution at 480 nm. Under appropriate UV irradiation, the 4f cerium electrons are promoted to the 5d orbitals. While electrons are in the lowest cerium 5dT 2g orbitals (5d xy , 5d xz and 5d yz ), they can be transferred to the neighboring indium orbitals which point in the adequate direction between the fluorine anions. For reminding, the DFT calculations showed that the rare earth ions localized in the indium sites are subject to oxidative stress (valence >3). The resulting reduction of indium at the monovalent oxidation state gives rise to an ns 2 monovalent cation with a highly polarizable electronic doublet. By consequence, a back electron transfer is highly probable. In this elpasolite matrix which possesses a redox metallic ions pair, luminescence can be monitored under appropriate UV excitation wavelength because of the succession of K-F-In/Ce along quasi-collinear axes. A new photochromic inorganic oxide is so obtained.

Luminescence and synthetic route [6,9]

ZnO is a II-IV wide gap semiconductor material with a large band gap energy of 3.37 eV at room temperature. Typical photoluminescence spectra exhibit two types of emission strongly related to the synthetic route. The UV excitonic recombination is located in the UV part [35][36][37][38][39] whereas the contribution of bulk and surface defects result in large emission band in the visible range [40][41][42][43]. The ratio between both contributions depends on the crystallinity rate and the capability to control surface states.

We studied optical properties of ZnO nanoparticles obtained through organometallic route [8,9] and supercritical fluid route [6,7] as a function of their size and shape. The particle shape control is developed by using different ligands (organometallic route) or precursor concentration and operating parameters as time residence (supercritical fluids route).

As an example, SEM pictures of spherical nanoparticles are reported on Figure 4. The luminescence properties of these materials were investigated. The corresponding emission graphs are reported on Figure 5. The luminescence of the nanomaterials obtained through organometallic synthesis present a visible emission. The recombination of the photogenerated electron-hole pairs leads to a wide band which is maximum at about 600 nm. The visible luminescence of ZnO is still subject of controversies. This "yellow orange" emission seems to be linked to the recombination of a shallowly trapped electron with a deeply trapped hole. In this process the photogenerated hole is first trapped by surface defects at the surface of the material. Then it migrates to vacancy levels located in the core of the particle leading to the stabilization of a deep hole trapped level. The emission takes place when the photogenerated hole trapped in the deep defect as an oxygen vacancy recombines with the photogenerated electron trapped in a shallow level located below the conduction band. The UV excitonic is generally observed for high temperature gas phase routes because these syntheses allow the generation of high crystalline quality ZnO materials. But the control of the shape and the size of the nanostructures is difficult. The supercritical fluids route makes possible the precipitation of nano-size particles. By controlling the pressure and temperature, it is possible to avoid the stabilization of bulk or surface defects. The lack of trapping levels within the gap makes possible a direct radiative recombination of the photogenerated electron hole pair. By consequence, the resulting emission spectrum is a unique sharp peak peaking in the UV range. Such optical properties announce significant difference of sensitivity regarding the environment of the nanoparticles and their surface activity.

CONCLUSION

As demonstrated, photoluminescence can provide a diagnostic on a material. One can play with the temperature to make a doping element migrate within a spinel structure. An evaluation of the thermal history of the compounds can be done. It is also possible to switch the response of a luminescent elpasolite under an appropriate UV excitation beam. Application such as UV detectors or data storages can be proposed by the monitoring of the emission. Finally, once the synthesis route is well-controlled, one can choose to favor UV or visible radiative de-excitation of the well-known ZnO matrix.  exc = 325nm a)

Figure 1 -

 1 Figure 1-a) X-ray diffraction patterns of the 1000 °C heat-treated Mn doped samples (Red arrows show the secondary phase: α-Al 2 O 3 corundum) and b) EDX-c he mical maps of the 15% Mn-content compound (corresponding Zn, Al and Mn molar percentages are reported in the figure)

Figure 2 .

 2 Figure 2. Room te mperature excitation and e mission spectra of 0.5% Mn-doped ZnAl 2 O 4 compound obtained after thermal treatment at (a) 1200 °C and (b) 1350 °C.

Figure 3 .

 3 Figure 3. Room temperature excitation and emission spectra of Ce 3+ doped Rb 2 KInF 6 matrix before irradiation, for a) Ce 3+ in indium site, and b) Ce 3+ in rubidium site.

Figure 4 .

 4 Figure 4. Spherical nanoparticles of ZnO obtained through organometallic synthetic route

Centered at 5 Figure 5 .

 55 Figure 5. Emission spectra of ZnO nanoparticles obtained through supercritical fluid route a) and organometallic synthetic route b) under a 325 nm excitation wavelength.

  table, Rietveld refinements were performed to determine the global divalent cation inversion rate (+') on the 0.5% Mn-doped ZnAl 2 O 4 compounds. Materials were heat-treated at various annealing temperatures ranging fro m 500°C to 1350°C. Results obtained considering the (Zn 0.995- Mn 0.005-' Al + ' )[Zn  Mn  ' Al 2--' ]O 4 formulae are listed in Table1.

	Table 1. Structural para meters of the 0.5% Mn-doped ZnAl 2 O 4 samples
		 + '*	a (Å)	Rp	R-Bragg
	500	0.159	8.09(5)	0.051	0.042
	600	0.087	8.08(9)	0.067	0.026
	800	0.059	8.09(5)	0.058	0.028
	1000	0.048	8.084(5)	0.061	0.023
	1200	0.047	8.089(4)	0.084	0.028
	1350	0.026	8.087(4)	0.121	0.039

*global inversion rate, **with all-non excluded model and not corrected for background.

ACKNOWLEDGMENTS

The work was funded by the French National Agency for Research ANR AMOS program (ANR-2010-BLANC-0820), the CNRS, the LaPhIA excellence cluster program (IDEX Université de Bordeaux), the pole de compétitivité Routes de lasers ® and the Conseil Régional d'Aquitaine.