Boundary elements method for microfluidic two-phase flows in shallow channels.

Boundary elements method for microfluidic two-phase flows in shallow channels.


In the following work we apply the boundary element method to two-phase flows in shallow microchannels, where one phase is dispersed and does not wet the channel walls. These kinds of flows are often encountered in microfluidic Lab-on-a-Chip devices and characterized by low Reynolds and low capillary numbers.

Assuming that these channels are homogeneous in height and have a large aspect ratio, we use depth-averaged equations to describe these two-phase flows using the Brinkman equation, which constitutes a refinement of Darcy’s law. These partial differential equations are discretized and solved numerically using the boundary element method, where a stabilization scheme is applied to the surface tension terms, allowing for a less restrictive time step at low capillary numbers. The convergence of the numerical algorithm is checked against a static analytical solution and on a dynamic test case. Finally the algorithm is applied to the non-linear development of the Saffman-Taylor instability and compared to experimental studies of droplet deformation in expanding flows.


I Introduction

Microhydrodynamics is a branch of fluid dynamics that deals with slow viscous flows at small length scales. In recent years the research field of microfluidics investigated the possibilities that microhydrodynamics offers to perform chemistry or biology on a micrometric scale. Such efforts have led to an increasing number of Lab-On-A-Chip applications in the last ten years Whitesides (). In this course droplet microfluidics has emerged seemann12 (), because it exploits the laminar flow in microchannels to precisely control and steer operation on droplets, which act as highly parallelizable reaction chambers.

Microfluidic length scales are in the order of tens to hundreds of micrometers. When microfluidic channels are filled with two immiscible liquids, for instance water and oil, the viscosities are in the order of and surface tension or interfacial tension in the order of , depending on the fluid mixture and surfactants. Due to the small length scale the flow resistance in these channels is high, which is one reason why flow rates usually range between a few to hundreds of with flow rates in the order of . The Reynolds number, Re, is small and therefore, it is often a reasonable approximation to discard the non-linear inertial terms and to consider Stokes flow, which is described in section II.

However, when considering two-phase flow even in the Stokes regime, the dynamics become non-linear due to the free interface between both liquids. The non-linearity stems from domains of different viscosity separated by a mobile interface under surface tension.

Two competing effects dominate the dynamics; one comes from viscous shear and the other from surface tension. The capillary number expresses the balance between viscosity and surface tension: , which is considered here to be between and .

Throughout the article we consider shallow channels that lie in a common plane. Instead of trying to resolve the full three-dimensional problem, we solve a depth-averaged problem, which is two-dimensional. For shallow channels the velocity profile in the thin direction (z-axis) is assumed to be parabolic, a hypothesis that is also used to derive Darcy’s law in two-dimensions (x-y plane). Darcy’s law states that the flow velocity is given by the pressure gradient divided by viscosity and a permeability coefficient , .

Although there have been propositions to account for tangential surface stresses in Darcy’s law Borhan (), the inability to impose tangential stresses and velocities renders this approach incomplete.

In this work we propose the use of the Brinkman equations instead, which include a correction to the Darcy’s law in form of the depth-averaged in-plane Laplacian, a reminiscence of the  Stokes equation. For droplet flows the Brinkman equation was to our knowledge first proposed by Boos et al. BoosThess () and Bush bush97 (), who treated the flow induced by a thermo-capillary effect.

The Brinkman equation is solved with a boundary element method (BEM), which eliminates one more dimension turning the problem from a  differential equation into an integral equation on a  line. While BEM approaches have been followed for  Stokes flows Bazhlekov (), for  Stokes flow wrobel () and Darcy flow Nadim (), recall that simulations of  Stokes flow cannot account for the confinement in the z-direction whereas Darcy’s law becomes invalid close to boundaries and interfaces. The use of the Brinkman equation requires high aspect ratios to justify depth-averaging but we will see that it might still accurately captures the dynamics for aspect ratios approaching . Close to boundaries or interfaces the Brinkman equation gives much better results than Darcy’s law, because it captures depth averaged boundary layers even if the averaged equations become inconsistent BoosThess ().

The derivation of the depth-averaged problem is presented in section II and the numerical method is described in section III together with a stabilization scheme for the surface tension on the interface and an acceleration using Gauss block pre-condensation and multi-core parallelism.

The method is applied to the non-linear development of the Saffman-Taylor instability of finger formation and to the numerical modeling of two recent experimental studies of droplet deformation in section IV. Section V concludes with a brief discussion of the method and its results.

Ii Governing equations

Throughout the article vectors and tensors are written in bold face unless they are represented by a Greek character. Scalars or components of vectors and matrices are written in normal face. All field variables are non-dimensionalized, using a characteristic length scale , the pressure scale and the velocity scale , which are build using the continuous fluids viscosity and surface tension .

Low Reynolds number flows are described by the  Stokes and continuity equation, where non-dimensional operators and variables in are denoted with a tilde.


The non-dimensional parameter compares the viscosity of the considered fluid phase against the viscosity of the carrier fluid, . For the dispersed phase , , and for the continuous phase , . Because of the small size and the horizontal alignment gravitational effects are neglected.

ii.1 Brinkman model for depth-averaged flow

The non-dimensional height is and is considered to be small, . As the flow is confined between two plates at a distance , one considers only fluid motion in the plane and neglects the vertical velocity component, which is equivalent to the assumption of constant pressure in the -direction.

Under this assumption the flow field writes . The two-dimensional velocity vector represents mean velocities, which demands . With these assumptions eq.(1) can be written in terms of two-dimensional variables and operators:


When the profile becomes a parabolic Poiseuille profile and its second derivative in eq.(2) is known. Using the parabolic profile in eq.(2) and depth-averaging over we get the amalgam equation of Darcy equation and  Stokes equation given in equation (3), which is called Brinkman equation and was first applied in granular media flows brinkman (),

Figure 1: Depth-averaged velocity profiles across a parallel channel of width for . The full line represents the solution of the  Stokes equation and the dashed line the solution of the Brinkman equation. The aspect ratio of the profiles shown from bottom up, .

We shall briefly illustrate the advantage of the Brinkman equation through a comparison of its solution for a flow in a rectangular duct of width with the  Stokes solution. The exact solution can be found by separation of variables and is given for instance in Langlois and Deville deville14 (). Depth-averaging the solution of the  Stokes equation gives the mean velocity across the channel.


In comparison, the mean velocity using of the depth-averaged Brinkman equation is:


Both solutions show at leading order a hyperbolic cosine with similar prefactors, for the Stokes equation instead of for the Brinkman equation and in the hyperbolic cosine a factor instead of . The depth-averaged velocity profiles are plotted in figure 1 for different aspect ratios. One observes that the solution from the Brinkman equation tends the solution of the  Stokes equation as the aspect ratio increases. Already for square channels, , the solutions are not too far from each other. Whereas a comparison with the Darcy equation, which is constant in , gives . Far away from the walls the Darcy equation gives correct results for high aspect ratios but fails near walls and for a moderate confinement.

In a more detailed analysis Gallaire et al. Gallaire () showed that even in the complex thermo-capillary flow around a droplet the averaged model agrees almost perfectly with  Stokes. Including the in-plane Laplacian yields two important improvements in comparison to Darcy’s law: 1) Tangential velocities and stress can be imposed on boundaries and 2) there appears a boundary layer near walls and interfaces that scales like , the non-dimensional height of the channel.

ii.2 In-flow and out-flow boundary conditions

Boundary conditions of the single-phase problem prescribe either the stress or the velocity. The typical no-slip boundary condition on channel walls is . In contrast to Darcy flow, the Brinkman model imposes normal and tangential velocities. The normal and tangent are given by a vector that contains their projections on the and axis, e.g. .

As typical inflow boundary condition the solution of the Brinkman equation in a straight channel flow is used. For a straight inflow boundary of length parameterized by , whose origin is in the middle of the boundary:


It is worth observing that the dimensionless inflow velocity is represented by the capillary number because the velocity is non-dimensionalized by . Integration of the velocity field along the inflow boundary relates the dimensional flow rate to the inflow capillary number.


The outflow boundary of parallel channel flow with constant pressure is imposed by zero tangential velocity, and constant normal stress .

ii.3 Droplet interface condition

Droplet interface conditions are continuity of velocity and discontinuity of interface stresses. The surface or interface stresses shall also be denoted by , whose components are surface stresses . Over the interface the normal stress is discontinuous due to surface tension and curvature. The tangential stress can also be discontinuous due to a varying surface tension. The varying surface tension is introduced by . The two principal curvatures are the in-plane curvature and the meniscus curvature in the thin direction .


In accordance with real droplet microfluidic systems we assume a non-wetting dispersed phase, thus the out-of-plane meniscus is approximately a half-circle, , and the in-plane curvature is corrected by a term, which was derived by Park and Homsy park84 (). With a perspective on discretization the stress jump is also expressed in terms of Frenet equations, which is proposed by Tryggvason and co-authors zaleski ().

In this formulation of the interface stresses jump we neglect the effect of dynamic film formation, which changes the out-of-plane curvature . This effect, which depends in a non-linear way on the capillary number, see Park et al. will be discussed in the conclusion.

Iii Numerical method

In order to solve the two-phase flow numerically the domain has to be discretized. It is desirable to apply an interface-tracking scheme, as the interface represents a localized force. This force is relatively high because it competes with viscous forces that scale like the capillary number, which is small.

Using a diffusive interface instead of a discrete interface can lead to problems at low capillary number. For instance Carlson et al. amberg10 () used a phase field method with a diffuse interface to study the dynamics of droplets in a bifurcation in at low Reynolds numbers. They observed a dependence of the solution on the thickness of the interface and had to use a very small time step. In order to avoid these difficulties we discretize the interface with mesh elements for the sake of precision and stability of the numerical algorithm.

A Boundary Element Method is implemented, where only boundaries are discretized. As a consequence it is unnecessary to remesh the whole domain as the interface evolves. For a complete description of the method and applications in viscous flow one may consult the book by Pozrikidis pozrikidis (). In this book the Brinkman equation is mentioned as a mean to compute unsteady flows, where a discretized time derivative due to the local acceleration is represented by . In the following section we elaborate the two-phase flow boundary integral form of the Brinkman equation and its boundary conditions. The procedure closely follows that for  Stokes flow wrobel () but we include it for a complete and comprehensive presentation of the method.

iii.1 Boundary integral formulation

Transformation of the Brinkman equation into a boundary integral form is done by integration by parts, which requires expressing the Brinkman equation from eq.(3) using the divergence of the in-plane stress tensor.

The in-plane stress tensor is,


Hence the Brinkman equation can be expressed as:


As a first step we integrate the Brinkman eq.(10) and continuity equation over a vector field and a scalar field on a domain one obtains:


Here and can be seen as test functions, as in the Finite Elements Method. The continuity equation has been multiplied by , which will be necessary later on for the droplet interface boundary condition. If the velocities and pressure fulfill this equation for whatever choice of test function they also solve the initial equation system eq.(3).

Performing integration by parts on eq.(11) gives,


The term that appears as a boundary integral on the boundary of the domain appears as a surface stresses or . Performing integration by parts once more makes and exchange roles.


After integrating by parts twice the Laplace operator that acted on the velocity has been transferred to . Since the vector field takes the role of a velocity in a reappearing Brinkman equation, can be interpreted as a pressure and the new boundary integral term can be interpreted as a surface stress :


When , and solve the Brinkman equation the domain integral disappears and only boundary integrals are left. Green’s functions of the Brinkman equation is used a test function because they have a Dirac like forcing, which ensure that in a numerical discretization of the problem the highest elements are usually the diagonal elements of the matrix, leading to a low condition number of the problem. A Dirac forcing on the boundary also allows expressing the domain integral solely based on the velocities at the location of the Dirac distribution.

iii.2 Green’s functions

In order to solve for the unknown velocities and stresses two Green’s functions are necessary. A third Green’s function can be used to compute the pressure in the domain. As a consequence of the Dirac forcing the velocity, pressure and stress fields have singularities at , where their respective values go to infinity.

In the following we present at first the two Green’s function of the Brinkman equation with a Dirac forcing of the stress equation in or direction. Taking two velocity vectors , two pressures and two stress tensor fields , with , and which verify:


The Dirac forcing occurs at the position with a vectorial weight . We are interested in two linear independent weight vectors and .

The solutions are given in indexed form, first index for the corresponding weight function and for the component of a vector or the component in a matrix. The Green’s function velocity field, pressure and stress tensor are taken from Pozrikidis pozrikidis ().


Introducing and , which are functions of aspect ratio times distance and are expressed by modified Bessel functions of second kind and . In the above and following definitions is the Kronecker delta and not to be confused the Dirac distribution .


The function behaves like a close to and is therefore weakly singular. Using this result in eq.(15) leads to the associated pressure field .


From velocity and pressure fields one derives the associated stress tensor, :


Inserting now and with or into eq.(13) gives two independent boundary integral equations with a forcing either in the or direction. Because of the Dirac forcing the domain integral does not disappear completely but is confined to a single point and eq.(13) becomes:


The domain integral over the Dirac selects the value of at , the position where the Green’s functions was forced. Depending on whether is in the domain, on the boundary or outside the prefactor is or . When it is desired to determine the pressure in the domain, one uses a third fundamental solution with a Dirac distribution forcing of the continuity equation, which corresponds to a point source in the continuity equation at . Instead of eq.(15) the Green’s function solves:


This point source solution is also solution of the Darcy equation and it’s fields are:


Inserted in eq.(13) this gives:


Again the prefactor depends whether is in the domain, on the boundary or outside, . Evaluating this integral equation allows to calculate the pressure at a given point in the domain once the stresses and velocities at the boundary are known.

iii.3 Boundary integral for two-phase flow

When dealing with two-phase flows the two fluid domains need to be coupled by their common boundaries. The channels internal boundary is the same as the droplets boundary, except for the sense of orientation. The orientation determines the direction of the normal and therefore changing the sign of one of the normals makes both integral paths identical. One includes the droplet into the boundary integral formulation by summing the contributions of the boundary integral of the droplet and the channel domain, where the channel domains internal boundary has a reversed integration path, i.e. inverted normals. Boundary integration is performed counter-clockwise around a domain, illustrated in Figure 2.

Figure 2: Schematic representation of how the two fluid domains are coupled. Arrows depict the boundary integral paths. a) Two separate domains. b) Changing the orientation of the channels internal boundary, . c) Combination of both integral contributions on the interface using the interface conditions in eq.(8). Piece-wise constant and piece-wise linear shape functions are sketched.

In the channel domain the viscosity ratio is and in the droplet . The dispersed phase boundary integral problem is:


Combining eq.(24) and (25) by adding the contributions under the integral along gives,


If the singularity at is located on the outer boundary the coefficient and if the singularity is located on the interface . Prescribing the jump in viscosity and jump in surface stresses , see eq.(8), specifies the fluid-fluid interface conditions and poses the combined problem.

A significant simplification is obtained for , where from eq.(26) only a diagonal matrix and right hand side remains,


In the absence of walls the velocity at the point depends only on known quantities and does not require the solution of a matrix. Even in the presence of outer boundaries the evolution problem can be reduced by one matrix inversion in a problem that involves only matrix vector multiplications, which will be shown in section III.6. We do not restrict our method to it is worth mentioning that Zhu et al. Power () have used this in the case of droplet motion.

iii.4 Spatial discretization

The boundaries are discretized by piece-wise straight line elements with local shape functions whose amplitudes are to be determined. Fixed boundary conditions like walls and inflow and outflow boundaries are discretized with piece-wise constant shape functions, whereas the droplet is discretized by piece-wise linear shape functions. The boundary integral form of the Brinkman equation, eq.(26) is discretized with a collocation method and integrated by Gauss-Legendre quadrature with nodes.

The quadrature uses weights , which are evaluated at base points . For instance applied to a generic real valued function integrated between and


Discretized normal and tangential vector , element length and local variable are given as:


The integral equation, eq.(26), turns into a double summation, the sum of the quadrature points, summed over all fixed elements and droplet elements.


In eq.(30), there are unknowns to be solved for; velocities or surface stresses on the channel boundaries and velocities on the drop interface. For each position of the test functions on one of the nodes, there are two discretized equations, eq.(30), for or . Summing up all known variables in the right hand side and forming a dense matrix with a vector of unknowns, yields a linear equation system.

Due to the weak and strong singularities of the Green’s functions, these functions diverge on the collocating element although the integral does not. In numerically integration the singularities need special attention. The test velocity is weakly singular like a logarithm; in fact its asymptotic development in terms of Bessel functions contains a . If the integration interval contains a singularity, the contribution is not evaluated in Gauss quadrature but integrated analytically.

The last member of eq.(19) for the stress contains a singularity but its divergent behavior is perpendicular to the normal and vanishes. To show this we consider an element collocating with a singularity, whose normals can be written as: .

Hence the strongly singular contribution cancels out, which allows setting them to zero on a collocating element and leaves an expression that can be numerically integrated. The numerical approximations of the Bessel functions with a subtractable singular term are given in the appendix A.

Velocities at the droplet interface are unknown and the provided interface conditions are the jumps in surface stresses, which we discretize from eq.(8).


where is defined on vertices and midpoints.

Validation: Marangoni flow of Boos and Thess

We demonstrate the convergence of the method by comparison to an analytical solution. A droplet placed in a Hele-Shaw cell and subject to a surface tension gradient was studied by Boos and Thess BoosThess (). The gradient in surface tension will lead to a motion that is tangential to the interface, which is named Marangoni effect. In their theoretical study they assumed a cylindrical droplet that is exposed to a linearly changing surface tension gradient. Their result is of particular interest because it provides a test case with a shear stress boundary condition, something impossible when using the Darcy equation. The length scale is equal to the droplet radius, hence we apply boundary conditions for a droplet of radius and surface tension gradient in the y-direction are radial velocity and tangential stress discontinuity .

The  Brinkman solution by Boos and Thess is given as:


The amplitude of the tangential velocity along the interface decreases with that is proportional to the channel height . A high inner viscosity , further reduces the tangential velocity.

Figure 3: Convergence for a cylindrical droplet exposed to a surface tension gradient. The error is given as and plotted against the number of elements on the circumference. L/H = 3, = 1/3, L/H = 6, , L/H= 3, , L/H=6, . The inset visualizes the flow field for and .

Comparing theoretical velocity against numerically obtained velocity shows a decrease of the maximum error with second order, shown in figure 3. Two aspect ratios and and two viscosity ratios and were used.

iii.5 Temporal discretization and stabilization of the fluid-fluid interface

Two schemes are implemented to integrate the solution in time, a one step explicit Euler scheme and a two-step Runge-Kutta scheme (Heun’s method).

Using the explicit Euler scheme the droplet is advanced in discrete time steps:


Where is the velocity field obtained by solving the boundary element problem eq.(30) using the nodes at .

Alternatively we use a two-step Runge-Kutta scheme:


With being the intermediate velocity field obtained using the nodes at .

During evolution, when the element size on the interface is twice as big or twice as small as the initial size, the points are redistributed equidistantly with the help of a rd order polynomial fit.

The time step is limited as a consequence of the mobile interface and non-linearities associated with surface tension and coupling domains of possibly different viscosities. In the limit of low capillary numbers when the surface tension dominates over viscous dissipation the problem evolves very slowly whereas the time step for stable integration is limited by the element size. For instance the circular interface with constant curvature of a droplet at rest, which is physically stable, develops nevertheless divergent oscillations when solving the evolution problem for a time step larger than some factor multiplied by the discrete element length . For a viscosity ratio , aspect ratio and a droplet of radius the maximal time step is about . This stiff constraint has been observed by others, e.g. Dai et al. shelley93 () when analyzing finger formation in a Hele-Shaw cell at low capillary numbers.

Figure 4: Relaxation of a droplet with eccentricity , area , aspect ratio and viscosity ratio . The black (resp. gray) line represents the stabilized scheme with elements and (resp. ). The non-stabilized scheme is represented by the green line , blue , red and green dashed line .

The numerical instability is illustrated here on a dynamic problem, where a droplet whose initial shape is an ellipse of eccentricity and area relaxes to a circular interface. A discretization of 100, 200 and 400 elements is tried with time steps on the threshold of instability of = 8 , so and . In figure 4 the difference between maximal and minimal in-plane curvature is plotted, which starts near two and goes to zero for the perfectly cylindrical droplet. In plotting the curvature, oscillations will appear much more pronounced. The time steps on the threshold lead to oscillations due to deformations on the evolving interface. As the time step is divided by two the scheme becomes stable, as shown by the green dashed line with and .

Hou et al. shelley94 () developed a linearization technique to remove the stiffness from surface tension but their approach requires a transformation of the equations in terms of tangential and angular coordinates (s,). We propose a different approach where the curvature is also linearized, however the variables remain in cartesian coordinates.

Interface stabilization technique

A semi-implicit scheme is used to stabilize these numerical capillary waves and therefore allow for larger time-steps. The scheme used here in two dimensions extends in a straight forward manner to three dimensions.

For a discrete interface the in-plane interface stress is derived from eq.(31) at the vertex at time step as:


Stabilization is achieved when the interface is ”made aware” of its displacement when providing a feedback loop. The position of the points a discrete time step after time is approximated by eq.(34) to be . Presuming the change in distance between two points to be negligible (), the linearized expression for the interface stress after a time increment is then given by:


A similar stabilization is applied to the second term in eq.(31). Using the stabilization scheme for the relaxing droplet allows for much larger time step. Computing the droplet relaxation with 400 nodes and a time step of and agrees with the non-stabilized solution that was computed at a time step that is 16 and 160 times smaller, respectively. For the stabilized scheme no oscillations are observed and stabilization is achieved independent of the discretization.

Convergence study: Deformable droplet in flow focusing

In lack of an analytical solution for the convergence study of deformable droplets we resort to comparison with a numerical result of increased spatial and temporal resolution. For this study a cross flow junction is simulated, where a fluid stream is focused by flow from two side channels. A droplet of diameter submitted to that flow is deformed and accelerated, figure 5 shows the numerical set-up with time lapsed droplet position. The left inflow condition is chosen to be and the inflow velocity from the side channels is .

Droplet viscosity is half the viscosity of the surrounding fluid, and droplet aspect ratio . Variables that are observed are the displacement of the center of mass of the droplet and its perimeter.

Displacement error is measured by the root-mean-square (rms) of the difference in position at time steps normalized by the total distance. The perimeter error is the rms of the difference of the droplet perimeter, sampled at time steps and normalized by the perimeter for an undeformed droplet .

Spatial discretization varies between as element size for the walls and for the droplet to account for the importance of the mobile interface. Temporal discretization ranges between . The reference solution uses and and , respectively.

Figure 5: The figure shows the computational domain with the droplet interface at several time steps. The three inset figures show the in-plane curvature at three selected time steps around the most deformed. The interface is free of any spurious oscillations.
Figure 6: Left figure a) convergence for varying element size , right figure b) convergence for varying time step . The blue line represents the displacement error and the red line the perimeter error . Full line for 1-step integration and dashed line for 2nd order Runge-Kutta scheme.

Figure 6 a) shows the errors for a fixed temporal resolution of and confirms second order convergence as expected from the previous study in section III.4. In fact the error decreases almost with 3rd order, which maybe due to the symmetries in the configuration. Figure 6 b) shows for a fixed spatial resolution of the temporal convergence. The convergence of the stabilized Euler scheme is of order as expected. However the second order Runge-Kutta scheme (Heun’s method) shows also first order convergence although with smaller error, which is likely due to the fact that an intermediate time step was performed. The semi-implicit interface stabilization scheme incorporates a first order Euler scheme and thus spoils any higher order scheme.

Without interface stabilization the spatial convergence, done with would fail for on the droplet, which follows from the empirical formula above . Likewise the temporal convergence with fixed spatial discretization of would require a time step as low as that is eight times lower than the time step in the refined solution.

Due to incompressible flow the area of the droplet, initially is theoretically conserved. In the worst case, the scheme with the lowest resolution finished the simulation loosing of the area, whereas the scheme with the highest resolution lost .

iii.6 Block pre-elimination and parallel scaling

The discretized problem results in a dense matrix with right-hand-side and vector of unknowns , containing velocities and stresses. The degrees-of-freedom (DOF) of the problem are . DOF associated to the static interface, the outer walls, and DOF associated to the dynamic interface, the droplet. If the droplet is surrounded by walls one generally finds that the DOF, , of the walls are larger than the DOF, , of the droplet. The problem without droplet is without evolution and independent of time. One therefore splits the matrix into four blocks:


Here is the influence of the static walls on themselves, like a resistance, is the influence of the droplet on the outer walls, like propulsion, is the influence of the outer walls on the droplet and the influence of the droplet on itself.

Matrix sizes are is , and are and is . The matrix is always the same because it’s boundary conditions and its element distribution does not change. We invert the matrix and save the inverse .

Applying the inverse to the upper part of the eq. 38:

Then replacing in the lower part of eq. 38:

And finally:


With called the Schur complement. These matrix-matrix and matrix-vector multiplications reduce the problem from a dense matrix to a dense matrix. In the case of equal viscosities , and is the identity matrix, so the velocity at the droplet interface is given explicitly by: .

The described boundary element code has been implemented in the C++ programming language. Independently of the block pre-elimination we have used shared memory parallelization. A direct matrix solver for general matrices implemented in LAPACK lapack () solves the matrix. Dense matrices in contrast to sparse matrices have stronger limitations for parallelism on distributed computers. A rather efficient way is to use OpenMPopenmp (), which allows shared memory parallelism on a multiprocessor and multicore environment. LAPACK routines automatically use these features, whereas the loop for matrix filling has been parallelized by OpenMP pragma.

At each iteration in time, matrix filling is done by nested loops. An outer loop over all boundary elements and a nested loop that integrates over all Greens functions at the collocation points, where the numerical integration is a further nested loop. This inner most loop over only six points is too small for efficient parallelization. Hence the loop over all Greens functions is parallelized.

Scaling for parallelism and pre-condensation The speed-up is demonstrated solving for a droplet advected in a rectangular channel using pre-condensation and multiple cores. The boundaries were discretized with 1440 DOF and its multiple by five and ten (7200 and 14400 DOF). The ratio between DOF on the droplet and on the fixed geometry, was changed between and .

Figure 7: Speed-up of the simulation by: a) OMP implementation, where the speed-up increases almost linearly to the number of cores, mean of and DOF shown. b) Pre-condensation, where the speed-up increases approximately linearly with the ratio of droplet to wall degrees-of-freedom.

The performance tests were performed on a Dell Server with 16 cores at 1.8 GHz. The computation time for a single iteration of 1440 DOF on a single core was with 7200 DOF and with 14400 DOF . Practically all the time is spent on matrix filling and matrix solving. There is an almost perfect scaling of computation time for matrix filling and for matrix solving.

In the presence of outer walls the fixed boundary problem can be pre-condensed and therefore reduces the time to fill and solve the matrix. In assessing the speed-up we ignore the time associated to pre-condensation, which includes one matrix inversion. The matrix inversion is more costly than direct solving of a matrix but it is done only once and becomes negligible to the overall time. The larger the fixed part of the matrix the higher the speed-up.

Increasing the number of cores shows a good scaling up to 16 cores for the problems with and unknowns, which achieve a speed-up of about compared to the single core case. The parallelism work about equally well for matrix filling and solving. The problem with unknowns shows more waiting time and achieves only a speed-up of 7. We shall exclude the results so the trends are uniform but apply only to the large scale problems, which are in fact the ones that have the highest need for acceleration. Figure 7 a) shows the speed-up time on a single core divided by time on multiple cores.

Comparing problems of varying ratio in figure 7 b) shows the mean speed-up out of four configurations, from DOF on single core and 16 cores and DOF on a single core and 16 cores. The error bars are computed from configurations different DOF with and without parallelism and indicate that the speed-up is quite uniform also when using OpenMP. The problem with DOF was stagnant at a speed-up by factor 3.

Iv Applications compared to theory and experiments

The physical significance of the presented method is demonstrated on three cases: Comparison to a linear stability analysis of the Saffman-Taylor instability and two experiments on droplet deformation.

iv.1 The Saffman-Taylor instability

Finger formation as a consequence of the Saffman-Taylor instability is an archetypical problem of dynamically developing interfaces. It has been studied numerically using Darcy’s law for instance by Dai et al. shelley93 () who looked at the amplification of small perturbations on the unstable interface.

The instability and resulting finger formation is due to a viscosity change over an interface, where a less viscous fluid drives a more viscous fluid. In opposition, surface tension damps surface perturbations that would cause highly curved regions. For a radial configuration the growth rate using the Brinkman equation was derived recently nagela () and showed better agreement than Darcy’s law for relatively large capillary numbers.

Introducing a sinusoidal perturbation of the interface, the growth of the perturbation is given as . The initial aspect ratio, radius over height, was set to , the viscosity ratio to and the capillary number was varied between and as to obtain unstable, neutral and stable modes. The instability is numerically simulated with a droplet interface of elements and a time step of .

Plotting perturbation amplitude against time for analytical and numerical results in figure 8 shows perfect agreement in the beginning of the simulation, then the simulation is slightly parting from the prediction as the radius evolves. The background illustration in the figure shows the interface evolution for every time step. Deviation for large at is due to non-linear saturation and at because the aspect ratio increases with the radius and therefore the solution shifts into the unstable regime.

Figure 8: Plotting the growth of a sinusoidal perturbation of wave number while injection of a liquid that is times less viscous than the surrounding liquid. The evolution of the interface for is depicted at several instances in time. The graph in the center compares prediction by linear stability analysis - - - against numerical simulation —.

When starting from stochastic initial conditions these simulations have allowed to observe tip splitting and opens perspectives to study the selection of number of fingers. When adding Marangoni stresses it is also possible to analyze the evolution of fingering in a medium of variable surface tension, like in a Hele-Shaw with thermal gradients.

iv.2 Comparison with experiments: Droplet stretching

Two different cases of droplet stretching in microchannels are presented; first at large and then at low aspect ratio.

Stretching in a hyperbolic flow Recently Ulloa et al. cordero14 () experimentally studied droplet deformation in diverging flow at low Reynolds numbers and aspect ratio of droplet radius to channel height . The viscosity ratio of inner fluid to outer fluid is in all their experiments as well as in our simulations. Their set-up is a microfluidic cross-junction with fluid injected from two opposite sides, where the fluid leaves the junction in two channels, oriented to the inlet channels.

From left inlet channel droplets enter the junction and get advected to the center of the junction, a stagnation point, and stretched by the flow. This stagnation point is a saddle point and due to imperfect alignment and the unstable nature of the equilibrium position the droplet gets streamed away after some time. Before getting streamed away the droplets have deformed in the flow field to an equilibrium shape. Ulloa et al. cordero14 () characterized the deformation by the difference in major and minor axis length and , normalized by their sum, . They investigated the influence of channel geometry, shear rate and droplet radius on the deformation.

Figure 9: Multiple snapshots of a droplet that enters a junction from the left. Gray arrows indicate the incoming flow. Capillary number , channel width and droplet radius . The inset b) shows droplet deformation, plotted against the distance between droplet center and the center of the junction.
Figure 10: Maximum deformation against a) droplet aspect ratio at capillary number and channel aspect ratio and b) against capillary number for droplet aspect ratio and channel aspect ratio . The black line is the trend extracted from the experimental data and the circles correspond to simulations.

Choosing that channel half width as a length scale , a junction of width , extending from to in the -direction and to in the -direction is simulated. For illustrative purposes figure 9 shows a simulation with one droplet at several instances. One sees the droplet coming from the left, stretching in the center and leaving the junction on the top. A maximum deformation of is reached at the center as shown in the inset.

In figure 10 we compare the scaling of maximum deformation with droplet radius normalized by the channel height in a channel of width , at constant capillary number . Furthermore we compare to the shear rate that has been non-dimensionalized by the viscosity, surface tension and channel width in a Channel of . Numerical simulation shows a dependence in and in . Experimentally Ulloa et al. found .

Stretching in extensional flow In a recent work Brosseau et al. barets () investigated the influence and adsorption of surfactants. They follow the idea of a microfluidic tensiometer, which was brought up by Cabral and Hudson cabral (), who proposed to measure the surface tension by observation of droplet deformation in an extensional flow created in a micro-channel.

Influence of surfactants is difficult to determine experimentally and microfluidics might help to investigate characteristics that are inaccessible by established means. When surfactants are absorbed on the fluid interface they change the surface tension and introduce surface tension gradients. Modeling surfactants in numerical simulations is a formidable task and the flow solver we developed is able to simulate deformable interfaces and surface tension gradients, while being computationally inexpensive as to run many simulations in order to fit surfactant models to experiments by reverse engineering. In the adsorption limited regime no convection-diffusion equation needs to be solved as the surfactants are spread homogeneously over the domain. In contrary in the diffusion limited regime an advection diffusion equation needs to be solved in addition to the Brinkman equation.

Figure 11: a) Droplet stretching at capillary number and b) at ; Droplet contours at several moments in time are plotted. Graph c) shows the evolution of the deformation at different positions of the droplet. Grey dots show the experimental data from Brosseau et al. for and . In figure d) the maximum deformation is plotted for different aspect ratios and . Experimental results by Brosseau et al. for are drawn by a dashed line.

Here we compare the first series of experiments of Brosseau et al. barets () without surfactants to our simulations. The channel they used has a constant height of and has wide section, wide and long. This wide section has an entry and an exit channel wide. For non-dimensionalization the length scale is the entry channel width , viscosity ratio is and channel aspect ratio , as in the experiments. Simulation of square channels are out of the nominal range of validity of the presented Brinkman model. Nevertheless the value of this comparison is to estimate the limitations of our method when approaching low aspect ratios.

The wide channel section is units long and units wide and joined by unit wide entry and exit channels. The edges where these channels join have been rounded with a radius of that can also be observed in the experimental setup. Three different droplet aspect ratios were used and .

The experimental study used the droplet velocity in the entry channel to set the capillary number, which we shall denote here . Illustration of the extensional flow together with the resulting deformation for and and radius is shown in figure 11 a) and b). At capillary number the droplet is more stretched in the entry channel, till and then also expands much wider than the droplet at . In figure 11 c) one can see this behavior in the deformation curve, where the initial deformation is negative, which corresponds to elongation. The deformation is measured by the difference in major and minor half axes and , . At the flow expands and shortly after the droplets reach their maximum deformation. Results are shown for the examples a) and b) as well as for a comparison with a given experimental data set, showing a good agreement.

The maximum deformation is shown in figure 11 d) against the capillary number for radii and . The experimentally obtained scaling law between deformation and capillary number is only approximately retrieved for and , with a root-mean-square error of about 40%.

It shall be stressed that although the channel aspect ratio is not in the regime of microchannels with large aspect ratios qualitative agreement could be achieved. Furthermore the Reynolds number was beyond the domain of validity of the Brinkman equation. The Reynolds number being estimated based on the properties of the carrier fluid (fluorinated oil: ), channel height and maximum velocity of the droplet , yields a maximum value of . Despite the relatively large values no dependence on the Reynolds number could be established from the experimental data.

A simulation of droplets with lateral confinement, as shown here, would not have been possible with the Darcy equation because without the in-plane Laplacian no coating boundary layers build up at the lateral boundaries of the droplet, resulting in collision of the droplet interface with the lateral walls.

V Conclusion

The boundary element method that is presented here solves for two-phase flows in geometries of moderate aspect ratio. Convergence of the numerical scheme has been verified as well as functionality of an interface stabilization scheme and acceleration due to Gauss block pre-condensation.

The dynamic evolution of the interface in time has been verified with the dispersion relation for the Saffman-Taylor instability based on the Brinkman equation. The perturbation amplitude of the linearized equations has been compared to the numerically obtained amplitude and results agree very well showing that the solver preserves stable, neutral and unstable modes in a non-linear problem

The deformation of a droplet in a linear flow in a wide channel at low Reynolds number agreed well with experimental findings, where the droplet aspect ratio and capillary number were varied. In a complementary study the validity of the  model equation was checked in the limit when channels have a square cross section, where droplets deform in a sudden expansion. Main results as the deformation over position diagram closely resembles the data obtained in experiments, as well as the scaling of the maximum deformation over capillary number are recovered. However the  model breaks down at low capillary numbers, where the surface tension forces the droplet interface very close to the walls.

Furthermore only flows in the absence of inertia can be solved by the boundary element method without resorting to domain integration, imposing a Reynolds number , whereas in some cases microfluidic experiments have Reynolds number larger than . Domain integration or dual reciprocity methods pozrikidis () can be used to include these terms, but they introduce a considerable effort compared to the presented procedure. In the experimental study at finite Reynolds number barets () to which we compared our simulations, the results seemed to be independent of the Reynolds number. An explanation could be that the Weber number, which measures the inertial terms compared to the surface tension remains small, at most in the experiments of Brosseau et al.. Hence the surface tension might be expected to dominate over inertial terms and the dynamics of the interface become not too much affected.

Given the restrictions of being a  approach with domain wise constant parameters it should be stressed that the method is rather fast. Simulations run on a desktop computer typically finish in a couple of minutes or hours depending on the problem size, which in return allows doing extensive parametric studies and simulating a large number of droplets.

The acceleration is due to the reduced number of degrees of freedom when using BEM but also to acceleration techniques and the interface stabilization. The latter is important since the droplet evolves on a capillary time scale that is , which for typical values like and gives . Considering for instance an element size of elements per unit length and a admissible time step it will take about iterations to advance in physical time. Applying stabilization allows for with no restriction on the spatial discretization. Without any kind of interface smoothing scheme the interface remains oscillation free. The relatively low computational effort comes to a price: Depth-averaging is based on a high aspect ratio and the soundness of the results becomes questionable when the droplets are not sufficiently confined.

The variety of problems that we solved demonstrates that the tool has the capability to simulate the flow of deformable droplets in complex channel geometries encountered in Lab On A Chip applications. In two recent publications we applied the solver to study droplet relaxation nagel2 () and droplet trapping nagel3 () and have obtained good agreement between simulation and experiments. Current work includes modified boundary conditions that take into account the effect of film formation by coupling the interface boundary condition to an asymptotic solution that modifies the out-of-plane curvature according the asymptotic developments of Park and Homsy park84 ().

The fact that the physical description and numerical algorithm is reduced to its essentials makes it relatively easy to implement and efficient to run. Therefore this work could provide a further step towards numerical simulation of droplet microfluidics and therefore enable experimentalists to use simulations to design or evaluate their Lab-on-a-Chip set-up.


We thank the European Research Council for funding this work through a starting grant (ERC SimCoMiCs – 280117).

Appendix A Numerical approximation of Bessel functions

The modified Bessel functions of second kind are developed as a series. The modified Bessel differential equation is:

Whose solutions are the modified Bessel functions of second kind and , where we are only interested in and . The differential equation has two singular points, one at and another at . Therefore the functions are developed in two zones around those points and meet at . An estimation of the maximum error around amounts to about . Note that in the definition of the Green’s function, eq.(17).

Expansion around

The numerical values of the terms obtained in the asymptotic development in powers of have been slightly modified as to minimize the error in the interval , such that the error is smaller than everywhere, sacrificing the convergence rate for large values of , where the and are almost zero anyway.

Expansion around

With , where is the Euler-Mascheroni constant: