U S Welzel

A Fréour

E J Kumar

Mittemeijer

Diffraction Stress Analysis Using Direction-Dependent Grain-Interaction Models

Keywords: diffraction stress analysis, grain interaction, surface anisotropy, grain-shape texture

Recently, so-called direction-dependent elastic grain-interaction models have been developed, which are capable of considering the effects of surface anisotropy or of a grain-shape (morphological) texture on the mechanical elastic constants and diffraction (X-ray) stress factors of polycrystals. The notion 'direction-dependent grain-interaction' signifies that different graininteraction constraints prevail along different directions in the polycrystal. This work summarizes recent developments and presents a comparison of direction-dependent elastic grain-interaction models.

Introduction

Elastic grain-interaction models describe the distribution of stresses and strains over the crystallographically differently oriented grains in a polycrystal. They can be employed for the calculation of mechanical elastic constants (relating mechanical strains to mechanical stresses) and diffraction elastic constants (diffraction (X-ray) stress factors for mechanically elastically anisotropic samples, relating (diffraction) lattice strains to mechanical stresses) from single-crystal elastic data [START_REF] Hauk | Structural and residual stress analysis by nondestructive methods[END_REF][START_REF] Welzel | [END_REF][3].

Usually, isotropic grain-interaction models, as for example the Reuss [4], Voigt [START_REF] Voigt | Lehrbuch der Kristallphysik[END_REF], Neerfeld-Hill [START_REF] Neerfeld | [END_REF]7] or Eshelby-Kröner [8,9] models, are considered. In such isotropic models, the same grain-interaction assumptions are adopted for all directions in the polycrystal. The Voigt and Reuss models involve extreme grain-interaction assumptions: either the strain tensor (Voigt model) or the stress tensor (Reuss model) for each crystallite is assumed to be equal to the mechanical strain or stress tensor, respectively; they are incompatible with physical reality because of the corresponding stress, respectively strain discontinuities at the grain boundaries. It was demonstrated by Hill [7] that the Voigt and Reuss models set bounds for the mechanical elastic constants. Further, Neerfeld [START_REF] Neerfeld | [END_REF] and Hill [7] found on an empirical basis that the arithmetic (or geometric, Hill, [7]) averages of X-ray [START_REF] Neerfeld | [END_REF] and macroscopic [START_REF] Neerfeld | [END_REF]7] elastic constants calculated according to the models of Voigt and Reuss are generally in good agreement with experimental data. Close to these values are also the results of the (mathematically more tedious) Eshelby-Kröner model.

Devised for bulk specimens with an isotropic microstructure, the above mentioned models imply mechanical elastic isotropy for crystallographically untextured polycrystals. However, a polycrystal cannot generally be considered as being mechanically elastically isotropic, even in the absence of crystallographic texture.

Thin films or surface layers of bulk-polycrystals, as examples, can exhibit mechanically elastically only transverse isotropy, owing to their microstructure and reduced dimensionality [10,11]. Thus, anisotropic grain interaction occurs. In a (columnar) thin film (or the surface layer of a bulk polycrystal), each crystallite is surrounded by neighbouring crystallites in only two dimensions. Thus, the grain-interaction perpendicular to the surface can be different from the grain interaction parallel to the surface [10,[12][13][START_REF] Welzel | Mittemeijer in: Diffraction analysis of the microstructure of materials[END_REF] (see also Stickforth [11], who introduced the notion 'surface-anisotropy' of bulk polycrystals).

Polycrystals with a grain-shape (morphological) texture also are macroscopically elastically anisotropic. The effect of a morphological texture on the elastic behaviour can be considered by employing an extension of the Eshelby-Kröner grain-interaction model [START_REF] Krier | [END_REF]16]; see, in particular, Koch et al. [17]). To this end, a calculation scheme involving ellipsoidal grains (inclusions) with their principal directions aligned along common axes in the specimen frame of reference, where the S 3 -axis is perpendicular to the sample surface, can be adopted (i.e., an ideal grain-shape texture occurs; for a detailed discussion, see [17]), whereas spherical grains (inclusions) are considered in the traditional Eshelby-Kröner model.

Mechanical Elastic Constants and X-Ray Stress Factors

In thin films, plane, rotationally symmetric states of mechanical stress/strain (mechanical averages of tensors are indicated by brackets ' < > ') are frequently met:

0 0 0 0 0 0 S S S S ε ε ε ε ⊥     =       , 0 0 0 0 0 0 0 S S S σ σ σ     =       . (1), (2)
The superscript S indicates that a tensor is expressed in the specimen frame of reference S. Under these loading conditions, which will be considered for the Vook-Witt and the inverse Vook-Witt models discussed below, two mechanical elastic constants A and B suffice, as only three independent non-zero stress and strain tensor components occur:

A S S ε σ = , (3)
B S S ε σ ⊥ = . (4)
For macroscopically elastically anisotropic specimens, the ϕ, ψ and hkl-dependent diffraction (Xray) stress factors ( , , ) ( ,, )

hkl S ij ij F hkl ϕ ψ ε ψ ϕ σ = . ( 5
)
Note that the ( , ,

) ij F hkl ψ ϕ
are not components of a tensor since they relate the lattice strain (a number) to the stress tensor (expressed in the S-system). For the case of transverse elastic isotropy and under a plane, rotationally symmetric state of stress, the diffraction strain is independent of the angle ϕ, thus, ϕ can be arbitrarily set to zero (or any other value) for the X-ray stress factors and ϕ as an index can be suppressed for the diffraction strain, so that Eq. ( 5) can be simplified:

( ) 11 22
( ,0, ) ( , 0, )

hkl S F hkl F hkl ψ ε ψ ψ σ = + . (6)
Usually, the lattice strain is plotted versus sin 2 ψ in a so-called sin 2 ψ-plot.

Thin Films / Surface Anisotropy

In (columnar) thin films and surface layers of bulk polycrystals, neighbouring grains surround a grain in only two dimensions. Grain-interaction assumptions taking this into account can be formulated as follows ( [10], see also Vook & Witt [18]; a mechanical loading state as given by Eqs.

(1) and ( 2) is considered): (i) in the plane of the film, the strain is rotationally symmetric and (ii) equal for all crystallites and (iii) the stresses perpendicular to the layer are zero for all crystallites, i.e. the crystallites can deform freely in this direction. The above formulated grain-interaction assumptions allow the calculation of the mechanical elastic constants A and B and the (sum of two, see Eq. ( 6)) diffraction stress factors [10,[12][13][START_REF] Welzel | Mittemeijer in: Diffraction analysis of the microstructure of materials[END_REF]. Extreme grain-interaction assumptions (e.g. the stresses perpendicular to the layer are zero for all crystallites) are adopted in the Vook-Witt (VW) model. In terms of the bulk models it could be said that, a Voigt behaviour is attributed to the in-plane directions, whereas a Reuss behaviour is attributed to the direction perpendicular to the film surface. Thus, the VW model will generally be incompatible with the true elastic behaviour of a polycrystal. To overcome this problem, the inverse Vook-Witt (iVW) model, required for the construction of an effective grain-interaction model, has been proposed by Welzel et al. [13,[START_REF] Welzel | Mittemeijer in: Diffraction analysis of the microstructure of materials[END_REF]. The grain-interaction assumptions for the iVW model are as follows: (i) the in-plane stress is rotationally symmetric and (ii) equal for all crystallites and (iii) the strain perpendicular to the film surface is equal for all crystallites. As in the VW model, the above formulated grain-interaction assumptions allow the calculation of the mechanical elastic constants A and B and the (sum of two, see Eq. ( 6)) diffraction stress factors [13,[START_REF] Welzel | Mittemeijer in: Diffraction analysis of the microstructure of materials[END_REF].

In an effective grain-interaction model, the mechanical and diffraction elastic constants are calculated as arithmetic averages of the elastic constants obtained from a set of extreme graininteraction models. For the following discussion such an average of (extreme) grain-interaction models will be called an effective grain-interaction model. A well known example for an effective grain-interaction is the Neerfeld-Hill model (cf. Introduction). The background of any averaging of (extreme) grain-interaction models to describe physical reality could be described as follows: A real sample is conceived to be constituted from separate volume fractions of crystallites, each of which obeys a certain type of grain-interaction.

The need for introducing the iVW model can be understood as follows: As compared to bulk materials where all directions for grain interaction can be equivalent, thin films posses two principal directions exhibiting possibly different types of grain interactions: the in-plane direction (all inplane directions are equivalent) and the direction perpendicular to the surface of the film. The number of grain-interaction models of extreme types of grain-interaction assumptions is two for bulk materials (Reuss and Voigt). Consequently, four types of extreme grain-interaction models then can be formulated for columnar thin films, as two principal directions, each with two extreme grain-interactions, occur. These extreme grain-interaction models are the Reuss, the Voigt, the Vook-Witt and the inverse Vook-Witt models.

The VW model and the iVW model imply that a polycrystal is macroscopically elastically anisotropic (transversely isotropic). Of course, this mechanical anisotropy cannot be revealed by the two elastic constants A and B (cf. Table 1 for an example). However, the occurrence of non-linear sin 2 ψ-plots for crystallographically untextured polycrystals subjected to a plane, rotationally symmetric state of stress (see Eq. ( 2)) indicates that indeed mechanical elastic anisotropy does occur; see Fig. 1, where pronounced curvature in the sin 2 ψ-plots can be observed. Experimental investigations of direction-dependent graininteraction in thin films have been conducted by van Leeuwen et al. [10] and Welzel et al. [13,[START_REF] Welzel | Mittemeijer in: Diffraction analysis of the microstructure of materials[END_REF].

Grain Shape (Morphological) Texture

In order to calculate the elastic constants of a polycrystal from single-crystal elastic data in the Eshelby-Kröner (EK) model, the crystallites surrounding an individual grain (inclusion) in a polycrystal are conceived as an elastically homogenous matrix with the elastic properties of the entire polycrystal [8,9]. Traditionally, a spherical shape of the inclusions is considered [9]. The effect of a grain-shape (morphological) texture on mechanical and diffraction elastic constants can be considered in the EK model by considering ellipsoidal inclusions with their principal axes aligned along common directions in the specimen frame of reference. The shape of the crystallites can be described by a shape parameter h , which is defined as the ratio of the principal axis of the ellipsoid in the zdirection ( 3 a ) of the specimen frame of reference and the principal axes of the ellipsoid in the x- direction ( 1a ) and the y-direction ( 2 a ) in the specimen frame of reference, respectively:

3 3 1 2 a a a a h = = (7)
Thus, the considered ellipsoids exhibit rotational symmetry with respect to the surface normal of the specimen. Only an ideal grain-shape, morphological texture is considered in the following, as only in this case unique mechanical elastic constants and X-ray stress factors can be calculated employing the EK model (for a more detailed discussion of the effect of a non-ideal morphological texture, see [17]).

Examples for mechanical and diffraction elastic constants calculated for grain aspect ratios 0.01 100 η ≤ ≤ are given in Fig. 2. For 1 h ¹ (i.e., non-spherical grains), macroscopic elastic anisotropy occurs. This is revealed by the occurrence of curved sin 2 ψ-plots (as this model allows the calculation of all components of the mechanical stiffness tensor, the mechanical elastic anisotropy, transverse isotropy in this case, could also be deduced from an inspection of this tensor; see [17] for examples). It has to be stressed that it is the preferential alignment of non-spherical grains that induces the macroscopic elastic anisotropy. Obviously, a polycrystal consisting of non-spherical grains with a random distribution of their (shape) orientations will be macroscopically elastically isotropic, as no preferred morphological direction occurs. The mechanical elastic constants exhibit a weak dependence on the grain-aspect ratio h whereas the effect on the diffraction elastic constants is more pronounced.

Comparison of the Direction-Dependent Elastic Grain-Interaction Models

At a first glance, the occurrence of surface anisotropy and the occurrence of a morphological texture bear no relation and thus, similarities among the elastic constants calculated with the VW/iVW and the EK models cannot be expected. However, upon comparing Fig. 1 with Fig. 2a, it is striking that similarities in the sin 2 ψ-plots do occur.

In Fig. 2b, the mechanical elastic constants A and B (cf. Eqs. ( 3) and ( 4)) calculated employing the EK model are shown as functions of the grain-aspect ratio h (line with open circles). For comparison, the values of A and B have been indicated also for the VW, the iVW and the Reuss and Voigt models. It can be observed that the values of A and B calculated according to the EK model fall between the corresponding VW and iVW values. As the grain shape approaches the limit of a flat disc ( )

0 h ®
, A and B tend towards the corresponding VW values, whereas in the limit of needle-like grains ( ) h ® ¥ , A and B tend towards the corresponding iVW values. Hence, the similarities among the different models revealed by the sin 2 ψ-plots also hold for the mechanical elastic constants A and B. These observations may be understood as follows.

In the case of flat-disc shaped grains (i.e., 0 η → ) with the disc faces parallel to the surface of the specimen, a polycrystal could be conceived as a stack of flat disc-shaped grains extending throughout the entire width of the polycrystal. If such a polycrystal is subjected to a plane, rotationally symmetric mechanical strain parallel to the disc interfaces, all discs (grains) will experience the same in-plane strain. Due to mechanical equilibrium conditions, they will experience the same stress ( 0 º ) perpendicular to the interfaces: the 'grain-interaction assumptions' are thus of the Vook-Witt type.

In the case of needle-shaped grains (i.e., η → ∞ ), a polycrystal could be conceived as an assembly of needle-shaped grains with the needle axes perpendicular to the surface of the specimen. If such a polycrystal is subjected to a plane, rotationally symmetric mechanical strain perpendicular to the needle axes and if one requires that the outer contour of the assembly does not change, all needles will experience the same in-plane stress in order to satisfy the mechanical equilibrium conditions. Due to the requirement of surface contour-coherence, the grains will experience the same strain along the needle axes: the 'grain-interaction assumptions' are thus of the inverse Vook-Witt type. It is possible to confirm the here observed similarities by analytical simplifications of certain tensors involved in the calculations in the Eshelby-Kröner model for the limits 0 h ® and h ® ¥ [START_REF] Welzel | [END_REF].

Conclusions

• Thin films can exhibit mechanical, macroscopic elastic anisotropy also in the absence of crystallographic texture. The elastic grain interaction of such thin films can be modelled employing an effective grain-interaction model combining the extreme Reuss, Voigt, Vook-Witt and inverse Vook-Witt models. • Polycrystals with a grain-shape (morphological) texture are mechanically elastically anisotropic.

Their elastic grain interaction can be modelled employing an extension of the traditional Eshelby-Kröner model; whereas spherical grains (inclusions) are considered in the traditional Eshelby-Kröner model, a calculation of mechanical elastic constants and diffraction stress factors is also possible considering ellipsoidal grains with their principal directions aligned along common axes in the specimen frame of reference. F F + versus sin 2 ψ (cf. Eq. ( 6)).

• For extreme grain morphologies (flat discs and needles), the grain-interaction in the Eshelby-Kröner model exhibits similarities with the grain interaction in the Vook-Witt and inverse Vook-Witt models, respectively. Even though the Vook-Witt and inverse Vook-Witt models have been developed in order to deal with the effect of elastic surface anisotropy of bulk polycrystals and thin films, they are thus capable of modelling the effect of an ideal morphological texture on the elastic grain interaction.

  for an hkl reflection at a specimen tilt angle ψ and rotation angle ϕ , to the mechanical stress tensor expressed in the specimen frame of reference S σ [1-3]: ,

Figure 2b :

 2b Figure 2b: Mechanical elastic constants A and B of a crystallographically untextured copper polycrystal according to the Vook-Witt (VW), inverse Vook-Witt (iVW) and the EK model as a function of the grain aspect ratio h (line with open circles). For comparison,the values for A and B according to the Reuss (R) and Voigt (V) models have also been indicated (singlecrystal compliances from[START_REF] Meyers | Mechanical metallurgy, principles and applications[END_REF]).

Table 1 :

 1

	Mechanical elastic constants A and B of a
	crystallographically untextured copper polycrystal (single-
	crystal elastic data from [19]) according to the Reuss (R),
	Voigt (V), Vook-Witt (VW) and inverse Vook-Witt (iVW)
	models.
	Cu	R	V	VW	iVW
	A (TPa -1 ) 5.76	4.65	5.19	4.91
	B (TPa -1 ) -6.72 -4.50 -5.59 -5.03


Alexandre Guion

Shahriar Afkhami

Stephane Zaleski

Jacopo Buongiorno

Simulations of microlayer formation in nucleate boiling

As a bubble grows outside of a cavity during nucleate boiling, viscous effects can be large enough compared to surface tension to impede liquid motion and trap a thin liquid layer, referred to as liquid microlayer, underneath the growing bubble. In practice, numerical simulations of nucleate boiling resolve the macroscopic liquid-vapor interface of the bubble, but resort to subgrid models to account for micro scale effects, such as the evaporation of the liquid microlayer. Evaporation models require initialization of the microlayer shape and extension, but models for microlayer formation are either physically incomplete or purely empirical. In this work, the Volume-Of-Fluid (VOF) method, implemented in the Gerris code, is used to numerically reproduce the hydrodynamics of hemispherical bubble growth at a wall, and resolve the formation of the liquid microlayer with an unprecedented resolution. Simulations are validated against the latest experimental data and compared to existing analytical models. Lastly, remaining gaps in building a generally applicable model for the formation of the microlayer are presented.

Introduction and motivation

The transport of latent heat makes boiling one of the most efficient mode of heat transfer, allowing a wide range of systems to improve their thermal performance, from nuclear power plants to microelectronic devices -see reviews by [START_REF] Yeh | Review of heat transfer technologies in electronic equipment[END_REF], [START_REF] Manglik | On the advancements in boiling, two-phase flow heat transfer, and interfacial phenomena[END_REF], and [START_REF] Kandlikar | History, advances, and challenges in liquid flow and flow boiling heat transfer in microchannels: a critical review[END_REF]. Subcooled boiling heat transfer for example is able to accommodate very high heat fluxes, thus cooling components very efficiently.

Predicting boiling heat transfer has therefore garnered significant attention, but remains complicated by the need to consider phenomena occurring over multiple scalessee Figure 1, from the adsorbed liquid layer at the wall at the nanometer scale -see [START_REF] Churaev | The effect of adsorbed layers on van der waals forces in thin liquid films[END_REF] and [START_REF] Chung | A review of recent progress on nano/micro scale nucleate boiling fundamentals[END_REF], up to the bubble diameter at the millimeter scale, see [START_REF] Dhir | Boiling heat transfer[END_REF], [START_REF] Stephan | Evaluation of heat and mass transfer phenomena in nucleate boiling[END_REF], [START_REF] Kim | Review of nucleate pool boiling bubble heat transfer mechanisms[END_REF], [START_REF] Kunkelmann | Modification and extension of a standard volume-of-fluid solver for simulating boiling heat transfer[END_REF], [START_REF] Guion | On the Liquid Microlayer Underneath a Vapor Bubble Growing at a Heated Wall[END_REF], [START_REF] Jiang | Dynamic modeling on bubble growth, detachment and heat transfer for hybrid-scheme computations of nucleate boiling[END_REF], [START_REF] Jung | An experimental method to simultaneously measure the dynamics and heat transfer associated with a single bubble during nucleate boiling on a horizontal surface[END_REF], [START_REF] Sato | A depletable micro-layer model for nucleate pool boiling[END_REF], [START_REF] Giustini | Evaporative thermal resistance and its influence on microscopic bubble growth[END_REF], Zou et al. (2016a), and Zou et al. (2016b).

Physical phenomena

During nucleate boiling a bubble grows outside of a surface imperfection, referred to as cavity, and from a preexisting embryo of vapor initially trapped in that surface imperfection. Boiling starts at the wall when the local temperature is sufficiently high to allow the vapor bubble to exist and grow outside of the cavity and into the surrounding liquid, see [START_REF] Bankoff | Entrapment of gas in the spreading of a liquid over a rough surface[END_REF]. At the inception of nucleate boiling at a cavity, the liquid in the vicinity of the heated wall is superheated and mass transfer occurs at the liquid-vapor interface of the bubble. The growth of the bubble is also fueled by the overpressure within the vapor, and impeded by the inertia of the liquid. During this initial phase of rapid growth, often referred to as the inertia-controlled phase of bubble growth, the bubble grows hemispherically, see Figure 2: bubble growth rates, U b , based on bubble cross section area, bubble height, and bubble width are identical. We define a capillary number Ca associated with bubble growth rate as:

Ca = µ l U b /σ
with µ l the liquid viscosity and σ the surface tension coefficient. As the bubble grows outside of the cavity, Ca can be large enough to impede liquid motion at the wall and trap a thin liquid layer underneath the bubble: the liquid microlayer -see experimental evidence by [START_REF] Moore | The measurement of rapid surface temperature fluctuations during nucleate boiling of water[END_REF], [START_REF] Hendricks | Initiation of cooling due to bubble growth on a heating surface[END_REF], [START_REF] Hospeti | Deposits formed beneath bubbles during nucleate boiling of radioactive calcium sulfate solutions[END_REF], [START_REF] Jawurek | Simultaneous determination of microlayer geometry and bubble growth in nucleate boiling[END_REF], [START_REF] Foltz | The measurement of surface temperatures with platinum films during nucleate boiling of water[END_REF], [START_REF] Judd | Laser interferometric investigation of the microlayer evaporation phenomenon[END_REF], [START_REF] Judd | A comprehensive model for nucleate pool boiling heat transfer including microlayer evaporation[END_REF], [START_REF] Koffman | Experimental study of microlayer formed during vapour bubble growth on the heated solid surface[END_REF], [START_REF] Moghaddam | Physical mechanisms of heat transfer during single bubble nucleate boiling of fc-72 under saturation conditions-i. experimental investigation[END_REF], [START_REF] Golobic | Experimental determination of transient wall temperature distributions close to growing vapor bubbles[END_REF], [START_REF] Kim | Detection of liquid-vapor-solid triple contact line in two-phase heat transfer phenomena using high-speed infrared thermometry[END_REF], [START_REF] Gao | An investigation of microlayer beneath nucleation bubble by laser interferometric method[END_REF], [START_REF] Jung | An experimental method to simultaneously measure the dynamics and heat transfer associated with a single bubble during nucleate boiling on a horizontal surface[END_REF], [START_REF] Jung | An experimental study on heat transfer mechanisms in the microlayer using integrated total reflection, laser interferometry and infrared thermometry technique[END_REF], [START_REF] Bigham | Role of bubble growth dynamics on microscale heat transfer events in microchannel flow boiling process[END_REF], [START_REF] Yabuki | Heat transfer mechanisms in isolated bubble boiling of water observed with mems sensor[END_REF], [START_REF] Yabuki | Microlayer formation characteristics in pool isolated bubble boiling of water[END_REF], and Zou et al. (2016a).

Once the bubble has reached a certain size, referred to as departure diameter, the bubble departs from its initial position at the heated wall and additional cooling mechanisms can occur.

From measurements of microlayer thickness over time by [START_REF] Kim | Detection of liquid-vapor-solid triple contact line in two-phase heat transfer phenomena using high-speed infrared thermometry[END_REF] in the case of water at atmospheric pressure, we identify two separate time scales τ 1 and τ 2 representative of microlayer formation and evaporation, respectively: τ 1 ∼ 10 -100µs, and τ 2 ∼ 1 -10ms. Hemispherical bubble growth was observed at both time scales -see Figure 2.

Modeling of boiling heat transfer

In the past, various empirical correlations and mechanistic models have been proposed to describe and predict boiling heat transfer -see reviews by [START_REF] Dhir | Boiling heat transfer[END_REF], [START_REF] Manglik | On the advancements in boiling, two-phase flow heat transfer, and interfacial phenomena[END_REF], and [START_REF] Warrier | Heat transfer and wall heat flux partitioning during subcooled flow nucleate boiling: a review[END_REF]. A wide range of boiling surface geometries, materials, working fluids, and operating conditions have been investigated. In particular, heat partitioning models by [START_REF] Judd | A comprehensive model for nucleate pool boiling heat transfer including microlayer evaporation[END_REF], [START_REF] Kurul | Multidimensional effects in forced convection subcooled boiling[END_REF], and Gilman (2014) provide a mechanistic description of the wall boiling phenomena: the total heat flux at the wall is divided into separate heat flux components that account for various modes of heat transfer at the wall, such as forced convection, evaporation, quenching, and sliding conduction. In this context, recent advanced diagnostics (Infrared Imagery, Particle Image Velocimetry, High Speed Video Imagery, and etc.) provide new and highlyresolved experimental data that can inform and validate models or correlations describing Through careful calibration and post-processing, it is possible to accurately measure the temperature and heat flux distributions on the boiling surface which reveals detailed features, including the effects of evaporation of the liquid microlayer as shown by [START_REF] Gerardi | Investigation of the pool boiling heat transfer enhancement of nano-engineered fluids by means of high-speed infrared thermography[END_REF] and [START_REF] Kim | Review of nucleate pool boiling bubble heat transfer mechanisms[END_REF]. Depending on the fluid, flow and nucleation temperature, the contribution from microlayer evaporation to overall heat transfer and bubble growth can be large in the case of water, see [START_REF] Gerardi | Investigation of the pool boiling heat transfer enhancement of nano-engineered fluids by means of high-speed infrared thermography[END_REF], or relatively small in the case of refrigerants, see [START_REF] Kim | Review of nucleate pool boiling bubble heat transfer mechanisms[END_REF]. Typical microlayer models focus on evaporation at the Triple-Phase-Line (TPL) region where the three phases (liquid, solid, and vapor) are in apparent contact, see [START_REF] Dhir | Simulation of boiling: How far have we come[END_REF] and [START_REF] Kunkelmann | Modification and extension of a standard volume-of-fluid solver for simulating boiling heat transfer[END_REF], or in the extended microlayer region, see [START_REF] Guion | On the Liquid Microlayer Underneath a Vapor Bubble Growing at a Heated Wall[END_REF], [START_REF] Jiang | Dynamic modeling on bubble growth, detachment and heat transfer for hybrid-scheme computations of nucleate boiling[END_REF][START_REF] Sato | A depletable micro-layer model for nucleate pool boiling[END_REF][START_REF] Sato | A depletable micro-layer model for nucleate pool boiling[END_REF]. In practice, numerical simulations of boiling resolve the macroscopic liquid vapor interface of boiling bubbles, but typically resort to subgrid models to include contributions at the microscopic scale such as from evaporation of the microlayer, which requires initialization of the microlayer shape and extension.

R b ∼ 1 mm δ ∼ 1 -10 µm δ r θ r min r out III IV r max I II III θ h ∼ 0.1µm ∼ 1 µm

Current gaps and proposed strategy

Existing models of microlayer formation are either physically incomplete, e.g. do not include the effect of surface tension or contact angle, see [START_REF] Van Stralen | Bubble growth rates in pure and binary systems: combined effect of relaxation and evaporation microlayers[END_REF], or are purely empirical, see [START_REF] Sato | A depletable micro-layer model for nucleate pool boiling[END_REF]. The objective of the present work is to leverage existing numerical methods to build a large numerical database and inform mirolayer formation, clarify remaining modeling challenges, and eventually build a generally applicable model for the formation of the microlayer. Existing experimental data that directly measures microlayer thickness profiles, see work from [START_REF] Judd | Laser interferometric investigation of the microlayer evaporation phenomenon[END_REF], [START_REF] Koffman | Experimental study of microlayer formed during vapour bubble growth on the heated solid surface[END_REF], [START_REF] Kim | Detection of liquid-vapor-solid triple contact line in two-phase heat transfer phenomena using high-speed infrared thermometry[END_REF], Jung &[START_REF] Jung | An experimental method to simultaneously measure the dynamics and heat transfer associated with a single bubble during nucleate boiling on a horizontal surface[END_REF]Jung &[START_REF] Jung | An experimental study on heat transfer mechanisms in the microlayer using integrated total reflection, laser interferometry and infrared thermometry technique[END_REF], does not cover a wide range of fluids and surfaces, which makes it difficult to generate or validate a general model on microlayer formation. In practice, measurements of microlayer thickness profiles are mostly obtained at times where the microlayer has already started to evaporate, with little insight on its formation dynamics. Recent advances by [START_REF] Yabuki | Microlayer formation characteristics in pool isolated bubble boiling of water[END_REF] in indirect measurements of microlayer thickness at short time scales can offer new insights and will be discussed. In the present work, we use Computational Fluid Dynamics with Interface Tracking Methods to numerically reproduce the hydrodynamics of bubble growth at a wall, and resolve the formation of the liquid microlayer at the wall for a wide range of conditions and fluids. For all simulated cases, the Volume-Of-Fluid Method implemented in the Gerris flow solver from [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF] is used to accurately track the liquid vapor interface over time, and predict the formation of the microlayer.

Assumptions made in this work

In this work, we focus on the mechanisms of microlayer formation. In particular, we only consider the hydrodynamics of bubble growth at the wall, and neglect thermal effects, separating microlayer formation (τ 1 ∼ 10 -100µs) from microlayer evaporation (τ 2 ∼ 1 -10ms τ 1 ). Part of the liquid layer that forms at the wall may also evaporate during such short time scales τ 1 , hence affecting its shape over time. However, we estimate that only 10% or less of the thickness of the liquid microlayer would be affected by evaporation during the formation phase, for a typical wall superheat at the wall ∆T (∆T = 1 to 10K) observed experimentally at atmospheric pressure, and considering evaporation from pure conduction through the liquid layer, for various growth times (1 to 100µs), and various layer thicknesses (0.1 to 10µm). Thermal effects, including the effect of mass transfer, on the motion of the contact line at short time scales are not included in the scope of the present work. Long range interaction forces of the liquid film with the substrate are also neglected as they may be relevant only at very short time scales where the film is very thin -see for example [START_REF] Rednikov | Singularity-free description of moving contact lines for volatile liquids[END_REF]. In this work, we implement a static contact angle at a wall for simplicity, which allows the contact line to move. However, one could replace it with a specific mobility law that would describe the behavior of the microscopic contact angle in the conditions of interest. Mesh dependence of moving contact lines simulations with grid size as previously shown by [START_REF] Moriarty | Effective slip in numerical calculations of movingcontact-line problems[END_REF], [START_REF] Weinstein | Scale dependence of contact line computations[END_REF], and [START_REF] Afkhami | A mesh-dependent model for applying dynamic contact angles to vof simulations[END_REF] is reported elsewhere, see [START_REF] Afkhami | Transition in a numerical model of contact line dynamics and forced dewetting[END_REF].

Contributions and structure of the paper

This work focuses on three main contributions: (i) identify the minimum set of dimensionless parameters that controls the hydrodynamics of microlayer formation at a wall during nucleate pool boiling (ii) generate a numerical database that resolves the dynamics of microlayer formation for all conditions of interest (e.g. water in lab experiment at atmospheric pressure, water at nuclear reactor pressure conditions), and (iii) identify remaining gaps in building a generally applicable model for microlayer formation, to be used during the initialization of boiling simulations. This paper specifically focuses on the central microlayer region, see region III in Figure 1. Future work includes modeling the outer-edge region, see region IV in Figure 1, using the present numerical database, and modeling the particular contact line dynamics using experimental data for specific wall conditions and for the fluid/solid pair of interest.

In section 2 we provide a dimensional analysis and general functional dependence for the microlayer formation dynamics. In section 3 we describe the computational framework used to simulate microlayer formation, and verify the correct implementation and convergence of the methods used. In section 4 we discuss the reduction of the parameter space, qualitative trends, and numerical results used to model microlayer formation, in particular we compare our results in the central region (region III, Figure 1) with the analytical model from [START_REF] Cooper | The microlayer in nucleate pool boiling[END_REF] and with latest experimental results measuring microlayer thickness at early times. In section 5 we conclude by discussing the contributions of this present work and the remaining opportunities for future work.

Problem statement and dimensional analysis

We consider the early growth of a steam bubble near the wall of a pool of initially stagnant liquid. We take advantage of the symmetry of the growth and consider the problem to be axisymmetric. Ten physical variables are involved in describing the hydrodynamics of microlayer formation: liquid and vapor viscosities, µ l and µ v , respectively; liquid and vapor densities, ρ l and ρ v , respectively; surface tension, σ; bubble growth rate, U b , the velocity at which the liquid/vapor interface moves inside the surrounding liquid; microscopic contact angle between the liquid vapor interface and the wall at a given reference length scale dx, the minimal mesh size, θ; radial distance from bubble root, r; time, t; and the unknown local thickness of the liquid microlayer forming at the wall, δ. Altogether, these ten variables involve three physical dimensions, and seven dimensionless Pi-groups can fully describe microlayer formation. The proposed dimensionless Pi-groups are:

Π 1 = δ/r c = δ * Π 2 = µ l /µ v = µ * Π 3 = ρ l /ρ v = ρ * Π 4 = r/r c = r * Π 5 = t/t c = t * Π 6 = θ Π 7 = Ca = µ l U b /σ
where the reference length and time scales r c and t c were defined as follows:

r c = µ l /(ρ l U b )
(2.1)

t c = r c /U b = µ l /(ρ l U 2 b ) (2.2)
with r c based on Reynolds number Re = 1, for which inertia of the liquid surrounding the bubble and viscous forces at the wall would be comparable. A unique functional form describes the microlayer shape δ * as a function of the six other Pi-groups, which can be written as follows:

δ * = f (µ * , ρ * , r * , t * , Ca, θ) (2.3)
In this study we have performed numerical simulations of microlayer formation to quantify f (Eq. 2.3). The relevant parameter space for liquid water boiling in lab and industrial conditions (Pressurized Water Reactors, or PWR) is as follows:

µ l /µ v ∈ [3.2; 23.3], ρ l /ρ v ∈ [5.9; 1650], Ca ∈ [0.001; 0.1], θ ∈ [0; 90 • ].
The range of capillary number, Ca, is obtained from water properties at saturation (µ l and σ) for both pressures (0.101MPa and 15.5MPa), and estimates of bubble growth rates U b from [START_REF] Mikic | On bubble growth rates[END_REF], including a range of cavity mouth radius r b :

U b = πh f g ρ v ∆T sat /(7ρ l T sat ) (2.4)
where ∆T sat is the so-called wall superheat at the inception of boiling:

∆T sat = 2σT sat /(r b ρ v h f g ) (2.5)
with T sat the saturation temperature, and h f g the latent heat of vaporization. In practice, bubble growth rate U b and wall superheat ∆T sat , estimated using Eq. 2.4 and 2.5, are experimentally measured.

Computational framework

A myriad of methods are available today to simulate two phase flows while tracking the interface between two fluids -see recent review by [START_REF] Kharangate | Review of computational studies on boiling and condensation[END_REF]. The Volume-Of-Fluid (VOF) method, see [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF], tracks the interface using a color function C that represents the volume fraction of one phase in each cell of the computational domain. Starting from an initial distribution, the color function C is then advected by the flow:

∂Cρ/∂t + ∇.(Cρ U ) = S (3.1)
The Gerris Flow Solver by [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF] (gfs.sourceforge.net/) solves the incompressible Navier-Stokes equations, and implements the VOF method to track the interface and compute the surface tension force:

ρ(∂ U /∂t + U .∇ U ) = -∇p + ∇.(2µ D) + σκδ s n + F external (3.2) ∂ρ/∂t + ∇.(ρ U ) = S (3.3) ∇. U = S (3.4)
where S is a volume source such that ∇. U = S; U is the velocity field; p is the pressure; κ is the curvature of the liquid vapor interface; n is the unit vector normal to the liquid vapor interface; D is the rate of strain D = ∇U + (∇U ) T ; ρ is the fluid density,

ρ = Cρ l + (1 -C)ρ v ; µ is the fluid viscosity, µ = Cµ l + (1 -C)µ v .
The computational domain in Gerris is spatially discretized using a quad mesh, which allows refining specific regions in space while coarsening others. In this work, we typically refine in the vicinity of the wall and at the interface to resolve microlayer formation, while coarsening away from the wall and the bubble interface where such refinement is not necessary. We include a time-dependent volumetric source of vapor within the vapor phase to reproduce the growth observed in lab experiments. The source of vapor fuels bubble growth uniformly inside the bubble, and pushes the liquid vapor interface at the desired rate U b , hence reproducing the hydrodynamics of hemispherical bubble growth. A timedependence of the source is necessary to reproduce the sustained growth of the bubble as its volume increases. For a given constant bubble growth rate U b (U b = dR b /dt) during the initial inertia-controlled phase of bubble growth, one can compute the volumetric

dV b /dt = 2πR 2 b (dR b /dt) = 2πR 2 b U b (3.5)
The rate of change of the equivalent hemispherical bubble dV b /dt is obtained by imposing a volumetric source term S in the vapor phase of volume V b :

SV b = dV b /dt (3.6)
Rearranging Eq. 3.6, we obtain an explicit expression for the volumetric source term S to impose in the vapor phase:

S = (dV b /dt)/V b = 3U b /(R b,0 + U b t) (3.7)
Note the volumetric source term S is uniform in space within the vapor phase at all times, yet time dependent, and function of three parameters: the bubble growth rate U b known from experiments or derived from Mikic's theory, see [START_REF] Mikic | On bubble growth rates[END_REF], the initial bubble radius R b,0 prescribed during the initialization of the numerical simulation, and the physical time t. We typically simulate the growth of the bubble outside of the cavity, and initialize the bubble radius R b,0 as the cavity mouth radius r b (see Eq. ??).

Method verification

We verify the implementation of the volumetric source term in axisymmetric configuration in Gerris in the absence of a wall, for the special case of a spherical bubble growth in a (semi)-infinite medium. We compare the volume and shape of the bubble to its simple analytical and spherical reference, and confirm convergence with grid size.

We also verify the convergence of the proposed method in the presence of a wallsee Figure 3. Two figures of merit are considered: the shape of the liquid microlayer forming at the wall, and the macroscopic bubble growth rate obtained from imposing the volumetric source of vapor inside the bubble. Four simulations of bubble growth at a wall are reported here, with increased refinement in the microlayer region: dx/R b,0 = 1/40, 1/80, 1/160, 1/320 (R b,0 = 10µm), specifically in a band of thickness 0.4R b,0 at the wall. Refinement elsewhere at the interface is kept identical: dx/R b,0 =1/10. All four simulations yield identical macroscopic bubbles shapes -see Figure 4(a). In Figure 4(b), the microlayer shape converges to a single profile as we refine the mesh in the microlayer region, from dx = 250 nm down to 125 nm, 62.5 nm, and 31.25 nm.

Numerical simulation results

Four regions (and scales) can be identified based on the curvature of the bubble interface, see Figure 5: (i) the contact line region where the interface meets the wall (negative curvature), (ii) the central region of the liquid microlayer (≈ zero curvature), (iii) the microlayer outer-edge region that connects the central region with the macroscopic bubble cap (up to the bubble nose where the curvature peaks), (iv) the macroscopic hemispherical bubble cap (uniform curvature along the interface).

Reducing the parameter space

We perform simulations with all four bounding values for both viscosity and density ratios, and we note the weak dependence of the microlayer shape within both ranges of interest. This weak dependance is supported by the argument that the properties of the liquid phase affect the hydrodynamics of this problem much more than those of the vapor phase, which in practice can be considered as an inviscid and massless fluid. Therefore, in the rest of the study, we will no longer consider the dependence of the microlayer formation on the viscosity and density ratios, which simplifies the function in Eq 2.3 for the conditions of interest.

Figure 5: Typical axisymmetric shape of a growing bubble at a wall. Interface is solid black line. Curvature κ is given by the colormap: blue in the contact line region is negative curvature (minimum reached at the contact point κ = -6.6 10 -1 [1/µm]), green in the central linear region is approaching zero curvature (κ ≈ 10 -3 [1/µm]), red at the bubble nose is maximum curvature (κ = 5.8 10 -2 [1/µm]), and yellow is the uniform curvature along the interface at the bubble cap as the bubble grows hemispherically (κ ≈ 1.5 10 -2 [1/µm]). The range used in the color map does not reflect the maximum positive value and minimum negative value taken by the curvature, but rather is centered around 0 to allow easier identification of all four regions, also denoted in the image with direct labeling.

Modeling the microlayer inner edge

A general model of moving contact lines is out of the scope for the present paper. However, we identify multiple scenarios that could take place at the inner edge of the microlayer, where the liquid/vapor interface meets the solid surface: i) no motion of the contact line, pinned either at the cavity mouth or elsewhere on the wall surface, ii) intermittent motion of the contact line, where pinning and depinning of the contact line occurs at the wall, and iii) continuous motion of the contact line, that could be described by a mobility law θ = g(Ca cl ), with Ca cl the capillary number associated with the motion of the contact line, and θ the contact angle. When the contact line is free to move (no pinning), its motion relative to the bubble growth dictates whether a microlayer forms at the wall or not.

Qualitatively, the radial component of the surface tension force applied to a liquid wedge decreases in magnitude for increasing contact angles. Hence the contact line is expected to move faster at higher contact angles and the wetted fraction underneath the bubble is expected to be reduced at higher contact angles. Quantitatively, physical models or experimental data on contact line motion are needed to conclude on whether a microlayer is expected to form or not, for a given fluid/solid pair and surface condition, which cannot be resolved by the Gerris flow solver alone. In this work, we simulate moving contact lines for a wide range of constant and uniform contact angles at the wall, ranging from 20 • to 90 • . An extension of this work could include a physical mobility law to apply the proposed methodology to a specific application for which the physical motion of the contact line is known.

The simulation of moving contact lines is made possible with Gerris through the use of a height-function based model of contact angle implemented by [START_REF] Afkhami | Height functions for applying contact angles to 2d vof simulations[END_REF], where a numerical contact angle θ is imposed in the cell at the wall containing the contact line. This microscopic contact angle is used to define the normal to the interface at the contact line, and therefore impacts both the computed curvature and the resulting surface tension force at the contact line. Cell normal velocities are used to advect volume fractions, resulting in an implicit slip along no-slip boundaries such as along the wall, which allows the contact line to move at the wall with the prescribed microscopic contact angle θ. As a consequence, the dynamics of the simulated moving contact line depend both on the microscopic contact angle prescribed and on grid size. Specific mesh convergence studies are reported in other papers, including most recently in [START_REF] Afkhami | Transition in a numerical model of contact line dynamics and forced dewetting[END_REF].

Typical simulation results are shown in Figure 6: we confirm that the wetted fraction underneath the bubble is reduced for higher contact angles, and we note the formation of a bulge of liquid at the inner edge of the microlayer. The surface tension force acts at the inner edge against the viscous force at the wall and pushes the contact line outwards, similar to the classic hole expansion problem, see [START_REF] Taylor | On making holes in a sheet of fluid[END_REF], [START_REF] Brenner | On the bursting of viscous films[END_REF], [START_REF] Savva | Viscous sheet retraction[END_REF], [START_REF] Gordillo | Asymptotic behavior of a retracting two-dimensional fluid sheet[END_REF][START_REF] Agbaglah | Longitudinal instability of a liquid rim[END_REF]. The Ohnesorge number Oh is defined as:

Oh = µ l / δρ l σ = Ca/δ *
Three regimes are identified based on the Ohnesorge number: Oh < 0.1, Oh ∈ [0.1; 10], and Oh > 10. As Oh increases, the capillary waves found at the tip disappear through the action of viscosity. In all simulations presented in this paper, we find that Oh ranges from 0.001to0.1, which corresponds to the first regime where capillary waves are expected in the vicinity of the tip, and damped away from the tip due to viscous effects, leaving the rest of the extended liquid layer unaffected by the motion of the rim, as confirmed in Figure 6.

The microscopic contact angle θ used in Gerris allows us to represent the behavior of the moving contact line as a boundary condition, or subgrid model, without attempting to resolve the complex physics at play at the scale of the contact line. For simplicity in the rest of the paper, we present results for low microscopic contact angles in order to observe microlayer formation. However, the informed reader would need to consider the specific motion of the contact line for the wall condition and solid/fluid pair of interest in order to determine whether a microlayer forms or not in such conditions, and quantify the relative motion of the contact line with respect to the macroscopic bubble.

Modeling the microlayer central region

An analytical model for microlayer formation has been proposed by [START_REF] Cooper | The microlayer in nucleate pool boiling[END_REF], using a simplified, single-phase boundary layer description of the two phase flow problem:

δ = C 0 ν l t g (4.1)
where t g is the growth time, C 0 a dimensionless constant, and ν l = µ l /ρ l . Under constant bubble growth rate U b , we obtain the dimensionless form:

δ * = C 0 (µ l /ρ l ) * (r/U b ) = C 0 √ r * (4.2)
With an initial bubble dimensionless radius R * b,0 , and assuming that the contact line is either pinned at its initial position (r * = R * b,0 ) or that its outward motion is negligible compared to the macroscopic bubble growth, then the square root model becomes:

δ * = C 0 r * -R * b,0 (4.3)
We present in Figure 7 the square root model along with our numerical simulation results, for the range of capillary numbers of interest, Ca ∈ [0.001; 0.1], and for a low contact angle θ = 10 • for which we observe the formation of the microlayer with negligible motion of the contact line compared to the radial extension of the bubble growth. Numerical simulations at higher contact angles were also performed for sake of generality, with contact angles ranging from 10 • to 90 • . As long as a microlayer forms, higher contact angle conditions yield identical microlayer profiles away from the contact line region, and do not affect the overall volume of liquid trapped underneath the bubble, as shown in the previous section and in Figure 6.

More specifically, we plot the simulated microlayer profiles δ * as a function of dimensionless radial distance r * -R * b,0 on the left panels of Figure 7, for three consecutive times, and for capillary numbers spanning the range of interest: Ca=0.001 (a), Ca=0.01 (c) and Ca=0.1 (e). We note that the microlayer continuously grows over time, and that the portion of the thickness profile that has already formed at a given time remains unchanged at later times, with the exception of the contact line region where the dewetting process increases liquid build-up over time at the inner-edge.

The volume of liquid V in associated with the build-up at the inner edge is typically negligible compared to the total volume of liquid V tot trapped under the bubble (V in /V tot < 0.1%), and will not be modeled in the later portion of the paper. This simplification only applies if the contact line's outward motion is negligible compared to the bubble growth, which depends on the physical behavior of the contact line for a given fluid/solid pair and surface condition. In our numerical experiments, and for the contact angle model described earlier, we note in [START_REF] Guion | Modeling and simulation of liquid microlayer formation and evaporation in nucleate boiling using computational fluid dynamics[END_REF] that the simplification breaks down for very low capillary numbers (Ca → 0) and large contact angles (θ → 90 • ), where the velocity of the contact line approaches the bubble growth rate, and no microlayer is left behind.

In Figures 7 (b) (d) (f), we plot the ratio δ * / r * -R * b,0 as a function of dimensionless radial distance r * -R * b,0 . We note that the ratio is constant in the central region of the microlayer, and that the magnitude of the ratio is close to 0.5 for all the considered capillary numbers. From this observation, we use C 0 = 0.5 in the square root model prediction, and compare it with our simulation results in Figures 8,9 and 10, for various times and all Ca of interests.

Modeling microlayer formation is particularly important to capture two key boiling characteristics: i) the wall heat flux, and ii) the total volume of liquid trapped underneath the bubble that evaporates and that fuels bubble growth up to its departure from the heated surface. The square root model is able to capture both thickness and volume of the microlayer in the central region, which corresponds to the thin part of the microlayer.

We also conclude from Figures 8,9 and 10 that the square root model does not capture the outer-edge region of the microlayer, and therefore the total volume of liquid trapped underneath the bubble is underestimated. This conclusion suggests the need for specific modeling of the transition from the central (thin) region of near-zero curvature, to the outer-edge (thick) region, of non-zero curvature which is left for future work.

Validating numerical results with experimental data and existing models

Initially in [START_REF] Cooper | The microlayer in nucleate pool boiling[END_REF], experimental data on organic liquids boiling at low pressures provided estimates of C 0 , in the range of 0.5 to 1.0, for which experimental data and square root model were found to agree within ±25%.

Recent advances in microlayer thickness measurements at different radial locations and short time scales by [START_REF] Yabuki | Microlayer formation characteristics in pool isolated bubble boiling of water[END_REF] uncovered two regimes based on a modified Bond number Bo: Bo < 13 : C exp = 0.13Bo 0.38 (4.4)

Bo 13 : C exp = 0.34 (4.5)

with Bo = ρ l r 2 U/tσ, and C exp the experimental measurement of C 0 . Rearranging Bo we obtain Bo = ρ l U 2 tU/σ = t * /Ca. The transition at Bo = 13 correspond to physical times t = Bo Ca t c in the range of 10 -3 µs to 10 -1 µs for the case of boiling water at 0.101M P a (Ca ∈ [0.001; 0.1], and t c ∼ 0.1µs). Therefore, only the second regime is relevant to microlayer formation, which occurs over tens of µs, and we note C exp = 0.34 the corresponding experimental measurement of C 0 from [START_REF] Yabuki | Microlayer formation characteristics in pool isolated bubble boiling of water[END_REF], which is of the same order of magnitude as the value of 0.5 we obtain in our numerical simulations.

Lastly, we use experimental measurements of microlayer thickness to validate our numerical results, for the case of water at saturation at 0.101M P a (data courtesy of Prof. Hyungdae Kim). In the experiment, a high-speed video camera was used to measure bubble radius over time, and an IR camera was used to measure the wall temperature distribution underneath the bubble. From the bubble height and width measurements shown in Figure 2, we conclude that the bubble grew hemispherically. An additional high-speed video camera was used to measure microlayer thickness. The measurement error was assessed using a solid convex lens of known curvature: the error in thickness measurement is found to be ∼ 10% for thicknesses greater than 0.5µm, and ∼ 30% otherwise. The earliest bubble growth rate was measured at t = 20µs, see in Figure 11, and drops dramatically at early times which indicates that the initial bubble growth rate during the first tens of µs may have been much higher than its first available measured value of 4.2m/s. We include this uncertainty in initial bubble growth rates in our validation below through a sensitivity study on the initial bubble growth rate: in blue, we assume the bubble growth rate to be constant during the first 20µs, and in red we assume its acceleration to be constant during the first 20µs. We compare our numerical results against experimental results in Figure 11 (b): black squares represent the earliest experimental measurement of the microlayer thickness, which took place at time t = 410µs; and grey triangles are obtained by deducing the microlayer thickness at t = 20µs, assuming that the heat transfer within the microlayer is purely by conduction, see [START_REF] Guion | On the Liquid Microlayer Underneath a Vapor Bubble Growing at a Heated Wall[END_REF]:

∂δ(r, t) ∂t = - k l ∆T sat (r, t) ρ l h f g δ (4.6)
with k l the liquid thermal conductivity, ∆T sat the local and time-dependent wall superheat, ρ l the density of the liquid and h f g the latent heat of vaporization. To simplify the calculation of this estimate, we assume constant wall superheat:

∆T sat (r, t) = ∆T sat,0
for which a closed form exist:

δ(r, t 0 ) = δ(r, t 1 ) 2 + 2k l ∆T sat,0 ρ l h f g (t 1 -t 0 ) (4.7)
with t 0 = 20µs, t 1 = 410µs and ∆T sat,0 = 11.7 • C in our case. Assuming constant wall superheat results in under estimating the thermal coupling with the wall surface, and therefore over estimating the volume of liquid evaporated, in particular in the region of small thickness δ (e.g. low radial distance r). The comparison in Figure 11 between modeled and measured microlayer thicknesses at t = 20µs show both qualitative and quantitative agreements. We also note in Figure 11 a departure from the model in the low thickness region. We expect this discrepancy to be reduced if we no longer assume constant wall superheat, but include the thermal coupling between the microlayer and the heated wall. Such coupling would require to solve for the heat equation in the substrate, while using Eq. 4.6 as a flux boundary condition at the wall surface, as shown in previous work by [START_REF] Guion | On the Liquid Microlayer Underneath a Vapor Bubble Growing at a Heated Wall[END_REF]. The initial evaporation of the microlayer would locally reduce the wall superheat and therefore slow down the rate of evaporation of the microlayer. This slow down effect of the microlayer evaporation due to conjugate heat transfer with the solid surface is currently not included in Figure 11, and could explain the over estimated values of the microlayer thickness in the low thickness region at t = 20µs.

The following properties of water at saturation at 0.101M P a were used: µ l = 3.5 10 -4 P a.s, ρ l = 950kg/m 3 . The first measured value of bubble growth rate was U b (t ∼ 20µs) ∼ 4.2m/s.

Summary and conclusions

In this work, we used Computational Fluid Dynamics with Interface Tracking Methods to numerically reproduce the hydrodynamics of bubble growth at a wall, and resolve the formation of the liquid microlayer at the wall with high resolution in space and time. We identified the minimum set of dimensionless parameters that controls microlayer formation, and reduced the parameter space for the conditions and fluid of interest, namely water at pressures ranging from 0.101M P a (lab conditions) and up to 15.5M P a with δ * the shape of the extended liquid microlayer, r c = µ l /(ρ l U b ) and t c = r c /U b the characteristic length and time scales used, respectively. The fluid analyzed in this study is water, however because of the use of dimensionless parameters, the results presented here are broadly applicable to other fluids with dimensionless parameters within the same ranges.

We found that three regions need to be modeled in order to accurately represent the thickness of the microlayer, and the total volume of liquid trapped underneath the bubble: (i) the inner region: our qualitative results confirm that the dynamics of the moving contact line at the inner edge of the microlayer dictate whether a microlayer forms or not. The proposed methodology can be used to include a physical mobility law in place of the current static microscopic contact angle used, and simulate the liquid build-up at the inner edge for a given fluid/solid pair and surface conditions for which the motion of the contact line is known.

(ii) the central region: our simulation results indicate that the square root model proposed by [START_REF] Cooper | The microlayer in nucleate pool boiling[END_REF] captures the central region of the microlayer accurately in the case of hemishperical bubble growth: δ * ≈ 0.5 r * -R * b,0

(5.2) as long as the motion of the contact line is negligible compared to the bubble growth rate. In addition, we benchmark our results against the earliest microlayer thickness measurements available (courtesy of Prof. Hyungdae Kim), and conclude on both qualitative and quantitative match.

(iii) the outer region: we note that the square root model breaks down at the microlayer outer-edge, where its curvature rapidly increases to meet the non-zero curvature of the bubble cap.

These results clarify the remaining steps needed to build a general model on microlayer formation that resolves the initial and rapid phase of bubble growth in nucleate boiling, called inertia-controlled phase, which lasts a tens of µs in the case of water at atmospheric pressure. The next phase of bubble growth, the thermal diffusion controlled phase, is much slower and lasts tens of milliseconds in the case of water at atmospheric pressure. In practice, the microlayer formation model allows initialization of boiling simulations at the end of the rapid inertia-controlled phase of bubble growth during which the microlayer forms.

Lastly, we conclude that the numerical methods used in this study have proved to be mature enough to investigate microlayer formation and its parameter space of interest, within reasonable simulation time while achieving sufficiently fine spatial and temporal resolution. However, additional experimental data on moving contact lines and microlayer thickness profiles at very short time and length scales relevant to microlayer formation are still needed to continue to benchmark numerical simulations and unlock potential predicting capabilities.

Figure 1 :

 1 Figure1: Multiple scales involved in nucleate boiling: from the millimeter scale at the bubble cap (top) to the nanometer scale at the adsorbed liquid layer (Region I, right insert). Region I: adsorbed layer region of thickness h 0 ∼ 1 -10nm. Region II: transition region between the non-evaporating adsorbed layer and the evaporating liquid microlayer (h h 0 and h 100nm). Region III: central microlayer region of typical thickness δ ∼ 1 -10µm. Note in the insert the black dotted line that represents the Triple Phase Line (TPL) where the three phases(liquid, vapor, and solid) are in apparent contact. The contact angle θ represents the angle between the contact line and the substrate at the microscopic scale. Region IV: microlayer outer-edge region of typical thickness δ 10µm, where the thickness rapidly increases from the central region to the bubble cap.

Figure 2 :

 2 Figure 2: Bubble growth rate U b as a function of time based on bubble cross section area, bubble height, and bubble width recorded by a High Speed Video camera on the side of the growing bubble (Courtesy of Prof. Hyungdae Kim). Fluid is water at saturation and atmospheric pressure. Typical time scales τ 1 ∼ 10 -100µs, and τ 2 ∼ 1 -10ms are representative of microlayer formation and evaporation, respectively. Hemispherical bubble growth was observed at both time scales.

Figure 3 :

 3 Figure 3: Computational domain for implementing a volumetric source of vapor to control the growth of the bubble in the presence of a wall.

Figure 4 :

 4 Figure 4: Axisymmetric shapes of the macroscopic bubble (a) and liquid microlayer (b) obtained in four simulations of bubble growth at a wall, with increased refinement dx * ∈ {1, 1/2, 1/4, 1/8} in the band: 0 δ * 160 at the wall, while the refinement elsewhere at the interface is kept identical: dx * = 4, with R * b,0 = 40 and r c = 0.025 µm.

Figure 6 :

 6 Figure 6: Microlayer thickness profiles δ * for two contact angles θ = 20 • (blue solid line) and θ = 80 • (black dashed line), as a function of dimensionless radial distance r * -R * b,0 . Ca=0.1, t * = 3.6 10 3 (r c = 10nm, t c = 1ns).

Figure 7 :

 7 Figure 7: Microlayer dimensionless profiles δ * plotted in left column: (a), (c), (e), and ratios δ * / r * -R * b,0 plotted in right column: (b), (d), (f), as a function of dimensionless radial distance r * -R * b,0 , with R * b,0 the initial dimensionless bubble radius. Each row corresponds to a different Ca number: Ca=0.001 in top row (a)(b), Ca=0.01 in middle row (c)(d), Ca=0.1 in bottom row (e)(f). Microscopic contact angle is θ = 10 • , and finest mesh size at the wall is 31.25nm. Same physical times are plotted in all panels: t = 10µs (blue solid line), t = 20µs (black dashed line), t = 30µs (red dot-dashed line). In each row, the following characteristic length and time scales are imposed: (r c = 1nm, t c = 1ns) in top row (a)(b), (r c = 10nm, t c = 10ns) in middle row (c)(d), (r c = 0.1µm, t c = 0.1µs) in bottom row (e)(f).

Figure 8 :

 8 Figure 8: Microlayer dimensionless profiles δ * (left column) and subsequent volume V * (right column) with V * = V /V tot , where V tot is the total volume of liquid trapped underneath the bubble), as a function of dimensionless radial extension (r * -R * b,0 )/t * , from simulations (black solid line), and from the square root model: δ * = 0.5 r * -R * b,0 (blue dashed line). Ca=0.001 (top row), Ca=0.01 (middle row), Ca=0.1 (bottom row). Physical time t = 10µs.

Figure 9 :

 9 Figure 9: Microlayer dimensionless profiles δ * (left column) and subsequent volume V * (right column) with V * = V /V tot , where V tot is the total volume of liquid trapped underneath the bubble), as a function of dimensionless radial extension (r * -R * b,0 )/t * , from simulations (black solid line), and from the square root model: δ * = 0.5 r * -R * b,0 (blue dashed line). Ca=0.001 (top row), Ca=0.01 (middle row), Ca=0.1 (bottom row). Physical time t = 20µs.

Figure 10 :

 10 Figure 10: Microlayer dimensionless profiles δ * (left column) and subsequent volume V * (right column) with V * = V /V tot , where V tot is the total volume of liquid trapped underneath the bubble), as a function of dimensionless radial extension (r * -R * b,0 )/t * , from simulations (black solid line), and from the square root model: δ * = 0.5 r * -R * b,0 (blue dashed line). Ca=0.001 (top row), Ca=0.01 (middle row), Ca=0.1 (bottom row). Physical time t = 30µs.

Figure 11 :

 11 Figure 11: Measured bubble growth rates (black squares) over time (a), and two assumed initial regimes: constant growth rate (blue dash-dotted line), or constant acceleration (red dashed line). Thickness δ (in µm) measured in lab experiment (black squares) of boiling water at saturation at 0.101M P a, 0.41ms after boiling inception. The deduced microlayer thickness at t = 20µs (grey triangles) is obtained assuming that the heat transfer within the microlayer is purely by conduction. The wall temperature at boiling inception is 111.7 • C, and the wall heat flux is 209kW/m 2 (courtesy of Prof. Hyungdae Kim). Predictions from the square root model are plotted for the two initial profiles of U b showed in (a). The first measured value of U b is U b (t ∼ 20µs) ∼ 4.2m/s.

† Email address for correspondence: jacopo@mit.edu


Optimal Location of Actuators and Sensors in Active Vibration Control

Keywords: Active Vibration Control, Optimal location and number, actuators and sensors, numerical results I

Several studies have been developed in order to find the optimal location of actuators and sensors in active control of structures. In this paper, a modified optimization criteria is proposed for these two optimization problems, ensuring good observability or good controllability of the structure, and considering residual modes to limit the spillover effects. Its efficiency is shown by comparison with classical criteria, illustrated for a simply supported beam and a rectangular plate. In these two applications, the number of active elements is discussed, using or neglecting the residual modes.

INTRODUCTION

In vibration suppression of structures, some parameters, such as location of actuators and sensors, have a major influence on the performance of the control system. It is well known that misplaced sensors and actuators lead to problems such as the lack of observability or controllability. In this way, different cost functions have been used to find the optimal locations of these active elements.

In the case of optimization of actuators locations, [START_REF] Arbel | Controllability measures and actuator placement in oscillatory systems[END_REF], [START_REF] Hac | Sensor and actuator location in motion control of flexible structures[END_REF] and [START_REF] Devasia | Piezoelectric actuator design for vibration suppression : plavement and sizing[END_REF] proposed to maximize a controllability criteria using a measure of the gramian matrix. This approach seeks to ensure active damping of all needed modes by minimizing the control energy. In a second approach [START_REF] Wang | A controllability index for optimal design of piezoelectric actuators in vibration control of beam structures[END_REF] suggested to maximize the control forces transmitted by the actuators to the structure. A third optimization cost function is a linear quadratic optimal framework. [START_REF] Dhingra | Multiobjective design of actively controlled structures using a hybrid optimization method[END_REF], [START_REF] Kondoh | The positioning of sensors and actuators in the vibration control of flexible systems[END_REF], [START_REF] Bruant | A methodology for determination of piezoelectric actuator and sensor location on beam structures[END_REF] and [START_REF] Yang | Optimization of noncollocated sensor/actuator location and feedback gain in control systems[END_REF] proposed a quadratic cost function taking into account the measurement error and control energy. In the case of sensors locations, the most usual performance uses the energy of the state output so as to maximize the information given by sensors. [START_REF] Baruh | Placement of sensors and actuators in structural control[END_REF] and [START_REF] Hac | Sensor and actuator location in motion control of flexible structures[END_REF] rather proposed to maximize measures of the gramian observability matrix in order to have optima without dependence on initial conditions.

In this paper, a methodology is presented to find the actuators and sensors location to increase control quality. In order to simplify the optimization problem, it was decided to search independently the optimal location of actuators and that of sensors.

In applications, the external disturbances are not allways wellknown. In addition, many regulator methods can be used. Then, it seems to be better to define an optimization methodology independent of these external work parameters. From these reasons, we propose here to obtain optimal actuators location by maximizing a modified criteria which can result from the two first different approaches. From a given number of eigenmodes, our objective is to find the actuators locations such that the actuators control is the largest possible using less possible electrical energy. In the same way, as the sensors have to inform the control system about the vibrationnal state of the structure for each eigenmode, we present a modified criteria in order to maintain the system output as large as possible.

The spillover effects appear as a strong problem of active control implementation on real structures. Very few papers take into account the residual modes in their developpements [START_REF] Gaudiller | Active control of flexible structures using a minimum number of components[END_REF], [START_REF] Collet | Shape optimization of piezoelectric sensors dealing with spill-over instabilities[END_REF]. Here, each criteria is presented considering them, in order to limit the actuation or observation by active elements on a certain subset of modes.

In section 2, we point out the active vibration control equations. In section 3 and 4, we present the optimization criteria used for sensors and actuators location. Comparisons with classical criteria are shown in section 5, illustrated for a simply supported plate with piezoelectric active elements. Discussion is developed about the optimal number of actuators and sensors, especially when residual modes are taken into account in the cost functions. Results are shown for a simply supported beam with point sensors.

EQUATIONS OF ACTIVE CONTROL

Consider a flexible structure with N a actuators and N s velocity sensors (all further developpements can be applied for displacement sensors). The motion's equations and the sensors' equations of the system in modal coordinates can be written as follows, when the N first eigenmodes are considered and when taking into account the N R residual eigenmodes:

αi + 2ζ i ω i αi + ω 2 i α i = Na l=1 b il Φ l i = 1..N (1) αi + 2ζ R i ω R i αi + (ω R i ) 2 α i = Na l=1 b R il Φ l i = 1..N R (2) y j = N l=1 c jl αl + N R l=1 c R jl αl j = 1.
.N s α i , αi and αi represent modal displacement, velocity and acceleration, ω i and ζ i are the natural frequency and damping ratio of the ith mode, and ω R i and ζ R i those of the residual modes; b il Φ l is the ith modal component of the control force due to Φ l applied to the actuator l, b R il Φ l is the ith residual modal component of the force appearing with the actuation of the actuator l. y j is the quantity measured from the jth sensor. c jl is the sensing constant of the jth sensor due to the motion of the lth mode and c R jl those due to the motion of the lth residual mode. b il , b R il , c jl , c R jl depend respectively of the lth actuator location and jth sensor location.

These equations can be written in a state-space form, using the state vector {x}(size

(N + N R ) + (N + N R )): {x} = {ω i α i αi } T (3) d dt {x} = [A]{x} + [B]{Φ} (4) {y} = [C]{x} (5) [A] (2N +2N R ,2N +2N R ) , [B] (2N +2N R ,Na) and [C] (Ns,2N +2N R )
are the state, control and output matrices given by:

[A] =     [0] [0] [ω i ] [ 0 ] [0] [0] [0] [ω R i ] -[ω i ] [0] -[2ζ i ω i ] [0] [0] -[ω R i ] [0] -[2ζ R i ω R i ]     (6) [B] =     [0] [b il ] [0] [b R il ]     (7) [C] = [0] [0] [c jl ] [c R jl ] (8)
Φ is the actuation vector: in the case of piezoelectric actuators, it contains electrical potentials.

From equation ( 4), several automatic tools can be used to actively control vibrations [START_REF] Preumont | Vibration control of active structures[END_REF] by applying a well chosen {Φ}. The actuation must be defined to bend the N first eigenmodes, but it also actuates the residual modes exciting them. This effect is called spillover. In fact, the best case would be having:

b il >> 0 ∀l ∈ 1, ...N a and ∀i ∈ 1...N b R il = 0 ∀l ∈ 1, ...N a and ∀i ∈ 1...N R ,
ensuring actuation (or controllability) for the N first modes and no influence (or non controllability) of the N R residual modes. This configuration should be approached by optimizing the actuators locations.

In the same way, the vibrational information given by the output sensors to the control system contains motions on residual modes. As the control system neglects these residual modes, the received information does not correspond to the needed information. The best case would be having:

c jl >> 0 ∀j ∈ 1, ...N s and ∀l ∈ 1...N c R jl = 0 ∀j ∈ 1, ...N s and ∀l ∈ 1...N R ,
in order to have each mode (l ∈ 1...N ) well observed (good observability) and each residual mode not observed (non observability). As before, this configuration should be approached by optimizing the sensors locations.

Hence, before setting up the regulator and observer system, the active elements' locations have to be defined.

OPTIMAL LOCATION OF ACTUATORS

Many studies have been developed on optimal location of actuators. Two usual approaches are used before setting up the regulator: minimizing the control energy required by maximizing a measure of the controllability gramian matrix, and maximizing in a certain way the control force.

In this work, we suggest considering these two approaches: finding the actuators' locations so that the actuators' forces are the largest possible (2nd approach) using electrical energy (1st approach) as possible.

In the next subsections, we briefly describe these approaches and explain our modified criteria.

Minimization of the control energy

This method is detailed in [START_REF] Hac | Sensor and actuator location in motion control of flexible structures[END_REF]. The objective here is to find actuators locations which minimize the control energy required to bring the modal system (considering the N first eigenmodes) to a desired state {x T } after some time T :

J = min {Φ} T 0 {Φ} T {Φ}dt (9)
The optimal solution gives the following optimal control energy :

J = [e [A]T {x 0 } -{x T }] T W -1 (T )[e [A]T {x 0 } -{x T } (10)
where W (T ) is the controllability gramian matrix defined by :

W (T ) = T 0 e [A]t [B][B] T e [A] T t dt (11)
Hence, as the actuators' effects are only in W (T ) from [B], minimizing J with respect to the actuators locations consists in minimizing W -1 (T ) or maximizing a measure of the controllability gramian matrix. [START_REF] Hac | Sensor and actuator location in motion control of flexible structures[END_REF] has shown that instead of using W (T ), a steady state W c can be considered to eliminate the dependency of the solution T . W c tends to a diagonal form with

(W c ) ii = (W c ) i+N +N R ,i+N +N R = (W U c ) ii (12) = Na j=1 b 2 ij 4ζ i ω i = 1 4ζ i ω i Na j=1 b 2 ij (13) (W c ) i+N,i+N = (W c ) i+2N +N R ,i+2N +N R = (W R c ) ii (14) = Na j=1 (b R ij ) 2 4ζ R i ω R i i = 1, ...N R (15) = 1 4ζ R i ω R i Na j=1 (b R ij ) 2 i = 1, ...N R (16) (W U c ) ii and (W R c )
ii equal to the energy transmitted from the actuators to the structure for the ith used and residual eigenmode.

Hence if one eigenvalue of (W U c ) ii is small there is an eigenmode that is difficult to control: there is no controllability for the system. Similarly, if one eigenvalue

(W c ) R
ii corresponding to a residual mode is high, the induced spillover effect can be important. In order to avoid these cases, several criteria can be maximized, as [START_REF] Arbel | Controllability measures and actuator placement in oscillatory systems[END_REF], [START_REF] Hac | Sensor and actuator location in motion control of flexible structures[END_REF], [START_REF] Collet | Shape optimization of piezoelectric sensors dealing with spill-over instabilities[END_REF]:

det(W U c ) -γdet(W R c ) (17) trace(W U c ) * (det(W U c )) 1/2N -γtrace(W R c ) * (det(W R c )) 1/2N (18) min i=1,N (W U c ) ii -γ max i=1,N R (W R c ) ii ( 19
)
where γ is a weighting constant. These three criteria ensure global controllability of the system for the N first eigenmodes and try to minimize the global excitation of the N R residual modes.

Maximization of the control force

In this approach, detailed in [START_REF] Wang | A controllability index for optimal design of piezoelectric actuators in vibration control of beam structures[END_REF], the actuators locations are obtained by maximizing the global control force

{f c } = [B]{f } from {f c } T {f c } = {Φ} T [B] T [B]
{Φ} (here residual modes are neglected).

Using singular value analysis, [B] can be written as:

[B] = [M ][S][N ] T where [M ] T [M ] = [I] = [N ] T [N ], [S] =     σ 1 0 0 0 σ i 0 0 0 σ Na 0 0 0     (20)
and

{f c } T {f c } = {Φ} T [N ][S] T [S][N ] T {Φ} (21) or ||{f c }|| 2 = ||{Φ}|| 2 ||S|| 2 .
σ 2 i ( 22
)
taking into account globally all modes. In this equation, σ i is refered to as the ith degree of controllability of the system. Its magnitude is a function of the location and size of the piezoelectric actuators.

When residual modes cannot be neglected, a similar criteria can be obtained, consisting in differentiating the N first modes from the N R residual modes and maximizing :

Na i=1 σ 2 i -γ Na i=1 (σ R i ) 2 (23)
where σ R i are the components of [S] R corresponding to b R , and γ is a weighting constant.

The modified criteria

The use of these criteria shows that the eigenmodes are taken into account globally. As our objective is to control each considered modes, without exciting residual modes, by transmitting a maximum control force with a minimum electrical energy, we suggest here modifying these criteria.

• First of all, instead of maximizing a global norm of W c which means minimizing the electrical energy, we suggest maximizing each diagonal term of W c (the energy transmitted from the actuators to the structure for each eigenmode (1...N ))and minimizing the N R other diagonal term of W c :

maximize {(W U c ) 11 , ..., (W U c ) NN } (24) minimize {(W R c ) 11 , ..., (W R c ) N R N R } (25)
This optimization problem has N a optimization parameters (the locations of the N a actuators) and N + N R optimization functions. It can be simplified by considering the new optimization criteria:

maximize { min i=1,...N (W U c ) ii -γ max i=1,...N R (W R c ) ii } (26)
where γ is a weighting constant. This criteria, with γ = 0 has been proposed by [START_REF] Arbel | Controllability measures and actuator placement in oscillatory systems[END_REF], but as the components of W c have not the same range, solving this problem can induce the study of particular modes instead of each of them, and then the obtained locations are not optimal.

Consequently, we decide to establish homogeneity between each term (W c ) ii by dividing each of them by its maximal value obtained if the ith mode is the specific mode to be controlled. This maximal value is the maximal energy which can be transmitted from the actuators for the ith eigenmode.

Hence, the optimization problem becomes, by using the homogeneous components: maximizing

{ min i=1,N (W U c ) ii max locations (W U c ) ii -γ max i=1,N R (W R c ) ii max locations (W R c ) ii } (27) = min i=1,N Na j=1 b 2 ij max locations Na j=1 b 2 ij -γ max i=1,N R Na j=1 (b R ij ) 2 max locations Na j=1 (b R ij ) 2 (28) and, ∀i = 1, ...N + N R 0 ≤ (W c ) ii max locations (W c ) ii ≤ 1 (29)
• In the same way, instead of maximizing the norm of the global control force for the considered modes, we suggest to maximize the applied force for each mode (1, ...N ) independently of {Φ} and to minimize the applied force for residual modes. For the ith mode, the modal equation can be written in the state space form :

ω i αi αi = 0 ω i -ω i -2ζ i ω i ω i α i αi + {0} {B i } {Φ} (30) and ||{f c }|| 2 = ||{Φ}|| 2 .||S i || 2 where S 2
i has only a non zero term :

σ2 i = Na j=1 b 2 ij (31)
Then, the optimization problem is:

maximize {||S 1 || 2 , ...||S N || 2 } = {σ 2 1 , ...σ 2 N } (32) and minimize {||S R 1 || 2 , ...||S R N || 2 } = {σ 2 1 , ...σ 2 N R } (33)
From the same arguments than before, this problem is thus transformed: maximizing

{ min i=1,N σ2 i max locations σ2 i -γ max i=1,N R σ2 i max locations σ2 i } (34) = min i=1,N Na j=1 b 2 ij max locations Na j=1 b 2 ij -γ max i=1,N R Na j=1 (b R ij ) 2 max locations Na j=1 (b R ij ) 2
(35) By comparing equations ( 28) and ( 35) it follows that the two different approaches give the same final criteria. Its greatest advantage is that all modes are studied with the same range. Residual modes are not neglected; their influence on the structure dynamic can be more or less important using γ. Furthermore, the expression inside ( 28) and ( 35) has a physical meaning: it is the mechanical energy transmitted for the ith mode divided by the maximal mechanical energy that could be received.

OPTIMAL LOCATION OF SENSORS

The optimal locations of sensors are determined in the same way than the optimal locations of actuators.

Maximization of the system output

When the system is released from the initial state {x(0)} = {x 0 }, as when it is subjected to a persistent disturbance, [START_REF] Hac | Sensor and actuator location in motion control of flexible structures[END_REF] has shown that maximizing the system output ∞ 0 {y} T {y}dt (for desired modes) yields maximizing the gramian observability matrix defined by :

W o = ∞ 0 e [A] T t [C] T [C]e [A]t dt (36)
and W o tends to a diagonal form (for displacement sensors

): i = 1, ...N (W o ) ii = (W o ) i+N,i+N = Ns j=1 c 2 ji 4ζ i ω i = 1 4ζ i ω i Ns j=1 c 2 ji (37)
Each diagonal term (W o ) ii corresponds to the maximization of the output energy J i for the ith mode, obtained if we consider the state equation reduced to the ith mode:

ω i αi αi = 0 ω i -ω i -2ζ i ω i ω i α i αi + [0] +[b il ] {Φ} (38) {y} = [C i ] 0 αi J = ∞ 0 {y} T {y}dt (39)
Consequently, if the ith eigenvalue of W o is small, it means that the ith mode will not be well observed. To avoid this case, [START_REF] Hac | Sensor and actuator location in motion control of flexible structures[END_REF], [START_REF] Baruh | Placement of sensors and actuators in structural control[END_REF] suggest to find the sensors locations by maximizing a measure of W o , like:

trace(W o ) * (det(W o ) 1/2N ) min i (W o ) ii (40)
These criteria take into account the eigenmodes globally.

They can be used with the residual modes as for actuators, in order to maximize the given information for the N first modes and minimize the information about the N R modes (in order to minimize the spillover effect).

The modified criteria

As we want to have a convenient information for the N first eigenmodes we suggest maximizing the output energy J i for each mode i (each diagonal term of W o denoted (W U o ) ii ) instead of maximizing the global system output, and minimizing the output energy of each residual mode

(called (W R o ) ii ): maximize {(W U o ) 11 , ..., (W U o ) NN } (41) minimize {(W R o ) 11 , ..., (W R o ) N R N R } (42)
This optimization problem has N s optimization parameters (the locations of the N s sensors) and N optimization functions. It can be simplified by considering the new optimization criteria:

maximize { min i=1,...N (W U o ) ii -γ max i=1,...N R (W R o ) ii } (43)
where γ is a weighting constant. As the components of W o have not the same range, solving this problem can induce the study of particular modes instead of each of them, and then the obtained locations will not be optimal. Hence, we decide to establish homogeneity between each term (W o ) ii by dividing each of them by its maximal value obtained if the ith mode is the specific mode to be measured. The optimization problem becomes to maximize the minimal value of the homogeneous components:

min i=1,N (W U o ) ii max locations (W U o ) ii -γ max i=1,N R (W R o ) ii max locations (W R o ) ii (44) = min
with ∀i = 1, ...N + N R 0 ≤ (W o ) ii max locations (W o ) ii ≤ 1 (46)
Then all modes are studied with the same range. The first quotient inside (44) represents the output energy measured for the ith mode divided by the maximal output energy which could be measured for the ith mode by the sensors.

APPLICATIONS

In this section, we aim to compare the different usual criteria of controllability and observability with one modified criteria. In order to make comparisons, we define two degrees.

A degree of controllability and a degree of observability

We define the following degree of controllability for each mode (i = 1, ..N + N R ):

DEGCi = Na j=1 b 2 ij max one actuator b 2 i1 (47)
It equals to the energy transmitted from the actuators to the structure for the ith mode divided by the maximal value energy obtained if the ith mode is optimally controlled by one actuator (if one actuator is located in order to optimally control this mode). When the N first degrees are over 100 % it means that each mode (i = 1, ...N ) is better controlled than when it is specifically controled by an actuator.

Inversely, for residual modes, the objective is to minimize DEGCi the higher possible.

In the same way, the following degree of observability is defined for each mode i:

DEGOi = Ns j=1 c 2 ji max one sensor c 2 1i (48)
It equals to the output energy measured from the sensors for the ith mode divided by the maximal value of output energy obtained if the ith mode is optimaly measured by one sensor (if one sensor is located to optimally measure this mode). When the N degrees are over 100%, each mode is better measured than when it is specifically measured by a sensor. Of course, the opposite is searched for residual modes.

In the next subsections, two applications are presented. The first one consists in comparing several criteria for actuators locations. The second one consists in studying the influence of residual modes in the optimization procedure.

First application: optimal location of piezoelectric actuators on a simply supported plate

A simply supported plate is considered, equiped with N a piezoelectric actuators to locate. In order to consider only pure bending motion, each actuator is made up of a pair of piezoelectric materials attached symetrically. They are assumed to be perfectly bonded to the surface of the plate, and their thickness is assumed to be small compared to the plate thickness. The geometrical and mechanical properties of the system are detailed in tables 1 and 2. The plate is divided in 400 rectangular elements which correspond to the possible locations of actuators. As the structure has two symmetry axes, we only consider its quarter. In this study, the first seven eigenmodes are used and all residual modes are neglected.

Simulations are presented for piezolectric actuators. Exactly the same results would be obtained for piezoelectric sensors measuring velocity [START_REF] Bruant | Optimal location and number of sensors and actuators in active control[END_REF].

The four following criteria are compared using DEGCi:

min i (W c ) ii max locations (W c ) ii called M IN L/LM AX (49) trace(W c ) * (det(W c )) 1/2N called T RACE * DET (50) min i (W c ) ii called L (51) Na i=1 σ 2 i called W AN G (52)
Results for different actuators numbers are shown figures 1 to 3. Figure 1.a represents, for each configuration, the smallest value DEGCi (i = 1...N ). The others plots represent the degree of controllability of each eigenmode.

By comparing results, the smallest value DEGi has the highest values for the criteria MINL/LMAX. This is because other criteria use matrice (gramian or singular) which have not homogeneous components. Consequently, some eigenmodes are favourite instead of others. For example, for TRACE*DET the 6th mode and also modes 3,4 and 5 have a degree of controllability higher than the modes 1, 2 and 7. This is also the case for L: the 1st mode is very well controlled instead of modes 4 and 6 which are very badly controlled even if seven actuators are used. In the optimization of WANG, the worst controlled mode is the fourth. Hence, these criteria ensure avoiding non controllability area, but do not ensure homogeneous con-trollability for considered modes. For MINL/LMAX, the variations of the degrees DEGCi are homogeneous, and they are more than 300% using seven actuators.

About the number of actuators

The controllability degree gives us information about the quality of the control for each mode. It can also be used to define the number of actuators needed. When the smallest value DEGCi, in figure 3, is over 100% it means that each mode is better controlled than when it is specifically controlled by an actuator. Consequently, we can choose the value 100% as criteria for the optimal number of actuators. From this point of view, the number of actuators needed is for MINL/LMAX 3, even though TRACE*DET and WANG need 4 actuators and L needs 7 actuators.

Furthermore, the smallest value DEGCi is over 200% for MINL/LMAX using 5 actuators. Then, when a robust control system which ensures good controlability (over 100% for each mode) is necessary even if one actuator is failing, we will use this configuration. For the other criteria, it is not possible to use only 5 actuators located optimaly.

In this application, we have not considered the residual modes. In fact, for better optimal configuration of actuators, it should be necessary to take into account them. In our modal approach, it would consist in adding a sufficient number of residual modes. The use of the modal basis gives us a phyiscal meaning of gramian matrix.

In the next subsection, the residual modes are used to locate the sensors.

Second application: optimal location of point sensors on a slender uniform beam

As it was explained previoulsy, residual modes can be dangerous in active control, leading to instabilities. The used output sensors contain informations on motion on all modes instead of only the needed motion of the N first modes. In the same way, the actuators actuate all modes of the real structure, and consequently can excite residual modes.

In order to illustrate the influence of residual modes in the optimization procedure, we consider here a simply supported uniform beam described in [START_REF] Hac | Sensor and actuator location in motion control of flexible structures[END_REF]. In this case, an analytical analysis can be easily developed, and as the structure is symmetric, only the first middle of the beam is studied and divided in 50 possible point sensors locations. The observability gramian is equal to:

(W o ) ii = 1 2ρAlζ i ω i ( NS j=1 sin 2 iπx j l )
where ρA is the mass density, x j the location of the jth sensor and l = 1 m , the beam length. We want to precisely measure the displacement of the first six modes while avoiding residual modes: modes 7, 8, 9 and 10.

The results are also available for point actuators.

The considered criteria is defined in (44). it was used for different sensors quantity, with γ = 0 and γ = 1. Results are shown in tables 3 and 4. For each case, the smallest degree of observability for needed modes (SDEGOi) and the highest degree of observability for residual modes (HDEGOi R ) are presented. To avoid measure of residual modes, HDEGOi R has to be minimal, and to have good information about the first 6 modes, SDEGOi has to be maximal, with HDEGOi R less than SDEGOi. First of all, it can be seen from these results that forgetting residual modes in the optimization procedure (γ = 0) can be dangerous: except when there are 3 sensors, the maximal value HDEGOi R is over the minimal value SDEGOi. It means that one residual mode is better observed than a needed mode. Consequently, the sensors information can be strongly disturbed, inducing spillover effect.

When residual modes are taken into account in the optimization criteria, HDEGOi R is smallest than SDEGOi except when there is one sensor. Comparing results for γ = 0 and for γ = 1 show that the value of each degree decreases between γ = 0 and γ = 1. When γ = 1, the optimization procedure has to adjust the sensors locations obtained for γ = 0 in order to avoid measurement of residual modes. These kind of simulations can also help us to define the number of needed sensors. A good compromise between needed modes and residual modes seems to be obtained by using 3 sensors: each needed modes is better measured than each residual modes (SDEGOi > HDEGOi R ). In an other part, the degrees of observability of the six first modes are over 100% and those of each residual modes is less than 100%.

CONCLUSION

In this paper, the problem of sensors and actuators locations is considered. We propose for each optimization problem a modified optimization criteria derived from usual approaches, ensuring good observability and controllability of each mode of the structure, and taking into account residual modes which have to be less observable and controllable.

Their efficiency are shown by comparing them with classical criteria, especially when the number of active elements is varying. The results show how important is the choice of the optimization criteria for the effectiveness of the control. They also give us some help to define the optimal sensors and actuators number. The great advantages of these criteria are their computational simplicity, their non-dependance with the external disturbances and with the applied control law. This methodology was used here considering actuators and sensors independently. It could be studied for collocated active elements. In the same way, same developpements could be done for active elements shape.

  Figure 1.

  Figure 3.

  Thus, maximizing this norm independently on the input {Φ} induces maximizing ||S|| 2 . Wang suggests to maximize

	Na
	i=1

table 1 :

 1 characteristics of the plate

	Length (m)	0.38
	width (m)	0.3
	thickness (m) 0.0031752 E (Pa) 20710 9
	ν	0.292
	ρ (kg/m 3 )	7870
	table 2: characteristics of the PZT 5A
	Length (m)	0.019
	width (m)	0.015
	thickness (m)	0.001
	33 (F/m) e 31 (C.m -2 ) e 32 (C.m -2 )	1.510 -8 -7.209 -7.209

table 3

 3

		: results for γ = 0
	sensors numbers SDEGOi HDEGOi R
	1	18%	96 %
	2	93%	100 %
	3	167%	144 %
	4	221%	325 %
	table 4: influence of the residual modes: γ = 1
	sensors number SDEGOi HDEGOi R
	1	13,5 %	34 %
	2	7 5 %	7 4 %
	3	116 %	82 %
	4	191 %	136 %