Asymptotic analysis of axisymmetric drop spreading

We study in this paper the time evolution of the spreading process of a small drop in contact with a ﬂat dry surface, using asymptotic techniques. We reduced the problem by solving a quasisteady self-similar macroscopic problem and matched with the precursor region solution, where the van der Waals forces are important. A ﬁnal nonlinear third-order ordinary differential equation has been solved numerically using shooting methods based on the fourth-order Runge-Kutta techniques. We obtained how the capillary number changes when the drop size decreases with time. The evolution process then diverges slightly from that obtained using the spherical cap approximation. The inﬂuence of gravity is also considered for both hanging and sitting drops. (cid:64)


I. INTRODUCTION
The spreading process of relatively small drops over a solid surface is very important from both the practical as well as the theoretical points of view ͓1͔. Many experimental works have been published in the literature in the past three decades. General reviews in this subject can be found elsewhere ͓2,3͔. Most experiments have been made with wetting fluids, i.e., with a positive spreading parameter Sϭ SG Ϫ SL ϪϾ0, where SG , SL , and are the solid-gas, solid-liquid, and liquid-gas free energy per unit area, respectively. is also called the surface tension. In most of the experiments, the drop radius R(t) is measured using photolithographic techniques, while in others the authors used interferometers. Experimentally it is found that the drop radius increases with time in the form RϳCt m , where different values for m ͑close to 0.1) have been deduced ͓4-6͔. Using a carefully designed experimental setup, Tanner ͓7͔ obtained for silicon oils, mϭ0.105. For small drops where gravity is not relevant, the macroscopic shape of the drop is found to be very close to a spherical cap ͓1͔. However, this small difference between the real shape and the spherical cap is very important in determining the time evolution of the drop spread. Even with small drops, gravity effects tend to be important as the drop radius increases with time.
Another controversial issue in this problem is related with the constant C. Using very simple arguments assuming a small drop shape as a spherical cap, it is possible to derive a very simple theory for the spreading process ͓8͔. The surface tension generated pressure gradient must be balanced by the viscous force, i.e., v/h 2 ϳh/R 3 , where is the viscosity of the fluid, v is the radial velocity of the fluid, and h is the thickness of the fluid drop. It follows that v/ϭCa ϳ(h/R) 3 ϳ 3 . Here is the angle of the surface profile assumed to be very small compared with unity and Ca is the usual capillary number. This is the so-called Tanner's law ͓7,9,10͔ written as ϭC T Ca 1/3 , where C T is the Tanner's constant. Replacing v by dR/dt and h by V/C V R 2 where V is the volume of the drop and C V a constant of order unity, we obtain from Tanner's law a first order ordinary differential equation for R(t), which gives the well known asymptotic behavior Rϳt 1/10 for t→ϱ. However, the constant C T is not universal and depends globally on the specific problem ͓11͔. Friz ͓12͔ obtained values for C T close to 3.4, Chen and Wada ͓13͔ obtained values for C T close to 3.6, while Diez et al. ͓5͔ obtained values around 3.8. Kalliadasis and Chang ͓14͔ obtained a value of 7.48 for the case of an advancing gas-liquid meniscus. Although Tanner's law is universal, the constant depends strongly on the physical configuration. Therefore the value of the constant obtained for the drop spread ͑close to 3͒ is very different from that obtained for an advancing meniscus. We anticipate here that our problem has an inflection point close to the edge of the drop, which does not appear in the configuration studied in ͓14͔. In the case of a spreading drop, the macroscopic region is governed by a nonlinear third order differential equation while in the meniscus problem the macroscopic problem is represented by an algebraic equation. In this paper we show how C T depends on the ratio of the van der Waals influence length a to be defined later to the actual drop size. Due to the changing drop size with time, C T also changes with time, modifying the drop size evolution.
From the analytical point of view, the spreading process including the disjoining pressure effects ͑van der Waals forces͒ has been studied by Starov et al. ͓15͔ and Chebbi and Selim ͓16͔. They used a quasisteady similarity solution. Starov employed the spherical cap approximation while *Author to whom correspondence should be addressed.
Chebbi and Selim obtained numerically the appropriate drop shape. In both cases, matching conditions are assumed to be at the inflection point. Gravity effects have not been included in the above-mentioned works. Brenner and Bertozzi ͓11͔ studied the stability of the quasisteady solution, deducing a simple criterion for the transition to gravity induced spread, using the experimental data obtained by ͓4͔. López et al. ͓17͔,Ehrhard ͓18͔,Voinov ͓20͔ and Chebbi and Selim ͓19͔ studied the spreading process with gravity. The disjoining pressure effects were not considered here.
The main purpose of this paper is to obtain the time evolution of the drop radius and the constant C T for the parametric set of interest, retaining the disjoining pressure and gravity effects.

II. FORMULATION
Using the lubrication approximation, the equation for the evolution of a free surface of a fluid under gravity and capillary forces is given by ͓8͔ where r is the radial coordinate for the assumed axisymmetrical drop, is the fluid density, g is the gravity acceleration, aϭͱ͉A͉/6 is the characteristic length of the nonretarded van der Waals force, which is only a fraction of an angstrom and A is the Hamaker constant, which is negative for a wetting liquid. A positive value of g indicates that the spreading process is taking place over the solid surface. Within this formulation, the averaged radial fluid velocity is then given by v͑ r,t ͒ϭ We introduce the following nondimensional variables: where HϭH 0 G() and RϭR 0 F() are the thickness at the center of the drop and the macroscopic radius of the drop, with H 0 and R 0 being the corresponding values at time t ϭ0. Using these nondimensional variables, the evolution equation ͑1͒ then transforms to where is the ratio of the van der Waals length to the size of the drop, ϭ 0 F/G 2 , which is assumed to be very small compared with unity, with 0 ϭaR 0 /H 0 2 being its value at time tϭ0. BϭgR 2 / corresponds to the Bond number, which relates the gravity to the surface tension forces. Here, BϭB 0 F 2 , with B 0 being the corresponding value at time t ϭ0. In nondimensional form, the radial averaged velocity is then At the macroscopic edge of the drop, the nondimensional radial velocity KϭK f is therefore where K is related to the usual capillary number by K ϭ3Ca(R 0 /H 0 ) 3 , with Caϭv f /. The total volume of the drop that remains invariant during the spreading process is written as To solve Eq. ͑4͒ with the corresponding boundary and initial conditions, we divide the problem in a macroscopic ͑surface tension-viscous-gravity͒ region where /ӷ1, and a thin region of order close to the edge of the drop, (1Ϫ)ϳ, where the effect of the van der Waals forces must be considered in the analysis. Due to the disparity in the two spatial scales, →0, the solution in both regions is to be obtained and properly matched.

A. Macroscopic region †"1؊…ӷ… ‡
Assuming a quasisteady self-similar solution to the macroscopic problem ͑where the effect of the van der Waals forces can be neglected͒, ϭ(), from the overall volume conservation ͑7͒ it follows that F 2 Gϭ1 and Vϭ2R 0 2 H 0 I. In this case, the macroscopic equation ͑4͒ reduces to to be solved with the boundary conditions to be solved with the first three boundary conditions in Eq. ͑9͒. From Eq. ͑5͒ it follows that the averaged radial velocity related to its value at the edge drop increases linearly with the radial coordinate, K r ϭK. There is a family of solutions for any value of K(B,). However, there is only one value of K that satisfies the additional boundary condition arising from matching with the inner solution of the precursor region. This in fact represents a well posed eigenvalue problem.

B. Inner or precursor region †"1؊…ϳ ‡
At the edge of the drop, where is of order , there is a very thin region with (1Ϫ)ϳ, where the nonretarded van der Waals forces cannot be neglected. We introduce for this region the following inner variables of order unity: The inner equation takes the universal form Here yЈϭdy/dx. Matching requires that ͓1͔ y→Ϫx͓3 ln͑Ϫx͔͒ 1/3 for x→Ϫϱ, ͑13͒ in order to assure that in the precursor film y→0 as x→ϱ. In terms of the outer variables, this matching condition transforms to for 1ӷ͑1Ϫ ͒ӷ. ͑14͒

III. RESULTS
Due to the appropriate matching condition, Eq. ͑14͒, we note that the eigenvalue K depends on and B and is therefore a function of time, because the thickness of the drop related to the van der Waals length decreases with time. However, this relationship is weak as we show below. We transform Eq. ͑10͒ to a parameter-free nonlinear equation given by where ϭK 1/4 and B* is a reduced Bond number given by B*ϭB/K 1/2 . The boundary conditions now take the form together with the matching condition ͑14͒. We integrate numerically Eqs. ͑15͒ and ͑16͒ using a fourth-order Runge-Kutta equation with an initial guess of d 2 /d 2 ͉ 0 until the matching condition is fulfilled. Instead of using matching condition ͑14͒, we continued the numerical integration into the inner region given by Eq. ͑12͒, as the macroscopic value of reached 100. The appropriate solution is obtained as we get the final condition in the precursor film, y→0 for x →ϱ. We used a step size of ⌬ϭ10 Ϫ10 for the macroscopic region and ⌬xϭ10 Ϫ3 for the inner region. Figure 1 shows the eigenvalue or reduced capillary number K, as a function of , for three different values of the reduced Bond number B*. Note the very weak dependence with , doubling the value of K in four decades in . Due to the fact that ϭ 0 F 5 () and BϭB 0 F 2 , then K is also a function of time. The lower curve, represented by solid squares, corresponds to the case of a spreading hanging drop with B*ϭϪ5, while the upper curve, represented by solid triangles, corresponds to the case of B*ϭ5. The middle curve neglects the effect of gravity. In this figure it is shown how gravity enhances the spreading rate over a solid. Figure  2 shows the volume integral I as a function of , obtained numerically. For Bϭ0, we can see that I is almost a constant, changing only in the third digit as changes in four decades. The value of I is very close to 0.25, which represents the shape of a spherical cap. A similar behavior for B* 0 is obtained. Therefore, in both cases, it shows that depends mainly on with ‫,0→ץ/ץ‬ thus justifying a quasisteady self-similar solution to the macroscopic Eq. ͑8͒ by neglecting ‫ץ/ץ‬ in Eq. ͑4͒. For B*Ӷ1, we use the fact that can be written by the spherical cap shape plus a small perturbation in the form ͓20͔  with g 2 ϳ1 and g 1 Ӷ1, but ͉d m g 1 /d m ͉ for mу1, could have values very large compared with unity. Here, g i (i ϭ1,2) vanishes at ϭ0 and ϭ1. Neglecting the effect of g i () in 2 in Eq. ͑10͒, we obtain a very good approximative solution for g i as However, with this approximation we are still unable to obtain in closed form the appropriate matching condition with the inner region, because we cannot neglect the effect of g 1 in 2 close to the edge of the drop. Close to →1, 1Ϫ 2 ϳ2(1Ϫ), while g 1 ϳ(1Ϫ)ln(1Ϫ). Nevertheless, we can extract useful information within this approximation. For Bϭ0, the integral I can be solved in closed form as which is also plotted in Fig. 2. This approximation for , given by Eqs. ͑17͒ and ͑18͒, also shows the existence of an inflection point close to the edge of the drop, where the second derivative of vanishes. Although does not change appreciably with time, the gradients change in a very important way with . Figure 3 shows the drop shape for ϭ10 Ϫ10 . For Bϭ0, it is impossible to note any difference for the drop shape for any other value of or that obtained using Eq. ͑17͒. The influence of gravity is weak regarding the drop shape. Figure 4 shows the gradients d/d for the macroscopic region, as a function of . This value is given at the inflection point and also close to the wall, indicating the existence of two different macroscopic apparent contact angles. Figure 4͑a͒ shows these gradients for the case of vanishing gravity. The solution for the gradient of at the inflection point obtained from Eq. ͑17͒ is also plotted in this figure, showing a high accuracy of this approximation. The macroscopic angle at the inflection point decreases with , while the angle close to the wall increases. This shows that the inflection point moves closer to the wall as increases. For B*ϭ5, Fig. 4͑b͒ shows that for increasing values of , the corresponding gradients at the inflection point decrease. This behavior is the opposite to the case with negative values of B*. In all cases, the inflection point moves closer to the wall as increases. For the axisymmetric geometry, our results are different from that obtained using the matching condition at the inflection point. We write the dependence of K on and B as K ϭK 0 ⍀(F 5 )⌫(B 0 F 2 ), with K 0 ϭK( 0 ,B 0 ) and ⍀(1) ϭ⌫(1)ϭ1. Introducing the following variables,  to be solved with the initial condition ⌿(0)ϭ1. Figure 5͑a͒ shows ⍀ as a function of ⌿ for two different values of 0 . In all cases a good correlation of the form ⍀Ӎ1ϩb ln ⌿ ϩc ln 2 ⌿ can be constructed. For 0 ϭ10 Ϫ10 , we obtained b ϭ0.01991 and cϭ0.008121. Figure 5͑b͒ shows similarly ⌫ as a function of ⌿ also for two different values of 0 . Surprisingly there is not a big difference between both curves, which can be very well correlated by ⌫Ӎ1ϩaB 0 ⌿ 2/5 , for all values of . In our case aӍ0.17. For comparison with previous results, for Bϭ0, Eq. ͑22͒ can be readily numerically solved, showing the results in Fig. 6. The well known solution, considering a constant ,⌿ϭ(1ϩ) 1/2 , is also plotted, showing a small but noticeable discrepancy as increases. For →ϱ, the asymptotic solution gives ⌿ϳ n with n plotted in Fig. 7. In physical units, the asymptotic behavior of the drop radius is given by When including the effect of gravity, the asymptotic behavior of ⌿, for large values of nondimensional time , is given by solving numerically the integral ͵ ⌿ s 3/5 ds ln 2 ͑ s ͒ ϳB 0 for →ϱ.
The transition from surface tension induced spread to gravity induced spread can be obtained by solving numerically Eq. ͑22͒. However, this transition can be well visualized by assuming a constant in the spreading process. In this case, Eq. ͑22͒ reduces to to be solved with the initial condition ⌿(0)ϭ1. The solution can be readily obtained by which is the well known behavior of gravity induced spread (Fϳ 1/8 ). One possible criterion for the transition could be that the second term at the right-hand side of Eq. ͑24͒ be of order unity, that is, ⌿ϳ(1/0.17B 0 ) 5/2 . This would correspond to an actual Bond number of BϭB 0 ⌿ 2/5 ϭ5.8824. However, due to the very low exponent of ⌿,2/5, the influence of gravity is felt long before. This means that the transition from surface tension induced spread to gravity induced spread is not sharp enough, as pointed out in ͓6͔. A very simple criterion is then when BϭB 0 ⌿ 2/5 ϭ1 or ⌿ϳB 0 Ϫ5/2 , that is, when the radius of the drop equals the capillary length. Figure 8 shows the transition for the case of B 0 ϭ0.2. In this figure, the actual Bond number BϭB 0 F 2 , which is the same as the ratio of the drop radius to the capillary length, is plotted as a function of the nondimensional time . The surface tension induced spread approximation given by Eq. ͑26͒ as well as the gravity induced spread asymptote, Eq. ͑27͒, are also plotted. Equation ͑26͒ gives a better approximation for values of Bр7, that is, after six decades in the time coordinate counted after reaching a linear profile in the log-log plot (ϳ10 2 ). Experimental results reported elsewhere ͓4,6͔ show the data for only three decades in time. Therefore, it is not possible to show the full transition from surface tension to gravity induced spread with the published data.
The quasisteady self-similar approach is also justified for small values compared with unity, of the ratio of the transit time, t t ϳR/v, to the evolution time, t e ϳK 0 /(dK/dt) or t e ϳB 0 /(dB/dt). Using the above results, these ratios can be written as Therefore, in both cases the quasisteady behavior is fully satisfied due to the fact that ⌿Ͼ1. Only at early times (⌿ ϳ1) the quasisteady approximation fails. It is also possible to obtain the constant of the Tanner's law, C T ϭ/Ca 1/3 , which in our notation can be written as C T ϭ(3/K) 1/3 d/d. This law is rather universal and relates the local macroscopic contact angle with its velocity of a fluid interface, except for the constant, which must be obtained from the global problem. Due to the existence of two different macroscopic angles of the drop, we also plotted Tanner's constant as a function of , using both angles and given in Fig. 9. When using the macroscopic angle close to the wall we note that the Tanner's constant does not change appreciably with , and is very close to 3.0, for all values of B. However, when using the macroscopic angle at the inflection point, it changes practically from 5.28 at ϭ10 Ϫ10 to 3.39 at ϭ10 Ϫ4 . The reported values are between the two curves, showing the difficulty associated with the measuring techniques of the appropriate macroscopic angle.
By way of illustration we computed the evolution Eq. ͑22͒ for a typical oil with the following data ͓21͔: Ϸ841 kg/m 3 , Ϸ0.0225 kg/ms and Ϸ0.035 kg/s 2 . The initial volume of the drop is Vϭ0.024 cm 3 . Under these conditions, the initial Bond number is close to 1.45 and the initial radius of the drop is 0.248 cm. Figure 10 shows the numerical results for the drop radius as a function of time. For comparison we also plot the classical solution Fϭ(1 ϩ) 1/10 in physical units. We can see in this figure that the influence of gravity is extremely important for estimating the evolution process.
In this paper we obtained, using asymptotic techniques, the evolution equation for the spread of a small drop of a wetting fluid over a dry surface. Using the disparity of the spatial scales between the drop size and the van der Waals length, it is possible to reduce the problem by solving a quasisteady self-similar macroscopic problem and matched with the precursor region solution, where the van der Waals forces are important. A final nonlinear third-order ordinary differential equation has been solved numerically using shooting methods based on the fourth-order Runge-Kutta techniques. We obtained that the averaged film velocity increases linearly with the radial coordinate and the radius of the drop increases asymptotically with time in the form R ϳt m with m slightly larger than 1/10 and changing very slowly with time.

ACKNOWLEDGMENTS
C.T. acknowledges the Comisión de Investigaciones Científicas de la Provincia de Buenos Aires and the Facultad de Ciencias Exactas de la Universidad de Buenos Aires for supporting a short visit to Argentina.