Computational Investigation of the Clustering of Droplets in Widening Pipe Geometries

Abstract

Experimentally, periodically released droplets in systems of widening pipes show clustering. This is surprising, as purely hydrodynamic interactions are repulsive so that agglomeration should be prevented. In the main part of this paper, we investigate the clustering of droplets under the influence of phenomenological hydrostatic forces and some hypothetical attraction. In two appendices, we explain why a direct numerical simulation for this system is rather more difficult (and probably not possible with current methods) than the “simple” geometry would suggest.

1 Introduction

The purpose of this work is the search for a theoretical understanding of the clustering of oil droplets (density 1.2 g/cm$^3$) suspended in water (density 1.0 g/cm$^3$) in an opening pipe as in the experiments by Li and Barrow [7]. Intuitively, one would expect that under the influence of purely hydrodynamic interactions, there should only dominate repulsion between droplets and no clustering should occur. We investigate the basic properties of the interaction between droplets, together with the effect of the simulation geometry with the phenomenological modeling of hydrodynamic forces and interactions. We desist from trying to implement a full-fledged direct numerical simulation (DNS) of droplets in flow, as the state of the art of direct numerical simulation is still far away from modelling particles, let alone droplets, in three dimensional fluid geometries with the exact boundary conditions. Some problems of the standard discretizations for DNS methods with freely moving bodies inside are explained in Appendix A, while the issues of “approximate” boundaries are discussed in Appendix B, to make it plain that we have selected our modeling approach from a higher insight of simulation methods, not due to ignorance of more elaborate methods.

Fig. 1.

Cylindrical pipe of the model system (above) with dimensions in meter and Hagen-Poiseuille flow profile (below) with the widening of the pipe to an “O”-shape at the right of $x=0,$ where the fanning out of the flow is computed based on the assumption of continuity, respectively incompressibility.

2 Our Modeling Approach

For simplification, we treat the experiment by Li and Barrow [7] as a widening pipe (see Fig. 1 above) without the butterfly-pipe mechanism which triggers the separation and release of the droplets.

Fig. 2.

Flow entering a pipe from a wider vessel with constant flow profile so that only after the “entrance length”, a trule parabolic Hagen-Poiseuille flow profile is reached.

2.1 Forces on Droplets in Pipes

We want to model small oil droplets with a radius of about 0.6 mm in water, so the surface tension should be large enough so that the deformation can be neglected and the droplets can be treated as rigid. We work with an experimental inside flowrate of 3.8 [ml/h]. Assuming a Hagen-Poiseuille flow with a parabolic profile, the maximal velocity in the center of the narrow inflow pipe (d=1.32 mm) will be $u_{in}=0.76$ [mm/s], while at the wider O-shaped outflow pipe it will be $u_{out}=0.34$ [mm/s], see Fig. 1. The Reynolds number Re based on the viscosity of water $\mu =10^{-3}$ [Pa s] with unit density $\rho$ and the total pipe diameter at inflow of 1.3 mm is

$$ Re_{max}=\frac{\rho u_{in} D_{\text {H}}}{\mu } \approx 4.4 \cdot 10^{-4}, $$

“deep” in the viscous regime of the Stokes flow ($Re<1$), and far away for inertia dominated ($Re>50$), let alone fully developed turbulent flow ($Re>1000$), so that vortices and turbulence should not play any role. For pipe flow, at inflows from a constant flow profile, it takes a certain distance, the entry (or entrance) length, until the flow profile changes to full developed Hagen-Poiseuille flow, see Fig. 2. For laminar flow, is typically taken “within 5%” [15] of the channel diameter, scaled by the Reynolds number, so for our case at the narrow inflow we have

$$ L_{h,n}= 0.05 \cdot Re D=5.2 \cdot 10^{-7} \text{[ }m], $$

and at the wider outflow, $L_{h,w}=1.460 \cdot 10^{-6} \text{[ }m],$ still considerably less than a droplet diameter. This means that for our setup, the flow assumes a parabolic profile practically instantaneously, even between droplets. The entrance length is a confirmation that viscous forces dominate and that the assumption of Hagen-Poiseuille flow is justified. (The assumption of a constant flow speed would reduce the Reynold number and therefore the entrance length, which would again confirm the prevalence of Hagen-Poiseuille flow.) We neglect any “elastic” deformations of the droplet which may lead to deformations of the shape, so for the force on the droplet we use Stokes law

$$\begin{aligned} F_d = 6 \pi \mu R v, \end{aligned}$$

(1)

Fig. 3.

(Dimensionless) wall correction factor and added mass over the reduced diameter from 0 to 1.

where v is the velocity difference between the center of mass of the droplet and the velocity of the flow in the pipe. The Stokes law is only valid for walls far away. For more narrow geometries, a (dimensionless) wall correction factor $F^*$ must be included, which takes the influence of the walls into account. Fitting the experimental data from Iwaoka and Ishii [5] for the dimensionless

$$\begin{aligned} r^* = \frac{r_{particle}}{r_{pipe}}, \end{aligned}$$

(2)

we found that at least for the experimental data range ($0 \le r^* \le 0.9$)

$$\begin{aligned} F^* = \frac{1}{(1-r^*)^{2.45}}. \end{aligned}$$

(3)

was a good parametrization, see Fig. 3. Another issue of the Stokes law is that it is valid only for equilibrium velocities. When the spheres undergo an acceleration in a dense liquid (with a density comparable to the density of the body), there is an added (or vitual) mass of the fluid around the particle which must be accelerated, too. [14]. For a sphere, that gives an additional inertia term, which in a cylindrical vessel depends on the vessel radius. For the boundaries at infinity, the correction is 0.5, i.e. for the acceleration of a spherical mass of volume V,  an additional mass of $0.5 V \rho$ must be taken into account. The correction coefficients for boundaries in finite distance have been computed theoretically by Smythe [13] and experimentally verified by Mellsen [8]. We have fitted the resulting graph to a (dimensionless) prefactor of

$$\begin{aligned} m^*=\frac{0.5}{(1-0.9243{r^*}^2)^{1.2579}}. \end{aligned}$$

(4)

The comparison to the wall correction factor can be seen in Fig. 3. Obviously, the added mass effect (for the accelerated motion) is considerably smaller than the wall correction factor for the motion at given velocity.

2.2 The Interaction of Droplets

It is possible to expand the repulsive force between two droplets which are approaching each other with velocity v as a correction to the Stokes law [2, 3]. For our simulation, we have parameterized the resulting prefactor of the interaction law for equal sized spheres, see Fig. 2a) of Goddard et al. [3]. This means that for distances $d>10r$ with r being the particle radius, the repulsive force decays so much that it corresponds to Stokes’ law alone while for distances $d< 0.1r$, the repulsive force increases as a prefactor to Stokes’s law with a 1/d dependence, while for the intermediate range $10r>d>0.1r$, there is a smooth transition between the repulsive (short range) and the neutral (long range) regime.

Fig. 4.

Release of the droplets with only hydrodynamic interactions and initial disorder.

3 Results: Gravity in X-Direction

For this preliminary investigation, we make the system more symmetric by setting the gravitation in x-direction: We first want to understand the effect of the interactions we have implemented without interference of an asymmetric effect from gravity between an upper and lower boundary. We integrate the equations of motion with MATLAB’s adaptive time ode15s-integrator and enforce the use of the BDF2 (Backward differentiation formula of second order) with a relative error tolerance of 0.1% and absolute errors of 1% of the droplet radius and 1% of the maximal velocity in the narrow channel. For this preliminary investigation, we work with a constant augmented mass of 50% only, without the corrections due to the closeness of the boundary. To take into account that the release of the droplet at the inlet will be affected by an oscillation of the droplet which may shift its center of mass, we add $\pm 3.75$ % equally distributed randomness. Our hope was that this disorder in the initial coordinates would trigger the clustering. The result for hydrodynamic forces only was disappointing and is shown in Fig. 4. The droplets are released like pearls on a string. While the initial disorder is amplified by the hydrodynamic interaction, so that the scattering of the vertical position around the symmetry axis is enhanced, no clustering occurs.

Fig. 5.

Force between the droplets from Eq. (5) drawn for $k_{mult}=1$ with the extension of the droplets.

3.1 Introducing an Additional Interaction

We next introduce an additional interaction. From comparison with the clustering in the experiment, we have to assume that there must be an additional attraction. On the other hand, once the droplets are in contact, they feel a kind of “elastic” repulsion, because they do not fuse. We construct the interaction for the droplet radius $r_{dr}$ so that

$$\begin{aligned} F^{int}=k_{mult} \left\{ \begin{array}{ll} \displaystyle -2.5 \cdot 10^{-4} \frac{r_{dr}}{r^2} &{} \text{ for } r>2.15 \, r_{dr}, \\ \displaystyle +10^{-4} \frac{(2 r_{dr}-r)}{r_{dr}}+ 10^{-5} \exp \left( 1.2 \frac{1.05 r_{dr}}{r-2r_{dr}} \right) &{} \text{ else, } \end{array} \right. \end{aligned}$$

(5)

so the upper term is a Coulomb-type attraction, while the lower term is a linear repulsion with an added exponential term to guarantee that no fusion (penetration) of droplets can occur. The Coulomb-type has been chosen because $1/r^2$ is the only algebraic form which allows to model a spatially decaying interaction between agglomerations by their centers of mass. While the coefficients in the interaction look a bit strange, they have been chosen so that repulsion changes continuously into attraction when the droplets come into contact, see Fig. 5.

Fig. 6.

Wordlines (x-coordinates over time) above and (three-dimensional) clusters for the final timestep of the worldlines below for $k_{mult}=0.8$. At $t=107$[s] (not shown), the first two clusters fuse to a seven-droplet cluster, but the three-droplet clusters behind seem to be stable.

3.2 Clustering for Different Values of $k_{mult}$

With attraction, the dynamics in the simulation develops two timescales: The slow timescale of the advection of the droplets through the pipe and the fast timescale due to the reordering of the droplet positions into clusters at close range under the influence of the attractive force. Depending on the prefactor $k_{mult},$ we obtain clusters with a different number of droplets. In Fig. 6, we have shown the result for $k_{mult}=0.8$ with three droplets per clusters in “equilibrium” and only the first cluster with four droplets. The release of three-droplet cluster was sustained also beyond the simulation time shown here. Interestingly, the number of droplets in a cluster is not proportional to the attraction strength, because the force equilibrium between attraction and hydrostatic repulsion is rather subtle: For $k_{mult}=0.4,$ a leading cluster of three droplets develops, which is then followed by pairs of droplets which over time fuse into clusters of four droplets, see Fig. 7. While the first four-droplet cluster fuses with the leading three-droplet cluster, the later four-droplet clusters stay intact when they are flushed along the system. Typically, the number of droplets in a cluster for a given value of $k_{mult}$ is constant, except the first cluster, which may have one droplet less or more and which then fuses with the second cluster.

Fig. 7.

Wordlines (x-coordinates over time) above and (three-dimensional) clusters for the final timestep of the worldlines below for $k_{mult}=0.4$. At the widening outlet, first two-droplet clusters form which then fuse into four-droplet clusters.

4 Summary, Conclusions and Future Work

For the current state of the art of fluid simulation techniques, direct numerical simulations of (hydrophobe) droplets in (hydrophile) fluids seem to be unfeasible. Therefore, to understand the experimentally observed clustering of droplets, we have proposed a phenomenological interaction of the droplets based on known hydrodynamic interactions with an additional Coulomb-type attraction. The phenomenology of this model is rather rich, as the number of droplets in a cluster is not proportional to the interaction strength: Due to a subtle balance of attraction and hydrodynamic repulsion, there can be clusters with more droplets for weaker attraction. In the next step, a better physical understanding of the possible attractive interactions and a matching with the experimental data have to be obtained.

Appendices

Appendix A: Problem of Direct Numerical Simulation

Ideally, one would investigate the flow of spherical droplets in another fluid with a “direct numerical simulation” of two-phase flow, with a geometry of droplets (themselves a continuum) inside a continuum of the outer fluids, with a suitable boundary in between. Unfortunately, such an approach is currently not feasible, neither from the computational nor from the algorithmic aspect. The success of the Finite Element Methods in structural dynamics (that means, problems which are predominantly linear) is so much taken for granted that in engineering tasks, there is a downright reflex to use (or at least recommend) Finite Elements. For flow problems of particles (droplets included) which come closer than several mesh sizes while additional refinement is not possible, there are indeed computational issues [9]. Already for two dimensions, the simulation of polygonal particles in fluids is rather difficult and expensive if the boundary conditions have to be taken into account exactly, see [9] for a suitable simulation method and [10, 11] for the resulting problems which have to be dealt with for the grid generation. If the boundary conditions are approximated, there is a danger of numerical instabilities which is treated in Appendix B. Curved boundaries of freely floating particles in fluids add an additional layer of computational complexity to the problem and have not yet been implemented according to our knowledge.

There are various packages which offer the option to simulate multi-phase flows via the Finite Volume method (FVM). In particular, the most striking advantage of the FVM is that it implicitly fulfils the incompressiblity condition, which for other continuum simulation methods must be ponderously enforced. Nevertheless, there are other drawbacks: First of all, the discretization into finite volumes makes a uniform connectivity for the underlying mesh necessary. As a result, in screenshots of simulations, small “steps” can be seen which represent the elementary volumes of the underlying uniformly connected grid. This step-like outlines leads to the problems with the zero-order approximation of boundaries which is discussed in Appendix B. A second drawback is that the non-linear Navier-stokes equations are solved in a linearized form: This gives a fast update speed, but leaves a remaining uncertainty about the validity of the result.

The least suitable method for two-phase flows is the Finite Difference Method: The finite-difference approximation of differential operators demands an underlying square grid, which severely limits the accuracy of the simulation. Worse, in spite of various approaches which claim to “interpolate” boundaries with “immersed”boundaries (interpolated and therefore of rather dubious validity), one can show the congruence of finite-difference methods with zero-order Finite Element methods [4]: This means that the boundaries in finite difference methods are always approximated in zero order, there are no boundary conditions “on” the boundary, just in the middle of the “element” next to the boundary. The issue of zero-order boundaries and the very likely occurence of numerical instabilities is treated in the next Appendix.

Appendix B: Little Known Hitches in the Approximation of Boundary Conditions to Zero Order

B. P. Leonard treated the snags in the application of low order finite difference methods in his highly readable, enligthening and amusing treatise “A survey of finite differences of opinion on numerical muddling of the incomprehensible defective confusion equation” [6]. The title refers to the confusion in the discussion of the error order when diffusion equations are augmented with advective terms, as well as the “remedies” which make the numerical solutions more “stable” but physically implausible. Nevertheless, Leonard’s treatise only focused on the discretization of the equations. There are additional issues for the treatment of boundaries, which affect in particular finite difference and finite volume methods, when smooth slanted or curved surfaces are implemented by a “stepped” profile outline. For a partial differential equation

$$\begin{aligned} \frac{d f(t)}{d t } = \hat{L} f(t), \end{aligned}$$

(6)

(with $\hat{L}$ as the sum of the spatial differential operators) the product

$$\begin{aligned} \tilde{f}(t) =\underbrace{ \exp \left( t (\hat{L} ) \right) }_{\text{ time } \text{ evolution } \text{ operator }} \underbrace{ \psi _0 }_{\text{ initial } \text{ condition }} \end{aligned}$$

(7)

is the formal solution (including additional spatial dependencies). When we solve the problem numerically, so that there is an additional error $\epsilon _L$ in the discretization of the spatial differential operator, our solution becomes

$$\begin{aligned} \tilde{f} = \exp \left( t (\hat{L} + \epsilon _L) \right) \psi _0. \end{aligned}$$

(8)

The solution of a numerical simulation is “stable” if we don’t generate any exponentially diverging noise via $\epsilon _L.$ If we add additional noise $\epsilon _0$ to the initial condition by not using the exact, but stepped we can write the time evolution formally as

$$\begin{aligned} \tilde{f} &= & \exp \left( t (\hat{L} + \epsilon _L) \right) \left( \psi _0 + \epsilon _0 \right) \nonumber \\ &= & \exp \left( t (\hat{L} + \epsilon _L) \right) \psi _0 + \exp \left( t (\hat{L} + \epsilon _L) \right) , \epsilon _0 \end{aligned}$$

(9)

and we can hope that the additional exponential term will not add much. As a simple test case of a higher-order solution method with zero order boundary, we can try to solve the harmonic oscillator where we replace the trapeze-method with the rectangle midpoint method. For quadrature over intervalls of length $\varDelta x$, the error for the midpoint method is

$$\begin{aligned} \epsilon _{mid}=\frac{1}{24} \varDelta x^3 f'' (x_i), \end{aligned}$$

(10)

while surprisingly, the error of the trapeze rule is twice as large,

$$\begin{aligned} \epsilon _{tra}=\frac{1}{12} \varDelta x^3 f'' (x_i), \end{aligned}$$

(11)

(both errors from [1]), but which can be understood graphically, see Fig. 8.

Fig. 8.

Approximation of the integral of function f in the interval $x_i$ to $x_i+1$ with the trapeze rule, which underestimates (red area) the integral more than the midpoint rule overestimates it (error compensation of red and green area). (Color figure online)

Fig. 9.

Solution of a damped harmonic oscillator with the exact solution in black and the ghost oscillations generated by the midpoint time integrator.

We can now shift from spatial coordinates x to time coordinates t and create time integrators for differential equations

$$\begin{aligned} \frac{d y}{ dt} =f(t) \end{aligned}$$

(12)

from quadrature formulae by integrating Eq. (12) to

$$\begin{aligned} y(t+\tau )=y(t) + \int _t^{t+\tau } f(t) d \tau , \end{aligned}$$

(13)

and insert our favorite quadrature scheme, which is the basis of the trapeze rule for time integration. When we replace the trapeze rule with the midpoint rule, we would expect better accuracy, because of the smaller local error Eq. (10) for the midpoint rule, compared to Eq. (11). The actual result can be seen in Fig. 9: While it starts well enough (at least more accurately than the Euler-Integrator would do), the numerical solution starts to develop high-frequency ghost oscillations with exponentially increasing amplitude, even for the (exponentially) decaying exact solution. The culprit is the midpoint value $f((x_i+x_{i+1} )/2)$ in the integration, which, when we look at the original quadrature in Fig. 8, clearly is not an initial or final value for the integral in Eq. 13 - quite in contrast to the values at the end of the interval which are used for the Trapeze rule. The random osciallations around the true values at the initial value of each timestep are then amplified to non-decaying oscillations of high frequency. The same can be expected to happen for the random assignments of zero-order boundary conditions in lattice methods. In fact, such oscillations have already been observed for interpolated boundaries of circles moving over square grids [12], and in truely zero-order fashion the noise amplitude did not decay even for increasing area size of the circles: A good reason to keep one’s fingers away from DNS simulations without the exact implementation of the boundaries. Such an exact implementation of spherical droplets in a fluid, on the other hand, would necessitate the use of curvilinear meshes near the droplet boundaries, and, as the positions change in every timestep, an adaptive remeshing in every timestep becomes necessary. From the standpoint of the first author of this article, who has experience with such issues for polygonal particles and two dimensions [10, 11], the issue of droplets in fluids in three dimensions are at present not yet tractable in direct numerical simulations.