Modelling and Numerical Simulation of Binary Droplet Collisions Under Extreme Conditions

Abstract

The complexity of binary droplet collisions strongly increases in case of immiscible liquids with the occurrence of triple lines or for high energetic collisions, where strong rim instabilities lead to the spattering of satellite droplets. To cope with such cases, the Volume of Fluid method is extended by an efficient interface reconstruction, also applicable to multi-material cells of arbitrary configuration, as well as an enhanced continuous surface stress model for accurate surface force computations, also applicable to thin films. For collisions of fully wetting liquids, excellent agreement to experimental data is achieved in different collision regimes. High-resolution simulations predict droplet collisions in the spattering regime and provide detailed insights into the evolution of the rim instability. Another challenge is the numerical prediction of the collision outcome in the bouncing or coalescence region, where the rarefied gas dynamics in the thin gas film determines the collision result. To this end, an important step forward became possible by modelling the pressure in the gas film. With the introduction of an interior collision plane within the flow domain, it is now possible to simulate droplet collisions with gas film thickness reaching the physically relevant length scale.

1 Introduction

Binary droplet collisions play an important role in nature, science and technology, ranging from cloud dynamic processes to technical spray applications. The collision can result in coalescence or fragmentation of the droplets. The collision outcome strongly affects the resulting droplet size distributions. Consequently, any sound description of droplet population dynamics is to be based on a fundamental understanding of (binary) droplet collisions as an elementary process. The main question here is which specific outcome a binary collision takes, depending on the specific dimensionless parameters of the process. Figure 1 displays some prototypical collision processes. Binary collision outcomes for identical droplets are mainly influenced by the Weber number (We) and the impact parameter (B; see Fig. 2) with

$$\begin{aligned} We=\frac{\rho _{\textrm{l}}U_{\textrm{rel}}^2 D_0}{\sigma }\quad \text {and}\quad B=\frac{b}{D_0},\end{aligned}$$

(1)

where $\rho _\textrm{l}$ is the liquid density, $U_{\textrm{rel}}$ the relative velocity, $D_0$ the initial droplet diameter, $\sigma$ the surface tension and b is the offset of the droplets’ trajectories. Further parameters of influence are the Ohnesorge number (Oh), as well as the ratios of density ($\psi$) and the viscosity ($\varphi$) with

$$\begin{aligned} Oh=\frac{\eta _{\textrm{l}}}{\sqrt{D_0\rho _{\textrm{l}}\sigma }},\quad \psi =\frac{\rho _{\textrm{l}}}{\rho _{\textrm{g}}}\quad \text {and}\quad \varphi =\frac{\eta _{\textrm{l}}}{\eta _{\textrm{g}}},\end{aligned}$$

(2)

where $\eta$ denotes the viscosity and the subscript g (l) refers to the ambient gas (liquid). The Ohnesorge number could be replaced by the Reynolds number $Re =\sqrt{We }/Oh$, but the use of Oh is advantageous since the collision outcome is then only dependent on We and B, as long as the material parameters and the droplet size are fixed. Another dimensionless parameter used for the description of the temporal evolution of a droplet collision is the dimensionless time $t^*=t U_{\textrm{rel}}/D_0$ with $t^*=0$ at the moment of initial contact.

Fig. 1

(reproduced from [34] with permission)

Prototypical collision processes (left to right, cf. Fig. 2 for the types I-VI): coalescence (I), bouncing (II), coalescence (III), near head-on separation (IV), off-center separation (V) (reproduced from [14] with permission, Copyright by Cambridge University Press) and spatter (VI)

Fig. 2

adapted from [28]

Regime map for binary droplet collision as a function of Weber number We and center offset $B=\frac{b}{D_0}$,

In previous work [14, 28, 31], mainly on droplet collisions at low Weber numbers, distinct regimes of collision outcomes were characterized and condensed into a collision regime map. An adapted map is shown in Fig. 2. In the regime of low Weber numbers, four different collision outcomes are identified depending on the combination of We and B. These are: coalescence (I, III), bouncing (II), near head-on separation (IV) and off-center separation (V); the latter two are also called reflexive separation and stretching separation, respectively. The numbers in brackets are the corresponding section numbers in the regime map. For larger We, the rim of the collision complex becomes unstable, fingers are formed and partly ejected as satellite droplets. This collision outcome is named spatter (VI) or shattering (also splashing, in particular for droplets colliding with a wall). Going from left to right in Fig. 1, we see examples for these different collision outcomes, from (I) to (VI). Highly relevant and current topics are binary collisions of drops of different liquids, either miscible or immiscible. Droplet collisions of different, miscible liquids appear for instance in spray drying processes, where different droplets’ flight history can lead to different liquid compositions [10]. This leads to Marangoni stress since the surface tension depends on the local composition [11]. This yields to changes of the boundaries in the collision regime map, cf. [1]. In the collision of droplets of immiscible liquids, which is important not only for pharmaceutical applications but also for modern combustion engines in which water is injected [3], the interfacial tensions determine the contact angle at the three-phase contact line. This can lead to a (partial) enclosure of the droplet of higher surface tension. Bouncing, coalescence, crossing and reflexive separation [3, 27] are possible outcomes of collisions of immiscible liquids. For partial wetting adhesive merging [40] and a retraction of the contact line after a full encapsulation [41] can occur additionally. Despite impressive experimental work such as, e.g., [14, 22, 23, 27, 28, 31, 42], a fundamental analysis of the local phenomena during (binary) droplet collisions requires full information of the instantaneous velocity and pressure fields within the liquid phase. This requires a computational approach, based on fully resolved, detailed numerical solution of the governing equations from continuum physics. This leads to deep insights into the local flow dynamics, a valuable input for model enhancement [26]. In the following, the results achieved in the sub-project A7 of the Collaborative Research Center SFB-TRR75 are presented.

2 Continuum Mechanical Model

The droplet collisions under consideration are modelled as incompressible (due to low Mach numbers), isothermal two-phase flows without phase change, where the two immiscible fluid phases are labelled by ± and occupy the time-dependent domains $\Omega ^\pm (t)$, such that $\Omega$ is disjointly decomposed into $\Omega =\Omega ^+(t)\cup \Omega ^-(t)\cup \Sigma (t)$, where $\Sigma (t)$ denotes the sharp interface. The velocity $\textbf{u}(t,\textbf{x})$ and pressure $p(t,\textbf{x})$ are governed by the Navier-Stokes equations

$$\begin{aligned} \rho \left( \partial _t\textbf{u}+\left( \textbf{u}\cdot \nabla \right) \textbf{u}\right) =-\nabla p+\eta \Delta \textbf{u}+\rho \textbf{g} \quad \text {and}\quad \nabla \cdot \textbf{u}=0,\end{aligned}$$

(3)

where $\textbf{g}$, $\rho$ and $\eta$ denote the acceleration due to gravity, the phase-specific mass density and dynamic viscosity, respectively. The interface is assumed to carry no mass and admits no tangential slip of the adjacent velocity fields $\textbf{u}$. With constant surface tension $\sigma$, one obtains the dynamic and kinematic jump conditions

$$\begin{aligned} \llbracket p\textbf{I}-\eta (\nabla \textbf{u}+\nabla \textbf{u}^\textsf{T})\rrbracket \textbf{n}_{\Sigma }=\sigma \kappa _\Sigma \textbf{n}_\Sigma \quad \text {and}\quad \llbracket \textbf{u}\rrbracket =\textbf{0},\end{aligned}$$

(4)

where $\textbf{n}_\Sigma$ denotes the unit normal and the so-called jump-bracket

$$\begin{aligned} \llbracket \phi \rrbracket (t,\textbf{x}):=\lim _{h\searrow 0}\left( \phi (t,\textbf{x}+h\textbf{n}_\Sigma )-\phi (t,\textbf{x}-h\textbf{n}_\Sigma )\right) \quad \text {for}\quad \textbf{x}\in \Sigma (t)\end{aligned}$$

(5)

expresses the jump of a phase-specific quantity $\phi$ across the interface $\Sigma$. Neglecting phase change, which is a valid assumption considering the small time-scales of droplet collision processes, the speed of normal displacement $V_\Sigma$ of the interface fulfils the kinematic condition $V_\Sigma =\langle \textbf{u},\textbf{n}_\Sigma \rangle$, i.e. the interface is passively advected by the flow.

3 Numerical Method

Due to the possible occurrence of topological changes during the droplet collision numerical methods with implicit treatment of the interface are advantageous. Since conservation of phase volume is also desired, the Volume of Fluid (VOF) method of [13] has been employed for direct numerical simulations within this project, using the in-house solver Free Surface 3D (FS3D), originally developed by [33]; cf. [7] for an overview. The volume indicator function

$$\begin{aligned} f(t,\textbf{x}):={\left\{ \begin{array}{ll}0&{}\text {for}\quad \textbf{x}\in \Omega ^+(t),\\ 1&{}\text {for}\quad \textbf{x}\in \Omega ^-(t),\end{array}\right. } \end{aligned}$$

(6)

encodes the assignment of a point $\textbf{x}$ to the respective phase, which is governed by

$$\begin{aligned} \partial _tf+\textbf{u}\cdot \nabla f=0\text {,}\end{aligned}$$

(7)

see Fig. 3 for illustration. One thus obtains a one-field formulation of the two-phase Navier-Stokes equations, namely

$$\begin{aligned} \rho \left( \partial _t\textbf{u}+\left( \textbf{u}\cdot \nabla \right) \textbf{u}\right) =-\nabla p+\eta \Delta \textbf{u}+\rho \textbf{g}+\textbf{f}_\Sigma . \end{aligned}$$

(8)

The volume specific forces due to surface tension $\textbf{f}_\Sigma$ enter the model either by using the continuous surface force (CSF) [2] or stress (CSS) [19] formulation, i.e.

$$\begin{aligned} \textbf{f}_{\Sigma }^{\textrm{CSF}}=\sigma \kappa _\Sigma \textbf{n}_\Sigma \delta _\Sigma \quad \text {or}\quad \textbf{f}_{\Sigma }^{\textrm{CSS}}=-\sigma \nabla \cdot (\delta _\Sigma (\textbf{I}-\textbf{n}_\Sigma \otimes \textbf{n}_\Sigma )),\end{aligned}$$

(9)

where $\delta _\Sigma$ denotes the surface Dirac distribution. The control volumes for the velocity and pressure, respectively, admit a staggered (MAC) arrangement in the Cartesian grid. Combining the pressure-projection of [4] with a first-order temporal discretization, the convective terms are semi-linearized by an alternating permutation of directionally split transport, which ensures second-order accuracy [39]. In order to obtain reliable results, it is crucial to maintain a sharp interface throughout the simulation. Thus, the passive advection of the volume fractions, i.e. the time-integration of Eq. (7), is carried out geometrically [32]. Since the accuracy of initial volume fractions crucially affects the simulation outcome, the algorithm of [16] is employed. The details of the interface reconstruction are deferred to Sect. 3.1. In the CSF formulation for the surface tension, cf. Eq. (9), the curvature $\kappa _\Sigma$ is obtained from height functions [29]. Figure 4 contains a flowchart of the procedure.

Fig. 3

Physical phases on a Cartesian grid represented by volume fractions f with Piecewise Linear Interface Calculation (PLIC, [32])

Fig. 4

Schematic flowchart of solver Free Surface 3D

3.1 Efficient Interface Reconstruction

In geometric VOF methods, the interface reconstruction $\Gamma$, which is an approximation of the true interface $\Sigma$ from volume fractions, constitutes a central element in terms of both accuracy of the method and consumption of computational resources. Since the physical phenomena under consideration feature highly dynamic deformations as well as a large amount of surface area, a robust and efficient method for the reconstruction of the interface from the volume fractions f is required. In the Piecewise Linear Interface Calculation (PLIC) approach, this task comprises the computation of the normal field and the positioning of a plane in each cell to match the associated volume fraction.

Fig. 5

Results from an experimental study of the wetting behaviour reported in [41]. This study was performed in a cooperation of sub-projects A7, B1 and C4 together with D. Baumgartner and C. Planchette from TU Graz. 73% Glycerol- water solution (G73) (transparent) colliding with an other droplet (blue) of Bromonaphtalene (B) as a partially wetting liquid or Silicon Oil M5 (SOM5) as a fully wetting liquid. From left to right: G73 + B just before collision; G73 + SOM5 rim formation; G73 + B disc formation; G73 + B retraction of the contact line after collision (dewetting). The collisions finally show separation into two or multiple droplets

Fig. 6

Three component topological configurations. From left to right: a anti-parallel, non-wetted interfaces, b parallel, fully-wetted interfaces, c independent non-wetted interfaces, d fully-wetted, layered interfaces, e interfaces with a contact line

As observed in experiments and shown in Fig. 5, the presence of a third phase enriches the topological spectrum by thin films and contact lines, which requires a robust and efficient method to capture all possible topologies. In analogy to a liquid interacting with a solid particle, for the three-fluid configurations under consideration, the liquid with the lower surface tension (denoted secondary, index 2) towards the gaseous face always covers the one with the higher surface tension (denoted primary, index 1). Thus, a sequential reconstruction can be justified, i.e. the primary normal $\textbf{n}_{\Gamma ,1}$ can be computed first and independently of the secondary normal $\textbf{n}_{\Gamma ,2}$, by resorting to the method of [43] employed for two-phase configurations. Figure 6 illustrates the topological configurations, each of which is associated to a different scheme for the computation of the secondary normal $\textbf{n}_{\Gamma ,2}$. While the normal $\textbf{n}_{\Gamma ,1}$ of the first component is obtained in the same way as in the two-phase case, the computation of $\textbf{n}_{\Gamma ,2}$ is performed in a case-sensitive manner as follows [30]:

• a. Anti-parallel interfaces with a thin gas film in between: $\textbf{n}_{\Gamma ,2} = -\textbf{n}_{\Gamma ,1}$

• b. Parallel interfaces, liquid film: $\textbf{n}_{\Gamma ,2} = \textbf{n}_{\Gamma ,1}$

• c. Non-wetted, independent interfaces (gas film): $\textbf{n}_{\Gamma ,2} = \nabla f_{2}$

• d. Fully wetted, layered interfaces (liquid film): $\textbf{n}_{\Gamma ,2} = \nabla (f_{1}+f_{2})$

• e. Interfaces form a contact line with contact angle $\theta$ (see [24]):

  $$\begin{aligned} \textbf{n}_{\Gamma ,2} = \textbf{n}_{\Gamma ,1} \cos (\theta ) + \frac{\nabla f_{2} - \left( \textbf{n}_{\Gamma ,1}\cdot \nabla f_{2}\right) \textbf{n}_{\Gamma ,1}}{||\nabla f_{2} - \left( \textbf{n}_{\Gamma ,1}\cdot \nabla f_2\right) \textbf{n}_{\Gamma ,1}||} \sin (\theta )\end{aligned}$$

• f. Two component cells with liquid 1 truncating the liquid 2’s reconstruction stencil: Treatment like three component cells (a-e).

A sequential positioning of the PLIC planes, which is described below, is performed such that the volume fractions are each enclosed. The correct orientation is determined in three steps:

• 1. The number of neighbouring cells are counted, in which an enclosed volume fraction as well as the true volume fraction $f_{2}$ are 1. The enclosed volume fraction results from an extension of the PLIC plane in cell (i, j, k) to its neighbours. The orientation with the maximum count is chosen.

• 2. If two orientations have the same count in step 1, the minimum of

  $$\begin{aligned} g_{i,j,k}(\textbf{n}_{\Gamma ,2}) = \sum \limits _{i^*=i-1}^{i+1}\sum \limits _{j^*=j-1}^{j+1}\sum \limits _{k^*=k-1}^{k+1} (f_{2,i^*,j^*,k^*}-f_{p,i^*,j^*,k^*}(\textbf{n}_{\Gamma ,2}))^2 \end{aligned}$$

  (10)

  is employed (similar to [25]) as a second criterion for the decision, which orientation yields the true topology. The predicted volume fraction $f_p$ is again computed in all neighbouring cells.

• 3. A topology check is performed. If the topological properties are not as assumed, the orientation is not chosen.

This method of computing the interfaces normals is explicit, thus relatively fast, and allows the reconstruction of thin liquid or gas films. Before Eq. (10) can be evaluated, the secondary PLIC interface must be positioned to match the associated volume fraction. With the discrete fields of approximative normals $\textbf{n}_{\Gamma ,i}$ at hand, positioning the PLIC planes in a cell translates to finding the unique root of a scalar monotonous function, namely the truncated cell volume. By exploiting the Gaussian divergence theorem in combination with a joint co-moving coordinate origin $\textbf{x}_0$, the truncated volume is conveniently parametrized in terms of the signed distance as a sum of face-based quantities. This is favourable in terms of general applicability and allows to restrain from the costly extraction of topological connectivity at runtime. By assigning to each face (index k) of the original cell an individual origin $\textbf{x}_{0,k}$ (co-moving with the intersection of the PLIC planes; cf. Fig. 6), the volume computation can be extended to the three-phase case, where the computational cells are truncated twice. The full mathematical details can be found in [18]. In both cases, the derivatives of the volume function, which are exploited in a higher-order root finding scheme, can be obtained at negligible cost. With the polyhedron truncation and volume computation being decisive for the computational effort, the developed approach was shown to be highly efficient. On average, only one to two truncations are required to position the plane, outperforming existing methods. Further information concerning the performance assessment for an extensive set of polyhedrons, volume fractions and normal orientations can be found in [17], along with comprehensive description of the algorithm. Ongoing research aims at embedding the developed volume computation into a minimization-based scheme to improve the estimation of the normal fields $\textbf{n}_\Gamma$ from volume fraction data f. The above described non-iterative topology-capturing algorithm for the orientation and the efficient positioning are the basis for simulations of immiscible liquids.

3.2 Surface Forces for the Collision of Immiscible Liquids

During the collision of immiscible liquids, thin films of the outer liquid 2 spread over the inner liquid 1 which has the higher surface tension $\sigma _{13}$ towards the continuous phase 3, cf. Fig. 5. The wetting behaviour in droplet collisions of two immiscible liquids 1 and 2 is governed by the interfacial tensions and is described by the spreading parameter $S = \sigma _{13}-\sigma _{23}-\sigma _{12}$ [27]. In case of $S > 0$, liquid 2 fully wets liquid 1, while partial wetting occurs if $S < 0$. In Sect. 3.1 it was already discussed that the reconstruction of various topologies of three-component situations are to be captured correctly in numerical simulations of collision processes of immiscible liquids. The modelling of the surface forces is equally affected, as the surface forces are modelled utilising the volume fractions’ gradients which are truncated for liquid 2. The CSS model, cf. Eq. (9), employs the interface normals

$$\begin{aligned} \tilde{\textbf{n}}_\Gamma = \frac{\nabla \tilde{f}}{||\nabla \tilde{f}||} \end{aligned}$$

(11)

and the approximation of the surface Dirac distribution

$$\begin{aligned} \tilde{\delta }_\Gamma =||\nabla \tilde{f}|| \text {.} \end{aligned}$$

(12)

The normals of the interfaces are approximated analogously to the two-phase PLIC algorithm, but applied to a quadratically smoothed volume fraction field $\tilde{f}$. Compared to calculations of the surface Dirac distributions from the PLIC planes’ areas, much smoother results are obtained by the approximation according to Eq. (12), as the discretisation effect is lower, cf. [19]. The acceleration due to the surface forces, i.e.

$$\begin{aligned} \left( \partial _t\textbf{u}\right) _\Gamma = \frac{2}{\rho _1 + \rho _3} \textbf{f}_\Gamma ^{\textrm{CSS}},\end{aligned}$$

(13)

results from a scaling with the average density at the interface, which again avoids an influence of the discretisation. The use of a local density proved insufficient in multiple studies with FS3D. A physical reasoning for this might be that the surface force is acting very locally. This becomes important for the extension of the CSS model to immiscible liquids, cf. [30]. This extension is facilitated by introducing partial surface tensions $\gamma _i$, following the idea of [15, 38]. The interfacial tensions between component i and j,

$$\begin{aligned} \sigma _{ij} = \gamma _i + \gamma _j, \end{aligned}$$

(14)

are described as sums of the partial surface tensions $\gamma _i$. For the further discussion, index 1 denotes the inner liquid, 2 the covering liquid and 3 the gaseous continuous phase. For the superposition of the interfaces’ contributions, the capillary pressure tensor, cf. [19],

$$\begin{aligned} \textbf{T}_i = -\nabla \cdot (\tilde{\delta }_\Gamma (\textbf{I}-\tilde{\textbf{n}}_\Gamma \otimes \tilde{\textbf{n}}_\Gamma )) \end{aligned}$$

(15)

is computed for each component to retrieve the resulting acceleration

$$\begin{aligned} \left( \partial _t \textbf{u}\right) _\Gamma =\frac{1}{\rho } \sum \limits _{i=1}^{3}\gamma _i \textbf{T}_i \text {.} \end{aligned}$$

(16)

A virtual volume fraction field of liquid 2 is constructed in the vicinity of three component cells to obtain a sufficient stencil for the smoothing of the volume fractions. It is extrapolated from the PLIC interfaces at the contact line. The orientations $\tilde{\textbf{n}}_{\Gamma ,i}$ and the approximated surface Dirac distributions $\tilde{\delta }_{\Gamma ,i}$ are computed with Eqs. (11) and (12) for the respective phase. The information of the interfaces present is necessary for the choice of the scaling density

$$\begin{aligned} \rho = {\left\{ \begin{array}{ll} 0.5(\rho _1 + \rho _3)&{}\text {in 3 component cells in the smoothed field},\\ 0.5(\rho _i+\rho _j)&{}\text {in 2 component cell with smoothed}\, i-j\text {-interface},\end{array}\right. } \end{aligned}$$

(17)

employed to calculate the accelerations. As the topology information is not available inside the smoothed field, the information from the unsmoothed field is used to identify cells with a thin film. For all other cells, the information of the smoothed volume fractions is used to identify the interfaces present. This leads to a non-smoothing approach in cells with a thin film and is identical to the original CSS model everywhere else. Like for the reconstruction, additional modelling of the surface forces at thin liquid films is necessary. In cells with a thin film of liquid 2 covering liquid 1 and cells with a contact line, which is surrounded by such liquid film cells, the accelerations are chosen as

$$\begin{aligned} \left( \partial _t \textbf{u}\right) _{\Gamma ,\textrm{film}} =\frac{2}{\rho _3+\rho _1} (\gamma _1 + 2\gamma _2 + \gamma _3 ) \textbf{T}_1,\end{aligned}$$

(18)

assuming that the interfaces are parallel to each other and the film is very thin. The interface of the covering liquid reduces the interfacial tension compared to the two component case. The thin film does not act like a liquid-liquid interface yet. Therefore, the surface force scales with the inner liquid’s density jump towards the gaseous phase 3. This choice of the density strongly affects the collision outcome of two immiscible liquid droplets. The choice of the density in three component cells according to Eq. (17) is motivated by the results presented in Sect. 4. The analysis for the collision of fully wetting liquids has shown that the modelling of the surface forces is very delicate. Simulations of the partial wetting behaviour qualitatively yield the expected behaviour with this choice of the scaling density. However, for partial wetting, no suitable experimental validation data could be found in the literature.

Fig. 7

Simulation of the interaction of 50% glycerol solution (dark) and the silicon oil M5 (light) compared to experimental results [27]. Top: Coalescence: $D_{1}=187~\mathrm {\upmu m}$, $B=0.089$, $U_{\textrm{rel}}=1.98~\mathrm {m/s}$; Bottom: Crossing separation: $D_{1}$ = $189~\mathrm {\upmu m}$, B = 0.014, $U_{\textrm{rel}}$ = $3.71~\mathrm {m/s}$. Both: $D_{2}$ = $195~\mathrm {\upmu m}$, $6.25~\upmu \mathrm {m/cell}$. Pictures of experiments by courtesy of C. Planchette

4 Droplet Collisions of Immiscible Liquids

In Sects. 3.1 and 3.2, a topology-capturing PLIC method as well as a film stabilisation CSS method were presented. With those new methods, numerical simulations of droplets of two immiscible liquids were enabled in FS3D. A comparison to experiments of near head-on coalescence as well as crossing separation of equally sized droplets of silicon oil M5 and a 50% glycerol-water solution are shown in Fig. 7. The visualisation of the PLIC interfaces was enabled by a Paraview plugin developed in sub-project A1, see Heinemann et al. in this volume. Gravitation is omitted in the simulations due to limitations of the computational domain size at an adequate resolution. Both cases shown are near head-on collision, but no symmetry is imposed in the algorithm. The film stabilisation in both the reconstruction as well as the CSS model enable the simulation of fully wetting liquids, where a thin film of the outer liquid covers the inner liquid. The choice of the density, which scales the interfacial forces, is essential to capture the collision outcomes. Also the smoothing has to be done carefully: In the simulations shown in Fig. 7, one smoothing step was sufficient. With half the resolution, no smoothing at all showed the best results, but at a higher resolution, the number of smoothing steps had to be increased to counteract the refinement. A detailed analysis of this preliminary finding is necessary with further increase of the resolution enabled by parallelisation of the algorithm. The successful comparison to experimental data shows that the methods utilised for the simulation of immiscible droplet collisions are valid. This enables future studies on the influence of the geometry parameters and fluid properties on the collision outcome of fully wetting immiscible liquid droplet collisions. According to [41], the influence of the wetting behaviour plays an important role for the collision morphology and outcome. Further investigations of the contact line modelling are therefore necessary for partially wetting droplets. This is planned as soon as sufficient experimental validation data is available across different collision regimes.

5 Bouncing Versus Coalescence

Under standard ambient conditions, the regime transition from bouncing to coalescence (II$\rightarrow$III, cf. Fig. 2) occurs around the critical Weber number $We _{\textrm{crit}}\approx 13.63$ [28]. The critical distance between the colliding droplets is in the order of 10 nm and the mean free path of gas molecules $\lambda$ is about 100 nm.

Fig. 8

Comparison of experiments (gray, reproduced from [23] with permission, Copyright by AIP Publishing) and direct numerical simulation (blue, [21]) for bouncing (top row, $We =9.33$, $Re =110.36$) and coalescence (bottom row, $We =13.63$, $Re =134.24$) of binary droplets for varying time instances

Since the head-on collision process can be considered symmetrical, cf. Fig. 8, the computational effort can be significantly reduced by numerically replacing one of the droplets by a mirror droplet via a collision plane. With gas film thicknesses $h\approx \lambda$ to be expected, however, a full spatial resolution of the problem remains infeasible due to the viscous time step limit. In order to correctly account for the effect of the thin gas film, a subgrid scale model was applied at the collison plane. Note that, while the thickness of the gas film between the droplets becomes very small during the collision process, the lateral film extension permits the application of a continuum model. Indeed, the so-called lubrication approximation is valid and was shown to accurately capture the relevant physics [12]. For a given film thickness h, the pressure p is governed by

$$\begin{aligned} \nabla \cdot \left( {\left( {\frac{h^3}{3}+C_1\lambda h^2+C_2\lambda ^2h}\right) \nabla p}\right) =\eta \left( {\partial _th+\nabla \cdot \left( {h\textbf{u}^\Sigma }\right) }\right) ,\end{aligned}$$

(19)

where $\textbf{u}^\Sigma =[u^\Sigma ,v^\Sigma ]^\textsf{T}$ denotes the lateral interface velocity and the coefficients $C_i$ account for the slip of the lateral velocity at the collision plane [6]. At the lateral boundary of $\Omega _{\textrm{gap}}$, Dirichlet boundary conditions for the pressure are applied. The local film thickness h is obtained from the volume fractions in the three cell layers above the collision plane, and determines the dynamic extension of the lubrication domain. E.g., for a collision plane with normal $\textbf{e}_z$ and base point $(0,0,z_{\textrm{lub}})$, one obtains

$$\begin{aligned} \Omega _{\textrm{gap}}(t;z_{\textrm{lub}}):=\{(x,y)\in \mathbb {R}^2:(x,y,z_{\textrm{lub}})\in \Omega \,\text {and}\,h(t;x,y)\le h_{\textrm{lub}}\},\end{aligned}$$

(20)

where $h_{\textrm{lub}}$ is chosen to cover two to three layers of cells. The numerical investigations were conducted resorting to two realisations of the collision plane.

5.1 Collision with a Symmetry Plane

Aiming at the numerical prediction of transition, we have conducted numerical investigations for a series of Weber numbers between 7.77 and 18.8 [21]. As stated above, the binary collision numerically corresponds to a single droplet colliding with the domain boundary, where symmetry conditions were applied for the velocity. This approach was found to be very well suited for capturing the physics for prescribed collision outcome, cf. again Fig. 8. The procedure works as follows: once the film thickness h reaches a certain threshold ($\approx 3$ cells), the surface forces exerted by the mirrored interface induce an artifical rupture of the gas film, yielding coalescence for all We. If the symmetry of the volume fraction field is replaced by a Dirichlet condition for the volume fractions, for analogous reasons, one always obtains only bouncing. Thus, switching between these two strategies at some user-defined time (if the gas film is sufficiently thin), allows to compute the desired outcome. For an unspecified collision outcome, however, no convergence with mesh refinement could be obtained. Above some threshold resolution, only bouncing was observed in the numerical simulations of configurations for which experiments show coalescence. It was found that the minimum gap height ($\approx 200$ nm) did not reach physically reasonable levels ($\approx 70$ nm), thus inhibiting coalescence. From a follow-up series of numerical experiments and comparison to [5], it was conjectured that this numerical artifact is caused by imposing symmetry boundary conditions at the collision plane, suggesting to shift the collision plane into the interior of the computational domain.

5.2 An Artificial Interior Collision Plane

For the collision with the domain boundary described above, the lubrication pressure was coupled to the flow solver in two steps: after the original time stepping scheme shown in Fig. 4 is performed, the bulk pressure at the boudary of $\Omega _{\textrm{gap}}$ serves as a Dirchlet boundary condition for Eq. (19). Since $\Omega _{\textrm{gap}}\subset \partial \Omega$, the resulting lubrication pressure is then in turn applied as a boundary condition for the pressure projection in the bulk, such that applying the resulting gradient ensures a divergence-free velocity. An internal collision plane, i.e. $\Omega _{\textrm{gap}}\not \subset \partial \Omega$ (left panel in Fig. 10), however, requires two modifications: firstly, the lubrication pressure cannot be prescribed as boundary condition for the bulk pressure. Instead, the lubrication pressure exerts a surface force (analogous to the CSF-model for the surface tension, cf. Sect. 3). Secondly, in preliminary numerical experiments, it was found that employing an intermediate bulk pressure as a local boundary condition for the lubrication equations induces numerical noise. In order to alleviate this influence, by taking into account the rotational symmetry, the local values are replaced by an average, say $p_{\textrm{BC}}$, over the boundary of $\Omega _{\textrm{gap}}$. Equation (19) translates into a linear system whose solution is computed by a standard GMRES scheme [35]. The drawing in Fig. 10 (top-left) illustrates that, as a result of the Cartesian domain decomposition of the flow solver FS3D, the collision plane only intersects a subset of the parallel processes. Furthermore, the domain $\Omega _{\textrm{gap}}$, within which the lubrication equation is solved, dynamically deforms along with the flattening of the droplet, implying potentially large spatial load imbalances. For an efficient solution, a load balancing scheme for deformable domain was designed and implemented. Figure 9 illustrates the concept. In comparison to the Gauss-Seidel-type approach of Sect. 5.1, the computational effort was reduced significantly.

Fig. 9

Initial distribution of degrees of freedom (left, involving only 4 of 12 processors), load balancing scheme and balanced distribution (involving all processors) (illustrated for the configuration of Fig. 10, where the colors correspond to the individual processors)

Fig. 10

Illustration of interior collision plane with processors (left) and minimum film thickness over time (right; $We =3.47$; FEM: courtesy of M. Chubynsky) with vertical lines corresponding to the instants shown in the bottom row

While the droplet deformation in Fig. 10 qualitatively coincides with the top row in Fig. 8 (collision with domain boundary), the modified subgrid-scale model allows to additionally obtain physically reasonable minimum film thicknesses for an appropriate choice of $p_{\textrm{BC}}$, even for a very coarse spatial resolution; see the diagram in Fig. 10 (top-right). However, recall that the boundary condition $p_{\textrm{BC}}$ resorts to the bulk pressure, whose absolute value does not necessarily carry any physical meaning for incompressible flows, since only the gradient is employed for the projection; cf. Fig. 4. Contrary, the lubrication pressure enters the model directly. Thus, any offset strongly influences the simulation outcome and a consistent intrinsic choice of $p_{\textrm{BC}}$ is subject of ongoing research.

6 Spattering Collisions at High Weber Numbers

In high energetic head-on collisions of two droplets with the same diameters and material properties ($We \ge \mathcal {O}(100)$), the rim of the collision complex formed on the collision plane becomes increasingly unstable with increasing Weber numbers [34]. If the Weber number is high enough, it yields substantial disintegration into secondary droplets. In order to prevent the thin liquid lamella emerging at high energy collisions from unphysical rapture, the stabilisation method of [9] has been further improved [20]. With this approach, interface tension forces are correctly calculated also at liquid interfaces with a neighbouring, additional interface. Since the numerical simulation of binary droplet collisions at high Weber numbers requires a high grid resolution, a domain adjustment technique based on [8] can be used to reduce the computational effort. With a small white noise disturbance (1% of the initial relative velocity of the droplets) exerted to the initial velocity field, simulation results are in excellent agreement with the corresponding experiments for $We =443$ [20]. A comparison of the results is reproduced in Fig. 11. These high-resolution calculations offer detailed insights into areas of the collision complex that are not accessible with experimental methods. Furthermore, the comparison of a Fast-Fourier Transform analysis (FFT) combined with statistical techniques applied to the VOF simulation’s results with theoretical predictions based on [44] revealed that the rim instability in binary droplet collisions is completely dominated by the Plateau-Rayleigh instability [20].

Fig. 11

Comparison of experiment (top, [22]. Reprinted with permission, Copyright by the American Physical Society) and direct numerical simulation (bottom, [20]) for head-on collision of binary droplet at $We =443$ and $Re =6207.3$

Simulations of binary droplet collisions at still higher Weber number ($We = 805$ in the following) require even higher grid resolutions. Grid independence studies need to be performed to determine the necessary resolution. The white noise signal imposed on the velocity field as initial disturbance, however, is grid-dependent as higher frequencies are resolved on a finer grid. Therefore, a function for the description of the initial disturbance which can be transferred to equidistant grids of different resolution is constructed. This function should have the non-regular structure of a white or coloured noise signal and consider all spatial directions equally. In [36], a first idea of this function has been presented and is given here in more detail. The noise function, written exemplarily for the first component,

$$\begin{aligned} u_{noise }\left( \tilde{x}_1, \tilde{x}_2, \tilde{x}_3 \right) = A \sum _{i=0}^{I_max /3}\sum _{k=0}^5 \prod _{l=1}^{3} \cos \left( \left( 3i+\pi _{k,l}\right) \Phi \tilde{x}_l\left( 1+aR_{6i+k}\right) +\omega _{6i+k} \right) \end{aligned}$$

(21)

is imposed on each component of the velocity field, where $\pi _{k,l}$ refers to the l-th component of the k-th permutation of (1,2,3). The spatial coordinates are normalised by multiplication with $2\pi$ and division by the maximum length of the computational domain. The imposed disturbance is normalised by $A:=\frac{1}{2}||\textbf{u}_{\textrm{noise}}||$. Due to the multiplication with a random number $R_k\in [-1,1]$, the summands are incommensurable with probability 1; thus, a regular structure of the noise signal is avoided. The factor $a\in ]0,1[$ prevents the suppression of individual frequencies, for example in case of $R_k=-1$. The addition with random numbers $\omega _k\in [0,2\pi ]$ in the argument of the cosine function prevents a local amplification of the signal. The factor $\Phi \in ]0,1]$ allows the representation of intermediate frequencies. The maximum frequency that can be represented is determined by the cell width. According to sampling theory, to reproduce a signal with the maximum frequency $f_{max }$, the sampling frequency must be more than twice as high [37]. Thus, the maximum representable frequency in the normalised domain $f_{max }\le N/2$, where N is the number of grid cells in the direction of the largest extension of the computational domain. The largest possible frequency in the argument of the cosine function is $\bar{f}_{max }=(I_{max }+2)\Phi (1+a)$. The upper limit of the sum (21) can thus be expressed as

$$\begin{aligned} I_{max } \le 3\left\lfloor \frac{1}{3}\left( \frac{N}{2\Phi \left( 1+a\right) }-2\right) \right\rfloor . \end{aligned}$$

(22)

To evaluate the $6 I_{max}/3$ summands per spatial direction, $12I_{max }$ random numbers are required. To validate the noise function, the comparison of the frequency spectrum of $u_{noise }\left( \tilde{x}_1, \tilde{x}_2, \tilde{x}_3 \right)$ with the spectrum of the previously used perturbation on a 256$\times$256$\times$64 grid is shown in Fig. 12 (left). The signal $u_noise$ is normalised here by the factor $u_{var}/\tilde{A}$, where $u_{var}$ is the variance of the original disturbance function, $\tilde{A}$ the arithmetically averaged amplitude of the disturbance function $u_{noise }$. Furthermore, the noise signal with a maximum frequency corresponding to the grid resolution 256$\times$256$\times$64 is transferred to disturb the initial velocity field on a grid with twice as many nodes in each direction (cf. Fig. 12 (right)). The agreement between the two disturbance signals is excellent.

Fig. 12

The noise signal $u_{noise }\left( \tilde{x}_1, \tilde{x}_2, \tilde{x}_3 \right)$ resolved on a 256$\times$256$\times$64 grid reveals a frequency spectrum which is comparable to the one of the white noise signal (left) and can be transferred to higher grid resolution, 512$\times$512$\times$128 (right); factor $a=0.95$, $\Phi =0.1$

A grid convergence study was performed for droplet collisions with $We =805$. Using symmetry conditions, the problem could be reduced to a quarter of a droplet. Its radius has been resolved by $0.\overline{11} N_x$ cells in a domain of $N_x\times N_x\times N_x/4$, where $N_x \in (512,1024,2048)$. However, the results do not yet show convergence. The collision complex initially spreads faster on grids with higher resolution. This finding indicates that there are additional physical effects on the smallest time scales that are not yet understood. One such effect might be the ageing of the interface. For further analysis we benefit from detailed insights into the evolution of the collision complex, the growth of liquid fingers on the rim and the detachment of secondary droplets. Figure 13 (left) shows the velocity field at time $t=0.54~ms$, where the first drops have already detached. Fluid is still flowing from the disc into the rim at high velocity. The detailed view of the rim region in Fig. 13 (center) illustrates that the viscous dissipation in the fluid is very high in a narrow transition region. The detailed view in Fig. 13 (right), which shows the pressure field in the area of the rim, confirms that fluid is flowing from the disc into the rim. The constrictions between two adjacent finger-forming structures are also characterised by a reducecd pressure. This causes them to grow together, whereas the high pressure in the area of the depicted single finger structure increases its growth in subsequent timesteps.

Fig. 13

Top view on spreading collision complex at $t=0.54~ms$. The rate of local viscous dissipation ($\mathrm {kg~m^2~s^{-3}}$, left) and the pressure difference to the ambience ($\mathrm {kg~m^{-1}~s^{-2}}$, right) are mapped onto the contour plot $f=0.5$

7 Conclusions

This project within the SFB-TRR75 was driven by the desire to understand small-scale phenomena in binary droplet collisions as a basis for new and improved multi-scale models of larger processes, which involve droplets under extreme conditions. The investigations’ concern were droplet collisions with immiscible liquids, the rim instability during spattering collisions at high Weber numbers and the regime transition between bouncing and coalescence in head-on collisions. In each of these different scenarios, the collision process is of multi-scale nature in itself: For collisions of immiscible liquids, triple lines and thin films of the two interacting liquids and ambient gas are present. The local forces in the vicinity of the triple line determine the overall collision outcome. Looking at collisions at high Weber numbers a liquid rim develops corrugations, which grow into liquid fingers and eventually pinch off to form satellite droplets. At the same time the rim is contracted by the surface tension forces of the extremely thin liquid lamella. Droplets approaching each other at low Weber numbers expel the thin gas film in the gap between the droplets. Coalescence can only occur if the Van-der-Waals forces come into action at gap thicknesses below 100 nm, a scale at which the flow in the gas film can be described by the rarefied flow equations and slip between gas and liquid becomes relevant. To allow for valid and accurate numerical simulations of such collision processes, the capabilities of the Volume-of-Fluid method were enhanced with new approaches like an efficient and topology-capturing PLIC algorithm for three immiscible volume fractions, an enhanced lamella detection and stabilisation algorithm for head-on collisions as well as a subgrid-scale model for computing the pressure in thin gas films based on lubrication approximation. However, due to the multi-scale character of the phenomena during droplet collisions, sufficient grid resolution remains a challenge: The collision outcome for immiscible droplets strongly depends on the modelling of the surface forces, where the mesh resolution and the choice of the scaling density proved to be important. However, the described methods enable a detailed calculation of droplet collisions of immiscible liquids also for asymmetric collisions. Also high energetic droplet collisions reveal a strong influence of the grid resolution. Direct Numerical Simulations were employed there to rigorously show that the Plateau-Rayleigh instability pattern dominates the rim instability - this is valuable for the development of analytical models for predicting the onset of spatter and the secondary droplets’ spectrum. Moreover, experimentally inaccessible data such as the local dissipation rates in the liquid throughout the collision process were obtained. Some challenges remain, caused by the multi-scale and multi-physical character of droplet collisions under extreme conditions. Predicting the transition between bouncing and coalescence, rather than computing the fine details retrospectively, still requires a more consistent coupling between the different solvers for the different governing equations, and a two-phase subgrid-scale modelling seems promising. A challenge of general nature is the simulation of flow instabilities, since the latter act as amplification systems, highly sensitive to disturbances in the fields. At this point, the introduction of carefully computed artificial noise is a key step enabling further research. Last, but not least, numerical simulations are restricted to the physics included in the underlying mathematical models. Here, we found indications for additional dissipative mechanics, being present in the experiments but not in the numerical simulations. Such discrepancies occur, if new surface area is generated at high local rates. This poses another extreme condition which deserves further investigation. To conclude, analysing binary droplet collisions on the small scales enhances the understanding and modelling applied on larger scales which enable predictions ranging from spray processes to weather forecast.