Keywords: modèles 3D, modèles polyédriques, alignement, normalisation, symmétrie de réflexion plane, analyse en composantes principales

Une nouvelle méthode d'alignement de modèles 3D est proposée. Cette approche est basée sur des propriétés de symétrie et utilise le fait que les techniques d'analyse en composantes principales (ACP) ont de bonnes propriétés vis à vis des plans de symétrie d'un objet tridimensionnel. La recherche rapide d'axes d'alignement optimaux parmi les vecteurs propres de l'ACP est une première étape du processus d'alignement. Les plans de symétrie du modèle 3D sont utilisés comme critère de sélection. Quand l'ACP ne fournit pas tous les axes, nous introduisons un coût d'invariance local de translation (CILT ) qui calcule une mesure de symétrie locale de translation d'une forme le long d'un axe. Les résultats expérimentaux sur la Base de données de Princeton montrent que la méthode proposée calcule des alignements cohérents à l'intérieur de chaque classe de modèles 3D.

Une nouvelle méthode d'alignement de modèles tridimensionnels 1 Introduction

Normalization of 3D models is a common pre-processing stage in many applications in computer graphics such as visualization, 3D Object recognition, 3D shape matching and retrieval [START_REF] Shilane | The Princeton shape benchmark[END_REF][START_REF] Zaharia | 3D versus 2D/3D shape descriptors: A comparative study[END_REF][START_REF] Bustos | An experimental comparison of feature-based 3D retrieval methods[END_REF]. 3D models are generally given in arbitrary scale, position and orientation in 3D-space and most of the methods do not satisfy geometrical invariance, then it is very important to normalize the models into a canonical coordinate frame before any processing. The normalization consists in two steps: the alignment to determine the pose-invariant and the scaling to make the scale-invariant. The alignment has always been the most difficult point in the normalization process.

Toperform an alignment, a concatenation of isometries in 3D-space (translation, rotation and reflection) must be selected to determine the canonical coordinate system. In most of the methods, the center of gravity of the model is chosen as the origin to secure the translation invariance. However, the choice of a good rotation is still a well discussed topic [PRM + 00, VSR01, CO02, SMKF04, ZP04, BKSV04, Kaz07]. The alignment problem addressed in this paper is different from the alignment approaches of [START_REF] Chen | A 3D model alignment and retrieval system[END_REF][START_REF] Kazhdan | An approximate and efficient method for optimal rotation alignment of 3D models[END_REF], where the purpose is to find the best alignment between two given 3D models. Here, we want to compute an intrinsic global coordinate system for each 3D object.

Figure 1: A model of chair. Left: before normalization. Right: well aligned model

When looking at a 3D model, we can say whether it is well aligned or not and we know, in most of the cases, how to find its good alignment (cf. Figure 1). When the 3D model has symmetries, the object is aligned with particular axes or symmetry planes. This is confirmed by Ferguson [Fer00] who noticed that symmetry detection is a key part of human perception and this fact has guided Podolak et al. [PSG + 06] when introducing principal symmetry axes. Our goal is to find a method that best aligns any 3D model (an alignment similar to what a human would select) and will consequently align two similar 3D models in the same way. In this paper, we show that by detecting the planar reflection symmetries we can select a set of good alignment axes. However this method is guaranteed to give the correct alignment for only some cases. Thus keeping only this type of symmetry is insufficient for computing the best alignment for any 3D model. An alternative method is to detect also the the local translational symmetry that has an interesting semantic meaning: the object has the same RR n ¡ 6408 geometrical properties in different parts along a given direction. To build our general alignment algorithm, we proceed as follows. We first focus on discrete detection of plane reflection symmetries and classify a model in terms of its symmetry group and the number of its mirror planes. This classification is used to select the good alignment axes among those found by the principal components analysis (PCA). Then we introduce (LT IC) cost that measures the invariance of a model with respect to local translation about a given direction. This measure is used to compute the remaining alignment axes when the model has at most one good alignment axis given by the PCA. This paper is structured as follows. Section 2 reviews related work on alignment and symmetry detection for 3D models. Section 3 presents our selection of the best alignment axes within the PCA-eigenvectors by analyzing the plane reflection symmetry, and section 4 describes our alignment method. Experimental results of are presented in section 5. Finally, we conclude and point out future work in section 6.

Related Work

The most well-known approach computing the alignment of 3D objects is the principle component analysis method (PCA) [PRM + 00, VSR01, SMKF04, ZP04, BKSV04] based on the computation of moments of 3D models. After a translation of the center of mass to the origin of the coordinate system, three principal axes computed with PCA are used to determine the orientation. The experiences show that PCA-alignment has two disadvantages: (i) It is often imprecise and can produce poor alignments; (ii) The principal axes are not always good at aligning orientations of different models within the same semantic class. Podolak et al. [PSG + 06] introduce a planar reflective symmetry transform (PRST) that computes a measure of the reflectional symmetry of a 3D shape with respect to all possible planes. They use it to define two new concepts for the global coordinate system, the center of symmetry and the principal symmetry axes: the principal symmetry axes are the normals of the orthogonal set of planes with maximal symmetry, and the center of symmetry is the intersection of those three planes. This approach has been improved by Rustamov with the augmented symmetry transform [START_REF] Rustamov | Augmented symmetry transforms[END_REF]. Other methods finding symmetries in 3D models have been presented. Minovic et al. [START_REF] Minovic | Symmetry identification of a 3-D object represented by octree[END_REF] compute symmetries of a 3D object represented by an octree. Their method is based on the computation of a principal octree aligned with the principal axes. Then they compute a measure of symmetry, the symmetry degree, reasoning with the number of distinct eigenvalues associated to the principal axes. Sun and Sherrah [START_REF] Sun | 3D symmetry detection using the extended gaussian image[END_REF] convert the symmetry detection problem to the correlation of the Gaussian image. Then rotational and reflectional symmetry directions are determined using the statistics of the orientation histogram. While Martinet et al. [START_REF] Martinet | Accurate detection of symmetries in 3D shapes[END_REF] use generalized moments to detect perfect symmetries in 3D shapes, Mitra et al. [START_REF] Mitra | Partial and approximate symmetry detection for 3D geometry[END_REF] compute partial and approximate symmetries in 3D objects.

Our goal is to align 3D models using their planar symmetry properties. Our method must be such that similar objects (i.e., objects belonging to a same semantic class) have similar alignments. As noticed in [START_REF] Minovic | Symmetry identification of a 3-D object represented by octree[END_REF], any plane of symmetry of a body is perpendicular to a principal axis. As a result, for models INRIA that have plane reflection symmetries, some PCA-coordinate planes coincide with some mirror planes. Therefore, we have chosen to use the PCA, not for global alignment, but for selection of robust partial alignment features of a model (i.e., only the principal axes that we consider good for a perfect alignment). Given a 3D model, the first key idea is to test the reflection symmetry of the PCA-coordinate planes. According to the result of this test, we select a set of principal axes and use them in our alignment method. When the model has at least two orthogonal mirror symmetries, the PCA gives the good alignment. In the other cases we use the local translational invariance cost along a direction to compute the good alignment axes. Before describing our alignment procedure, let us classify the 3D polygonal models with respect to their plane reflection symmetry and select classes of objects where PCA gives a good alignment.

Symmetry & 3D Objects

In the following M will denote a 3D polygonal model represented by its surface S composed of a set of triangular facets T = {T 1 , ..., T NT }, T i ⊂ R 3 , given by a set of vertices P = {p 1 , ..., p NP },

p i = (x i , y i , z i ) ∈ R 3 .
We study the reflection planes in the symmetry groups [START_REF] Dubrovin | Modern geometry, methods and applications[END_REF], and use them to discriminate different classes of mirror symmetry. Then, we discuss for each class when the PCA alignment has good properties with respect to the planar reflective symmetry.

Plane Reflection Symmetry Analysis

A plane reflection symmetry is defined by a mirror plane π that can be parameterized by its unit normal n and its scalar distance δ from the origin. This symmetry associates to each point p of S a mirror reflection point q on S defined by: q = p -2 (n T • pδ) n.

According to [START_REF] Dubrovin | Modern geometry, methods and applications[END_REF], studying the plane reflection symmetries of a 3D polyhedral model and the types of symmetry groups, we can distinguish five classes of 3D polyhedral models:

1. G C : 3D models that have cyclic symmetry. They have n mirror planes (n > 1) that pass through a fixed axis, such as a regular n-pyramid, a simple rectangular table (n = 2) and a simple square table (n = 4).

2. G D : 3D models that have dihedral symmetry. They have n mirror planes (n > 1) that pass through a particular axis with one mirror plane perpendicular to the axis, such as a regular n-prism or regular n-bipyramid.

3. G R : 3D models that have rotation symmetry such as the five convex regular polyhedra called platonic solids. It contains three sub-groups: G T of tetrahedral symmetry (6 mirror planes), G O of octahedral symmetry (9 mirror planes) and G I of icosahedron symmetry (15 mirror planes).

4. G 1 : 3D models that have only one plane reflection symmetry. This is the case for many natural objects such as the airplanes, the animals, the humans, the chairs, the cars, etc.

5. G 0 : 3D models that don't have any plane reflection symmetry, such as the plants and the trees.

RR n

¡ 6408

This classification is valid for perfect plane reflection symmetries. We will extend it to approximate mirror reflections (cf. section 4.1).

Figure 2: Models from different classes of mirror symmetry:

Tetrahedron ∈ G T , Octahedron ∈ G O , Icosahedron ∈ G I , Table ∈ G even C , Streetlight lamp ∈ G odd C , Hourglass ∈ G even D , 5-prism ∈ G odd D , Chair ∈ G 1 , Plant ∈ G 0 .

Principal Components & Plane Reflection Symmetry Analysis

In this section, we explore the relation between the principal components analysis (PCA) and the plane reflection symmetry analysis. In our proofs, we have retained the "Continuous" PCA (CPCA) [START_REF] Vranic | Tools for 3D-object retrieval: Karhunen-Loeve transform and spherical harmonics[END_REF] because it appears to be more complete and the most stable of all the approaches we have studied. CPCA computes three orthogonal eigenvectors of the covariance matrix C. As noticed in [START_REF] Minovic | Symmetry identification of a 3-D object represented by octree[END_REF], when π is a mirror plane of S and n is the unit normal of π, then π passes through the center of gravity of S and n is an eigenvector of the covariance matrix C that is a principal component axis of S (cf. in appendix, the proof in the continuous case). If S has n mirror planes that pass through a fixed axis (in the cases G C , G D , G R of section 3.1), then n different eigenvectors are associated to the same eigenvalue: in this case, S has a discrete rotational symmetry of order n (n > 1) with respect to the same axis. Besides, we note that if S has a set of dual orthogonal reflection planes, the CPCA detects at least two orthogonal normals associated to one dual orthogonal INRIA mirror plane of this set. In what follows, for each class described in section 3.1 we discuss the position of these vectors with the reflection symmetries:

-If M ∈ G C and n is even (M ∈ G even C
), then the CPCA detects two orthogonal normals associated to two orthogonal reflection planes and the axis of the axial symmetry (the intersection of the mirror planes). When M ∈ G C and n is odd (M ∈ G odd C ), the CPCA gives only one normal associated to one mirror plane.

-If M ∈ G D , then the CPCA gives at least two orthogonal normals; the first is associated to one of the n mirror planes and the second is the axis of the axial symmetry. If, furthermore, n is even, then the CPCA detects the third axis associated to the mirror plane that is orthogonal to the first given mirror.

-If M ∈ G O , then the CPCA detects three orthogonal normals associated to three orthogonal reflection planes, contrarily to the cases of G T and G I , where the CPCA gives only one normal associated to one mirror plane.

-If M ∈ G 1 , then the CPCA gives only one normal associated to its mirror plane.

Thus, when M ∈ G even C G D G O , the CPCA detects at least two good alignment axes and when M ∈ G odd C G T G I G 1 , the CPCA gives only one good alignment axis. Finally, when M ∈ G 0 , the CPCA doesn't detect any good alignment axis. We summarize our discussion using N GA (M), which accounts the number of the good alignment axis computed by the CPCA. Given the symmetry class of the model M, N GA (M) is defined as follows:

N GA (M) =    2 + if M ∈ G even C G D G O 1 if M ∈ G odd C G T G I G 1 0 if M ∈ G 0 4 Alignment of 3D Objects
Given a 3D model M, we aim to develop a general algorithm that will compute N GA (M) and will select the set of the good alignment axes given by the CPCA and if necessary, to compute the rest of alignment axes in order to complete the pose coordinate system. Algorithm: Compute good alignment axes 1. Translate the input 3D model M from its center of gravity to the origin of world coordinate system.

2. Compute the three CPCA eigenvectors v 1 , v 2 , v 3 of the covariance matrix C and rotate the translated model in the new CPCA-coordinate system R (v 1 ; v 2 ; v 3 ) which has the eigenvectors as rows.

3. Test the reflection symmetry for each coordinate plane normal to a CPCAaxis, (xy-, yz-, zx-coordinate plane) and deduce N GA (M). This step is detailed in section 4.1.1.

4. Select the good alignment axis/axes according to the value of N GA (M):

(if N GA (M) = 2 + ) Return the three good alignment axes R ga (n 1 ; n 2 ; n 3 ) = (v 1 ; v 2 ; v 3 ). RR n ¡ 6408 (if N GA (M) = 1) Return n ∈ {v 1 , v 2 , v 3
} the normal of the unique mirror plane as the first good alignment axis and rotate the 3D model in the new coordinate system R (n; v

2 ; v 3 ) if n = v 1 , in R (n; v 3 ; v 1 ) if n = v 2 , in R (n; v 1 ; v 2 ) if n = v 3 . (if N GA (M) = 0)
Return n the normal of the plane with maximal reflection symmetry (cf. section 4.1.2) as the first good alignment axis and rotate the 3D model in a new coordinate system R (n; v 2 ; v 3 ).

5. If N GA (M) ∈ {0, 1}, Compute the direction vector with maximal local translational invariance cost as will be shown in the algorithm of section 4.2.3 and return the three good alignment axes R ga (n; n 2 ; n 3 ).

Plane Reflection Symmetry

There have been two approaches for measuring imperfect symmetry:

-The symmetry distance of a shape with respect to a given symmetry is the minimum mean squared distance from the given shape to its perfectly symmetric shape. This measure has the advantage to estimate the symmetry in 3D surface points. While this distance is precise and robust for measuring symmetry, it is expensive for large models.

-The symmetry descriptor similarity of a shape with respect to a given symmetry is the distance between a shape descriptor of the given shape and that of its perfectly symmetric shape. This measure has been well useful in order to approximate the symmetry distance. The efficiency of the symmetry description in 3D space enables a fast comparison of the amount of reflection symmetries with respect to several planes.

Continuous Symmetry Distance

Let us first define S γ as the reflective surface of S with respect to a plane γ. It is represented by a set of triangles T γ = {T 1 , ..., T NT }, given by a set of vertices P γ = {p 1 , ..., p NP }. Following previous works on distance estimation between 3D surfaces [CRS98, ASCE02] and on symmetry distance [START_REF] Zabrodsky | Symmetry as a continuous feature[END_REF], we define the continuous symmetry distance CSD γ of S with respect to a plane reflection γ as:

CSD γ (S) = 1 A p∈S d(p, S γ ) ds, (1)
where A denotes the area of S and d is the distance between a point p of S and S γ , such that:

d(p, S γ ) = min p ∈Sγ p -p 2 ,
with . 2 is the usual Euclidean norm. The integral of the symmetry error over the whole surface is computed by summing the contributions of all the triangles in T. We obtain a more precise result by taking into account all points of S. The computation of these integrals is only slightly more expensive than the discrete case [START_REF] Zabrodsky | Symmetry as a continuous feature[END_REF]. However, it is necessary to sample the surface S uniformly in order to obtain correct point-S γ distances: each triangular facet is sampled uniformly. The integral over each triangle T i ⊂ T is then approximately done with sums of integrals over triangles obtained by sampling T i .

INRIA

In fact, for each vertices of each sample triangle T ⊂ T i it would be necessary to calculate the distance to all triangles of S in order to find the minimum distance. This leads to a complexity O(N T N S ) that is very expensive for large models. This complexity have been reduced in [START_REF] Cignoni | Metro: Measuring error on simplified surfaces[END_REF] by using a local search processing in order to reduce the number of point-triangle distance evaluations.

The idea is to partition the bounding box into cubic cells and use them in an indexing scheme for the fast search of the nearest triangle of S γ to the sampling point.

If γ is a perfect mirror plane of S, then CSD γ (S) is null. As we want to retain the quasi-perfect mirror planes, we will approximate this definition so that we say that γ is a mirror plane of S when SD γ (S) < ( 0). This test will be used in step 3 of the algorithm described in section 4 in order to select the mirror planes among the coordinate plane normal to a CPCA-axis.

Symmetry Descriptors

The symmetry descriptor represents the symmetries of a given model with respect to several planes in 3D space. It is generally associated to a given shape descriptor that represents a model with a spherical function or a 3D function that rotates with the model. Kazhdan et al. [START_REF] Kazhdan | Symmetry descriptors and 3D shape matching[END_REF] define a symmetry descriptor using the planes through its center of gravity. Podolak et al. [PSG + 06] extended this work by considering symmetries with respect to all possible planes thought a model's bounding volume.

Following [START_REF] Kazhdan | Symmetry descriptors and 3D shape matching[END_REF], we consider a symmetry descriptor that represents the symmetries of a 3D model with respect to only planes through its center of gravity and in the angular neighborhood to the planes normal to the CPCA-axes. We apply our symmetry descriptor to the spherical shape descriptor computed by the Gaussian Euclidean Distance Transform [START_REF] Kazhdan | Symmetry descriptors and 3D shape matching[END_REF] of the surface. Measuring imperfect symmetry is used in step 4 (N GA (M) = 0) of the algorithm described in section 4. Specifically, given the symmetry descriptor values we select the good axis by finding the plane with maximal symmetry.

Local Translational Invariance

Traditionally, in geometry, the translational symmetry is the invariance of an infinite object after a particular translation. We extend this definition to a finite object, in particular, to a 3D model. We define here the local translational symmetry that will be used in this section. This symmetry implies that a 3D model has the same geometrical properties in different parts along a given direction.

Finding the direction that maximizes the local translational invariance is the last step in our general alignment scheme. More precisely, we look for local translational symmetries with respect to all directions perpendicular to the first good alignment axis (cf. section 4). To do this we need to compute a shape descriptor f defined over a one-dimensional interval, that represents a 3D model along a given direction. We need also to define a measure of symmetry for f with respect to local translation along this direction. For this purpose, we describe a method which selects the direction with maximal translational invariance.

Shape Description over 1-D Space

Let d ∈ R 3 be a unit direction vector and π d (ρ), ρ ∈ R, be a family of planes perpendicular to d and at the signed distance ρ from the center of the coordinate system. By taking I d the interval defined by the limits of the 3D surface S in the direction d, we represent S as follows:

S = ρ∈I d S d (ρ),
where S d (ρ) is the 3D sub-shape of S limited by the planes π d (ρδ) and

π d (ρ + δ), δ ∈ R.
In what follows, we consider a function f d defined on the interval I d and having values on a scalar or vector space such that f d (ρ) is a shape descriptor of S d (ρ) for any ρ ∈ I d . Methods computing the shape descriptor f d are given in section 4.2.4.

Local Translational Invariance

LTIC for Alignment

In this section, we investigate the use of the LTIC in 3D to compute a good alignment axis with respect to translational symmetry. More precisely, we want to select the second alignment axis by finding the direction with maximal local translational invariance cost among the directions perpendicular to the first axis n computed in section 4. In order to evaluate the LTIC, we generate a set of direction vectors perpendicular to the first good alignment axis n, obtained by rotating the coordinate system about n: Let R K n be the set generated by the transformation R k which is the rotation about n by the angle θ k = πk K where 0 ≤ k < K:

R k = R θ k (1,0,0) • R , RR n ¡ 6408 with R θ k (1,0,0) =   1 0 0 0 cos(θ k ) -sin(θ k ) 0 sin(θ k ) cos(θ k )  
and R (cf. section 4) is the matrix that contains n in the first row. In what follows, we associate to each R k ∈ R K n one unit direction vector d k equal to the second row of R k . For each d k , a shape descriptor f k is introduced. Now, the problem of computing the good alignment axis is to find the direction d ga or its associated rotation R ga , that maximizes the LT IC(f k ):

R ga = argmax R k ∈R K n LT IC(f k ).
(3)

The algorithm for computing the good alignment axes (given the first one) can be summarized as follows:

Algorithm: Compute direction with maximal LT IC 1. Translate the input 3D model M from its center of gravity to the origin, and scale the translated model such that the average distance of a point on the surface to the new coordinate origin is 1. 3. Return R ga associated to f ga with maximal LT IC.

2. Given a matrix R , for each θ k , (a)
Given a matrix R , this algorithm finds the direction vector with maximal local translational invariance cost. The second good alignment axis n 2 is the direction vector d ga that is the second row of R ga and is perpendicular to the first axis n. The third good alignment axis n 3 is naturally the third row of R ga .

Three shape descriptor models for f k

Suppose the surface S k is positioned in the coordinate system defined by (n, d k , n ∧ d k ), and S k (ρ) and I k are defined as in section 4.2.1 with d = d k . Three shape descriptors models G k , E k and F k are introduced to represent S k . They use only one coordinate (along the axis n ∧ d k ) as the axis n is already selected in the good coordinate system and d k -coordinate is fixed in S k (ρ).

-Global average description G k :

G k (ρ) = 1 A k (ρ) p∈S k (ρ) (n ∧ d k ) T • p ds, INRIA where A k (ρ) = RR p∈S k (ρ) ds denotes the area of S k (ρ).
-Global extremum description E k :

E k (ρ) = " min p∈S k (ρ) ((n ∧ d k ) T • p) , max p∈S k (ρ) ((n ∧ d k ) T • p)
« .

-Vector shape description F k :

Let J k be the interval defined by the limits of the 3D surface S k in the direction n ∧ d k . The bounding box of mesh S k is partitioned into M k cells along the direction n ∧ d k following a regular sampling of J k . Then, S j k (ρ) is the intersection of the shape S k (ρ) and the j th cell.

S k (ρ) = M k j=1 S j k (ρ).
This descriptor represents S k (ρ) with a collection of areas and averages associated to the shapes

S j k (ρ), 1 ≤ j ≤ M k : F k (ρ) = a j k (ρ), g j k (ρ) M k where a j k (ρ) = p∈S j k (ρ) ds if S j k (ρ) = ∅, ∅
otherwise.

g j k (ρ) = 1 a j k (ρ) p∈S j k (ρ) (n ∧ d k ) T • p ds if S j k (ρ) = ∅, ∅ otherwise.

Discrete computation

With the introduced definitions, we deduce a discrete version of the function f k represented on N k points regularly sampled on I k . To define f k at the same scale in any direction d k , the number of samples N k is such that the interval L(I k )

N k has a fixed length 2δ (cf. section 4.2.1) for any orientation k. A unit of measurement N = scale 2δ should be fixed for all 3D models. In our case, N = 1 2δ (N = 32, 64, 128) because M is scaled (see the first step of the algorithm in section 4.2.3).

N k = N L(I k ) . (4)
Similarly, we take M k = N L(J k ) when computing the vector shape description F k .

Finally, for each shape descriptor proposed here, we use a distance dist(f k (i), f k (i )) (where f k (i), 1 ≤ i ≤ N k ) in order to apply definition 1 and choose a normalized error N 0 fixed for all 3D models. We consider that f k

(i) = f k (i ) if dist(f k (i), f k (i )) < N .
In our implementation, dist(, ) is the usual Euclidean

INRIA

norm . 1 for the global average description and the global extremum description; it is defined as follows for the vector shape description:

dist(F k (i), F k (i )) = d F k (i, i ) + d F k (i , i) A k (i) + A k (i ) , ù where d F k (i, i ) = g j k (i) =∅ 1≤j≤M k a j k (i) g j k (i ) =∅ min 1≤j ≤M k g j k (i) -g j k (i ) 1 .
The distance d F k (i, i ) leads to a complexity O(M k M k ). In order to reduce the computing time, we reduce the number of distance evaluations g j k (i)-g j k (i ) 1 , 1 ≤ j, j ≤ M k . We make an a priori coherence assumption: we suppose that the index j m = argmin 1≤j ≤M k g j k (i)g j k (i ) 1 is not far from the index j. Thus, we first test the sparsity of g j k (i ) and we stop if g j k (i ) = ∅. Otherwise, the adjacent points are processed, in order of increasing distance from g j k (i) and we stop when we find g jm k (i ) = ∅. We stress that all not tested points are farther than the found point.

Experimental Results

In order to evaluate our alignment algorithm, we ran the experiments with the Test Princeton 3D Shape Benchmark database [START_REF] Shilane | The Princeton shape benchmark[END_REF] consisted of 907 polygonal models categorized into 92 distinct classes. We found that our approach produces coordinate frames that are robust and semantically correct for most of the models. Figure 5 shows a number of models from different classes aligned by our method. Moreover, our approach provides similar alignments for models belonging to the same class, see for example the alignments of the mailbox class in Figure 6. We evaluated the performance of our alignment method by generating a test set of 33 distinct classes which generally are not well aligned by CPCA. Table 2 gives, for each class, the percentage of perfect alignment (i.e., accurately similar to what a human would select) and compares the results of the CPCA method to our method using the shape descriptors G, E and F introduced in section 4.2.4. As Table 1 indicates, we note that for the three descriptors (G, E and F ), our general scheme provides better alignment performance, with perfect-alignment percentages that are generally close to 100%. When using the shape descriptor F , our method provides more accurate alignment results than using the descriptors G and E.

To measure the efficacy of our alignment algorithm in shape retrieval task, we apply it as a normalization step in a general retrieval process. We use the shape descriptor DLA [START_REF] Chaouch | 3D model retrieval based on depth line descriptor[END_REF] that represents each model by a set of depth lines transformed into sequences and the dynamic programming distance DPD that measures the similarity between the depth line descriptors. To compare objectively the retrieval effectiveness, for both types of alignment methods, we compute Precision-Recall diagrams commonly used in information search (the query is not counted in the answer as in [START_REF] Vranic | 3D Model Retrieval[END_REF]) and four quantitative measures for evaluating query results (see [START_REF] Shilane | The Princeton shape benchmark[END_REF] for a description of this measures):

(1) The Nearest Neighbor (NN), (2) The First Tier (FT), (3) the Second Tier (ST), (4) the Discounted Cumulative Gain (DCG). Comparing the curves as well as the NN, FT, ST and DCG values in Figure 7, we conclude that our alignment method clearly outperforms the CPCA.

Efficiency:

The O(N T ) complexity of the CPCA algorithm makes our approach clearly faster than the existing alignment approaches based on symmetry in 3D rotation space. As you can see in axes. The most time-consuming stage is the symmetry descriptor algorithm for finding the plane with maximal symmetry. This descriptor was computed on 20% of the models that don't have any good alignment axis within CPCA-axes.

v = v π + d v n, v = v π -d v n.
Suppose g is the the center of gravity of S.

We first construct the covariance matrix of S.

C = 1 S v∈S (v -g) • (v -g) T ds = 1 2S v∈S (v -g) • (v -g) T ds + 1 2S v ∈S (v -g) • (v -g) T ds = 1 2S v∈S (v π -g + d v n) • (v π -g + d v n) T ds + 1 2S v∈S (v π -g -d v n) • (v π -g -d v n) T ds = 1 S v∈S (v π -g) • (v π -g) T ds + v∈S d 2 v n • n T ds C • n = 1 S »Z Z v∈S (vπ -g) • (vπ -g) T ds - • n + 1 S »Z Z v∈S d 2 v n • n T ds - • n = 1 S Z Z v∈S (vπ -g) • (vπ -g) T • n ds + 1 S Z Z v∈S d 2 v n • n T • n ds
By previous Lemma g ∈ π and by orthogonal projection for all v onto π i.e., v π ∈ π, we get

(v π -g) T • n = n T • (v π -g) = n T • v π -n T • g = δ -δ = 0
here n is unit vector, thus we have

n T • n = 1
Combining the three equations, we obtain

C • n = 1 S v∈S d 2 v ds n = λ n
Therefore, the unit normal n of the mirror plane π is the eigenvector of S and 1 S v d 2 v is the corresponding eigenvalue.

RR n

¡ 6408

Figure 3 :

 3 Figure 3: Alignments of models with different N GA using our method: Hot air balloon and Hourglass models N GA = 2 + , Chair model N GA = 1, Tree model N GA = 0. Row 1 : CPCA Alignments, Row 2 : Testing the reflection symmetry for CPCA-coordinate planes, Row 3 : Computing the direction with maximal local translational invariance, Row 4 : Our Alignment results.

INRIA

  Cost Definition 1. Given a shape descriptor f d defined over an interval I d and a unit vector d, we say that f d has local translational invariance along d in an interval I ⊂ I d if for all ρ, ρ ∈ I, f d (ρ) = f d (ρ ). In order to measure the local translational symmetry of a shape descriptor f d , we detect the maximal sub-intervals I i of I d such that f d has local translational invariance along d in I i . The cost of this symmetry is defined as follows: Definition 2. Given a shape descriptor f d defined over I d , the local translational invariance cost (LT IC) of f d along d is the sum of the lengths of the maximal intervals I i of I d where f d has local translational invariance along d in I i : LT IC(f d ) = Ii∈I L(I i ) (2) where L(I i ) is the length of I i and I = {I i ⊂ I d |I i maximum; f d has local translational invariance along d in I i }.

  Compute the transformation R k and the associated direction vector d k . (b) Rotate the transformed model (step 1) in the coordinate system R k in order to obtain S k . (c) Compute the interval I k of length L I k and the shape descriptor f k defined over I k . (d) Measure the LT IC(f k ).

Figure 4 :

 4 Figure 4: The three shape descriptors G, E and F corresponding to the chair model are computed along the vertical direction. Surface positions associated to intervals of local translational invariance are colored in blue while the others are colored in red.

Figure 5 :Figure 6 :

 56 Figure 5: Alignments of models from different classes using our method.

  C • n = λ n. Let π = {u ∈ R 3 |n T • u = δ} be the mirror plane of S. Then, ∀ v ∈ S, ∃ (v , v π , d v ) ∈ (S, π, R) such that

Table 1 :

 1 Repartition of 3D models of the Test PSB database with respect of the number of retained CPCA-axes in our alignment method .

Table 2 :

 2 table 1, the CPCA provides, in our general algorithm, a quick alignment for 28.5% of the models that have at least two good alignment Perfect alignment percentages for some classes (311 models). Comparing the Accuracy of CPCA and our method with the shape descriptors G, E and F .

	INRIA

Conclusion

We have presented a new alignment method for 3D models. It retains the principal axes of the CPCA with respect to approximate reflection plane symmetry.

We have introduced a new notion of cost (LT IC) that measures the invariance of a model with respect to local translation about a given direction. This measure is used to compute the remaining alignment axes.

Our experiments show that our approach consistently aligns the 3D objects: we obtain 100% in 24 classes among the 33 classes tested and the others never exceed less than 75% of correct alignment. Moreover, our alignment method provides more accurate results than the CPCA when it is used as a normalization step in a 3D shape retrieval method.

Appendix

Lemma 3. Let π be a mirror plane of S and g be the center of gravity of S. Then g ∈ π .

Lemma 4. Let π be a mirror plane of S and n be the unit normal of π. Then n is an eigenvector of S.

Proof.

The vector n is an eigenvector of the covariance matrix C of S if ∃ λ = 0 such


Shu Yang

Boles Law Mielniczuk

Moulay Saïd

El Youssoufi

Tomasz Hueckel

A note on evolution of pressure and flow within an evaporating capillary bridge

Experiments with evaporation of capillary bridges between two glass spheres show that the bridge gorge radius decreases much faster than the contact radius, distorting the original constant mean curvature bridge shape. In addition, the Laplace pressure calculated from local principal curvatures exhibits high gradients along the bridge moving external surface, most commonly with a high suction near the triple phase contact and positive pressure near the gorge. The high suction results from a negative external curvature at contact. Numerical dynamic simulations with a moving evaporating interface do not currently allow for reproducing a negative external curvature at contact. A series of static simulations are shown based on a representation of an experimentally observed interface, which does include the negative curvature at contact. The resulting Laplace pressure distribution is close to the experimental ones. Most importantly, the pressure gradients induce a consistent flow of liquid from the central area of the bridge, axially toward the solid contact, and then along the solid interface toward the contact area. The flow is believed to contribute to contact pinning. Pinning is viewed as one of the precursors of capillary bridge rupture.

Introduction

The capillarity phenomena are encountered in various areas of mechanical engineering, chemical engineering, geoengineering and energy engineering. They are encountered in people's everyday life, such as rain drop and floating ants. Most of them are too familiar to be given special notice. Among them are capillary bridges forming between wetted solid grains.

The effect of drying on capillary bridges is of critical importance in many new technologies such as laser printing, Atomic Force Microscopy, CO 2 sequestration, to mention a few, but also in traditional technologies, as petroleum engineering, soil, agriculture engineering, for its role in the onset of drying cracking [1][2][3][4][5]. The basic theory of capillary pressure, or surface pressure, or Laplace pressure, that is the pressure difference between the gas phase and the liquid phase separated by a fixed, axially symmetric, curved interface is due to Young [6] and Laplace [START_REF] Laplace | Traité de Mécanique Céleste[END_REF], who separately postulated that the pressure difference is proportional to the mean principal curvature a e-mail: yangshu0813@foxmail.com of the gas/liquid interface. It is commonly expressed (see, e.g., [8]) as follows:

p a -p l = γ 1 r g + 1 r ext , (1)
where p a and p l are air and liquid pressures, r g and r ext are the gorge and the external profile curvature radii, respectively, positive for a convex curve, while γ is surface tension for water/air interface.

In the case of a water bridge between two spherical grains the principal curvatures that are easily measurable are the gorge curvature (always positive) and the external curvature (either negative, positive or zero) at the gorge. Alternatively, pressure difference can be calculated at the triple contact [9,10], where the values of external local principal curvature, as well as the principal curvature of the cross-section of the bridge normal to the bridge tangent plane can be measured, under some reasonable assumptions (as of radial symmetry).

The Young-Laplace equation is valid for static surfaces, with a curved, infinitely thin and immobile interface and hence without account for evaporation [8]. For the former case, and for a weightless liquid, the pressures in water and in the air are uniform, and their difference is constant, and hence through eq. ( 1) the interface mean curvature is deemed to be constant. That implies that the interface surface is one of the following constant mean curvature surfaces identified originally by Delaunay [11]: nodoid, catenoid, unduloid, segment of a sphere, or cylinder [12]. Equation (1) also implies for all the three media forming the bridge being in equilibrium, that the interfaces of separation at their intersection subject to surface tension forces between the media are in equilibrium and thus meet at contact angles that depend only on the surface tension between the pairs of the media. Most commonly the surface tensions are constant and hence contact angles (called hence equilibrium contact angles) are constant. In reality, for sessile drops they are known to depend on temperature, contaminants, and direction of the motion of the contact [13][14][15][16].

Evaporation of the capillary bridge between two spheres brings additional experimentally observed factors into the picture. These are: i) evaporation flux occurring at the interface between the liquid and surrounding gas; ii) the consequent loss of total liquid volume ranging from its initial value to zero (in most cases); iii) possible change of temperature, as evaporation is a weakly endothermic process; however, temperature effects are considered as minor in our considerations; iv) evolution of the shape and position of the liquid/gas interface; stick-slip behavior of the triple interface line between solid, liquid and gas; v) ensuing intermittent changes in the surface area of liquid/solid contact and of its perimeter length; vi) evolution of contact angle (often linked to the contact line behavior); vii) internal liquid flow within the bridge; viii) evolution of the original total attractive force between the spheres; ix) possible mechanical rupture of the bridge (for certain geometric proportions) as a single water body prior to the transition of the whole liquid volume into vapor.

The majority of these phenomena were measured for evaporating bridges and reported (except for iii) and vii) by Mielniczuk et al. [START_REF] Mielniczuk | Geo-Congress 2013[END_REF][START_REF] Mielniczuk | [END_REF][19] (see also Maeda et al. [20]). Temperature changes (iii)), and internal flows (vii)) were either surmised from the general character of vaporization (iii)), or observed, but not measured (vii)). Surprisingly, very little attention has been given to evaporation of liquid bridges, as opposed to extension of liquid bridges [10,21], which was intensely studied since 1970. Luckily, quite a lot of effort has been devoted to the evaporation of sessile droplets [22,23], which brought a substantial progress in our understanding of evaporation of small liquid bodies.

This paper focuses on several features of the behavior of the capillary bridges during evaporation observed and reported in Mielniczuk et al. [19] that differ from the expected patterns of the processes involved. These include:

-initial differences in and evolution of Laplace pressure at different locations of the evolving liquid/gas interface as calculated from the evolving curvatures; -pinning-depinning-repinning behavior of the triple contact line during the process of evaporation; -possible effect of the above on the liquid pressure within the bridge; -the distribution and evolution of evaporation flux at the interface; -the role of the vapor flux and pressure gradients in generating internal flows within the bridge.

Pinning-depinning-repinning have been widely reported, discussed and theorized in the context of both evaporation and displacement of sessile drops [13,24], with two major types of solid surface roughness, nicknamed for their shape as Great Wall of China [25] and fakir [26] /\/\/\, with the liquid in contact with either the solid or gas filling the surface microcavities [27][28][29]. Pinning consists in arresting of, otherwise smooth, re-positioning of the triple phase interface along the solid/liquid boundary. Flows within evaporating drops [23,[30][31][32], and capillary bridges in extension [33] were postulated, simulated and measured for some time. However, the dynamics of capillary bridges induced by evaporation somehow escaped a closer scrutiny.

This paper is prompted by prior observations by Mielniczuk et al. [19] of pinning, depinning and repinning during evaporation of capillary bridges between two smooth glass spheres together with the associated contact angle evolution. Especially, repinning of capillary bridges, when they reach tall bridge proportions calls for attention as a precursor of the spontaneous rupture of the drying bridges, while still containing a substantial volume of liquid. We hence present data from the same experiments regarding the evolution of principal curvatures at different locations of the moving gas/liquid interface, and hence of the calculated Laplace pressures, showing substantial gradients along the interface. We then use these experimental Laplace pressure distributions to validate an approximate numerical procedure to calculate the pressure within the liquid bridge and the associated flow velocity field. The approximation (employing a commercial software) consists of a series of static solutions based on the input of experimental shape and position of the bridge free boundary. At the boundary the interface equilibrium conditions are assumed, and through imposed curvatures also the Laplace pressure at the boundary, which serve as boundary conditions to the pressure field over a boundary segment. The solution hence includes incipient values of Laplace pressure, approximate pressure field across the body of the bridge, and the resulting flow velocity patterns within it. The data set used for the simulations refers to the constant rate of evaporation period of the experiment as documented by Mielniczuk et al. [START_REF] Mielniczuk | [END_REF]. Consequences of the observed patterns are then discussed. Temperature effects, variable surface tension, and specifically Marangoni force effects [32,34], as well as ionized water effect [35] or thin film effects [36] are not considered. And neither are specific micro-scale solid surface roughness models.

Experimentally observed evolution of contact diameter, pinning, curvatures, and Laplace pressure during evaporation

The experiments were conducted on pairs of smooth glass spheres of 8 mm diameter linked by demineralized water bridges at various constant separations exposed to evaporation at a constant temperature of 21 • C and constant environmental relative humidity of 25%. In what follows we analyze the data of the behavior of the water bridge with a grain separation of 1.3 mm and initial water volume of 4 μL. During the test measurements of the total capillary force change between the grains were made using a precision balance, against the volume of water evaporated. Materials and methods are described in detail in Mielniczuk et al. [START_REF] Mielniczuk | [END_REF]. In the previous papers we reported only the changes of curvature at the bridge gorge, as well as the external curvature (of lateral profile), also at the bridge gorge on the basis of what we were able to calculate with the Laplace pressure using eq. (1). In what follows, we report the values of the principal gorge and external radii at four points at different distances from the horizontal symmetry plane of the bridge for 18 states of decreasing volume of water. The principal radii and curvature are located at mutually orthogonal planes crossing the surface of the interface, as shown in fig. 1.

The principal gorge radius of curvature at a point of a bridge is that of a conic cross-section of a cone tangent to the bridge at a point at a given height h, r hp (which may be a circle, ellipse, parabola or hyperbola depending on the considered point height, and hence conic section slope) and is calculated from the corresponding horizontal cross-section radius of a circle, r h . The values of the radii, the corresponding local bridge opening angle, β, and the local value of Laplace pressure are shown for selected instances in fig. 2. In the rest of the paper the assumption is made that axial symmetry is valid as well as that Fig. 2. Profile and geometrical and pressure characteristics of capillary bridges at constant separation of 1.3 mm during evaporation at the relative water content of the original volume of 100%, 69% and 33%, respectively. The values of the radii of curvature and the pressure are calculated as averages for two upper quadrants (left and right) of the projection shown for all other replicas conducted. Radii are in mm, angles are in degrees, and pressure is in Pa. The figures also show, how much faster the motion of the gorge is compared to an almost immobile triple contact line (see also fig. 3). gravity effect does not affect the symmetry with respect to the equatorial plane.

The most important consequence of the occurrence of the evaporation flux at the free boundary, which is the only external change affecting the bridge, is the evolution of the very same free boundary, i.e. the shape of the bridge, of its curvatures, and thus the Laplace pressures at the evaporation interface. By inspection of the evolution in fig. 2 it can be seen that the bridge surface appears to evolve from nodoid through catenoid (not observed in our experiments) to unduloid, to eventually cylinder (not shown here because occurring during an unstable phase), It is visible that the gorge shrinks more than three times faster than the contact. The rupture range shows a relatively wide scatter in terms of the volume evaporated pointing to a possible role of randomness of glass surface imperfections in generating the bridge rupture. The criterion of bridge rapture remains elusive; (b) Laplace pressure vs. relative loss of water volume. At gorge, it is mostly positive and growing as rupture is approached, while at contact, it is predominantly negative, slightly fading near rupture. Also shown is the difference between the gorge and contact pressure, always positive and growing suggesting a push from the gorge toward the contact.

which are all surfaces of constant mean curvatures of Delaunay [19,37,38]. However, there is a visible localized departure from such ideal transformation, consisting in the intermittent evolution of the triple contact line, i.e. the locus of contact of the three phases, marked with red crosses in fig. 2.

Figure 3(a) presents the relative change of contact radius, compared to the relative change in gorge radius. The contact radius decreases altogether by 15% only, whereas the gorge radius decreases by 60%, by the moment the bridge ruptures. In addition, the displacement of contact radius is not continuous, it rather qualifies as stick-slip motion. It needs to be remembered, it is not the displacement of the physical particles that the contact radius measures, rather it is the liquid/gas interface translation, marking the position of particles that are to become gas. Notably, for nearly 1/3 of the volume loss, up to -dV /V 0 = 0.177, the contact radius does not change more than 3.4%, in fact from 1.22 mm to 1.18 mm, what may qualify as pinning, then slips for another 5%, to pin again, while losing next 25% of volume, to slip for about 4% of the original size, and finally to gently come to a halt, after which the bridge ruptures.

Spontaneous rupture of capillary bridges is an interesting phenomenon, as for taller bridges it occurs when only 1/4 to 1/3 of the bridge volume has evaporated. It appears that rupture is preceded by a mechanical instability, linked to a drop in suction, which eventually becomes a positive pressure associated with other forms of geometric transformations [START_REF] Mielniczuk | Geo-Congress 2013[END_REF][START_REF] Mielniczuk | [END_REF].

The evolution of the gorge radius is quite different. Its rate is nearly constant for the first 2/3 of the process, but while the contact slows down toward the end, the gorge accelerates. As the bridge volume is affected by the square of the radii, the aforementioned difference is exacerbated. This is the first indication that the loss of mass is highly non-uniform throughout the bridge, with the center being substantially squeezed, while the contact areas barely changing.

The pinning behavior has been discussed broadly for decades in the context of sessile drops both under evaporation and motion on an incline [27,[START_REF] De Gennes | Gouttes, bulles, perles et ondes[END_REF]. Most of the effort is centered on the contact angle hysteresis, or differences between advancing and retreating contact angles, and explicit microscopic models for the occurrence of the above, including surface roughness or surface tension variations, see e.g. Bormashenko [16]. In the context of capillary bridges, the contact angle has been considered for bridges in extension or compression [21]. Mielniczuk et al. [19] found experimentally that repinning of contact in evaporating bridges correlated with the onset of instability with a substantial drop in suction followed by an increase of Laplace pressure at gorge and a very fast conversion of the unduloidal bridge into a slim cylindrical water stick, and eventual rupture. For tall bridges, which exhibit positive pressure at gorge from the onset, the instability corresponds to a fast increase in gorge pressure and formation of a cylindrical water stick.

The evolution of Laplace pressures at the gorge and at the contact points to a substantial difference in pressure values between the two locations at the initial state. Note that the initial state corresponds to an instant exposure to the relative humidity of the environment, hence instant generation of the evaporation flux. Thus, it can be expected that the initial state does not correspond to the equilibrium situation, which implies no evaporation, hence the uniformly saturated gas environment. In other terms, also the initial Laplace pressure and the interface shape including contact angle are in principle different from those at equilibrium.

Over the course of evaporation, the gorge pressure keeps being positive and increases in value by almost 15 times. However, this is in general agreement with gorge Laplace pressures for other so called slim bridges. Notably, stout bridges initiate their evolution with a negative pressure, which slightly increases, arrives at a maximum of pressure, but shortly before that initialize a departure of the pattern of a continuous growth to then become unstable. At the same time the evolution of the initially stout bridges consists mainly in a reduction of the gorge radius accompanied by a stick-slip behavior of the contact radius, through which the bridge becomes geometrically a slim bridge with the proportions 1:1 of the bridge halfheight to contact radius. In the latter phase, the Laplace negative pressure drastically drops and turns into a positive pressure before either losing all its water through evaporation, or via a water-wire instability [19,[START_REF] Mielniczuk | [END_REF].

The pressure at the contact area is negative (suction) from the very beginning and relatively high (about -20 Pa) and it drops to a minimum of -28 Pa shortly before 30% of water volume loss. It then suffers a sharp turn to a prolonged segment of growth (pressurization of the contact area). In other terms, the initial suction subsides nearly to zero (effectively -4 Pa) prior to rupture. Notably, the end of the initial increase of suction, and the onset of re-pressurization of the contact area coincides with the depinning or the onset of the slip phase. So, the difference between the Laplace pressure value at the gorge and at the contact along the evaporation interface is consistently positive, and growing from about 20 Pa up to nearly 50 Pa. There are few contributions to such developments. First, as reported, the gorge radius is decreasing much faster than the contact radius; second, a continuous pressurizing the center of the bridge and thus finally squeezing of water from the center toward the solid interface and toward the contact. A similar pattern is seen from the numerical simulation of evaporation of sessile drops [30,34].

It must be emphasized that the high suction values are caused by relatively high external curvatures at contact, which are the only possible source of the negative component of Laplace pressure, and at the same time, a substantially decreasing curvature of the circumference at the highly inclined and gradually rotating normal crosssection of the bridge "trunk". It is to be restated that while quantitatively the experimental bridge realizations showed some scatter, mainly because of the way how much symmetrically water would get distributed when placed manually between the grains, the qualitative repeatability was remarkable. Details of the procedure to determine the curvature by fitting Delauney curves are given in Mielniczuk et al. [37].

There are two main possible consequences of the observed (or truly, calculated) significant Laplace pressure gradients: water pressure gradients within the bridge body, and caused by the latter, an internal water flow within the bridge. However, measurements of such quantities were not possible in our experiments. To assess characteristics of both pressures and flow velocities, we have undertaken their approximate numerical simulations, under a set of hypotheses.

Approximate numerical simulations

Theoretical and numerical framework

The behavior of an evolving capillary bridge exposed to a drying atmosphere is quite challenging to simulate numerically. Most of the relevant numerical work has been focused on sessile drops or rise in capillary tubes. Most of work on capillary bridges is analytical (e.g. [12,21,37]) and refers to static bridges.

For an evaporating capillary bridge the problem is stated as a Navier-Stokes flow with a constant viscosity, bounded on one side by a frictional solid surface with a Navier slip law, and a free and moving, evaporating liquid/gas surface, on another, as well as symmetry conditions on the remaining boundaries, see fig. 4(a). In two hallmark problems of dynamic interface, namely moving or drying sessile drop and moving capillary meniscus, the central problem is the initial position, shape and evolution of the free surface of the liquid/gas interface (see, e.g., a comprehensive review of Scardovelli and Zaleski [41]) and its relationship to the triple contact line. The interface is considered as a discontinuity of density, viscosity, pressure and because of the latter, it results to be curvilinear. Physically, it appears to be extremely thin, but numerically, it often is represented as of certain finite thickness. Clearly, the evolution of the interface requires a continuous remeshing of the field. There is a number of problems with numerical simulations of the interface and its evolution. The first of them is a common assumption of a non-slip condition at the liquid/solid interface, which is clearly an untenable assumption for an evaporating liquid surface, with a moving (often in a stick and slip fashion) triple contact line, as pointed out in [42]. Otherwise, using a Volume of Fluid technique implying an explicit slip length results in a convergence failure and mesh size dependence. Diffuse interface algorithm seems to be a good solution for some problems [43], but for our purpose based on interface curvature would require further modifications.

An epic story is that of the contact angle. Clearly, a constant equilibrium value of that entity is not supported experimentally, or to some, not even measurable experimentally (see Bormashenko [16]). For de Gennes et al. [START_REF] De Gennes | Gouttes, bulles, perles et ondes[END_REF] and Ramé and Garoff [44], its evolution is linked to either receding or advancing motion of the triple line. Roughly speaking, a change in contact angle takes place, when contact radius is pinned, and vice versa, the angle does not change when contact is de-pinned. Numerically, a number of contact angle models have been proposed, each with its own deficiency. Namely, the "slip models" lead to a pressure singularity at the contact line, or to a flow kinematics substantially different from those seen in experiments, or both [45], while at no-slip boundary in an analytical solution (no evaporation) of Huh and Scriven [46] shear stress becomes singular approaching the contact. Their solution implies a significant pressure on the free surface, to be balanced by the Laplace pressure at the interface, which requires the interface to be strongly curved (see also Dussan [47]). Concepts of formation/disappearance of new interfaces [48] and/or inclusion of surface heterogeneities at several scales [49,50] require additional experimental data not accessible without additional modeling at molecularand microscales. Among others, that includes data on surface tension in the free surface near the contact line that deviates from the equilibrium ones.

Importantly, it is a practice in dealing with the contact angle to use a piecewise linear representation of the interface near the triple line. That likely is a suitable approximation when the main concern is the flow of the liquid within the corner (see Ramé and Garoff [44], Afkhami et al. [42]). However, our focus is on the Laplace pressure and its evolution toward rupture, and possible precursors of that occurrence, for which a most faithful representation of the interface curvature is crucial. So, piecewise linear approximations with the zero interface-curvature are not suitable for that purpose.

In view of the above-advanced complications, and our limited experimental information, we opt for an approximate approach to simulate the static representation of the deionized water behavior within the capillary bridge. We have used a commercial software COMSOL Multiphysics for simulation. As a characteristic dimension of the bridge we take the separation of 1.3 mm, while surface tension of the standard liquid mixture used in University of Montpellier laboratory is taken as 0.00496 N/m (see [START_REF] Mielniczuk | [END_REF], and [51]). Hence the Bond number Bo 1 and the effects of gravity can be ignored. The flow velocity of liquid is slow enough (Ma 0.3), which means that the gas and liquid can be regarded as incompressible fluid. Temperature is assumed as constant at 294.15 K.

We therefore consider the incompressible fluid flow of the liquid and gas phases as described by Navier-Stokes equations, continuity equations in both phases and the conditions at the liquid/gas interface

ρ ∂v ∂t + (v c • ∇)v = ∇ • ϑ ∇v + (∇v) T -pI + F , (2) ∇ • v = 0, (3)
v c = v(X, t) - ∂x(X m , t) ∂X m ∂X m (X, t) ∂t , ( 4
)
where ρ is the density of fluid, ϑ is the fluid velocity, t is time, I is identity matrix, ν is dynamic viscosity, and F is the volume force set here equal to zero as the gravity is neglected. The extra velocity variable v c is a convective variable, while X, x and X m (X, t), are Lagrangean, Eulerian and Lagrangean moving mesh point coordinates at time t. The latter one is a characteristic of the Arbitrary Lagrange-Euler (ALE) method to describe coordinate points of the mesh to avoid mesh distortion (see e.g. Donea et al. [52]). Boundary and initial conditions are specified as follows: symmetry condition is applied at r = 0, and at z = 0, including zero flow velocity and zero shear stress components in r and z directions, respectively. The linear momentum balance at a generic liquid/gas interface implies that

n 1 T 1 -n 1 T 2 = γ(∇ t • n 1 )n 1 -∇ t γ, ( 5
)
where n 1 is the unit normal vector pointing toward gas phase, T 1 and T 2 are the stress tensors, respectively, in liquid and gas, γ is the surface tension associated with the interface, and ∇ t is the surface gradient. ∇ t • n 1 is twice the mean curvature of the surface, and ∇ t = (I -n 1 n T 1 )∇. The relationship between evaporating mass flux and velocity in liquid and gas phase is described by the Rankine -Hugoniot condition:

v 1 = v 2 + M f 1 ρ 1 - 1 ρ 2 n 1 , (6)
where v 1 and v 2 are the flow velocities in the liquid phase and gas phase at the liquid-gas interface, and ρ 1 and ρ 2 are the density of the liquid and gas, respectively. M f is the evaporating liquid mass flux (see, e.g., Scardovelli and Zaleski [41]). The liquid/solid interface is characterized by the Navier slip, with the resulting friction force

F fr = - ϑ δ v ( 7
)
with δ being the slip-length, taken usually as 1/10 of the grid size. Finally, at the triple line of solid/liquid/gas interface the surface tension force is applied at the contact angle θ c away from the tangent to the solid

F c = γ cos(θ c ). (8)
In the control algorithm of moving mesh, the easiest way to control the mesh motion is to pre-specify the grid motion at the boundary and automatically solve for the internal grid point velocities by using Laplace equation. The influence of the boundary conditions of the mesh on the grid movement is gradually permeating into the interior to obtain grid point variables and updates.

The mesh boundary conditions are shown in fig. 4(b). At the boundaries of the considered field at r = 0 and r = 4 mm, the prescribed mesh displacements are zero in z-direction, while at z = 0 and z = 2.5 mm, the prescribed mesh displacements are zero in r-direction. The prescribed mesh displacement is zero in both r-and z-direction at the surface of the solid sphere. Winslow smoothing is employed for the Laplace equation to update the internal grid (see, e.g., Knupp [53]). Therefore, the effect of mesh boundary conditions can be gently introduced into the internal grid, so that the mesh is smoothly deformed. The velocity of mesh at the liquid-gas interface is then calculated by the equation

v mesh = v 1 • n 1 - M f ρ 1 n 1 . (9)
A separate issue is that of modeling evaporation. Once underwent the phase change, vapor moves via diffusion in gas phases (convection is disregarded), which is described by the equation

∂c ∂t = ∇ • (D∇c), ( 10
)
where c is the concentration of vapor, and D is the diffusion coefficient. The symmetry condition is applied at r = 0 and the horizontal axis of symmetry at z = 0. The transport of diluted species is zero in r-direction and in z-direction, respectively:

-n • (-D∇c) = 0.

At the liquid-gas interface, vapor concentration is assumed to be saturated and an empirical relationship between temperature near the interface T in and saturated pressure P sat can be approximated as [START_REF] Irvine | Steam and Gas Tables with Computer Equations[END_REF] P sat = exp 9.487 -

3.893 × 10 3 T in + 230.47 . ( 12
)
The saturated vapor pressure can be converted to concentration by using the ideal gas theory

c sat = P sat RT in , ( 13
)
where R is the ideal gas constant. At the external boundary of gas phase, we impose vapor concentration as zero, controlled by relative humidity, there, H = 0, hence,

c ∞ = c sat H = 0.
The initial conditions are as follows: the value of the velocity is zero throughout both liquid and gas phases:

V = 0. (14)
The initial values of displacement of the mesh coordinates, R and Z, are zero throughout the liquid body,

dR t=0 = 0, dZ t=0 = 0. ( 15
)
The initial pressure is zero. In this situation, the pressure is the relative pressure of the liquid phase and gas phase:

P t=0 = 0. ( 16
)
The initial value of vapor concentration is set as the vapor concentration of the environment.

c t=0 = c sat H = 0. ( 17
)
The free triangular mesh was chosen for the gas phase as easier to adapt to the boundary than free quad mesh, and the solution is easier to converge. The mapped mesh is built inside the water bridge to avoid the singularity and deformation of the free interface near the liquid-gas interface. The maximum element size of mapped mesh is 0.02 mm. Discretization type is P2+P1, which implies that the velocity components are described by the second-order function, while the pressure field is linear. Figure 5(a) shows the schematic of the mesh for the initial configuration.

Computational strategy

Our effort to employ the COMSOL code to simulate Laplace pressure evolution and flow in a capillary bridge during dynamic motion of the interface coupled with the evaporation flux has been unsuccessful in generating reasonable results. The main reason appears to be a commonly used approximation of the interface near the triple line as piecewise linear in order to simulate well the contact angle. However, that makes it impossible to simulate the negative Laplace pressure component, which can only come from a negative curvature of the interface line, which is excluded by the piecewise linear approximation. Finally, the resultant pressure across the bridge is postulated in COMSOL as constant. This is far from what we measure experimentally. One complication arises that molecular scale physics of the contact area is not well known, with a variety of mechanisms proposed, including formation of a thin continuous film left after receding contact [START_REF] Keita | [END_REF], molecular layering in the liquid phase near the contact, or formation of a single molecular layer [56,57]. Instead of following the dynamic process explicitly, we convert the process in a sequence of static equilibrium states defined by a successive approximated position and shape of the interface as obtained in the experiments. For the given boundary, we obtain a FE solution through COMSOL code in terms of Laplace pressure resulting from the condition of the interface equilibrium (5), the corresponding pressure field in the liquid, and the corresponding flow velocity field. All the above fields correspond to an instant of the onset of flow in response to placing the interface at particular positions. We do not allow the software to proceed and reconstruct the interface as composed of a circle and a straight line near the contact, which we consider un-physical, as it eliminates the possibility of a strong suction at contact (see Dussan [47]).

As evident from the above outline, evaporation, or specifically evaporation flux, is postulated to play no direct role in the process. Clearly, the evaporating flux can be calculated from the diffusion process of vapor from a known fully saturated interface toward completely dry air environment at a remote boundary. Thus we are in the position to associate the evaporation flux with a particular position of the interface. Such a strategy of employing experimentally observed shape and position of the interface has been successfully used in other contexts, such as drying of water in a capillary channel [START_REF] Keita | [END_REF].

Because in the formulation we impose the position and shape of the interface based on experimental observations, our effort loses the meaning of prediction of class A. Rather, it is oriented at explaining certain mechanisms involved. In particular, we are interested in the pressure field and the flow pattern and amounts, which were not measured in our experiments.

Results

Six different configurations, labelled as A-F , of phase transition interface have been examined, for a proper representation, linked to the first phase of de-pinning, middle phase of re-pinning, sliding and final repining, shown in fig. 5(b). They correspond to specific instances characterized by the particular volume of liquid as shown in table 1.

Note that the interface configuration is based on the experimental values of the contact diameter, 2r h and gorge diameter 2r g , between these points a horizontal symmetry axis centered circle was fit. The former of the two points is corrected for a gravity-induced departure from symmetry with respect to the r-axis, until a control on the half-volume is verified. The resulting shape of the interface is not exactly the same as the experimental one, but very close. The driving intention was to allow a negative curvature at the contact point, necessary to allow negative pressures at contact. Except for such a procedure, COMSOL software converts an assumed interface into a (usually smaller) circle with a tangent straight line toward the contact point. As mentioned earlier, the latter approximation arbitrarily eliminates a possible negative pressure at the contact point.

Prior to examining the simulation of the pressures and flow it is helpful to realize how highly non-uniform and non-linear the vapor flux is. Clearly, there is no experimental data to corroborate these results. The classical approach to free evaporation is to assume it as driven by the diffusion of vapor in the air, between the interface with water, where we assume saturation concentration and an ambient level at an assumed distance, as by Sobac and Brutin [START_REF] Sobac | Droplet Wetting and Evaporation: From Pure to Complex Fluids[END_REF]. Deegan et al. [22] and Hu and Larson [30,34] concerning evaporation of a droplet assumed its shape as a spherical cap and found a high concentration of flux near the droplet periphery. In this study, the "Transport of Diluted Species" (TDS) module of COM-SOL Multiphysics is used to obtain the vapor concentration distribution with the form of the liquid/gas interface pre-imposed as the experimental one. At each new configuration the interface is adjusted. As shown in fig. 6(a), the value on the onset near the contact is higher by about 33% compared to average values in the middle of the bridge. The overall flux increases as the process continues, with the concentration near the contact approaching 40% of the values in the middle of the bridge. The Ushaped evaporation flux distribution has been seen both for the evaporating menisci [START_REF] Keita | [END_REF][START_REF] Morris | [END_REF], and for evaporating sessile drops [30,34,57,60]. The "edge enhancement" of the evaporation flux is attributed to a greater probability of an evaporating molecule's volatization at the thin periphery, than when leaving from the center of a drop, or a bridge (see Deegan et al. [60]). Notably, the pinned configurations B, C, D and F yield the highest evaporation flux values at the contact. This coincidence invites more investigations.

Figure 6(b) shows the results of simulations of pressure evolution at gorge and at contact. Notably, the former are positive and the latter ones are negative, as in the experiments. Clearly, they are derivatives of the evolution of the radii and hence, the mean interface curvatures, whereas the interface, while based on the experimental findings, is actually fitted. Figure 7 presents a comparison of the experimental values of both r g and r h , as well as r hext and the corresponding values adopted in simulations. It needs to be remembered that at each configuration r hn is at a different coordinate z and that is projected over a different angle β, so it is not an evolution of the same field variable. The values of the Laplace pressure simulated at contact points as they evolve in the course of evaporation have to be compared to the values obtained experimen-tally shown in fig. 3(b). The simulation values are generally lower (smaller suction) than in the experiments, but relatively levelled, except for the last point, which is a clear outlier (see below). The conclusion: contact is under consistent suction. Although, the two calculations are decoupled in our approach, it is intuitive that water is under the highest suction at the contact area, where the flux of evaporation is highly concentrated, see fig. 6(a).

Figure 7 presents the comparison of the experimental values of both r g and r h , as well as r hext and the corresponding values adopted in simulations. It needs to be remembered that at each configuration r h is at a different coordinate z and that it is projected over a different angle β, so it is not an evolution of the same field variable.

The full comparison of the simulation of Laplace pressure vs. experimental data along the entire interface is presented in fig. 8 for three of the considered configurations: A, C and F. As can be seen for most of the process, the Laplace pressure values obtained through simulations follow the experimental trend and to some degree even the order of magnitude of the experimental data. First of all, in both cases we have a positive pressure in the central region, and negative pressures near the contact line. The simulated data usually exhibit visibly lower values at contact than the experimental ones, while the data closer to the gorge level show a better agreement. The Laplace pressures at the last successful configuration (F, at 1.8 μl) do not follow the previous trend. First, the values of suction near contact (> 50 Pa) are much higher than the experimental values, while near gorge the positive pressures are underestimated by the prediction. The static solution near the contact point has been known for long time to produce infinitely high suction values at the triple line (see Huh and Scriven [46]) and requires a cutoff near the contact (see e.g. Dussan and Davis [61]). Incidentally, configura- tion F, at 1.8 μl was the last one for which the solution could be obtained through COMSOL, and with doubtful results. For further configurations there were no convergent solutions. On the physical side, it needs to be noted that among several realizations of the 4 μl bridge tested, the one shown in fig. 2 [19], is quite significant (∼ 20% of the initial volume). The average volume loss at rupture for 1.3 mm separation bridges is only 40% [START_REF] Mielniczuk | Geo-Congress 2013[END_REF].

The Laplace pressure at the interface can be viewed as boundary conditions activating the dynamics. Using the prescribed experimental shape and position of the interface COMSOL provides an instantaneous pressure field for the entire bridge. Figure 9 presents incipient pressure fields for two of all 6 configurations tested numerically. The most important observation from these results is that the pressure field evolves from an initial configuration (A) where only very central part of the bridge is under positive pressure and the entire part in contact with the solid is under suction, to E, at 2.8 μl (last considered reliable) where almost the entire central part of the bridge is under positive (but smaller) pressure and only a limited ring near the triple contact where the water body is really thin is under suction. This is consistent with our previously published results for the evaporation of the 4 μl bridge at 1.3 mm separation [19].

Notably, for a different realization of the same experiments (there is a clear dispersion of the results due to uneven water volume distribution) the resultant intergranular force component due to pressure acting across the liquid/solid contact is initially repulsive, and small, whereas the resultant intergranular force component due to surface tension acting along the triple gas/liquid/solid contact line is initially attractive and decreases to near zero (ref. [19], fig. 7). In principle, one should not be tempted to directly compare the finite element stress (pressure alone based) representations with the representation which employs a concept of surface tension, which is a line force. It suffices to note that in the obtained FE solution the positive pressures initially decrease to increase toward the second half of the process, and suction follows the same trend: initial decrease followed by an increase. Notably, the geometrical parameters used in the macro-scale theory refer to the total surface area of contact and perimeter length as opposed to a partial (central) area for pressure in FE solution and to a peripheral ring where suction acts for surface tension.

For completeness it is of interest to see the results concerning shear stress (or viscous, for their proportionality to shear strain, with viscosity coefficient) distribution along the most relevant boundaries. For configuration C, fig. 10 presents components of shear stress along the liquid-solid interface, along which Navier frictional conditions (eq. ( 6)) are imposed. Shear stresses are very small over most of the contact.

Visible is a peak of shear stress near the triple line, with a singularity in the r stress component. Notably, a singularity is inherent in the solutions near the contact line (see Huh and Scriven [46]; Snoeijer and Andreotti [50]), unless the notion of Navier condition is revised (see e.g. Gerbeau and Lelievre [49]).

Finally, the flow pattern is shown for both liquids: water and gas, again for the advanced configurations C and E, in fig. 11. For other configurations the pattern is very similar. The flow pattern consists of three areas: an intense, mainly inward radial flow at the gorge of the bridge; a mid-section flow, that is almost entirely axial, i.e. parallel to the bridge axis of rotation and directed to- Laplace pressure calculation is based on an approximation of the experimentally observed position and shape of the interface, and using the liquid/gas interface linear momentum balance condition (eq. ( 5)).

ward the solid contact; and finally, near the solid grain, a predominantly radial outward flow parallel to the substrate surface. The arrows refer to the points positioned at their start. First, it needs to be noticed that the velocities are very small, indeed. Again, it helps to realize that the flow is not induced by evaporation, rather, it is forced by the evolving shape/position of the interface (taken from the experiment) induced by evaporation and/or by the Laplace pressure.

Notably, the flow near the contact line has an opposite direction to the triple line motion.

Discussion and conclusions

The results concerning flow and flow patterns within the bridge, clearly confirm the previous postulates that evaporation may lead to induction of a visible flow within the bridge. While our qualitative experimental attempts were supporting that notion, a comprehensive quantitative experimental study is still missing, for the capillary bridges, while abounds for drying droplets.

The obtained flow velocities are small, but reassuringly of the same order of magnitude, about 0.002 mm/s, as in evaporating sessile drop in the classical result obtained by Hu and Larson [30,34]. Also the pattern of the flow is similar, at least in the near-solid contact area, away from the central evaporating area and toward and along the liquid/solid interface and ending near the triple contact line (fig. 12(a) and (b)). In our case the flow pattern seems to confirm the pressure-driven mechanism. Our approach is to consider the actual interface shape and position, of which mean curvatures control the Laplace pressure. Then the flow within the evaporating bridge can be seen as follows: away from the site of positive Laplace pressure at the central part of the interface, parallel to the horizontal plane of bridge symmetry toward the vertical axis of symmetry, turning vertically toward the grains through a nearzero pressure segment of the bridge, and finally, turning along the grain boundary toward the suction-dominated interface segment near triple contact. Still in simpler terms: the positive pressure caused by the faster shrinking gorge pumps the liquid through the bridge center toward the triple contact perimeter held back by suction.

The above seems also to confirm the leading hypotheses of our numerical approach: that evaporation flux does not directly or strongly affect the flow (as in hypotheses discussed by Sobac and Brutin [START_REF] Sobac | Droplet Wetting and Evaporation: From Pure to Complex Fluids[END_REF] for sessile drops), which is dominated by pressure at the interface. That hypothesis in other terms is equivalent to considering the entire process through a series of snapshots taken for subsequent configurations of the interface taken as incipient dynamic process, without allowing for reaching any steady state. Notably, the simulations conducted this way were convergent until less than half of the liquid evaporated. However, that is quite close to the situation when capillary bridges actually undergo spontaneous rupture [START_REF] Mielniczuk | [END_REF]19,[START_REF] Mielniczuk | [END_REF].

Clearly, all the flow is a result of non-uniform pressure within the bridge, and thus a departure of the bridge shape from the Laplace-Young equilibrium, i.e. constant mean curvature interface surface. We did measure as rigorously as possible the mean curvatures and found them ranging from positive to negative, the latter ones are predominant in the triple contact area (especially for larger separations and taller bridges). The gorge principal curvature and circumferential curvatures above (and below) are always positive. The negative component of the curvature necessarily comes from the external (profile) curvature.

Apart from other intricacies and singularities of the FE modeling of flow with a moving liquid/gas interface, including the moving triple contact point, the introduction of a straight (or actually conical, in 3D) segment of the interface near the contact makes it impossible to re-Fig. 11. Pattern of flow of both liquid within the bridge and gas obtained in simulations at configurations C and E. It can be seen that at configuration C the axial flow in the central part of the bridge is much smaller than in configuration E. This can be related to the fact that in the simulations in configuration E, both pressure near the gorge and suction near the triple contact reach higher values than in configuration C. However, flow velocities for both configurations are of comparable magnitudes near the triple contact, and as a matter of fact, near the entire interface. pressure interface segment at the gorge, parallel to symmetry plane at the bridge mid-height, toward bridge interior, turning vertical to approach the grains, across a near-zero pressure segment of the liquid body, and finally, turning along the grain boundaries toward the suction-dominated interface near the triple contact. In other terms: the faster shrinking gorge exerts the positive pressure, which pushes the liquid through the bridge center toward the triple contact perimeter held back by suction. (b) Schematic of flow in an evaporating sessile drop (from Hu and Larson [30,34]). Similarly to the flow in a lower half of the capillary bridge, the flow is central, downward followed by spreading along the solid/liquid interface toward the air/liquid solid/interface. produce a negative curvature at contact, and hence the negative pressure at that point. This is what lead us to adopt an approach through a sequence of snapshots of imposed interface shapes and positions, with an experiment lead simulation of the negative profile curvature at the contact line. An experiment lead simulation implies that the shape and position of the interface are taken directly from the experiment. This technique appeared to return Laplace pressure values quite close to the experimental values, in particular, realistic values of negative pressure at the contact point. At the advancement of the process quite close to the configuration corresponding to the experimentally observed first instability of the bridge (see Mielniczuk et al. [19]) the numerically obtained suctions became unrealistically high.

Our last comment refers to the role of internal flow in pinning understood as an arrest in the translation of the triple contact line. Despite the huge effort that has been devoted to explaining mechanisms of anchoring of water at the triple contact via solid surface effects in sessile drops, there are no established analytical ways of expressing the role of flow in the positioning of the interface. Hence, the comment is essentially a speculation, hopefully a useful one. In brief, we do observe in simulations a massive displacement of liquid from the central area of the bridge to the contact area. Liquid is being pushed out by a much faster decrease of the gorge radius accompanied by an increasing and accelerating positive pressure in the gorge area. At the same time, water in the contact area is being subject to a high negative pressure.

Notably, in our calculations the pressures in our calculations result from an (as close as possible) approximation of the experimental interface line, in particular, its curvatures. Equally important is to remember that our results showing a high concentration of vapor flux near the contact area, are effectively decoupled from the calculations of flow, and come from a different boundary value problem (vapor diffusion), while still based on the same approximation of the experimental instantaneous position of the interface. In consequence, we are not in the position to use the mass balance of water over an incremental time-elapse, as we do not know what an instantaneous displacement of interface could be over the given amount of time. Nevertheless, we can calculate the total volume of liquid entering the outer ring of the bridge body across a vertical cross-section cut, say at 1/4 of the current contact radius from the contact (which is a somewhat arbitrary measure of liquid flow) which yields for configuration B, an amount of 1.02 * 10 -11 m 3 /s, while the total corresponding volume of water converted into vapor along the interface between the section cut and the entire contact perimeter is 1.19 * 10 -12 m 3 /s, which is almost ten times smaller volume than the entering fluid volume. Thus, the evaporated volume from the considered area is barely 10% of the massive inflow of water pushed into that area due to configuration changes. In other terms, it is concluded that the mass displacements and flow due to Laplace pressure changes are by far the dominant effect. However, at a smaller scale of the tip of the contact, say at 0.0007 mm, that is 0.088% of the current contact radius away from the contact, it is the other way round, i.e. the total inflow of 8.71 * 10 -16 m 3 /s, is smaller than the evaporated volume, which is 3.46 * 10 -14 m 3 /s. At this scale we are near the singularity of the evaporation flux from practically a film. These findings are in agreement with the differences at different scales discussed by Snoeijer and Andreotti [50] for drops. It is useful to realize that the motion of water particles near contact is consistently in the opposite direction of the displacement of the triple contact line. For all practical purposes, at the macroscale, the substantial influx of water into the contact zone may to a large degree contributing to what is perceived as the pinning of the perimeter.

In summary, on the basis of approximate simulation results we propose the following assessment of the understanding of the mechanics of evaporating capillary bridges:

-experimental observations indicate that evaporating water in the capillary bridge is subject to highly differentiated Laplace pressure (calculated from curvatures): with a suction near the contact and positive pressure near the gorge; -the gorge radius decreases much faster than the triple contact radius, the latter one occasionally being so slow, that it appears as pinned; -an accurate simulation of the dynamics of evaporating bridge with a moving evaporating interface requires a precise representation of the (negative) interface curvature near the contact, which is not available in com-mercial numerical codes that use the piecewise linear interface representation; -an approximate approach based on the simulation of a series of patterns of an incipient flow induced by placing the interface using an experimentally observed configuration allows us to calculate a well-predicted Laplace pressure, pressure in the liquid, and liquid flow, as follows: -the observed geometrical disparity between tightening of the gorge, and low mobility of the contact is concomitant with a displacement of mass of water from the central part of the bridge toward the perimeter of the contact. This flow pattern appears to be affiliated with the pinning (or at least substantial slowing) of the contact during drying; -the above is in agreement with the highly non-uniform distribution of the evaporation flux estimated from a vapor diffusion problem along the liquid/gas interface, with a high evaporation flux near the contact and much smaller flux near the center; -the simulations of the evaporation flux on the one hand, and pressures and flow on the other are decoupled; however, both take into account the experimental position of the interface; -the simulations were not possible for advanced states of evaporation known from experiments to approach unstable behavior.

A final conclusion is that to arrive at the ability to predict the flow process induced in a capillary bridge, the computational capability is needed to conduct a rigorous dynamic simulation of the flow linked to the evaporation flux across the freely moving interface, with a better understood pinning kinematics criterion, sensitive to the flow velocity near (or) at the triple point.

TH acknowledges MIST Laboratory (IRSN -CNRS -University of Montpellier) for financial support, which made his several research visits in Montpellier as a visiting professor possible.

Fig. 1 .

 1 Fig.1. Liquid/gas interface surface principal curvatures, r hn and rext at an arbitrary height h of the capillary bridge between two spheres. The principal radii of curvature and curvatures themselves can be determined from 2D profile photos. The principal internal radius is that of a tangent cone crosssecting ellipsoid.

Fig. 3 .

 3 Fig. 3. (a) Evolution of (upper) contact radius and gorge radius vs. relative volume change of water for a bridge at 1.3 mm separation. It is visible that the gorge shrinks more than three times faster than the contact. The rupture range shows a relatively wide scatter in terms of the volume evaporated pointing to a possible role of randomness of glass surface imperfections in generating the bridge rupture. The criterion of bridge rapture remains elusive; (b) Laplace pressure vs. relative loss of water volume.At gorge, it is mostly positive and growing as rupture is approached, while at contact, it is predominantly negative, slightly fading near rupture. Also shown is the difference between the gorge and contact pressure, always positive and growing suggesting a push from the gorge toward the contact.

Fig. 4 .

 4 Fig. 4. (a) Schematic of the capillary water bridge. The solid/gas interface is immobile. The liquid/gas interface is an initially undetermined moving boundary; (b) Mesh boundary conditions. Zero displacement components at the fixed boundaries; satisfying Laplace equation at the mesh points, R and Z in both fluid bodies.

Fig. 5 .

 5 Fig. 5. (a) An example of the initial finite element mesh configuration for the liquid bridge and the surrounding gas, neighboring the undeformable solid grain; (b) Six considered configurations of the evolving phase transition interface (see table 1). The configurations correspond to the experimental ones, as determined from fig. 2. Notably, configurations B, C and D correspond to a pinned contact point. Configuration F is the last one, for which the computations would yield a solution.

  Fig. 5. (a) An example of the initial finite element mesh configuration for the liquid bridge and the surrounding gas, neighboring the undeformable solid grain; (b) Six considered configurations of the evolving phase transition interface (see table 1). The configurations correspond to the experimental ones, as determined from fig. 2. Notably, configurations B, C and D correspond to a pinned contact point. Configuration F is the last one, for which the computations would yield a solution.

Fig. 6 .

 6 Fig. 6. (a) Evaporation flux distribution shown along the vertical coordinates of the free interface of the capillary bridge for the decreasing liquid volumes during evaporation. Note that the flux is consistently about 2.5 times higher near the triple contact line compared to the gorge region. Note also an acceleration of evaporation with time; (b) computed and experimental Laplace pressure values at gorge, p g , and at contact, pc. The frame marks the range within which the solution is considered reliable. Notably, gorge values are positive as in experiments and remarkably close. The contact suction values are also of the order of the experimental ones, but not as close.

Fig. 7 .

 7 Fig. 7. (a) Computed and experimental radii of curvature at the bridge gorge, rg, and at contact, r h evolving with the decreasing volume. The frame shows the range of the simulation validity. The simulation result is better at the gorge than at the contact, the latter is mainly affected by the initial point; (b) Computed and experimental pressures at the bridge gorge, p g , and at contact, p c, evolving at decreasing volume. The frame shows the range of simulations.

  reached the lowest volume of water at rupture. The scatter range of the volumes of water at rupture in the experiments shown in figs. 3(a) and (b), see also Mielniczuk et al.

Fig. 8 .

 8 Fig. 8. Simulated and experimental Laplace pressure distribution along the interface for configuration A, B, . . . , F tested. Both the trend and the values of Laplace pressure are well captured for configurations A, B, C, D, and E. The solution for configuration F exhibits unrealistically high values of Laplace pressure at the contact point and hence is considered as invalid.Laplace pressure calculation is based on an approximation of the experimentally observed position and shape of the interface, and using the liquid/gas interface linear momentum balance condition (eq. (5)).

Fig. 9 .

 9 Fig. 9. Pressure contours for simulated configurations A, and E. Notably, in the initial configuration, A, positive pressure affects almost the entire central half of the water body, while the portions in contact with the solid are under suction.In configuration E the central part of the contact of the grain with liquid is under pressure while only the ring near the triple contact is under suction. Note that the solution for configuration F exhibits irrealistically high values of suction near the triple contact point and hence is considered as invalid.

Fig. 10 .

 10 Fig. 10. Radial component of shear (or viscous) stress along the bridge boundary with the grain for configuration C, with the advanced flow.

Fig. 12 .

 12 Fig. 12.(a) Schematic of flow in an evaporating capillary bridge, summarized as follows: away from the positive Laplace pressure interface segment at the gorge, parallel to symmetry plane at the bridge mid-height, toward bridge interior, turning vertical to approach the grains, across a near-zero pressure segment of the liquid body, and finally, turning along the grain boundaries toward the suction-dominated interface near the triple contact. In other terms: the faster shrinking gorge exerts the positive pressure, which pushes the liquid through the bridge center toward the triple contact perimeter held back by suction. (b) Schematic of flow in an evaporating sessile drop (from Hu and Larson[30,34]). Similarly to the flow in a lower half of the capillary bridge, the flow is central, downward followed by spreading along the solid/liquid interface toward the air/liquid solid/interface.

Table 1 .

 1 Selected configurations of gas/liquid interface for the indicated current water volume and evaporated fraction of the original water volume.

	Configuration code	A	B	C	D	E	F
	Current water volume	4 μl	3.3 μl	3.1 μl	3.0 μl	2.8 μl	1.8 μl
	Fraction of the volume evaporated	0	0.175	0.225	0.25	0.3	0.55

Author contribution statement

SY performed all analyses, calculations and produced numerical results and their theoretical assessment. BM and MSE-Y contributed to building the experimental setup. BM performed all the experiments and processed experimental data. BM and MSE-Y contributed to the experimental data processing and analyses. TH contributed to the test design, assessment and analyses of the experimental data, calculations and the numerical results. All contributed to the writing of the manuscript.