In this text we will formalise a novel solution, the Bayesian Volumetric Map (BVM), as a framework for a metric, short-term, egocentric spatial memory for multimodal perception of 3D structure and motion. This solution will enable the implementation of top-down mechanisms of attention guidance of perception towards areas of high entropy/uncertainty, so as to promote active exploration of the environment by the robotic perceptual system. In the process, we will to try address the inherent challenges of visual, auditory and vestibular sensor fusion through the BVM. In fact, it is our belief that perceptual systems are unable to yield truly useful descriptions of their environment without resorting to a temporal series of sensory fusion processed on a short-term memory such as the BVM.

I. INTRODUCTION

Perception has been regarded as a computational process of unconscious, probabilistic inference. Aided by developments in statistics and artificial intelligence, researchers have begun to apply the concepts of probability theory rigorously to problems in biological perception and action. One striking observation from this work is the myriad ways in which human observers behave as near-optimal Bayesian observers, which has fundamental implications for neuroscience, particularly in how we conceive of neural computations and the nature of neural representations of perceptual variables [START_REF] Knill | The Bayesian brain: the role of uncertainty in neural coding and computation[END_REF].

Consider the following scenario -an observer is presented with a non-static 3D scene containing several moving entities, probably generating some kind of sound: how does this observer perceive the 3D structure of all entities in the scene and the 3D trajectory and velocity of moving objects, given the ambiguities and conflicts inherent to the perceptual process? Given these considerations, the research presented on this text regards a Bayesian framework for artificial perception models.

In this text we will formalise a novel framework for a metric, short-term, egocentric spatial memory that will enable the implementation of top-down mechanisms of attention guidance of perception towards areas of high uncertainty, so as to promote active exploration of the environment by the robotic perceptual system.

To support our research work, an artificial multimodal perception system (IMPEP -Integrated Multimodal Perception Experimental Platform) has been constructed at the This publication has been supported by EC-contract number FP6-IST-027140, Action line: Cognitive Systems. The contents of this text reflect only the author's views. The European Community is not liable for any use that may be made of the information contained herein. ISR/FCT-UC consisting of a stereovision, binaural and inertial measuring unit (IMU) setup mounted on a motorised head, with gaze control capabilities for image stabilisation and perceptual attention purposes -see Fig. 1.

II. BACKGROUND AND DEFINITIONS

The perceptual brain is known to be divided into two separate, albeit interdependent, pathways of sensory processing: the ventral pathway, popularly denominated as the "what" pathway, which is concerned primarily with perceptual classification and recognition tasks, and the dorsal pathway, popularly known as the "where" pathway, which is dedicated to fast processing of sensory information, with the sole purpose of yielding spatial representations (e.g. positioning, structure and motion), whatever the nature of the entity it is analysing. Given the perceptual problem exposed earlier on, the latter is of particular interest to our work. These spatial representations are believed to be metric and egocentric in lower-level areas of the dorsal pathway so as to promote fast and accurate interaction with the surrounding environment.

Given these facts, the framework for spatial representation that will be presented in the rest of this section satisfies the following criteria: it is egocentric and metric in nature; it allows for the representation of dynamical spatial occupation of the environment, avoiding any need for any assumptions on the nature of those objects, or in other words, for data association. Data association is thus effectively postponed to higher-level processing. Metric maps are very intuitive, yield a rigorous model of the environment and help to register measurements taken from different locations. Grid-based maps are the most popular metric maps in mobile robotics applications. One of the most popular grid-based maps is the occupancy grid, which is a discretised random field where the probability of occupancy of each cell is kept, and the probability values of occupancy of all cells are considered independent between each other [START_REF] Elfes | Using occupancy grids for mobile robot perception and navigation[END_REF]. The absence of an object based representation permits the ease of fusing low level descriptive sensory information onto the grids without necessarily implicating data association.

In our specific application domain, common occupancy grid configurations which assume regularly partitioned euclidean space to build the cell lattice are not appropriate. Hence, we chose a log-spherical coordinate system spatial configuration (see Fig. 2) for our Bayesian Volumetric Maps (BVM), thus promoting an egocentric trait and yielding more precision for objects closer to the observer, which seems to agree with biological perception.

The BVM is primarily defined by its range of azimuth and elevation angles, and by its maximum reach in distance ρ Max , which in turn determines its log-distance base through b = a log a (ρ Max -ρ Min ) N

, ∀a ∈ R, where ρ Min defines the egocentric gap, for a given number of partitions N , chosen according to application requirements. The BVM space is therefore effectively defined by

Y ≡ ] log b ρ Min ; log b ρ Max ] × ]θ Min ; θ Max ] × ]φ Min ; φ Max ] (1)
In practice, the BVM is parametrised so as to cover the full angular range for azimuth and elevation.

Each BVM cell is defined by two limiting log-distances, log b ρ min and log b ρ max , two limiting azimuth angles, θ min and θ max , and two limiting elevation angles, φ min and φ max , through:

Y ⊃ C ≡ ] log b ρ min ; log b ρ max ] × ]θ min ; θ max ] × ]φ min ; φ max ]
(2) where constant values for log-distance base b, and angular ranges ∆θ = θ max -θ min and ∆φ = φ max -φ min , chosen according to application resolution requirements, ensure BVM grid regularity. Finally, each BVM cell is formally indexed by the coordinates of its far corner, defined as C = (log b ρ max , θ max , φ max ).

More recently, Coué et al. [START_REF] Coué | Bayesian occupancy filtering for multitarget tracking: an automotive application[END_REF] and Tay et al. [START_REF] Tay | An efficient formulation of the bayesian occupation filter for target tracking in dynamic environments[END_REF] expanded on the occupancy grid by explicitly introducing Bayesian filtering. A two-step mechanism estimates, at each time step, the state of the occupancy grid by combining a prediction step (history) and an estimation step (incorporating new measurements). This approach is derived from the Bayesian filtering approach [START_REF] Jazwinsky | Stochastic Processes and Filtering Theory[END_REF], and is thus named the Bayesian Occupancy Filter (BOF).

To compute the probability distributions for the current states of each cell, the Bayesian Program (BP) formalism, as first defined by Lebeltel [START_REF] Lebeltel | Programmation Bayésienne des Robots[END_REF], will be used throughout this text.

III. BAYESIAN VOLUMETRIC MAPS FOR MULTIMODAL PERCEPTION

A. Sensor fusion advantages and challenges

The use of more than one sensor promotes a robustness increase on the observation and characterisation of a physical phenomenon. In fact, using different types of sensors allows for the dilution of each sensor's individual weaknesses through the use of the strengths of the remainder.

There is evidence that humans fuse perceptual cue information following mainly two general strategies [START_REF] Ernst | Merging the senses into a robust percept[END_REF]: combination, that expresses interactions between sensory signals that are not redundant, and integration, that expresses interactions between sensory signals that are redundant. Combination has the purpose of maximising information coming from different cues, whilst the goal of integration is to minimise variance in the sensory estimate to increase its reliability. For several estimates resulting from combination to be integrated into a single estimate, they must be in the same units and referred to the same coordinate system, and hence must undergo a process called promotion [START_REF] Ernst | Merging the senses into a robust percept[END_REF].

We will try to explicitly or implicitly address each of the challenges of sensor fusion as described in [START_REF] Ernst | Merging the senses into a robust percept[END_REF] using the BVM, for vision, audition and vestibular sensing. It is our belief that perceptual systems are unable to yield useful descriptions of their environment without resorting to a temporal series of sensory fusion processed on a shortterm memory such as the BVM. We propose to use vestibular sensing as ancillary information to promote visual and auditory sensing to satisfy the requirements for integration, enumerated above.

B. Using Bayesian filtering for visuoauditory integration

The Bayesian Program presented in Fig. 3 is based on the solution presented by Tay et al. [START_REF] Tay | An efficient formulation of the bayesian occupation filter for target tracking in dynamic environments[END_REF], adapted so as to conform to the BVM egocentric, three-dimensional and log-spherical nature.

The estimation of the joint state of occupancy and velocity of a cell is answered through Bayesian inference on the decomposition equation given in Fig. 3. This inference effectively leads to the Bayesian filtering formulation as used in the BOF grids -see Fig. 4. In this context, prediction propagates cell occupancy probabilities for each velocity and cell in the grid -P (O C V C |C). During estimation, P (O C V C |C) is updated by taking into account the observations yielded by the sensors Questions: from the Bayesian filter estimation will then be used for the prediction step in the next iteration.

S i=1 P (Z i |V C O C C) to obtain the final state estimate P (O C V C |Z 1 • • • Z S C). The result Program                                                                                                                                  Description                                                                                                                  Specification                                                                                                       
P (C A C O C O -1 C V C Z1 • • • Z S ) = P (A C )P (V C |A C )P (C|V C A C )P (O -1 C |A C )P (O C |O -1 C ) S i=1 P (Zi|V C O C C)
P (Oc Vc|z1 • • • z S c) → P (Oc|z1 • • • z S c) P (Vc|z1 • • • z S c)

C. Using the BVM for sensory combination of vision and audition with vestibular sensing

Consider the simplest case, where the sensors may only rotate around the egocentric origin and the whole perceptual system is not allowed to perform any translation. In this case, the vestibular sensor models will yield measurements of angular velocity and position, which can then be easily used to manipulate the BVM, which is, by definition, in spherical coordinates.

To maintain a head-centred coordinate system for the BVM, which would obviously shift in accordance to head turns, instead of rotating the whole map, the most effective solution is to perform the equivalent index shift. This process is described by redefining C: C ∈ Y indexes a cell in the BVM by its far corner, defined as C = (log b ρ max , θ max -

θ inertial , φ max -φ inertial ) ∈ C ⊂ Y.
This process obviously relies on the assumption that inertial precision on angular measurements is greater than the chosen resolution parameters for the BVM.

D. Dealing with sensory synchronisation

The BVM model presented earlier assumes that the state of a cell C, given by (O C , V C ), and the observation by any sensor i, given by Z i , correspond to the same time instant t.

In accordance with the wide multisensory integration temporal window theory for human perception reviewed in [START_REF] Spence | Multisensory integration: maintaining the perception of synchrony[END_REF], the BVM may be used safely to integrate auditory and vision measurements as soon they become available; preliminary tests using the BVM update model show that this, in fact, promotes an effect similar to the well-known temporal ventriloquism, given the inherent auditory measurement frequency as opposed to vision. Spatial ventriloquism, on the other hand, is implicitly ensured due to the inherent properties of the Bayesian formulation of visuoauditory integration (i.e. modality reliability expressed in terms of uncertainty). Promotion through vestibular sensing is also not a problem, since inertial readings are available at a much faster rate than visuoauditory perception.

IV. SENSOR MODELS

A. Vision sensor model

Our motivations suggest for the vision sensor model a tentative data structure analogous to neuronal population activity patterns to represent uncertainty in the form of probability distributions -a spatially organised 2D grid has each cell associated to a population code simulation, a set of probability values of a neuronal population encoding a probability distribution [START_REF] Pouget | Information processing with population codes[END_REF]. The stereovision algorithm used for visual depth sensing is an adaptation of the fast and simple coherence detection approach by Henkel [START_REF] Henkel | A Simple and Fast Neural Network Approach to Stereovision[END_REF], yielding an estimated disparity map δ(k, i) and a corresponding confidence map λ(k, i). For visual perception of occupancy, this stereovision sensor described can be decomposed into simpler linear (1D) depth ρ(k, i) measuring sensors per projection line/pixel (k, i), each oriented in space with spherical angles (θ(k, i), φ(k, i)).

This algorithm is then easily converted from its deterministic nature into a probabilistic implementation simulating the population code-type data structure. This results in probability distributions on sensor measurements made available as likelihood functions taken from sensor readings -soft evidence, or "Jeffrey's evidence" in reference to Jeffrey's rule [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF]; the relation between vision sensor measurements Z and the corresponding readings δ and λ is thus described by the calibrated expected value ρ( δ) and standard deviation σ ρ (λ) for each sensor.

We have decided to model these sensors in terms of their contribution to the estimation of cell occupancy in a similar fashion to the solution proposed by Yguel et al. [START_REF] Yguel | Efficient GPU-based Construction of Occupancy Grids Using several Laser Range-finders[END_REF].

In the spirit of Bayesian programming, we again start by defining the relevant variables:

• C ≡ (log b ρ max , θ max , φ max ), O C and Z have the same meaning as before. However, once a projection line (θ, φ), with θ min ≤ θ ≤ θ max ∧ φ min ≤ φ ≤ φ max , is established for a sensor, only log b ρ max varies throughout the respective line-of-sight, thus effectively indexing each cell. Therefore, by abuse of notation and in order to simplify references to cells in the line-of-sight, these will be referred to using the abstraction C ∈ N, 1 ≤ C ≤ N , where N = log b (ρ M ax -ρ M in ) denotes the total number of cells in the line-of-sight.

• G C ∈ G C ≡ O N -1 represents the state of all cells in the line-of-sight except for C. Each g C is, thus, an (N -1)-tuple of the form ([O 1 = o 1 ], • • • , [O c-1 = o c-1 ], [O c+1 = o c+1 ], • • • , [O N = o N ]) given a specific cell [C = c].
The following expression gives the decomposition of the joint distribution of the relevant variables according to Bayes' rule and dependency assumptions:

P (Z C O C G C ) = P (C)P (O C |C)P (G C |O C C)P (Z|G C O C C) (3)
The parametric form and semantics of each component of the joint decomposition are then as follows:

• P (C) and P (O C |C) represent a priori information on the environment. Note that

P (C)P (O C |C) is, in fact, formally equivalent to P (A C )P (V C |A)P (C|V C A C )P (O -1 |A)P (O|O -1 )
when considering scene dynamics. The probability of a cell being empty is

P Empty = P ([O C = 0]|C). • P (G C |O C C) ≡ P (G C |C) represents the probability
that, knowing a state of a cell, the whole line-of-sight is in a particular state [START_REF] Yguel | Efficient GPU-based Construction of Occupancy Grids Using several Laser Range-finders[END_REF].

• P (Z|G C O C C) is sensor-dependent but, in any case, for all (O C , G C ) ∈ O × G C ,
P k (Z) = L k (Z, µ ρ (k), σ ρ (k)), µ ρ (k) = ρ( δ) σ ρ (k) = 1 λ σ min (4)
with σ min and ρ( δ) taken from calibration, the former as the estimate of the smallest error in depth yielded by the stereovision system and the latter from the intrinsic camera geometry. The likelihood function constitutes, in fact, the elementary sensor model as defined above for each vision sensor. Equation ( 4) only partially defines the resulting probability distribution by specifying the random variable over which it is defined and an expected value plus a standard deviation/variance -a full definition requires the choice of a type of distribution that best fits the noisy pdfs taken from the population code data structure. The traditional choice, mainly due to the central limit theorem, favours normal distributions N (Z, µ ρ (k), σ ρ (k)). Considering what happens in the mammalian brain, this choice appears to be naturally justified -biological population codes often yield bellshaped distributions around a preferred reading [START_REF] Treue | Seeing multiple directions of motion -physiology and psychophysics[END_REF], [START_REF] Born | Structure and Function of Visual Area MT[END_REF], [START_REF] Knill | The Bayesian brain: the role of uncertainty in neural coding and computation[END_REF], [START_REF] Pouget | Information processing with population codes[END_REF].

However, the fact that depth sensors always yield positive readings may be contradicted by the circumstance that normal distributions assign non-zero probabilities to negative depth values; even worse, close to the origin (Z = 0) this distribution assigns a high probability to negative depth values! With this purpose, we have adapted Yguel et al.'s Gaussian elementary sensor model so as to additionally perform the transformation to distance log-space, as follows

P k ([Z = z]) =        ]-∞;1] N (µ(k -0.5), σ(σ ρ ))(u)du, z ∈ [0; 1] z +1 z N (µ(k -0.5), σ(σ ρ ))(u)du, z ∈ ]1; N ] ]N ;+∞] N (µ(k -0.5), σ(σ ρ ))(u)du, z = "No Detection" (5)
where µ(•) and σ(•) are the operators that perform the required spatial coordinate transformations, and k = µ ρ is assumed to be the log-space index of the only occupied cell significant distance to the origin would produce an error in the angular position, like when a test pilot is in a centrifuge. In this model, X, Y and Z refer to the three axes of the robotic vision head in egocentric coordinates. The orientation of the system in space is encoded using a rotation matrix Θ. Angular velocity of the head is encoded using the yaw y, pitch p and roll r conventions. Yaw rotations are rotations around the Z axis; pitch around the Y axis and roll around X. When a rotation consists of a combination of yaw, pitch and roll rotation, the three rotations are applied successively and in this order. The rotation update is given by

Θ t+δt = Θ t .R(δy, δp, δr) (10)
where R(y, p, r) = c(y).c(p) c(y).s(p).s(r)s(y).c(r) c(y).s(p).c(r) + s(y).s(r) s(y).c(p) c(y).c(r) + s(y).s(p).s(r) -c(y).s(r) + s(y).s(p).c(r) -s(p) c(p).s(r) c(p).c(r) and the instantaneous angular velocity is defined as:

Ω =   δy/δt δp/δt δr/δt  
The calibrated inertial sensors in the IMU provide direct egocentric measurements of body angular velocity and linear acceleration (including gravity G). Given the motion of the system, we can define the probability distribution of the sensory inputs. The gyros will measure Ω t with added Gaussian noise, i.e. Φ t = Ω t + η t Φ , where η t Φ is a three-dimensional vector, the elements of which follow independent Gaussian distributions with mean 0 and standard deviation σ Φ . The accelerometers will measure the gravito-inertial acceleration F with added Gaussian noise, i.e. Υ t = F t + η t Υ , where η t Υ is a three-dimensional vector, the elements of which follow independent Gaussian distributions with mean 0 and standard deviation σ Υ . F is the resultant acceleration due to linear acceleration and gravity. Given the geocentric body linear acceleration A and the system orientation Θ, we can compute F . Transforming to the egocentric frame of reference we have

F = Θ -1 .(G -A) (11)
The sensor data at time t is therefore defined by (Φ t , Υ t ), and the state of our system at time t by (Θ t , Ω t , A t ). Estimates for spherical angles (θ inertial , φ inertial ) are then easily derived from the pitch-roll-yaw configuration of Θ.

As suggested in [START_REF] Laurens | Bayesian processing of vestibular information[END_REF], even in the absence of any sensory information, motion estimates for which the rotational velocity and acceleration are low are more probable. This can be described in a simple way using a Gaussian distribution. Having N (x, µ, σ) = e -(x-µ) 2 /(2.σ 2 ) √ 2.π.σ 2 the probability distribution for acceleration is given by P (A t ) ∝ N (|A t |, 0, σ A ); similarly for angular velocity Ω we have P (Ω t ) ∝ N (|Ω t |, 0, σ Ω ).

V. CONCLUSIONS

We have formalised herewith a novel solution, the Bayesian Volumetric Map, a framework for a metric, shortterm, egocentric spatial memory for multimodal perception of 3D structure and motion. This solution allows: the estimation of 3D structure and local motion states through perceptual fusion, involving vision, audition and inertial sensing, effectively postponing data association to higherlevel perceptual processing; the implementation of top-down mechanisms of attention guidance of perception towards areas of high entropy/uncertainty, so as to promote active exploration of the environment by the robotic perceptual system.

Further details on the calibration and implementation of these models can be found at http://paloma.isr.uc. pt/~jfilipe/BayesianMultimodalPerception.

Fig. 1 .

 1 Fig. 1. View of the current version of the Integrated Multimodal PerceptionExperimental Platform (IMPEP), on the left. On the right, the IMPEP perceptual geometry is shown: {E} is the main reference frame for the IMPEP robotic head, representing the egocentric coordinate system;{C l,r } are the stereovision (respectively left and right) camera referentials; {M l,r } are the binaural system (respectively left and right) microphone referentials; and finally {I} is the inertial measuring unit's coordinate system.

  log

Fig. 2 .

 2 Fig. 2. Egocentric, log-spherical configuration of the Bayesian Volumetric Maps.

1 C

 1 Relevant variables: C ∈ Y: indexes a cell on the BVM; A C : identifier of the antecedents of cell C (stored as with C); Z1, • • • , Z S ∈ {"No Detection"} ∪ Z: independent measurements taken by S sensors; O C , O -: binary values describing the occupancy of cell C, for current and preceding instants, respectively; V C : velocity of cell C, discretised into n possible cases ∈ V ≡ {v1, • • • , vn}. Decomposition:

1 C 1 C

 11 Parametric forms:P (A C ): uniform;P (V C |A C ): histogram; P (C|V C A C ): Dirac, 1 iff c log b ρ = a log b ρ + v log b ρ δt, c θ = a θ + v θ δt and c φ = a φ + v φ δt (constantvelocity assumption); P (O -|A C ): probability of preceding state of occupancy given set of antecedents; P (O C |O -): defined through transition matrix T = 1- 1-, where represents the probability of non-constant velocity; P (Zi|V C O C C): direct measurement model for each sensor i, given by respective sub-BP. Identification: None.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Bayesian Program for the estimation of Bayesian Volumetric Map current cell state.

  the probability distribution over Z depends only on the first occupied cell. Knowing the position of the first occupied cell in the projection line, which will be denoted as [C = k], P (Z|G C O C [C = k]) gives the probability of a measurement if [C = k] would be the only occupied cell in the line-of-sight. This particular distribution over Z is called the elementary sensor model, denoted by P k (Z). Given the first occupied cell [C = k] on the line-of-sight, the likelihood functions yielded by the population code data structure become

Proceedings of the 2008 International Conference on Cognitive Systems University of Karlsruhe, Karlsruhe, Germany, April 2-4, 2008

in the line-of-sight, which represents the coordinate interval ]k -1; k].

The answer to the Bayesian Program question in order to determine the sensor model P (Z|O C C) for vision, which is in fact related to the decomposition of interest P (O C Z C) = P (C)P (O C |C)P (Z|O C C), is answered through Bayesian inference on the decomposition equation given in [START_REF] Coué | Bayesian occupancy filtering for multitarget tracking: an automotive application[END_REF]; the inference process will dilute the effect of the unknown probability distribution P (G C |O C C) through marginalisation over all possible states of G C . In other words, the resulting direct model for vision sensors is based solely on knowing which is the first occupied cell on the line-of-sight and its relative position to a given cell of interest C.

To correctly formalise the Bayesian inference process, a formal auxiliary definition with respective properties follow.

Definition 1: A k c ∈ G C is the set of all tuples for which the first occupied cell is

), which are empty, and the kth,

), which are empty, and the (k -1)th, (o k ), which is occupied. Then,

Empty (1 -P Empty ). It now becomes possible to determine P (Z|O C C) in order to express the desired joint distribution P (Z O C C). This process leads to four distinct possible cases, that will be described next.

In the case of detection given an occupied cell [C = c], the sensor measurement can only be due to the occupancy of this cell or a cell before it in terms of visibility.

Thus [START_REF] Yguel | Efficient GPU-based Construction of Occupancy Grids Using several Laser Range-finders[END_REF],

Equation ( 6) has two terms: the left term that represents the case where [C = c] is occupied and the right term that comes from the aggregation of all the remaining probabilities around the last possible cell that might produce a detection: [C = c] itself. The "No Detection" case ensures that the distribution is normalised.

In the case of no detection given an occupied cell [C = c], which would correspond most probably to the effects of occlusion from earlier cells,

In the case of a measurement from detection knowing that [C = c] is empty, where a erroneous detection is yielded by the sensor (the so-called false alarm),

There are three terms in the empty cell, from left to right, corresponding respectively to before the detection, after the detection and no detection at all. Again, the "No Detection" case ensures that the distribution is normalised.

In the case of no detection knowing that [C = c] is empty, which will either be due to a miss-detection or a completely empty line-of-sight corresponding to G ∅ ,

B. Audition sensor model

The audition sensor model used as a source of observations for BVM cell state updates is fully described in [START_REF] Pinho | A Bayesian Binaural System for 3D Sound-Source Localisation[END_REF].

C. Vestibular sensor model

To process the inertial data, we adapted the Bayesian model proposed by Laurens and Droulez [START_REF] Laurens | Bayesian processing of vestibular information[END_REF] for the human vestibular system. The aim is to provide an estimate for the current angular position and angular velocity of the system, that mimics human vestibular perception. Since we only consider the simplest case, where sensors may only rotate around the egocentric origin, the angular rotation measurements may be safely assumed to be independent; linear acceleration might have a centripetal component that depends on the distance to the origin, but since the model is only detecting gravity, only a sustained rotation and a


Cátia Pinho

João Filipe Ferreira

Pierre Bessière

Jorge Dias

A Bayesian Binaural System for 3D Sound-Source Localisation

In this text we present a Bayesian system of auditory localisation in distance, azimuth and elevation using binaural cues only. We describe its supporting sensor model and calibration procedure. The binaural system is also integrated in a spatial representation framework for multimodal perception of 3D structure and motion -the Bayesian Volumetric Map (BVM). This solution will enable the implementation of an active perception system with great potential in applications as diverse as social robots or even robotic navigation.

I. INTRODUCTION

Although vision might be the dominant sense in humans, we rely on hearing as our only panoramic, long-range sensory system. The ability not only to detect and identify a sound, but also to pinpoint swiftly and accurately the location of its source can bring substantial advantages. This applies equally to a predator stalking its prey in the wild [START_REF] King | The shape of ears to come: dynamic coding of auditory space[END_REF] and to robotic applications such as [START_REF] Okuno | Sound and Visual Tracking for Humanoid Robot[END_REF], [START_REF] Schauer | A Computational Model of Early Auditory-Visual Integration[END_REF], [START_REF] Andersson | Robot Phonotaxis with Dynamic Sound-source Localization[END_REF], [START_REF] Hörnstein | Sound localization for humanoid robots -building audio-motor maps based on the HRTF[END_REF], [START_REF] Koene | Gaze shift reflex in a humanoid active vision system[END_REF], [START_REF] Rodemann | Purely Auditory Online-Adaptation of Auditory-Motor Maps[END_REF]. Moreover, auditory stimulus localisation is also an important component driving attention and gaze shifts, especially when the target is not in sight [START_REF] Koene | Gaze shift reflex in a humanoid active vision system[END_REF].

In this text we present a Bayesian system of auditory localisation in distance, azimuth and elevation using binaural cues only. We describe its supporting sensor model and calibration procedure. The binaural system is also integrated in a spatial representation framework for multimodal perception of 3D structure and motion, the Bayesian Volumetric Map (BVM) -for more details, please refer to [START_REF] Ferreira | Bayesian Models for Multimodal Perception of 3D Structure and Motion[END_REF].

To support our research work, an artificial multimodal perception system (IMPEP -Integrated Multimodal Perception Experimental Platform) has been constructed at the ISR/FCT-UC consisting of a stereovision, binaural and inertial measuring unit (IMU) setup mounted on a motorised head, with gaze control capabilities for image stabilisation and perceptual attention purposes -see Fig. 1. This solution will enable the implementation of an active perception system with great potential in applications as diverse as social robots or even robotic navigation (Fig. 2).

The Bayesian Program (BP) formalism, as first defined by Lebeltel [START_REF] Lebeltel | Programmation Bayésienne des Robots[END_REF], will be used to define the sensor model presented herewith.

II. BAYESIAN BINAURAL SENSOR MODEL

The Bayesian binaural system presented herewith is composed of three distinct and consecutive processors (Fig. 3):

This publication has been supported by EC-contract number FP6-IST-027140, Action line: Cognitive Systems. The contents of this text reflect only the author's views. The European Community is not liable for any use that may be made of the information contained herein. the monaural cochlear unit, which processes the pair of monaural signals {x 1 , x 2 } coming from the binaural audio transducer system by simulating the human cochlea, so as to achieve a tonotopic representation (i.e. a frequency band decomposition) of the left and right audio streams; the binaural unit, which correlates these signals and consequently estimates the binaural cues and segments each sound-source; and, finally, the Bayesian 3D sound-source localisation unit, which applies a Bayesian sensor model so as to perform localisation of sound-sources in 3D space.

A. Cochlear and auditory periphery processing

The first stages of auditory processing consist of cochlear and auditory periphery processing, which produces what is called an auditory image model (AIM) [START_REF] Patterson | Time-domain modeling of peripheral auditory processing: A modular architecture and a software platform[END_REF]. The AIM processor implements a functional model of a cochlea that simulates the phase-locked activity that complex sounds produce in the auditory nerve. Spectral analysis, the first stage of the AIM, is performed by a bank of auditory filters which converts each digitised wave that composes the stereo signal into an array of filtered waves. This processing is done using gammatone filters [START_REF] De Boer | Synthetic whole-nerve action potentials for the cat[END_REF], [START_REF] Immerseel | Digital implementation of linear gammatone filters: Comparison of design methods[END_REF], linearly distributed along a frequency scale measured in equivalent rectangular bandwidths (ERBs), as defined by [START_REF] Patterson | SVOS Final Report: The auditory filterbank[END_REF] for simulating the cochlea, obtaining a model of basilar membrane motion (BMM) through frequency band decomposition.

The second stage of the AIM simulates the mechanical/neural transduction process performed by the inner haircells. It converts the BMM into a neural activity pattern (NAP), which is the AIM's representation of the afferent activity in the auditory nerve [START_REF] Patterson | Time-domain modeling of peripheral auditory processing: A modular architecture and a software platform[END_REF]. In this stage the envelopes of the signals are first compressed, and then subjected to halfwave rectification followed by a squaring and lowpass filtering, resulting in m stereo audio signal pairs corresponding to m frequency channels with respective frequency centre

f k c , x ′ 1 (n), x ′ 2 (n) f k c , k = 1 • • • m.

B. Binaural cue processing

Sound waves arising from a source on our left will arrive at the left ear first. This small, but perceptible, difference in arrival time (known as an ITD, interaural time difference) is an important localisation cue and is detected by the inferior colliculus in primates, which acts as a temporal correlation detector array, after the auditory signals have been processed by the cochlea. Similarly, for intensity, the far ear lies in the head's "sound shadow", giving rise to interaural level differences (ILDs) [START_REF] King | The shape of ears to come: dynamic coding of auditory space[END_REF], [START_REF] Kapralos | Auditory Perception and Spatial (3D) Auditory Systems[END_REF]. ITDs vary systematically with the angle of incidence of the sound wave relative to the interaural axis, and are virtually independent of frequency, representing the most important localisation cue for low frequency signals (< 1500 Hz in humans). ILDs are more complex than ITDs in that they vary much more with sound frequency. Low-frequency sounds easily travel around the head and produce negligible ILDs. ILD values produced at higher frequencies are larger, and are increasingly influenced by the filter properties of each external ear, which imposes peaks and notches on the sound spectrum reaching the eardrum. Instead of being centred on the interaural axis, cones of confusion associated with particular ILD values take a different shape for each sound frequency.

Moreover, when considering sound sources within 1 -2 meters of the listener, binaural cues alone can even be used to fully localise the source in 3D space (i.e. azimuth, elevation and distance). Iso-ITD surfaces form hollow cones of confusion with a specific thickness extending from each ear in a symmetrical configuration relatively to the medial plane. On the contrary, iso-ILD surfaces, which are spherical surfaces delimit hollow spherical volumes, symmetrically placed about the medial plane and centred on a point on the interaural axis [START_REF] Shinn-Cunningham | Tori of confusion: Binaural localization cues for sources within reach of a listener[END_REF]. Thus, for sources within 2 meters range, the intersection of the ILD and ITD volumes is a torus-shaped volume [START_REF] Shinn-Cunningham | Tori of confusion: Binaural localization cues for sources within reach of a listener[END_REF]. If the source is more than 2 meters away, the change in ILD with source position is too gradual to provide spatial information (at least for an acoustically transparent head), and the source can only be localised to a volume around the correct cone of confusion [START_REF] Shinn-Cunningham | Tori of confusion: Binaural localization cues for sources within reach of a listener[END_REF].

Given this background, we have decided to adapt the solution by Faller and Merimaa [START_REF] Faller | Source localization in complex listening situations: Selection of binaural cues based on interaural coherence[END_REF] to implement the binaural processor. Using this algorithm, interaural time difference and interaural level difference cues are only considered at time instants when only the direct sound of a specific source has nonnegligible energy in the critical band and, thus, when the evoked ITD and ILD represent the direction of that source (corresponding to the process involving the superior olivary complex (SOC) and the central nucleus of the inferior colliculus (ICc) in mammals). They show how to identify such time instants as a function of the interaural coherence (IC). The source localisation suggested by the selected ITD and ILD cues are shown to imply the results of a number of published psychophysical studies related to source localisation in the presence of distractors, as well as in precedence effect conditions [START_REF] Zurek | The precedence effect[END_REF]. This algorithm thus amplifies the signal-to-noise ratio and facilitates auditory scene analysis for multiple auditory object tracking, and is briefly summarised in the following paragraphs -for more details, please refer to [START_REF] Faller | Source localization in complex listening situations: Selection of binaural cues based on interaural coherence[END_REF].

The ITD and IC, denoted respectively by τ (n) and c 12 (n), where n indexes the sample currently being processed, are estimated from the normalised cross-correlation functions of the signals from left and right ear for each centre frequency f c , respectively x ′ 1 and x ′ 2 . The normalisation of the cross-correlation function is introduced in order to get an estimate of the IC, defined as the maximum value of the instantaneous normalised cross-correlation function. This estimate describes the coherence of the left and right ear input signals. In principle, it has a range of [0; 1], where 1 occurs for perfectly coherent x ′ 1 and x ′ 2 . However, due to the DC offset of the halfwave rectified signals, the values of c 12 are typically higher than 0 even for independent (nonzero) x ′ 1 and x ′ 2 . Thus, the effective range of the interaural coherence c 12 is compressed to [a; 1] by the neural transduction. The compression is more pronounced (larger a) at high frequencies, where the low pass filtering of the half-wave rectified critical band signals yields signal envelopes with a higher DC offset than in the signal wave forms [START_REF] Faller | Source localization in complex listening situations: Selection of binaural cues based on interaural coherence[END_REF].

The ILD, denoted as ∆L(n), is then computed using the signal levels at the corresponding offsets [START_REF] Faller | Source localization in complex listening situations: Selection of binaural cues based on interaural coherence[END_REF]. Note that due to the envelope compression the resulting ILD estimates will be smaller than the level differences between the ear input signals. For coherent ear input signals with a constant level difference, the estimated ILD (in dB) will be 0.23 times that of the physical signals [START_REF] Faller | Source localization in complex listening situations: Selection of binaural cues based on interaural coherence[END_REF].

When several independent sources are concurrently active in free field, the resulting cue triplets {∆L(n), τ (n), c 12 (n)} can be classified into two groups [START_REF] Faller | Source localization in complex listening situations: Selection of binaural cues based on interaural coherence[END_REF]: (1) Cues arising at time instants when only one of the sources has power in that critical band. These cues are similar to the free-field cues -localisation is represented in {∆L(n), τ (n)}, and c 12 (n) ≈ 1. (2) Cues arising when multiple sources have non-negligible power in a critical band. In such a case, the pair {∆L(n), τ (n)} does not represent the direction of any single source, unless the superposition of the source signals at the ears of the listener incidentally produces similar cues. Furthermore, when the two sources are assumed to be independent, the cues are fluctuating and c 12 (n) < 1. These considerations motivate the following method for selecting ITD and ILD cues. Given the set of all cue pairs, {∆L(n), τ (n)}, only the subset of pairs is considered which occurs simultaneously with an IC larger than a certain threshold, c 12 (n) > c 0 . This subset is denoted

{∆L(n), τ (n)|c 12 (n) > c 0 } (1)
The same cue selection method is applicable for deriving the direction of a source while suppressing the directions of one or more reflections. When the "first wave front" arrives at the ears of a listener, the evoked ITD and ILD cues are similar to the free-field cues of the source, and c 12 (n) ≈ 1. As soon as the first reflection from a different direction arrives, the superposition of the source signal and the reflection results in cues that do not resemble the freefield cues of either the source or the reflection. At the same time IC reduces to c 12 (n) < 1, since the direct sound and the reflection superimpose as two signal pairs with different ITD and ILD. Thus, IC can be used as an indicator for whether ITD and ILD cues are similar to free-field cues of sources or not, while ignoring cues related to reflections.

Faller and Merimaa's cue selection method, as the authors point out, can be seen as a "multiple looks" approach for localisation, which provides the motivation for our implementation. Multiple looks have been previously proposed to explain monaural detection and discrimination performance with increasing signal duration [START_REF] Viemeister | Temporal integration and multiple looks[END_REF]. The idea is that the auditory system has a short-term memory of "looks" at the signal, which can be accessed and processed selectively. In the context of localisation, the looks would consist of momentary ITD, ILD, and IC cues. With an overview of a set of recent cues, ITDs and ILDs corresponding to high IC Fig. 4. Example of the use of an adaptation of the cue selection method proposed by [16] using a 1 s "multiple looks" buffer. Represented in the figure is a histogram of collected ITD cues corresponding to high IC levels for a particular frequency channel of a 1 s audio snippet. This histogram is interpreted as a distribution corresponding to the probability of the occurrence of ITD readings, which is then used as a conspicuity map in order to perform a summary cross-correlogram over all frequencies (see main text for more details).

values are adaptively selected and used to build a histogram that provides a statistical description of gathered cues (see Fig. 4).

Finally, the binaural processor capitalises on the multiple looks configuration and implements a simple auditory scene analysis algorithm for detection and extraction of important auditory features to build conspicuity maps and ultimately a saliency map, thus providing a functionality similar to the role of the external nucleus of the inferior colliculus (ICx) in the mammalian brain. The first stage of this algorithm deals with figure-ground (i.e. foreground-background) segregation and signal-to-noise ratio. In signal processing, the energy of a discrete-time signal x(n) is given by [

19] E = ∞ -∞ |x(n)| 2
Using this notion, a simple strategy can be followed to selectively apply the multiple looks approach to a binaural audio signal buffer so that only relevant audio snippets are analysed. This strategy goes as follows: given a binaural signal buffer of N samples represented by the tuple {x ′ 1 (n), x ′ 2 (n)}, the average of the energies of the component signals

x ′ 1 (n) and x ′ 2 (n) is E avg = N 1 |x ′ 1 (n)| 2 + N 1 |x ′ 2 (n)| 2 2 (2)
and can be used as a noise gate so that only when E avg > E 0 ITDs, ILDs and ICs triplets are collected for the buffer, yielding multiple looks values only for relevant signals (just the ITD-ILD pairs corresponding to high IC values are kept in conspicuity maps per frequency channel), while every other buffer instantiation is labelled as irrelevant noise. E 0 can be fixed to a reasonable empirical value or be adaptive, as seems to happen with human hearing. A set of results exemplifying this algorithm is presented on Fig 5.

Once the multiple looks information is gathered, since ITDs are proven to be stable across frequencies for a specific sound source at a given azimuth regardless of range or elevation, the ITD conspicuity maps may be summed over all frequencies, in a process similar to what is believed to occur in the ICx, in computational terms known as a summary cross-correlogram (again see Fig. 4). From the resulting one-dimensional signal, the largest peaks may be Fig. 5. Binaural processing results of an approximately 30 second-long audio snippet of a typical "cocktail party" scenario, with the main voice calling out "Nicole, look at me", while other voices can be heard coming from sites close to the robotic head, elsewhere in the lab. The active perception head was moved while the main voice speaker was kept still, first keeping the speaker to the right and slowly travelling towards the centre, then keeping the speaker to the left and again slowly moving towards the centre. Top -the effect of the signal power-based figure-ground segregation noise gate is shown (dashed line represents gate threshold); Middle -ITD estimates for the most salient sound; Bottom -corresponding azimuth estimates. These results show the performance of the binaural processor under difficult conditions, the only "failure" being the estimates corresponding to the 14 s instant: for a signal power above the interest threshold, the background noise (i.e., some other voice in the lab) was more salient than the main voice.

taken as having been effected by the most important soundsources represented in the auditory image. Then, a search is made across each frequency band to find the closest ITD and its ILD pair, for each reference ITD, thus building n-sized vectors (for m = n-1 frequency channels) for each relevant sound source of the form

A = [τ, ∆L(f 1 c ) • • • ∆L(f m c )] (3)
C. Bayesian sensor model Finally, regarding the Bayesian 3D sound-source localisation unit, auditory sensor space is defined as a log-spherical volumetric occupancy grid Y ′ , with each cell being indexed by its far corner C ≡ (log b ′ ρ max , θ max , φ max ) ∈ C ′ ⊂ Y ′ -this configuration follows the same formalism as the Bayesian Volumetric Map (BVM) framework, described in [START_REF] Ferreira | Bayesian Models for Multimodal Perception of 3D Structure and Motion[END_REF].

The set of cell indices C ′ in auditory space is a subset of cell indices C in BVM space, and thus the indexing used for the auditory sensor space Y ′ also corresponds to cells in the BVM space Y; however, different resolutions are assumed to accommodate the difference between visual and auditory accuracy ratings. Hence the need of a related base for logspace given by b ′ = a k log a b , ∀a ∈ R. This relation ensures that log-space in the auditory sensor spatial domain and logspace in the BVM are related by factor k, in the sense that one cell in b ′ log-space corresponds to k cells in b log-space. Consequently, the important spatial relation C ′ = N i=1 C is also ensured through the additional relations ∆θ ′ = i × ∆θ and ∆φ ′ = j × ∆φ and allows measurements from one cell in auditory space to precisely correspond to N = i × j × k cells in the BVM (for more details please refer to [START_REF] Ferreira | Bayesian Models for Multimodal Perception of 3D Structure and Motion[END_REF]).

The direct audition sensor model is formulated as the first question of the Bayesian Program in Fig. 6, where all relevant variables and distributions and the decomposition of the corresponding joint distribution, according to Bayes' rule and dependency assumptions, are defined. The use of the auxiliary binary random variable S C , which signals the presence or absence of a sound-source in cell C, and the corresponding family of probability distributions P (S C |O C C) ≡ P (S C |O C ) promotes the assignment of probabilities of occupancy close to 1 for cells for which the binaural cue readings seem to indicate a presence of a sound-source and close to .5 otherwise (i.e. the absence of a detected sound-source in a cell doesn't mean that the cell is empty). The family of distributions is given in tabular fashion: obviously, P ([

S C = 1]|[O C = 0]) = 0 and P ([S C = 0]|[O C = 0]) = 1, while P ([S C = 1]|[O C = 1]) = P SS and P ([S C = 0]|[O C = 0]) = 1 -P SS
are dependent on the probability assigned to an occupied cell corresponding to the position of a sound-source, denoted by P SS , which can be empirically chosen or statistically learned through the analysis of several typical perceptual scenarios.

The second question corresponds to the estimation of the position of cells most probably occupied by sound sources, through the inversion of the direct model through Bayesian inference on the joint distribution decomposition equation. The former is used as a sub-BP for the BVM, while the answer to the latter yields a gaze direction of interest in terms of auditory features for the multimodal attention system, using a maximum a posteriori (MAP) method.

Program                                                                                                                                          Description                                                                                                                        Specification                                                                                                         

III. BAYESIAN BINAURAL SYSTEM CALIBRATION

As can be seen on the BP in Fig. 6, calibration of the binaural system involves the characterisation of the families of normal distributions P (τ

|S C O C θ max ) and P (∆L(f k c )|τ S C O C C) ≈ P (∆L(f k c )|S C O C C
) through descriptive statistical learning of their central tendency and statistical variability. This is done in an equivalent manner as with commonly used head-related transfer function (HRTF) calibration processes (see, for example, [START_REF] Calamia | Three-dimensional localization of a close-range acoustic source using binaural cues[END_REF]) and is described in the following paragraphs.

A set M C of n-dimensional measurement vectors such as defined in equation ( 3) is consequently collected per cell C ∈ C ′ . The full set of collected measurement vectors for all cells in auditory sensor space Y ′ is expressed as M = M C . Denoting M C = M \ M C as the set of measurements for all cells other than C, the statistical characterisation process of each family of distributions is effected through

P (τ |[S C = 1] O C θ max ) ≡ N (τ, µ τ (M C ), σ τ (M C )) (4a) P (τ |[S C = 0] O C θ max ) ≡ N (τ, µ τ (M C ), σ τ (M C )) (4b) P (∆L(f k c )|[S C = 1] O C C) ≡ N (∆L(f k c ), µ ∆L(f k c ) (M C ), σ ∆L(f k c ) (M C )) (4c)
P (∆L(f k c )|[S C = 0] O C C) ≡ N (∆L(f k c ), µ ∆L(f k c ) (M C ), σ ∆L(f k c ) (M C )) (4d)
Auditory calibration is performed by presenting a broadband audio stimulus through a loudspeaker positioned in well-known spatial coordinates corresponding to the geometric centre of each cell C ∈ C ′ so as to sample space according to the auditory sensor space Y ′ . The experimental setup used for this purpose is described in Fig. 7.

The acquisition method may be simplified by a factor of 4 by taking into account the spatial redundancies of auditory sensing, namely the symmetry enforced by the backto-front ambiguity and the left-to-right antisymmetry for both ITDs and ILDs, to reduce calibration space to the front-left quadrant.

A further simplification of the procedure consists in positioning the loudspeaker, for each of the N d considered distances from the binaural system, precisely in front of the active perception head (i.e. (θ, φ) = (0, 0)) and to rotate the active head so that the whole range of azimuths and elevations of the auditory sensor space is covered. This replaces the several minutes taken to reposition the loudspeaker by hand (now only happening N d times) by a few seconds of head motions for each cell. The full procedure is depicted in Fig. 8 for a typical calibration context.

IV. CONCLUSIONS

Full 3D auditory localisation has rarely been explored in robotic applications (see, for example, [START_REF] Calamia | Three-dimensional localization of a close-range acoustic source using binaural cues[END_REF] for a review on this subject); this work contributes with a novel probabilistic solution that produces localisation estimates based on binaural cues yielded by a robust binaural processing unit. This solution has been designed so as to provide a sensor model to be used by a multimodal perception framework, the Bayesian Volumtric Map, described in [START_REF] Ferreira | Bayesian Models for Multimodal Perception of 3D Structure and Motion[END_REF]. Further details on the implementation of these models can be found at http://paloma.isr.uc.pt/~jfilipe/ BayesianMultimodalPerception.

Fig. 1 .

 1 Fig. 1. View of the current version of the Integrated Multimodal PerceptionExperimental Platform (IMPEP), on the left. On the right, the IMPEP perceptual geometry is shown: {E} is the main reference frame for the IMPEP robotic head, representing the egocentric coordinate system;{C l,r } are the stereovision (respectively left and right) camera referentials; {M l,r } are the binaural system (respectively left and right) microphone referentials; and finally {I} is the inertial measuring unit's coordinate system.

Fig. 2 .

 2 Fig. 2. Typical application context of the IMPEP active perception system.

Fig. 3 .

 3 Fig. 3. The IMPEP Bayesian binaural system.

Fig. 6 .

 6 Fig. 6. Bayesian Program for binaural sensor model.

Fig. 7 .

 7 Fig. 7. Experimental setup for the binaural system calibration procedure.

Fig. 8 .

 8 Fig. 8. Schematic of the experimental acquisition method of the auditory calibration procedure. A typical auditory sensor space is depicted, roughly based both on known human auditory precision ratings [20] and the specifications of the pan and tilt unit model PTU-46-17.5 from Directed Perception, characterised by an azimuthal range of [0 o ; 90 o ] with ∆θ = 2 o resolution and an elevation range of [-30 o ; 30 o ] with ∆φ = 10 o resolution, and a number of partitions N = 4, of which only the furthest N d = 2 partitions would be used -the number of cells to sample amount to [(90/2) × (60/10) × 2] = 540. Thus, for 20 readings of a 1 s broadband stimulus per cell and approximately 1 s to set up each position for the loudspeaker, calibration would take about 10800×2 s = 6 h. This procedure allows partitioning calibration into one session per loudspeaker distance, in this case of ≈ 3 h.

Proceedings of the 2008 International Conference on Cognitive Systems University of Karlsruhe, Karlsruhe, Germany, April 2-4, 2008


Xiaolong Zhang

Liquid

Liquid film dynamics in

come

Presented by

Xiaolong ZHANG 06-08 Oct. 2021, Toulouse, France Email: xiaolog.zhang@cea.fr 01 METHODOLOGY

LIQUID FILM DYNAMICS

Pulsating heat pipe (PHP, Caloduc Oscillant in French)

• Promising, high-efficient, passive heat transfer devices, particularly in thermal control system for aerospace applications. • A capillary meandering between evaporators (heat sources) and condensers (heat sinks), with optional adiabatic sections. • Partially filled with working fluid: vapor bubbles and liquid plugs.

• Self-sustained oscillation of liquid plugs as a result of pressure imbalance among the bubbles. • Liquid films are deposited by receding plugs.

Motivation

• Liquid film dynamics is crucial, because heat exchange through the films is its dominant mechanism. • Modeling of PHPs (Ph.D. thesis of X. ZHANG).

PHASE 1: Film Dynamics

• Two-dimensional axisymmetric hydrodynamics of thin liquid film, when the meniscus oscillates from pinned contact line. Key phenomena related to the liquid film:

• Length variation: meniscus oscillation and contact line motion (receding) determined by local superheating and contact angle. • Thickness variation: liquid evaporation and condensation of the film.

PHP MODELING AND SPACE EXPERIMENTS

Single Branch PHP Modeling

CONCLUSION

Attempts of establishing empirical correlations for predicting PHP performance have failed, because of the large number of system parameters and the complexity of PHP behavior. Therefore, the thermo-hydrodynamic simulation that reflects the physical phenomena remains the most promising way for developing PHP design tools in the future. Compared to more complicated numerical models, the one-dimensional description of PHP appears as an optimal choice that can produce accurate results without oversimplification.

• Film dynamics: The future PHP models should account for the liquid film phenomena presented above: spatial and temporal variation of the film thickness during meniscus oscillation and film thinning due to liquid evaporation; contact line recedes in fluid channel as a result of capillary action and liquid evaporation, which may be simplified as a function of local contact angle and superheating and be integrated into PHP system modeling. • In the next stage of CASCO (Code Avance de Simulation de Caloduc Oscillant) development, the constant liquid film thickness, as an input parameter, will be replaced with the liquid film model presented in single branch PHP. To validate the model, extensive comparisons with experimental data will be conducted. • When the meniscus advances over the liquid film, the interface has a wavy appearance (ripples). The crest and the trough are truthfully captured by the numerical calculation.

Contact angle variation

• Variation of θ is periodic, because θ is imposed by the film shape (controlled by the meniscus motion).

• The initial θ is the maximum value achieved during oscillation. θ remains nearly constant for a large part of a period, because the meniscus is away from the contact line.

Film dynamics with pinned contact line and meniscus oscillation

• The lubrication approximation and Kelvin effect.

• local contact angle θ and meniscus curvature R m are oscillating.

• meniscus velocity U is low; shape determined by surface tension (hemispherical). • Films: mostly flat and thin (~100 μm), low Re number and strong viscous forces.

Contact line motion: dewetting problem

• contact line receding because of capillary action and liquid evaporation. • Minimal contact angle decreases with ΔT, which implies that it will reach zero for certain ΔT. For larger ΔT, the solution does not exist.

• ΔT here is merely around 1mK, which is very small, compared to the value we encounter in PHP. A growing dewetting ridge has been observed. • Ca cl monotonically increases with both ΔT and θ.

• For large ΔT, the dependence of Ca cl on ΔT is more significant than that on θ , which signifies that the evaporation effect prevails. x cl

• Curvatures are replaced by the model lumped shape.

• Vapor is superheated, described as the ideal gas.

• Liquid-vapor interface has the saturation temperature corresponding to vapor pressure p. • Film deposited by the receding meniscus is a function of meniscus speed (from film dynamics). • Contact line recedes as a function of local superheating and contact angle (from film dynamics). • Film thickness 𝛿 is uniform, but evolves with time.

• Results will be compared to experimental data. t

  X., & Nikolayev, V. (2021). Liquid film dynamics with immobile contact line during meniscus oscillation. Journal of Fluid Mechanics, 923, A4. doi:10.1017/jfm.2021.540 • Zhang, X., & Nikolayev, V. Liquid film dynamics during meniscus oscillation. Proceedings of IHPC (International Heat Pipe Conference) / IHPS (International Heat Pipe Symposium), Gelendzhik, Russia, September, 7-10, 2021

••

  A capillary tube with an end sealed and the other open connect to pressure reservoir of constant p e . • It consists of an evaporator (E), an adiabatic section (A) and a condenser (C). Inside the tube a bubble, a meniscus (film) and a plug. One dimensional PHP simulation code, developed by V. Nikolayev in CEA. • Complex PHP loops are projected on a straight tube from a starting point, applied with repeating wall conditions: evaporator (E), adiabatic section (A) and condenser (C)… • Currently, liquid films are represented by a constant thickness (input parameter). PHP experiments in space PHP development for the satellites • ESA MAP (Microgravity Application Program) project, TOPDESS (Twophase Passive Thermal Devices for Deployable Space Systems). • ESA Space PHP experiment on HTH (Heat Transfer Host) on board of International Space Station (foreseen for 2026). • CNES R&T DSO/TB/TH-2020,0026119: collaboration of Airbus DS and CEA.

  to experiments of Lips et al. (2010) • An oscillating liquid plug in a cylindrical capillary tube.

  speed as a function of ΔT and θ

  • dewetting ridge forms by collecting the liquid. Matter Physics (SPEC), Alternative Energies and Atomic Energy Commission of France (CEA) Xiaolong ZHANG Thesis director: Vadim NIKOLAYEV Space PHP prototype (left) and schematic (right)

  • Contact line receding must be included in calculation. Experiments of L. Fourgeaud et al. (2016) • An ethanol film is first deposited on solid substrate by a receding liquid meniscus. • Dewetting occurs and the contact line recedes.

  53°C, T c =10°C, L e =15cm, L a =5cm, L c =20cm. Thesis structure Objectives • PHP performance in space: PHP with large inner diameter (to reduce the viscous resistance); can only function under microgravity condition. • Test of PHP prototype for cooling electronic devices in satellites. • Obtain and use experimental data for validation and improvement of numerical design tools. Numerical calculations demonstrate • The growth of the dewetting ridge • Contact line receding speed logarithmically decreases with the ridge size: nearly constant under most experimental conditions Results are verified by comparison to analytical formula proposed by J. Snoeijer & J. Eggers(2010)

•

  Effect of contact line receding on the liquid film.

	PHASE 2: Single Branch PHP

Acknowledgements

• EU-CEA Numerics program.

• CNES co-finding for the Ph.D. thesis.

• ESA (European Space Agency) TOPDESS project.


domain: phys spi spi.meca.mefl

text:
Xiaolong Zhang
email: xiaolong.zhang@cea.fr

Vadim S Nikolayev

Liquid film dynamics during meniscus oscillation

come

Introduction

Pulsating Heat Pipe (PHP), a heat transfer device that is made out of a simple capillary tube, usually contains a sequence of liquid plugs and vapor bubbles of the same substance. These bubbles have a length much larger than the tube radius. In a properly functioning PHP, the liquid plugs oscillate between the evaporator (heated section) and the condenser (cooled section). Liquid films that are deposited when the menisci recede, provide the principal channel for the heat and mass exchange in PHP [1].

Previous studies have established a relationship between the film thickness and the velocity of constant receding meniscus, under the adiabatic conditions [2][3][4]. Additionally, the film thinning due to evaporation has been investigated in the case of the continuous liquid film on tube wall [5][6][7].

However, in practice, because of the evaporation, the liquid film rupture occurs in the evaporator of PHP, creating the triple contact lines where the heat and mass transfers are intense. It leads to a singularity in a narrow vicinity of the contact line, where the film thickness approaches zero and the evaporation rate in extremely high. Janeček and Nikolayev [8] present a theoretical model that relaxes this singularity by accounting for the Kelvin effect. Nonetheless, the impact of meniscus oscillation was not included.

In this article, we consider the simplest oscillation case, where the contact line is immobile (pinned on a surface defect). In order to investigate the film evolution during meniscus oscillation and with wall superheating, a 2D model that describes the dynamics of liquid-vapor interface has been developed. This work is twofold. First, we focus on the vicinity of the immobile contact line under evaporation by ignoring the meniscus motion. Next, we analyze the contact angle response to the meniscus oscillation under adiabatic conditions. Finally, we simulate both effects together.

Description of the two-dimensional physical model

Consider a circular capillary tube of an inner radius R, containing a liquid plug and a long vapor bubble of the same fluid. Due to the film evaporation, a contact line appears and is pinned. The physical processes and the vapor-liquid interface are assumed to be symmetric with respect to the tube axis, and film thickness h is small as compared with the radius, h R. Therefore, the problem can be described in 2D, in the cross-section plane of the tube, shown in Fig. 1.

liquid vapor film 2R U contact line x y h(x,t)
tube wall:

𝑇𝑇 w = 𝑇𝑇 sat + Δ𝑇𝑇 x m h s x s R m parabola (x f ,h f ) Figure 1.
Meniscus moving through a capillary tube filled with liquid; the contact line is pinned at x = 0.

The tube is assumed to be thin enough to make the gravity force negligible. Because the vapor density and viscosity are small compared to the liquid, one can safely assume that the vapor pressure p is spatially homogeneous in the bubble and neglect the vapor-side viscous stress on the liquid-vapor interface.

Adopting the classical Bretherton approach [3], the interface is divided into the film and meniscus regions. In the film region, the interface is nearly flat (|∂h/∂x| 1), it can thus be described with the lubrication approximation of the hydrodynamics equation. The meniscus region is the part of interface situating in the tube center, where the interface is controlled by the surface tension alone, thus being of constant curvature (shown in Fig. 1 with a circle of radius R m and center at x s ).

With aforementioned assumptions, the lubrication theory [9] for the film region results in the equation that describes the interface dynamics,

∂h ∂t + ∂ ∂x h 3 3µ ∂∆p ∂x = - J ρ , (1)
where h = h(x, t) is the local film thickness, and J is the mass flux across the interface defined to be positive at evaporation. Here, µ and ρ are the liquid shear viscosity and density, respectively. ∆p is the interface pressure jump. Conventional hypotheses concerning the liquid film mass exchange [9] are applied. A linear temperature profile in the radial direction is assumed in the thin liquid film, thus the energy balance at the interface results in the mass flux

J = k(T w -T int ) hL , (2)
where T int is the temperature of the vapor-liquid interface, k is the liquid heat conductivity, and L is the latent heat. The problem is formulated here for a general case where the tube wall can be superheated or subcooled with the respect to the saturation temperature T sat corresponding to the imposed vapor pressure p. The wall superheating is denoted ∆T . The tube wall temperature thus is T w = T sat + ∆T . The pressure jump ∆p = p -p l (with p l , the liquid pressure) across the interface obeys the Laplace equation

∆p = σ K + 1 R -h σ ∂ 2 h ∂x 2 + 1 R m , ( 3
)
where σ is the surface tension, K is the 2D interface curvature in the axial cross-section shown in Fig. 1; in the small-slope approximation, K ∂ 2 h/∂x 2 . Because of this limitation, such a film theory is not able to describe the meniscus region. The radial contribution to the curvature R -1 m is assumed to be independent of x, because h is much smaller than R. At large x, Eq. ( 1) results in an increasingly larger h, where the surface tension gradually predominates over the viscous forces and controls the interface, so its curvature ∂ 2 h/∂x 2 becomes constant. Such a condition corresponds to a parabolic shape in the axial plane. This parabola matches the circular meniscus, which results in the condition

∂ 2 h/∂x 2 | x=x f = R -1 m , (4)
defined at the ending point (x f , h f ) of the film region.

Relaxing the contact line singularity with the Kelvin effect

If the vapor-liquid interface was at a fixed saturation temperature (T int = T sat ), the mass flux J (Eq. 2) would diverge at the contact line h = 0 as J ∼ ∆T /h, which is nonphysical because total evaporated mass (the integral of J) would be infinite. The Kelvin effect, i.e. the dependence of T int on the interfacial pressure jump ∆p

T int = T sat [1 + ∆p/(ρL)] (5)
can relax the singularity [8], because it allows T int to vary along the interface. Equations (2, 5) result in

J = [∆T -∆pT sat /(ρL)]k/(hL). (6)
T int now can attain the wall temperature T w at the contact line, which results in J(x → 0) = 0 and ∆p(x → 0) = ∆p cl , where a constant pressure jump at the contact line is introduced: ∆p cl = Lρ∆T /T sat . By combining Eqs. (1,6), one obtains the governing equation.

Unlike other microscopic approaches [10], the problem is now regular (i.e. ∆p is not divergent). As the Kelvin effect alone is capable of relaxing the contact line singularity, the other microscopic effects such as hydrodynamic slip, Marangoni effect, and interfacial kinetic resistance [8] are not crucial. They are not included in our model for the sake of clarity.

Contact line model with fixed contact angle and substrate superheating

Consider now the steady case, where the meniscus motion is ignored and the contact angle is fixed. When T w = T sat , the liquid forms a straight wedge near the contact line, as shown in Fig. 2a. When the substrate is heated, a flow inside the wedge brings the liquid towards the contact line to compensate the mass loss due to evaporation, thus creates a viscous pressure drop. Interfacial curvature forms near the contact line [1], which creates a difference between the microscopic contact angle θ micro , the actual slope at the contact line, and the interface slope θ farther away from the contact line, cf. Fig. 2b. The characteristic length for this effect is K ∼ 10 -100 nm [11] from the contact line. At the length scale x meso K but much smaller than the film thickness ∼ 10 µm, one can define the experimentally measurable interface slope θ called the apparent contact angle. Their relationship, cf. Fig. 6a, can be expressed as

θ = θ(∆T, θ micro ), (7)
The region x < 0 < x meso is often referred to as microregion. Numerical calculation of θ has been completed by Janeček and Nikolayev [8]. It is based on the steady version of Eq. ( 12) (i.e. with ∂h/∂t = 0) solved with the boundary conditions (8a). Instead of using ∆p cl = Lρ∆T /T sat , the weaker form (the second equation of (8a)) was employed for numerical proposes.

Other two boundary conditions for this fourth order differential equation are the imposed slope θ micro at the contact line and the condition of zero (microregion scale) curvature (8b):

h (x = 0) = 0, ∂∆p/∂x| x→0 = 0, (8a) ∂h/∂x| x=0 = θ micro , ∂ 2 h/∂x 2 | x=xmeso = 0. ( 8b
)
Note that the interface slope saturates at x K , θ is thus independent of x meso . Figure 3 demonstrates an example of θ as a function of θ micro and ∆T for pentane at 1 bar, cf, Eq. (7). θ(∆T, θ micro ) monotonously grows with both θ micro and the evaporation intensity controlled by ∆T ; evidently, θ(∆T = 0) = θ micro .

One now introduces θ evp-min (∆T ) ≡ θ(∆T, θ micro = 0), (9) which is a lower bound for θ(∆T, θ micro ), i.e. its value for the complete wetting case [11]. Therefore, θ micro = 0 curve in Fig. 3 represents θ evp-min (∆T ). For example, for ∆T =25 mK, one obtains θ evp-min (∆T = 25 mK) 7.8 • , which signifies that the apparent angle θ will always be larger than 7.8 • at this ∆T value, whatever are the wetting conditions.

Interface profile and contact angle during meniscus oscillation

In this section, an adiabatic problem is considered, ∆T = 0 and ∆p cl = 0, which means θ = θ micro . The mass exchange appearing at the macroscopic scale is insignificant and can be safely neglected. A meniscus oscillates from the pinned contact line. The center of meniscus x m , which is the experimentally measurable quantity, travels periodically with a period P and an amplitude A. One can assume its harmonic oscillation

x m (t) = x i + A[1 -cos(2πt/P )], (10)
where x i is the initial meniscus center position. Initially, the equilibrium profile satisfying this condition ∂h/∂t = 0 is used in Eq. ( 1). The contact line is pinned at x = 0. The contact angle θ micro is expected to vary between the static advancing θ adv and the static receding θ rec angles (the difference of which defines the wetting hysteresis) while the contact line remains immobile. The governing equation is solved for x ∈ [0, x f ]. The boundary conditions are Eqs. (4, 8a) and the condition h(x = x f ) = h f where h f is yet to be defined. As explained in sec. 2, the condition (4) means that the solution matches a parabola at x = x f . The parabola is in fact a small-angle approximation of the circular meniscus shape, so its geometrical parameters can be defined from those of the meniscus (cf. Fig. 1), in particular from its curvature radius and the position of circle's lowest point (x s , h s ). The conditions of matching are both the continuity and the smoothness:

h f = (x f -x s ) 2 /(2R m ) + h s , ∂h/∂x| x=x f = (x f -x s )/R m , (11)
where h s is defined with the condition R m + h s = R [4]; x s satisfies the equation R m + x s = x m (cf. Fig. 1).

In the present algorithm, x f is chosen in such a way that the difference x m -x f remains constant during oscillations. This way, h f = h(x f ) remains large with respect to the film thickness and more or less constant. As discussed in [4], the solution is nearly independent of the specific choice of h f (and thus of x f ). The film shape and the unknown parameters (R m , x s , h s , h f ) can be determined at each time moment. If the oscillation is sufficiently slow, the meniscus radius R m is close to R throughout the oscillation but deviates from it for larger meniscus velocities.

Numerical implementation

For the harmonic oscillation case, the meniscus velocity is U (t) = U 0 sin(2πt/P ). The velocity amplitude U 0 = 2πA/P is convenient to be chosen as a characteristic velocity to define the capillary number Ca 0 = µU 0 /σ. Other variables are converted into dimensionless form as (dimensionless quantities are denoted with a tilde): x = αRCa

1/3 0 x, h = αRCa 2/3 0 h, t = αRCa 1/3 0 U -1 0 t, p = σα -1 R -1 p
and T = T sat T , where α 1.34. αRCa 2/3 0 is the Bretherton film thickness [3] corresponding to the maximal receding velocity U 0 .

With such a "natural" scaling four main dimensionless parameters are left: θ i , P , Ca 0 , and N . All the quantities will be studied in this parametric space. The dimensionless amplitude is linked to the period as à = P /(2π). The governing equation in dimensionless form reads

∂ h ∂ t + ∂ ∂ x h3 3 ∂∆p ∂ x = N ∆p -pcl h ( 12
)
where,

N = µkT sat (LραR Ca 0 ) -2 , (13)
which describes the magnitude of the Kelvin effect in the microregion. Equation ( 12) is discretized using the finite volume method (FVM), which is more stable numerically [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF] than the more conventional finite difference method. The variables such as h and their even order derivatives are defined at center of each grid segment (called node), while the odd order derivatives are defined at the segment ends. The FVM has the advantage that the liquid flux is continuous at the segment ends. Nonlinear terms are managed by iteration: they include values from the previous iteration. The numerical algorithm is similar to that used by Nikolayev [9]. The mesh size is exponentially refined near the contact line (as x → 0) to capture the contact angle variation without considerably increasing the total number of nodes [9]. 4.2. Interface profile during oscillation Fig. 4 illustrates numerical results of the interface profiles at several time moments during oscillation. The meniscus motion follows the harmonic law (10). The liquid film is deposited until t = P/2. For t > P/2, the meniscus advances over the deposited film. The interface has the wavy appearance near the end of film, resembling the ripples that appears during the steady meniscus advancing discussed in Nikolayev and Sundararaj [5]. The interface profiles h(x, 0) and h(x, P ) are indistinguishable, which confirms the periodicity of oscillation. The numerical results are compared with the experiments of Lips et al. [START_REF] Lips | [END_REF]. In the experiments, a capillary tube (inner radius 1.2mm) contains two vapor bubbles and a short liquid plug of pentane. The variation of pressure difference between the vapor bubbles forces the plug oscillate, with the frequency 3.7 Hz and the velocity amplitude U 0 = 0.24 m/s approximately. The plug motion is recorded by the high-resolution camera and the velocity and the curvature of meniscus are then recovered from the image analysis. Because the motion is not perfectly harmonic, these experimental data were used in the numerical calculations in Eqs. (10) and (4).

In their experiments, the Weber number W e 0 = 2RρU 2 0 /σ = 5.66 is larger than unity. The Reynolds number Re 0 = 2RρU 0 /µ = 1521 is quite high and the impact of inertia on the shape of the central meniscus part must be taken into consideration. However, the thin film can still be considered as controlled by the viscosity only, because Ca 0 ∼ 3.74 × 10 -3 remains low. Under these conditions, the film ripples in the transition region close to the meniscus, which appear during the plug advancing, can be compared to the calculations.

Figure 5 presents several snapshots of plug oscillation. The left column shows the original images of Lips et al. [START_REF] Lips | [END_REF] and the corresponding simulation results. The liquid film in transition region (in the block) is enlarged in the right column. The wavy appearance of the interface is truthfully captured by the numerical calculation.

Unfortunately, the quantitative comparison of the film thickness (i.e. the vertical coordinate) is hardly possible since the cylindrical aberration of the glass capillary is not corrected and the contact line cannot be distinguished. One can nonetheless compare the axial lengths. The size of one pixel in mm can be obtained from the known outer tube diameter (4 mm) that is visible in the original images. From the experimental images, the axial distance between the crest and the trough (Fig. 5) of the film ripple is 0.50 ± 0.02 mm, while in the simulation, it is 0.51 ± 0.01 mm. This distance is almost constant in time. Evidently, the agreement is excellent.

Contact angle variation

A typical variation of θ during oscillation is plotted in Fig. 6a. The initial contact angle θ i is the maximum contact angle achieved during the periodical motion.

In the beginning of a period, the capillary forces lead to the fast contact angle reduction until the meniscus recedes far enough so the curvature gradient reduces and the contact angle becomes nearly constant for a large part of a period. This nearly constant value is quite insensitive to both θ i and P . θ slowly decreases until the ripples in the near-meniscus region approach the contact line during the backward stroke (Fig. 6b) at the end of a period. This causes the contact angle oscillations, during which its minimal value θ min is attained:

θ(t) ≥ θ min , (14)
which depends quite weakly on θ i (Fig. 6a). This minimal value of contact angle is of importance (cf. sec. 5 below). Since the motion is periodical, the contact angle θ i is attained at t = P .

Combining film evaporation and meniscus oscillation

In this section we consider the meniscus oscillations together with the evaporation effect (i.e., for a positive ∆T ). Two approaches are considered. First, by applying a multi-scale reasoning, we combine the results of sections 3 and 4. Second, the problem of meniscus oscillation with evaporation is solved numerically.

Multi-scale approach

In this section we combine the above results for the case of heating conditions. Consider a case of small ∆T . The evaporation in the film and meniscus regions can thus be neglected during an oscillation period and the macroscopic results of sec. 4 still apply. However, because of the contact line singularity, the evaporation is not negligible in the microregion and leads to a difference between the apparent contact angle θ and the microscopic contact angle θ micro , as defined in sec. 3. First, instead of solving Eq. ( 12) for this case, we apply here the multi-scale reasoning introduced in sec. 3. When a solution of the full problem exists, the micro and macro regions can be matched for a given ∆T . They are connected through the relationship (7), which means the oscillatory film shape eventually imposes the value of θ micro .

In the sprite of the multi-scale reasoning, θ must simultaneously satisfy the conditions:

1. Similarly to the case where ∆T = 0, during oscillations with a pinned contact line, the moving film imposes the θ value. This implies that the inequality (14) should hold for ∆T > 0, cf. Fig. 6a.

2. On the other hand, θ is defined by Eq. ( 7), which is bounded from below: θ ≥ θ evp-min (∆T ).

The condition (ii) suggests that θ(t), throughout oscillation, should remain larger than θ evp-min (∆T ). From the inequality (14), one obtains

θ min ≥ θ evp-min (∆T ), (15)
which presents a necessary condition for matching of two regions. With the equality sign θ min = θ evp-min (∆T max ), this equation defines a superheating limit ∆T max . Since θ evp-min (∆T ) is an increasing function (cf. Fig. 3), the superheating limit ∆T max is an upper bound. Thus the inequality (15) is satisfied when ∆T < ∆T max .

To determine ∆T max graphically, Fig. 7 presents an example where θ(t) (bottom and vertical axes, extracted from Fig. 6a) is plotted together with θ evp-min (∆T ) (top and vertical axes, extracted from Fig. 3). During the oscillation, the minimum value θ min 2 • is attained. From the dependence θ evp-min (∆T ), one can deduce that ∆T max 1 mK. A solution of oscillation problem with evaporation should be nonexistent if ∆T > 1mK.

Numerical approach

In this section, we describe the results of the numerical solution of Eq. ( 12). In Fig. 8, variations of θ micro during meniscus oscillation are plotted for small ∆T of 0.2mK, 0.5mK, 0.7mK and 0.9mK. θ micro attains the minimal value near the time moment t = 0.968P . The minimum value of θ micro decreases with superheating, which implies that θ micro reaches zero for certain ∆T . The minimal values for ∆T = 0.2mK, 0.5mK, 0.7mK and 0.9mK, are 1.2 • , 0.82 • , 0.53 • and 0.2 • respectively. The simulations confirm that the reasoning of the previous section is correct, because the minimal value of θ micro tends to zero when ∆T approaches 1mK. Moreover, the variation of θ micro is strong near the minimum. It corresponds to the short time lapse when the largest ripple formed near the meniscus (cf. Fig. 4) approaches the contact line, cf. Fig. 9a. The interface curvature exhibits the rapid and strong changes near this time moment, cf. Fig. 9b. After the ripple disappearance (for t > 0.968P ), the pressure jump quickly relaxes to the steady-state profile, in which ∆p saturates towards both the contact line and the meniscus, and decreases monotonically in between, for example see the curve for t = 0.974P .

Paradox and its solution

Both above approaches show that the solution of the oscillation problem is nonexistent when the capillary superheating exceeds a tiny value, much smaller than the realistic superheating. The reason for this paradox is the pinned contact line: if the contact line receded, another degree of liberty would appear so the contact angle would not be constrained any more. The contact line is necessarily depinned when θ attains θ (∆T, 0) (actually, a larger value θ (∆T, θ rec ) but θ (∆T, 0) gives a lower bound).

As this maximum superheating is considerably smaller than that encountered in practice (where it is rather of several degrees K, see e.g. Fourgeaud et al. [7]), one can deduce that the contact line receding at evaporation is inevitable and necessary to be accounted for in PHP modeling.

While the calculation has been completed here only for the pentane case, one can safely state that ∆T max is substantially smaller than realistic superheating used in thermal engineering applications for many other fluids.

Conclusions

Contact line effects have been investigated during meniscus oscillation and film evaporation for the simplest case, where the contact line is pinned. Apparent contact angle varies during oscillation and encounters a minimal value of several degrees. This minimal value exhibits a weak dependence on the initial value of contact angle. Because the apparent contact angle monotonically grows with the superheating for a fixed microscopic contact angle, the minimum places an upper bound for the tube superheating. This upper bound of the order of 1mK is considerably smaller than a typical experimental superheating, which means that the contact line receding is inevitable at oscillations so the dry spots should grow during the meniscus receding in the PHP evaporator. This also means the necessity to include the contact line receding when modeling the oscillation and film evaporation effects simultaneously.

Figure 2 .

 2 Figure 2. Sketch of the straight liquid wedge with no phase exchange (a), and curved wedge at evaporation (b).

Figure 3 .

 3 Figure 3. θ as a function of ∆T for different θ micro computed for pentane at 1 bar. The curve for θ micro = 0 corresponds to θ evp-min (∆T ).

Figure 4 .

 4 Figure 4. Periodic interface shape variation at oscillation for P = 50, Ca 0 = 1 × 10 -3 and θ i = 20 • . The labels give the times corresponding to each profile, and the arrows indicate the meniscus motion direction.

Figure 5 .

 5 Figure 5. Liquid film shape during meniscus oscillation: comparison with experimental results by Lips et al. [13].

Figure 6 .

 6 Figure 6. Contact angle and interface variation for Ca 0 = 10 -3 , P = 50.

Figure 7 .

 7 Figure 7. Function θ(∆T, 0) (the green line, for pentane at 1bar) plotted together with θ variation during meniscus oscillation (blue dashed line, θ i = 20 • ).

Figure 9 .

 9 Figure 9. Film thickness and pressure profiles in the micro-region around the moment t = 0.968P when the minimal θ micro is achieved: ∆T =0.7mK, θ i = 10 • , P = 50.

  Figure 8. Variations of θ micro from t = 0.9P to 0.98P , for ∆T of 0.2mK, 0.5mK, 0.7mK and 0.9mK. θ i = 10 • , P = 50.

		4
		3
	θ micro (°)	2
				ΔΤ =0.2mΚ
		1		ΔΤ =0.5mΚ
				ΔΤ =0.7mΚ
				ΔΤ =0.9mΚ
		0 0.9	0.92	0.94	0.96	0.98
				t/P

Acknowledgments

The present work is supported by the project TOPDESS, financed through the Microgravity Application Program by the European Space Agency. This article is also a part of the PhD thesis of X. Z. co-financed by the CNES and the CEA NUMERICS program, which has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 800945. An additional financial support of CNES awarded through GdR MFA is acknowledged. We are indebted to S. Lips for providing us access to his experimental results.