# A multipole based Boundary Element Method for moving boundary problem in axisymmetric geometry

Jitendra Singh, Alain Glière, Jean-Luc Achard

► To cite this version:

Jitendra Singh, Alain Glière, Jean-Luc Achard. A multipole based Boundary Element Method for moving boundary problem in axisymmetric geometry. 30 th International Conference on Boundary Elements and other Mesh Reduction Methods., Jul 2009, Maribor, Slovenia. pp.43-54, 10.2495/BE080051. hal-00392329

HAL Id: hal-00392329

https://hal.science/hal-00392329v1

Submitted on 13 Feb 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A multipole based Boundary Element Method for moving boundary problems in axisymmetric geometry

J. Singh<fnref n="1" />, A. Glière<fnref n="1" /> & J.-L. Achard<fnref n="1" />,<fnref n="2" />

## Abstract

The Boundary Element Method (BEM) has been an effective method for modelling a number of particular problems which may be described by potential flows or, in a quite different vein, by creeping motions. However when it comes to moving boundary flows the BEM requires a significant amount of computational time. In the context of accelerating the BEM, the Fast Multipole Method (FMM) has been used for two and three-dimensional problems but they are not yet reported for axisymmetric ones. Our purpose is to present a multipole based boundary element method (MM-BEM) for moving boundary problems in the axisymmetric case. The proposed method takes advantage of grouping the ring sources in the axisymmetric domain, in order to reduce the amount of direct computations. Direct computations are only performed when the ring sources are located close to the evaluation point. Here, MM-BEM is implemented to simulate the impact of a drop onto a liquid surface, modelled assuming potential flow. It is shown that the proposed method provides the same results as the conventional BEM.

Keywords: multipole expansion, multipole coefficient, Boundary Element Method, potential flow, axisymmetric problem, moving boundary.

## 1 Introduction

A handful of practical applications characterized by moving boundaries can be modelled in axisymmetric environment, specifically the applications involving drops and bubbles [1–3]. The objective is to keep track of the moving boundary at different instants of time. The interest of the axisymmetric formulation lies in

<figure><img src="image_2.png" /></figure>

<footnote n="1">¹CEA-LETI MINATEC, France</footnote>

<footnote n="2">²LEGI, Microfluidics, Interfaces & Particles Team, France</footnote>
the fact that, while keeping the whole three-dimensional physical aspect of the problem, the solution is carried out in a two-dimensional domain. Thus, the implementation becomes easier and less computational resources are required.

The BEM has been an efficient technique for tracking moving boundaries, especially when dealing with potential and viscous fluid flow assumptions [4]. The method has tremendous advantage that the boundary condition is treated neatly and resolved very accurately. However, simulating highly distorted interfaces require distributing a large number of nodes, in which case the BEM requires a significant amount of computational time. Moreover, if the problem is solved at several time sequences the aforesaid issue turns demanding. The iterative solver has particular advantage for moving boundary problems, since the solution of the previous time step provides a good initial guess for the following step. However, the free surfaces are Dirichlet type boundary, for which the integral equation is Fredholm integral equation of the first kind, which is not sufficiently well conditioned. The numerical solution requires a large number of iterations.

In the context of accelerating the BEM, the fast iterative solver techniques, like FMM, tackles efficiently the issues of time and memory in two and three-dimensional problems. However, FMM is not applicable for axisymmetric problems, due to the presence of ring sources.

In this work we implement a multipole method coupled with BEM (MM-BEM) for moving boundary flow in the extensively studied application of drop splash, wherein a droplet is subjected to impinge on a liquid free surface. The application represents an interesting case of highly distorted moving interface and is modelled by assuming the potential flow of a fluid in an axisymmetric environment. The multipole expansion is first used to separate the ring sources from the evaluation point, which in turn enables to make a grouping structure for the ring sources. The contribution of a group is then calculated at the axis of symmetry in terms of multipole coefficients and this information is efficiently used at the evaluation point. In this manner we obtain a fast matrix vector multiplication which is coupled with an iterative solver to accelerate the solution procedure.

Outline of the paper is as follows. In the second and third sections we present the gradient and the multipole formulation of the Dirichlet problem and the general features of the MM-BEM technique. The application to liquid drop impact is presented in the fourth section. The two last sections consist of numerical results and a brief conclusion.

## 2 Axisymmetric gradient formulation

In the context of the potential flow theory, Dirichlet boundary conditions are applied to moving free surfaces. The integral equation formulation for the Dirichlet problem results in a Fredholm integral equation of the first kind, which is not well conditioned in terms of numerical solution. However, this Fredholm integral equation can be converted to a second kind equation for the boundary distribution of the normal derivative, which is better conditioned. The gradient formulation for the Laplace equation is well established [3], and, for

<figure><img src="image_3.png" /></figure>
completeness, we present it here for the axisymmetric case. Starting with the integral equation for the three dimensions Laplace equation defined over the surface S which bounds the domain of interest [4]

$$c(\mathbf{x})\varphi(\mathbf{x}) = \int_{S} \frac{\partial\varphi}{\partial n'}(\mathbf{x}')G(\mathbf{x}, \mathbf{x}')dS' - \int_{S} \varphi(\mathbf{x}')\frac{\partial G}{\partial n'}(\mathbf{x}, \mathbf{x}')dS', \quad (1)$$

where $\varphi$ is the velocity potential, $G(\mathbf{x},\mathbf{x}')=1/4\pi r$ is the free space Green's function and $r = |\mathbf{x} - \mathbf{x}'|$ is the distance between the evaluation point $\mathbf{x}$ and the source point $\mathbf{x}'$. Moreover, $\partial\varphi/\partial n = \mathbf{n} \cdot \nabla\varphi$ denotes the normal derivative with $\mathbf{n}$ the unit outward normal vector on S, and $c(\mathbf{x})$ depends on the local geometry of S at $\mathbf{x}$ ($c=1/2$ if S is regular at $\mathbf{x}$). From eqn (1) an integral equation for the normal derivative of the potential at the evaluation point can be derived as

$$c(\mathbf{x}) \frac{\partial \varphi}{\partial n}(\mathbf{x}) = \int_{S} \frac{\partial \varphi}{\partial n'}(\mathbf{x}') [\mathbf{n} \cdot \nabla G(\mathbf{x}, \mathbf{x}')] dS' - \int_{S} \varphi(\mathbf{x}') \left[ \mathbf{n} \cdot \nabla \left( \frac{\partial G}{\partial n'}(\mathbf{x}, \mathbf{x}') \right) \right] dS'. \quad (2)$$

The first and second surface integrals in eqn (2) express boundary distribution of point sources and point sources dipoles respectively. This equation is hypersingular due to the kernel in the second integral and can get regularized by following equivalence between doublet and vortex distributions (see e.g., reference [5])

$$\int \varphi(\mathbf{x}') \left[ \mathbf{n} \cdot \nabla \left( \frac{\partial G}{\partial n'}(\mathbf{x}, \mathbf{x}') \right) \right] dS' = \int \mathbf{n} \cdot \left[ (\nabla' \varphi \times \mathbf{n}') \times \nabla' G(\mathbf{x}, \mathbf{x}') \right] dS'. \quad (3)$$

Let the z-axis in cylindrical polar coordinate system $(\rho, \phi, z)$ coincide with the axis of revolution of the domain of interest. We assume that the velocity potential is rotationally symmetric about the z-axis. The azimuthal component of the outer normal vanishes over the boundary S. The velocity potential and its gradient depend only on $\rho$ and $z$. The meridian plane where the evaluation point $\mathbf{x}$ lies is arbitrary. For sake of simplicity it corresponds here to $\phi=0$. The differential area $dS$ is expressed in the form $\rho d\phi d\Gamma$, where $d\Gamma$ is the differential arc length along the contour $\Gamma$ of the surface in a meridian plane. An important consequence of axial symmetry is to simplify the integrand in the right hand side of eqn (3), since

$$[(\nabla' \varphi \times \mathbf{n}') \times \nabla' G] = - \frac{\partial \varphi}{\partial S'}(\mathbf{x}') [\mathbf{i}_{\phi'} \times \nabla' G], \quad (4)$$

with $\partial\varphi/\partial S'$ being the tangential component of the velocity, acting here as a vortex strength and $\mathbf{i}_{\phi'}$ is the unit vector on S at $\mathbf{x}'$ perpendicular to the meridian plane. The sense of $\mathbf{s}$, the unit tangent vector, is such that $(\mathbf{s}, \mathbf{n}, \mathbf{i}_{\phi'})$ is right-handed in that order. Taking into account all the notations and axisymmetric simplifications, we recast eqn (2) into the form

$$c(\mathbf{x}) \frac{\partial \varphi}{\partial n}(\mathbf{x}) = \int_{\Gamma} \frac{\partial \varphi}{\partial n'}(\mathbf{x}') \left( \mathbf{n}_{\rho} \frac{\partial G^{AX}}{\partial \rho} + \mathbf{n}_{z} \frac{\partial G^{AX}}{\partial z} \right) \rho(\mathbf{x}') d\Gamma' \\ - \int_{\Gamma} \frac{\partial \varphi}{\partial S'}(\mathbf{x}') \left( \mathbf{n}_{\rho} \frac{\partial \tilde{G}^{AX}}{\partial z} + \mathbf{n}_{z} \frac{\partial \tilde{G}^{AX}}{\partial \rho'} \right) \rho(\mathbf{x}') d\Gamma', \quad (5)$$

<figure><img src="image_4.png" /></figure>
where

$$G^{AX} = \frac{1}{4\pi} \int_{0}^{2\pi} \frac{d\phi'}{|\mathbf{x} - \mathbf{x}'|} \quad \text{and} \quad \tilde{G}^{AX} = \frac{1}{4\pi} \int_{0}^{2\pi} \frac{\cos\phi'}{|\mathbf{x} - \mathbf{x}'|} d\phi' \quad (6)$$

are free space axisymmetric Green's function. In the conventional BEM formulation, explicit relations for $G^{AX}$ and $\tilde{G}^{AX}$ are expressed in terms of complete elliptic integrals of the first and second kinds. We rather adopt a multipole expansion based grouping technique.

### 3 Axisymmetric multipole formulation

In this section we describe the multipole techniques for the solution of eqn (5) in order to accelerate the matrix vector multiplication for iterative solver. The evaluation and source points are first separated and then an azimuthal integration is carried out to obtain the axial multipole moments. We begin with the expansion of the Green's function associated with the three dimension Laplace equation in spherical coordinates system $(r, \theta, \phi)$ [6]

$$\frac{1}{4\pi} \frac{1}{|\mathbf{x} - \mathbf{x}'|} = \frac{1}{4\pi} \sum_{n=0}^{\infty} \sum_{m=0}^{n} \varepsilon_m \frac{(n-m)!}{(n+m)!} \frac{r'^n}{r^{n+1}} P_n^m(\cos\theta) P_n^m(\cos\theta') \cos m(\phi - \phi'), \quad (7)$$

where $\varepsilon_m$ is the Neumann factor, $\varepsilon_0=1$ and $\varepsilon_i=2$ ($i=1,2,\text{etc.}$) and $P_n^m$ are the associated Legendre functions. To obtain the series expansion of the free space axisymmetric Green function $G^{AX}$ we integrate the above series in the azimuthal direction, keeping the origin on the axis of symmetry

$$G^{AX} = \frac{1}{4\pi} \sum_{n=0}^{\infty} \sum_{m=0}^{n} \varepsilon_m \frac{(n-m)!}{(n+m)!} \frac{r'^n}{r^{n+1}} P_n^m(\cos\theta) P_n^m(\cos\theta') \int_{0}^{2\pi} \cos m(\phi-\phi') d\phi'. \quad (8)$$

Since the evaluation point lies in the meridian plane, $\phi=0$ and

$$\int_{0}^{2\pi} \cos m(\phi - \phi') d\phi' = \begin{cases} 0, & m \neq 0 \\ 2\pi, & m \to 0 \end{cases}. \quad (9)$$

Therefore, only the $m=0$ terms survive and we obtain

$$G^{AX} = \frac{1}{2} \sum_{n=0}^{\infty} \frac{r'^n}{r^{n+1}} P_n(\cos\theta) P_n(\cos\theta'). \quad (10)$$

The multipole coefficients for the first integral in eqn (5) are obtained by substituting eqn (10), written as

$$\int_{\Gamma} \frac{\partial \varphi}{\partial n'}(\mathbf{x}') \left( n_{\rho} \frac{\partial G^{AX}}{\partial \rho} + n_{z} \frac{\partial G^{AX}}{\partial z} \right) \rho(\mathbf{x}') d\Gamma' \\ = \sum_{n=0}^{\infty} \left( n_{\rho} \frac{\partial}{\partial \rho} \left( \frac{P_n(\cos\theta)}{r^{n+1}} \right) + n_{z} \frac{\partial}{\partial z} \left( \frac{P_n(\cos\theta)}{r^{n+1}} \right) \right) M_n(O) \quad (11)$$

with $M_n(O)$ being the multipole coefficients at point $O$ due to the ring sources in $\Gamma$, given by

$$M_n(O) = \int_{\Gamma} \frac{\partial \varphi}{\partial n'}(\mathbf{x}') r'^n P_n(\cos\theta') \rho(\mathbf{x}') d\Gamma'. \quad (12)$$

<figure><img src="image_5.png" /></figure>
The series expansion of the second free space axisymmetric Green's function $\tilde{G}^{AX}$ is obtained as follows:

$$ \tilde{G}^{AX} = \frac{1}{4\pi} \sum_{n=0}^{\infty} \sum_{m=0}^{n} \varepsilon_m \frac{(n-m)!}{(n+m)!} r'^{n} P_n^m(\cos\theta) P_n^m(\cos\theta') \int_0^{2\pi} \cos\phi' \cos m(\phi-\phi') d\phi'. \quad (13) $$

Integrating the above series in the azimuthal direction, while keeping the origin on the axis of symmetry and the evaluation point in the meridian plane, we get

$$ \int_{0}^{2\pi} \cos \phi' \cos m(\phi - \phi') d\phi' = \begin{cases} 0, & m \neq 1 \\ \pi, & m \to 1 \end{cases}. \quad (14) $$

Only the terms $m=1$ survive and the eqn (13) becomes

$$ \tilde{G}^{AX} = \frac{1}{2} \sum_{n=0}^{\infty} \frac{(n-1)!}{(n+1)!} \frac{r'^n}{r^{n+1}} P_n^1(\cos\theta) P_n^1(\cos\theta'). \quad (15) $$

The multipole coefficients for the second integral in eqn (5) are obtained by substituting eqn (15), written as

$$ \begin{aligned} & \int_{\Gamma} \frac{\partial \varphi}{\partial s'}(\mathbf{x}') \left( n_{\rho} \frac{\partial \tilde{G}^{AX}}{\partial z} + n_z \frac{\partial G^{AX}}{\partial \rho'} \right) \rho(\mathbf{x}') d\Gamma' \\ & = \sum_{n=0}^{\infty} n_{\rho} \frac{\partial}{\partial z} \left( \frac{(n-1)! P_n^1(\cos\theta)}{(n+1)! r^{n+1}} \right) \tilde{M}_n(O) + \sum_{n=0}^{\infty} n_z \left( \frac{P_n(\cos\theta)}{r^{n+1}} \right) \hat{M}_n(O) \end{aligned} \quad (16) $$

where the multipole coefficients are given by

$$ \tilde{M}_n(O) = \int_{\Gamma} \frac{\partial \varphi}{\partial s'}(\mathbf{x}') r'^n P_n^1(\cos\theta) \rho(\mathbf{x}') d\Gamma', \quad (17) $$

$$ \hat{M}_n(O) = \int_{\Gamma} \frac{\partial \varphi}{\partial s'}(\mathbf{x}') r'^{n-1} \left[ n P_n(\cos\theta) \sin\theta + P_n^1(\cos\theta) \cos\theta \right] \rho(\mathbf{x}') d\Gamma'. \quad (18) $$

From the above expansions it is clear that the ring sources and the evaluation point are separated. Therefore, a grouping structure for the ring sources can be employed. The contribution at the evaluation point due to the ring sources in a group is calculated using the multipole coefficients corresponding to the group.

The above multipole coefficients are derived by assuming that the source ring is nearer to the expansion centre $O$ than the evaluation point. In fact in the literature they are termed as exterior multipole moments. Since the interest is towards maximizing the multipole calculation we use also the interior multipole coefficients, which can be derived on the same lines [7].

Admissibility criterion for the multipole calculations depends upon the ratio of distance of the group from the expansion centre to that of distance of the evaluation point. The distance of the group from the expansion centre is taken as the average of the distances of the elements which are contained in the group from the expansion centre. The contribution at the evaluation point in group D (Figure 1) from the ring sources contained in group A uses the exterior

<figure><img src="image_6.png" /></figure>
multipoles coefficients for group D computed at expansion centre *O*. In the vice-versa case, (ring sources in group D and evaluation point in group A) the contribution is calculated by interior multipoles at expansion centre *O*.

The description of the method based on multipole technique is as follows:

- First, the ring sources are grouped. This has been done by defining an arc length parameter along the contour, as shown in Figure 1.

- Second, the positions of the multipole expansion centres are determined, based on the geometry of the problem. For example the points *O* and *O'* are sufficiently good choice for the cylindrical geometry (Figure 1).

- Third, the calculation of the interior and the exterior multipoles for the grouped ring sources at the multipole expansion centres are performed.

- Finally, the contribution at the evaluation point due to a particular group is computed by multipole coefficients if the admissibility criterion is satisfied; otherwise the direct calculations are used for each element of the group.

<figure><img src="image_7.png" /><figcaption>Figure 1: Grouping structure and expansion centres for a cylindrical geometry.</figcaption></figure>

# 4 Application to drop impact

The impact of a liquid drop onto a liquid surface, which is an example of a highly distorted moving interface, is chosen as an application case for the multipole based BEM. Numerous theoretical, experimental and numerical studies have been performed during more than one century in order to better understand this phenomenon which, apart from being interesting and beautiful, is of considerable practical importance. Among other numerical techniques, the Volume Of Fluid [8], the level-set [9] and the Boundary Element Method [1] have been used.

## 4.1 Problem statement

An axisymmetric liquid domain, bounded in the meridian plane (*ρ*,*z*) by the liquid-gas interface Γ is considered (Figure 2). Due to the impulsive nature of the

<figure><img src="image_8.png" /></figure>
drop impact, an irrotational flow model is appropriate. A harmonic velocity potential subject to Dirichlet boundary conditions, can thus be defined. The local time derivative of the velocity potential is calculated by combining the normal momentum balance with the Laplace pressure jump to obtain the Bernoulli's equation. A dimensionless formulation is used:

$$ \frac{\partial \varphi}{\partial t} = -\frac{v^2}{2} - \frac{H}{We} - \frac{z}{Fr} - \frac{2}{Re} \frac{\partial^2 \varphi}{\partial n^2} \quad (19) $$

<figure><img src="image_9.png" /><figcaption>Figure 2: Initial geometry and potential. The neck region is smoothed.</figcaption></figure>

In this equation $v$ is the dimensionless velocity, $H$ is the total curvature and $Re$, $Fr$ and $We$ are respectively the Reynolds, Froude and Weber numbers. The drop radius $R$ is adopted as length scale and $\sqrt{\sigma/R\rho_L}$ as velocity scale, where $\sigma$ and $\rho_L$ are the surface tension and the liquid density. The ratio between the length and velocity scale gives the time scale. Viscous effects are partially taken into account as the normal viscous stress appears in eqn (19). The complete description of this model can be found in Georgescu et al. [2].

In order to avoid topological change, the falling spherical drop is initially placed in contact with the liquid surface. A negligibly small volume of the sphere is cut-out and, to avoid a singularity of the capillary pressure, the intersection of the sphere with the initially flat liquid surface is smoothed in the meridian plane by a circular element (Figure 2). The drop initially falls with a constant velocity onto an initially still liquid surface. The discontinuity of the velocity potential in the neck region is eliminated by means of the quadratic smoothing procedure described by Weiss and Yarin [10].

## 4.2 Solution procedure

The application represents a transient free-boundary problem that repeatedly involves two types of calculations: (a) solution of Laplace equation for normal components of the velocity, (b) updating of the potential at the forthcoming time step by using time marching scheme of Runge-Kutta and displacement of the interface.

Within the first type of calculation, the Laplace equation is solved by the Boundary Element Method using gradient formulation. First order elements are

<figure><img src="image_10.png" /></figure>
used to approximate the unknowns. The vortex strength is calculated by interpolating the given velocity potential over the interface by the cubic splines. The solution is obtained by grouping source/vortex rings using multipole techniques discussed in previous section coupled with the GMRES iterative solver without preconditioning. In order to maximize the number of multipole calculations the multipole coefficients are calculated at two different expansion centres. Two expansion centres are used and their positions are evaluated at each time step, following the deformation of the interface. The normal velocity obtained at previous time step is used as an initial guess for the iterative solver for the subsequent time step.

Within the second type of calculation, the Bernoulli's equation is solved by 4^th^ order explicit Runge-Kutta scheme for the local time derivative of the velocity potential. The time steps are chosen according to the stability criterion described by Georgescu et al. [2]. The interface position is then updated according to the normal velocity component.

**5 Results**

The presented MM-BEM has been implemented for the axisymmetric static potential problems. With 18 terms in the expansion, this method provides accurate results when compared to analytical solution of Dirichlet or Neumann problems. The computational time is almost reduced by half when compared to the conventional BEM [7].

<figure><img src="image_11.png" /><figcaption>Figure 3: Comparison of numerical and experimental profiles (Fr=90, We=43) obtained by MM-BEM and Liow [11]. Bottom image reproduced with permission from Cambridge University Press.</figcaption></figure>

<figure><img src="image_12.png" /></figure>
The method is adapted for free surface flows and specifically to the chosen application of a drop impinging on a liquid surface. The number of nodes on the interface varies from 20 to 150 and the multipole series expansion is truncated after 21 terms. The comparison of the free surface profiles obtained by simulation through MM-BEM and experiments [11] at different instants are shown in Figure 3 (Fr=90, We=43). The qualitative agreement is correct.

Furthermore, a 10% of total CPU time saving is obtained when using MM-BEM instead of the conventional BEM. Simulated profiles are identical in both methods.

# 6 Conclusion

A multipole based BEM has been implemented in an axisymmetric environment for the simulation of moving boundary problem of drop impact onto a liquid surface. By multipole techniques the interaction between individual nodes and ring sources in far field has been replaced by the interaction between node and group of ring sources. It has been shown that results obtained through MM-BEM matches well with available experimental data.

Although the time reduction is not significant yet, there is a wide scope of improvement in MM-BEM. The tuning of the parameters, for example, the number of terms in multipole series expansion, efficient choice of multipole centre in order to increase the validity of the far field and optimizing the admissibility criterion can further lead to more CPU time saving. To the best of our knowledge, this is the first attempt to accelerate the axisymmetric moving boundary value problem using a multipole expansion based technique.

# References

[1] Oguz, H.N. & Prosperetti, A., Bubble entrainment by the impact of drops on liquid surfaces. *Journal of Fluid Mechanics*, 219, pp. 143–179, 1990.

[2] Georgescu, S.C., Achard, J.-L. & Canot, E., Jet drops ejection in bursting gas bubble processes. *European Journal of Mechanics B-Fluids*, 21(2), pp. 265–280, 2002.

[3] Pozrikidis, C., Three-dimensional oscillations of inviscid drop induced by surface tension. *Computer and Fluids*, 30, pp. 417–444, 2001.

[4] Brebbia, C.A., Telles, J.C.F. & Wrobel, L.C., *Boundary Element Techniques*, Springer-Verlag: Berlin and New York, 1984.

[5] Brockett, T.E., Kim, M-H. & Park, J-H., Limiting forms for surface singularity distributions when the field point is on the surface. *Journal of Engineering Mathematics*, 23, pp. 201–31, 1989.

[6] Morse, P.M. & Feshbach, H., *Methods of theoretical physics*, McGraw-Hill: New York, 1953.

[7] J. Singh, A. Glière & J.-L. Achard, Multipole accelerated BEM for axisymmetric potential problem. (in preparation).

<figure><img src="image_13.png" /></figure>
[8] Morton, D., Rudman, M. & Liow, J. L., An investigation of the flow regimes resulting from splashing drops. *Physics of Fluids, 12*(4), pp. 747–763, 2000.

[9] Watanabe, Y., Saruwatari, A., & Ingram, D.M., Free-surface flows under impacting droplets. *Journal of Computational Physics*, 227(4), pp. 2344–2365, 2008.

[10] Weiss, D.A. & Yarin A.L., Single drop impact onto liquid films: neck distortion, jetting, tiny bubble entrainment, and crown formation. *Journal of Fluid Mechanics*, 385, pp. 229–254, 1999.

[11] Liow, J.L., Splash formation by spherical drops. *Journal of Fluid Mechanics*, 427, pp. 73–105, 2001.

<figure><img src="image_14.png" /></figure>