Computers & Fluids 46 (2011) 306–311
Contents lists available at ScienceDirect
Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d
Three-dimensional thin film and droplet flows over and past surface features with complex physics Y.C. Lee a, H.M. Thompson b, P.H. Gaskell b,⇑ a b
Department of Mechanical Engineering, Heriot-Watt University, Edinburgh EH14 4AS, UK School of Mechanical Engineering, University of Leeds, Leeds LS2 9JT, UK
a r t i c l e
i n f o
Article history: Received 30 April 2010 Received in revised form 10 August 2010 Accepted 11 August 2010 Available online 17 August 2010 Keywords: Thin film flow Droplets Mixing Flexible surface Multigrid Topography
a b s t r a c t This paper describes, in the main, the benefits of using a general Newton globally convergent solver within an adaptive multigrid framework for solving discretised forms of lubrication models of threedimensional, free-surface flow over and/or past substrate topography, involving complex physics. The two illustrative gravity-driven problems considered address solute mixing in a continuous thin film flow and droplet migration down an inclined substrate. The computational price paid for the flexibility offered by the solver is investigated alongside the overall benefits of adaptive local mesh refinement and multigridding. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction Three-dimensional thin film and droplet flow on heterogeneous substrates is of enormous significance in engineering and industrial applications, such as heat exchangers, tribological components and precision manufacturing applications [1]. Reliable film thickness/droplet control in such products and processes is often very difficult to achieve in practice since the free-surface disturbances induced by small-scale substrate chemical heterogeneity, patterning or occlusions result in non-uniformities that can persist over length scales several orders of magnitude greater than the size of the features that cause them [2]. The present lack of reliable experimental data is testament to the difficulties of designing such systems experimentally and flow simulation is likely to play an increasingly important role in their interrogation for the foreseeable future. The present work focuses on the accurate and efficient numerical solution of the equations governing such flows which feature complex physics. The numerical prediction of thin film and droplet flows is still at a relative early stage of development since the computational resources required to solve the associated system of Navier–Stokes equations to the required accuracy in three-dimensions are currently lacking. Consequently, numerical studies to date have employed a range of simplifying assumptions to make the problems ⇑ Corresponding author. E-mail address:
[email protected] (P.H. Gaskell). 0045-7930/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2010.08.008
tractable. For example, the neglect of fluid inertia has enabled the Boundary Element method to provide an effective means of investigating film flows over topography, see e.g. [3,4]; however, the majority of studies have employed the long wave approximation. The latter enables three-dimensional film flows to be represented by two-dimensional lubrication equations [5]. In the case of droplet spreading/migration great reliance is placed on the use of a disjoining pressure model to alleviate the stress singularity that arises at the three-phase contact line(s) present in the flow [6,7]. Lubrication equations have modelled successfully a wide range of different thin film and droplet flows and the underlying equations have been solved using a variety of numerical methods. The Finite Element method, for example, is increasingly being used to solve weak forms of the lubrication equations [8]. However by far the most common approach to date has been to solve Finite Difference/volume analogues of them. Until quite recently timesplitting schemes proved popular in this respect but over the last 5 years or so the advantages of utilising an alternative Full Approximation Storage (FAS), Full Multigrid (FMG) strategy have become increasingly apparent [9]. This paper describes the benefits of using a general Newton globally convergent flow solver within a flexible and adaptive multigrid framework for solving the increasingly complex system and number of equations that arise when thin film and droplet flows are subjected to additional physical effects. Its advantages are demonstrated by solving two new three-dimensional, gravity-driven
307
Y.C. Lee et al. / Computers & Fluids 46 (2011) 306–311
free-surface flows, involving the mixing of solute species on rigid and flexible substrates and droplet flow past an occlusion. 2. Problem definition 2.1. Mathematical formulation Consider, for illustrative purposes, Fig. 1(a) and (b) which show the motion of a thin liquid film/droplet, of general height H(X, Y, T), over a substrate which is inclined at an angle h to the horizontal and contains a single splitter occlusion. The liquid is assumed Newtonian and incompressible, with constant density, q, viscosity, l, and surface tension, r. Proceeding as in [7,10] with = H0/L0 small [5], where L0 is the characteristic in-plane length scale and H0 the undisturbed fully-developed asymptotic film thickness or initial droplet height, the resulting generic lubrication equation for h = H/H0 is:
" 3 # ( 3" b #) @h @ h @p @ h @p a Bo ; þ ¼ 2 sin h @t @x 3 @x @y 3 @y
ð1Þ
with the nondimensional pressure, p = P/P0, given by
a 6 2 p ¼ 3 r2 h þ a 61=3 Nh þ b½Bo h cos h PðhÞ; b b
ð2Þ
where the pressure datum is set to zero; a and b are constants, the values of which specify the nature of the problem of interest. Spatial and temporal terms are scaled via (x, y) = (X, Y)/L0 and t = T/T0, respectively. Following several previous studies, see for example [6,7], the disjoining pressure term, P(h), in Eq. (2) is given by
" n m # ðn 1Þðm 1Þ h h PðhÞ ¼ rð1 cos hc Þ ; h h h ðn mÞ
ð3Þ
and is used to alleviate the stress singularity that arises at the associated three-phase wetting line. This disjoining pressure model assumes that a thin, energetically stable wetting layer, of thickness h* = H*/H0, is present across the entire flow domain; the values (n, m) = (3 ,2) are employed. For the case of thin film flow (a = 1, b = 0), Fig. 1(a), the characteristic film height, H0 = (3lQ0/qgsin h)1/3; g is the gravitational constant, Q0 is the volumetric flow rate per unit width, pressure
(a)
P0 = qgL0sin h/2, time T0 = L0/U0, velocity U0 = 3Q0/2H0, b = L0/Lc, where Lc = (rH0/3qgsin h)1/3 is the capillary length; N = Ca1/3cot h measures the influence of the normal component of gravity on the flow profile with the Capillary number, Ca = lU0/r. In the case of droplet flow (a = 0, b = 1), Fig. 1(b), the characteristic droplet height, H0 = hcR0/2 (hc being the droplet’s equilibrium contact angle), pressure, P0 = r/L0, time T0 = lL0/r3, velocity U0 = L0/T0, and length scale L0 = dR0 (d is a ratio relating the size of the domain to the size of the drop); Bo ¼ qgL20 =r is the Bond number representing the ratio of gravity to surface tension forces. 2.2. Boundary conditions The boundary conditions in the case of thin film flow arise from the assumption that the flow is fully developed both upstream and downstream and is tangential at the spanwise boundaries. In the case of droplet spreading/migration, simple symmetrical boundary conditions are employed, as the solution domain is taken to be sufficiently large to ensure that the droplet never reaches one of its edges. When occlusions are present within the flow the following boundary condition is applied at the associated static occlusion wetting line [11], rhw n ¼ 1 tan p2 hs , where hw(x, y) is the static wetting curve, n is the outward pointing unit normal and hs the static contact angle. The remaining condition applied at the occlusion boundary is that of no-flux. 2.3. Incorporating additional physics Thin film flow can be subject to several competing physical influences [1]. For example, for a liquid whose viscosity can be assumed independent of solute concentration, c, and which obeys the well mixed approximation [12], the governing advection–diffusion concentration equation [13] is given by:
" 2 " 2 # # @c h @p @c h @p @c ¼ 2 þ : @t @x 3l @x 3l @y @y
ð4Þ
Substrate flexibility can be accounted for by relating its deflection, e = E/H0, and f = h + e, to the corresponding hydrostatic and capillary pressures exerted by the liquid film via a direct application of Newton’s second law [14]. If the flexible substrate is homogeneous,
(b)
Fig. 1. Schematic of gravity-driven film and droplet flow past an occlusion, or blockage, of length, Lb, and width, Wb, showing the associated geometry defined on a domain of size, LS WS, for the case of: (a) thin film flow over a substrate containing a flexible patch (length, Lf and width, Wf); (b) droplet migration and division. Direction of flow is from upper left to lower right.
308
Y.C. Lee et al. / Computers & Fluids 46 (2011) 306–311
infinitely long and thin (non-dimensional thickness k) with uniform tension, f, in the longitudinal and transverse directions and to have constant density, qm, then within the limits of the lubrication approximation ( 1) the equation to be solved is:
@2e @2e ! þ @x2 @y2
!
q @2f @2f ¼ f þ m k Bo cos h þ 2 þ 2 ; @x @y q
ð5Þ
where = f/r is the ratio of tension in the flexible substrate to the surface tension of the liquid. The additional boundary conditions required to close the problem following incorporation of either, or both, of the above physical effects are: the inlet concentration profile c(x = 0, y) = c0 and e = 0 in regions where the substrate is rigid.
using a combined Full Approximation Storage (FAS) and Full Multigrid (FMG) technique. Relaxation/smoothing on grid Gk is performed using a fixed number of pre- and post- Red-Black GaussSeidel Newton iterations. 3.3. General Newton globally convergent solver With u representing the vector of unknowns of the equation set to be solved, that is in the present work Eqs. (1) and (2) together with either one or both of Eqs. (4) and (5), the same can be expressed more generally as: n Nk ðunþ1 k Þ ¼ Fk ðuk Þ:
ð6Þ
The system of residual equations to be minimised can then be written as:
3. Method of solution Below only the new, key features of the proposed general numerical procedure are given, with just an overview provided of the other salient points since they are reported extensively elsewhere. 3.1. Spatial and temporal discretisation
n gk ðuk Þ ¼ Nk ðunþ1 k Þ Fk ðuk Þ ¼ 0;
ð7Þ
in which uk denote the vector {ui,k} of unknowns. Employing a Newton approach, based on forward-differences, and neglecting terms of order (Duk)2 and higher, leads to:
gk ðuk þ Duk Þ ¼ gk ðuk Þ þ Jk Duk ;
The above equations, written in discrete form, are solved at each node of a rectangular computational domain. Use of the Crank–Nicholson method to approximate the time-derivative terms appearing in Eqs. (1) and (4) together with central differencing for spatial derivatives, leads to discrete analogues that are second order accurate in both space and time [9,13,14]. 3.2. Multigrid strategy Following the multigrid approach employed in [9], a sequence of progressively finer grids ðGk : k ¼ 0; 1; . . . ; KÞ, with uniform grid spacing Dk, is defined. Each grid level, Gk , has nk ¼ 2kþkc þ1 þ 1 nodes per unit length in both co-ordinate directions where kc is a constant defining the resolution of the coarsest grid level, such that the mesh size associated with it and that on finer grid levels is Dk ¼ 2ðkþkc þ1Þ . The corresponding time-dependent, nonlinear, coupled set of governing lubrication equations together with any extra equations embodying additional physics are solved efficiently
ð8Þ
where gk(uk + Duk) = 0. This yields the following system of linear equations for the Newton step:
Duk ¼ J1 k gk ðuk Þ;
ð9Þ
which is solved by LU decomposition and the corrections Duk are then added to the solution vector. The associated Jacobians, Jk , are obtained numerically, thus circumventing: (i) the need to derive them analytically [15] followed by hard coding the same each time an equation representing a new physical effect is added to the system, or an exiting equation modified; (ii) situations for which analytic forms of a Jacobian cannot be easily derived or do not exist. Accordingly, each derivative term of the particular Jacobian is computed numerically using forward differencing; the Newton approach employed relies on the assumption that the solution is relatively smooth. The strategy adopted to ensure convergence is to require that the Newton step decreases the product:
c f
y
x
(a)
y
x
(b)
Fig. 2. Three-dimensional thin film flow past a splitter occlusion of length, lb = 0.5, and width, wb = 0.1, showing the associated automatic adaptive mesh refinement at steady-state for: (a) free-surface disturbance; (b) concentration distribution (c = 1 and 0 defined at the inlet boundary along x = 0). Direction of flow is from bottom left to top right.
309
Y.C. Lee et al. / Computers & Fluids 46 (2011) 306–311
nk ¼
1 2 1 J ðDuk Þ2 ¼ gk ðuk Þ gk ðuk Þ; 2 k 2
MG (h, p) MG (h, p, c) MG (h, p, c, e)
ð10Þ
CPU time, t
which at the same time provides a descent direction for nk, such that its gradient satisfies:
rnk Duk ¼ ðgk ðuk Þ Jk Þ ðJ1 k gk ðuk ÞÞ ¼ gk ðuk Þ gk ðuk Þ < 0: ð11Þ
1000
100
When Duk is small, the system of equations is solved using a quadratically convergent full Newton step. 10 10000
3.4. Spatial and temporal adaptivity
100000
10000000
Number of unknowns, N Another key component of the solution strategy is the inclusion of both error controlled spatial and temporal adaptivity, the former automatically determines where fine grids are needed to capture details of a rapidly evolving flow. Adaptive grid refinement is implemented via a relative local truncation error skk1 P e, where e is a user-specified tolerance; large values of skk1 indicate regions of significant error between solutions on successive grid levels so that further local mesh refinement is required. For further details the reader is referred to [9]. Automatic adaptive time-stepping, on the other hand, is achieved by employing a temporal error control algorithm based on predictor–corrector stages [7]. The method employed provides an implicit and second order accurate approach by using time-stepping based on local error estimates, obtained from the difference between the current and predicted solutions and which act as an indicator of whether to increase or decrease the time step in a controlled manner. This provides an efficient means of minimising the computational expense associated with repeated time step failure. 4. Results and discussion The first problem explores gravity-driven (h = 30°) thin film flow involving the mixing of solute species on rigid and flexible substrates containing a blunt (semi-circular) nosed splitter occlusion. In all cases, the solution is integrated forward in time to steady-state from an initial planar free-surface, where f = 1, and hs ¼ p2 .
Fig. 4. Comparison of CPU time (s) dependence on the total number of freedoms for flow past a splitter occlusion for the case of: (i) fluid flow only, MG (h, p); (ii) fluid flow with concentration gradients, MG (h, p, c); and (iii) fluid flow with concentration gradients and substrate flexibility, MG (h, p, c, e).
The liquid involved (water: viscosity l = 0.001 Pa s, density q = 1000 kg m3, surface tension r = 0.07 N m1) with Q0 = 1.635 106 m2 s1, has an asymptotic film thickness H0 = 100 lm such that N = 0.12, indicating that the normal component of gravity has negligible effect on the free-surface shape. Unless stated otherwise a total of k = 4, W(2,2)-multigrid levels (kc = 4) are used, the coarsest and finest levels corresponding to D0 = 1/32 and D4 = 1/ 512, respectively, with automatic mesh adaptivity carried out at the final three finest grid levels (e = 0.1jgk(uk)j). Fig. 2(a) shows the free-surface disturbance resulting from flow past a splitter occlusion of length lb = 0.5 and width wb = 0.1, centered at (xs,ys) = (0.75, 0.5), placed in a square (ls = ws = 1) domain having side walls (no-flow boundaries) along y = 0 and y = 1. In this problem Lc = 0.78 mm and b = 39.15. A characteristic horse-shoe bow-wave is formed upstream of and surrounding the splitter occlusion, together with a greatly elevated free-surface at the static wetting line. The maximum free-surface disturbance occurs at the front of the occlusion and is 34% greater than H0; it decreases along the side of the occlusion such that at the exit of the domain the film thickness adjacent to the occlusion has reduced significantly to be 13.5% greater than H0, inline with the results observed in [11].
f
1.35
film height, f
1.3
e
y
x
1.25 1.2
Υ = 2.5 Υ=5 Υ = 10 Υ = 20 Υ=∞
1.15 0.44
0.46
0.48
0.5
0.52
0.54
0.56
y
(a)
(b)
Fig. 3. Three-dimensional thin film flow past the splitter occlusion of Fig. 2 with a square flexible patch located upstream of it (lf = wf = 0.25 centered at (0.25, 0.5)) at steadystate showing: (a) the free-surface disturbance and associated substrate deflection; (b) the effect on film height at the static wetting curve around the splitter occlusion when the flexible patch has different values of . Direction of flow is from bottom left to top right.
310
Y.C. Lee et al. / Computers & Fluids 46 (2011) 306–311
Fig. 5. (a) Three-dimensional droplet flow down an inclined plane encountering a splitter occlusion of length, lb = 1.0, and width, wb = 0.1, centered at (1.5, 0.5) and the associated adaptive mesh refinement at a particular instant. Direction of flow is from top left to bottom right. Top view, (b) and (c), of droplet flow down an inclined plane past a splitter occlusion, shown at progressive time intervals: h = 30° (left); h = 60° (right). Direction of flow is from left to right.
Y.C. Lee et al. / Computers & Fluids 46 (2011) 306–311
The addition of solute transport effects, via Eq. (4), with the solute concentration at the inlet (x = 0) set as c = 0 for 0 6 y 6 0.5 and c = 1 for 0.5 < y 6 1.0, reveals that mixing occurs upstream of the splitter occlusion, see Fig. 2(b). It shows that solute mixing is greatest when the gradient in c is large and diminishes as it flows downstream until it encounters the splitter to produce uniform concentration levels on either side of the channel outlet (x = 1.0). In the example considered, where the splitter is placed at the centre of the domain, the uniform concentration levels at outlet are 0.655 (0 6 y 6 0.45) and 0.345 (0.55 6 y 6 1) respectively, and are equally displaced from the average concentration c = 0.5. Relocating the splitter occlusion in the spanwise direction leads to significant changes in the outlet concentration levels on each side of it. Note that in order to aid visualisation of the adaptive mesh structures, only two mesh refinement levels are shown. Both plots in Fig. 2 show that the mesh refinement strategy has successfully identified the regions with sharp gradients in film thickness and/ or solute concentration that require additional mesh refinement. The complexity of the problem is now further increased by coupling in Eq. (5) and replacing a section of the rigid substrate upstream of the occlusion with a square flexible patch of size lf = wf = 0.25, centered at (xf, yf) = (0.25, 0.5), with k = 0.1, qm = 1000 kg m3 and = {2.5, 5, 10, 20, 1}. The resultant free-surface disturbance and associated substrate deflection for = 2.5 is illustrated in Fig. 3(a) and shows a significant reduction in the freesurface level above the patch. In this case the automatic mesh adaptivity scheme introduces additional mesh refinement over and around the patch, compared with the rigid substrate case shown in Fig. 2(a). Fig. 3(b) shows that varying the tension ratio has a significant effect on the resultant free-surface profile around the splitter occlusion. For example, decreasing from 1 (rigid substrate) to 2.5 produces a 3% increase in the maximum free-surface level and a free-surface depression 13% lower than the asymptotic film thickness on the top of the patch (emax = 0.14). The flexible patch has a comparatively smaller influence on the distribution of solute throughout the domain. For example, for = 2.5 the solute concentration at the outlet x = 1 for 0 6 y 6 0.45 and 0.55 6 y 6 1 is 0.65 and 0.35, respectively, indicating only a very small improvement in mixing from the rigid substrate case reported in Fig. 2, with corresponding values 0.655 and 0.345, respectively. Fig. 4 illustrates the performance of the general Newton globally convergent solver. It shows the CPU time needed to solve the above problem with: (i) the fluid flow Eqs. (1) and (2) only; (ii) the fluid flow and concentrations Eqs. (1), (2) and (4); and (iii) the fluid flow, concentration and flexible substrate Eqs. (1), (2), (4) and (5). It is clear that the extra cost from solving equations embodying additional physics, without the necessity of having to derive a priori and hard-code Jacobians [13,15], is small and that the solver retains the same attractive O(N) efficiency in each case, where N is the number of unknowns. The final problem considers gravity-driven droplet flow down a flat, inclined, rigid substrate containing a splitter occlusion. The droplet properties are taken from the experimental study [16] and are those of a paraboloid silicon oil drop (46V10) of volume 8.3 mm3 with q = 924 kg m3, l = 0.00915 Pa s and r = 0.0205 N m1. An equilibrium contact angle, hc = 45°, is employed which results in H0 = 0.934 mm and R0 = 2.378 mm and the droplet flows down a rectangular domain (ls = 3ws; Lc = 4 mm; b = 1) with the splitter occlusion, lb = 1.0, wb = 0.1, centered at (xs, ys) = (1.5, 0.5). Solutions are obtained using a 25 9 (kc = 2) coarse grid, a total of k = 7, W(2,2)-multigrid levels and automatic adaptive mesh refinement performed on the final five finest grid levels (e = 0.1jgk(uk)j), as shown in Fig. 5(a). The precursor film thickness is given a value H* = 9.43 lm (h* = 0.01); the effect of the value of h* on droplet spreading/migration is discussed in detail in [17].
311
When h = 30°, Fig. 5(b) shows that the droplet, which is round, divides into two equal parts on either side of the occlusion. A small amount of fluid is retained at its front and the droplet halves eventually re-combine, beyond the splitter, again leaving a small amount of fluid on the downstream end of the latter. Increasing h to 60°, Fig. 5(c), leads to the formation of a cusp as the ovalshaped droplet encounters the upstream end of the splitter. In the corresponding case without an occlusion the emergence of a series of sub-droplets, or pearlings, has been observed experimentally [16] and predicted computationally [17], and similar behaviour is observed along the length of the splitter: the droplet divides into five smaller droplets on either side. The latter eventually flow past the occlusion and re-combine with the result that two droplets migrate down the slope and a pearling transition begins once again [16]. 5. Conclusion A general Newton globally convergent lubrication flow solver, implemented within an efficient adaptive multigrid framework, that is capable of solving supplementary coupled equations embodying additional physics without requiring the derivation of awkward analytical Jacobians, is presented and used successfully to solve gravity-driven flow problems of varying complexity. The additional cost involved from incorporating the solution of additional lubrication-like equations embodying additional physics is insignificant and the attractive O(N) efficiency of multigrid schemes in general is preserved. The coupled Multigrid scheme as presented offers the accuracy, efficiency and convenience desired to explore a wide range of practical micro-fluidic film and droplet spreading flows. References [1] Craster RV, Matar OK. Dynamics and stability of thin liquid films. Rev Mod Phys 2009;81(3):1131–97. [2] Dècre MMJ, Baret JC. Gravity-driven flows of viscous liquids over twodimensional topographies. J Fluid Mech 2003;487:147–66. [3] Blyth MG, Pozrikidis C. Film flow down an inclined plane over a threedimensional obstacle. Phys Fluids 2006;18(5):052104. [4] Baxter SJ, Power H, Cliffe KA, Hibberd S. Three-dimensional thin film flow over and around an obstacle on an inclined plane. Phys Fluids 2009;21(3):032102. [5] Oron A, Davis SH, Bankov SG. Long-scale evolution of thin liquid films. Rev Mod Phys 1997;69:931–80. [6] Schwartz LW, Eley RR. Simulation of droplet motion on low-energy and heterogeneous surfaces. J Colloid Interf Sci 1998;202:173–88. [7] Gaskell PH, Jimack PK, Sellier M, Thompson HM. Efficient and accurate time adaptive multigrid simulation of droplet spreading. Int J Num Meth Fluids 2004;45(11):1161–86. [8] Becker J, Grun G, Seemann R, Mantz H, Jacobs K, Mecke KR, Blossey R. Complex dewetting scenarios captured by thin-film models. Nat Mater 2003;2(1): 59–63. [9] Lee YC, Thompson HM, Gaskell PH. An efficient adaptive multigrid algorithm for predicting thin film flow on surfaces containing localised topographic features. Comput Fluids 2007;36(5):838–55. [10] Gaskell PH, Jimack PK, Sellier M, Thompson HM, Wilson MCT. Gravity-driven flow of continuous thin liquid films on non-porous substrates with topography. J Fluid Mech 2004;509:253–80. [11] Sellier M, Lee YC, Thompson HM, Gaskell PH. Thin film flow on surfaces containing localized arbitrary occlusions. Comput Fluids 2009;38(1):171–82. [12] Howison SD, Moriarty JA, Ockendon JR, Terril EL, Wilson SK. A mathematical model for drying paint layers. J Eng Math 1997;32:377–94. [13] Gaskell PH, Jimack PK, Sellier M, Thompson HM. Flow of evaporating gravitydriven thin liquid films over topography. Phys Fluids 2006;18(1):013601. [14] Lee YC, Thompson HM, Gaskell PH. Thin film flow over flexible membranes containing surface texturing: bio-inspired solutions. J Eng Tribo IMechE Part J 2009;223(3):337–45. [15] Lee YC, Thompson HM, Gaskell PH. FILMPAR: a parallel algorithm designed for the efficient and accurate computation of thin film flow on functional surfaces containing micro-structure. Comput Phys Commun 2009;180(12):2634–49. [16] Podgorski T, Flesselles J-M, Limat L. Corners, cusps and pearls in running drops. Phys Rev 2001;87(3):036102. [17] Koh YY, Lee YC, Gaskell PH, Jimack PK, Thompson HM. Droplet migration: quantitative comparisons with experiment. Euro Phys J – Special Topics 2009;166:117–20.