Search in book...
Toggle Font Controls
Create new playlist

Name your new playlist

Playlist description (optional)
Sign In

Email address

Password

Forgot Password?

or

Continue with Facebook

Continue with Google
Sign Up

Full Name

Email address

Confirm Email Address

Password

or

Continue with Facebook

Continue with Google

Chapter 24

The Shape of Drops

24.1 Introduction

Liquids form drops of distinct shapes when deposited on surfaces. This is a phenomenon we know from everyday life, e.g., when observing the shape of raindrops on a window or on the roof of a car. The contour of the drops formed is defined primarily by the size of the drop, the surface tension of the liquid, and the wetting behavior of the liquid on the surface. The theory for deriving the contour was originally described by Bashforth¹ and Adams² in 1883 [1]. It is interesting to note that the equation describing the contour of droplets was the first practical equation to be tackled by a group of numerical schemes referred to as the Adams-Bashforth methods (see section 27.2.6). In this section we will derive this contour equation for both sessile and pendant droplets.

24.2 Derivation

The drop is assumed to have a circular footprint and to be rotationally symmetric about the z-axis (see Fig. 24.1). The coordinate system in this cut is chosen at the apex of the drop. The position variable s follows the curvature. At any given point, the tangent to the shape of the drop forms the angle θ. The drop is assumed to have a circular contact area and wet the surface with the contact angle Θ₀.

f24-01-9781455731411 — Fig. 24.1 Cut view through a sessile drop with the important geometric relations used for deriving the differential equation describing the drop contour.

Therefore each cut through the drop parallel to the surface on which the drop is deposited will be a circle. The drop is doubly curved with one radius (r_in-plane) being described by the curvature in-plane and one radius r_θ being described by the curvature out-of-plane. The apex of the drop is coinciding with the origin of the coordinate system (x = 0 and z = 0). Here the drop will assume the form of a sphere and the two radii will be identical ${r_{θ} |}_{x = 0, z = 0} = {r_{in-plane} |}_{x = 0, z = 0} = r_{max}$ ${r_{θ} |}_{x = 0, z = 0} = {r_{in-plane} |}_{x = 0, z = 0} = r_{max}$ .

From the Young-Laplace equation (see Eq. 20.11), follows

$\frac{1}{r_{in-plane}} + \frac{1}{r_{θ}} = \frac{Δ p}{γ}$ $\frac{1}{r_{in-plane}} + \frac{1}{r_{θ}} = \frac{Δ p}{γ}$

where the pressure difference is given by the gravitational force, and thus

$\frac{1}{r_{in-plane}} + \frac{1}{r_{θ}} = \frac{Δ ρ g z + p_{0}}{γ}$ $\frac{1}{r_{in-plane}} + \frac{1}{r_{θ}} = \frac{Δ ρ g z + p_{0}}{γ}$

(Eq. 24.1)

where we have introduced the density difference ∆ρ between the density of the drop and the ambiance. This generalization allows using the derived formula also for cases where a drop is suspended in a second liquid; such is the case, e.g., in droplet microfluidics. The constant p₀ is the ambiance pressure at the origin of the coordinate system, i.e., at z = 0. As stated, the two radii become identical at the apex of the drop, in which case Eq. 24.1 becomes

$\frac{2}{r_{\max}} = \frac{p_{0}}{γ}$ $\frac{2}{r_{\max}} = \frac{p_{0}}{γ}$

(Eq. 24.2)

which allows us to remove the constant p₀ by substituting Eq. 24.2 into Eq. 24.1, which yields

$\frac{1}{r_{in-plane}} + \frac{1}{r_{θ}} = \frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}}$ $\frac{1}{r_{in-plane}} + \frac{1}{r_{θ}} = \frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}}$

(Eq. 24.3)

Ultimately, the contour of the drop is to be found as a function z (x). Therefore the two radii have to be expressed as functions of x and z. For r_θ we find

$\begin{array}{l} x = \sin θ r_{θ} \\ r_{θ} = \frac{x}{\sin θ} \end{array}$ $\begin{array}{l} x = \sin θ r_{θ} \\ r_{θ} = \frac{x}{\sin θ} \end{array}$

si5_e (Eq. 24.4)

Using Eq. 24.4 we can rewrite Eq. 24.3 as

$\frac{1}{r_{in-plane}} + \frac{\sin θ}{x} = \frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}}$ $\frac{1}{r_{in-plane}} + \frac{\sin θ}{x} = \frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}}$

(Eq. 24.5)

Eq. 24.5 is complete and accounts for all the physical effects that must be considered. However, it cannot be solved in this form. Depending on this equation, we may also want to express the two variables x and z differently because we only have one equation and two variables. There have been several attempts in the literature to proceed from here on. We will begin by the classical derivation performed by Bashforth and Adams.

24.3 Bashforth and Adams: Curvature Expressed as Z (X)

Bashforth and Adams simply used the inverse of the curvature according to Eq. 3.102 for describing the in-plane curvature, from which we find

$κ_{in-plane} = \frac{\frac{d^{2} z}{d x^{2}}}{{(1 + {(\frac{d z}{d x})}^{2})}^{\frac{3}{2}}} = \frac{1}{r_{in-plane}}$ $κ_{in-plane} = \frac{\frac{d^{2} z}{d x^{2}}}{{(1 + {(\frac{d z}{d x})}^{2})}^{\frac{3}{2}}} = \frac{1}{r_{in-plane}}$

si7_e (Eq. 24.6)

Using Eq. 24.6 we can rewrite Eq. 24.5 as

$\frac{\frac{d^{2} z}{d x^{2}}}{{(1 + {(\frac{d z}{d x})}^{2})}^{\frac{3}{2}}} + \frac{\sin θ}{x} = \frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}}$ $\frac{\frac{d^{2} z}{d x^{2}}}{{(1 + {(\frac{d z}{d x})}^{2})}^{\frac{3}{2}}} + \frac{\sin θ}{x} = \frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}}$

si8_e (Eq. 24.7)

Obviously, we also need a way to express sin θ as a function of z since the latter will be the independent variable. The function connecting z and x is

$\tan θ = \frac{d z}{d x}$ $\tan θ = \frac{d z}{d x}$

(Eq. 24.8)

Here, we use Eq. 3.24, which allows us to express sin θ as a function of tan θ, into which we insert Eq. 24.8. This results in

$\sin θ = \frac{\tan θ}{\sqrt{1 + {(\tan θ)}^{2}}} \sin θ = \frac{\frac{d z}{d x}}{\sqrt{1 + {(\frac{d z}{d x})}^{2}}}$ $\sin θ = \frac{\tan θ}{\sqrt{1 + {(\tan θ)}^{2}}} \sin θ = \frac{\frac{d z}{d x}}{\sqrt{1 + {(\frac{d z}{d x})}^{2}}}$

si10_e (Eq. 24.9)

Using Eq. 24.9 in Eq. 24.7 results in

$\begin{array}{l} \frac{\frac{d^{2} z}{d x^{2}}}{{(1 + {(\frac{d z}{d x})}^{2})}^{\frac{3}{2}}} + \frac{\frac{\frac{d z}{d x}}{\sqrt{1 + {(\frac{d z}{d x})}^{2}}}}{x} = \frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}} \\ \frac{d z}{d x} \frac{1}{x} (1 + {(\frac{d z}{d x})}^{2}) + \frac{d^{2} z}{d x^{2}} = (\frac{Δ ρ g}{γ} z + \frac{2}{r_{\max}}) {(1 + {(\frac{d z}{d x})}^{2})}^{\frac{3}{2}} \end{array}$ $\begin{array}{l} \frac{\frac{d^{2} z}{d x^{2}}}{{(1 + {(\frac{d z}{d x})}^{2})}^{\frac{3}{2}}} + \frac{\frac{\frac{d z}{d x}}{\sqrt{1 + {(\frac{d z}{d x})}^{2}}}}{x} = \frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}} \\ \frac{d z}{d x} \frac{1}{x} (1 + {(\frac{d z}{d x})}^{2}) + \frac{d^{2} z}{d x^{2}} = (\frac{Δ ρ g}{γ} z + \frac{2}{r_{\max}}) {(1 + {(\frac{d z}{d x})}^{2})}^{\frac{3}{2}} \end{array}$

si11_e (Eq. 24.10)

Eq. 24.10 is a nonlinear, second-order ODE that gives the contour of the drop as a function z (x). As it is a second-order differential equation, we need two boundary conditions that are given by

$\begin{array}{l} z (0) & = 0 \\ {\frac{d z}{d x} |}_{x = 0} & = 0 \end{array}$ $\begin{array}{l} z (0) & = 0 \\ {\frac{d z}{d x} |}_{x = 0} & = 0 \end{array}$

si12_e

The first boundary condition is due to the fact that the origin is selected to be at z = 0. As this is a maximum, the first derivative at this point is 0, which gives rise to the second boundary condition. However, it is virtually impossible to find an analytical solution to this equation directly. In their original work, Bashforth and Adams used a Taylor series expansion to approximate the equation numerically. Even today, this equation is difficult to solve, which is why several different approaches have been sought.

24.4 Brien, Ben, and Van den Brule: Curvature Expressed as Function of θ (Sessile Drops)

24.4.1 Derivation

In 1991 Brien, Ben, and van den Brule suggested a solution for expressing Eq. 24.5 as a function of θ[3]. This case is especially suitable for sessile drops, as in this case θ is a monotonously increasing independent variable that increases from θ = 0 at the apex to θ = θ₀ at the contact point.

In this case, the curvature is expressed using Eq. 3.98 as the change of the angle θ, which is a function of the path s. We thereby obtain

$κ_{in-plane} = \frac{1}{r_{in-plane}} = \frac{d θ}{d s}$ $κ_{in-plane} = \frac{1}{r_{in-plane}} = \frac{d θ}{d s}$

in which case Eq. 24.5 becomes

$\begin{array}{l} \frac{d θ}{d s} + \frac{\sin θ}{x} = \frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}} \\ \frac{d θ}{d s} = \frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}} - \frac{\sin θ}{x} \end{array}$ $\begin{array}{l} \frac{d θ}{d s} + \frac{\sin θ}{x} = \frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}} \\ \frac{d θ}{d s} = \frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}} - \frac{\sin θ}{x} \end{array}$

si14_e (Eq. 24.11)

We also need means of expressing x and z as functions of θ, which are given by

$\frac{d x}{d s} = \cos θ$ $\frac{d x}{d s} = \cos θ$

(Eq. 24.12)

$\frac{d z}{d s} = \sin θ$ $\frac{d z}{d s} = \sin θ$

(Eq. 24.13)

Dividing Eq. 24.12 and Eq. 24.13 by Eq. 24.11, respectively, results in

$\frac{d x}{d θ} = \frac{\cos θ}{\frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}} - \frac{\sin θ}{x}}$ $\frac{d x}{d θ} = \frac{\cos θ}{\frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}} - \frac{\sin θ}{x}}$

si17_e (Eq. 24.14)

$\frac{d z}{d θ} = \frac{\sin θ}{\frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}} - \frac{\sin θ}{x}}$ $\frac{d z}{d θ} = \frac{\sin θ}{\frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}} - \frac{\sin θ}{x}}$

si18_e (Eq. 24.15)

As can be seen, instead of one second-order ODE (see Eq. 24.10) we now have two coupled first-order ODEs (see Eq. 24.14 and Eq. 24.15). This system is significantly simpler to solve numerically.

24.4.2 Numerical Solution

24.4.2.1 Equations and Boundary Values

We will perform the numerical solution in Maple. In addition to Eq. 24.14 and Eq. 24.15 we will add a third ODE, which will make our solution more practical for application. This equation is the function of the drop volume, which is given by

$\begin{array}{l} d V = π x^{2} d z \\ \frac{d V}{d s} = π x^{2} \sin θ \end{array}$ $\begin{array}{l} d V = π x^{2} d z \\ \frac{d V}{d s} = π x^{2} \sin θ \end{array}$

si19_e (Eq. 24.16)

where we have used Eq. 24.12 to express dz as a function of ds. We now divide Eq. 24.16 by Eq. 24.11, resulting in

$\frac{d V}{d θ} = \frac{π x^{2} \sin θ}{\frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}} - \frac{\sin θ}{x}}$ $\frac{d V}{d θ} = \frac{π x^{2} \sin θ}{\frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}} - \frac{\sin θ}{x}}$

si20_e (Eq. 24.17)

The system of first-order ordinary differential equations we want to solve is thus given by Eq. 24.14, Eq. 24.15, and Eq. 24.17. The independent variable of this system is θ ranging from 0 to Θ₀. The unknown variables are x (θ), z (θ), V (θ), and r_max. Since we have four first-order ODEs we need four boundary conditions that are given by

$\begin{array}{l} x (θ = 0) = 0 \\ z (θ = 0) = 0 \\ V (θ = 0) = 0 \\ V (θ = Θ_{0}) = V_{0} \end{array}$ $\begin{array}{l} x (θ = 0) = 0 \\ z (θ = 0) = 0 \\ V (θ = 0) = 0 \\ V (θ = Θ_{0}) = V_{0} \end{array}$

si21_e

where both Θ₀ and V₀ are user-supplied values. In the following, we will solve these equations numerically using Maple.

The Maple listing we will be using the solve the contour of the drop is shown in listing 24.1. We will briefly walk through it. In line 1 we restart the Maple server and include Core.txt. In line 4 we define the gravitational constant g and the variable maxpoints, which defines the number of supporting points for the numerical solver. 500 is a good starting point but this value may need to be increased, especially if the drop spreads very far, which is the case, as we will see, for very small values of Θ₀ and very high values for V₀. The following lines define the three ODEs that are given by Eq. 24.14, Eq. 24.15, and Eq. 24.17. As you may have noticed, there is a fourth ODE that simply states that the solver should consider r_max as a constant, i.e., not dependent on the independent variable θ. As you can see, the ODEs have been written as functions of the density (called _rho) and the surface tension (called _gamma).

Listing 24.1

[Maple] Listing for calculating the contour of a sessile drop. A digital version of the listing can be found under the name SessileDrop.mw in the supplementary material.

restart:read "Core.txt":
#define the ODEs and the boundary conditions
g:=9.81:maxpoints:=500;
ode1:=(_rho,_gamma)->diff(x(theta),theta)=cos(theta)/(-sin(theta)/x(theta)+_rho*g/_gamma*z(theta)+2/rmax(theta));
ode2:=(_rho,_gamma)->diff(z(theta),theta)=sin(theta)/(-sin(theta)/x(theta)+_rho*g/_gamma*z(theta)+2/rmax(theta));
ode3:=(_rho,_gamma)->diff(V(theta),theta)=Pi*(x(theta))^2*sin(theta)/(-sin(theta)/x(theta)+_rho*g/_gamma*z(theta)+2/rmax(theta));
ode4:=diff(rmax(theta),theta)=0;
bc1:=x(0)=0;
bc2:=z(0)=0;
bc3:=V(0)=0;
bc4:=(V0,theta0)->V(theta0)=V0;
#this function extracts the plot contour from the numerical result for a nicer plot
extractshape:=proc(data) local show,i;
 show:=Matrix(2*maxpoints,2);
 for i to maxpoints do
   show[i+maxpoints,1]:=data[2,1][i,4]:
   show[i+maxpoints,2]:=data[2,1][maxpoints,5]-data[2,1][i,5]:
   show[maxpoints+1-i,1]:=-data[2,1][i,4]:
   show[maxpoints+1-i,2]:=data[2,1][maxpoints,5]-data[2,1][i,5]:
 end do:
 return show;
end proc:
#this function performs the numerical solution
solution:=(V0,theta0,_rho,_gamma)->dsolve({ode1(_rho,_gamma),ode2(_rho,_gamma),ode3(_rho,_gamma),ode4,bc1,bc2,bc3,bc4(V0,evalf[5](theta0/180*Pi))},{x(theta),z(theta),V(theta),rmax(theta)},numeric,method=bvp[midrich,buildestimation(0,evalf[5](theta0/180*Pi),[theta,x(theta),z(theta),V(theta),rmax(theta)],[theta->20*Lc*sin(theta),theta->(1-cos(theta)),theta->theta/theta0*V0,s->1],false,maxpoints)):
#water
rho:=1:local gamma:=72:Lc:=sqrt(gamma/(rho*g)):
#water on glass
quickplot([extractshape(solution(0.1,10,rho,gamma)),extractshape(solution(1,10,rho,gamma)),extractshape(solution(5,10,rho,gamma)),extractshape(solution(10,10,rho,gamma)),extractshape(solution(50,10,rho,gamma)),extractshape(solution(100,10,rho,gamma))],-10..10,0..5,"x [mm]","z[mm]");
#water
rho:=1:local gamma:=72:Lc:=sqrt(gamma/(rho*g)):
#water on PTFE
quickplot([extractshape(solution(1,115,rho,gamma)),extractshape(solution(5,115,rho,gamma)),extractshape(solution(10,115,rho,gamma)),extractshape(solution(50,115,rho,gamma)),extractshape(solution(100,115,rho,gamma)),extractshape(solution(300,115,rho,gamma)),extractshape(solution(500,115,rho,gamma)),extractshape(solution(1000,115,rho,gamma))],-15..15,0..15,"x [mm]","z [mm]");
#mercury
rho:=13.5:local gamma:=485.5:Lc:=sqrt(gamma/(rho*g)):
#mercury on glass
quickplot([extractshape(solution(1,140,rho,gamma)),extractshape(solution(5,140,rho,gamma)),extractshape(solution(10,140,rho,gamma)),extractshape(solution(50,140,rho,gamma)),extractshape(solution(100,140,rho,gamma)),extractshape(solution(300,140,rho,gamma)),extractshape(solution(500,140,rho,gamma)),extractshape(solution(1000,140,rho,gamma))],-15..15,0..15,"x [mm]","z [mm]");

24.4.2.2 The Function extractshape

In line 15 we define the function extractshape, which takes the output of the numerical solver as input. It simply extracts the contour and turns it into an array with the first column containing the x-values and the second column containing the z-values. Maple delivers the output of the numerical solver as an array where the entry data[2,1] contains the actual output. data[1,1] contains the heading explaining the layout of the data. In this case, it contains the variable names theta, V(theta), rmax(theta), x(theta), z(theta). It is always helpful to look at this entry first in order to know which values to extract. In our case, we need the entries in column 4 (x-values) and 5 (z-values), which the function extractshape simply copies into a new symmetrical array so that we may plot a symmetrical plot contour.

24.4.2.3 The Function solution

The function solution in line 27 performs the actual numerical solving. It uses the Maple command dsolve and accepts arguments for the drop volume V₀, the contact angle Θ₀ (expected in degrees), as well as the density and the surface tension, which have to be entered in g cm⁻³ and mN m⁻¹, respectively. The contact angle is first converted to radians before being passed on to the fourth boundary condition function bc4. We then pass the argument numeric and the desired method to be a boundary value problem solver method=bvp[midrich] to dsolve.

However, we also need to provide a choice for an algorithm. If we simply pass method=bvp, Maple will abort the numerical solver stating that

   Error, (in dsolve/numeric/bvp) system is singular at left endpoint, use midpoint method instead.

By indicating that Maple should use a midpoint method, we can avoid this problem. Midpoint methods are often useful to avoid endpoint singularities. In this case, we indicate a Richardson interpolation midpoint algorithm (midrich).

24.4.2.4 The Function buildestimation

We also have to provide an initial solution to the solver. This is why we use our function buildestimation, which we have introduced in section 2.3.5. We pass the boundaries (again converting the upper limit for θ to radians) and indicate the order in which we will pass functions for calculating the initial solution in a list: theta, x(theta), z(theta), V(theta), rmax(theta). The next argument is a list in form of the functions that Maple will use to calculate the initial solution. Remember, we only need to supply four functions here since the values for the independent variable θ will be automatically calculated by buildestimation and passed as an argument. For these functions, x was assumed to increase with the sine of θ, whereas z was assumed to increase with the cosine of θ. Note that z (θ = 0) = 0 at the apex of the drop, which is why the function reads (1 − cos θ). You may have noticed the scaling factor 20 Lc for the function x (θ). Here, Lc is the capillary length, and the value 20 is simply by experimenting. It turns out that this estimation for x is a rather critical number that strongly influences whether or not the solver will actually converge. This is especially critical for small contact angles Θ₀ and high volumes V₀. The given values provide sufficient stability for most cases, but they may be worth optimizing further if required. For V(theta) we simply assume a linear increase and rmax(theta) is primed with 1. We don’t requested outputting plots for the estimations (false) and we indicate that the estimation is to be made for the number of points given by maxpoints.

24.4.2.5 Contours of Water Drops

We are now ready to solve the drop contours. We start with the case of water on a glass surface for which we assume a contact angle Θ₀ = 10°. This is the code section starting from line 30. Here we have first defined the values for the density, the surface tension, and the capillary length. One thing to note here is that gamma is protected in Maple since it is used as an internal variable. This is why we need to declare our gamma as local. Otherwise, Maple would throw an error. Then we call the function quickplot that was described in section 2.3.2. As arguments we pass the outputs of several numerical solutions using the same contact angle, density, and surface tension but different volumes: 0.1 µl, 1 µl, 5 µl, 10 µl, 50 µl, and 100 µl. These outputs are first passed to the function extractshape, which converts the numerical data to an array that can be plotted. We plot the results in the range of − 10 mm to 10 mm for the x-axis and 0 mm to 5 mm for the z-axis with the according axes labels. The output of this plot is shown in Fig. 24.2a. As can be seen, the drop spreads very quickly, i.e., the values for x become very large, very quickly. However, the values for z seem to converge at a given point. We will elucidate this fact in a moment.

f24-02-9781455731411 — Fig. 24.2 Numerically calculated drop contours of water drops on surfaces. The calculations have been made using the given Maple code. The volumes of the drops are 0.1 µl, 1 µl, 5 µl, 10 µl, 50 µl, and 100 µl for a), and 0 µl, 5 µl, 10 µl, 50 µl, 100 µl, 300 µl, 500 µl, and 1000 µl for b).

In the next step, we will calculate the drop contour for a water drop on a PTFE surface. This is the code section starting from line 35. Here we assume a contact angle of Θ₀ = 115°. We first define the properties of water and proceed likewise. Here, we use the following drop volumes: 0 µl, 5 µl, 10 µl, 50 µl, 100 µl, 300 µl, 500 µl, and 1000 µl. The output of the calculation can be seen in Fig. 24.2b. As expected, the drop stays significantly more compact due to the high contact angle and the high surface tension of water. Again, the values for z seem to converge.

24.4.2.6 Contours of Mercury Drops

This becomes even more obvious in the last case we want to consider; a mercury drop on a glass surface. Here we assume an even bigger contact angle of Θ₀ = 140°. This is the code section starting from line 40. The density of mercury is obviously higher than that of water, as is the surface tension. We used the following volumes for calculating the drop contour: 1 µl, 5 µl, 10 µl, 50 µl, 100 µl, 300 µl, 500 µl, and 1000 µl. As can be seen, the drop starts to flatten out at higher volumes (see Fig. 24.3). The convergence of the z-values we previously observed are even more strongly expressed here.

f24-03-9781455731411 — Fig. 24.3 Numerically calculated drop contour of mercury drop on a glass surface. The calculations have been made using the given Maple code. The volumes of the drops are 1 µl, 5 µl, 10 µl, 50 µl, 100 µl, 300 µl, 500 µl, and 1000 µl.

24.4.2.7 Convergence Values for z

As already observed, for large volumes V₀ the values for z (θ = 0) converge to the value h_max, which is independent of the drop volume. A rough geometric approximation of this value can be made looking at Fig. 24.4. Here, we assume a drop of very large radius for which case the drop becomes approximately a circle with radius r_max. We know that at the apex the two curvatures are identical, in which case the Young-Laplace equation (see Eq. 20.11) becomes

f24-04-9781455731411 — Fig. 24.4 Geometric derivation of the height convergence h_max for sessile drops. In this case, the drop is assumed to have a spherical shape and a very large diameter footprint.

$\begin{array}{l} \frac{2}{r_{\max}} = \frac{Δ ρ g h_{\max}}{γ} \\ h_{\max} = \frac{2}{r_{\max}} L_{c}^{2} \end{array}$ $\begin{array}{l} \frac{2}{r_{\max}} = \frac{Δ ρ g h_{\max}}{γ} \\ h_{\max} = \frac{2}{r_{\max}} L_{c}^{2} \end{array}$

si22_e (Eq. 24.18)

On the other hand, h_max can be derived geometrically from Fig. 24.4 to be

$\begin{array}{l} h_{\max} = (1 - \cos (Θ_{0})) r_{\max} \\ r_{\max} = \frac{h_{\max}}{1 - \cos (Θ_{0})} \end{array}$ $\begin{array}{l} h_{\max} = (1 - \cos (Θ_{0})) r_{\max} \\ r_{\max} = \frac{h_{\max}}{1 - \cos (Θ_{0})} \end{array}$

si23_e (Eq. 24.19)

Substituting Eq. 24.19 in Eq. 24.18 yields

$\begin{array}{l} h_{\max} = 2 \frac{L_{c}^{2}}{h_{\max}} (1 - \cos (Θ_{0})) \\ h_{\max} = L_{c} \sqrt{2 (1 - \cos (Θ_{0}))} \end{array}$ $\begin{array}{l} h_{\max} = 2 \frac{L_{c}^{2}}{h_{\max}} (1 - \cos (Θ_{0})) \\ h_{\max} = L_{c} \sqrt{2 (1 - \cos (Θ_{0}))} \end{array}$

si24_e (Eq. 24.20)

For the three cases discussed, these values would amount to 0.3 mm for the water drop on glass, 3.2 mm for the water drop on PTFE, and 3.6 mm for the mercury drop on glass. The correct values found from the numerical solutions are 0.49 mm for the water drop on glass, 4.74 mm for the water drop on PTFE, and 3.8 mm for the mercury drop on glass. As can be seen, the estimations fit pretty nicely for the first and the third case, whereas in the second case we have not reached the convergence value. This rough estimation shows why water on glass will form very shallow and very widely spreading pools whereas it will form compact drops on a nonwetting surface. The same can be seen for liquids with high surface tension, e.g., mercury, which will also form compact drops. Whether or not a liquid will form puddles or compact drops solely depends on its wetting contact angle Θ₀.

24.4.2.8 Working With the Maple File

As indicated, you can find the Maple file containing the listing in the supplementary material. Feel free to experiment with it and change the values at will. When changing the input values you will most likely run across the following Maple error eventually:

   Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging.

This error indicates that the numerical solver fails to find a converging solution. This problem can sometimes be resolved by changing the number of support points, i.e., by either increasing or decreasing the value of maxpoints. If the solve still fails to find a solution, it may be worth adapting the estimation functions. During the preparation of the file, the x (θ) and z (θ) functions were found to be most critical. For small contact angles Θ₀ try increasing the scaling factor. Often, it is possible to track the solver’s problems by starting with a scenario for which a solution was found and adapt the estimation functions such that they closely match the found solution. Once the estimation functions are adapted, try changing the input values slowly and observe in which direction the solution deviates from the given initial solutions. In these cases, it is advisable to turn on the plot output of the estimation function by passing true in the function buildestimation. See also section 2.3.5 for details on this function.

24.5 Del Río and Neumann: Curvature Expressed as Function of S (Pendant Drop)

24.5.1 Derivation

The derivation described so far is suitable for the derivation of sessile drops. However, it fails to give correct results for the case of drops suspended, e.g., at the tips of the tap, a pipette, or a leaf. The teardrop shape of pendant drops is one of the most natural feeling geometries, and has been used by mankind for centuries as template for the design of jewellery or ceramics. The reason why the algorithms outlined so far fail in the case of the pendant drop is due to the fact that the angle θ is not unique in this case. Drops with high volumes tend to show constrictions close to the suspension point, which results in a (local) decrease of the angle θ (see Fig. 24.5). Therefore a different independent variable should be chosen.

f24-05-9781455731411 — Fig. 24.5 Difficulties using θ as independent variable for pendant drop analysis. θ increases from 0° to θ = 90° at the drop’s maximum diameter. From there on, θ will increase until the infliction point, from where on it will decrease until it reaches again θ = 90° . After this point, it will decrease again.

A very simple and elegant algorithmic implementation to solve this problem was suggested by del Río and Neumann in 1997 [4]. The authors opted for using the path variable s, which follows the contour (see Fig. 24.1). Referring to Eq. 24.5, the in-plane radius rin-plane can be expressed according to Eq. 3.98 as

$κ_{in-plane} = \frac{1}{r_{in-plane}} = \frac{d θ}{d s}$ $κ_{in-plane} = \frac{1}{r_{in-plane}} = \frac{d θ}{d s}$

(Eq. 24.21)

Inserting Eq. 24.21 in Eq. 24.5 yields

$\begin{array}{l} \frac{d θ}{d s} + \frac{\sin θ}{x} & = \frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}} \\ \frac{d θ}{d s} & = - \frac{\sin θ}{x} + \frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}} \end{array}$ $\begin{array}{l} \frac{d θ}{d s} + \frac{\sin θ}{x} & = \frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}} \\ \frac{d θ}{d s} & = - \frac{\sin θ}{x} + \frac{Δ ρ g z}{γ} + \frac{2}{r_{\max}} \end{array}$

si26_e (Eq. 24.22)

Furthermore we know that

$d x = \cos θ d s \to \frac{d x}{d s} = \cos θ$ $d x = \cos θ d s \to \frac{d x}{d s} = \cos θ$

(Eq. 24.23)

$d z = \sin θ d s \to \frac{d z}{d s} = \sin θ$ $d z = \sin θ d s \to \frac{d z}{d s} = \sin θ$

(Eq. 24.24)

which can be directly derived from (see Fig. 24.1). Eq. 24.21, Eq. 24.22, and Eq. 24.23 are already a set of three first-order ODEs that will give us x (s), z (s), and θ (s) as functions of the independent variable s. We will add Eq. 24.16 as a fourth ODE, which allows us to calculate the volume of the drop. This completes the set of ODEs used to determine the shape of the pendant drop.

Originally, del Río and Neumann used the COLSYS library originally developed by Ascher et al.[5]. However, it turns out it is actually not very difficult to solve this set of ODEs. It can be conveniently done using Maple and even a custom-written solver, as we will see in section 27.4.3.4.

24.5.2 Numerical Solution

24.5.2.1 Equations and Boundary Values

We will use the four first-order ODEs, Eq. 24.21, Eq. 24.22, Eq. 24.23, and Eq. 24.16. We will add a fifth ODE that simply indicates that r_max is to be considered a constant:

$\frac{d r_{\max}}{d s} = 0$ $\frac{d r_{\max}}{d s} = 0$

(Eq. 24.25)

The independent variable of our ODE system is the path variable s ranging from 0 to the maximum value s_max, which we do not know. In order to make our problem more applicable for practical applications we will chose the boundary conditions such that we can define the width of the suspension point of the drop. In practical applications a drop would most likely be suspended, e.g., on a pipette tip or a needle for which we know the diameter. Thus when reaching s_max we want the x-value to reach a user-definable value r₀ such that x (s = s_max) = r₀.

As we have four first-order ODEs we need a total of five boundary conditions that are given by

$\begin{array}{l} θ (s = 0) = 0 \\ x (s = 0) = 0 \\ z (s = 0) = 0 \\ V (s = 0) = 0 \\ x (s = s_{\max}) = r_{0} \end{array}$ $\begin{array}{l} θ (s = 0) = 0 \\ x (s = 0) = 0 \\ z (s = 0) = 0 \\ V (s = 0) = 0 \\ x (s = s_{\max}) = r_{0} \end{array}$

si30_e

where both s_max and r₀ are user-supplied values.

In a practical experiment we would suspend a drop to a capillary through which we can add liquid to the drop or remove liquid from it. By doing so, we increase the circumference of the drop contour. Theoretically, we could rewrite the ODE system as functions of the volume and thereby indicate the total volume of liquid we want the drop to have. However, it will be challenging to find good estimation functions that take the volume as independent variable. It is significantly simpler to use s as the independent variable and output the total drop volume once the contour is found. Therefore, in order to emulate our experimental setup we would supply increasing values of s_max to the ODE system and observe how the volume of the calculated drop contour increases. You may have noticed that we do not account for the contact angle in the given equations. This is due to the fact that mostly tubes and capillaries with semi-circular capillary walls are used for pendant drop experiments (see Fig. 24.6). These can accommodate for any contact angle, which is why we do not need to consider it in the equations. In cases where the contact angle is desired, we may exchange the last boundary condition for

f24-06-9781455731411 — Fig. 24.6 A semi-circular capillary wall can accommodate for any contact angle: around 0° (a), close to 180° (c), as well as values in between (b).

$θ (s - s_{\max}) = Θ_{0}$ $θ (s - s_{\max}) = Θ_{0}$

which accounts for the contact angle leaving the value x (s = s_max) variable. In the following, we will solve these equations numerically using Maple.

The Maple listing we will be using to solve the contour of the drop is shown in listing 24.2. We already know most of the function. In line 1 we restart the Maple server and include the file Core.txt. In line 4 we define the gravitational constant g which is negative in this case. Therefore, gravity will pull the drop rather than compressing it, as is the case for sessile drops. Furthermore, we define the variable maxpoints, which defines the number of supporting points for the numerical solver. The value of 200 was found to be sufficient for most calculations on pendant drops. The lines that follow define the ODEs and the boundary values keeping the radius of the suspension point r_max and the upper limit of the independent variable s_max. The ODEs have again been written as functions of the density (called _rho) and the surface tension (called _gamma).

Listing 24.2

[Maple] Listing for calculating the contour of a pendant drop. A digital version of the listing can be found under the name PendantDrop.mw in the supplementary material.

restart:read "Core.txt":
#define the ODEs and the boundary conditions
g:=-9.81:maxpoints:=200;
ode1:=(_rho,_gamma)->diff(phi(s),s)=-sin(phi(s))/x(s)+_rho*g/_gamma*z(s)+2/rmax(s);
ode2:=diff(x(s),s)=cos(phi(s));
ode3:=diff(z(s),s)=sin(phi(s));
ode4:=diff(V(s),s)=Pi*x(s)^2*sin(phi(s));
ode5:=diff(rmax(s),s)=0;
#these are the boundary conditions
bc1:=phi(0)=0;
bc2:=x(0)=0;
bc3:=z(0)=0;
bc4:=V(0)=0;
bc5:=(smax,r0)->x(smax)=r0;
#this function extracts the plot contour from the numerical result for a nicer plot extractshape:=proc(data) local show,i;
 show:=Matrix(2*maxpoints,2);
 for i to maxpoints do
   show[i+maxpoints,1]:=data[2,1][i,5]:
   show[i+maxpoints,2]:=data[2,1][i,6]-data[2,1][maxpoints,6]:
   show[maxpoints+1-i,1]:=-data[2,1][i,5]:
   show[maxpoints+1-i,2]:=data[2,1][i,6]-data[2,1][maxpoints,6]:
 end do:
 #print some more infos
 printf("V=%4.2f ul
", data[2,1][maxpoints,2]);
 return show;
end proc:
#this function performs the numerical solution
solution:=(smax,r0,_rho,_gamma)->dsolve({ode1(_rho,_gamma),ode2,ode3,ode4,ode5,bc1,bc2,bc3,bc4,bc5(smax,r0)},{phi(s),x(s),z(s),V(s),rmax(s)},numeric,method=bvp[midrich],buildestimation(0,smax,[s,phi(s),x(s),z(s),V(s),rmax(s)],[s->s/smax*Pi,s->s,s->s,s->s,s->1],false,maxpoints)):
#water
rho:=1:local gamma:=72:
#from a 1 mm diameter tube
quickplot([extractshape(solution(0.6,0.5,rho,gamma)),extractshape(solution(1,0.5,rho,gamma)),extractshape(solution(2,0.5,rho,gamma)),extractshape(solution(3,0.5,rho,gamma)),extractshape(solution(4,0.5,rho,gamma)),extractshape(solution(5,0.5,rho,gamma)),extractshape(solution(5.4,0.5,rho,gamma))],-4..4,-5..0,"x [mm]","z [mm]");
#from a 0.5 mm diameter tube
ickplot([extractshape(solution(0.6,0.25,rho,gamma)),extractshape(solution(1,0.25,rho,gamma)),extractshape(solution(2,0.25,rho,gamma)),extractshape(solution(3,0.25,rho,gamma)),extractshape(solution(4,0.25,rho,gamma)),extractshape(solution(4.3,0.25,rho,gamma))],-4..4,-5..0,"x [mm]","z [mm]");
#mercury
rho:=13.5:local gamma:=485.5:
#from a 1 mm diameter tube
ickplot([extractshape(solution(0.6,0.5,rho,gamma)),extractshape(solution(1,0.5,rho,gamma)),extractshape(solution(2,0.5,rho,gamma)),extractshape(solution(3,0.5,rho,gamma)),extractshape(solution(4,0.5,rho,gamma)),extractshape(solution(4.2,0.5,rho,gamma))],-4..4,-5..0,"x [mm]","z[mm]");
#from a 0.5 mm diameter tube
quickplot([extractshape(solution(0.6,0.25,rho,gamma)),extractshape(solution(1,0.25,rho,gamma)),extractshape(solution(2,0.25,rho,gamma)),extractshape(solution(3,0.25,rho,gamma)),extractshape(solution(3.4,0.25,rho,gamma))],-4..4,-5..0,"x [mm]","z [mm]");

24.5.2.2 The Function extractshape

We already know the function extraxtshape in line 18 which simply extracts the drop contour from the numerical solutions. Again, the output data are provided by Maple in array form with the entry data[1,1] containing the headings of the output table and the entry data[2,1] containing the actual data. In this case, column 4 contains the x-values whereas column 5 contains the z-values. We have made a small addition to extractshape, which is shown in line 27. This line outputs the volume of the calculated drop contour at the endpoint of the independent variable V (s_max) which is simply the last entry of the second column of the output data. The second column contains the data of V (s) (check the heading entry in data[1,1] to verify).

24.5.2.3 The Function solution

The function solution in line 32 is similar to the implementation used in listing 24.1. Again, we use a midpoint algorithm due to the discontinuity of the singularity at the lower bound. We indicate that the following dependent variables are sought: theta(s), x(s), z(s), v(s), rmax(s). Again, we provide estimations of a first solution for these variables. We assume θ to be linearly increasing with increasing values of s. These initial values are required for the algorithm to converge. If the values are not provided in a way that they increase the solver will typically not converge. Interestingly, the estimations for the other variables can be as simple as that. By indicating to the solver that we assume these values to increase with increasing s, we will obtain stable solutions. Also it is sufficient to prime r_max with the value of 1. This is somewhat surprising as the original implementation of del Río and Neumann emphasized that a good estimation for r_max was key to obtaining stable solutions. In addition the authors opted for priming the values of x and z with elliptical fits, which was not found necessary for the given Maple code.

24.5.2.4 Contours of Water Drops

In line 37 we consider a drop of water suspended on a tube of capillary of 1 mm diameter. Thus we use a value of 0.5 mm for r_max. We proceed by slowly increasing the value of s_max using the values of 0.6 mm, 1 mm, 2 mm, 3 mm, 4 mm, 5 mm, and 5.4 mm. When increasing to values above 5.4 mm, the solver fails to achieve a solution and Maple will throw the following error:

Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging.

It may be possible to further optimize the estimation functions or change the total number of points, i.e., the variable maxpoints, but in most cases this will allow only a slight increase of this value. The output of line 37 can be seen in Fig. 24.7a. For values of s_max> 5.4 mm the solver will fail to find a stable solution as there is no stable drop configuration that can accommodate this s_max value. In practical applications, drops above this value will detach. The function extractshape also outputs the volumes of the drops corresponding to the indicated values of s_max. In this case these values are 0.12 µl, 0.45 µl, 2.10 µl, 5.54 µl, 10.95 µl, 17.20 µl, and 18.24 µl. Therefore we can see that drops with volumes larger than 18.24 µl cannot be suspended on capillaries of 1 mm diameter.

f24-07-9781455731411 — Fig. 24.7 Contours calculated for pendant water drops suspended on capillaries of different diameters. For higher values of s_max the solver fails to find stable solutions, indicating that such drops would detach and not form a stable suspended drop.

In line 40 we repeat this calculation for drops suspended to a capillary of 0.5 mm diameter. Here we use the following values for s_max: 0.6 mm, 1 mm, 2 mm, 3 mm, 4 mm, and 4.3 mm, resulting in drops of volumes 0.08 µl, 0.27 µl, 1.51 µl, 4.34 µl, 8.70 µl, and 9.81 µl. This is again the limit for which the solver finds a stable solution. The calculated contours are shown in Fig. 24.7b.

24.5.2.5 Contours of Mercury Drops

In line 45 we change the density and the surface tension values to the values of mercury and repeat the calculation for the capillary of 1 mm diameter (see Fig. 24.8a) and 1 mm diameter (see Fig. 24.8b). Despite the higher surface tension, mercury can only be suspended in significantly smaller drop volumes to capillaries due to its higher density.

f24-08-9781455731411 — Fig. 24.8 Contours calculated for pendant mercury drops suspended on capillaries of different diameters. For higher values of s_max the solver fails to find stable solutions, indicating that such drops would detach and not form a stable suspended drop.

24.5.2.6 Determining Stable Drops

Please also note one interesting aspect when looking at the output produced from line 45. The input values for s_max are 0.6 mm, 1 mm, 2 mm, 3 mm, 4 mm, and 4.3 mm, resulting in drops of volumes 0.12 µl, 0.45 µl, 2.06 µl, 5.16 µl, 8.70 µl, and 8.39 µl. Notice that the last volume value is actually smaller than the second-to-last value. By increasing the upper limit of the path variable we have actually resulted in a drop with smaller volume. Returning to our experimental setup, this would mean that we add liquid to the drop, thereby obtaining a drop of smaller contour length. This is physically not possible. Therefore, it is likely to assume that such a drop will actually not reach this contour even though it is numerically stable. Another example of this phenomenon is displayed in Fig. 24.9. Here a mercury drop was suspended to a capillary of 10 mm diameter. The values for s_max used were 5.1 mm, 6 mm, 7 mm, 8 mm, 9 mm, 10 mm, and 10.4 mm. The resulting volumes found were 30.5 µl, 95.38 µl, 119.39 µl, 110.04 µl, 76.86 µl, 47.77 µl, and 81.58 µl.

f24-09-9781455731411 — Fig. 24.9 Discontinuity observed on a pendant mercury drop suspended on a capillary of diameter, 10 mm. As the contour length s_max grows the volume first increases and then decreases again.

Again, one can see that the volume first increases and then decreases again. If you look closely you can see that the last two profiles even range above the z = 0 line toward the suspension point. It is likely to assume that albeit a numerical solution can be found for these cases, these drops will not occur physically. Therefore it seems reasonable to only slowly increase the values for s_max and check continuously for decreasing total volume values, which indicate the end of the physically stable regime.

24.5.2.7 Working With the Maple File

As stated, you can find the Maple file containing the listing in the supplementary material. Keeping in mind the stated constraints with respect to drop stability, this sheet will form a solid basis on which you can experiment with drop contours.

24.6 Comparison With Experimental Data

Until now, we have merely studied the results of the drop shapes derived numerically. In Fig. 24.10 you can see the comparison of the numerically predicted and the actual physical shapes of a sessile and a pendant drop of water with a volume of 5 µl each. The sessile drop was deposited on a poly(methyl methacrylate) substrate. The contact angle of water on poly(methyl methacrylate) is around 72° and was measured with the image acquisition software of the instrument used. The pendant drop was suspended on a capillary with a 0.31 mm outer diameter. As you can see in both cases, our numerical predictions fit the actual drop size really well. It is therefore safe to assume that the numerical solutions obtained from listing 24.1 and listing 24.2 will yield physically correct shapes.

f24-10-9781455731411 — Fig. 24.10 Comparison of the physical drop contour with the numerically derived drop contours (red). a) A 5 µl sessile drop of water deposited on a poly(methyl methacrylate) substrate. The water contact angle is around 72° and was determined using the image acquisition software SCA20 on a DSA30 Drop Shape Analyzer by Krüss (www.kruss.de). Overlaid on the physical drop contour is the contour predicted using listing 24.1. As you can see, the drop contour is correctly predicted. b) A 5 µl pendant drop suspended to a capillary with an outer diameter of 0.31 mm. The physical drop contour is overlaid by the drop contour predicted using listing 24.2. As you can see, the contour is also correctly predicted as well.

24.7 Drops Inside of a Fluid

Until now we have only considered the case of the sessile drop on the surface surrounded by ambiance (usually air). However, Eq. 24.5 can also be used to determine the shape of drop contours with a drop inside of a second liquid or a gas bubble adhering to a surface. In all of these cases, we only need to employ the difference in densities and the interfacial surface tensions as input values for ρ and γ. These values can be negative (as is the case of a gas bubble immersed in water). In the case of the sessile drop, the contact angle Θ₀ may also need to be changed. Apart from these modifications the given listings can be used directly to solve these cases. It may be required to adapt the estimation functions to account for sign changes, etc.

24.8 Historical Development of Drop-Shape Analysis

The theory of drop shape analysis was first described by Bashforth and Adams is still the underlying model to this day. Obviously in the late 19^th century developing solutions to nonlinear second-order ODEs was a tedious work that was only possible by usage of reference tables. Nowadays, instruments for surface tension measurements use algorithms based on a complete fit of contour points to solve Eq. 24.5. However, several simplified algorithms and techniques have been developed over the years that allowed measurement of surface tension at significantly reduced computational cost.

In 1976 Hartland and Hartley proposed algorithms based on formula translating system (FORTRAN) code for solving the equation for which they provided data [6]. The first account of this technique being used for the determination of surface tension (rather than simply predicting the shape of a droplet) dates back to the work of Andreas et al. in 1938 [7]. The authors introduced a slightly simplified version of the method of Bashforth and Adams using only two diameters acquired by means of a camera. This method is often referred to as d_e/d_smethod and was explained in section 22.2.4.2. It is still one of the most widely used techniques for measuring surface tension. Numerous methods have been described over the years that use slightly different diameters as suggested for the d_e/d_s method, e.g., the minimum and maximum diameter of the droplet as suggested by Winkel [8]. In all cases, the underlying fluid mechanics is still the same.

The first to use a digital image recording system were Girault et al. in 1984 [9]. In this work Eq. 24.5 was fitted polynomially to the measured contour. Numerous contributions to the literature have been made over the years with the aim of optimizing the speed and the precision of the fitting algorithms by using different sampling points of the measured droplet contour.

Maybe the most important improvement in algorithm developed was provided by Rothenberg et al., who developed a method referred to as axisymmetric drop shape analysis (ADSA) [10]. This method fits the measured contour to Eq. 24.5 using a nonlinear fitting algorithm using the sum of the squares between the measured and the calculated curve. Most modern instruments use modifications of the ADSA algorithms.

24.9 Summary

In this section we have derived the equation of the drop contour for sessile and pendant drops. We encounter droplets of liquids in almost all aspects of microfluidics. As stated, this equation (given by Eq. 24.7) has historical significance not only for the field of fluid mechanics, but also for the field of numerics as it spurred the development of a very important class of numerical methods for solving differential equations. As we have seen, the shape of both sessile and pendant drops can be derived conveniently with algebra tools, e.g., Maple. The contours found can be matched very closely to experimental data obtained from sessile and pendent water droplets. The derived Maple worksheets form a very good basis for further experimentation on drop contours. In section 27.4.3.4 we will use a custom-written solver for solving this equation using a slightly different algorithm.

..................Content has been hidden....................

You can't read the all page of ebook, please click here login for view all page.

Table of Contents for Chapter 24: The Shape of Drops

Create new playlist

Sign In

Sign Up