# Generalized Navier Boundary Condition for a Volume Of Fluid approach using a Finite-Volume method.

###### Abstract

In this work, an analytical Volume Of Fluid (VOF) implementation of the Generalized Navier Boundary Condition is presented based on the Brackbill surface tension model. The model is validated by simulations of droplets on a smooth surface in a planar geometry. Looking at the static behavior of the droplets, it is found that there is a good match between the droplet shape resolved in the simulations and the theoretically predicted shape for various values of the Young’s angle. Evaluating the spreading of a droplet on a completely wetting surface, the Voinov-Tanner-Cox law () can be observed. At later times scaling follows , suggesting spreading is limited by inertia. These observations are made without any fitting parameters except the slip length.

## I Introduction

For any multi-phase flow over a wall where dissipation at the contact line becomes of the same order as the bulk dissipation, a good understanding of contact line behavior is essential (Squires and Quake, 2005). However, despite many investigations (Dussan, 1979; Gennes, 1985; Coninck and Blake, 2008; Eral et al., 2013), for applications including adhesion of liquids to a solid surface (Furmidge, 1962), transport of liquid water in fuel cells (Zhu et al., 2008), liquid infused surfaces (Wexler et al., 2015), and coating (Vandre et al., 2013), contact line behavior is still not well understood.

The physics of a static contact angle between a liquid and a gas on a smooth solid surface is well established (Young, 1805; Gauss, 1830). However, real surfaces are not completely smooth. They are not chemically homogeneous, and/or have roughness. This causes static contact angle hysteresis and contact line pinning, both of which are difficult to model. When looking at a moving contact line instead of a static contact angle things get even more complicated. Dynamic contact angle behavior is not well understood, even on a completely smooth solid surface. The origin of this poor understanding of the moving contact line is twofold: there is the contact line singularity problem, and the question of how the contact angle depends on contact line velocity.

The contact line singularity problem was first identified by Huh and Scriven (1971). While normally it is a good approximation to use the no-slip condition as boundary condition on the wall, for corner flow this assumption causes the viscous stress and pressure to scale as and thus to diverge as at the contact line. Numerous methods have been proposed to resolve this discontinuity or work around it. Hocking (1977) showed that, using domain perturbation method in cylindrical coordinates, any slip-velocity model (Navier, 1823; Thompson and Troian, 1997) resolves the velocity singularity. Another method to circumvent the contact line singularity has been to use the Cahn-Hilliard-van der Waals model (Cahn, 1977) as the basis for either diffusive interface models (Seppecher, 1996; Jacqmin, 2000), or for precursor film models (Hervet and de Gennes, 1984; Gennes, 1985). Precursor film models use a disjoining pressure (Derjaguin, 1955; Teletzke et al., 1987; Herring and Henderson, 2010) to model the van der Waals forces that cause the formation of a precursor film ahead of the interface, removing the singularity. More exotic models have suggested local shear-thinning (Cox, 1998) and non-constant surface tension (Shikhmurzaev, 1993, 1997) as possible solutions.

For the question of how the contact angle depends on the contact line velocity there are also various models. Typically a distinction is made between the local microscopic contact angle at the contact line and the macroscopic apparent contact angle, which is observed in experiments. Due to experimental limitations to access sufficiently small length scales, the apparent contact angle is measured away from the contact line, and curvature of the interface causes this angle to be different from the microscopic contact angle (Chen and Wada, 1989). Arguably, the easiest method to define the dynamic microscopic contact angle is to simply assume it is fixed and the same as the static contact angle (Eggers, 2004). Instead of a fixed dynamic contact angle, Molecular Kinetic Theory (MKT) (Blake and Haynes, 1969; Blake, 2006) predicts a microscopic contact angle which changes with the contact line velocity. The Voinov-Tanner-Cox (Voinov, 1976; Tanner, 1979; Cox, 1986) law describes the relation between the apparent contact angle and the microscopic contact angle. This law is based on the assumption that it is possible to choose a length scale arbitrary close to the contact line. This makes it impossible to identify a characteristic length scale of the contact line geometry, and reduces its physics to a balance between capillary and viscous stresses (Kafka and Dussan, 1979). Using the lubrication approximation one can now derive the Voinov (1976) equation for some specific asymptotic limits, and matching solutions of the Voinov equation with the mesoscopic hydrodynamic solution further away from the wall then gives the Cox-Tanner-Voinov law (Voinov, 1976; Tanner, 1979; Cox, 1986).

Although there are a lot of different models to describe contact lines, there is no consensus on what is the correct description of the physics. Many models have multiple fitting parameters which can be tuned to give the same results (Sibley et al., 2014). Because impurities and surface heterogeneities have a large effect on measurements, experiments are very difficult to reproduce. On top of this, one needs access to microscopic length scales to get to the details and the outcome of experiments on a macroscopic level only depends very weakly on these small length scales (Bonn et al., 2009). With the advent of Molecular Dynamics (MD) (Koplik et al., 1988, 1989; Thompson and Robbins, 1989, 1990) simulations, and new experimental techniques, such as Atomic Force Microscopy (AFM) (Checco et al., 2006), it is now possible to access length and time scales at the contact line that were not accessible before. While MD simulations still can only probe small systems for short times, a couple of fundamental discoveries have been made. Although it had been argued that continuum models break down all together at the contact line (Dussan, 1979), it was found that both the Navier-Stokes equations and Young’s equation hold up even down to the nanometer scale (Bocquet and Charlaix, 2010; Seveno et al., 2013). Furthermore, support was found for contact line slip (Thompson and Robbins, 1989), precursor films (He and Hadjiconstantinou, 2003), and non-static dynamic microscopic contact angles (Shen and Ruth, 1998). Apart from simulations, precursor films have also been found experimentally (Hardy, 1919; Kavehpour et al., 2003).

In addition to the experimental difficulties and often contradicting findings, another reason that there is no consensus is that all of the above mentioned models have some fundamental short comings and/or work in different regimes. Further analysis of contact line slip models by Huh and Mason (1977) showed that, even though stresses are not diverging anymore, the pressure still shows a weak singularity i.e. the pressure diverges, but becomes finite when integrated. While using a slip length model is sufficient from a modeling perspective, a divergent pressure cannot be right from a physics perspective. The precursor film on the other hand does successfully work around the contact line singularity for both viscous stresses and pressure. However, they are typically not seen under partial wetting conditions (Snoeijer and Andreotti, 2013). The different models for microscopic and macroscopic contact angles have some limitations too. While a model which describes the contact angle as a function of the capillary number might describe dynamic contact angle hysteresis correctly, there is also static contact angle hysteresis, which causes contact line pinning of non-moving contact lines. Any model that lets the contact angle only depend on the local capillary number will not be able to capture static contact angle hysteresis. Assuming a static angle ,on the other hand, does not capture static contact angle hysteresis, is not able to properly describe the flow of a liquid over a chemically patterned surface with different wetting properties, and does not accurately predict contact line velocity Yamamoto et al. (2013).

In this work, a validation study of a Volume Of Fluid (VOF) implementation of a contact line model called the Generalized Navier Boundary Condition (GNBC) is presented. While the non-sharp interface of the VOF method implicitly resolves the contact line singularity problem, even with a no-slip condition at the wall (Seppecher, 1996; Chen et al., 2000; Sibley et al., 2013), the question of what is the right contact angle is still valid for this method. In addition to applying the Navier slip condition at the wall, the GNBC uses the reduced Young’s stress as a restoring force when the contact angle deviates from its equilibrium value, and is informed by Molecular Dynamics (Qian et al., 2003; Gentner et al., 2004). While previous implementations used a friction factor to link the reduced Young’s stress to a contact line velocity (Gerbeau and Lelievre, 2009), this approach does not work for the VOF model, because of the implicit slip of the interface (Saha and Mitra, 2009). Instead, in this work, the reduced Young’s stress is used directly as a source term in the navier stokes equations (Mahady et al., 2015). A consequence of using the GNBC is that the contact angle no longer is a constraint imposed onto the system, but that the contact angle is a self-selecting variable. Using this approach has a couple of advantages over existing models. Because there is no enforcement of a model that relates the contact angle to contact line speed, this model can reproduce static contact angle hysteresis, without artificially fixing the position of the contact line (Fang et al., 2008; Dupont and Legendre, 2010). Additionally, a variable contact angle model can describe flow over a chemically patterned surface (Dupuis and Yeomans, 2004; Wang et al., 2008). While the Generalized Navier Boundary Condition has been implemented for continuum simulations using a diffuse interface Cahn Hillard method (Qian et al., 2006; Ren and Weinan, 2007; Wang et al., 2008; Cai et al., 2015), using a Volume Of Fluid approach has the advantage that many less grid points are needed to resolve the interface (Ding et al., 2007). In contrast to the body force term derived by Mahady et al. (2015), the model presented in this work describes wetting of a dry surface. Opposed to Deganello et al. (2011), the model presented here is a analytical derivation of the line tension force.

In our work, we uncover the following findings: when looking at the static behavior of droplets, it is found there is a good match between the droplet shape resolved in the simulations and the theoretically predicted shape for various values of the Young’s angle. Investigating the spreading of a droplet on a completely wetting surface, the Voinov-Tanner-Cox law (Voinov, 1976; Tanner, 1979; Cox, 1986) () can be observed. Late time scaling follows , suggesting spreading is limited by inertia. These observations are made without any fitting parameters except the slip length.

## Ii Theory

### ii.1 Body force reformulation of uncompensated Young’s stress

The traditional form of the Generalized Navier Boundary Condition (GNBC) used in, for example, Arbitrary Lagrangian-Eulerian simulations (Gerbeau and Lelievre, 2009), looks like:

(1) |

Here the first two terms on the left hand side represent the Navier slip boundary condition (Navier, 1823), and the third term represents the unbalanced Young’s stress. is the slip coefficient, is velocity, is the velocity of the wall, is any vector tangent to the wall, is the dynamic viscosity, is the normal of the wall pointing outward, is the surface tension coefficient, is the dynamic contact angle, is the equilibrium Young’s contact angle, is the vector tangent to the wall and normal to the contact line, and the distribution is defined as (Gerbeau and Lelievre, 2009):

(2) |

where is any smooth function, and denotes the Lebesgue measure (i.e. the length measure) on the contact line. The slip length in equation 1 is equal to .

Analogous to the methods of Brackbill et al. (1992) a Volume Of Fluid expression is derived for the above equation for the Generalized Navier Boundary Condition. However, instead of relating the uncompensated Young’s stress to a velocity using a slip coefficient, the uncompensated Young’s stress is modeled as an extra body force acting at the contact line in the Navier-Stokes equation. The reason for this approach is that the intrinsic contact line slip in a VOF code is large enough that converting the contact line tension directly to a velocity on the wall does not move the contact line. The Navier-slip boundary condition component of the above equation is left unchanged in our implementation of the Generalized Navier Boundary Condition.

The goal of this derivation is to find an expression that rewrites the contact line force as a body or volume force, so it can be treated as a momentum source term in the Navier-Stokes equations. The first step is to find a relation, that rewrites the line force as surface force acting on the wall:

(3) |

The points, , form the contact line or triple point, is the length of the line segment in a the small volume of integration , and is the side of that is part of the wall, and in which the points lay. An additional constraint for is that it is zero outside of the interface region:

(4) |

where is the normal to the contact line in the plane of the wall, and describes the plane .

Consider a system of two fluids, fluid , and fluid , separated by an interface, and define a discontinuous function, , to distinguish between the two phases:

(5) |

An example of such a function would be the density of two different incompressible liquids, , and . In this case the position of the contact line can be found with:

(6) |

In order to change this problem from a boundary value problem on the contact line to an approximate continuous model, one can define a continuous function , which varies smoothly over thickness going from to , . is a length scale of the order of the grid size , and defines the width of the transition region from to . The two functions and are related through the interpolation function :

(7) |

which is normalized as:

(8) |

is bounded as:

(9) |

is differentiable, and decreases monotonically with . The continuous function is defined such that:

(10) |

i.e. the function approaches as the interface thickness goes to zero. is differentiable because is, and:

(11) |

where is the two-dimensional gradient in the plane of the wall. Using Gauss’ theorem and the realization that is constant within each fluid, the above integral can be written as:

(12) |

where , thus converting the surface integral to a line integral. To pull the normal out of the integral, its weighted mean is calculated. Since is bounded, its maximum contribution to the line integral is . Integral 12 can thus be approximated as:

(13) |

where is the point on closest to , and is the radius of the contact line.

The integral in equation 13 can be bounded by:

(14) |

where in the limit , is zero everywhere except for . Taking the corresponding limit of over the interface gives:

(15) |

As goes to , is thus equivalent to:

(16) |

Because the Brackbill surface tension already takes care of the component of the reduced Young’s stress, only the component needs to be modeled. Using the delta function, the contact line force can now be written as a surface force as follows:

(17) |

converting the line integral over the contact line into an integral over the wall surface. Equation 16 can be used as an approximation for the delta function when the interface has a finite thickness. Substitution gives:

(18) |

and comparing equation 3 with 18, the surface force can be identified as:

(19) |

As a last step, the surface integral needs to be converted to a volume integral. Since is independent of the distance away from the wall, integrating over is the same as multiplying with mesh size, i.e. :

(20) |

where is the surface area of the wall in volume .

While the above equation has the correct limiting behavior, it was found that, due to the diffuse nature of the interface, at too low resolution the interface gets pulled apart. To counter this phenomenon, a modification is proposed to localize the contact line force more to the interface:

(21) |

Where is the Heaviside step function. Depending on the value of the Young’s angle, this function truncates the contact line force. For this means that the contact line does not get pulled apart, while for this prevents the gas phase from being pulled into the droplet at the contact line. The value is to normalize the function. The new function for the contact line restoring force at low resolution thus becomes:

(22) |

where the same definition is used for as in equation 20.

### ii.2 Numerical implementation

This section focuses on the numerical implementation of the above derived equation in the Volume Of Fluid (VOF) solver that comes with OpenFOAM (Jasak, 1996; Ltd, 2011). This involves implementing the uncompensated Young’s stress and the Navier slip boundary condition. In this code the general phase parameter is called , and has the following properties:

(23) |

Phase parameter is stored as a separate field, just like velocity and pressure, and its evolution is calculated using the following transport equation:

(24) |

where is the phase averaged velocity, and is the velocity difference between the liquid and the gas phase. This equation is equivalent to a material derivative, but rewritten to minimize numerical diffusion Rusche (2002).

The volume fraction is used to calculate phase-averaged densities, velocities, and viscosities, which are used in the momentum balance

(25) |

and the continuity equation:

(26) |

In the above equations is the density, is the velocity, time, is pressure is the viscosity, is gravity, is the surface tension force, is the contact line tension force, and is the dyadic product. The density , velocity , and viscosity , are all phase averaged using .

The surface tension force is calculated using the expression:

(27) |

where is the surface tension coefficient,

(28) |

is the curvature of the interface, and

(29) |

is the normal of the interface Brackbill et al. (1992).

Using , the line tension force defined in equation 22 is rewritten as:

(30) |

and is implemented as:

(31) |

where again is the normal of the wall pointing outward. The equilibrium Young’s angle can be defined uniquely for any grid cell along the wall, so an arbitrary wettability pattern can be created. and are properties of the mesh that can be accessed directly in OpenFOAM. The contact line tension source term is solved explicitly along with the surface tension source term.

The Navier slip condition is implemented as:

(32) |

where is the slip length is the identity matrix, and is the dyadic product of the wall normal, , with itself. Using this formulation, the velocity perpendicular to the wall is always set to zero.

## Iii Results

The first section of the results focuses on the static behavior of the implementation of the Generalized Navier Boundary Condition presented above. The second section covers the dynamic contact angle behavior. To speed up simulation times all results are for 2D droplets in a planar geometry, and using the symmetry of the system only half of the droplet is simulated. All simulations are performed without gravity and represent water droplets in air at room temperature.

### iii.1 Static

Figure 1 shows the interface () of various droplets with different Young’s angles at a resolution of in a box of . The half droplets in the figure have a surface area of about (i.e. a radius of when the Young’s angle is ). As initial condition these droplet where given their equilibrium shape, and it can be seen that after a simulation time of the droplets have maintained their shape.

To further quantify these droplets, figure 2 shows the height and radius of the base of the droplets as a function of the Young’s angle. The radius of the base of the droplet is calculated as:

(33) |

and the height of the droplet is equal to:

(34) |

In the above equations:

(35) |

is the radius of the droplet, is the Young’s angle, and is the surface area of the droplet. While there are small deviations between the theoretically predicted droplet shapes and the simulations, both match well.

The convergence of the error as a function of resolution is shown in figure 3. The x axis shows the resolution. The values of the time steps at these resolutions are: , , , and , to keep the Courant number of the same order between simulations. The y axis shows the absolute error, , between the theoretical value for the radius of the base of the droplet and the simulated value for a Young’s angle of . The graph shows that the code is close to 2nd order accurate.

Because the interface gets thinner as resolution increases the Volume Of Fluid method is mesh dependent. As can be seen in figure 4, this results in the pressure peak on the surface getting sharper with increasing resolution. The increased pressure inside the drop, left of the pressure peak, is the Laplace pressure.

The reason there is a pressure peak at the contact line for a stationary droplet in the first place can be seen in figure 5. This figure shows the value of the pressure peak as a function of the Young’s angle. From the figure it is clear that there is a strong dependence of the pressure peak on the Young’s angle. Since the line tension force only acts parallel to the surface, this suggest that the pressure peak is the result of the surface tension force calculated by the Brackbill et al. (1992) model. If the vertical component of the reduced Young’s stress also were implemented, an extra term, proportional to , would have been present in the plot as an additional contribution to the pressure peak. Since we are only concerned about solid surfaces, this term was not incorporated in the model.

Due to the large pressure peak there is also a small residual slip velocity at the contact line, which can be observed in figure 6. As was the case with the pressure, this velocity peak gets sharper at higher resolutions. As can be seen in figure 2 the residual slip velocity does not negatively affect the shape of the drop. However, spurious currents are a well known issue of Volume Of Fluid solvers Deshpande et al. (2012) and, if needed, can be controlled by making the time step sufficiently small.

To investigate to what extent the pressure and velocity peaks affect the solution they are both integrated over the wall for various resolutions. Figure 7 shows this integral for the pressure, and figure 8 for the slip velocity. The integration limits in both plots are from to . For the integral over the pressure peak this approach makes sure that the integral is not affected by the Laplace pressure inside the droplet, and the velocity integral uses the same limits to be consistent with the pressure. It can be appreciated how the pressure integral converges to a constant value and the velocity integral approaches zero as resolution increases, showing a convergent solution for both the pressure and velocity.

### iii.2 Dynamic

For the validation of the dynamic case, the starting point is again a 2D droplet in planar geometry. However, in this case, the equilibrium Young’s contact angle is set to , and the spreading behavior as a function of time is investigated.

Figure 9 shows the interface () at various times for a box of and with a resolution of . As expected the droplet keeps spreading until the edge of the simulation box is reached and the simulation is stopped.

Figure 10 shows how the droplet spreads by showing the radius of the droplet as a function of time for two different resolutions of and . For lower resolutions the contact line was pulled apart and the results are not shown in the graph. For reference, the power law is also shown. The power law describes the late time spreading behavior for low viscosity axisymmetric droplets (Biance et al., 2004; Ding and Spelt, 2007). At a resolution of initially the contact line hardly moves, but then the spreading radius as a function of time becomes proportional to: . At the larger resolution of it can be seen that the scaling also converges to .

Figures 11 and 12 show the pressure profile and velocity at the wall for various times. The pressure shows a sharp peak already, but especially at there is quite some noise in both the curves for pressure and velocity. The peaks are expected to become even sharper with increasing resolution (Afkhami et al., 2009). While it is known that for the Navier-slip boundary condition pressure is divergent (Devauchelle et al., 2007), whether the reduced Young’s stress provides a cut-off mechanism for the pressure has not yet been analytically determined.

Because of the difference in the spreading curves in figure 10 an additional simulation was performed at a much larger resolution. The simulation domain for this simulation is with a resolution of , but the mesh is refined at the wall to better capture the curvature at the contact line. The grid cells at the wall are about cubed. The droplet has an initial radius of and makes an angle to with the surface. The function which was used in the above simulations to keep the contact line tension localized was omitted in this simulation because it is not needed at larger resolutions.

Figure 13 shows the interface () of the droplet for various times. Figure 14 shows the corresponding radial position of the contact line as a function of time. As was observed in Figure 10, this smaller droplet also shows inertia dominated spreading.

Because of the larger resolution in these simulations the apparent contact angle can accurately be determined. Figure 15 shows the apparent contact angle as a function of the capillary number. It can be appreciated that the simulation recovers the Voinov-Tanner-Cox (Voinov, 1976; Tanner, 1979; Cox, 1986) law: ().

Figures 16 and 17 show the pressure peak and slip velocity at the wall at various times. Because of the high resolution at the wall both are very sharp peaks confined to the interface. The fluctuations of the pressure and velocity that was observed in figures 11 and 12 is no longer there, suggesting the simulations have fully converged.

## Iv Conclusions & Discussion

An implementation of the Generalized Navier boundary condition for the Volume Of Fluid method is presented in this work. In analogy with the Brackbill surface tension model, a body force representation is developed for the contact line tension, while the Navier slip condition is applied on the wall. A validation of the code is presented for both a static case of a droplet maintaining its equilibrium shape, and a dynamic case of a spreading droplet.

It is shown how, on a completely smooth solid surface, in a system without gravity, the shape of the simulated droplets matches with their theoretically predicted shape for various Young’s angles. In addition, it is shown how the pressure peak and corresponding velocity peak at the interface converge with increasing resolution.

For the dynamic case it is found that the spreading of the droplet scales as . This suggests that spreading is limited by inertia. Also the Voinov-Tanner-Cox law is observed (). This behavior is observed without using any fitting parameters.

An interesting future application of this model is flow over pattered surfaces (Dupuis and Yeomans, 2004; Wang et al., 2008). Another topic of interest is the possibility to investigate both static and dynamic contact angle pinning. One can apply any pattern of equilibrium contact angles on a surface, and study how different patterns pin the contact line, either in a static or in a dynamic system.

###### Acknowledgements.

Many thanks to Sidney Nagel and the Nagel group at the Department of Physics of the University of Chicago for stimulating discussions and feedback.## References

- Squires and Quake (2005) T. M. Squires and S. R. Quake, Reviews of modern physics 77, 977 (2005).
- Dussan (1979) E. Dussan, Annual Review of Fluid Mechanics 11, 371 (1979).
- Gennes (1985) P.-G. D. Gennes, Reviews of modern physics 57, 827 (1985).
- Coninck and Blake (2008) J. D. Coninck and T. Blake, Annu. Rev. Mater. Res. 38, 1 (2008).
- Eral et al. (2013) H. B. Eral, J. M. Oh, et al., Colloid and polymer science 291, 247 (2013).
- Furmidge (1962) C. G. L. Furmidge, Journal of colloid science 17, 309 (1962).
- Zhu et al. (2008) X. Zhu, P. Sui, and N. Djilali, Journal of Power Sources 181, 101 (2008).
- Wexler et al. (2015) J. S. Wexler, A. Grosskopf, M. Chow, Y. Fan, I. Jacobi, and H. A. Stone, Soft matter (2015).
- Vandre et al. (2013) E. Vandre, M. Carvalho, and S. Kumar, Physics of Fluids 25, 102103 (2013).
- Young (1805) T. Young, Philosophical Transactions of the Royal Society of London pp. 65–87 (1805).
- Gauss (1830) C. F. Gauss, Principia generalia theoriae figurae fluidorum in statu aequilibrii. (Gottingae: Dieterichs, 1830).
- Huh and Scriven (1971) C. Huh and L. Scriven, Journal of Colloid and Interface Science 35, 85 (1971).
- Hocking (1977) L. Hocking, Journal of Fluid Mechanics 79, 209 (1977).
- Navier (1823) C. Navier, Mémoires de l’Académie Royale des Sciences de l’Institut de France 6, 389 (1823).
- Thompson and Troian (1997) P. A. Thompson and S. M. Troian, Nature 389, 360 (1997).
- Cahn (1977) J. W. Cahn, The Journal of Chemical Physics 66, 3667 (1977).
- Seppecher (1996) P. Seppecher, International journal of engineering science 34, 977 (1996).
- Jacqmin (2000) D. Jacqmin, Journal of Fluid Mechanics 402, 57 (2000).
- Hervet and de Gennes (1984) H. Hervet and P.-G. de Gennes, Comptes-rendus des séances de l’Académie des sciences. Série 2, Mécanique-physique, chimie, sciences de l’univers, sciences de la terre 299, 499 (1984).
- Derjaguin (1955) B. Derjaguin, Kolloid Zhurnal 17, 205 (1955).
- Teletzke et al. (1987) G. F. Teletzke, H. T. Davis, and L. Scriven, Chemical Engineering Communications 55, 41 (1987).
- Herring and Henderson (2010) A. Herring and J. Henderson, The Journal of chemical physics 132, 084702 (2010).
- Cox (1998) R. Cox, Journal of Fluid Mechanics 357, 249 (1998).
- Shikhmurzaev (1993) Y. D. Shikhmurzaev, International Journal of Multiphase Flow 19, 589 (1993).
- Shikhmurzaev (1997) Y. D. Shikhmurzaev, Journal of Fluid Mechanics 334, 211 (1997).
- Chen and Wada (1989) J.-D. Chen and N. Wada, Physical review letters 62, 3050 (1989).
- Eggers (2004) J. Eggers, Physical review letters 93, 094502 (2004).
- Blake and Haynes (1969) T. Blake and J. Haynes, Journal of colloid and interface science 30, 421 (1969).
- Blake (2006) T. D. Blake, Journal of Colloid and Interface Science 299, 1 (2006).
- Voinov (1976) O. Voinov, Fluid Dynamics 11, 714 (1976).
- Tanner (1979) L. Tanner, Journal of Physics D: Applied Physics 12, 1473 (1979).
- Cox (1986) R. Cox, Journal of Fluid Mechanics 168, 169 (1986).
- Kafka and Dussan (1979) F. Y. Kafka and E. Dussan, Journal of Fluid Mechanics 95, 539 (1979).
- Sibley et al. (2014) D. N. Sibley, A. Nold, N. Savva, and S. Kalliadasis, Journal of Engineering Mathematics pp. 1–23 (2014).
- Bonn et al. (2009) D. Bonn, J. Eggers, J. Indekeu, J. Meunier, and E. Rolley, Reviews of modern physics 81, 739 (2009).
- Koplik et al. (1988) J. Koplik, J. R. Banavar, and J. F. Willemsen, Physical review letters 60, 1282 (1988).
- Koplik et al. (1989) J. Koplik, J. R. Banavar, and J. F. Willemsen, Physics of Fluids A: Fluid Dynamics (1989-1993) 1, 781 (1989).
- Thompson and Robbins (1989) P. A. Thompson and M. O. Robbins, Physical Review Letters 63, 766 (1989).
- Thompson and Robbins (1990) P. A. Thompson and M. O. Robbins, Physical Review A 41, 6830 (1990).
- Checco et al. (2006) A. Checco, Y. Cai, O. Gang, and B. M. Ocko, Ultramicroscopy 106, 703 (2006).
- Bocquet and Charlaix (2010) L. Bocquet and E. Charlaix, Chemical Society Reviews 39, 1073 (2010).
- Seveno et al. (2013) D. Seveno, T. D. Blake, and J. D. Coninck, Physical review letters 111, 096101 (2013).
- He and Hadjiconstantinou (2003) G. He and N. Hadjiconstantinou, Journal of Fluid Mechanics 497, 123 (2003).
- Shen and Ruth (1998) C. Shen and D. W. Ruth, Physics of Fluids (1994-present) 10, 789 (1998).
- Hardy (1919) W. B. Hardy, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 38, 49 (1919).
- Kavehpour et al. (2003) H. P. Kavehpour, B. Ovryn, and G. H. McKinley, Physical review letters 91, 196104 (2003).
- Huh and Mason (1977) C. Huh and S. Mason, Journal of fluid mechanics 81, 401 (1977).
- Snoeijer and Andreotti (2013) J. H. Snoeijer and B. Andreotti, Annual review of fluid mechanics 45, 269 (2013).
- Yamamoto et al. (2013) Y. Yamamoto, T. Ito, T. Wakimoto, and K. Katoh, International Journal of Multiphase Flow 51, 22 (2013).
- Chen et al. (2000) H.-Y. Chen, D. Jasnow, and J. Viñals, Physical review letters 85, 1686 (2000).
- Sibley et al. (2013) D. N. Sibley, A. Nold, N. Savva, and S. Kalliadasis, The European Physical Journal E 36, 26 (2013).
- Qian et al. (2003) T. Qian, X.-P. Wang, and P. Sheng, Physical Review E 68, 016306 (2003).
- Gentner et al. (2004) F. Gentner, R. Rioboo, J. Baland, and J. D. Coninck, Langmuir 20, 4748 (2004).
- Gerbeau and Lelievre (2009) J.-F. Gerbeau and T. Lelievre, Computer Methods in Applied Mechanics and Engineering 198, 644 (2009).
- Saha and Mitra (2009) A. A. Saha and S. K. Mitra, Journal of colloid and interface science 339, 461 (2009).
- Mahady et al. (2015) K. Mahady, S. Afkhami, and L. Kondic, Journal of Computational Physics 294, 243 (2015).
- Fang et al. (2008) C. Fang, C. Hidrovo, F. min Wang, J. Eaton, and K. Goodson, International Journal of Multiphase Flow 34, 690 (2008).
- Dupont and Legendre (2010) J.-B. Dupont and D. Legendre, Journal of Computational Physics 229, 2453 (2010).
- Dupuis and Yeomans (2004) A. Dupuis and J. M. Yeomans, Future Generation Computer Systems 20, 993 (2004).
- Wang et al. (2008) X.-P. Wang, T. Qian, and P. Sheng, Journal of fluid mechanics 605, 59 (2008).
- Qian et al. (2006) T. Qian, X.-P. Wang, and P. Sheng, Journal of Fluid Mechanics 564, 333 (2006).
- Ren and Weinan (2007) W. Ren and E. Weinan, Physics of Fluids (1994-present) 19, 022101 (2007).
- Cai et al. (2015) X. Cai, H. Marschall, M. Wörner, and O. Deutschmann, Chemical Engineering & Technology 38, 1985 (2015).
- Ding et al. (2007) H. Ding, P. D. Spelt, and C. Shu, Journal of Computational Physics 226, 2078 (2007).
- Deganello et al. (2011) D. Deganello, T. Croft, A. Williams, A. Lubansky, D. Gethin, and T. Claypole, Journal of Non-Newtonian Fluid Mechanics 166, 900 (2011).
- Brackbill et al. (1992) J. Brackbill, D. Kothe, and C. Zemach, Journal of computational physics 100, 335 (1992).
- Jasak (1996) H. Jasak, Ph.D. thesis, Imperial College London (University of London) (1996).
- Ltd (2011) O. Ltd, OpenFOAM (Version 2.1.1) [Computer software] (2011), URL http://www.openfoam.org.
- Rusche (2002) H. Rusche, Ph.D. thesis, Imperial College (2002).
- Deshpande et al. (2012) S. S. Deshpande, L. Anumolu, and M. F. Trujillo, Computational Science & Discovery 5, 014016 (2012).
- Biance et al. (2004) A.-L. Biance, C. Clanet, and D. Quéré, Physical Review E 69, 016301 (2004).
- Ding and Spelt (2007) H. Ding and P. D. M. Spelt, Journal of fluid mechanics 576, 287 (2007).
- Afkhami et al. (2009) S. Afkhami, S. Zaleski, and M. Bussmann, Journal of computational physics 228, 5370 (2009).
- Devauchelle et al. (2007) O. Devauchelle, C. Josserand, and S. Zaleski, Journal of Fluid Mechanics 574, 343 (2007).