in $L_\infty$-norm. On the other hand for the LBO-Essad. method, the $L_2$-order of convergence is much smaller than 1 and the $L_\infty$-norm saturates for $d \times \max_{j \in [1,N]} (\kappa(x_j)) \le 0.2$. This first comparison indicates that the discretization of the mixed area $A_{mixed}$ for evaluating the discrete Laplace-Beltrami Operator is of paramount importance when the triangles are obtuse. Whereas the curvature calculation is very sensitive to the chosen discrete method, the evaluation of the unit normal vector is not (Fig. 3.6b). The $L_2$-order of convergence is nearly 3, like, in a lesser extend, the convergence rate in $L_\infty$-norm.

With disturbances of 20% the initial edge size $d$, the shape of the triangles wander strongly away from the equilateral structure. Figure 3.7a shows that both methods saturate at the same level with the grid refinement when the curvature is considered. Despite this lack of convergence,

<figure><img src="image_41.png" /><figcaption>FIGURE 3.7 – LBO methods - Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$ and $\mathbf{n}$, $p = 20\%$</figcaption></figure>

the unit normal vector continues to converge, for both methods, but at the first order instead of the third order measured for regular meshes.

## IF-methods

For the Integral Formulation methods, namely the IF-Trygg. [123] and IF-triang. methods, the same analysis as the one proposed for the LBO-methods, is performed for a regular mesh ($p = 0$) and a disturbed one ($p = 0.2$).

For regular meshes ($p = 0$), the numerical curvatures converge with the decrease in the dimensionless grid size $d \times \max_{j \in [1,N]} (\kappa(x_j))$ (Fig. (3.8a)). The accuracy of the IF-Trygg. method is second order for both norms. The results are worse for the IF-triang. method. Although the $L_2$-order of convergence is equal to 2 for coarse grids, it decreases to 1.5 asymptotically. The behaviour of the $L_\infty$-norm, firstly second order, becomes first order for the fine meshes. Thus, for the finest grid the maximal relative error for the IF-triang. method is larger than 2 decades to that of the IF-Trygg. one. The numerical approximation of the normal unit vector is slightly better for the IF-Trygg. method than for the IF-triang. method (Fig. (3.8b)). However, they converge both with a second order accuracy for the $L_2$ and $L_\infty$-norms. This first set of tests seems to indicate that the IF-Trygg. method is more efficient than the IF-triang. method.

As for the LBO-methods, the curvature approximation is deeply worsen as soon as 20% disturbances are applied to the mesh (Fig. 3.9a). Again, the relative errors remains constant, whatever the grid size is. In that case, the IF-triang. method seems to be slightly better, about twice smaller than the IF-Trygg. method. The unit normal vector approximations keep convergent behaviours, with a first order slope, whatever the norm and the method (Fig. 3.9b). As for the
<figure><img src="image_42.png" /><figcaption>FIGURE 3.8 – Integral Formulation - Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$ and $\mathbf{n}$, $p = 0\%$</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>FIGURE 3.9 – Integral Formulation - Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$ and $\mathbf{n}$, $p = 20\%$</figcaption></figure>

curvature, the IF-triang. method is better, but the difference is not really significant.

## SR-methods

The last family of methods concerns the Surface Reconstruction techniques. Several approaches are compared in this section by using the notations presented in Sec. 3.2.4. Some are based on the first vertex neighbourhood, the SR-$\mathcal{N}_1$-$\mathbf{n}_{iter}-0$ [138], SR-$\mathcal{N}_1$-$\mathbf{n}_m$-0 [26], SR-$\mathcal{N}_1$-$\mathbf{n}_c$-0, SR-$\mathcal{N}_1$-$\mathbf{n}_m$-$a_{00}$ and SR-$\mathcal{N}_1$-$\mathbf{n}_c$-$a_{00}$ methods, and others are expressed with the second neighbourhood, the SR-$\mathcal{N}_2$-$\mathbf{n}_m$-0, SR-$\mathcal{N}_2$-$\mathbf{n}_c$-0, SR-$\mathcal{N}_2$-$\mathbf{n}_m$-$a_{00}$ and SR-$\mathcal{N}_2$-$\mathbf{n}_c$-$a_{00}$ methods. Following the same methodology as before, the accuracy of the approximations for the curvature and the unit normal vector are carried out as a function of the grid size, in a first step for a regular mesh, and then for a 20% disturbed one.

Figure 3.10 illustrates the convergence of the curvature $\kappa_i$, for a regular grid ($p=0$), of the different methods based on the first neighbourhood $\mathcal{N}_1$ (Fig. 3.10a) and second neighbourhood $\mathcal{N}_2$ (Fig. 3.10b). Whatever the norm, the methods are second order accurate. However, the best results are obtained with methods that enforce the vertex $x_i$ to belong to the surface (SR-$\mathcal{N}_k$-$\mathbf{n}_\alpha$-0, with $k \in \{1, 2\}$, and $\alpha \in \{c, m, iter\}$). Likewise, the increase in the stencil from $\mathcal{N}_1$ to $\mathcal{N}_2$ for the reconstructed surface seems to cause a slight deterioration in the approximation. The same remarks apply when the relative errors on the normal vector are studied (Fig. 3.11).
<figure><img src="image_44.png" /><figcaption>FIGURE 3.10 – Surface Reconstruction - Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$, $p = 0\%$</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>FIGURE 3.11 – Surface Reconstruction - Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\mathbf{n}$, $p = 0\%$</figcaption></figure>

When the triangular mesh is no longer regular ($p = 20\%$), the curvature convergence is maintained, but with a reduced convergence order equal to one (Fig. 3.12). Again, enforcing the

<figure><img src="image_46.png" /><figcaption>FIGURE 3.12 – Surface Reconstruction - Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$, $p = 20\%$</figcaption></figure>

vertex $x_i$ to lay on the surface provides slightly better results. On the other hand, the relative errors on the normal vectors are not sensitive to the perturbation magnitude (Fig. 3.13). Indeed, both the error levels and convergence orders are preserved for disturbed meshes.
<figure><img src="image_47.png" /><figcaption>FIGURE 3.13 – Surface Reconstruction - Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\mathbf{n}$, $p = 20\%$</figcaption></figure>

### 3.3.2 Comparisons

In this section, the best methods of each family of discretization for the curvature and normal vector are now compared. This comparison step is performed for $p = 20\%$ which leads to the largest errors, but also for $p = 0\%$. Indeed, one may use (or imagine) smoothing procedures to recover a nearly regular mesh which would be close to what is obtained for $p = 0\%$. Thus, for the Laplace-Beltrami-Operator discretization, only the LBO-Meyer method [76] is considered. Concerning the discrete Integral Formulation, both the IF-Trygg. [123] and IF-Triang. are retained. For the last family of methods, the best ones are those satisfying that the vertex $x_i$, where the curvature and normal vector are calculated, belongs to the reconstructed surface. Thus, both the SR-$\mathcal{N}_1$-$\mathbf{n}_m$-0 [26] and SR-$\mathcal{N}_2$-$\mathbf{n}_m$-0 methods are kept, depending on the size of the neighbourhood for the calculation of the approximate surface.

#### Errors on curvature and normal vector

For the curvature calculation on a regular grid, although all the approaches exhibit convergence behaviours, the SR-$\mathcal{N}_1$-$\mathbf{n}_m$-0 method provides the best results up to nearly two order smaller than those obtained with IF-Triang. when the $L_2$ norm is considered. This gap is increased up to three decades when the $L_\infty$ norm is chosen (Fig. 3.14a).

<figure><img src="image_48.png" /><figcaption>FIGURE 3.14 – Synthesis - Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$ and $\mathbf{n}$, $p = 0\%$</figcaption></figure>

The differences between the methods are much smaller when the analysis is focused on the unit normal vector (Fig. 3.14b). The most interesting result is probably the third order conver-
gence of the LBO-Meyer method, which provides, on the one hand a quite bad approximation for coarse meshes and on the other hand best results for the finest meshes.

For the curvature calculation with 20%-magnitude perturbations, only the reconstruction methods SR-$\mathcal{N}_1$-$\mathbf{n}_m$-0 and SR-$\mathcal{N}_2$-$\mathbf{n}_m$-0 keep going to converge with the mesh refinement, but with a rate falling from second to first order for both norms (Fig. 3.15a); the other methods saturate. It is worth noticing that, even for small random perturbations, the quality of the cur-

<figure><img src="image_49.png" /><figcaption>FIGURE 3.15 – Synthesis - Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$ and $\mathbf{n}$, $p = 20\%$</figcaption></figure>

vature approximations is highly deteriorated (Fig. 3.16). Indeed, for a given dimensionless mesh

<figure><img src="image_50.png" /><figcaption>FIGURE 3.16 – Variation of the relative errors as a function of the maximum amplitude $p$ of the perturbations</figcaption></figure>

size $d \times \max(\kappa) = 1.22 \times 10^{-2}$, adding tiny random disturbances of maximal magnitude $p = 5\%$ increases the error induced by Reconstruction Surface methods by at most one decade, nearly two decades for the IF-triang. method and three decades for LBO-Meyer and IF-Trygg. methods The use of irregular meshes for the normal evaluation reduces the order of convergence to one, except for the Surface Reconstruction methods which are quite insensitive to the perturbations (Fig. 3.15b).

A partial conclusion on the accuracy of the different numerical approximations can now be drawn. Contrary to the LBO and Integral Formulation methods, the Surface Reconstruction seems promising for the curvature calculation since it preserves convergence properties, even for disturbed meshes. Indeed, such meshes are very common : they can for example be met in front tracking methods for which the interface between two immiscible fluids are discretized by a set of triangles whose vertices are advected by the flow field. Although the accuracy of the normal
vector approximation is better for the Surface Reconstruction methods, the LBO and Integral Formulation remains attractive, despite much larger errors, especially with the $L_\infty$-norm.

## CPU-cost and efficiency

Although the Surface Reconstruction methods provide very often the best approximations, their computational times are 3 times larger than that stemming from the LBO-Meyer approximation, and up to 10 times larger when compared to the IF-Trygg. method for the finest meshes. Thus, to perform pertinent comparisons, the accuracy of the different numerical methods must be studied in regards of the CPU time required to get them. The machine used is an Intel(R) Xeon (R) CPU E7-8890 v4 @ 2,20Ghz.

Whereas the Laplace Beltrami Operator and Integral Formulation methods compute the surface tension, from which the curvature and the unit normal vector are directly deduced, the Surface Reconstruction method needs at first the evaluation of an analytical approximation of the discrete surface, then the independent calculation of the curvature and the normal vector to finally evaluate the surface tension. The expensive part of these algorithms lies in the surface tension calculation for the LBO and IF methods, and in the Surface Reconstruction for the SR method. Therefore, the computation times for the evaluations of the curvature, the unit normal vector, the surface tension or the three together are similar.

Figures 3.17 illustrates the relative errors committed on the curvature approximations as a function of the CPU time spent to evaluate $\kappa_i$, $\mathbf{n}_i$ and $\mathbf{f}_i$ over the whole vertices of the mesh. It must be noticed that the increase in the CPU time corresponds obviously to the decrease in the mesh size. For $p = 0\%$, the most efficient methods are IF-Trygg. and SR-$\mathcal{N}_1$-$\mathbf{n}_m$-$0$ (Fig. 3.17a). Assuming a relative error of 1%, the aforementioned methods take about 0.02 s and require a

<figure><img src="image_51.png" /><figcaption>FIGURE 3.17 – Relative errors on $\kappa$ in function of the time in $L_2$ (continuous line) and $L_\infty$ (dash line) norms, with (a) $p = 0\%$ and (b) $p = 20\%$.</figcaption></figure>

mesh size $d \times \max_i(\kappa(x_i)) = 0.3$ for IF-Trygg and 0.5 for SR-$\mathcal{N}_1$-$\mathbf{n}_m$. But for $p = 20\%$, this accuracy can only be reached by the Surface Reconstruction methods, for a CPU time about 5 times larger than with $p = 0\%$ and a grid size $d \times \max_i(\kappa(x_i)) = 0.2$ (Fig. 3.17b).

### 3.3.3 Summary

For stationary surfaces, the evaluation of the curvature, which behaves as the surface tension, is the most sensitive step. Only the Surface Reconstruction methods prove to be convergent and accurate with the mesh refinement, whatever the regularity the surface mesh and the applied norm. The Laplace-Beltrami Operator and Integral Formulation methods are not robust
regarding the mesh regularity. And even for regular grids where they provide good results, their efficiencies are never better than the Surface Reconstruction methods.

These tests for static surfaces give us a first glimpse of the accuracy and robustness of the different approaches to get the surface properties. The next section is focused on the coupling between the surface advection and the computation of the curvature and normal vector.

## 3.4 Analysis of advected surfaces

In the previous section, the calculation of the curvature and normal vector estimates are carried out for a fixed discrete surface: the resulting errors are then only governed by the mesh size, the regularity of the triangles and the accuracy of the approximation methods. Such an approach has already been used in literature [67, 128, 132]. The only drawback of this static approach is that it does not take into account all the processes involved in the simulation of moving or deforming interfaces.

The aim of this part is to study, in what extend the advection of the interface mesh with all associated numerical procedures (advection, reseeding, swapping) may alter the results, both in the convergence process and in the accuracy.

### 3.4.1 Problem presentation

To address this problem, an analytical flow field acting on a moving interface, with deformations, is proposed in the first part of this section. Thereafter, the mechanisms used to advect and manage the interface mesh deformation are presented. These supplementary processes must also be taken into account for the definition of the reference value, which forms the third part. Later on, the methodology of the numerical tests is presented.

#### Velocity field

In order to design a dynamic interface with deformation, a divergence-free flow field defined in spherical coordinate frame is first considered

$$ \mathbf{v} = \frac{Q}{r^2} \mathbf{e}_r, \quad r > 0 \qquad (3.24) $$

with $4\pi Q$ the algebraic volume flow rate through the sphere of radius $r$. Starting from a discrete surface located at $z = z_0$, each vertex $x_i$, initially at coordinates $\mathbf{x}(t=0) = (r(t=0), \theta, \varphi)$, is advected in the radial direction according to equation:

$$ r(t) = \sqrt[3]{3Qt + \frac{z_0^3}{\cos^3\theta}} \qquad (3.25) $$

Depending on the sign of parameter $Q$, the origin of coordinates corresponds to a sink or a source. If $Q < 0$, the surface is stretched to the singularity, forming an elongated ligament characterized by two principal curvatures with opposite signs. On the contrary, for $Q > 0$, the surface inflates, producing an increasing bubble with a sharp lateral edge, similar to an impact on a liquid film. The interface topology in this configuration also looks like a droplet impinging a flat liquid film. These two cases are treated in separated subsections, as they clearly define quite different interface deformation dynamics, with curvature singularities of liquid finger type for the case $Q < 0$ whereas the other case provides an abrupt curvature evolution of liquid impact type in the most varying interface geometry zones.
## Surface advection

The mesh of the initial surface is made up of equilateral triangles with an edge size $d$. The vertices of these elements are then randomly displaced with a maximal disturbance $p$, according to relation (3.20). Advection from $t = n\Delta t$ to $(n+1)\Delta t$ of each vertex $x_i$ of the mesh is the result of the numerical integration of equation $dx_i/dt = V_i$ by a second order Runge-Kutta scheme (Heun's method). The vertex velocity $V_i$ is spatially interpolated by the PERM method [74] which relies on the knowledge of the exact velocity components (Eq. 3.24) on the staggered Eulerian grids. The fix Eulerian grid size is always twice the edge length $d$ of the equilateral triangles, before applying any perturbation. The time step $\Delta t$ is based on the Courant number equal to 0.5. Using the maximal instantaneous velocity magnitude over all the vertices $x_i$, namely $\max_i(Q/r_i^2) = Q(3Qt + z_0^3)^{-2/3}$, and the Eulerian grid size $2d$, the instantaneous time step writes $\Delta t = (3Qt + z_0^3)^{2/3} \times d/Q$.

During the transport, the edges of the triangles are stretched or compressed, leading eventually to undesirable elements. To overcome this issue and in order to keep triangle edges in the range $[d/2, 2d]$, several procedures are implemented [99]. A new vertex may be inserted at the middle of a triangle edge, if its size exceeds $2d$. Conversely, if the triangle edge is less than $d/2$, it is suppressed by merging its two vertices at the center of the former edge. At last, to improve the quality of the triangles, edges may be flipped between two flatten triangles sharing their longest edge. All these procedures, which aim to keep a good approximation of the real interface, are performed at each time step (see Chap. 2 for more details). Concerning the transport of the interface mesh, both the exact advection (the edge velocity is not interpolated but uses the analytical flow field Eq. 3.25) and the discrete advection scheme described before, as in real multiphase flow simulations, are performed. Indeed, considering an exact advection is interesting here in order to separate several numerical effects linked to interface management, advection or geometrical properties.

The error calculations for the moving surface are restricted to $\theta \le \theta_0$, with $\theta_0 = \arctan(L_0/z_0)$, $L_0 = 1$, $z_0 = 1$ for $Q < 0$ and $z_0 = 0.01$ for $Q > 0$. The extension of the initial flat surface is large enough to avoid edge effects during curvature and normal vector computations.

## Reference surface

In Sec. 3.3, the mesh vertices were exactly located on the stationary surface. For advected surfaces, the time integration, as well as the different procedure developed to improve the mesh quality, do not ensure that the vertices belong to the analytical surface any more after exact or discrete advection. Therefore, to perform comparisons and error measurements (3.21) (3.22), a reference surface must be defined; this issue is addressed below.

Two different cases are met, depending on the exact or approximative nature of the vertex transport.

For an exact transport, the advected vertices are located on the analytical surface (3.25). Due to the mesh improvement procedures, the modified or new nodes are projected on the reference surface in order that all the vertices keep an exact position as a function of time. Therefore, reference values are simply calculated with relations (3.18) and (3.19).

On the other hand, the advection of the vertices by a numerical scheme (Sec. 3.4.1) does not allow to stay on the analytical surface at any time. To compare the curvature and the normal vector to reference solutions, each vertex $x_i(t^{(n)})$ of the mesh must be associated to a point $x_i^{\text{ref}}(n\Delta t)$ located on the analytical surface $S^{\text{ref}}(n\Delta t)$. This point, defined as the intersection between this surface and the streamline stemming from the considered vertex $x_i(t^{(n)})$, is obtained as follows (Fig. 3.18). The vertex $x_i(t^{(n)})$ is first advected with the exact transport equation
<figure><img src="image_52.png" /><figcaption>FIGURE 3.18 – Choice of the reference value to compute the relative errors.</figcaption></figure>

backward in time to give a point $x_i^{\text{ref}}(0)$ belonging to the initial flat surface at $z = z_0$. Since the exact and approximative trajectories differ, the coordinates $\mathbf{x}_i^{\text{ref}}(0)$ and $\mathbf{x}_i(t^{(0)})$ (assuming the node $x_i$ existed at $t=0$ in the mesh) are different from each other. Therefore, the time required to move the vertex $x_i$ from $\mathbf{x}_i(t^{(n)})$ to $\mathbf{x}_i^{\text{ref}}(0)$, is not $n\Delta t$. The reference point $x_i^{\text{ref}}(n\Delta t)$ is then simply obtained by advecting forward the point $x_i^{\text{ref}}(0)$ with the exact transport equation over the time interval $n\Delta t$. Finally, the reference curvature and the normal vector are computed at coordinates $\mathbf{x}_i^{\text{ref}}(n\Delta t)$ with relations (3.18) and (3.19) applied to the analytical surface $S^{\text{ref}}(n\Delta t)$.

## Methodology of the numerical tests

Following the same methodology as in the previous section, the LBO-Meyer [76], IF-Trygg. [123], IF-Triang., and the SR-$\mathcal{N}_1$-$\mathbf{n}_m$-0 [26] methods are considered. In addition, variants of the Surface Reconstruction method are studied for larger stencils ($\mathcal{N}_2$ instead of $\mathcal{N}_1$) or by keeping the degree of freedom $a_{00}$ in the evaluation of the approximated surface $Z = f(X, Y)$ (Eq. 3.14). For all these methods, the tracking of the interface as a function of time is carried out for an initial mesh made up of, either perfect equilateral triangles, or irregular patterns (see. Eq. 3.20).

To point out the role played by the advection scheme on the curvature and normal vector evaluation, simulations are carried out into two steps.

First, the displacement of the vertices is performed in an exact way in order to get rid of errors associated to the transport scheme. Thus, the difference between the numerical and exact solutions can only be due to the mesh structure, namely the spatial distribution of the triangles and their departure from the ideal equilateral shape. Contrary to the stationary surfaces, whose vertices always belong to six triangles, the addition or suppression of vertices/triangles during the advection process modify the number of neighbours.

In a second step, the transport scheme is accounted for, in order to estimate how much it reduces the accuracy of the numerical approximations for the curvature and normal vector.

### 3.4.2 Ligament formation – Case $Q < 0$

For a sink located at the origin $O(r = 0)$, the parameter $Q$ is negative. The front tracking surface, initially at $z = z_0$, is sucked all the more that the vertices are close to the singularity. The resulting surface exhibits a thin tail coming closer and closer to $O$ as time runs out. It is interesting to notice that, as can be observed in a two-phase flow, the evolution of the surface
looks like a long liquid ligament which becomes thinner and thinner as a function of time. In real two-phase flow, the surface tension instabilities would potentially break up the structure to produce droplets.

In the next subsections, the initial flat surface is located at $z_0 = 1$ and the volume flow rate is set to $4\pi Q = -0.4\pi/3$. The calculations are carried out until time $t = 9$ s. Figure 3.19 illustrates the shape of the surface at the final instant, for vertices advected in an exact way.

<figure><img src="image_53.png" /><figcaption>FIGURE 3.19 – Bottom and side views of the discrete mesh transported in a spherical velocity field, $Q < 0$</figcaption></figure>

## Exact transport – Ideal case

In the ideal case, vertices added by the remeshing process are projected back on the analytical surface and simulations are carried out with the exact transport model, starting from the surface at $z = z_0$ discretized by equilateral triangles (see Eq. 3.25).

<figure><img src="image_54.png" /><figcaption>FIGURE 3.20 – Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms for $\kappa$ and $\mathbf{n}$, $p = 0\%$ with exact transport - $Q < 0$</figcaption></figure>

The behaviors of the curvature error as a function of the reduced mesh size (Fig. 3.20) show convergences in $L_2$-norm, but with very different orders depending whether IF or LBO methods are used (0.3, 0.34, 0.47 respectively for LBO-Meyer, IF-Trygg, and IF-Triang.) or the SR approximations are considered (nearly second order). This distinction is also true for the $L_\infty$-norm : whereas the asymptotic convergence is close to 2 for the SR-methods, a saturation is observed for the two other classes of methods. Concerning the normal vector, all the methods
are convergent with a second order slope in $L_2$-norm, order that is preserved in $L_\infty$-norm with the SR methods and reduced to first order for the IF and LBO-methods.

Introducing 20%-perturbations in the initial mesh destroys the convergence of the curvature (Fig. C.1 in C.1) for the LBO and IF approximations in $L_2$-norm, and even exhibits a divergence behavior in $L_\infty$-norm. This deterioration of accuracy is also measured for the SR solutions where the order of convergence falls from 2 to 1, whatever the norm considered. The normal vector approximations are also getting worst for the LBO and IF methods, decreasing from 2 to 1, whereas their orders are unchanged for the SR approximations.

## Numerical scheme for the advection

Unlike the ideal case, the new vertices from the remeshing process are not projected back on the analytical surface and the advection is a second order Runge-Kutta scheme (Heun's method). The vertex velocity is interpolated by the PERM method [74].

Figure 3.21 gives the relative errors for the curvature and the normal vector, when the initial flat surface is made of equilateral triangles and an advection of the vertices is carried out by the numerical scheme. For the curvature evaluation, the substitution of the exact transport

<figure><img src="image_55.png" /><figcaption>FIGURE 3.21 – Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms for $\kappa$ and $\mathbf{n}$, $p = 0\%$ with Lagrangian transport - $Q < 0$</figcaption></figure>

by a numerical scheme increases the level of the errors for all methods (Figs. 3.20a and 3.21a). Moreover, the SR-$\mathcal{N}_1$-$\mathbf{n}_m$-$0$ [26] approximation loses precision to become as accurate as the LBO-Meyer [76] solution, which is much cheaper in computation time than SR-$\mathcal{N}_1$-$\mathbf{n}_m$-$0$. To recover more accurate solutions, the number of neighbouring triangles utilized to define the reconstructed surface is increased from $\mathcal{N}_1$ to $\mathcal{N}_2$, and the degree of freedom $a_{00}$ of the analytical surface $Z = f(X, Y)$ (Eq. 3.14) is let free in the fitting process. These changes lead to the SR-$\mathcal{N}_2$-$\mathbf{n}_m$-$a_{00}$ method. The use of the numerical scheme for the vertex advection also deteriorates the order of convergence. As discussed here-above, the improvement in the quality of the approximations is the result of the increase in the stencil size with, for the Surface Reconstruction methods, the addition of one degree of freedom to compute the surface fitting. These two conditions ensure a smoothing effect of the surface irregularities, which originates from local remeshing performed to preserve a good mesh quality (see Sec. 3.4.1 for details).

For the normal vector calculation, the level of the errors is higher than the one obtained with the exact transport and the order of convergence is 1 (Figs. 3.20b and 3.21b).

A quick insight on results initialized with the flat surface made of equilateral triangles, whose vertices are randomly perturbed with maximum magnitude equal to a 20% edge length (see Fig. C.2 in C.1), shows that no significant difference appears between solutions obtained with $p = 0\%$ and $p = 20\%$. Thus, contrary to the exact advection case where the orders of
convergence were substantially reduced with the loss of mesh regularity, this is not observed
here because the transport scheme introduces errors that dominate those produced by the lack
of mesh regularity.

CPU-cost and efficiency

The choice of a numerical method results from an equilibrium between the accuracy it
provides and the computation time required to get the solution. Figure 3.22 shows the relative

<figure><img src="image_56.png" /><figcaption>FIGURE 3.22 – Relative errors on $κ$ (a) and $n$ (b) according to CPU time in $L_2$ (continuous line) and $L_∞$ (dash line) norms, $p = 0\% - Q < 0$ and numerical transport.</figcaption></figure>

errors for the curvature and normal vector as a function of the CPU-cost, for an initial flat
surface made up of undisturbed ($p = 0\%$) triangles advected by the numerical scheme. It has to
be noticed that, for a given method and mesh, the computation time is identical between the
curvature and the normal vector because both approximations are highly linked together (see
eq. 3.10 for LBO- and IF-methods, and Eqs. 3.18 and 3.19 for the SR-methods).

For the curvature evaluation, the SR-$N_2-n_m-a_{00}$ and IF-Triang. methods are the most CPU-
time efficient as soon as the relative error in L2-norm is less than 20%. In that case, the afore-
mentioned methods are 2, even 3 times more accurate for the same CPU cost. Their superiority
is less clear when the normal vector has to be approximated (Fig. 3.22b). Indeed, unlike for
the curvature, no method is really distinguishable from the others, except maybe the IF-Trygg.
method [123] which goes a little bit faster, mainly when a coarse approximation is required. The
computation times for disturbed initial meshes with $p = 20\%$ provide similar conclusions (C.1,
Fig. C.3).

3.4.3 Bubble shape – Case $Q > 0$

When the flow emerges from the singularity $O(r = 0)$, the initial surface at $z = z_0$ is inflated
in the radial direction, giving rise to a hemispherical shape connected on its boundary to a
nearly flat surface. Like for $Q < 0$, this peculiar shape may be related to interfaces observed in
two-phase fluid flows like droplet impinging onto a liquid film on a wall. A last interpretation
could be the rising of a gas bubble which emerges through a gas/liquid interface. However, this
case may be more critical than for $Q < 0$ because the region of rapid variation of curvature is
much wider.

In this part, the initial surface is located at $z_0 = 0.01$ and the volume flow rate is equal
to $4\pi Q = +0.4\pi/3$. Again, the surface advection is carried out until time $t = 9$ s. Figure 3.23
illustrates the shape of the surface at the final instant, for vertices advected in an exact way.
<figure><img src="image_57.png" /><figcaption>FIGURE 3.23 – Top and side views of the discrete mesh transported in a spherical velocity field, $Q > 0$</figcaption></figure>

### Exact transport – Ideal case

Like in Sec. 3.4.2 for $Q < 0$, the exact advection of the initial surface meshed by equilateral triangles is first studied, where the additional nodes, stemming from the refinement step, are projected on the analytical surface to ensure that all vertices lie on this surface. Figure 3.24a exhibits that the errors on the curvature evaluation, calculated with the Laplace-

<figure><img src="image_58.png" /><figcaption>FIGURE 3.24 – Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$ and $\mathbf{n}$, $p = 0\%$ with exact transport – $Q > 0$</figcaption></figure>

Beltrami Operator and Integral Formulations methods, stagnate when the mesh size decreases. The Surface Reconstruction methods SR-$\mathcal{N}_1$-$\mathbf{n}_m$-0 and SR-$\mathcal{N}_2$-$\mathbf{n}_m$-0 converge with a slope slightly better than 1.

Applying 20%-perturbations to the initial mesh does not change the behaviors of the errors with the mesh refinement : the levels of the errors are only slightly increased for the curvature and the normal vector approximations (Fig. C.4 in C.1).

### Numerical scheme for the advection

As a reminder, the advection is performed with a second order Runge-Kutta scheme (Heun's method) and the vertex velocity is interpolated by the PERM method of [74].

The accounting of the advection scheme for the initial surface meshed by equilateral triangles increases unsurprisingly the error levels, especially for the curvature variable (Fig. 3.25a versus Fig. 3.24a). On the finest mesh, the saturation of the LBO- and IF-methods is found again in
<figure><img src="image_59.png" /><figcaption>FIGURE 3.25 – Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$ and $\mathbf{n}$, $p = 0\%$ with Lagrangian transport - $Q > 0$</figcaption></figure>

$L_2$-norm and $L_\infty$-norm. Only the SR-$\mathcal{N}_2\mathbf{n}_m$-$a_{00}$ method sticks out of the others and gives the smallest level of errors on the curvature. The order of convergence, beginning from 1 on coarse meshes, becomes smaller than 1 when refining the mesh.

As usual, the error levels for the normal vector calculation are quite low (Fig. 3.25b), compared to those obtained for the curvature. The approximations based on the smallest stencil provide the best results (LBO-Meyer, IF-Trygg, and SR-$\mathcal{N}_1\mathbf{n}_m$-$0$). The two other methods, namely IF-Triang. and SR-$\mathcal{N}_2\mathbf{n}_m$-$a_{00}$, give similar errors, slightly larger than the previous ones.

The addition of 20% random perturbations to the initial discretization of the flat surface increases slightly the error levels (Fig. C.5 in C.1). For the curvature evaluation (Fig. C.5a), the SR-$\mathcal{N}_2\mathbf{n}_m$-$a_{00}$ method seems converging with a first order accuracy, except for the two finest meshes for which the error tends to stagnate. This increase in the relative error is found for all the methods, since it is related to a bad approximation of the surface by the discrete mesh, as discussed here-above for $p = 0\%$, and not to the evaluation method of the curvature itself.

Concerning the approximation of the normal vector (Fig. C.5b), the same comments given for $p = 0\%$ apply : except the finest meshes, for which SR-$\mathcal{N}_2\mathbf{n}_m$-$a_{00}$ is the most accurate, the methods based on the smallest stencils (LBO-Meyer, IF-Trygg. and SR-$\mathcal{N}_1\mathbf{n}_m$-$0$) provide usually the lower errors.

## CPU-cost and efficiency

The same analysis, as the one carried out for $Q < 0$ on the relative accuracy versus the CPU-time (Sec. 3.4.2), is now considered (Fig. 3.26). The data of the relative error are extracted from Fig. 3.25. The SR-$\mathcal{N}_2\mathbf{n}_m$-$0$ method, which was clearly better than the others for estimating the curvature when $Q < 0$, turns out to be also the best choice, even if its efficiency becomes less obvious (Fig. 3.26a). However, this is the only one capable of providing approximations lower than 10% in $L_2$-norm, with a $L_\infty$-norm often very much lower than those given by the other methods.

On the whole, the IF-Trygg. method [123] is the most efficient for the calculation of the normal vector (Fig. 3.26b). It is more than ten times less expensive than SR-$\mathcal{N}_2\mathbf{n}_m$-$a_{00}$ method for coarse meshes, but this gap is closed for the finest surface discretizations. To ensure that the normal vectors are well approximated on all the mesh vertices, the $L_\infty$-norm must be considered. In that case, only the SR-$\mathcal{N}_2\mathbf{n}_m$-$a_{00}$ method allows getting relative errors lower than 1%. The computation times for disturbed initial meshes with $p = 20\%$ provide similar conclusions (C.1, Fig. C.6).
<figure><img src="image_60.png" /><figcaption>FIGURE 3.26 – Relative errors on $κ$ (a) and $n$ (b) according to CPU time in $L_2$ (continuous line) and $L_∞$ (dash line) norms, $p = 0\% - Q > 0$.</figcaption></figure>

### 3.4.4 Summary

The robustness of the Surface Reconstruction methods has already been proved for stationary surfaces, regarding the initial disturbances. The key point to get a good approximation remains the evaluation of the curvature of the surface which requires the most restrictive condition on the mesh size. The accounting of the mesh vertex advection by a numerical scheme deteriorates the relative accuracy of the solution whatever the methods to estimate the surface properties. However, the SR-$N_2-n_m-a_{00}$ method still remains the most competitive. In a more general way, methods using a larger stencil (SR-$N_2-n_m-a_{00}$ and IF-Triang. methods) are the most accurate in the dynamic case with Lagrangian transport.

## 3.5 Conclusion

Three families of numerical methods have been investigated to evaluate the local curvature and the unit normal vector to a surface. They were compared on steady and unsteady interfaces, the later being newly proposed in the present work. To this end, an original analytical solution of a time dependent surface has been developed. Starting from a flat surface and depending on the direction of the flow rate, a long thin ligament or a hemispherical shape bubble joined at its boundary to a flat surface grows with time. These dynamic surface evaluation problems are interesting as they allow to recover geometrical trends found in real two-phase flows.

A first general conclusion, valid for both static and dynamic surfaces, is that there is no major difficulty to get an accurate approximation of the normal vector. The main issue lies in the curvature calculation. Indeed, its analytical expression is based on second order derivatives and then requires regularity. The Laplace-Beltrami-Operator approximations and the Integral Formulation methods are very accurate and efficient for regular meshes. However, as soon as perturbations are introduced, for example on the equilateral triangle mesh for the stationary surface or through the remeshing procedures (dynamic surfaces) that are used during the time integration to track the surface, the curvature evaluation does not converge any more with the decrease of the mesh size. The third approach, the Surface Reconstruction method, is basically more expensive than the previous popular methods. However, it is more robust and provides better results as soon as the surface mesh evolves or is perturbed. Nonetheless, special cares must be taken for SR-$N_1$ methods as it is done in this paper. Their neighbourhood must be extended to the second one when the degree of freedom for reconstructing the surface is not reached. Generally, even if the SR methods are less accurate on coarse meshes, they become more robust with mesh refinement. So, the Surface Reconstruction methods using the first neighbouring
rhood excel in the dynamic case with exact transport, that is the initial and additional vertices brought by the mesh management process are projected back on the exact surface. Whereas the SR-$\mathcal{N}_2$ methods shine in a more segregate case with Lagrangian transport and mesh perturbation. Although expensive, the Surface Reconstruction methods turn out to be the most efficient methods when the surface is advected : for a given relative accuracy, the CPU-time is often the lowest compared to other methods.

This thoroughly study on the surface geometric properties has emphasized the fact that the remeshing process plays a fundamental role to get a correct evaluation of the curvature. An expected improvement could be to make use of the fitted surface to preserve a good representation of the theoretical one. The new or moved vertices, stemming from the refinement or coarsening step of the discrete mesh, could be projected on the local reconstructed surface in order to preserve a high quality approximation of the analytical surface. Moreover, the approaches used in the Surface Reconstruction methods are the most used and also among the simplest. Other techniques exist (especially in the domain of computer aided graphics that could give better results [14, 95, 100]). However, it has to be noted that all these approaches have to be coupled with Navier-Stokes resolution. Hence, the precision brought by these methods must be compared in light of those used in the Navier-Stokes discretisation.
# Bibliographie

## Accuracy and convergence of the curvature and normal vector discretizations for 3D static and dynamic front-tracking interfaces

[10] J. BRACKBILL, D. KOTHE et C. ZEMACH. "A continuum method for modeling surface tension". In: *J. Comput. Phys.* 100.2 (1992), p. 335-354 (cf. p. 31, 32, 73).

[13] R. CANU, S. PUGGELLI, M. ESSADKI, B. DURET, T. MÉNARD, M. MASSOT, J. REVEILLON et F. DEMOULIN. "Where does the droplet size distribution come from?" In: *Int. J. Multiphase Flow* 107 (2018), p. 230-245 (cf. p. 36, 72, 73).

[14] F. CAZALS et M. POUGET. "Estimating differential quantities using polynomial fitting of osculating jets". In: *Comput. Aided Geom. Des.* 22.2 (2005), p. 121-146 (cf. p. 97).

[21] J. DELHAYE. "Jump conditions and entropy sources in two-phase systems. Local instant formulation". In: *Int. J. Multiphase Flow* 1.3 (1974), p. 395-409 (cf. p. 23, 72).

[22] F. DEMOULIN, P. BEAU, G. BLOKKEEL, A. MURA et R. BORGHI. "A new model for turbulent flows with large density fluctuations: Application to liquid atomization". In: *Atomization and Sprays* 17 (2007), p. 315-345 (cf. p. 73).

[23] M. DESBRUN, M. MEYER, P. SCHRODER et A. BARR. "Implicit fairing of irregular meshes using diffusion and curvature flow". In: *SIGGRAPH99*. 1999, p. 317-324 (cf. p. 75).

[26] J. DU, B. FIX, J. GLIMM, X. JIA, X. LI, Y. LI et L. WU. "A simple package for front tracking". In: *J. Comput. Phys.* 213 (2006), p. 613-628 (cf. p. 29, 33, 73, 78, 79, 83, 85, 90, 92, 107, 161, 163).

[31] M. ESSADKI. "Contribution to a unified Eulerian modeling of fuel injection : from dense liquid to polydisperse evaporating spray." Thèse de doct. Universite Paris-Saclay, 2018 (cf. p. 75, 76, 81).

[34] R. FEDKIW, T. ASLAM, B. MERRIMAN et S. OSHER. "A non-oscillatory Eulerian approach to interfaces in multi-material flows (the Ghost Fluid method)". In: *J. Comput. Phys.* 152 (1999). Sous la dir. de J. of COMPUTATIONAL PHYSICS, p. 457-492 (cf. p. 30, 73).

[45] R. GOLDMAN. "Curvature formulas for implicit curves and surfaces". In: *Comput. Aided Geom. Des.* 22 (2005), p. 632-658 (cf. p. 34, 78, 79).

[67] X. LI, G. XU et Y. J. ZHANG. "Localized discrete Laplace-Beltrami operator over triangular mesh". In: *Computer Aided Geometric Design*. T. 39. 2015, p. 67-82 (cf. p. 73, 80, 88).

[73] N. MAX. "Weights for computing vertex normals from facets normals." In: *Journal of Graphics Tools* 4.2 (1999), p. 1-6 (cf. p. 79).

[74] R. MCDERMOTT et S. POPE. "The parabolic edge reconstruction method (PERM) for Lagrangian particle advection". In: *J. Comput. Phys.* 227.11 (2008), p. 5447-5491 (cf. p. 49, 89, 92, 94).
[76] M. MEYER, M. DESBRUN, P. SCHRODER et A. H. BARR. "Discrete Differential-Geometry Operators for Triangulated 2-Manifolds". In: *Visualization and Mathematics III*. Sous la dir. de H.-C. HEGE et K. POLTHIER. Springer, 2003, p. 35-57 (cf. p. 33, 36, 75, 76, 81, 85, 90, 92, 107, 152, 154).

[77] E. MOHAMED, D. FLORENCE, de CHAISEMARTIN STÉPHANE, L. ADAM, M. THIBAULT et M. MARC. "Statistical modeling of the gas-liquid interface using geometrical variables : Toward a unified description of the disperse and separated phase flows". In : *Int. J. Multiphase Flow* 120 (2019), p. 103084 (cf. p. 36, 73).

[91] K. POLTHIER, J. HASS, D. HOFFMAN, A. JAFFE, H. ROSENBERG, R. SCHOEN et M. WOLF. "Computational aspects of discrete minimal surfaces." In : *Proc. of the Clay Summer school on Global Theory of Minimal Surfaces*. 2002 (cf. p. 75).

[94] S. POPINET et S. ZALESKI. "A front-tracking algorithm for accurate representation of surface tension." In : *Int. J. Numer. Meth. Fluids* 30 (1999), p. 775-793 (cf. p. 76, 78).

[95] H. POTTMANN, J. WALLNER, Y.-L. YANG, Y.-K. LAI et S.-M. HU. "Principal curvatures from the integral invariant viewpoint". In : *Comput. Aided Geom. Des.* 24.8 (2007). Discrete Differential Geometry, p. 428-442 (cf. p. 97).

[99] I. ROGHAIR, V. S. ANNALAND et J. KUIPERS. "An improved Front-Tracking technique for the simulation of mass transfer in dense bubbly flows". In : *Chem. Eng. Sci.* 152 (oct. 2016), p. 351-369 (cf. p. 50, 55, 59, 89).

[100] S. RUSINKIEWICZ. "Estimating curvatures and their derivatives on triangle meshes". In : *Proceedings. 2nd International Symposium on 3D Data Processing, Visualization and Transmission, 2004. 3DPVT 2004.* 2004, p. 486-493 (cf. p. 97).

[108] S. SHIN et D. JURIC. "Modeling Three-Dimensional Multiphase Flow Using a Level Contour Reconstruction Method for Front Tracking without Connectivity". In : *J. Comput. Phys.* 180.2 (2002), p. 427-470 (cf. p. 29, 73).

[116] G. TAUBIN. "A signal processing approach to fair surface design." In : *SIGGRAPH '95*. 1995, p. 351-358 (cf. p. 75).

[121] P. TRONTIN, S. VINCENT, J. ESTIVALEZES et J. CALTAGIRONE. "A subgrid computation of the curvature by a particle/level-set method. Application to a front-tracking/ghost-fluid method for incompressible flows". In : *J. Comput. Phys.* 231.20 (2012), p. 6990-7010 (cf. p. 78).

[122] G. TRYGGVASON, B. BUNNER, A. ESMAEELI, D. JURIC, N. AL-RAWAHI, W. TAUBER, J. HAN, S. NAS et Y.-J. JAN. "A Front-Tracking Method for the Computations of Multiphase Flow". In : *J. Comput. Phys.* 169.2 (2001), p. 708-759 (cf. p. 29-31, 76, 77).

[123] G. TRYGGVASON, R. SCARDOVELLI et S. ZALESKI. *Direct Numerical Simulations of Gas-Liquid Multiphase Flows*. 1re éd. Cambridge University Press, 2011 (cf. p. 31, 77, 82, 85, 90, 93, 95, 107).

[125] G. VAUDOR, T. MÉNARD, W. ANISZEWSKI, M. DORING et A. BERLEMONT. "A consistent mass and momentum flux computation method for two phase flows. Application to atomization process". In : *Computers and Fluids* 152 (2017), p. 204-216 (cf. p. 73).

[128] M. WARDETZKY. "Convergence of the Cotangent Formula : An Overview". In : *Discrete Differential Geometry* 38 (2008), p. 275-286 (cf. p. 88).

[132] G. XU. "Discrete Laplace-Beltrami operators and their convergence". In : *Comput. Aided Geom. Des.* 21.8 (2004). Geometric Modeling and Processing 2004, p. 767-784 (cf. p. 88).

[136] D.-h. YU. *Natural Boundary Integral Method and Its Applications*. T. 539. Mathematics and Its Applications. Springer Netherlands, 2002 (cf. p. 73).
[138] A. Z. ZINCHENKO, M. A. ROTHER et R. H. DAVIS. "A novel boundary-integral algorithm for viscous interaction of deformable drops." In : *Phys. Fluids* 9 (1997) (cf. p. 78, 79, 83).

[139] D. ZUZIO, J.-L. ESTIVALEZES, P. VILLEDIEU et G. BLANCHARD. "Numerical simulation of primary and secondary atomization". In : *Comptes Rendus Mécanique* 341 (2013), p. 15-25 (cf. p. 73).
# Chapitre 4

## Identification d'objets grâce à leurs propriétés surfaciques

### Sommaire

<table><tr><td>4.1</td><td>Précision des méth. de calc. de κ sur un nuage de gouttes de forme sphé.</td><td>103</td></tr><tr><td>4.1.1</td><td>Mise en œuvre de l'étude</td><td>104</td></tr><tr><td>4.1.2</td><td>Grandeurs enregistrées</td><td>106</td></tr><tr><td>4.1.3</td><td>Résultats et analyse</td><td>106</td></tr><tr><td>4.2</td><td>Echelles intermédiaires - Caractérisation des ligaments</td><td>112</td></tr><tr><td>4.2.1</td><td>Distribution surfacique de courbure</td><td>112</td></tr><tr><td>4.2.2</td><td>Méthodologie de l'étude numérique</td><td>113</td></tr><tr><td>4.2.3</td><td>Résultats et analyses</td><td>119</td></tr><tr><td>4.3</td><td>Ellipsoïdes déformées par des fonctions sinueuses et variqueuses</td><td>126</td></tr><tr><td>4.3.1</td><td>Déformation sinueuse</td><td>130</td></tr><tr><td>4.3.2</td><td>Déformation variqueuse</td><td>134</td></tr></table>

Nous avons montré dans le chapitre précédent la précision et la robustesse de différentes méthodes de calcul des propriétés géométriques de surface. Ainsi, parmi les trois familles étudiées (LBO, IF, SR), les méthodes de Reconstruction de Surface (SR) s'avèrent les plus fiables pour évaluer les caractéristiques géométriques de surfaces (courbure, normale), et par conséquent pour avoir une bonne estimation des forces de tension de surface. Toutefois, leur rôle ne se limite pas qu'à cela. Ainsi, dans ce chapitre, nous présentons un travail préliminaire dont l'objectif est l'utilisation des méthodes de calcul des propriétés géométriques de surface dans l'identification de formes. La finalité serait d'avoir un système d'apprentissage utilisant une base de données de différentes structures (sphères, ligaments) plus ou moins déformées, qui permettrait d'identifier l'ensemble de formes apparaissant dans un jet atomisé par exemple. Pour ce faire, nous introduisons une quantité basée sur la répartition des courbures de la structure, que nous évaluons sur des formes simples.

### 4.1 Précision des méthodes de calcul de courbure sur un nuage de gouttes de forme sphérique

En prenant l'exemple d'un jet atomisé, plusieurs échelles sont mises en jeu : les grandes échelles représentées par le jet principal, les échelles intermédiaires qui comprennent des ligaments et des gouttes plus ou moins grosses de forme sphérique ou non. Quant aux petites échelles, elles sont constituées de gouttes de forme sphérique ou quasi-sphérique. L'attention
dans cette partie est portée sur cette dernière classe d'échelles où les forces capillaires dominent, ce qui engendre des structures de forme sphérique. Le but ici est de déterminer la précision des méthodes de calcul de courbure et de vecteur normal pour identifier les gouttes à travers leur courbure moyenne sur l'ensemble des marqueurs constituant le maillage. Ainsi, dans un premier temps, est décrite la mise en œuvre numérique de cette étude, dans laquelle deux paramètres interviennent : la taille et la régularité du maillage de l'interface. Par la suite, les résultats des simulations sont exposés et commentés.

### 4.1.1 Mise en œuvre de l'étude

<figure><img src="image_61.png" /><figcaption>FIGURE 4.1 – Gouttes sphériques avec différentes résolutions de courbure $d_r = 1/2r$, $r = 1, 10$.</figcaption></figure>

Soit une sphère de rayon $r$. Son interface est initialement discrétisée par des triangles quasi-équilatéraux. Quel que soit le rayon $r$ de la sphère, la taille moyenne $d$ des arêtes des triangles est définie par $d \approx 1/2$. Le paramètre caractéristique le plus important n'est pas la taille $d$ des arêtes ou le rayon $r$ de la sphère mais la résolution $d_r$ de la courbure qui est le rapport entre $d$ et $r$ ($d_r = d/r$), dans notre cas $d_r = 1/2r$. La figure 4.1 présente différentes sphères pour $r \in [1, 10]$.

Afin d'étudier la sensibilité du calcul numérique des courbures, le maillage initialement régulier est perturbé. Ce processus est détaillé dans la section suivante, où on parlera indistinctement de gouttes, bulles ou sphères.

<figure><img src="image_62.png" /><figcaption>FIGURE 4.2 – Perturbation maximale du maillage de l'interface de 5% et 40%.</figcaption></figure>

#### Perturbation du maillage

Soit $m_0 \in \mathcal{M}$ un marqueur de la discrétisation $\mathcal{D}$ de la sphère de résolution de courbure $d_r$ et de surface $\partial\Omega$. Ses coordonnées sont désignées par $\bar{x}_{m_0}$ dans la base canonique ($e_1, e_2, e_3$)
<figure><img src="image_63.png" /><figcaption>FIGURE 4.3 – Cône dans lequel est déplacé le marqueur $m_0$.</figcaption></figure>

de $\mathbb{R}^3$, avec $e_\alpha^\beta = \delta_{\alpha\beta}$. La perturbation est effectuée en se plaçant dans une nouvelle base orthonormée directe $(\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3)$, avec comme origine du repère le centre $O$ de la sphère et tel que $\mathbf{b}_3 = \bar{x}_{m_0}/\|\bar{x}_{m_0}\|$ (Fig. 4.3). Le marqueur $m_0$ est potentiellement bougé dans une calotte sphérique de contour $\Gamma$, définie par l'intersection entre la surface $\partial\Omega$ et un cône de sommet $O$. Il est construit tel que pour tout point $A \in \Gamma$ et de coordonnées $\bar{x}_A$ :

$$ \|\bar{x}_{m_0} - \bar{x}_A\|/r = p \cdot d_r $$

avec $p$ l'amplitude maximale de perturbation fixée et $d_r$ la résolution de la courbure de la sphère. Après perturbation, la nouvelle position aléatoire de $m_0$ (4.1) dans la base $(\mathbf{b}_1, \mathbf{b}_2, \mathbf{b}_3)$ est repérée par $(\theta, \varphi)$.

$$ \bar{x}_{m_0}^{\text{new}} = \begin{cases} x = r \sin(\theta) \cos(\varphi) \\ y = r \sin(\theta) \sin(\varphi) \\ z = r \cos(\theta) \end{cases} \qquad (4.1) $$

Les deux angles sont fonction d'un couple de réels $(\alpha_m, \beta_m) \in [0, 1]^2$ issus d'une distribution uniforme standard et du demi-angle $\theta_{max}$ définissant l'ouverture du cône.

$$ \begin{aligned} \theta &= \alpha_m \theta_{max} \\ \varphi &= \beta_m 2\pi \end{aligned} \qquad (4.2) $$

Un dernier changement de base est effectué pour revenir dans la base canonique de $\mathbb{R}^3$.

**Calcul de $\theta_{max}$** : il est évalué en se servant du triangle formé par les points $O, A$ et de la position initiale du marqueur $m_0$.

$$ \begin{aligned} \cos(\theta_{max}/2) &= \frac{[r^2 - p^2 \cdot d_r^2/4]^{1/2}}{r} \\ &= \left[1 - \frac{p^2 \cdot d_r^2}{4}\right]^{1/2} \\ \sin(\theta_{max}/2) &= \frac{p \cdot d_r}{2} \end{aligned} $$
La formule trigonométrique $\cos(a + b) = \cos(a) \cos(b) - \sin(a) \sin(b)$ permet d'obtenir :

$$ 
\begin{aligned}
\cos(\theta_{max}) &= \left[ 1 - \frac{p^2 \cdot d_r^2}{2} \right] \\
\theta_{max} &= \arccos \left[ 1 - \frac{p^2 \cdot d_r^2}{2} \right]
\end{aligned}
\qquad (4.3) $$

### 4.1.2 Grandeurs enregistrées

Soit une sphère de résolution de courbure $d_r$ fixée. Nous rappelons que $d_r = d/r$, avec $d$ la taille des arêtes de triangles du maillage de l'interface avant perturbation et $r$ le rayon de la sphère. Pour une résolution de courbure donnée, $10^5$ sphères sont construites, chacune avec un maillage aléatoire issu d'une distribution uniforme standard et une perturbation maximale fixée ($p \in [0, 0.4]$), comme décrit dans la section 4.1.1. Nous notons $M_r^{(i)}$ l'ensemble des marqueurs $m$ de la i-ème sphère de résolution de courbure $d_r$, $i \in [1, 10^5]$. Ensuite, pour chaque sphère de rang $i$ et de résolution $d_r$ est définie la quantité $\langle \kappa \rangle_r^{(i)}$, moyenne arithmétique des courbures en chaque marqueur de l'interface de la sphère (4.4) $i$.

$$ \langle \kappa \rangle_r^{(i)} = \sum_{m \in M_r^{(i)}} \kappa_m / \text{card}(M_r^{(i)}) \qquad (4.4) $$

L'écart relatif à la courbure exacte $d_r/d = 1/r$ est alors :

$$ Err(\langle \kappa \rangle_r^{(i)}) = (\langle \kappa \rangle_r^{(i)} d/d_r - 1) \qquad (4.5) $$

L'objectif est de voir la répartition (ou PDF) de l'écart entre cette moyenne arithmétique des courbures sur une sphère et sa courbure exacte sur une population de $10^5$ sphères. L'impact sur cette distribution de l'amplitude maximale $p$ des perturbations du maillage, de la résolution $d_r$ de la courbure et du choix des méthodes de calculs des courbures en chaque marqueur est également discuté.

### 4.1.3 Résultats et analyse

Les graphes présentés dans cette section donnent les distributions des erreurs relatives (4.5) sur les $10^5$ sphères pour les différents niveaux $d_r$ de raffinement du maillage et de l'amplitude $p$ de la perturbation.

#### Influence de la résolution $d_r$ de courbure

La résolution de la courbure est définie par le paramètre $d_r$ qui représente le rapport de la taille des arêtes d'un triangle sur le rayon de la sphère. Les simulations ont été menées pour $d_r = 1/2r$, $r = 1, 10$ et $p \in \{5\%, 20\%, 40\%\}$. Les résultats pour $r \in \{2, 4, 8\}$ sont présentés par la suite, les autres suivant la même tendance.

Dans un cadre idéal, à savoir un maillage de l'interface constitué de triangles quasi-équilatéraux dans le cas $p = 5\%$, les méthodes les plus adaptées pour identifier des sphères sont LBO-Meyer, IF-Triangle, et IF-Trygg, avec une erreur relative inférieure à 1%, quelle que soit la résolution $d_r$ (Fig. 4.4). Ensuite viennent les méthodes SR-$\mathcal{N}_1$-$\mathbf{n}_m$-0 et SR-$\mathcal{N}_2$-$\mathbf{n}_m$-0, qui pour de faibles résolutions ($d_r = 1/4$) ont une erreur comprise entre 2% et 9% (Fig. 4.4a). Cette erreur diminue avec le raffinement du maillage de l'interface, pour descendre en dessous de 0.6% quand
$d_r = 1/16$ (Fig. 4.4c).

L'introduction de perturbations ($p = 20\%$) baisse la précision des méthodes LBO-Meyer, IF-Triangle et IF-Trygg. Pour de faibles résolutions ($d_r = 1/4$), elles restent les mieux appropriées pour identifier une goutte sphérique (Fig. 4.5a). Par contre, elles ne convergent pas avec les maillages, voire divergent pour certains. Ainsi, pour les différents maillages $d_r \in \{1/4, 1/8, 1/16\}$, les valeurs moyennes de la distribution pour LBO-Meyer restent autour de 0.1%. Elles sont respectivement de 0.5%, 0.8% et 0.9% pour IF-Triangle, et de 2%, 2.3% et 2.4% pour IF-Trygg (Fig. 4.5). Les méthodes $\text{SR-}\mathcal{N}_1\text{-n}_m\text{-0}$ et $\text{SR-}\mathcal{N}_2\text{-n}_m\text{-0}$, quant à elles, convergent avec le maillage, avec un ordre de convergence de 2 : pour $\text{SR-}\mathcal{N}_1\text{-n}_m\text{-0}$, les valeurs moyennes de distribution de l'erreur relative sont d'environ 2%, 0.5% et 0.12%; elles sont de 8%, 2% et 0.5% pour $\text{SR-}\mathcal{N}_2\text{-n}_m\text{-0}$. Ces écarts entre $\text{SR-}\mathcal{N}_1$ et $\text{SR-}\mathcal{N}_2$ se justifient par le fait que $\text{SR-}\mathcal{N}_2$ possède un stencil plus étendu. Cela a pour avantage de régulariser d'éventuels bruits sur l'interface, mais comporte la contrepartie d'être moins local et de ce fait être moins précis.

L'introduction de perturbations nous a montré que la hiérarchie que nous avions pour de faibles perturbations (LBO et IF prédominants à $p = 5\%$, Fig. 4.4) a été chamboulée pour des perturbations plus importantes (SR et LBO prédominants à $p = 20\%$, Fig. 4.5). D'où l'importance d'étudier l'influence des perturbations sur ces tests.

## Influence de la perturbation $p$

Cette partie est consacrée à l'impact de la perturbation $p$ du maillage de l'interface sur la distribution de l'écart relatif à la courbure exacte. Trois amplitudes de perturbations sont effectuées : $p \in \{0.05, 0.2, 0.4\}$. La résolution initiale de la courbure est $d_r = 1/20$. Les résultats sont présentés sur la figure 4.6.

Pour de faibles perturbations $p = 5\%$, toutes les méthodes sont équivalentes, et ont une erreur relative $Err(<\kappa>) \le 0.4\%$ (Fig. 4.6a). En effet, la forme des triangles équilatéraux initiaux est peu perturbée, ce qui est avantageux pour les méthodes LBO-Meyer [76], IF-Triangle et IF-Trygg [123]. Ces dernières dépendent fortement de la qualité du maillage. Quant aux méthodes de Reconstruction de Surface $\text{SR-}\mathcal{N}_1\text{-n}_m\text{-0}$ [26] et $\text{SR-}\mathcal{N}_2\text{-n}_m\text{-0}$, la densité de marqueurs pour représenter la surface est l'élément prépondérant dans la qualité de l'évaluation des propriétés géométriques de la surface. Ainsi, la résolution $d_r = 1/20$, à savoir la taille des arêtes des triangles vaut $1/20$ du rayon $r$ de la sphère, est suffisante pour avoir une erreur inférieure à 0.4%.

L'augmentation de la perturbation $p$ a une double influence sur les distributions.

- Les méthodes LBO-Meyer, IF-Triangle et IF-Trygg voient leur erreur moyenne de la distribution s'écarter de 0%. En d'autres termes, elles deviennent moins précises. Ainsi, au maximum de perturbation $p = 40\%$, les moyennes de distributions de LBO-Meyer, IF-Triangle et IF-Trygg deviennent respectivement 1%, 2% et 5.5% (Fig. 4.6c). En effet, les perturbations brisent la forme équilatérale des triangles, à laquelle sont sensibles ces méthodes. En ce qui concerne les méthodes de Reconstruction de Surface, elles sont peu, voire pas impactées.

- En plus de l'augmentation de la valeur moyenne de la distribution, les perturbations étalent la distribution pour les méthodes LBO-Meyer, IF-Triangle et IF-Trygg (Fig. 4.6b et 4.6c). Comme précédemment, les méthodes $\text{SR-}\mathcal{N}_1\text{-n}_m\text{-0}$ et $\text{SR-}\mathcal{N}_2\text{-n}_m\text{-0}$ ne sont pas impactées.

Cela vient confirmer les conclusions du chapitre 3 sur la sensibilité des méthodes LBO et IF à la qualité des triangles du maillage de l'interface. Ces dernières dépendent fortement de la forme des triangles reliés au marqueur considéré. Elles donnent de bons résultats lorsque les triangles sont presque équilatéraux, mais ont des résultats dégradés dès que l'on s'écarte de cette configuration idéale.
Conclusion

L'étude de la précision des méthodes quant à l'identification d'une population de $10^5$ gouttes sphériques montre que pour des résolutions de courbure $d_r$ grossières, la méthode LBO-Meyer est celle qui donne de meilleurs résultats. Par contre, si on veut s'affranchir de l'influence des perturbations du maillage, qui est incontournable lors de la simulation d'écoulements diphasiques, et contrôler les erreurs, le choix est porté sur les méthodes de Reconstruction de Surface, dans ce cas-ci $SR-\mathcal{N}_1-n_m-0$ et $SR-\mathcal{N}_2-n_m-0$.
<figure><img src="image_64.png" /><figcaption>(a) Perturbation $p = 5\%$ et $d_r = 1/4$</figcaption></figure>

<figure><img src="image_65.png" /><figcaption>(b) Perturbation $p = 5\%$ et $d_r = 1/8$</figcaption></figure>

<figure><img src="image_66.png" /><figcaption>(c) Perturbation $p = 5\%$ et $d_r = 1/16$</figcaption></figure>

FIGURE 4.4 – Influence de la résolution de la courbure : PDF de l'écart relatif $Err(<\kappa>)$, pour différentes résolutions $d_r$ et une perturbation $p = 5\%$.
4.1. Précision des méth. de calc. de κ sur un nuage de gouttes de forme sphé.

<figure><img src="image_67.png" /><figcaption>FIGURE 4.5 – Influence de la résolution de la courbure : PDF de l'écart relatif Err(<κ>)
pour différentes résolutions d<sub>r</sub> et une perturbation p = 20%.</figcaption></figure>
<figure><img src="image_68.png" /><figcaption>(a) Perturbation $p = 5\%$ et $d_r = 1/20$</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>(b) Perturbation $p = 20\%$ et $d_r = 1/20$</figcaption></figure>

<figure><img src="image_70.png" /><figcaption>(c) Perturbation $p = 40\%$ et $d_r = 1/20$</figcaption></figure>

FIGURE 4.6 – Influence de la perturbation $p$ : PDF de l'écart relatif $Err(<\kappa>)$, pour différentes perturbations $p$ et une résolution de courbure $d_r = 1/20$.
## 4.2 Echelles intermédiaires - Caractérisation des ligaments

Précédemment, nous nous sommes intéressés aux petites échelles constituées de gouttes sphériques. Cette partie, centrée sur les ligaments, est une approche préliminaire dans l'identification à travers la caractérisation de ces structures. Nous présenterons, dans un premier temps, la distribution surfacique de la courbure moyenne, élément que nous utiliserons pour caractériser les formes rencontrées. Cette distribution est par la suite testée sur des formes simples, à savoir des sphères déformées progressivement en ellipsoïdes. Un degré de complexité est apporté en introduisant des déformations sinueuses et variqueuses aux structures ellipsoïdales formées.

### 4.2.1 Distribution surfacique de courbure

Soit un domaine $\Omega$ de surface $\partial\Omega$. La distribution surfacique de courbure représente la fraction de la surface $\partial\Omega$ pour laquelle la courbure moyenne est la même. Avant de présenter cette grandeur, nous introduisons une quantité beaucoup plus générale notée $DF_G(\eta)$ ($DF$ pour density function). Elle représente la distribution associée à la fonction $G$ et la quantité $\eta$. Elle est définie par la relation (4.6).

$$ \forall \eta \in \mathbb{R}, \quad DF_G(\eta) = \left[ \int_{\mathbb{R}} G(\hat{\eta}) \delta(\hat{\eta} - \eta) d\hat{\eta} \right] / \int_{\mathbb{R}} G(\hat{\eta}) d\hat{\eta} \quad (4.6) $$

avec $\eta \in \{\kappa, \kappa_G, \kappa_{min}, \kappa_{max}\}$, $G$ une fonction travaillant sur la quantité $\eta$. Nous choisissons de travailler par la suite avec $\eta = \kappa$. Aussi, selon les valeurs de la fonction $G$, nous pouvons étudier différents cas de figure :

- $G(\kappa) = 1$ : En utilisant la constante 1, la distribution $DF_1(\kappa)$ donne l'occurrence ou le nombre de sommets possédant cette courbure $\kappa$ sur le domaine étudié.

- $G(\kappa) = S(\kappa)$ : Nous retrouvons ici la distribution surfacique de courbure. Avec $S(\kappa)$ la surface de courbure $\kappa$, $DF_S(\kappa)$ représente la fraction de surface de courbure moyenne $\kappa$.

Dans notre étude, la surface $\partial\Omega$ possède une discrétisation $\mathcal{D}$, avec $\mathcal{M}$ l'ensemble des marqueurs $m$. La relation (4.6) a besoin d'être adaptée à cette discrétisation. Les différentes grandeurs continues et leurs contreparties discrètes sont rassemblées dans le tableau 4.1.

<table><caption>TABLE 4.1 – Discrétisation des termes de la distribution $DF_G(\kappa)$.</caption><thead><tr><th>Titre</th><th>Continue</th><th>Discret</th></tr></thead><tbody><tr><td>Distribution</td><td>DF<sub>G</sub>(κ)</td><td>df<sub>G</sub>(κ)</td></tr><tr><td>Courbure</td><td>κ</td><td>κ<sub>m</sub>, m ∈ M</td></tr><tr><td>Fonction de κ</td><td>G(κ)</td><td>G<sub>m</sub>, m ∈ M</td></tr><tr><td rowspan="3">Domaine de définition</td><td rowspan="3">R, dκ</td><td>I<sub>κ</sub> = {κ<sub>m</sub>, m ∈ M}</td></tr><tr><td>N<sub>c</sub>, classes (ou subdivisions) de I<sub>κ</sub></td></tr><tr><td>Δκ = [max(I<sub>κ</sub>) - min(I<sub>κ</sub>)] / N<sub>c</sub></td></tr><tr><td>Dirac</td><td>δ(κ - κ̂)</td><td>Π[(κ<sub>m</sub> - κ̂)/Δκ]/Δκ</td></tr></tbody></table>

Ainsi, les grandeurs continues $(\kappa, G(\kappa))$ sur la surface $\partial\Omega$ deviennent des valeurs discrètes $(\kappa_m, G_m)$ sur les marqueurs $m \in \mathcal{M}$. L'intervalle $I_\kappa$ des valeurs possibles de $\kappa$ est subdivisé
en $N_c$ classes de courbure $\hat{\kappa}$ de taille $\Delta\kappa$ constante. Les classes de courbures $\hat{\kappa}$ constituent donc des sous-ensembles de $I_\kappa$ ($\hat{\kappa} \subset I_\kappa$). En dehors de cet intervalle, la fonction $G(\kappa)$ est nulle. La classe de courbure $\hat{\kappa}$ représente l'ensemble des courbures comprises dans l'intervalle $[\hat{\kappa} - \Delta\kappa/2, \hat{\kappa} + \Delta\kappa/2]$. Par la suite, $\hat{\kappa}$ désignera indistinctement cet intervalle et le milieu de cet intervalle. La distribution de Dirac $\delta$ est approximée par le rapport de la fonction porte $\Pi$ et de la taille $\Delta\kappa$ des classes de courbures. La fonction porte est définie par la relation (4.7).

$$ \forall t \in \mathbb{R}, \quad \Pi(t) = \begin{cases} 1 & \text{pour } t \in [-0.5, 0.5] \\ 0 & \text{sinon} \end{cases} \qquad (4.7) $$

Dans le cas de la distribution surfacique de courbure où $G(\kappa) = S(\kappa)$, nous choisissons de construire la surface discrète $S_m$ autour du marqueur $m$ comme étant l'aire composite $A_{mixed}$ de la méthode LBO-Meyer, définie dans la section 3.2.2, page 75.

$$ S_m = A_{mixed} \qquad (4.8) $$

Cette surface est construite en se servant du marqueur $m$, du milieu des arêtes reliées à $m$ et d'un dernier point dépendant de la forme des triangles. Ainsi, sur un triangle aigu, il correspond au centre du cercle circonscrit (Fig. 3.2a, page 76). Dans le cas de triangles obtus, c'est le milieu de l'arête opposée au sommet d'angle obtus qui est pris en compte (Fig. 3.2b, page 76). Pour tout $m \in \mathcal{M}$, les surfaces $S_m$ sont disjointes et leur somme est égale à la surface totale discrète $S_{tot}$.

$$ S_{tot} = \sum_{m \in \mathcal{M}} S_m \qquad (4.9) $$

L'ensemble de ces notations nous permet de poser la discrétisation de la relation (4.6) selon le choix de $G(\kappa)$.

$$ df_1(\hat{\kappa}) = \frac{\sum_{m \in \mathcal{M}} \Pi [(\kappa_m - \hat{\kappa})/\Delta\kappa]}{\sum_{m \in \mathcal{M}} \Delta\kappa} \qquad (4.10) $$

$$ df_S(\hat{\kappa}) = \frac{\sum_{m \in \mathcal{M}} S_m \Pi [(\kappa_m - \hat{\kappa})/\Delta\kappa]}{\sum_{m \in \mathcal{M}} S_m \Delta\kappa} \qquad (4.11) $$

Ainsi, pendant que la distribution $df_1(\hat{\kappa})$ décompte la proportion du nombre de marqueurs ayant une courbure $k_m \in \hat{\kappa}$, la distribution $df_S(\hat{\kappa})$ relève la portion de surface dont la courbure est comprise dans la classe de courbure $\hat{\kappa}$.

### 4.2.2 Méthodologie de l'étude numérique

La distribution $df_G$ définie précédemment (4.10-4.11) est utilisée pour caractériser des ligaments. Ces derniers pouvant être de formes diverses et variées, la densité $df_G$ est d'abord testée sur des formes simples, à savoir une famille d'ellipsoïdes.
<table><caption>TABLE 4.2 – Évolution du nombre nbsom de sommets, nbele d'éléments, de la surface S<sub>tot</sub> et des courbures κ<sub>min</sub> et κ<sub>max</sub> pour les ellipsoïdes d'étirement α ∈ {20%, 80%, 140%} en fonction de la résolution de courbure d<sub>r0</sub> de la sphère initiale servant à les engendrer.</caption><thead><tr><th></th><th colspan="3">α = 20%</th><th colspan="3">α = 80%</th><th colspan="3">α = 140%</th></tr></thead><tbody><tr><th>κ<sub>min</sub></th><td colspan="3">0.8647</td><td colspan="3">0.7858</td><td colspan="3">0.8306</td></tr><tr><th>κ<sub>max</sub></th><td colspan="3">1.44</td><td colspan="3">3.24</td><td colspan="3">5.76</td></tr><tr><th>1/d<sub>r0</sub></th><th>nbsom</th><th>nbele</th><th>S<sub>tot</sub></th><th>nbsom</th><th>nbele</th><th>S<sub>tot</sub></th><th>nbsom</th><th>nbele</th><th>S<sub>tot</sub></th></tr><tr><th>4</th><td>162</td><td>320</td><td>12.47</td><td>408</td><td>812</td><td>13.86</td><td>2110</td><td>4216</td><td>18.81</td></tr><tr><th>8</th><td>642</td><td>1280</td><td>12.66</td><td>1025</td><td>2046</td><td>14.01</td><td>3251</td><td>6498</td><td>15.72</td></tr><tr><th>16</th><td>2562</td><td>5120</td><td>12.71</td><td>3505</td><td>7006</td><td>14.06</td><td>6750</td><td>13496</td><td>15.73</td></tr><tr><th>32</th><td>10242</td><td>20480</td><td>12.72</td><td>13875</td><td>27746</td><td>14.08</td><td>21669</td><td>43334</td><td>15.74</td></tr><tr><th>64</th><td>40962</td><td>81920</td><td>12.73</td><td>55547</td><td>111090</td><td>14.08</td><td>86597</td><td>173190</td><td>15.74</td></tr><tr><th>128</th><td>163842</td><td>327680</td><td>12.73</td><td>223832</td><td>447660</td><td>14.08</td><td>340374</td><td>680742</td><td>15.74</td></tr><tr><th>256</th><td>655362</td><td>1310720</td><td>12.73</td><td>901619</td><td>1803234</td><td>14.08</td><td>1363294</td><td>2726584</td><td>15.74</td></tr><tr><th>512</th><td>2621442</td><td>5242880</td><td>12.73</td><td>3610158</td><td>7220312</td><td>14.08</td><td>5441908</td><td>10883812</td><td>15.74</td></tr></tbody></table>

### Construction de la famille d'ellipsoïdes

La famille d'ellipsoïdes est formée en partant d'une sphère de rayon $r_0$. La sphère est étirée dans la direction $x$ et contractée dans les deux autres directions $y$ et $z$. Le volume est conservé entre la sphère initiale et les ellipsoïdes résultantes dans l'optique de pouvoir les comparer. Ainsi, considérons $m \in \mathcal{M}$, un marqueur de la discrétisation $\mathcal{D}$ de la sphère de rayon $r_0$, et de coordonnées $\bar{x}_m$. La résolution de courbure est notée $d_{r0}$. Nous rappelons que $d_{r0} = d/r_0$, avec $d$ la taille des arêtes des triangles discrétisant la surface. Des perturbations d'amplitude maximale $p \in [0, 0.4]$ sont appliquées sur le maillage de la sphère initiale, comme décrit dans la section 4.1.1. Un étirement du maillage d'amplitude $\alpha \in [0.05, 1.4]$ est opéré dans la direction $x$, de telle sorte que les coordonnées $\bar{x}_m$ de $m$ vérifient la relation (4.12) :

$$ \bar{x}_m/r_0 = \begin{cases} x = a \sin(\theta) \cos(\varphi)/r_0 \\ y = b \sin(\theta) \sin(\varphi)/r_0 \\ z = b \cos(\theta)/r_0 \end{cases} \qquad (4.12) $$

avec $a$ et $b$ respectivement le demi-grand axe et le demi-petit axe de l'ellipsoïde, et :

$$ a/r_0 = (1 + \alpha) \quad , \quad b/r_0 = 1/(1 + \alpha)^{1/2} \qquad (4.13) $$

Le jeu de paramètres définissant les ellipsoïdes se résume donc au triplet $\{d_{r0}, p, \alpha\}$, avec $d_{r0}$ pour la résolution de la courbure et le calcul du volume, $p$ la perturbation du maillage et $\alpha$ l'étirement. Les différentes ellipsoïdes de la famille sont construites à partir de la précédente (Fig. 4.7). Ainsi, l'ellipsoïde d'étirement $\alpha_n$ est générée en se servant de celle d'étirement précédent $\alpha_{n-1}$ ($\alpha_{n-1} < \alpha_n$). A chaque étape, les procédures de gestion du maillage sont appliquées (enrichissement, désenrichissement, lissage, échange d'arêtes), comme présentées dans les sections 2.3-2.5, pp 49-57.

### Choix du nombre $N_c$ de classes de courbures

Comme évoqué dans le tableau 4.1, le passage de la version continue à la formulation discrète de la distribution $df_G$ nous amène à subdiviser l'intervalle $I_k$ des courbures des ellipsoïdes en $N_c$ classes de courbures de même largeur $\Delta\kappa$. Afin de déterminer le nombre de classes pour lequel on a une convergence de la distribution $df_G$, nous nous intéressons à l'évolution de la distribution
<figure><img src="image_71.png" /><figcaption>FIGURE 4.7 – Formation de la famille d'ellipsoïdes à partir d'une sphère de résolution $d_{r0} = 1/16$.</figcaption></figure>

$df_G$ en fonction du nombre de classes $N_c$. Cette convergence est étudiée pour la distribution surfacique ($df_S$) et la distribution de décompte ($df_1$), et l'évolution de $df_G$ pour chaque classe $\hat{\kappa}$ de courbure est regardée.

Ainsi, soit un ellipsoïde de paramètres $\{d_{r0}, p, \alpha\}$ donnés. Un nombre initial de classes $N_c^0$ de courbures est fixé, et la distribution $df_G$ pour chacune des classes de courbures est évaluée. A l'étape suivante, chaque classe $\hat{\kappa}$ de courbure est subdivisée en 3. Rappelons que $\hat{\kappa}$ désigne indistinctement l'ensemble des courbures dans l'intervalle $[\hat{\kappa} - \Delta\kappa/2, \hat{\kappa} + \Delta\kappa/2]$ et le milieu de cet intervalle. Subdiviser cet intervalle en 3 permet donc de garder $\hat{\kappa}$ centré dans le nouvel intervalle à l'itération suivante et de suivre la même classe de courbure au fur et à mesure que le nombre $N_c$ de classes augmente (ou que la taille de classe $\Delta\kappa$ diminue sachant que $\Delta\kappa = [\max(I_\kappa) - \min(I_\kappa)]/N_c$). Les tests présentés ont été effectués sur un ellipsoïde de paramètres $\{d_{r0} = 1/256, p = 0\%, \alpha = 140\%\}$. La méthode utilisée pour évaluer les courbures est SR-$\mathcal{N}_{1-m}-0$. La résolution de courbure $d_{r0} = 1/256$ permet de s'affranchir des effets du maillage. Quant aux étirements $\alpha$ inférieures, les conclusions sont les mêmes.

Trois classes de courbure $\hat{\kappa} \in \{1.65, 3.3, 4.94\}$ sont étudiées. Ces classes sont caractéristiques des faibles courbures, des courbures intermédiaires et des fortes courbures rencontrées sur cet ellipsoïde. Les conclusions sont les mêmes pour les autres classes $\hat{\kappa}$ de courbure. L'évolution de la distribution $df_G$ pour chaque classe montre une droite constante caractéristique d'une indépendance en fonction du nombre de classe $N_c$ (ou de la taille de classe $\Delta\kappa$). Nous présentons ces résultats sur la figure 4.8. Nous remarquons que, quelle que soit la classe concernée, pour un nombre de classes compris entre 10 et $10^4$, $df_G$ reste constant. Toutefois, des variations apparaissent pour un faible nombre de classes ($N_c = 3$) de courbures et également pour un nombre élevé ($N_c > 10^4$). Dans le premier cas, cela s'explique par une discrétisation grossière de l'intervalle des courbures. Dans le second cas, cela est dû à la trop grande quantité de classe par rapport à la finesse du maillage de l'interface. Le nombre de points dans une classe tend alors
<figure><img src="image_72.png" /><figcaption>FIGURE 4.8 – Évolution de $df_G$ en fonction du nombre $N_c$ de classes de courbure pour un ellipsoïde de paramètres $d_{r_0} = 1/256$, $p = 0\%$, $\alpha = 140\%$.</figcaption></figure>

vers l'unité, ce qui engendre de forte variation dans l'évolution de $df_G$ et une sous-représentation statistique des classes considérées.

L'évolution de $df_G$ étant constante pour des valeurs raisonnables de $N_c$ (entre 10 et $10^4$), nous choisissons de fixer $N_c = 20$ pour la suite de l'étude. Un exemple de ce découpage est donné sur la figure 4.9. Nous aurons à travailler avec des ellipsoïdes de résolution beaucoup plus faible que celle présentée actuellement. Ainsi, la valeur de $N_c = 20$ qui fait partie des limites inférieures pour la résolution de courbure $d_{r_0} = 1/256$ deviendra une limite supérieure au fur et à mesure que la résolution du maillage baisse. Prendre une valeur plus grande dégraderait les résultats de la $df_G$ sur des maillages grossiers.

<figure><img src="image_73.png" /><figcaption>FIGURE 4.9 – Représentation des $N_c = 20$ classes de courbures sur un ellipsoïde d'étirement $\alpha = 20\%$ et de résolution de courbure. Dans la figure de droite, un zoom sur une classe de courbure (petit rectangle noir) permet de voir le maillage de l'interface.</figcaption></figure>

## Choix de la distribution et de la représentation

Comme mentionné dans la section 4.2.1, le choix de la fonction $G$ permet de distinguer les distributions surfaciques ($df_S$) des distributions s'apparentant à un dénombrement ($df_1$). Toutefois, la question est de savoir laquelle de ces distributions est la plus appropriée pour caractériser des objets et reste indépendante du traitement opéré sur le maillage de l'interface.

Pour représenter les distributions $df_G$, nous les traçons en fonction de $\hat{\kappa} \times r_0$, avec $r_0$ le rayon de la sphère de volume équivalent et $\hat{\kappa}$ les classes de courbure. Cela permet de comparer des ellipsoïdes de volume différent et de déterminer leur écart à la sphère de volume équivalent. En
<figure><img src="image_74.png" /><figcaption>FIGURE 4.10 – Sur-enrichissement de la zone de forte courbure pour un ellipsoïde de paramètre $d_{r_0} = 1/16$, $p = 0\%$ et $\alpha = 20\%$.</figcaption></figure>

effet, il est à remarquer que $\kappa_{max} \times r_0 = (1 + \alpha)^2$, avec $\kappa_{max} = \hat{\kappa}_{max} + \Delta\kappa/2$. Ce choix a donc l'avantage de déterminer l'étirement $\alpha$ de l'ellipsoïde. Par ailleurs, intégrer $df_G$ sur $\hat{\kappa} \times r_0$ nous donne $r_0$ et permet ainsi de remonter au volume de l'ellipsoïde.

<figure><img src="image_75.png" /><figcaption>FIGURE 4.11 – Distribution surfacique $df_S$ (trait continu) et distribution de décompte $df_1$ (trait discontinu) pour deux stratégies de maillage de l'interface sur un ellipsoïde de paramètres $d_{r_0} = 1/16$, $p = 0\%$, $\alpha = 20\%$.</figcaption></figure>

Quant au choix de la distribution, sachant que l'une des distributions se base sur des surfaces élémentaires construites autour des marqueurs de l'interface, et l'autre sur la quantité de marqueurs, un test pour les discriminer doit jouer sur ces facteurs. Le test mis en place consiste préalablement à sur-enrichir une zone quelconque d'une sphère. Cette sphère servira ensuite à générer la famille d'ellipsoïdes d'après la méthode décrite dans la section 4.2.2. La figure 4.10 montre un sur-enrichissement sur la zone de forte courbure pour un ellipsoïde de résolution de courbure $d_{r_0} = 1/16$, de perturbation de maillage $p = 0\%$ et d'étirement $\alpha = 20\%$. Les résultats de ce test sont présentés sur la figure 4.11. La distribution surfacique $df_S$ est en traits continu et la distribution de décompte $df_1$ en traits discontinus. En bleu sont les résultats obtenus sans la zone de sur-enrichissement et en rouge avec cette zone. Dans les deux cas de figure, les distributions surfaciques $df_S$ sont identiques. Une différence apparaît cependant lorsqu'on compare les distributions de décompte $df_1$. En effet, le sur-enrichissement a été effectué dans
une zone de forte courbure. Ainsi, le nombre de marqueurs dans les classes $\hat{\kappa}$ de forte courbure est plus important que lorsqu'il n'y a pas ce sur-enrichissement initial. La distribution $df_1$ pour les trois dernières classes de courbure $\hat{\kappa} \in \{1.37, 1.40, 1.43\}$ dans le cas du sur-enrichissement est largement supérieure à son équivalent dans le cas "classique" (Fig. 4.11). En contrepartie, pour les autres classes de courbures, la distribution $df_1$ du cas "sur-enrichi" est en dessous de la distribution du cas "classique".

Ce résultat montre la dépendance de l'allure de la distribution $df_1$ relativement à la stratégie de gestion du maillage de l'interface. Concernant les autres étirements $\alpha$, les conclusions restent inchangées. Par conséquent, nous utiliserons par la suite que les distributions surfaciques $df_S$ car contrairement aux distributions $df_1$, elles ne sont pas sensibles à la stratégie de gestion du maillage. La distribution $df_1$ n'est donc pas pertinente pour caractériser des structures.

## Calcul d'erreurs

Pour quantifier les erreurs, les courbures moyennes $\kappa_m$ des marqueurs $m$ utilisées dans les relations (4.10-4.11) sont évaluées grâce aux méthodes de Reconstruction de Surface. Le terme $\kappa_m^{ref}$ est employé lorsque la courbure est calculée de manière exacte. Quant au terme $df_S^{ref}$, il désigne la distribution calculée de manière exacte. Les étapes du calcul sont répertoriés dans l'annexe F. Afin de quantifier la précision des méthodes de Reconstruction de Surface sur l'évaluation de $\kappa$ et de la distribution surfacique $df_S$, nous calculons les erreurs relatives en norme $L_2$ et $L_\infty$. Deux familles d'erreurs sont introduites :

- les erreurs dites "globales" qui s'intéressent à l'ensemble des classes de courbures $\hat{\kappa}$ et à l'ensemble des courbures $\kappa_m$ des marqueurs (4.14-4.15). Elles nous informent du comportement des distributions et de l'évaluation des courbures sur l'ensemble du maillage des ellipsoïdes.

$$ \text{Err}_2^{\text{rel}}(df_S) = \left[ \frac{\sum_{\hat{\kappa} \subset I_{\kappa}} \left[ df_S(\hat{\kappa}) - df_S^{\text{ref}}(\hat{\kappa}) \right]^2}{\sum_{\hat{\kappa} \subset I_{\kappa}} \left( df_S^{\text{ref}}(\hat{\kappa}) \right)^2} \right]^{1/2}, \quad \text{Err}_{\infty}^{\text{rel}}(df_S) = \frac{\max_{\hat{\kappa} \subset I_{\kappa}} |df_S(\hat{\kappa}) - df_S^{\text{ref}}(\hat{\kappa})|}{\max_{\hat{\kappa} \subset I_{\kappa}} |df_S^{\text{ref}}(\hat{\kappa})|} \quad (4.14) $$

$$ \text{Err}_2^{\text{rel}}(\kappa) = \left[ \frac{\sum_{\kappa_m \in I_\kappa} \left[ \kappa_m - \kappa_m^{\text{ref}} \right]^2}{\sum_{\kappa_m \in I_\kappa} \left( \kappa_m^{\text{ref}} \right)^2} \right]^{1/2}, \quad \text{Err}_{\infty}^{\text{rel}}(\kappa) = \frac{\max_{\kappa_m \in I_\kappa} |\kappa_m - \kappa_m^{\text{ref}}|}{\max_{\kappa_m \in I_\kappa} |\kappa_m^{\text{ref}}|} \quad (4.15) $$

- les erreurs dites "locales". Dans ce cas de figure, l'attention est portée à une classe $\hat{\kappa}$ de courbure donnée et à l'ensemble des courbures $\kappa_m$ des marqueurs de cette classe (4.16-4.17). Elles nous permettent de voir la précision des méthodes d'évaluation des courbures pour chaque niveau de courbures rencontrés (faibles, intermédiaires, fortes) mais également la précision des distributions surfaciques sur chaque classe $\hat{\kappa}$ de courbure.

$$ \text{Err}^{\text{rel}}(df_S)|_{\hat{\kappa}} = \frac{|df_S(\hat{\kappa}) - df_S^{\text{ref}}(\hat{\kappa})|}{|df_S^{\text{ref}}(\hat{\kappa})|} \quad (4.16) $$

$$ \text{Err}_2^{\text{rel}}(\hat{\kappa}) = \left[ \frac{\sum_{\kappa_m \in \hat{\kappa}} \left[ \kappa_m - \kappa_m^{\text{ref}} \right]^2}{\sum_{\kappa_m \in \hat{\kappa}} \left( \kappa_m^{\text{ref}} \right)^2} \right]^{1/2}, \quad \text{Err}_{\infty}^{\text{rel}}(\hat{\kappa}) = \frac{\max_{\kappa_m \in \hat{\kappa}} |\kappa_m - \kappa_m^{\text{ref}}|}{\max_{\kappa_m \in \hat{\kappa}} |\kappa_m^{\text{ref}}|} \quad (4.17) $$
### 4.2.3 Résultats et analyses

#### Distribution $df_S$

Nous rappelons que nous travaillons sur des familles d'ellipsoïdes construites en nous servant des procédures décrites par les relations (4.12) et (4.13). Ainsi pour l'ensemble des données qui seront présentées par la suite, le rayon initial est $r_0 = 1$, la perturbation initiale du maillage est $p = 0\%$ et l'étirement est $\alpha \in \{20\%, 80\%, 140\%\}$. La figure 4.12 présente l'évolution de la distribution surfacique $df_S$ en fonction du paramètre adimensionnel $\hat{\kappa} \times r_0$. La résolution de la courbure est $d_{r0} = 1/256$.

<figure><img src="image_76.png" /><figcaption>FIGURE 4.12 – Distribution surfacique $df_S$ pour des ellipsoïdes de paramètres $d_{r0} = 1/256$, $p = 0\%$, $\alpha \in \{20\%, 80\%, 140\%\}$.</figcaption></figure>

Nous pouvons observer que pour un faible étirement $\alpha = 20\%$, la distribution surfacique reste assez compacte, avec une plage de courbure comprise entre 0.9 et 1.5. En effet, pour de faibles étirements (voir fig. 4.7b), l'ellipsoïde reste assez proche de la sphère dont la distribution est un pic à 1. Les courbures dans le plan du demi-petit axe et autour de cette région commencent à diminuer (d'où l'apparition de classes de courbure inférieures à 1) mais restent prédominantes ; la fraction de surface qui représente cette classe de courbure est élevée (pour $\hat{\kappa} \times r_0$, $\alpha = 20\%$, $df_S = 10$). En revanche, des courbures de plus en plus importantes sont engendrées sur les extrémités. Toutefois elles restent minoritaires en quantité de surface. Plus on étire la structure ($\alpha \in \{80\%, 140\%\}$), plus la queue de la distribution s'étend vers les fortes courbures, dénotant leur apparition aux extrémités de l'ellipsoïde. En terme de surface, les fortes courbures demeurent minoritaires par rapport aux faibles courbures qui restent autour de la valeur $\hat{\kappa} = 1$.

Maintenant que nous connaissons l'évolution de la distribution surfacique $df_S$ en fonction de l'étirement $\alpha$ de l'ellipsoïde, nous allons nous focaliser sur la précision des méthodes de Reconstruction de Surface SR pour évaluer les courbures et les distributions tout au long du processus de formation des ellipsoïdes. Les méthodes utilisées sont $\text{SR-}\mathcal{N}_1\text{-}\mathbf{n}_m\text{-}a_{00}$, $\text{SR-}\mathcal{N}_1\text{-}\mathbf{n}_m\text{-}0$ (Glimm), $\text{SR-}\mathcal{N}_2\text{-}\mathbf{n}_m\text{-}a_{00}$ et $\text{SR-}\mathcal{N}_2\text{-}\mathbf{n}_m\text{-}0$. Seule la méthode $\text{SR-}\mathcal{N}_1\text{-}\mathbf{n}_m\text{-}0$ (Glimm) est présentée, les conclusions étant les mêmes pour les autres méthodes. Nous nous intéressons dans un premier temps aux erreurs "globales" décrites par les relations (4.14-4.15).
## Erreurs "globales"

Pour l'étude de la convergence en maillage des courbures et des distributions surfaciques, nous devons caractériser la finesse de chaque maillage. Pour ce faire, le paramètre $d_{min} \times \kappa_{max}$ est utilisé plutôt que la résolution de la courbure $d_{r0}$. En effet, $d_{min} \times \kappa_{max}$ prend en compte l'impact de l'ensemble des procédures de gestion de maillage dans la création des différentes ellipsoïdes. Il correspond à la taille minimale des triangles de l'interface ramenée au rayon de courbure le plus petit. Concernant $d_{min}$, il est obtenu en se servant de la classe de courbure la plus élevée. Elle regroupe les plus fortes courbures rencontrées sur l'ellipsoïde. Connaissant la surface occupée par cette classe de courbure et en supposant que tous les triangles y sont équilatéraux, on en déduit la taille minimale $d_{min}$. Quant à $\kappa_{max}$, elle correspond à la courbure analytique la plus grande sur l'ellipsoïde.

<figure><img src="image_77.png" /><figcaption>FIGURE 4.13 – Erreurs relatives en norme $L_2$ (trait continu) et $L_\infty$ (trait discontinu) sur la courbure $\kappa$ en fonction de $d_{min} \times \kappa_{max}$ pour des ellipsoïdes de paramètres $p = 0\%$, $\alpha \in \{20\%, 80\%, 140\%\}$.</figcaption></figure>

— Sur la courbure $\kappa$: Le premier point sur lequel nous nous concentrons est la convergence de la méthode SR-$\mathcal{N}_1$-$\mathbf{n}_m$-0 (Glimm) en fonction du maillage, et ceci pour des ellipsoïdes de différents étirements $\alpha$. La figure 4.13 présente les erreurs relatives sur la courbure $\kappa$ en norme $L_2$ (traits continus) et $L_\infty$ (traits discontinus) en fonction du paramètre $d_{min} * \kappa_{max}$. Pour un faible étirement $\alpha = 20\%$, l'ordre de convergence en norme $L_2$ est de 2. En effet, l'ellipsoïde avec un étirement $\alpha = 20\%$ reste semblable à une sphère et le maillage initial est ainsi peu perturbé par les procédures de gestion du maillage (Fig. 4.7b). Nous retrouvons l'ordre de convergence pour les méthodes SR dans le cas d'un maillage non perturbé ($p = 0\%$) de la section 3.3.2, page 85. L'augmentation de l'étirement des ellipsoïdes ($\alpha \in \{80\%, 140\%\}$) induit une plus grande fréquence d'application des méthodes de gestion du maillage et par conséquent un impact plus important de la qualité du maillage sur les courbures du maillage (Fig. 4.7c et 4.7d). De cette modification du maillage résulte une baisse de l'ordre de convergence de l'erreur relative sur la courbure : elle converge avec un ordre de 1. Ce résultat est conforme avec les conclusions de convergence des méthodes SR pour une perturbation du maillage $p = 40\%$ (voir Sec. 3.3.2, page 85). Par ailleurs, les courbes pour les ellipsoïdes d'étirement $\alpha \in \{80\%, 140\%\}$ présentent une cassure dans leur allure (à $d_{min} \times \kappa_{max} = 3.5 \times 10^{-2}$ pour $\alpha = 140\%$, et à $d_{min} \times \kappa_{max} = 5.5 \times 10^{-2}$ pour $\alpha = 80\%$). Cet incident peut être corrélaté à l'application
des procédures de gestion de maillage. En effet, lorsqu'on analyse l'évolution du nombre
de marqueurs sur les ellipsoïdes (Tab. 4.2, page 114), il y a une consistance dans leur
évolution pour l'ellipsoïde de paramètre α = 20%; le nombre de marqueurs est multiplié
par 4 lorsque la résolution de courbure $d_{r0}$ est divisée par 2. On retrouve ce facteur 4
pour $\alpha \in \{80\%, 140\%\}$ et les 5 résolutions les plus fines. Ce coefficient baisse cependant
sur les résolutions grossières ($1/d_{r0} \in \{4, 8, 16\}$). Ainsi pour résumer, $\forall d_{r0}$ et $\alpha = 20\%$,
les procédures de gestion du maillage se sont toutes bien enclenchées. Tandis que pour
$\alpha \in \{80\%, 140\%\}$, deux phases sont observables. Pour $1/d_{r0} \in \{4, 8, 16\}$, les critères pour
la mise en œuvre des procédures de gestion de maillage n'étaient pas tous satisfaits. Pour
$1/d_{r0} \in \{32, 64, 128, 256, 512\}$, les critères étaient pleinement remplis pour un fonction-
nement optimal des procédures. Ce comportement est une explication de la rupture dans
l'allure des courbes d'erreurs sur la figure 4.13. On retrouvera cette cassure dans certaines
courbes à venir.

<figure><img src="image_78.png" /><figcaption>FIGURE 4.14 – Erreurs relatives en norme $L_2$ (trait continu) et $L_\infty$ (trait discontinu) sur la distribution $df_S$ en fonction de $d_{min} \times \kappa_{max}$ pour des ellipsoïdes de paramètres $p = 0\%$, $\alpha \in \{20\%, 80\%, 140\%\}$.</figcaption></figure>

— Sur la distribution surfacique *df*<sub>S</sub>: Les erreurs relatives sur la distribution surfacique *df*<sub>S</sub> sont représentées sur la figure 4.14. Quelle que soit l'étirement α de l'ellipsoïde, toutes les erreurs relatives convergent avec un ordre global de 2. Toutefois, ces erreurs diminuent avec l'augmentation de l'étirement, avec respectivement pour α ∈ {20%, 80%, 140%}, l'erreur relative Err<sup>rel</sup>(*df*<sub>S</sub>) ∈ {10<sup>-3</sup>, 4×10<sup>-4</sup>, 8×10<sup>-5</sup>} en norme L<sub>2</sub> pour d<sub>min</sub>×κ<sub>max</sub> = 10<sup>-2</sup>. Une explication de cette baisse de l'erreur avec l'augmentation de l'étirement α est la variation de la taille Δκ des classes de courbure κ̂. En effet, le nombre de classe N<sub>c</sub> de courbures restant constant, la variation de l'étirement α des ellipsoïdes fait changer les courbures maximale max(*I*<sub>κ</sub>) et minimale min(*I*<sub>κ</sub>) sur ces dernières, avec Δκ<sub>α=20%</sub> < Δκ<sub>α=80%</sub> < Δκ<sub>α=140%</sub>. Ainsi, pour un marqueur *m* de l'interface et pour un même niveau d'erreur commis sur l'évaluation de sa courbure κ<sub>m</sub> quel que soit l'étirement α, il est plus susceptible de changer de classe de courbure κ̂ sur de faibles étirements α que sur des étirements beaucoup plus élevés. Par conséquent, le numérateur des relations (4.14) devient beaucoup plus faible avec l'augmentation de l'étirement des ellipsoïdes. Cette baisse s'observe également au dénominateur de la relation de calcul des erreurs relatives sur la distribution surfacique, mais reste cependant de moins grande ampleur que celle
du numérateur. Nous pouvons voir cette baisse du numérateur avec l'étirement $\alpha$ sur la figure 4.15. Elle représente la contribution de chaque classe $\hat{\kappa}$ de courbure dans le calcul de l'erreur relative, en d'autres termes $[df_S(\hat{\kappa}) - df_S^{ref}(\hat{\kappa})]^2 / \sum_{\hat{\kappa} \subset I_{\kappa}} df_S^{ref}(\hat{\kappa})^2 = f(\hat{\kappa})$. Elle permet également de voir que les contributions ont une tendance globale descendante lorsqu'on se déplace des faibles classes de courbures vers les plus fortes. En effet, plus on avance vers les fortes classes de courbures, plus les fractions de surface qu'elles occupent deviennent négligeables. Par conséquent, même s'il pourrait y avoir des variations relatives importantes de surface sur les fortes classes de courbures, les variations absolues sont trop petites pour peser face à celles des faibles classes de courbures. Ainsi, dans le cas de l'ellipsoïde d'étirement $\alpha = 80\%$, la classe de courbure $\hat{\kappa} \times r_0 = 1$ a une contribution de $10^{-11}$ alors que celle à $\hat{\kappa} \times r_0 = 3.2$ a une contribution de $2 \times 10^{-16}$ (Fig. 4.15).

<figure><img src="image_79.png" /><figcaption>FIGURE 4.15 – Contribution de chaque classe $\hat{\kappa}$ de courbure dans le calcul de $Err^{rel}(df_S)$ pour un ellipsoïde de paramètres $d_{r_0} = 1/256$, $p = 0\%$</figcaption></figure>

L'étude des erreurs relatives sur la courbure $\kappa$ et sur la distribution surfacique $df_S$ nous montre que nous sommes globalement précis sur l'évaluation de ces deux quantités et que nous pouvons contrôler l'erreur commise lors de leur calcul. Ces résultats, somme toute encourageants, masquent la contribution de chaque classe de courbure, d'où l'idée d'étudier l'évolution des erreurs "locales" présentées dans les relations (4.16) et (4.17).

### Erreurs "locales"

L'attention est portée ici sur la précision de la distribution surfacique $df_S$ et également des méthodes d'évaluation des courbures, ceci sur chaque classe $\hat{\kappa}$ de courbure. Les tests ont été menés sur l'ensemble des $N_c = 20$ classes de courbures pour tous les ellipsoïdes. Nous présentons les résultats que sur trois classes de courbures. Ces classes de courbures, pour chaque ellipsoïde d'étirement $\alpha$ sont représentatives des faibles, intermédiaires et fortes classes de courbures rencontrées sur chaque ellipsoïde.

Les erreurs relatives "locales" sur les courbures $Err^{rel}(\hat{\kappa})$ et sur la distribution surfacique $Err^{rel}(df_S)|_{\hat{\kappa}}$ pour un ellipsoïde de paramètre $\alpha = 20\%$ sont présentées sur la figure 4.16. L'ensemble des classes de courbures sélectionnées $\hat{\kappa} \in \{0.8791, 1.08, 1.282\}$ converge à l'ordre 2 (Fig. 4.16a). Nous retrouvons ainsi l'ordre de convergence de 2 illustré dans le cas de l'étude
<figure><img src="image_80.png" /><figcaption>FIGURE 4.16 – Erreurs locales sur la courbure $\kappa$ (4.16a) et sur la distribution surfacique $df_S$ (4.16b) pour un ellipsoïde de paramètre $p = 0\%$, $\alpha = 20\%$.</figcaption></figure>

sur les erreurs globales (Fig. 4.13). Ainsi quel que soit le type de classe de courbure faible, in-
termédiaire ou forte rencontré sur cet ellipsoïde, les méthodes SR de Reconstruction de Surface
ont une bonne évaluation des courbures. Il faut toutefois mentionner que les courbures sont
globalement faibles et restent toutes du même ordre, centrées autour de 1.

Quant aux erreurs " locales " sur la distribution surfacique dfS (Fig. 4.16b), leurs ordres de
convergence de 2 viennent confirmer la tendance observée sur l'erreur globale (Fig. 4.14). Pour
dmin × κmax > 1.5 × 10⁻², l'erreur croît lorsqu'on passe des faibles classes de courbures aux fortes
classes de courbures. Cependant, lorsque dmin × κmax < 1.5 × 10⁻², une inversion s'opère entre
les classes ˆκ = 1.08 et ˆκ = 1.282. Cela peut s'expliquer par un effet dit de "vases communicants".
En effet, un marqueur m avec une courbure κm mal-évaluée peut basculer ainsi d'une classe de
courbure vers une autre classe de courbure adjacente. Cet effet combiné entre toutes les classes
de courbures peut ainsi engendrer une augmentation ou une baisse des erreurs locales commises
sur le calcul des distributions surfaciques.

Concernant les ellipsoïdes de paramètre $\alpha \in \{80\%, 140\%\}$, les courbes d'erreurs "locales" sont présentées respectivement sur les figures 4.17 et 4.18. Leur allure est en accord avec les conclusions de la section précédente.

<figure><img src="image_81.png" /><figcaption>FIGURE 4.17 – Erreurs locales sur la courbure $\kappa$ (4.17a) et sur la distribution surfacique $df_S$ (4.17b) pour un ellipsoïde de paramètre $p = 0\%$, $\alpha = 80\%$.</figcaption></figure>

Conclusion

Deux types de distributions ont été présentées : la distribution de décompte $df_1$ donnant la fraction de population de marqueurs compris dans une classe de courbure $\hat{\kappa}$, et la distribution $df_S$ qui permet d'avoir la fraction de surface associée à une classe de courbure $\hat{\kappa}$ donnée. Il s'avère cependant que la distribution de décompte $df_1$ est tributaire de la stratégie de gestion du maillage de l'interface des objets étudiés, alors que la distribution surfacique $df_S$ y est indépendante. Cette dernière est donc la seule appropriée pour caractériser des structures (ici des ellipsoïdes) car,
<figure><img src="image_82.png" /><figcaption>FIGURE 4.18 – Erreurs locales sur la courbure $κ$ (4.18a) et sur la distribution surfacique $df_S$ (4.18b) pour un ellipsoïde de paramètre $p = 0\%$, $α = 140\%$.</figcaption></figure>

en plus d'être indépendante de la stratégie de maillage, elle possède un ordre de convergence global de 2, donnant ainsi un contrôle sur les niveaux d'erreurs en fonction de la finesse du maillage. Par ailleurs, les conclusions sur la précision des méthodes de Reconstruction de Surface sont confortées dans cette partie. Disposant donc d'un outil précis pour la caractérisation de structures, il serait intéressant de se pencher sur des objets beaucoup plus complexes.
## 4.3 Ellipsoïdes déformées par des fonctions sinueuses et variqueuses

<figure><img src="image_83.png" /><figcaption>FIGURE 4.19 – Diagramme de Clift [17]</figcaption></figure>

Les structures ligamentaires que l'on peut rencontrer dans un écoulement diphasique ne se limitent pas qu'à des formes simples. En effet, les travaux de CLIFT et al. [17] et BHAGA et al. [9] mettent en lumière la variété de formes prises par des bulles ou gouttes en fonction des nombres de Reynolds Re, Eötvös Eo et Morton Mo dans leur mouvement dans des liquides uniquement sous l'action de la gravité (Fig. 4.19). Dans cette partie, nous nous intéressons plus particulièrement au régime sphérique, ellipsoïdal et wobbling. Ce dernier sera construit en introduisant des perturbations sinueuses et variqueuses de l'interface sur des ellipsoïdes. L'optique est de pouvoir distinguer ces structures variqueuses et sinueuses de leurs équivalents ellipsoïdales par le truchement de leur distribution surfacique de courbure.

Les déformations sinueuses et variqueuses sont appliquées sur des ellipsoïdes construits selon les relations (4.12) et (4.13). Les paramètres des ellipsoïdes sont {$p = 0\%$, $\alpha = 200\%$}. Soit $\bar{x}_m$ les coordonnées d'un marqueur $m$ de l'interface des ellipsoïdes. Les nouvelles coordonnées $\bar{x}_m^{new}$ des marqueurs $m$ de l'interface pour une déformation sinueuse sont données par la relation (4.18). Quant à la déformation variqueuse, elle est résumée dans la relation (4.19).

$$ \bar{x}_m^{new} = \begin{cases} x = a \sin(\theta) \cos(\varphi) \\ y = b \sin(\theta) \sin(\varphi) + q(x) \\ z = b \cos(\theta) \end{cases} \quad \text{(sinueux)} \quad (4.18) $$

$$ \bar{x}_m^{new} = \begin{cases} x = a \sin(\theta) \cos(\varphi) \\ y = b \sin(\theta) \sin(\varphi)[1 + q(x)] \\ z = b \cos(\theta)[1 + q(x)] \end{cases} \quad \text{(variqueux)} \quad (4.19) $$

avec

$$ \forall s \in \mathbb{R}, \quad q(s) = A \sin[2\pi(s - \phi)/\lambda] \quad (4.20) $$
$\lambda$ est la longueur d'onde de la déformation, $A$ est son amplitude et $\phi$ le déphasage. Tous ces paramètres sont fonction du demi-grand axe $a$ et du demi-petit axe $b$. Ainsi soient $\{p_{\lambda}, p_{A}, p_{\phi}\} \subset \mathbb{R}$. Ils caractérisent respectivement l'intensité des paramètres $\lambda, A, \phi$ comme indiqué dans la relation (4.21).

$$ \lambda = p_{\lambda} \times a, \quad A = \begin{cases} p_{A} \times b & (\text{sinueux}) \\ p_{A} & (\text{variqueux}) \end{cases}, \quad \phi = p_{\phi} \times \lambda \quad (4.21) $$

Par la suite, pour plus de simplicité, seuls les paramètres $\{p_{\lambda}, p_{A}, p_{\phi}\}$ sont employés pour caractériser les déformations sinueuses/variqueuses. Pour remonter à la longueur d'onde, à l'intensité et au déphasage, on se réfère à la relation (4.21).

Les figures 4.20 à 4.23 illustrent les différentes déformations issues du choix des paramètres $\{p_{\lambda}, p_{A}, p_{\phi}\}$. Les paramètres des déformations prennent les valeurs $p_{\lambda} \in \{50\%, 100\%, 200\%\}$, $p_{A} \in \{5\%, 20\%\}$ et $p_{\phi} \in \{0\%, 25\%\}$. Ainsi $p_{\lambda} = 100\%$ implique que la périodicité spatiale de la déformation équivaut au demi-grand axe $a$. Si $p_{A} = 20\%$, l'amplitude de la déformation vaut 20% du demi-petit axe $b$ dans le cas sinueux. Dans le cas variqueux, c'est 20% de $r(x)$, avec $r(x) = b[1 - x^2/a^2]^{1/2}$. Cela correspond, à $x$ donné, au rayon du cercle compris dans le plan $yz$. Pour $p_{\phi} = 25\%$, le déphasage spatial de la déformation est égal à 25% de la longueur d'onde. Pour de grandes longueurs d'onde ($p_{\lambda} \in \{100\%, 200\%\}$) et de faibles amplitudes de déformations ($p_{A} = 5\%$), les structures déformées restent assez semblables à des ellipsoïdes (Fig. 4.20b, 4.20c, 4.21b, 4.21c, 4.22b, 4.22c, 4.23b, 4.23c). Dans le cas de faibles amplitudes $p_{A}$ de déformations, des différences apparaissent lorsque la longueur d'onde $p_{\lambda}$ est une fraction du demi-grand axe $a$ ($p_{\lambda} = 50\%$). Ces différences sont notables avec l'augmentation de l'amplitude $p_{A}$ de la déformation (Fig. 4.20d, 4.20e, 4.20f, 4.21d, 4.21e, 4.21f, 4.22d, 4.22e, 4.22f, 4.23d, 4.23e, 4.23f).

C'est dans ce cadre que nous choisissons de travailler davantage sur les structures sinueuses/variqueuses de paramètres $\{p_{\lambda} = 50\%, p_{A} = 20\%\}$. En effet, elles présentent le plus de disparité avec un ellipsoïde et peuvent être rattachées à des phénomènes physiques comme la perturbation transverse et sinusoïdale d'un jet, l'écoulement d'un fluide à forte viscosité ou encore la propagation de perturbations le long d'un jet. Les distributions surfaciques $df_S$ pour le cas sinueux sont premièrement abordées avant de passer à celles de la structure variqueuse.
<figure><img src="image_84.png" /><figcaption>FIGURE 4.20 – Déformations sinueuses ($p_{\phi} = 0\%$) appliquées à un ellipsoïde de paramètres {$d_{r0} = 1/32, p = 0\%, \alpha = 200\%$}.</figcaption></figure>

<figure><img src="image_85.png" /><figcaption>FIGURE 4.21 – Déformation sinueuse ($p_{\phi} = 25\%$) appliquée à un ellipsoïde de paramètres {$d_{r0} = 1/32, p = 0\%, \alpha = 200\%$}.</figcaption></figure>
<figure><img src="image_86.png" /><figcaption>FIGURE 4.22 – Déformation variqueuse ($p_{\phi} = 0\%$) appliquée à un ellipsoïde de paramètres $\{d_{r0} = 1/32, p = 0\%, \alpha = 200\%\}$.</figcaption></figure>

<figure><img src="image_87.png" /><figcaption>FIGURE 4.23 – Déformation variqueuse ($p_{\phi} = 25\%$) appliquée à un ellipsoïde de paramètres $\{d_{r0} = 1/32, p = 0\%, \alpha = 200\%\}$.</figcaption></figure>
4.3.1 Déformation sinueuse

La déformation sinueuse consiste en une translation des plans yz dans la direction y, en fonction de la position x (Fig. 4.20d, 4.21d). Dans un premier temps, la distribution surfacique de ces structures pour p<sub>φ</sub> = 0% est présentée et analysée. Par la suite, l'effet du déphasage φ est discuté.

Distribution surfacique *df*<sub>S</sub> pour *p*<sub>φ</sub> = 0%

<figure><img src="image_88.png" /><figcaption>FIGURE 4.24 – Distribution surfacique *df*<sub>S</sub> pour une structure sinueuse de paramètres {p<sub>λ</sub> = 50%, p<sub>A</sub> = 20%, p<sub>φ</sub> = 0%}. L'ellipsoïde de départ est définie par {d<sub>r<sub>0</sub></sub> = 1/256, p = 0%, α = 200%}.</figcaption></figure>

La distribution surfacique *df*<sub>S</sub> d'une structure sinueuse de paramètres {p<sub>λ</sub> = 50%, p<sub>A</sub> = 20%, p<sub>φ</sub> = 0%}, issue d'un ellipsoïde de paramètre {d<sub>r<sub>0</sub></sub> = 1/256, p = 0%, α = 200%} est illustrée sur la figure 4.24. Elle est représentée en couleur rouge, tandis que la couleur bleue représente la distribution surfacique d'un ellipsoïde de même volume et de même longueur. La distribution surfacique de la structure sinueuse présente 3 parties.

— La première zone, pour $\hat{\kappa} \times r_0 \in [-0.1, 3]$ (Fig. 4.24), ressemble à une parabole et regroupe la plus grande fraction de surface de la structure. Le maximum en échelle log est atteint pour $\hat{\kappa} \times r_0 = 1$ avec $df_S(\hat{\kappa}) = 8 \times 10^{-1}$. Elle correspond aux régions où se forment des points selles dus à la déformation sinueuse. Des courbures négatives, dues à des courbures principales de signe opposé, commencent également à être générées et leur intensité est fonction de l'amplitude $p_A$ de la déformation. Plus $p_A$ est élevé, plus les courbures dans les vallées et les bosses sont grandes en valeur absolue. Les régions comprises entre ces vallées et ces bosses gagnent également en surface.

— La deuxième zone couvre les classes de courbures $\hat{\kappa} \times r_0 \in [3, 9]$. On retrouve l'allure de la distribution surfacique associée à un ellipsoïde, à savoir une décroissance progressive de la fraction de surface de chaque classe $\hat{\kappa}$ de courbure au fur et à mesure qu'on se déplace vers les fortes courbures. Cette zone, proche des extrémités de la structure sinueuse, comprend encore quelques cols induits par la déformation sinueuse.
— La dernière zone ($\hat{\kappa} \times r_0 \in [9, 11]$) correspond aux extrémités de la structure sinueuse où la déformation accentue les fortes courbures par rapport à l'ellipsoïde de référence.

<figure><img src="image_89.png" /><figcaption>FIGURE 4.25 – Distributions surfaciques $df_S(\hat{\kappa}_+)$ et $df_S(\hat{\kappa}_-)$ pour une structure sinueuse de paramètres {$p_\lambda = 50\%, p_A = 20\%, p_\phi = 0\%$}. L'ellipsoïde de départ est définie par {$d_{r_0} = 1/256, p = 0\%, \alpha = 200\%$}.</figcaption></figure>

L'ensemble de ces régions de la distribution surfacique de la structure sinueuse permet de la différencier de celle d'un ellipsoïde. Afin de pousser plus loin l'analyse de la distribution surfacique de la structure sinueuse, un regard est porté aux distributions surfaciques associées aux courbures principales. En notant $\kappa_G$ la courbure de Gauss, les courbures principales s'obtiennent par la relation $\kappa_{\pm} = \kappa \pm [\kappa^2 - \kappa_G]^{1/2}$. Les courbures principales, notées $\kappa_+$ et $\kappa_-$ sont alors définies par l'équation (4.22).

$$ \kappa_{+} = \kappa + [\kappa^{2} - \kappa_{G}]^{1/2} \quad , \quad \kappa_{-} = \kappa - [\kappa^{2} - \kappa_{G}]^{1/2} \qquad (4.22) $$
Les distributions surfaciques $df_S(\hat{\kappa}_+)$ et $df_S(\hat{\kappa}_-)$ qui leur sont associées sont répertoriées en rouge sur la figure 4.25. En bleu, on retrouve les mêmes distributions $df_S(\hat{\kappa}_+)$ et $df_S(\hat{\kappa}_-)$ mais pour un ellipsoïde de référence de paramètres $\{d_{r_0} = 1/256, p = 0\%, \alpha = 200\%\}$.

- La distribution $df_S(\hat{\kappa}_+)$ de la structure sinueuse (Fig. 4.25b) est comparable à sa contrepartie ellipsoïdale sur l'intervalle $\hat{\kappa}_+ \times r_0 \in [2, 9]$, malgré un artefact en $\hat{\kappa}_+ \times r_0 = 3.5$. Hormis cet intervalle, elle présente une queue de distribution allongée qui dénote de l'augmentation des fortes courbures à l'extrémité de la structure sinueuse ($\hat{\kappa}_+ \times r_0 \in [9, 12]$). Il y a également la présence de classes de courbures plus faibles, avec $\hat{\kappa}_+ \times r_0 = 1.5$. Il est a noter que toutes les classes $\hat{\kappa}_+$ de courbure sont positives et non-nulles. Cela a pour implication que la structure ne présente pas de surface plane (les deux courbures principales devant être simultanément nulles). Ainsi, les déformations locales de l'intervalle sont pilotés par le signe de la courbure principale $\kappa_-.

- Sur la distribution $df_S(\hat{\kappa}_-)$ (Fig. 4.25b), les classes de courbures $\hat{\kappa}_- \times r_0 \in [-2, 2.5]$ regroupent la plus grande fraction de surface de la structure. La présence de classes de courbures négatives implique que localement des régions de la structure sinueuse sont des points-selles. Sur cette plage de classes de courbures, ce sont celles proches de 0 qui dominent. Par la suite, l'allure de la distribution $df_S(\hat{\kappa}_-)$ est semblable à celle de la référence ellipsoïdale, même si elle demeure légèrement inférieure, sur les classes de courbures $\hat{\kappa}_- \times r_0 \in [2.5, 9]$. Enfin, la queue de la distribution est légèrement allongée, avec l'apparition de plus de classes de fortes courbures sur $\hat{\kappa}_+ \times r_0 \in [9, 10]$.

L'impact de la perturbation sinueuse sur la structure ellipsoïdale de référence peut se résumer en trois parties. La zone 1 concerne les faibles classes de courbures. La distribution surfacique est fortement modifiée et il peut y avoir la création de classes de courbures plus faibles. La zone 2 regroupe les fortes classes de courbures et celles intermédiaires. Sur cette zone, la distribution surfacique de la structure sinueuse reste semblable à celle de l'ellipsoïde, avec peu voire pas de modification. Enfin la zone 3 rend compte de l'impact de la déformation sinueuse sur les fortes courbures de l'ellipsoïde. Distinguer ces zones et rendre compte de l'allure de la distribution surfacique sont des éléments qui permettent de caractériser la structure étudiée.

**Influence du déphasage $\phi$ : cas $p_{\phi} = 25\%$**

La section précédente nous a présenté l'allure des distributions surfaciques de courbure $\kappa$ et des courbures principales $\kappa_+$ et $\kappa_-$ ($df_S(\hat{\kappa})$, $df_S(\hat{\kappa}_+)$ et $df_S(\hat{\kappa}_-)$). Nous nous intéressons dans cette partie aux modifications induites sur ces distributions par le déphasage $\phi$. La figure 4.26 compare les distributions sus-mentionnées pour des déphasages $p_{\phi} = 0\%$ (bleu) et $p_{\phi} = 25\%$ (rouge).

Dans la zone 1 de faibles courbures, les distributions sont identiques. En effet, pour $\hat{\kappa} \times r_0 \in [-0.1, 2.2]$, nous retrouvons la forme "parabolique" (en échelle log) de la distribution $df_S(\hat{\kappa})$ (Fig. 4.26a). Sur la distribution $df_S(\hat{\kappa}_+)$, les pointes sont quasi-identiques sur $\hat{\kappa}_- \times r_0 \in [1.5, 2.2]$ avec une légère différence sur la fin de cet intervalle (Fig. 4.26b). Quant à la distribution $df_S(\hat{\kappa}_-)$, sur $\hat{\kappa}_- \times r_0 \in [-2, 2.5]$ nous avons la même forme de montagne avec une bosse pour les classes de courbures autour de 0 (Fig. 4.26c). Pour la zone 2 qui, globalement pour l'ensemble des distributions, couvre les classes de courbures sur $[3, 9]$, on observe la décroissance de la fraction de surface avec l'augmentation de la courbure. Les distributions pour le cas $p_{\phi} = 25\%$ sont supérieures à celles du cas $p_{\phi} = 0\%$. Elles restent toutefois du même ordre grandeur. Cette différence peut s'expliquer en examinant la zone 3 qui recouvre l'intervalle $[9, 12]$. En effet, la déformation sinueuse avec un déphasage $p_{\phi} = 25\%$ n'a que peu, voire pas d'impact sur les fortes courbures aux extrémités de l'ellipsoïde. La zone 3 se trouve alors vide, ce qui constitue un gain de surface pour les autres classes de courbures contrairement au cas $p_{\phi} = 0\%$.
<figure><img src="image_90.png" /><figcaption>(a) $df_S(\hat{\kappa})$</figcaption></figure>

<figure><img src="image_91.png" /><figcaption>(b) $df_S(\hat{\kappa}_+)$</figcaption></figure>

<figure><img src="image_92.png" /><figcaption>(c) $df_S(\hat{\kappa}_-)$</figcaption></figure>

FIGURE 4.26 – Comparaison des distributions surfaciques $df_S(\hat{\kappa})$, $df_S(\hat{\kappa}_+)$ et $df_S(\hat{\kappa}_-)$ pour $p_\phi \in \{0\%, 25\%\}$. La structure sinueuse a pour paramètres $\{p_\lambda = 50\%, p_A = 20\%\}$. L'ellipsoïde de départ est définie par $\{d_{r_0} = 1/256, p = 0\%, \alpha = 200\%\}$.
Conclusion

L'étude de la distribution surfacique de la structure sinueuse nous a conduit à distinguer 3 zones caractéristiques permettant de différencier la structure sinueuse d'un ellipsoïde. Ainsi la zone la plus caractéristique des modifications engendrées par la déformation et qui impacte une majeure partie de la surface de la structure est celle localisée aux niveaux des faibles classes de courbures. Sur la gamme de paramètres employés pour l'étude ($p_{\phi} \in \{0\%, 25\%\}$), cette zone est peu modifiée et constitue donc un repère important pour décrire la structure.

4.3.2 Déformation variqueuse

Dans cette partie, nous nous intéressons maintenant à la distribution surfacique d'une struc-
ture ellipsoïdale perturbée par une déformation variqueuse (4.19). Les paramètres de la déforma-
tion sont {p<sub>λ</sub> = 50%, p<sub>A</sub> = 20%, p<sub>φ</sub> = 0%}. Une image de la structure variqueuse est répertoriée
sur la figure 4.22d.

Distribution surfacique *df*<sub>S</sub> pour *p*<sub>*φ*</sub> = 0%

<figure><img src="image_93.png" /><figcaption>FIGURE 4.27 – Distribution surfacique *df*<sub>S</sub> pour une structure variqueuse de paramètres {p<sub>λ</sub> = 50%, p<sub>A</sub> = 20%, p<sub>φ</sub> = 0%}. L'ellipsoïde de départ est définie par {d<sub>r<sub>0</sub></sub> = 1/256, p = 0%, α = 200%}.</figcaption></figure>

La figure 4.27 présente en rouge la distribution surfacique *df*<sub>S</sub> pour une structure variqueuse de paramètres {p<sub>λ</sub> = 50%, p<sub>A</sub> = 20%, p<sub>φ</sub> = 0%}. En bleu, se trouve la distribution d'un ellipsoïde de référence et de paramètres {d<sub>r<sub>0</sub></sub> = 1/256, p = 0%, α = 200%}. La distribution surfacique de la structure variqueuse présente 2 parties.

— La première région couvre $\hat{\kappa} \times r_0 \in [0.2, 2.5]$. On distingue d'abord un plateau allant des faibles classes de courbures jusqu'à $\hat{\kappa} \times r_0 = 2$, qui montre une certaine régularité dans la déformation de la structure. Ensuite on observe une chute assez importante de la distribution $df_S$, passant de $5 \times 10^{-1}$ à $1.5 \times 10^{-2}$ sur $\hat{\kappa} \times r_0 \in [2, 2.5]$.

— La deuxième zone est sur l’intervalle $\hat{\kappa} \times r_0 \in [2.5, 8.75]$. La distribution surfacique de la structure variqueuse décroît en se déplaçant vers les fortes courbures, mais légèrement moins lentement que celle associée à l’ellipsoïde.
<figure><img src="image_94.png" /><figcaption>FIGURE 4.28 – Distributions surfaciques $df_S(\hat{\kappa}_+)$ et $df_S(\hat{\kappa}_-)$ pour une structure variqueuse de paramètres $\{p_\lambda = 50\%, p_A = 20\%, p_\phi = 0\%\}$. L’ellipsoïde de départ est définie par $\{d_{r_0} = 1/256, p = 0\%, \alpha = 200\%\}$.</figcaption></figure>

En ce qui concerne les distributions surfaciques $df_S(\hat{\kappa}_+)$ et $df_S(\hat{\kappa}_-)$ sur les courbures principales, elles sont regroupées sur la figure 4.28.

— Pour $\hat{\kappa}_+ \times r_0 \in [1.5, 4.5]$ (Fig. 4.28a), la distribution $df_S(\hat{\kappa}_+)$ oscille autour de celle de l'ellipsoïde. On a l'apparition d'une classe de courbure plus faible comparé au cas ellipsoïdal $(\hat{\kappa}_+ \times r_0 = 1.5)$. Sur le reste de l'intervalle $[4.5, 9]$, la distribution $df_S(\hat{\kappa}_+)$ suit la tendance de celle de l'ellipsoïde tout en restant supérieure à cette dernière. Cela est du aux oscillations observées dans le premier intervalle. Une partie des surfaces a été reportée sur les fortes classes de courbures. Par ailleurs, toutes les classes $\hat{\kappa}_+$ de courbures sont positives, d'où les types de déformations locales seront pilotés par le signe de la courbure principale $\kappa_-$.
— Contrairement à la précédente, la distribution $df_S(\hat{\kappa}_–)$ de la structure variqueuse est fortement modifiée sur les faibles classes de courbures, avec la création de classes de courbures négatives (Fig. 4.28b). Ainsi pour $\hat{\kappa}_– \times r_0 \in [-2, 2]$, on observe un plateau légèrement ascendant jusqu’à $\hat{\kappa}_– \times r_0 = 1.7$. Ce plateau traduit encore la régularité des déformations sur la structure variqueuse. A partir de $\hat{\kappa}_– \times r_0 = 1.5$ la distribution $df_S(\hat{\kappa}_–)$ de la structure variqueuse continue sa décroissance en suivant l’allure de celle associée à l’ellipsoïde.

Influence du déphasage $\phi$ : cas $p_{\phi} = 25\%$

Comme dans le cas sinueux, l'idée est de voir l'impact sur la distribution de la modification des paramètres de la déformation, et par conséquent d'identifier les zones de la distribution qui sont caractéristiques de la perturbation apportée. Ainsi, la figure 4.29 compare les distributions pour différents déphasages $p_{\phi} \in \{0\%, 25\%\}$.

Un premier point est l'amplification des classes de fortes courbures avec $p_{\phi} = 25\%$. Tout un ensemble de classes de courbures apparaissent ainsi sur l'intervalle $[9, 14]$ quelle que soit la distribution considérée. Ces classes supplémentaires par rapport au cas $p_{\phi} = 0\%$ constituent des surfaces en moins pour les autres classes. Ainsi, sur l'intervalle $[3, 9]$, toutes les distributions du cas $p_{\phi} = 0\%$ sont supérieures à celles du cas $p_{\phi} = 25\%$. Seuls sur les intervalles $\hat{\kappa} \times r_0 \in [0, 2]$ pour $df_S(\hat{\kappa})$ (Fig. 4.29a) et $\hat{\kappa} \times r_0 \in [-2, 2]$ pour $df_S(\hat{\kappa}_–)$ (Fig. 4.29c) on retrouve une similarité entre les deux cas $p_{\phi} \in \{0\%, 25\%\}$. Une légère différence se manifeste sur la distribution $df_S(\hat{\kappa})$ avec l'apparition de classes de courbures négatives à $\hat{\kappa} \times r_0 = -0.1$) (Fig. 4.29a). Cela implique qu'il y a des régions où les courbures principales sont du même ordre mais de signe opposé.
<figure><img src="image_95.png" /><figcaption>(a) $df_S(\hat{\kappa})$</figcaption></figure>

<figure><img src="image_96.png" /><figcaption>(b) $df_S(\hat{\kappa}_+)$</figcaption></figure>

<figure><img src="image_97.png" /><figcaption>(c) $df_S(\hat{\kappa}_-)$</figcaption></figure>

FIGURE 4.29 – Comparaison des distributions surfaciques $df_S(\hat{\kappa})$, $df_S(\hat{\kappa}_+)$ et $df_S(\hat{\kappa}_-)$ pour $p_\phi \in \{0\%, 25\%\}$. La structure variqueuse a pour paramètres $\{p_\lambda = 50\%, p_A = 20\%\}$. L'ellipsoïde de départ est définie par $\{d_{r_0} = 1/256, p = 0\%, \alpha = 200\%\}$.
Conclusion

En comparaison au cas sinueux, la déformation variqueuse avec $p_\phi = 0\%$ n'induit que la distinction de 2 zones sur les distributions surfaciques $df_S(\hat{\kappa})$ et $df_S(\hat{\kappa}_−)$, la distribution $df_S(\hat{\kappa}_+) $ se détachant très peu de celle de l'ellipsoïde de référence (Fig. 4.30). Ainsi, la première zone localisée sur les faibles classes de courbures présente une distribution assez plane, dénotant de la régularité de la déformation et de l'alternance entre monts et vallées sur la structure variqueuse. La deuxième zone cumule les classes intermédiaires de courbure et les fortes classes de courbures. Dans cette zone, la distribution suit l'évolution de celle associée à l'ellipsoïde. Il y a une décroissance progressive des fractions de surfaces lorsqu'on avance vers les fortes courbures. Sur les extrémités de la structure, des classes de courbures plus grandes ne sont pas engendrées pas la déformation, d'où l'inexistence d'une troisième zone. L'impact du déphasage $p_\phi = 25\%$ est en majeure partie répercutée sur la zone 3, celle des fortes classes de courbures. Ainsi, la structure variqueuse voit ses fortes classes de courbures augmenter tandis que la bulle sinueuse perd complètement cette zone de fortes classes de courbures (Fig. 4.31). La zone 1 des faibles classes de courbures demeure celle présentant le plus d'information pour différencier les bulles sinueuses et variqueuses. Quand à la zone 2 des classes de courbures intermédiaires, même s'il peut y avoir de légères différences par rapport à la référence ellipsoïdale, l'allure des distributions reste du même ordre.

Au vu des ces résultats préliminaires, nous sommes confortés dans notre choix d'utiliser les distributions surfaciques de courbures pour caractériser les interfaces. Pour pousser notre analyse, elles sont évaluées sur des bulles ascensionnelles issues des cas tests 1 et 2 de HYSING et al. [51] en 2D, avec une adaptation 3D par ADELSBERGER et al. [2]. Les configurations des tests sont rappelées brièvement. Le domaine de simulation est une boîte rectangulaire de dimensions $1 \times 2(\times 2 \text{ en 3D})$ comprenant une sphère de rayon $r_0 = 0.25$ centrée en $[0.5, 0.5]$ (ou $[0.5, 0.5, 0.5]$ en 3D). Des conditions de no-glissement sont imposées aux parois ($\mathbf{u} = 0$). Les paramètres physiques et adimensionnels des cas tests sont rassemblés dans le tableau 4.3.

<table><caption>TABLE 4.3 – Paramètres physiques et adimensionnels pour les cas tests de [51].</caption><thead><tr><th>Test case</th><th>ρ<sub>1</sub></th><th>ρ<sub>2</sub></th><th>μ<sub>1</sub></th><th>μ<sub>2</sub></th><th>g</th><th>σ</th><th>Re</th><th>Eo</th></tr></thead><tbody><tr><td>1</td><td>1000</td><td>100</td><td>10</td><td>1</td><td>0.98</td><td>24.5</td><td>35</td><td>10</td></tr><tr><td>2</td><td>1000</td><td>1</td><td>10</td><td>0.1</td><td>0.98</td><td>24.5</td><td>35</td><td>125</td></tr></tbody></table>

L'indice 1 est utilisé pour le fluide lourd environnant et l'indice 2 pour la bulle, avec $\rho$ la masse volumique, $\mu$ la viscosité dynamique, $g$ la pesanteur, $\sigma$ le coefficient de tension de surface et :

$$Re = \frac{\rho_1 U_g L}{\mu_1} \quad Eo = \frac{\rho_1 U_g^2 L}{\sigma} \quad U_g = \sqrt{gL} \quad L = 2r_0$$

Les figures 4.32 et 4.33 présentent respectivement les bulles du cas test 1 à $t = 3$ et du cas 2 à $t = 1.5$. Elles ont été générées par symétrie axiale à partir des résultats issus des simulations 2D de [51]. D'après le diagramme de Clift (Fig 4.19), la première bulle est de type ellipsoïdale. La deuxième bulle quant à elle se situe entre les régimes "skirted" et "dimpled ellipsoidal-cap". Les distributions surfaciques de courbure $df_S(\hat{\kappa})$ des deux bulles sont présentées sur la figure 4.34. Premier point, elles diffèrent des distributions surfaciques que nous avons obtenues jusqu'à présent. Concernant la distribution $df_S$ du cas test 1 (couleur bleue sur Fig. 4.34), nous avons une prédominance des classes de faibles courbures qui correspondent à la calotte et la base aplatie de la bulle. Lorsqu'on se déplace vers les fortes classes de courbures, la décroissance de leur distribution s'accélère. Les fortes classes de courbures ont une fraction de surface minoritaire
<figure><img src="image_98.png" /><figcaption>(a) $df_S(\hat{\kappa})$</figcaption></figure>

<figure><img src="image_99.png" /><figcaption>(b) $df_S(\hat{\kappa}_+)$</figcaption></figure>

<figure><img src="image_100.png" /><figcaption>(c) $df_S(\hat{\kappa}_-)$</figcaption></figure>

FIGURE 4.30 – Comparaison des distributions surfaciques $df_S(\hat{\kappa})$, $df_S(\hat{\kappa}_+)$ et $df_S(\hat{\kappa}_-)$ entre les structures sinueuse et variqueuse avec pour paramètres {$p_\lambda = 50\%, p_A = 20\%, p_\phi = 0\%$}. L'el-lipsoïde de départ est définie par {$d_{r_0} = 1/256, p = 0\%, \alpha = 200\%$}.
<figure><img src="image_101.png" /><figcaption>FIGURE 4.31 – Comparaison des distributions surfaciques $df_S(\hat{\kappa})$, $df_S(\hat{\kappa}_+)$ et $df_S(\hat{\kappa}_-)$ entre les structures sinueuse et variqueuse avec pour paramètres $\{p_\lambda = 50\%, p_A = 20\%, p_\phi = 25\%\}$. L'ellipsoïde de départ est définie par $\{d_{r_0} = 1/256, p = 0\%, \alpha = 200\%\}$.</figcaption></figure>
<figure><img src="image_102.png" /><figcaption>FIGURE 4.32 – Bulle 3D issue du cas test 1, à $t = 3$ [2, 51]</figcaption></figure>

<figure><img src="image_103.png" /><figcaption>FIGURE 4.33 – Bulle 3D issue du cas test 2, à $t = 1.5$ [2, 51]</figcaption></figure>

<figure><img src="image_104.png" /><figcaption>FIGURE 4.34 – Distributions surfaciques de bulles issues des cas tests 1 et 2 de [51].</figcaption></figure>

et sont localisés au niveau du renflement de la bulle. La distribution $df_S$ du cas test 2 (couleur rouge sur Fig. 4.34) est beaucoup plus étendue en terme de classes de courbures $\hat{\kappa}$. Les classes
dominantes en terme de surface sont centrées autour des valeurs $\hat{\kappa} = -1$ et $\hat{\kappa} = 1$. La bulle présente donc des régions convexes et concaves en proportion égale. Le passage entre ces deux régions nécessite l'annulation de la courbure d'où la présence de classes de courbures centrées autour de $\hat{\kappa} = 0$. Elles constituent la deuxième force en terme d'occupation de surface, place qu'elles partagent avec les courbures intermédiaires. Les fortes classes de courbures ferment naturellement le podium.

Ces tests opérés sur différentes structures ont montré la diversité des distributions surfaciques. C'est un premier pas dans l'élaboration d'une base de données pour caractériser et identifier les formes de bulles/gouttes dans des écoulements diphasiques.
# Bibliographie

## Identification d'objets grâce à leurs propriétés surfaciques

[2] J. ADELSBERGER, S. GROSS, M. KLITZ, A. RÜTTGERS et M. GRIEBEL. "3D Incompressible two-phase flow benchmark computations for rising droplets". In: *11th World Congress on Computational Mechanics (WCCM XI), 5th European Conference on Computational Mechanics (ECCM V), 6th European Conference on Computational Fluid Dynamics (ECFD VI)*. Sous la dir. d'E. ÕNATE, J. OLIVER et A. HUERTA. 2013 (cf. p. 138, 141).

[9] D. BHAGA et M. E. WEBER. "Bubbles in viscous liquids: shapes, wakes and velocities". In: *J. Fluid Mech.* 105 (1981), p. 61-85 (cf. p. 126).

[17] R. CLIFT, J. R. GRACE et M. E. WEBER. *Bubbles, Drops, and Particles*. Academic Press, 1978 (cf. p. 126, 146).

[26] J. DU, B. FIX, J. GLIMM, X. JIA, X. LI, Y. LI et L. WU. "A simple package for front tracking". In: *J. Comput. Phys.* 213 (2006), p. 613-628 (cf. p. 29, 33, 73, 78, 79, 83, 85, 90, 92, 107, 161, 163).

[51] S. HYSING, S. TUREK, D. KUZMIN, N. PAROLINI, E. BURMAN, S. GANESAN et L. TOBISKA. "Quantitative benchmark computations of two-dimensional bubble dynamics". In: *Int. J. Numer. Methods Fluids* 60.11 (2009), p. 1259-1288 (cf. p. 138, 141, 146).

[76] M. MEYER, M. DESBRUN, P. SCHRODER et A. H. BARR. "Discrete Differential-Geometry Operators for Triangulated 2-Manifolds". In: *Visualization and Mathematics III*. Sous la dir. de H.-C. HEGE et K. POLTHIER. Springer, 2003, p. 35-57 (cf. p. 33, 36, 75, 76, 81, 85, 90, 92, 107, 152, 154).

[123] G. TRYGGVASON, R. SCARDOVELLI et S. ZALESKI. *Direct Numerical Simulations of Gas-Liquid Multiphase Flows*. 1re éd. Cambridge University Press, 2011 (cf. p. 31, 77, 82, 85, 90, 93, 95, 107).
# Chapitre 5

## Conclusion générale et perspectives

Ce travail a eu pour objectif de développer des outils afin de caractériser les interfaces présentes au sein d'écoulements diphasiques à phases séparées, au travers de leurs courbures et normales. Outre la bonne évaluation des forces de tension de surfaces, caractériser les interfaces diphasiques lors de leur complexe évolution topologique peut permettre de remonter à des propriétés physiques de l'écoulement.

Ainsi, un des premiers éléments fut la détermination de la stratégie pour décrire les changements de l'interface. Le choix s'est porté sur la méthode du front-tracking. Elle a l'avantage de donner explicitement la position de l'interface et d'accéder aisément aux propriétés de la surface. Afin de suivre correctement tous les changements topologiques de l'interface, nous avons entièrement implémenté dans le code maison FUGU développé au sein de l'équipe TCM de MSME un ensemble de procédures de gestion de maillage. Elles sont réparties en deux groupes. Le premier vise a garder un maillage de l'interface adéquat grâce aux procédures d'enrichissement et de désenrichissement. Le second vise l'amélioration de la qualité du maillage pour avoir des triangles les plus équilatéraux possibles. Il regroupe les procédures de lissage, d'échanges d'arêtes et de suppression de triangles de très mauvaise qualité. Nous avons également mis en place une procédure permettant le contrôle du volume en usant des normales aux marqueurs de l'interface et de leurs vitesse. D'un point de vue qualitatif, nos procédures fournissent d'excellents résultats vis à vis du test proposé par LEVEQUE [65] et repris par ENRIGHT et al. [30], qui consiste à déformer une sphère dans un champ tourbillonnaire. Qualitativement, nous avons aussi de bons résultats car nous observons une convergence d'ordre 2 avec le maillage de l'erreur spatiale entre la sphère déformée et une de référence.

Le second volet de notre travail a concerné les méthodes utilisées dans le cadre du front-tracking pour évaluer les courbures et les normales aux marqueurs de l'interface. Ce travail a fait l'objet d'un article accepté dans Journal of Computational Physics. Trois familles de méthodes ont été comparées entre elles. La première se base sur la discrétisation de l'opérateur de Laplace-Beltrami (LBO). La deuxième s'appuie sur la formulation intégrale de la force de tension de surface (IF). Et la dernière correspond aux méthodes de reconstruction de surface (SR). Deux tests ont été menés. Le premier test dit statique est effectué sur une surface présentant une large variété de courbures, donnant lieu localement à la présence de surfaces de type parabolique, point-selle et cylindrique. Le second test dit dynamique est un test original qui prend en compte les schémas de transport lagrangien des marqueurs de l'interface. L'évolution de la surface à chaque instant est connue analytiquement et ses propriétés géométriques sont déduites aisément. Ce dernier test offre également la possibilité, selon le choix d'un paramètre, de se trouver dans des situations où nous avons soit une évolution "smooth" de la surface avec une variation brusque et localisée de la courbure (cas de l'impact d'une goutte), soit le développement d'un pincement (cas après la rupture d'une goutte). Dans le cas statique avec une surface perturbée ou pas, les méthodes SR s'avèrent les plus précises avec un ordre 2 de convergence en maillage pour l'option sans perturbation $p = 0\%$ et un ordre 1 lorsque $p = 20\%$. Elles fournissent également le meilleur
rapport erreur relative/temps CPU. Dans le cas dynamique, pendant que les méthodes LBO et IF peinent, voire ne convergent pas, les méthodes SR sont les seules à converger avec un ordre proche de 1. Elles ont donc été sélectionnées dans notre étude pour la dernière partie de notre travail.

Ce dernier volet a été consacré à l'introduction de la distribution surfacique $df_S$, à son étude et la caractérisation de différentes interfaces rencontrées dans des écoulements diphasiques (cas des bulles/gouttes de type sphérique, ellipsoïdale, wobbling [17]). Contrairement à son pendant $df_1$ orienté sur le dénombrement, elle est indépendante de la stratégie de gestion de maillage employée par l'utilisateur. Indépendante d'une discrétisation raisonnable de la plage de courbures (ou par rapport au nombre de classes de courbure $\hat{k}$), elle peut s'évaluer analytiquement sur des ellipsoïdes. Elle y converge à l'ordre 2 sur chaque classe de courbure $\hat{k}$ et aussi globalement avec le raffinement du maillage de l'interface. Pour une sphère, elle présente un pic localisé à la courbure de celle-ci. Pour des ellipsoïdes, elle a une allure en $1/x$, qui est d'autant plus étirée vers les fortes classes de courbures que l'ellipsoïde est allongé. Dans le cas de bulles de type wobbling, nous avons étudié deux types de déformation (sinueux et variqueux) dont les paramètres sont la longueur d'onde, l'amplitude et le déphasage de la perturbation. Ainsi pour de grandes longueurs d'onde et une faible amplitude, aucune démarcation n'a été repérée entre les distributions surfaciques de courbures des bulles déformées et leur contrepartie ellipsoïdale de volume et d'étirement identiques. Des différences apparaissent pour de petites longueurs d'onde et une grande amplitude de déformation. Trois parties se distinguent alors. La première zone, qui regroupe les petites classes de courbures, présente une forme "parabolique" dans le cas sinueux et un plateau dans le cas variqueux. Ces formes se retrouvent également en se penchant sur les distributions surfaciques des courbures principales, notamment la distribution $df_S(\hat{k}_-)$.
La seconde partie couvre les classes de courbures intermédiaires, avec des distributions qui ont du même ordre que ce qu'on observe sur les références ellipsoïdales. Quant à la troisième partie, elle concerne les fortes courbures. Elle est plus ou moins présente et prononcée en fonction du déphasage de la déformation. Entre les cas sinueux et variqueux, les différences les plus notables se situent sur les zones 1 et 3 des distributions, à savoir sur les petites classes de courbures et sur les plus grandes. Ces distinctions nous permettent d'identifier les deux cas de figure. Pour aller plus loin, les distributions de bulles issues des cas tests de HYSING et al. [51] ont aussi été effectuées et présentent des allures uniques et différentes de objets précédemment étudiés.

Ce travail se veut une étape préliminaire dans l'identification automatique d'interfaces. Au vu des études menées, une première perspective est l'expansion de l'étude paramétrique sur structures wobbling en prenant une gamme beaucoup plus étendue pour les paramètres que sont la longueur d'onde, l'amplitude et le déphasage de la déformation. Il serait intéressant de voir comment les distributions que nous avons obtenues et présentées jusqu'à maintenant évoluent avec ces paramètres. Par la suite, le périmètre des analyses peut être étendu, dans un premier temps, aux autres formes rencontrées dans le diagramme de Clift [17], puis à des formes plus générales en utilisant les harmoniques sphériques. Sur ce dernier point, nous avons déjà entamé un travail sur les calculs analytiques impliquant les harmoniques sphériques et la génération de structures à partir d'elles [24, 25, 71, 72, 97, 101, 105, 133, 134]. Tous ces résultats constitueront une base de données, qui servira de socle pour la mise en place d'un système d'apprentissage.
# Bibliographie

## Conclusion générale et perspectives

[17] R. CLIFT, J. R. GRACE et M. E. WEBER. *Bubbles, Drops, and Particles*. Academic Press, 1978 (cf. p. 126, 146).

[24] S.-H. DONG. "Comment on 'A single-sum expression for the overlap integral of two associated Legendre polynomials'". In: *J. Phys. A : Math. Gen.* 35 (2002), p. 4187-4188 (cf. p. 146).

[25] S.-H. DONG et R. LEMUS. "The Overlap Integral of Three Associated Legendre Polynomials". In: *Appl. Math. Lett.* 15 (2002), p. 541-546 (cf. p. 146).

[30] D. ENRIGHT, F. LOSASSO et R. FEDKIW. "A Fast and Accurate Semi-Lagrangian Particle Level Set Method". In: *Computers & Structures* 83.6 (fév. 2005), p. 479-490 (cf. p. 50, 63, 64, 145).

[51] S. HYSING, S. TUREK, D. KUZMIN, N. PAROLINI, E. BURMAN, S. GANESAN et L. TOBISKA. "Quantitative benchmark computations of two-dimensional bubble dynamics". In: *Int. J. Numer. Methods Fluids* 60.11 (2009), p. 1259-1288 (cf. p. 138, 141, 146).

[65] R. J. LEVEQUE. "High-resolution conservative algorithms for advection in incompressible flows". In: *SIAM J. Numer. Anal.* 33.2 (avr. 1996), p. 627-665 (cf. p. 50, 63, 145).

[71] H. A. MAVROMATIS. "A single-sum expression for the overlap of two associated Legendre polynomials". In: *J. Phys. A : Math. Gen.* 32 (1999), p. 2601-2603 (cf. p. 146).

[72] H. A. MAVROMATIS et R. S. ALASSAR. "A Generalized Formula for the Integral of Three Associated Legendre Polynomials". In: *Appl. Math. Lett.* 12 (1999), p. 101-105 (cf. p. 146).

[97] G. RACAH. "Theory of Complex Spectra. II". In: *Phys. Rev.* 62 (nov. 1942), p. 438-462 (cf. p. 146).

[101] M. SATO. "General formula of the Racah coefficient". In: *Progress of Theoretical Physics* 13.4 (avr. 1955) (cf. p. 146).

[105] D. SÉBILLEAU. "On the computation of the integrated products of three spherical harmonics". In: *J. Phys. A : Math. Gen.* 31.34 (août 1998), p. 7157-7168 (cf. p. 146).

[133] Y.-L. XU. "Fast evaluation of the Gaunt coefficients". In: *Math. Comput.* 65.216 (oct. 1996), p. 1601-1612 (cf. p. 146).

[134] Y.-L. XU. "Fast evaluation of Gaunt coefficients : recursive approach". In: *J. Comput. Appl. Math.* 85 (1997), p. 53-65 (cf. p. 146).
# Annexe A

## La méthode de Cardan

La méthode de Cardan est une méthode de résolution d'équations polynomiales de degré 3. Ainsi, nous voulons résoudre l'équation suivante :

$$P(x) = ax^3 + bx^2 + cx + d = 0$$

Pour cela, $P(x)$ est transformé par un changement de variable en :

$$P(z) = z^3 + pz + q \text{ avec } x = z - \frac{b}{3a}$$

En effet :

$$\begin{aligned} P(z) &= a \left(z - \frac{b}{3a}\right)^3 + b \left(z - \frac{b}{3a}\right)^3 + c \left(z - \frac{b}{3a}\right) + d \\ &= az^3 + \left(c - \frac{b^2}{3a}\right) z - \frac{b^3}{27a^2} + \frac{b^3}{9a^2} - \frac{bc}{3a^2} + d \end{aligned}$$

pour $a \neq 0$ et en regroupant les monômes en $z$, nous avons :

$$P(z) = 0 \Leftrightarrow P(z) = z^3 + \left(\frac{c}{a} - \frac{b^2}{3a^2}\right)z - \frac{b}{27a}\left(\frac{2b^2}{a^2} - \frac{9c}{a}\right) + \frac{d}{a} = 0$$

Nous posons ainsi $p$ et $q$ tels que :

$$\begin{aligned} p &= \frac{c}{a} - \frac{b^2}{3a^2} \\ q &= \frac{b}{27a} \left( \frac{2b^2}{a^2} - \frac{9c}{a} \right) + \frac{d}{a} \end{aligned}$$

Nous avons alors $\Delta$, le discriminant de cette équation de la manière suivante :

$$\Delta = -(4p^3 + 27q^2)$$

Selon le signe de $\Delta$, nous nous trouvons dans trois cas possibles :

### 1. Discriminant $\Delta$ négatif

L'équation $P(z) = 0$ admet une solution réelle et deux complexes conjugués. En posant $u = \left[ \left(-q + \sqrt{-\Delta/27}\right) / 2 \right]^{1/3}$ et $v = \left[ \left(-q - \sqrt{-\Delta/27}\right) / 2 \right]^{1/3}$, la solution réelle $z_1 = u + v$, d'où $x_1 = z_1 - b/(3a)$. Quant aux racines complexes, nous les obtenons ainsi :

$$\begin{cases} z_2 = ju + j\bar{v} \\ z_3 = j^2u + \bar{j}^2v \end{cases} \text{ avec } j = e^{i2\pi/3}$$
2. Discriminant $\Delta$ positif

Dans ce cas, nous avons trois solutions réelles. La forme trigonométrique de ces solutions est la suivante :

$$z_k = 2\sqrt{\frac{-p}{3}}\cos\left(\frac{1}{3}arccos\left(\frac{-q}{2}\sqrt{\frac{27}{-p^3}}\right) + \frac{2k\pi}{3}\right), \text{ avec } k \in 0, 1, 2$$

Et nous retrouvons $x_k$ avec la relation $x_k = z_k - b/(3a)$.

3. Discriminant $\Delta$ nul

Si $\Delta = 0$, l'équation admet deux solutions réelles, une simple et une double.

$$\begin{cases} z_1 = 3q/p \\ z_2 = z_3 = -3q/2p \end{cases}$$

Dans le cas où $p$ et $q$ sont nuls, l'équation admet une solution triple triviale.
# Annexe B

## Calcul de la courbure par la méthode LBO - Cas d'un sommet à 3 triangles

Afin d'éclaircir notre lanterne, nous exposons notre démarche pour effectuer les calculs analytiques de l'opérateur de Laplace-Beltrami. Soit la sphère $\Omega$ de rayon $R$, centrée au point $O$ de coordonnées $(0, 0, 0)$. Sa surface est notée $\partial\Omega$. L'opérateur de Laplace-Beltrami est calculé au point $x_0$. Le stencil est d'abord plan, régulier, avec $[x_0\tilde{x}_i] = L, i = 1, 2, 3$ et $x_0 \in \partial\Omega$. Les points $\tilde{x}_i, i = 1, 2, 3$ sont projetés sur $\partial\Omega$, obtenant ainsi les points $x_i, i = 1, 2, 3$.

La figure suivante résume la procédure :

<figure><img src="image_105.png" /><figcaption>FIGURE B.1 – Vue de profil du stencil 3 et de sa projection sur la sphère de rayon $R$</figcaption></figure>

$P$ est le projeté de $x_i$ sur l'axe $(Ox_0)$. Nous établissons les relations entre les différents paramètres $L, r = [x_iP], h = [x_0P]$.

$$ \frac{L}{R} = \tan(\beta) \quad \text{or} \quad \tan(\beta) = \frac{2 \tan(\frac{\beta}{2})}{1 - \tan^2(\frac{\beta}{2})} \\ \text{et} \quad \tan(\frac{\beta}{2}) = \frac{d}{2\sqrt{R^2 - \frac{d^2}{4}}} $$

En remplaçant et après simplification, nous obtenons :

$$ L = \frac{d\sqrt{1 - \frac{d^2}{4R^2}}}{1 - \frac{d^2}{2R^2}} $$
Par le théorème de Thalès, nous avons :

$$ \frac{R}{\sqrt{L^2 + R^2}} = \frac{R-h}{R} $$

En utilisant l'expression de L, et après simplification, nous aboutissons à :

$$ h = \frac{d^2}{2R} $$

Pour finir, avec le théorème de Pythagore, nous trouvons r :

$$ r^2 + h^2 = d^2 \Leftrightarrow r^2 = d^2 - h^2 = d^2 - \frac{d^4}{4R^2} $$

Ainsi, l'expression de r est :

$$ r = d\sqrt{1 - \frac{d^2}{4R^2}} $$

Les points $(x_i)_{i \in \mathbb{N}}$ ont pour coordonnées :

$$ x_0 = \begin{cases} 0 \\ 0 \\ R \end{cases} \text{ et pour } i \in [1, \mathcal{N}_{x_0}] \ , \quad x_i = \begin{cases} r \cos\left(\frac{2\pi(i-1)}{\mathcal{N}_{x_0}}\right) \\ r \sin\left(\frac{2\pi(i-1)}{\mathcal{N}_{x_0}}\right) \\ R-h \end{cases} $$

Nous rappelons l'expression de l'opérateur de Laplace-Beltrami tel que décrite dans [76] :

$$ \mathbf{K}(\mathbf{x}_i) = \frac{1}{2\mathcal{A}_{\text{mixed}}} \sum_{j=1}^{\mathcal{N}_i} (\cot \alpha_{ij} + \cot \beta_{ij}) (\mathbf{x}_i - \mathbf{x}_j) $$

Ayant fait l'hypothèse d'un stencil régulier, nous avons que :

$$ \cot \alpha_{ij} = \cot \beta_{ij} = \cot \alpha $$

D'où dans le cas de stencils réguliers, l'opérateur de Laplace-Beltrami se réécrit :

$$ \mathbf{K}(\mathbf{x}_i) = \frac{\cot \alpha}{\mathcal{A}_{\text{mixed}}} \sum_{j=1}^{\mathcal{N}_i} (\mathbf{x}_i - \mathbf{x}_j) = \frac{\cot \alpha}{\mathcal{A}_{\text{mixed}}} \sum_{j=1}^{\mathcal{N}_i} \mathbf{u}_{ij} $$

Étant calculé au point $x_0$, nous avons :

$$ \mathbf{u}_{0i} = x_0 - x_i = \begin{cases} -r \cos\left(\frac{2\pi(i-1)}{\mathcal{N}_{x_0}}\right) \\ -r \sin\left(\frac{2\pi(i-1)}{\mathcal{N}_{x_0}}\right) \\ h \end{cases} $$

Nous pouvons remarquer, qu'à un facteur $-r$ près, les deux premières coordonnées de $\mathbf{u}_{0i}$ constituent les parties réelles et imaginaires des racines $n$-ièmes de l'unité (ici avec $n = \mathcal{N}_{x_0}$). Or :

$$ \sum_{k=1}^{\mathcal{N}_{x_0}} \exp\left(i \frac{2\pi(k-1)}{\mathcal{N}_{x_0}}\right) = 0 $$

Ainsi, nous pouvons simplifier l'expression de $\mathbf{K}(\mathbf{x}_0)$ :

$$ \mathbf{K}(\mathbf{x}_0) = \frac{\cot \alpha}{\mathcal{A}_{\text{mixed}}} \mathbf{u}_{\mathcal{N}_{x_0}} \quad \text{avec} \quad \mathbf{u}_{\mathcal{N}_{x_0}} = \begin{cases} 0 \\ 0 \\ \mathcal{N}_{x_0}h \end{cases} $$
<figure><img src="image_106.png" /><figcaption>FIGURE B.2 – Aire utilisée pour construire $\mathcal{A}_{mixed}$ dans le cas d'un triangle obtus.</figcaption></figure>

Reste plus qu’à expliciter cot α et $\mathcal{A}_{mixed}$. Commençons par $\mathcal{A}_{mixed}$. Dans le cas du stencil 3, elle s’exprime comme suit :

$$ \mathcal{A}_{mixed} = \frac{3}{2} A_T \quad \text{avec} \quad A_T = \text{aire du triangle } (x_0 x_1 x_2) $$

$$ \mathbf{a} = x_1 \vec{x}_2 = \begin{Bmatrix} -\frac{3}{2}r \\ r \frac{\sqrt{3}}{2} \\ 0 \end{Bmatrix}, \quad \mathbf{b} = x_1 \vec{x}_0 = \begin{Bmatrix} -r \\ 0 \\ h \end{Bmatrix} $$

Nous calculons $A_T$ en passant par le produit vectoriel de $\mathbf{a}$ et $\mathbf{b}$ :

$$ 4A_T^2 = \|\mathbf{a} \times \mathbf{b}\|^2 $$

$$ \mathbf{a} \times \mathbf{b} = \begin{Bmatrix} \frac{rh\sqrt{3}}{2} \\ \frac{3rh}{2} \\ \frac{r^2\sqrt{3}}{2} \end{Bmatrix} $$

$$ \begin{aligned} \text{D'où } 4A_T^2 &= 3(rh)^2 + \frac{3}{4}r^2 \\ &= \frac{3}{4}r^2(r^2 + 4h^2) \\ &= \frac{3}{4}d^2\left(1 - \frac{d^2}{4R^2}\right)\left[d^2\left(1 - \frac{d^2}{4R^2}\right) + \frac{d^4}{R^2}\right] \end{aligned} $$

$$ 4A_T^2 = \frac{3}{4}d^4\left(1 - \frac{d^2}{4R^2}\right)\left(1 + \frac{3d^2}{4R^2}\right) $$

$$ \text{Donc } A_T = \frac{d^2 \sqrt{3}}{4} \sqrt{\left(1 - \frac{d^2}{4R^2}\right) \left(1 + \frac{3d^2}{4R^2}\right)} $$

Nous en tirons donc $\mathcal{A}_{mixed}$ dans le cas d'un stencil 3 :

$$ \mathcal{A}_{mixed} = \frac{3\sqrt{3}}{8} d^2 \sqrt{\left(1 - \frac{d^2}{4R^2}\right) \left(1 + \frac{3d^2}{4R^2}\right)} $$
Ensuite, nous calculons cot α avec la formule :

$$ \cot \alpha = \frac{\mathbf{a} \cdot \mathbf{b}}{\sqrt{\|\mathbf{a}\|^2 \|\mathbf{b}\|^2 - (\mathbf{a} \cdot \mathbf{b})^2}} $$

$$ \|\mathbf{a}\|^2 = 3r^2 = 3d^2 \left(1 - \frac{d^2}{4R^2}\right) $$

$$ \|\mathbf{b}\|^2 = r^2 + h^2 = d^2 $$

$$ \mathbf{a} \cdot \mathbf{b} = \frac{3}{2}r^2 = \frac{3}{2}d^2 \left(1 - \frac{d^2}{4R^2}\right) $$

$$ \text{D'où } \cot \alpha = \frac{\frac{3}{2} d^2 \left(1 - \frac{d^2}{4R^2}\right)}{\sqrt{3d^4 \left(1 - \frac{d^2}{4R^2}\right) - \frac{9}{4} d^4 \left(1 - \frac{d^2}{4R^2}\right)^2}} \\ = \frac{\sqrt{1 - \frac{d^2}{4R^2}}}{\sqrt{\frac{4}{3} - \left(1 - \frac{d^2}{4R^2}\right)}} $$

$$ \text{Donc } \cot \alpha = \sqrt{3} \frac{\sqrt{1 - \frac{d^2}{4R^2}}}{\sqrt{1 + \frac{3d^2}{4R^2}}} $$

Nous réécrivons donc le rapport $\frac{\cot \alpha}{A_{mixed}}$ :

$$ \frac{\cot \alpha}{A_{mixed}} = \frac{8}{3d^2 \left(1 + \frac{3d^2}{4R^2}\right)} $$

Nous pouvons ainsi explicitier l'opérateur de Laplace-Beltrami pour un stencil régulier de 3 :

$$ \mathbf{K}(\mathbf{x}_0) = \frac{8}{3d^2 \left(1 + \frac{3d^2}{4R^2}\right)} \mathbf{u}_{\mathcal{N}_{x_0}} $$

$$ \mathcal{N}_{x_0} = 3 $$

$$ \mathbf{K}(\mathbf{x}_0) = \frac{8h}{d^2 \left(1 + \frac{3d^2}{4R^2}\right)} \mathbf{n} \text{, avec } \mathbf{n} = \begin{cases} 0 \\ 0 \\ 1 \end{cases} \text{, et } h = \frac{d^2}{2R} $$

$$ \text{Ainsi, } \mathbf{K}(\mathbf{x}_0) = \frac{\frac{4}{R}}{1 + \frac{3d^2}{4R^2}} \mathbf{n} $$

Nous en déduisons la courbure $\kappa_H = \frac{1}{2}(\frac{1}{R_1} + \frac{1}{R_2})$ :

$$ \begin{aligned} \kappa_H &= \frac{1}{2} \| \mathbf{K}(\mathbf{x}_0) \| \\ &= \frac{\frac{2}{R}}{1 + \frac{3d^2}{4R^2}} \end{aligned} $$

Pour une sphère de rayon $R$, la courbure est $\kappa_H^{exact} = \frac{1}{R}$. Le calcul analytique avec la formule de [76] donne une courbure $\kappa_H^{an.} = \frac{2}{1 + \frac{3d^2}{4R^2}}$. Ce que nous retrouvons avec les simulations numériques :
<figure><img src="image_107.png" /><figcaption>Courbure en fonction de la taille relative du stencil d/R</figcaption></figure>
# Annexe C

## Complément sur les comparaisons de méthodes

### C.1 Errors for 20%-disturbed meshes

As described in the body of the present article, when the initial mesh was made up of equilateral triangles, the remeshing process has destroyed the regular structure of the discrete surface. As a consequence, the measured errors were much larger than those recorded for the stationary surface, especially for the curvature. Therefore, introducing random disturbances on the initial location of the vertices does not significantly modify, neither the behaviour, nor the level of the errors. This remark is valid whatever the *Q*-value, as well for the exact advection (Eq. 3.25) :

— see Fig. C.1 versus Fig. 3.20 for $Q < 0$,

— see Fig. C.4 versus Fig. 3.24 for $Q > 0$

as for the transport numerical scheme :

— see Fig. C.2 versus Fig. 3.21 for $Q < 0$,

— see Fig. C.5 versus Fig. 3.25 for $Q > 0$.

and for the CPU-cost and efficiency :

— see Fig. C.3 versus Fig. 3.22 for $Q < 0$,

— see Fig. C.6 versus Fig. 3.26 for $Q > 0$.

<figure><img src="image_108.png" /><figcaption>FIGURE C.1 – Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$ and $\mathbf{n}$, $p = 20\%$ with exact transport - $Q < 0$</figcaption></figure>
<figure><img src="image_109.png" /><figcaption>FIGURE C.2 – Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$ and $\mathbf{n}$, $p = 20\%$ with Lagrangian transport - $Q < 0$</figcaption></figure>

<figure><img src="image_110.png" /><figcaption>FIGURE C.3 – Relative errors on $\kappa$ (a) and $\mathbf{n}$ (b) according to CPU time in $L_2$ (continuous line) and $L_\infty$ (dash line) norms, $p = 20\%$ - $Q < 0$ and numerical transport.</figcaption></figure>

<figure><img src="image_111.png" /><figcaption>FIGURE C.4 – Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$ and $\mathbf{n}$, $p = 20\%$ with exact transport - $Q > 0$</figcaption></figure>
<figure><img src="image_112.png" /><figcaption>FIGURE C.5 – Relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$ and $\mathbf{n}$, $p = 20\%$ with Lagrangian transport - $Q > 0$</figcaption></figure>

<figure><img src="image_113.png" /><figcaption>FIGURE C.6 – Relative errors on $\kappa$ (a) and $\mathbf{n}$ (b) according to CPU time in $L_2$ (continuous line) and $L_\infty$ (dash line) norms, $p = 20\%$ - $Q > 0$ and numerical transport.</figcaption></figure>
# Annexe D

## Complément sur les méthodes de Reconstruction de Surface

### D.1 Peculiar stencil for the calculation of the SR-$\mathcal{N}_1$-$\mathbf{n}_m$-0 method

The calculation of the curvature with the Surface Reconstruction SR-$\mathcal{N}_1$-$\mathbf{n}_m$-0 [26] may lead to very bad approximations, as shown in Fig. 3.24a for the $L_\infty$-norm. To better understand the reason of such a result, the vertices at the root of these large errors have to be identified. To this end, the surface shape is subdivided into three sub-domains $\mathcal{D}_i$, characterized by very different principal curvature values (Fig. D.1) : $\mathcal{D}_1$ corresponds to the nearly flat annular part, $\mathcal{D}_2$ is the

<figure><img src="image_114.png" /><figcaption>FIGURE D.1 – Subdivision of the discrete mesh in three zones.</figcaption></figure>

hemispherical surface and $\mathcal{D}_3$ is referred to as the narrow strip which connects the two previous domains.

Figure D.2 summarizes the relative errors on the curvature calculation as a function of the mesh size $d \times \max_{x_i \in \mathcal{D}_j}(\kappa(x_i))$, for $1 \le j \le 3$ and the usual numerical methods. The comparison of these curves with Fig. 3.24a indicates that the increase in the relative error is due to vertices located in the sub-domain $\mathcal{D}_3$ where the two principal curvatures are large and have opposite signs (Fig. D.2c).

A thoroughly examination of the vertices with large errors shows that the nodes used to define the fitted surface $Z = f(X, Y)$ (Eq. 3.14) are distributed in a peculiar way. Due to the refinement procedure, which consists in adding a vertex at the middle of an edge, three nodes are compulsorily aligned so that their projections onto the local plane $Z = 0$ belong to an unique
<figure><img src="image_115.png" /><figcaption>FIGURE D.2 – Decomposition per zone of the relative errors in $L_2$ (continuous line) and $L_\infty$ (dash line) norms on $\kappa$, $p = 0\%$ with exact transport - $Q > 0$</figcaption></figure>

straight line noted $\Delta$. Furthermore, this situation often produces a first neighbourhood $\mathcal{N}_1$ with 4 triangles sharing the same vertex (Fig. D.3), what does not provide enough points to determine in a unique way the fitted surface. Thus, an arbitrary additional vertex must also be taken into

<figure><img src="image_116.png" /><figcaption>FIGURE D.3 – Stencil of the marker where the curvature and normal vector are not well evaluated.</figcaption></figure>

account in the second neighbourhood $\mathcal{N}_2 \setminus \mathcal{N}_1$ to calculate the function $Z = f(X, Y)$. Based on our meshes, the bad curvature evaluations seem resulting from configurations where the fifth point, projected over the $Z = 0$ plane, is very close to the straight line $\Delta$ (four points nearly aligned). To confirm this intuition, the four nodes belonging to the first neighbourhood $\mathcal{N}_1$ of a vertex labelled $M_0$, located in the $D_2$ domain, have been extracted from the mesh, whereas the fifth node $M_5$ is able to move all around $M_0$, while remaining on the analytical surface. Let the angle $\theta = M_5'M_0'M_1'$, with $M_i'$ the projection of $M_i$ onto the plane $Z = 0$ and the length $M_5'M_0'$
kept constant, then the relative error on the curvature is maximal when $\theta = 0$ [$\pi$] (Fig. D.4). The

<figure><img src="image_117.png" /><figcaption>FIGURE D.4 – Relative error on $\kappa$ in function of the angle of rotation $\theta$ of the markers $M_i \in (\mathcal{N}_2 \setminus \mathcal{N}_1)$, for the variant "Glimm - $\mathcal{N}_{1-1}$". The vertical dash lines are the principal directions.</figcaption></figure>

level difference measured on the curvature error may reach three decades between the average value of order 0.01 and the maximum $\approx 10$. Before ending this section, it is worth noticing that such a nearly alignment of vertices happens also in the meshes of the other domains $\mathcal{D}_1$ and $\mathcal{D}_2$ (Fig. D.1). However, in those cases, the two principal curvatures are almost identical and therefore, the loss of accuracy in one direction is probably counterbalanced by the other nodes taking part in the surface fitting. On the other hand, when the principal curvatures are significantly different, like in $\mathcal{D}_3$ domain, the loss of accuracy cannot be compensated by the other vertices.

To conclude this section, the use of the Surface Reconstruction SR-$\mathcal{N}_1$-$\mathbf{n}_m$-0 by Du et al. [26] must be used carefully when the first neighbourhood does not have enough degree of freedom to calculate the fitted surface. One way to overcome this issue is to always consider the two first neighbourhoods $\mathcal{N}_2$ to avoid degenerated situations as described here-above.
# Annexe E

## Effet de la pondération dans les méthodes SR

Nous profitons également de cette étude numérique pour tester l'introduction d'une pondération lors de l'établissement des méthodes de Reconstruction de surface. En effet, jusqu'ici, nous avons utilisé une pondération de 1 dans la méthode des moindres carrés servant à établir l'équation locale de la surface discrète. L'obtention de l'équation locale au marqueur m, de coordonnées $\bar{x}_m$, s'obtient en minimisant le résidu RMS évoqué dans l'équation (E.1).

$$ RMS(m) = \sum_{m_i \in \mathcal{N}_k^*(m)} \omega_i [f(X_i, Y_i) - Z_i]^2 \quad (E.1) $$

avec $\omega_i$ la pondération associée au marqueur $m_i$, de coordonnées locales $X_i, Y_i, Z_i$, f la fonction polynomiale définie par l'équation (3.14), page 78, et $k \in \{1, 2\}$ selon le choix du voisinage pour les méthodes SR-$\mathcal{N}_1^*$ ou du second voisinage pour les méthodes SR-$\mathcal{N}_2^*$. Concernant la pondération $\omega_i$, nous faisons le choix de l'exprimer sous la forme $1/(1+x^2)$, comme présentée dans l'équation (E.2).

$$ \omega_i = \frac{1}{1 + \left[ \frac{\|\bar{x}_i - \bar{x}_m\| \times N}{\min_{m_i \in \mathcal{N}_k^*(m)} \|\bar{x}_i - \bar{x}_m\|} \right]^2} \quad (E.2) $$

avec, $N \in \{0, 1, 2, 4, 8\}$, un entier permettant de jouer sur la localité de la reconstruction. De ce fait, plus l'on s'éloigne du marqueur central m, moins les marqueurs de son voisinage ont de l'importance dans l'établissement de l'équation locale de la surface.

### E.1 Analyse de l'influence de la pondération et des perturbations sur la distribution surfacique

Nous nous intéressons dans cette partie à la caractérisation de la densité surfacique de courbure sur une famille d'ellipsoïdes de volume fixé. Deux paramètres rentrent en ligne de compte dans notre étude, à savoir les perturbations du maillage de l'interface, d'amplitude maximale $p \in [0, 0.4]$, et l'introduction ou non d'une pondération $\omega_i$ lors de l'établissement de l'équation locale de la surface par la méthode des moindres carrés.
<figure><img src="image_118.png" /><figcaption>FIGURE E.1 – Evolution de la densité surfacique de courbure pour une famille d'ellipsoïdes de divers étirements.</figcaption></figure>

Distribution surfacique de courbure pour une famille d'ellipsoïdes de volume constant

Avant toute chose, nous présentons dans la figure E.1 la densité surfacique de courbure théorique pour une famille d'ellipsoïdes, construite selon les relations (4.12) et (4.13), avec $\alpha \in \{20\%, 40\%, 60\%, 80\%, 100\%, 120\%, 140\%\}$. Ces ellipsoïdes sont générées à partir d'une sphère, de maillage initial $r/d = 256$. Les courbures sont calculées de manière analytique. La plage de courbure présente sur chaque ellipsoïde est divisée en 20 classes de tailles identiques.

Comme nous pouvons l'observer, pour un faible étirement ($\alpha = 20\%$), la $df_S$ est très proche de la forme d'un pic. En effet, pour de faibles valeurs de $\alpha$, l'ellipsoïde est presque assimilable à une sphère. Cette dernière ne présente qu'une seule courbure, et sa distribution $df_S$ est un pic localisé à la valeur de sa courbure. Plus la valeur de $\alpha$ s'éloigne de 0, plus la distribution s'étale, avec l'apparition de zones de fortes courbures. Toutefois les zones de faible courbure, situées au centre de l'ellipsoïde, continuent d'occuper la majeure partie de la surface. Les régions de forte courbure sont focalisées aux extrémités allongées, couvrant par conséquent une surface bien moindre (de l'ordre de $10^{-3}$ pour les fortes courbures sur l'ellipsoïde de paramètre $\alpha = 140\%$).

Restriction des tests aux méthodes de reconstruction de surface

L'idée des tests menés dans cette partie est de voir la justesse avec laquelle nos différentes méthodes arrivent à retrouver ces distributions surfaciques de courbures. Pour cela, nous nous servons de la mesure introduite en (4.14) qui nous donnent l'écart en norme $L_2$ et $L_\infty$ entre nos méthodes et la solution analytique, et ceci pour différentes finesse de maillages de l'interface, représentées par la quantité $d_{moy} * \kappa_{max}$. Afin de justifier notre choix de nous restreindre qu'aux méthodes de reconstruction de surface SR-$\mathcal{N}_1^*$ et SR-$\mathcal{N}_2^*$, nous montrons, dans un premier temps, les résultats obtenus en se servant de la méthode LBO, avec différentes perturbations initiales du maillage (Fig. E.2).

Dans des configurations idéales, à savoir des ellipsoïdes proches de la sphère ($\alpha = 20\%$) et
aucune perturbation du maillage (p = 0), la méthode LBO présente un ordre 1 de convergence en norme $L_2$ et des erreurs relatives très faibles pour la résolution la plus fine ($10^{-3}$ pour $d_{moy} * \kappa_{max} = 8 \times 10^{-3}$, Fig. E.2a). Cette convergence en norme $L_2$ est perdue dès l'introduction de perturbations dans le maillage de l'interface, avec des erreurs relatives restant inférieures à $10^{-2}$ (Fig. E.2b et E.2c). Cependant, les erreurs en norme $L_\infty$ ne convergent pas, voire divergent avec les perturbations du maillage. De surcroît, plus on s'éloigne de la sphère (α prenant des valeurs grandes), plus le comportement des erreurs devient erratique, avec des erreurs relatives pouvant dépasser $10^1$.

Tous ces éléments nous confortent à mettre de côté les méthodes LBO. Reste à tester les méthodes de reconstruction de surface en mettant en lumière leur sensibilité aux perturbations p du maillage et l'utilisation de moindres carrés pondérés.

## Influence de la pondération $\omega_i$

Nous nous plaçons dans la configuration où aucune perturbation n'est introduite sur le maillage de la sphère initiale servant à générer la famille d'ellipsoïdes. Le poids $\omega_i$ (E.2) introduit dans l'obtention de l'équation locale de la surface par minimisation du résidu (E.1) prend en compte l'éloignement des points par rapport au marqueur central, et donc leur influence de l'établissement de l'équation locale. Le facteur $N \in \{0, 1, 2, 4, 8\}$ inséré dans l'expression de $\omega_i$ permet d'avoir différents valeurs de poids pour le même poids, tout en gardant la forme générale de la pondération choisie en $1/(1+x^2)$. Nous présentons dans la figure (cite) les erreurs en norme $L_2$ et $L_\infty$ sur les distributions surfaciques de courbure pour les valeurs de $N \in \{0, 1, 8\}$, les autres étant identiques. Les erreurs sont calculées selon (4.14).

Contrairement à ce que nous avions avec les méthodes LBO, les méthodes SR (ici SR-$\mathcal{N}_{1-m}-0$) convergent. Quelle que soit la pondération utilisée, nous observons une convergence de 2 en norme $L_2$ sur des structures proches de la sphère, à $\alpha = 20\%$. Cet ordre diminue lorsque la structure s'éloigne de plus en plus de la sphère, pour se stabiliser à 1 sur la gamme d'étirements employés. La norme $L_\infty$ converge également à l'ordre 1 sauf pour les étirements correspondant à $\alpha = 60\%$ et $\alpha = 80\%$, où la convergence est de 0.5. A résolution de maillage la plus fine (de l'ordre de $d_{moy} * \kappa_{max} = 10^{-2}$), aucune différence notable ne ressort entre la version sans pondération $N=0$ (Fig. E.3a) et les versions avec pondération $N=1, N=8$ (Fig. E.3b-E.3c). Cela s'explique par le fait que la distance entre les marqueurs de l'interface est suffisamment petite pour capter toutes les courbures de l'interface. Ainsi, l'amélioration apportée par la prise en compte de la localité des marqueurs devient négligeable.

A contrario, l'introduction de la pondération améliore les résultats sur les maillages grossiers. Dans notre cas, en particulier pour les étirements $\alpha = 120\%$ et $\alpha = 140\%$ à $d_{moy} * \kappa_{max} = 1$ (Fig. E.3b-E.3c). Bien que ces erreurs restent relativement élevées, il n'en disconvient pas que la pondération les réduit comparativement au cas $N=0$. En effet, la structure étant maillé grossièrement, les modifications topologiques entre deux marqueurs risquent de ne pas être captées. Donner un poids qui décroît avec la distance mitige ce manque d'information inhérent au maillage. Toutefois, nous ne remarquons pas d'écarts significatifs entre $N=1$ et $N=8$.

Ainsi, au vu des résultats de cette partie mettant en avant l'amélioration des erreurs pour des maillages grossiers par l'introduction de pondération non uniforme, nous présentons l'impact des perturbations p du maillage initial que pour $N=1$. Il est à noter que les calculs ont également été fait pour $N \in \{0, 1, 2, 4, 8\}$.

## Influence de la perturbation p

Un perturbation d'amplitude maximale p du maillage initial de la sphère génératrice de la famille d'ellipsoïdes est introduite. La figure E.4 synthétise les résultats pour $p \in \{0\%, 20\%, 40\%\}$

A p = 20% de perturbation (Fig. E.4b), la convergence en norme $L_2$ de la distribution $df_S$ commence à s'infléchir pour des structures sphéroïdales ($\alpha = 20\%$ et $\alpha = 40\%$) mais reste entre
1.5 et 2. Les structures ellipsoïdales présentent quant à elles sur leur maillage une convergence en norme $L_2$ et $L_\infty$ de 1. A $p = 40\%$ de perturbation (Fig. E.4c), toutes les structures sont à la même enseigne et continuent de converger à l'ordre 1 en norme $L_2$ et $L_\infty$. Bien que les erreurs soient pour grandes pour $\alpha = 20\%$ et $\alpha = 40\%$ entre le cas $p = 0\%$ et $p = 40\%$ (Fig. E.4a et E.4c), celles pour $\alpha \ge 60\%$ restent inchangées.
<figure><img src="image_119.png" /><figcaption>FIGURE E.2 – Erreurs relatives en norme L2 (trait continu) et L∞ (tiret) sur la densité surfacique de courbure d'une famille d'ellipsoïdes, avec la méthode LBO.</figcaption></figure>
<figure><img src="image_120.png" /><figcaption>FIGURE E.3 – Influence de la pondération - Erreurs relatives en norme L2 (trait continu) et L∞ (tiret) sur la densité surfacique de courbure d'une famille d'ellipsoïdes, à perturbation p = 0%, avec la méthode SR-N1-nm-0.</figcaption></figure>
<figure><img src="image_121.png" /><figcaption>(a) Perturbation $p = 0\%$</figcaption></figure>

<figure><img src="image_122.png" /><figcaption>(b) Perturbation $p = 20\%$</figcaption></figure>

<figure><img src="image_123.png" /><figcaption>(c) Perturbation $p = 40\%$</figcaption></figure>

FIGURE E.4 – Influence de la perturbation - Erreurs relatives en norme L2 (trait continu) et L∞ (tiret) sur la densité surfacique de courbure d'une famille d'ellipsoïdes, à paramètre $N = 1$, avec la méthode SR-$\mathcal{N}_1$-$\mathbf{n}_m$-0.
# Annexe F

## Calcul exact de la distribution surfacique $df_S$ pour des ellipsoïdes

Cette section présente les différentes étapes pour le calcul exact de la distribution surfacique $df_S$ pour des ellipsoïdes d'étirement $\alpha$. La construction de ces ellipsoïdes est mentionné dans la section 4.2.2, page 113. Pour des facilités de calcul, l'étirement est effectué ici dans la direction $z$, avec pour conséquence le nouveau système de coordonnées (F.1) :

$$ \bar{x}/r_0 = \begin{cases} x = b \sin(\theta) \cos(\varphi)/r_0 \\ y = b \sin(\theta) \sin(\varphi)/r_0 \\ z = a \cos(\theta)/r_0 \end{cases} \qquad (F.1) $$

avec $a$ et $b$ le demi-grand axe et le demi-petit axe de l'ellipsoïde vérifiant la relation (F.2)

$$ a/r_0 = (1 + \alpha) \quad , \quad b/r_0 = 1/(1 + \alpha)^{1/2} \qquad (F.2) $$

Une première étape est la détermination de l'élément différentiel surfacique $dA$ permettant de calculer des surfaces ou portions de surfaces sur une ellipsoïde.

### Élément différentiel surfacique $dA$

L'équation paramétrique de l'ellipsoïde donnée par la relation (F.1) permet grâce aux dérivés partielles en $\theta$ et $\varphi$ de $\bar{x}$ d'aboutir à l'élément différentiel surfacique $dA$. Nous notons :

$$ \bar{x}_{\theta} = \frac{\partial \bar{x}}{\partial \theta}, \quad \bar{x}_{\varphi} = \frac{\partial \bar{x}}{\partial \varphi} $$

Ensuite, nous calculons les coefficients $E$, $F$ et $G$ qui sont les coefficients de la 1ère forme fondamentale associée à l'équation paramétrique (F.1).

$$ E = \bar{x}_{\theta} \cdot \bar{x}_{\theta}, \quad F = \bar{x}_{\theta} \cdot \bar{x}_{\varphi}, \quad G = \bar{x}_{\varphi} \cdot \bar{x}_{\varphi} $$

L'élément différentiel surfacique $dA$ est défini par l'équation (F.3).

$$ dA = \sqrt{EG - F^2} d\theta d\varphi \qquad (F.3) $$

Ce qui nous donne dans ce cas de figure la relation (F.4).
$$
dA = b|\sin\theta|\sqrt{a^2\sin^2\theta + b^2\cos^2\theta}d\theta d\varphi \quad (F.4)
$$

Ainsi pour calculer l'aire A(θ₁, θ₂, φ₁, φ₂) d'une surface comprise entre [θ₁, θ₂] et [φ₁, φ₂], il suffit
d'intégrer dA sur ce domaine, comme mentionné dans la relation (F.5).

$$
A(\theta_1, \theta_2, \varphi_1, \varphi_2) = \int_{\theta_1}^{\theta_2} \int_{\varphi_1}^{\varphi_2} b |\sin \theta| \sqrt{a^2 \sin^2 \theta + b^2 \cos^2 \theta} d\theta d\varphi \quad (\text{F.5})
$$

Si cette aire peut être évaluée avec un logiciel de calcul formel, il n’en demeure pas moins
qu’il faut au préalable déterminer le domaine $[\theta_1, \theta_2] \times [\varphi_1, \varphi_2]$ pour le calcul de la distribution
df<sub>S</sub>.

Détermination du domaine $[\theta_1, \theta_2] \times [\varphi_1, \varphi_2]$

Dans le calcul de la distribution surfacique df<sub>S</sub>, le domaine [θ₁, θ₂] × [φ₁, φ₂] sur lequel une classe ̂κ de courbure est définie est connu par les courbures κ au bord de ce domaine. Comme les ellipsoïdes sur lesquels nous travaillons sont de révolution autour de l'axe z, l'intervalle [φ₁, φ₂] = [0, 2π]. Seuls θ₁ et θ₂ servent à délimiter la classe ̂κ de courbure. Ainsi le problème est le suivant : connaissant les courbures κ(θ₁) et κ(θ₂), quelles sont les valeurs de θ₁ et θ₂ ?

Dans un premier temps, voyons comment est calculée la courbure κ. Avec les conventions
utilisées dans la relation (F.1), la courbure moyenne κ et la courbure de Gauss κ_G sont évaluées
avec le système d'équation (F.6).

$$
\left\{
\begin{aligned}
R &= (x^2 + y^2 + z^2)^{1/2} \\
H &= (x^2/b^4 + y^2/b^4 + z^2/a^4)^{-1/2} \\
\kappa &= H^3 (2b^2 + a^2 - R^2) / (2b^4 a^2) \\
\kappa_G &= H^4 / (b^4 a^2)
\end{aligned}
\right.
\qquad (\text{F.6})
$$

A partir de la courbure de Gauss $\kappa_G$

Soit $\kappa_G^0$ une courbure de Gauss donnée. Nous cherchons $\theta \in [0, \pi]$ tel que $\kappa_G(\theta) = \kappa_G^0$, avec :

$$
\kappa_G(\theta) = \frac{a^2}{[(a^2 - b^2) \sin^2 \theta + b^2]^2} \quad (F.7)
$$

Les solutions qui nous intéressent sont celles comprises dans [0, π], plus précisément [0, π/2] car
celles dans l'intervalle [π/2, π] s'obtiennent par symétrie. La solution de ce problème est donnée
la relation (F.8).

$$
\theta = \arcsin \left( \left[ \frac{\sqrt{a^2/\kappa_G^0 - b^2}}{a^2 - b^2} \right]^{1/2} \right) \qquad (F.8)
$$
A partir de la courbure moyenne $\kappa$

Soit $\kappa^0$ une courbure moyenne donnée. Pour les mêmes raisons que précédemment, nous cherchons $\theta \in [0, \pi/2]$ tel que $\kappa(\theta) = \kappa^0$, avec :

$$ \kappa(\theta) = \frac{a[(a^2 - b^2) \sin^2 \theta + 2b^2]}{2b[(a^2 - b^2) \sin^2 \theta + b^2]^{3/2}} \qquad (F.9) $$

Rappelons que $ab^2 = r_0^3$ (voir Eq. F.2). Cela revient à résoudre l'équation (F.10).

$$ -2b\kappa^0 \left[ (a^2 - b^2) \sin^2 \theta + b^2 \right]^{3/2} + (a^3 - r_0^3) \sin^2 \theta + 2r_0^3 = 0 \qquad (F.10) $$

Notons $X = \sin^2 \theta$, $P(X) = -2b\kappa^0 \left[ (a^2 - b^2)X + b^2 \right]^{3/2} + (a^3 - r_0^3)X + 2r_0^3$ et $X_0$ la solution de $P(X) = 0$. Alors la solution de (F.10) est :

$$ \theta = \arcsin(X_0^{1/2}) \qquad (F.11) $$
Bibliographie

[1] T. ABADIE, J. AUBIN et D. LEGENDRE. "On the combined effects of surface tension force calculation and interface advection on spurious currents within Volume of Fluid and Level Set frameworks". In : *J. Comput. Phys.* 297 (2015), p. 611-636 (cf. p. 34, 35).

[2] J. ADELSBERGER, S. GROSS, M. KLITZ, A. RÜTTGERS et M. GRIEBEL. "3D Incompressible two-phase flow benchmark computations for rising droplets". In : *11th World Congress on Computational Mechanics (WCCM XI), 5th European Conference on Computational Mechanics (ECCM V), 6th European Conference on Computational Fluid Dynamics (ECFD VI)*. Sous la dir. d'E. ÕNATE, J. OLIVER et A. HUERTA. 2013 (cf. p. 138, 141).

[3] H. T. AHN et M. SHASHKOV. "Multi-material interface reconstruction on generalized polyhedral meshes". In : *J. Comput. Phys.* 226 (2007), p. 2096-2132 (cf. p. 27).

[4] S. L. ANNA. "Droplets and Bubbles in Microfluidic Devices". In : *Annu. Rev. Fluid Mech.* 48 (2016), p. 285-309 (cf. p. 20, 21).

[5] M. ATAEI, M. BUSSMANN, V. SHAAYEGAN, F. COSTA, S. HAN et C. B. PARK. "NPLIC : A machine learning approach to piecewise linear interface construction". In : *Computers and Fluids* 223 (2021), p. 104950 (cf. p. 27).

[6] E. AULISA, S. M. and R. SCARDOVELLI et S. ZALESKI. "Interface reconstruction with least-squares fit and split advection in three-dimensional Cartesian geometry". In : *J. Comput. Phys.* 225 (2007), p. 2301-2319 (cf. p. 27).

[7] E. BABINSKY et P. SOJKA. "Modeling drop size distributions". In : *Prog. Energy Combust. Sci.* 28 (2002), p. 303-329 (cf. p. 20).

[8] S. BALACHANDAR et J. K. EATON. "Turbulent Dispersed Multiphase Flow". In : *Annu. Rev. Fluid Mech.* 42 (2010), p. 111-133 (cf. p. 20).

[9] D. BHAGA et M. E. WEBER. "Bubbles in viscous liquids : shapes, wakes and velocities". In : *J. Fluid Mech.* 105 (1981), p. 61-85 (cf. p. 126).

[10] J. BRACKBILL, D. KOTHE et C. ZEMACH. "A continuum method for modeling surface tension". In : *J. Comput. Phys.* 100.2 (1992), p. 335-354 (cf. p. 31, 32, 73).

[11] N. BREMOND et E. VILLERMAUX. "Atomization by jet impact". In : *J. Fluid Mech.* 549 (2006), p. 273-306 (cf. p. 20).

[12] C. E. BRENNEN. *Fundamentals of Multiphase Flows*. Cambridge University Press, 2005 (cf. p. 21).

[13] R. CANU, S. PUGGELLI, M. ESSADKI, B. DURET, T. MÉNARD, M. MASSOT, J. REVEILLON et F. DEMOULIN. "Where does the droplet size distribution come from ?" In : *Int. J. Multiphase Flow* 107 (2018), p. 230-245 (cf. p. 36, 72, 73).

[14] F. CAZALS et M. POUGET. "Estimating differential quantities using polynomial fitting of osculating jets". In : *Comput. Aided Geom. Des.* 22.2 (2005), p. 121-146 (cf. p. 97).
[15] S. CHANDRA et C. T. AVEDISIAN. "On the Collision of a Droplet with a Solid Surface". In: *Proceedings : Mathematical and Physical Sciences* 432.1884 (jan. 1991), p. 13-41 (cf. p. 21).

[16] P.-C. CHEN, W.-C. WANG, W. L. ROBERTS et T. FANG. "Spray and atomization of diesel fuel and its alternatives from a single-hole injector using a common rail fuel injection system". In: *Fuel* 103 (2013), p. 850-861 (cf. p. 21).

[17] R. CLIFT, J. R. GRACE et M. E. WEBER. *Bubbles, Drops, and Particles*. Academic Press, 1978 (cf. p. 126, 146).

[18] G. COMPÈRE, E. MARCHANDISE et J.-F. REMACLE. "Transient adaptivity applied to two-phase incompressible flows". In: *J. Comput. Phys.* 227 (2008), p. 1923-1942 (cf. p. 23).

[19] M. CUTFORTH, P. T. BARTON et N. NIKIFORAKIS. "An efficient moment-of-fluid interface tracking method". In: *Computers and Fluids* 224 (2021), p. 104964 (cf. p. 27).

[20] M. DAI et D. P. SCHMIDT. "Adaptive tetrahedral meshing in free-surface flow". In: *J. Comput. Phys.* 208 (2005), p. 228-252 (cf. p. 23).

[21] J. DELHAYE. "Jump conditions and entropy sources in two-phase systems. Local instant formulation". In: *Int. J. Multiphase Flow* 1.3 (1974), p. 395-409 (cf. p. 23, 72).

[22] F. DEMOULIN, P. BEAU, G. BLOKKEEL, A. MURA et R. BORGHI. "A new model for turbulent flows with large density fluctuations : Application to liquid atomization". In : *Atomization and Sprays* 17 (2007), p. 315-345 (cf. p. 73).

[23] M. DESBRUN, M. MEYER, P. SCHRODER et A. BARR. "Implicit fairing of irregular meshes using diffusion and curvature flow." In: *SIGGRAPH99*. 1999, p. 317-324 (cf. p. 75).

[24] S.-H. DONG. "Comment on 'A single-sum expression for the overlap integral of two associated Legendre polynomials'". In: *J. Phys. A : Math. Gen.* 35 (2002), p. 4187-4188 (cf. p. 146).

[25] S.-H. DONG et R. LEMUS. "The Overlap Integral of Three Associated Legendre Polynomials". In: *Appl. Math. Lett.* 15 (2002), p. 541-546 (cf. p. 146).

[26] J. DU, B. FIX, J. GLIMM, X. JIA, X. LI, Y. LI et L. WU. "A simple package for front tracking". In: *J. Comput. Phys.* 213 (2006), p. 613-628 (cf. p. 29, 33, 73, 78, 79, 83, 85, 90, 92, 107, 161, 163).

[27] L. DUCHEMIN, S. POPINET, C. JOSSERAND et S. ZALESKI. "Jet formation in bubbles bursting at a free surface". In: *Phys. Fluids* 14.9 (sept. 2002), p. 3000 (cf. p. 19).

[28] C. DUMOUCHEL. "On the experimental investigation on primary atomization of liquid streams". In: *Exp. Fluids* 45 (2008), p. 371-422 (cf. p. 20).

[29] V. DYADECHKO et M. SHASHKOV. "Moment-of-fluid interface reconstruction". In: *Tech. Rep. LA-UR-05-7571* (2005) (cf. p. 27).

[30] D. ENRIGHT, F. LOSASSO et R. FEDKIW. "A Fast and Accurate Semi-Lagrangian Particle Level Set Method". In: *Computers & Structures* 83.6 (fév. 2005), p. 479-490 (cf. p. 50, 63, 64, 145).

[31] M. ESSADKI. "Contribution to a unified Eulerian modeling of fuel injection : from dense liquid to polydisperse evaporating spray." Thèse de doct. Universite Paris-Saclay, 2018 (cf. p. 75, 76, 81).

[32] M. FEBRES et D. LEGENDRE. "Enhancement of a 2D front-tracking algorithm with a non-uniform distribution of Lagrangian markers". In: *J. Comput. Phys.* 358 (2018), p. 173-200 (cf. p. 34, 36).
[33] M. FEBRES et D. LEGENDRE. "Enhancement of a 2D front-tracking algorithm with a non-uniform distribution of Lagrangian markers". In: *J. Comput. Phys.* 358 (2018), p. 173-200 (cf. p. 55).

[34] R. FEDKIW, T. ASLAM, B. MERRIMAN et S. OSHER. "A non-oscillatory Eulerian approach to interfaces in multi-material flows (the Ghost Fluid method)". In: *J. Comput. Phys.* 152 (1999). Sous la dir. de J. of COMPUTATIONAL PHYSICS, p. 457-492 (cf. p. 30, 73).

[35] D. FUSTER, G. AGBAGLAH, C. JOSSERAND, S. POPINET et S. ZALESKI. "Numerical simulation of droplets, bubbles and waves : state of the art". In : *Fluid Dyn. Res.* 41 (2009), p. 065001 (cf. p. 20).

[36] D. E. FYFE, E. S. ORAN et M. J. FRITTS. "Surface Tension and Viscosity with Lagrangian Hydrodynamics on a Triangular Mesh". In : *J. Comput. Phys.* 76 (1988), p. 349-384 (cf. p. 23).

[37] S. GANESAN et L. TOBISKA. "An accurate finite element scheme with moving meshes for computing 3D-axisymmetric interface flows". In : *Int. J. Numer. Methods Fluids* 57 (2008), p. 119-138 (cf. p. 23).

[38] S. GANESAN et L. TOBISKA. "Arbitrary Lagrangian-Eulerian finite-element method for computation of two-phase flows with soluble surfactants". In : *J. Comput. Phys.* 231 (2012), p. 3685-3702 (cf. p. 23).

[39] F. GIBOU, R. FEDKIW et S. OSHER. "A review of level-set methods and some recent applications". In : *J. Comput. Phys.* 353 (2018), p. 82-109 (cf. p. 28).

[40] J. GLIMM, M. GRAHAM, J. GROVE, X. LI, T. SMITH, D. TAN, F. TANGERMAN et Q. ZHANG. "Front tracking in two and three dimensions". In : *Computers & Mathematics with Applications* 35.7 (1998). Advanced Computing on Intel Architectures, p. 1-11 (cf. p. 28).

[41] J. GLIMM, X. LI, Y. LIU, Z. XU et N. ZHAO. "Conservative Front Tracking with Improved Accuracy". In : *SIAM J. Numer. Anal.* 41.5 (2003), p. 1926-1947 (cf. p. 29).

[42] J. GLIMM, J. W. GROVE, X. L. LI et D. C. TAN. "Robust Computational Algorithms for Dynamic Interface Tracking in Three Dimensions". In : *SIAM J. Sci. Comput.* 21.6 (2000), p. 2240-2256 (cf. p. 29).

[43] J. GLIMM, J. W. GROVE, X. L. LI, K.-M. SHYUE, Y. ZENG et Q. ZHANG. "Three-dimensional front tracking". In : *SIAM J. Sci. Comput.* 19.3 (1998), p. 703-727 (cf. p. 28).

[44] J. P. GOIS, A. NAKANO, L. G. NONATO et G. C. BUSCAGLIA. "Front tracking with moving-least-squares surfaces". In : *J. Comput. Phys.* 227 (2008), p. 9643-9669 (cf. p. 63).

[45] R. GOLDMAN. "Curvature formulas for implicit curves and surfaces". In : *Comput. Aided Geom. Des.* 22 (2005), p. 632-658 (cf. p. 34, 78, 79).

[46] S. GONG, W. MA et H. GU. "An experimental investigation on bubble dynamics and boiling crisis in liquid films". In : *Int. J. Heat Mass Transfer* 79 (2014), p. 694-703 (cf. p. 21).

[47] M. GOROKHOVSKI et M. HERRMANN. "Modeling Primary Atomization". In : *Annual Review of Fluid Mechanics* 40 (2008), p. 343-366 (cf. p. 20).

[48] R. GUILLAUMENT, S. VINCENT et J.-P. CALTAGIRONE. "An original algorithm for VOF based method to handle wetting effect in multiphase flow simulation". In : *Mech. Res. Commun.* 63 (2015), p. 26-32 (cf. p. 27).

[49] M. HERRMANN. "Detailed Numerical Simulations of the Primary Atomization of a Turbulent Liquid Jet in Crossflow". In : *J. Eng. Gas Turbines Power* 132 (juin 2010), p. 061506 (cf. p. 20).
[50] C. W. HIRT et B. D. NICHOLS. "Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries". In: *J. Comput. Phys.* 39 (1981), p. 201-225 (cf. p. 26).

[51] S. HYSING, S. TUREK, D. KUZMIN, N. PAROLINI, E. BURMAN, S. GANESAN et L. TOBISKA. "Quantitative benchmark computations of two-dimensional bubble dynamics". In : *Int. J. Numer. Methods Fluids* 60.11 (2009), p. 1259-1288 (cf. p. 138, 141, 146).

[52] M. JEMISON, M. SUSSMAN et M. SHASHKOV. "Filament capturing with the Multimaterial Moment-of-Fluid method". In : *J. Comput. Phys.* 28 (2015), p. 149-172 (cf. p. 27).

[53] D. JURIC et G. TRYGGVASON. "Computations of boiling flows". In : *Int. J. Multiphase Flow* 24.3 (1998), p. 387-410 (cf. p. 28).

[54] I. KATAOKA. "Local instant formulation of two-phase flow". In : *Int. J. Multiphase Flow* 12.5 (1986), p. 745-758 (cf. p. 24).

[55] S. KESHANI, W. R. W. DAUD, M. NOUROUZI, F. NAMVAR et M. GHASEMI. "Spray drying: An overview on wall deposition, process and modeling". In : *J. Food Eng.* 146 (2015), p. 152-162 (cf. p. 21).

[56] D.-A. KOFFI BI, M. TAVARES, E. CHENIER et S. VINCENT. "Accuracy and convergence of the curvature and normal vector discretizations for 3D static and dynamic front-tracking interfaces". In : *submitted in J. Comput. Phys.* (2021) (cf. p. 37).

[57] D.-A. KOFFI BI, M. TAVARES, É. CHÉNIER et S. VINCENT. "A Review of Geometrical Interface Properties for 3D Front-Tracking Methods". In : *Turbulence and Interactions*. Sous la dir. de M. DEVILLE, C. CALVIN, V. COUAILLIER, M. DE LA LLAVE PLATA, J.-L. ESTIVALÈZES, T. H. LÊ et S. VINCENT. Cham : Springer International Publishing, 2021, p. 144-149 (cf. p. 37).

[58] M. KUCHARIK, R. V. GARIMELLA, S. P. SCHOFIELD et M. J. SHASHKOV. "A comparative study of interface reconstruction methods for multi-material ALE simulations". In : *J. Comput. Phys.* 229 (2010), p. 2432-2452 (cf. p. 27).

[59] KUPRAT, ANDREW, KHAMAYSEH, AHMED, GEORGE, DENISE, LARKEY et LEVI. "Volume Conserving Smoothing for Piecewise Linear Curves, Surfaces, and Triple Lines". In : *J. Comput. Phys.* 172.1 (sept. 2001), p. 99-118 (cf. p. 60).

[60] T. LANGRISH et D. FLETCHER. "Spray drying of food ingredients and applications of CFD in spray drying". In : *Chem. Eng. Process.* 40 (2001), p. 345-354 (cf. p. 21).

[61] J. C. LASHERAS, E. VILLERMAUX et E. J. HOPFINGER. "Break-up and atomization of a round water jet by a high-speed annular air jet". In : *J. Fluid Mech.* 357 (1998), p. 351-379 (cf. p. 20).

[62] R. LEBAS, T. MENARD, P. BEAU, A. BERLEMONT et F. DEMOULIN. "Numerical simulation of primary break-up and atomization : DNS and modelling study". In : *Int. J. Multiphase Flow* 35 (2009), p. 247-260 (cf. p. 20, 30).

[63] A. H. LEFEBVRE et V. G. McDONELL. *Atomization and Sprays*. Second. Combustion : An International Series. CRC Press, 2017 (cf. p. 21).

[64] D. LEGENDRE et J. MAGNAUDET. "The lift force on a spherical bubble in a viscous linear shear flow". In : *J. Fluid Mech.* 368 (1998), p. 81-126 (cf. p. 23).

[65] R. J. LEVEQUE. "High-resolution conservative algorithms for advection in incompressible flows". In : *SIAM J. Numer. Anal.* 33.2 (avr. 1996), p. 627-665 (cf. p. 50, 63, 145).

[66] H. LHUISSIER et E. VILLERMAUX. "Bursting bubble aerosols". In : *J. Fluid Mech.* 696 (2012), p. 5-44 (cf. p. 19).

[67] X. LI, G. XU et Y. J. ZHANG. "Localized discrete Laplace-Beltrami operator over triangular mesh". In : *Computer Aided Geometric Design*. T. 39. 2015, p. 67-82 (cf. p. 73, 80, 88).
[68] K. LUO, H. PITSCH, M. PAI et O. DESJARDINS. "Direct numerical simulations and analysis of three-dimensional n-heptane spray flames in a model swirl combustor". In : *Proc. Combust. Inst.* 33 (2011), p. 2143-2152 (cf. p. 21).

[69] J. MAGNAUDET. *Méthodes numériques pour les problèmes d'inter interfaces*. Lalonde les Maures : Ecole de Printemps, mai 1997 (cf. p. 24).

[70] P. MARMOTTANT et E. VILLERMAUX. "On spray formation". In : *J. Fluid Mech.* 498 (2004), p. 73-111 (cf. p. 20).

[71] H. A. MAVROMATIS. "A single-sum expression for the overlap of two associated Legendre polynomials". In : *J. Phys. A : Math. Gen.* 32 (1999), p. 2601-2603 (cf. p. 146).

[72] H. A. MAVROMATIS et R. S. ALASSAR. "A Generalized Formula for the Integral of Three Associated Legendre Polynomials". In : *Appl. Math. Lett.* 12 (1999), p. 101-105 (cf. p. 146).

[73] N. MAX. "Weights for computing vertex normals from facets normals." In : *Journal of Graphics Tools* 4.2 (1999), p. 1-6 (cf. p. 79).

[74] R. MCDERMOTT et S. POPE. "The parabolic edge reconstruction method (PERM) for Lagrangian particle advection". In : *J. Comput. Phys.* 227.11 (2008), p. 5447-5491 (cf. p. 49, 89, 92, 94).

[75] T. MÉNARD, S. TANGUY et A. BERLEMONT. "Coupling level set/VOF/ghost fluid methods : Validation and application to 3D simulation of the primary break-up of a liquid jet". In : *Int. J. Multiphase Flow* 33 (2007), p. 510-524 (cf. p. 30).

[76] M. MEYER, M. DESBRUN, P. SCHRODER et A. H. BARR. "Discrete Differential-Geometry Operators for Triangulated 2-Manifolds". In : *Visualization and Mathematics III*. Sous la dir. de H.-C. HEGE et K. POLTHIER. Springer, 2003, p. 35-57 (cf. p. 33, 36, 75, 76, 81, 85, 90, 92, 107, 152, 154).

[77] E. MOHAMED, D. FLORENCE, de CHAISEMARTIN STÉPHANE, L. ADAM, M. THIBAULT et M. MARC. "Statistical modeling of the gas-liquid interface using geometrical variables : Toward a unified description of the disperse and separated phase flows". In : *Int. J. Multiphase Flow* 120 (2019), p. 103084 (cf. p. 36, 73).

[78] L. MORAWSKA, G. JOHNSON, Z. RISTOVSKI, M. HARGREAVES, K. MENGERSEN, S. C. B. C. CHAO, Y. LI et D. KATOSHEVSKI. "Size distribution and sites of origin of droplets expelled from the human respiratory tract during expiratory activities". In : *Aerosol Science* 40.3 (2009), p. 256-269 (cf. p. 21).

[79] C. MUNDO, M. SOMMERFELD et C. TROPEA. "Droplet-wall collisions : Experimental studies of the deformation and breakup process". In : *Int. J. Multiphase Flow* 21.2 (1995), p. 151-173 (cf. p. 21).

[80] N. NIKOLOPOULOS, A. THEODORAKAKOS et G. BERGELES. "A numerical investigation of the evaporation process of a liquid droplet impinging onto a hot substrate". In : *Int. J. Heat Mass Transfer* 50 (2007), p. 303-319 (cf. p. 21).

[81] W. F. NOH et P. WOODWARD. "SLIC (Simple line interface calculation)". In : *Proceedings of the Fifth International Conference on Numerical Methods in Fluid Dynamics June 28 – July 2, 1976 Twente University, Enschede*. Sous la dir. d'A. I. van de VOOREN et P. J. ZANDBERGEN. T. 59. Lecture Notes in Physics. Berlin, Heidelberg : Springer, Berlin, Heidelberg, juil. 1976 (cf. p. 26).

[82] S. NUKIYAMA. "The maximum and minimum values of the heat Q transmitted from metal to boiling water under atmospheric pressure." In : *Int. J. Heat Mass Transfer* 9 (1966), p. 1419-1433 (cf. p. 21).

[83] E. OLSSON et G. KREISS. "A conservative level set method for two phase flow". In : *J. Comput. Phys.* 210 (2005), p. 225-246 (cf. p. 28).
[84] E. OLSSON, G. KREISS et S. ZAHEDI. "A conservative level set method for two phase flow II". In : *J. Comput. Phys.* 225 (2007), p. 785-807 (cf. p. 28).

[85] S. OSHER et J. A. SETHIAN. "Front Propagating with curvature dependent speed : Algorithms based on Hamilton-Jacobi formulations". In : *J. Comput. Phys.* 79.1 (nov. 1988), p. 12-49 (cf. p. 27).

[86] A. PAUDEL, Z. A. WORKU, J. MEEUS, S. GUNS et G. V. den MOOTER. "Manufacturing of solid dispersions of poorly water soluble drugs by spray drying : Formulation and process considerations". In : *Int. J. Pharm.* (2012), p. 127-32 (cf. p. 21).

[87] C. S. PESKIN. "Numerical analysis of blood flow in the heart". In : *J. Comput. Phys.* 25.3 (1977), p. 220-252 (cf. p. 31).

[88] M. PEZO et V. STEVANOVIC. "Numerical prediction of critical heat flux in pool boiling with the two-fluid model". In : *Int. J. Heat Mass Transfer* 54 (2011), p. 3296-3303 (cf. p. 21).

[89] J. E. PILLIOD. *An Analysis of Piecewise Linear Interface Reconstruction Algorithms for Volume-of-fluid Methods*. Davis : University of California, sept. 1992 (cf. p. 27).

[90] J. E. PILLIOD et E. G. PUCKETT. "Second-order accurate volume-of-fluid algorithms for tracking material interfaces". In : *J. Comput. Phys.* 199 (2004), p. 465-502 (cf. p. 27).

[91] K. POLTHIER, J. HASS, D. HOFFMAN, A. JAFFE, H. ROSENBERG, R. SCHOEN et M. WOLF. "Computational aspects of discrete minimal surfaces." In : *Proc. of the Clay Summer school on Global Theory of Minimal Surfaces*. 2002 (cf. p. 75).

[92] S. POPINET. "An accurate adaptive solver for surface-tension-driven interfacial flows". In : *J. Comput. Phys.* 228.16 (2009), p. 5838-5866 (cf. p. 34-36).

[93] S. POPINET. "Numerical Models of Surface Tension". In : *Annu. Rev. Fluid Mech.* 50.1 (2018), p. 49-75 (cf. p. 31, 32, 35).

[94] S. POPINET et S. ZALESKI. "A front-tracking algorithm for accurate representation of surface tension." In : *Int. J. Numer. Meth. Fluids* 30 (1999), p. 775-793 (cf. p. 76, 78).

[95] H. POTTMANN, J. WALLNER, Y.-L. YANG, Y.-K. LAI et S.-M. HU. "Principal curvatures from the integral invariant viewpoint". In : *Comput. Aided Geom. Des.* 24.8 (2007). Discrete Differential Geometry, p. 428-442 (cf. p. 97).

[96] E. G. PUCKETT. "A volume-of-fluid interface tracking algorithm with applications to computing shock wave refraction". In : *Proceedings of the Fourth International Symposium on Computational Fluid Dynamics*. Sous la dir. de H. DWYER. Davis, CA, 1991, p. 933-938 (cf. p. 27).

[97] G. RACAH. "Theory of Complex Spectra. II". In : *Phys. Rev.* 62 (nov. 1942), p. 438-462 (cf. p. 146).

[98] Y. RENARDY et M. RENARDY. "PROST : A Parabolic Reconstruction of Surface Tension for the Volume-of-Fluid Method". In : *J. Comput. Phys.* 183 (2002), p. 400-421 (cf. p. 34).

[99] I. ROGHAIR, V. S. ANNALAND et J. KUIPERS. "An improved Front-Tracking technique for the simulation of mass transfer in dense bubbly flows". In : *Chem. Eng. Sci.* 152 (oct. 2016), p. 351-369 (cf. p. 50, 55, 59, 89).

[100] S. RUSINKIEWICZ. "Estimating curvatures and their derivatives on triangle meshes". In : *Proceedings. 2nd International Symposium on 3D Data Processing, Visualization and Transmission, 2004. 3DPVT 2004.* 2004, p. 486-493 (cf. p. 97).

[101] M. SATO. "General formula of the Racah coefficient". In : *Progress of Theoretical Physics* 13.4 (avr. 1955) (cf. p. 146).

[102] S. SAZHIN. *Droplets and Sprays*. Springer, 2014 (cf. p. 20).
[103] R. SCARDOVELLI et S. ZALESKI. "Interface reconstruction with least-square fit and split Eulerian-Lagrangian advection". In : *Int. J. Numer. Methods Fluids* 41 (2003), p. 251-274 (cf. p. 27).

[104] SCARDOVELLI et ZALESKI. "Direct Numerical Simulation of free-surface and interfacial flow". In : *Annu. Rev. Fluid Mech.* 31 (jan. 1999), p. 567-603 (cf. p. 32).

[105] D. SÉBILLEAU. "On the computation of the integrated products of three spherical harmonics". In : *J. Phys. A : Math. Gen.* 31.34 (août 1998), p. 7157-7168 (cf. p. 146).

[106] S. SHIN. "Computation of the curvature field in numerical simulation of multiphase flow". In : *J. Comput. Phys.* 222.2 (2007), p. 872-878 (cf. p. 32).

[107] S. SHIN, S. I. ABDEL-KHALIK, V. DARU et D. JURIC. "Accurate representation of surface tension using the level contour reconstruction method". In : *J. Comput. Phys.* 203 (2005), p. 493-516 (cf. p. 29, 32).

[108] S. SHIN et D. JURIC. "Modeling Three-Dimensional Multiphase Flow Using a Level Contour Reconstruction Method for Front Tracking without Connectivity". In : *J. Comput. Phys.* 180.2 (2002), p. 427-470 (cf. p. 29, 73).

[109] S. SHIN et D. JURIC. "A hybrid interface method for three-dimensional multiphase flows based on front tracking and level set techniques". In : *Int. J. Numer. Methods Fluids* 60.7 (2009), p. 753-778 (cf. p. 32).

[110] S. SHIN, I. YOON et D. JURIC. "The Local Front Reconstruction Method for direct simulation of two- and three-dimensional multiphase flows". In : *J. Comput. Phys.* 230.17 (2011), p. 6605-6646 (cf. p. 29).

[111] L. SHUI, J. C. EIJKEL et A. van den BERG. "Multiphase flow in microfluidic systems – Control and applications of droplets and interfaces". In : *Adv. Colloid Interface Sci.* 133 (2007), p. 35-49 (cf. p. 21).

[112] S. SOM et S. AGGARWAL. "Effects of primary breakup modeling on spray and combustion characteristics of compression ignition engines". In : *Combust. Flame* 157 (2010), p. 1179-1193 (cf. p. 21).

[113] D. SUN et W. TAO. "A coupled volume-of-fluid and level set (VOSET) method for computing incompressible two-phase flows". In : *Int. J. Heat Mass Transfer* 53 (2010), p. 645-655 (cf. p. 30).

[114] M. SUSSMAN et E. G. PUCKETT. "A Coupled Level Set and Volume-of-Fluid Method for Computing 3D and Axisymmetric Incompressible Two-Phase Flows". In : *J. Comput. Phys.* 162 (2000), p. 301-337 (cf. p. 29).

[115] M. SUSSMAN, P. SMEREKA et S. OSHER. "A level set approach for computing solutions to incompressible two-phase flow". In : *J. Comput. Phys.* 114 (1994), p. 146-159 (cf. p. 28).

[116] G. TAUBIN. "A signal processing approach to fair surface design." In : *SIGGRAPH '95*. 1995, p. 351-358 (cf. p. 75).

[117] M. TAVARES. "Simulation et modélisation multi-échelle d'écoulements diphasiques". Thèse de doct. Université Paris-Est, juin 2019 (cf. p. 30, 34, 36).

[118] M. TAVARES, D.-A. KOFFI BI, E. CHENIER et S. VINCENT. "A Two-Dimensional Second Order Conservative Front-Tracking Method with an Original Marker Advection Approach Based on Jump Relations". In : *Communications in Computational Physics* 27.5 (2020), p. 1550-1589 (cf. p. 37, 45).
[119] M. TAVARES, D.-A. KOFFI BI, E. CHÉNIER et S. VINCENT. "A Front-Tracking Method for Multiphase Flows with a Sharp Interface Representation". In: *Turbulence and Interactions*. Sous la dir. de M. DEVILLE, C. CALVIN, V. COUAILLIER, M. DE LA LLAVE PLATA, J.-L. ESTIVALÈZES, T. H. LÊ et S. VINCENT. Cham : Springer International Publishing, 2021, p. 189-195 (cf. p. 37).

[120] D. J. TORRES et J. U. BRACKBILL. "The Point-Set Method : Front-Tracking without Connectivity". In : *J. Comput. Phys.* 165 (2000), p. 620-644 (cf. p. 29, 30, 47).

[121] P. TRONTIN, S. VINCENT, J. ESTIVALEZES et J. CALTAGIRONE. "A subgrid computation of the curvature by a particle/level-set method. Application to a front-tracking/ghost-fluid method for incompressible flows". In : *J. Comput. Phys.* 231.20 (2012), p. 6990-7010 (cf. p. 78).

[122] G. TRYGGVASON, B. BUNNER, A. ESMAEELI, D. JURIC, N. AL-RAWAHI, W. TAUBER, J. HAN, S. NAS et Y.-J. JAN. "A Front-Tracking Method for the Computations of Multiphase Flow". In : *J. Comput. Phys.* 169.2 (2001), p. 708-759 (cf. p. 29-31, 76, 77).

[123] G. TRYGGVASON, R. SCARDOVELLI et S. ZALESKI. *Direct Numerical Simulations of Gas-Liquid Multiphase Flows*. 1re éd. Cambridge University Press, 2011 (cf. p. 31, 77, 82, 85, 90, 93, 95, 107).

[124] S. O. UNVERDI et G. TRYGGVASON. "A front-tracking method for viscous, incompressible, multi-fluid flows". In : *J. Comput. Phys.* 100.1 (1992), p. 25-37 (cf. p. 28).

[125] G. VAUDOR, T. MÉNARD, W. ANISZEWSKI, M. DORING et A. BERLEMONT. "A consistent mass and momentum flux computation method for two phase flows. Application to atomization process". In : *Computers and Fluids* 152 (2017), p. 204-216 (cf. p. 73).

[126] E. VILLERMAUX. "Fragmentation". In : *Annu. Rev. Fluid Mech.* 39 (2007), p. 419-446 (cf. p. 19).

[127] E. VILLERMAUX et B. BOSSA. "Drop fragmentation on impact". In : *J. Fluid Mech.* 668 (2011), p. 412-435 (cf. p. 21).

[128] M. WARDETZKY. "Convergence of the Cotangent Formula : An Overview". In : *Discrete Differential Geometry* 38 (2008), p. 275-286 (cf. p. 88).

[129] S. WILDEMAN, C. W. VISSER, C. SUN et D. LOHSE. "On the spreading of impacting drops". In : *J. Fluid Mech.* 805 (2016), p. 636-655 (cf. p. 21).

[130] F. XIAO, Y. HONMA et T. KONO. "A simple algebraic interface capturing scheme using hyperbolic tangent function". In : *Int. J. Numer. Methods Fluids* 48 (2005), p. 1023-1040 (cf. p. 27).

[131] F. XIAO, S. LI et C. CHEN. "Revisit to the THINC scheme : A simple algebraic VOF algorithm". In : *J. Comput. Phys.* 230 (2011), p. 7086-7092 (cf. p. 27).

[132] G. XU. "Discrete Laplace-Beltrami operators and their convergence". In : *Comput. Aided Geom. Des.* 21.8 (2004). Geometric Modeling and Processing 2004, p. 767-784 (cf. p. 88).

[133] Y.-L. XU. "Fast evaluation of the Gaunt coefficients". In : *Math. Comput.* 65.216 (oct. 1996), p. 1601-1612 (cf. p. 146).

[134] Y.-L. XU. "Fast evaluation of Gaunt coefficients : recursive approach". In : *J. Comput. Appl. Math.* 85 (1997), p. 53-65 (cf. p. 146).

[135] D. L. YOUNGS. "Time-dependent multi-material flow with large fluid distortion." In : *Numerical methods for Fluid Dynamics* 24(2) (1982), p. 273-285 (cf. p. 26).

[136] D.-h. YU. *Natural Boundary Integral Method and Its Applications*. T. 539. Mathematics and Its Applications. Springer Netherlands, 2002 (cf. p. 73).
[137] D. ZHAO, Y. TOBA, Y. SUZUKI et S. KOMORI. "Effect of wind waves on air-sea gas exchange : proposal of an overall CO2 transfer velocity formula as a function of breaking-wave parameter". In : *Tellus* (2003), p. 478-487 (cf. p. 19).

[138] A. Z. ZINCHENKO, M. A. ROTHER et R. H. DAVIS. "A novel boundary-integral algorithm for viscous interaction of deformable drops." In : *Phys. Fluids* 9 (1997) (cf. p. 78, 79, 83).

[139] D. ZUZIO, J.-L. ESTIVALEZES, P. VILLEDIEU et G. BLANCHARD. "Numerical simulation of primary and secondary atomization". In : *Comptes Rendus Mécanique* 341 (2013), p. 15-25 (cf. p. 73).