NISTIR 7272

Thermodynamics of Sessile Drops on a Rigid

Substrate: A Comparison of Two Theories

Jeffrey W. Bullard
NISTIR 7272

Thermodynamics of Sessile Drops on a Rigid

Substrate: A Comparison of Two Theories

Jeffrey W. Bullard

Building and Fire Research Laboratory

November 2005

U.S. DEPARTMENT OF COMMERCE

Carlos M. Gutierrez, Secretary

TECHNOLOGY ADMINISTRATION

William Jeffrey, Acting Under Secretary of Commerce for Technology

NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY

William Jeffrey, Director

AdLBtJDit

Abstract

Minimization of free energy is used to calculate the equilibrium shapes of sessile liquid

drops under different wetting conditions and under the influence of a gravitational

body force. The total system free energy is assumed to obey an equation recently

derived by Searcy [1], and the shapes of sessile drops of different sizes are calculated

and compared with experimental measurements. Stationary states in the free energy are

found using standard variational principles, which lead not only to an Euler-Lagrange

differential equation for the drop shape, but also to the boundary conditions at the

three-phase junction where the liquid, solid, and surrounding vapor meet. The analysis

shows that the Searcy equation leads to the classical Young-Dupr ́e equation for the

thermodynamic contact angle. Drop shapes calculated by this approach are in excellent

agreement with experimental observations of drops of pure water resting on a clean

polytetrafluoroethylene (PTFE) substrate.

Keywords: thermodynamics, equilibrium, liquids, wetting

i
Contents

Abstract i

List of Figures iii

List of Tables iv

1 Introduction 1

2 Theory 3

2.1 Derivation of Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2.2 Evaluation of Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.2.1 Evaluation of Boundary Condition at θ = 0 . . . . . . . . . . . . . . . . . . . 9

2.2.2 Evaluation of Boundary Condition at θ = π/2 . . . . . . . . . . . . . . . . . . 9

2.2.3 Evaluation of the Euler-Lagrange Equation . . . . . . . . . . . . . . . . . . . 11

3 Numerical Results 12

4 Discussion and Conclusions 14

A Appendix 16

A.1 First Fundamental Form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

A.2 Surface Normal Vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

A.3 Second Fundamental Form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

A.4 Mean Curvature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

B Appendix 19

ii
List of Figures

1 (a) Schematic drawing of a sessile liquid drop partially wetting a flat solid substrate

with thermodynamic contact angle Θw. ~g indicates the gravitational acceleration

vector. (b) Axisymmetric liquid drop, showing the spherical polar coordinate system

that is used in the derivation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

2 Enlarged view of region near the three-phase junction between solid, liquid, and vapor. 10

3 Cross-sectional profiles of equilibrium shapes of liquid sessile drops with volume of

1.0 mL. Profiles are shown for Θw = 45◦, 90◦, and 135◦. . . . . . . . . . . . . . . . . 13

4 Cross-sectional profiles of equilibrium shapes of liquid sessile drops with Θw = 90◦.

Profiles are shown drop volumes of 0.1 mL, 0.5 mL, 1.0 mL, and 2.0 mL. . . . . . . . 13

5 Photograph of a 2 mL drop of pure water resting on a clean, flat surface of PTFE.

The red dashed line shows the profile calculated using the same physical parameters

as in the experiment. Photograph courtesy of F. Barberis [2]. . . . . . . . . . . . . . 14

6 Pressure difference across the liquid-vapor surface as a function of volume of sessile

liquid drops when γlv = 72.2 mJ/m2, ρ = 1000 kg/m3, and g = 9.8 m/s2. . . . . . . . 15

7 Cylindrical polar coordinate system used in the numerical solution of the EulerLagrange equation. The liquid-vapor surface is parameterized by its arc length s.

ψ(s) is the tangent angle to the drop (ψ ≤ 0). . . . . . . . . . . . . . . . . . . . . . . 20

iii
List of Tables

1 Dimensionless groupings of variables. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

iv
1 Introduction

Wetting and adhesion phenomena are the basis of a diverse range of natural phenomena and

engineering processes. The tendency of a liquid to wet a solid is determined by the free energy

reduction that accompanies the covering of the solid by the liquid, relative to the combined free

energies of the same mass of liquid and solid in isolation.

A large number of theoretical and experimental investigations of wetting have been published

in the past 200 years. Many of these have focused on the equilibrium shape of a liquid drop or

a solid particle resting atop a planar solid surface. A schematic of a sessile liquid drop partially

wetting a flat solid is shown in Fig. 1(a). If the drop is sufficiently small, or if the influence of

gravity is negligible, then the drop is observed to have a shape that approximates a spherical cap,

although the particular shape is not of concern at this point. What is important is that the angle,

Θw, at which the liquid-vapor surface meets the solid-vapor surface appears to be a thermodynamic

property depending only on the three different interfaces that are present (liquid-vapor, solid-vapor,

and liquid-solid). This angle often is called the thermodynamic contact angle for the system.

Figure 1: (a) Schematic drawing of a sessile liquid drop partially wetting a flat solid substrate with

thermodynamic contact angle Θw. ~g indicates the gravitational acceleration vector. (b) Axisymmetric liquid drop, showing the spherical polar coordinate system that is used in the derivation.

The accepted relation between the thermodynamic contact angle and the thermodynamic properties of the surfaces is the Young-Dupr ́e equation:

cos Θw =

γsv − γls

γlv

(1)

where γsv, γls, and γlv are the surface free energy densities of the solid-vapor surface, the liquid-solid

surface, and the liquid-vapor surface, respectively.

At constant temperature and volumes of the three phases in Fig. 1(a), the thermodynamic

potential, G, which governs equilibrium and stability of the system, is [3]

G = Gb +

X

i

γiAi (2)

1
where Gb is the contribution to G of the bulk phases, Ai

, is the area of the i-th type of surface, and

γi ≡



∂G

∂Ai



Aj6=Ai

is the surface energy density of the i-th type of surface. The summation in Eq. (2) represents the

contribution to the thermodynamic potential of the excess free energy stored in the surfaces.

At thermodynamic equilibrium any infinitesimal change in the configuration of the system

cannot decrease G. For changes at constant temperature and constant volume of the liquid, solid,

and vapor, Gb is constant and the condition of equilibrium is

δG = δ (γlvAlv + γsvAsv + γlsAls) ≥ 0 (3)

which is a particular form of the more general Gibbs criterion [3, 4].

Recently, Searcy has proposed an alternative to Eq. (3) for the condition of equilibrium,

δGl = δ [γlvAlv + (γlv + λl) Als] ≥ 0 (4)

where now the potential Glis the free energy of the liquid, and λlis the partial free energy (energy

per mole of liquid) of bonding of the liquid to the solid, per unit area of the liquid-solid surface

that is formed. The parameter λl can be interpreted in the following way. Imagine that a liquid

drop is somehow deformed into exactly the shape it will have when it wets the substrate, but that

it is not yet brought into contact with that substrate. The free energy of the drop in this deformed

state, relative to the bulk liquid, is well-defined and depends only on its surface area and γlv. Now

the deformed drop is brought into contact with the solid substrate. The shape of the drop does not

change, so the only change in free energy is due to the chemical bonding of liquid and solid across

the newly formed surface. The parameter λlis equal to the difference between the free energy of

the drop in the latter case and the free energy of the identically deformed drop without wetting,

per unit area of contact between liquid and solid. The free energy change is assigned solely to the

liquid, while the free energy of the substrate is assumed to remain constant.

Both Eq. (3) and Eq. (4) ignore the influence of gravity on the thermodynamics of the system,

and it is well known that the equilibrium shape of a macroscopic sessile drop is sensibly affected

by gravity. Most notably, the equilibrium shape, which is size-invariant in the absence of gravity,

is dependent on the volume of the drop when gravity is present.

The purpose of this paper is to provide a thermodynamic basis for calculating the equilibrium

shape of an axisymmetric sessile liquid drop under the influence of gravity. That is, the approach

used here relies on a formal minimization of the appropriate thermodynamic potential, using variational calculus. The analysis leads to an Euler-Lagrange differential equation describing the shape

of the drop. In addition, the approach also naturally derives the boundary condition that must

prevail at the three-phase junction. Prior investigations similar to this have usually assumed the

boundary condition according to Eq. (1). Therefore, the approach used here seems more satisfying

from a thermodynamic perspective. Furthermore, this approach provides a means to directly test

the form of the thermodynamic potential proposed by Searcy in Eq. (4). It will be shown that the

implications of Searcy’s equation, at least for the shape and properties of macroscopic sessile drops,

is indistinguishable from those of the accepted thermodynamic potential originally formulated by

Gibbs, a special case of which is given in Eq. (2).

2
2 Theory

2.1 Derivation of Governing Equations

The system being considered is depicted in Fig. 1(a). It consists of a liquid drop having constant

volume V immersed in a vapor phase and resting on a planar solid substrate. The drop is assumed

to have an axis of rotational symmetry represented by the vertical dashed line in Fig. 1(b). A

gravitational body force ~g is assumed to act in the −x3 direction.

For purposes of being explicit, the following assumptions are made about the system before

proceeding with the analysis.

Assumption 1 The solid is rigid and has a planar interface with the vapor and liquid.

Assumption 2 Each phase has fixed composition and can undergo no chemical reactions or segregation.

Assumption 3 The surface energy density of each surface is constant. In particular, there is no

crystalline anisotropy in the surface energy density of any of the interfaces.

Assumption 4 Each phase has constant mass and volume.

This assumption is reasonable for the kinds of wetting phenomena that are being examined here.

It should be mentioned that this assumption prohibits the formation or annihilation of vacancies

in the solid or liquid that can change the volume at constant mass. However, the contribution to

the configurational component of the thermodynamic potential that these process would allow is

usually small and insensitive to changes in shape at constant mass.

Assumption 5 The vapor is a thermal reservoir that fixes the temperature of the liquid and solid

phases.

Assumption 6 The liquid drop is rotationally symmetric about an axis which intersects the centroid

of the drop.

Assumption 7 The gravitational field is constant over the volume of the drop.

The thermodynamic potential is composed of a contribution Gb from the bulk, a net contribution

Gi from the interfaces, and a contribution Gg due to gravitational potential energy:

G = Gb + Gi + Gg (5)

For the system in Fig. 1, Gi can be written according to Gibbs [3],

3
Gi = γlvAlv + γsvAsv + γlsAls (6)

As stated in the Introduction, Searcy [1] proposed the following modified equation for the system

in Fig. 1

Gi = γlvAlv + (γlv + λl) Als (7)

For the remainder of this section, we will use Eq. (7) to represent Gi, and will compare the results

to those that can be obtained when Eq. (2) is used instead.

The contribution to G of gravitational potential energy is given by an integral over the volume

of the drop:

Gg = −

Z

V

ρ~g · ~r(x1, x2, x3)dV (8)

where ~r is the position vector relative to the origin O, ~g is the acceleration due to gravity, and ρ is

the density of the liquid. Assumptions 4 and 7 allow us to simplify Eq. (8) to

Gg = ρg Z

V

xˆ3 · ~r(x1, x2, x3)dV (9)

where now g = |~g| is the magnitude of the acceleration due to gravity. The integral in Eq. (9) will

be recognized as the expression for the elevation of the centroid of the drop above the plane of the

substrate.

Because of Assumption 1, we may conveniently use a spherical polar coordinate system, as

shown in Fig. 1(b). The rotational symmetry of the drop means that its shape is independent of

the azimuthal angle φ. Thus the liquid-vapor interface may parametrized by the polar angle θ, so

that the ordered pair of parametric functions (X1(θ), X3(θ)) specifies the location of liquid-vapor

interface for any value of θ. X1 and X3 both are one-to-one functions of θ unless some portion of the

liquid-vapor interface lies parallel to the position vector R(θ) in Fig. 1(b) for θ ∈ (0, π/2). Later,

in Section 3, a different parameterization will be helpful for numerical solution of the problem.

However, for now the parameterization in θ is convenient for the variational derivation because it

provides fixed minimum and maximum values of the parameter (0 and π/2, respectively).

With the coordinates system defined, for now, as in Fig. 1(b), we can rewrite Eq. (7) as

Gi = π (γlv + λl) X2

1

π

2



+ 2πγlv Z π/2

0

X1(θ)

q

X ̇

1

2

(θ) + X ̇

3

2

(θ) dθ (10)

where X ̇

1 and X ̇3 represent the first derivatives of X1(θ) and X3(θ), respectively. Also, X1 and X3

both can be written in terms of R in Fig. 1(b):

X1 = R sin θ X ̇

1 = R ̇ sin θ + R cos θ (11)

X3 = R cos θ X ̇

3 = R ̇ cos θ − R sin θ (12)

4
Variable Dimensionless Variable Dimensionless

λl α =

λl

γlv

R R =

R

L

G G =

G

πγlvL2 ρg β =

ρgL2

γlv

V V =

V

πL3

Table 1: Dimensionless groupings of variables.

Therefore, Eq. (10) becomes

Gi = π (γlv + λl) R

2

π

2



+ 2πγlv Z π/2

0

R(θ)

q

R ̇ 2(θ) + R2(θ) sin θ dθ (13)

The gravitational term Gg in Eq. (9) can be explicitly written in integral form as

Gg = 2πρg Z π/2

0

sin θ dθ

Z R(θ)

0

(r cos θ) r

2 dr

=

π

2

ρg Z π/2

0

R

4

(θ) cos θ sin θ dθ (14)

The thermodynamic potential in Eq. (5) can now be written in integral form as

G = G ̃



R, R ̇



= π (γlv + λl) R

2

(π/2) + π

Z π/2

0

h

2γlvR

p

R ̇ 2 + R2 +

ρ

2

gR4cos θ

i

sin θ dθ (15)

We seek to minimize F subject to the constraint that the volume V of the liquid is fixed:

L

3 ≡ V =

2π

3

Z π/2

0

R

3

sin θ dθ = constant (16)

where we define L = V

1/3

to be a characteristic dimension of the drop.

The analysis is expedited by defining dimensionless variables as shown in Table 2.1. With these

dimensionless groupings, Eqs. (15) and (16) become

G



R, R ̇



= (1 + α) R2

π

2



+

Z π/2

0



2R

p

R ̇ 2 + R2 +

β

2

R4cos θ



sin θ dθ (17)

V =

2

3

Z π/2

0

R3sin θ dθ (18)

The variational approach used here is adapted from Segel [5]. Assume that R(θ) is the extremizing function, and examine small perturbations of this function having the form

5
R(θ) = R(θ) + 1ζ1(θ) + 2ζ2(θ) (19)

R ̇ (θ) = R ̇ (θ) + 1

 ̇ζ1(θ) + 2 ̇ζ2(θ) (20)

where ζ1 and ζ2 are arbitrary C1 smooth functions of θ and 1, 2  1.

We can substitute the trial solution, Eqs. (19) and (20) into Eqs. (17) and (18) to obtain

G (1, 2) = (1 + α) R2

π

2



+

Z π/2

0



2R

p

R ̇ 2 + R2 +

β

2

R4cos θ



sin θ dθ (21)

V (1, 2) = 2

3

Z π/2

0

R3sin θ dθ (22)

Using the method of Lagrange undetermined multipliers, the minimzation of Eq. (21) subject to

the constraint in Eq. (22) is equivalent to the minimzation of a new function:

Γ (1, 2) = G (1, 2) − η V (1, 2) (23)

where η is an undetermined constant. The function Γ is a function only of the fluctuations 1 and

2. A stationary value of Γ must satisfy

∂Γ

∂1









2

=

∂Γ

∂2









1

= 0 (24)

For convenience in subsequent equation, we adopt the shorthand notation

Γ1 =



2R

p

R ̇ 2 + R2 +

β

2

R4cos θ



sin θ (25)

Γ2 =

2

3

R3sin θ (26)

Then, evaluating the left derivative in Eq. (24),



∂Γ

∂i



1,2=0

= 2 (1 + α) R

π

2

 dR



π

2



di











1,2=0

+

"Z π/2

0



∂Γ1

∂i

− η

∂Γ2

∂i



dθ

#

1,2=0

= 0 i = 1, 2 (27)

The derivatives inside the integral can be expanded by the chain rule,

6
∂Γi

∂j

=

∂Γi

∂R

dR

dj

+

∂Γi

∂R ̇

dR ̇

dj

+

∂Γi

∂θ

dθ

dj

j = 1, 2;i = 1, 2 (28)

Because θ is the independent variable, dθ/dj = 0, so the third term in this equation is zero.

At the endpoints θ = 0 and θ = π/2, the R varies with 1 and 2 according to

dR(0)

d1

= ζ1(0) (29)

dR(0)

d2

= ζ2(0) (30)

dR



π

2



d1

= ζ1

π

2



(31)

dR



π

2



d2

= ζ2

π

2



(32)

In addition, Eq. (19) implies that

R

π

2







1,2=0

= R

π

2



;

∂Γi

∂R









1,2=0

=

∂Γi

∂R

;

and finally

∂Γi

∂R ̇









1,2=0

=

∂Γi

∂R ̇

i = 1, 2 (33)

Substituting Eqs. (29)–(33) into Eqs. (23), (24) and (28) gives

∂Γ

∂i

= 0 = 2 (1 + α) R

π

2



ζi

π

2



+

Z π/2

0



∂Γ1

∂R

− η

∂Γ2

∂R



ζi dθ

+

Z π/2

0



∂Γ1

∂R ̇

− η

∂Γ2

∂R ̇



 ̇ζi dθ (34)

The second integral can be rewritten using integration by parts:

Z π/2

0



∂Γ1

∂R ̇

− η

∂Γ2

∂R ̇



 ̇ζi dθ = ζi

π

2





∂Γ1

∂R ̇

− η

∂Γ2

∂R ̇



θ=π/2

− ζi(0) 

∂Γ1

∂R ̇

− η

∂Γ2

∂R ̇



θ=0

−

Z π/2

0

d

dθ



∂Γ1

∂R ̇

− η

∂Γ2

∂R ̇



ζi dθ (35)

7
Also, from Eq. (26), it is clear that ∂Γ2/∂R ≡ ̇ 0. Therefore, substituting and grouping like terms

in Eq. (34) gives

0 = ζi

π

2





2 (1 + α) R

π

2



+

∂Γ1

∂R ̇



θ=π/2

− ζi(0) 

∂Γ1

∂R ̇



θ=0

+

Z π/2

0



∂Γ1

∂R

−

d

dθ



∂Γ1

∂R ̇



− η

∂Γ2

∂R



ζ1dθ i = 1, 2 (36)

Because ζ1 and ζ2 are arbitrary smooth functions, the only way to guarantee that Eq. (36) vanishes

is if each term separately is zero:

2 (1 + α) R

π

2



+

∂Γ1

∂R ̇









θ=π/2

= 0 (37)

∂Γ1

∂R ̇









θ=0

= 0 (38)

Z π/2

0



∂Γ1

∂R

−

d

dθ



∂Γ1

∂R ̇



− η

∂Γ2

∂R



ζ1 dθ = 0 i = 1, 2 (39)

Eq. (39) was derived under the assumption that ζ1 and ζ2, while arbitrary, are constant. Therefore,

in principle, η could depend on them. However, it may be shown [5] that η must be independent of

ζ1 or ζ2 or both. Therefore, the term in square brackets of the integrand in Eq. (39) must be zero:

∂Γ1

∂R

−

d

dθ



∂Γ1

∂R ̇



− η

∂Γ2

∂R

= 0 (40)

To make the notation more compact, let

f



R, R ̇



≡ Γ1 − η Γ2 (41)

Then Eqs. (37), (38), and (40) take the form

∂f

∂R ̇









θ=0

= 0 (42)

∂f

∂R ̇









θ=π/2

= −2 (1 + α) R

π

2



(43)

∂f

∂R

−

d

dθ



∂f

∂R ̇



= 0 (44)

8
2.2 Evaluation of Governing Equations

Before solving Eq. (44), we can gain insight about the form the solution must have just by further

evaluating the two boundary conditions, Eqs. (42) and (43), and the Euler-Lagrange equation (44),

in terms of R. In particular we will see that the boundary conditions specify the thermodynamic

contact angle at the three-phase junction (θ = π/2), and that the Euler-Lagrange equation provides

a physical meaning to the Lagrange multiplier, η.

2.2.1 Evaluation of Boundary Condition at θ = 0

Substituting for f in Eq. (42) gives

∂f

∂R ̇









θ=0

=

"

2RR ̇

p

R ̇ 2 + R2

#

θ=0

=

2R(0)R ̇ (0)

q

R ̇ 2(0) + R2(0)

= 0 (45)

From Fig. 1(b) we require that R(0) > 0 for finite drop volumes. Therefore, Eq. (45) implies that

R ̇ (0) = 0 Boundary Condition 1 (46)

which means that the drop profile is flat at its apex. This is a natural consequence of the rotational

symmetry that has been imposed in the statement of the problem. Under rotational symmetry,

any nonzero value of R ̇ at θ = 0 would correspond to a singularity in mean curvature which is

unphysical, as we will see later.

2.2.2 Evaluation of Boundary Condition at θ = π/2

Substituting for f in Eq. (43) gives

∂f

∂R ̇









θ=π/2

=

2R



π

2



R ̇



π

2



q

R ̇ 2



π

2



+ R2



π

2

 = −2 (1 + α) R

π

2



(47)

which can be simplified to

R ̇

π

2



= − (1 + α)

r

R ̇ 2

π

2



+ R2

π

2



(48)

Squaring both sides and rearranging leads to a further simplification

R ̇

π

2



= ±

(1 + α) R



π

2



p

− (2α + α2)

(49)

But in Eq. (48), the radical on the right must be non-negative, so we expect that

R ̇

π

2



< 0 ⇐⇒ α > −1

Therefore, we choose the negative root in Eq. (49).

9
Figure 2: Enlarged view of region near the three-phase junction between solid, liquid, and vapor.

R ̇

π

2



= −

(1 + α) R



π

2



p

− (2α + α2)

(50)

This equation places limits on the possible values of α, because R ̇ (π/2) is required to be a real

value. Therefore, a condition on α is

−2 < α ≤ 0 (51)

or, equivalently

−2γlv < λl ≤ 0 (52)

Furthermore, R ̇ (π/2) can be related to the thermodynamic contact angle, Θw, in Fig. 1(a). An

enlarged view of the region near the three-phase junction is shown in Fig. 2. referring to the small

triangle outside the drop in this figure, basic trigonometry requires that

R (θ) cos (∆θ) − R (θ + ∆θ) = p

R ̇ 2 + R2 cos (π − Θw) ∆θ (53)

Dividing both sides of Eq. (53) by ∆θ, taking the limit as ∆θ approaches zero, and evaluating

at θ = π/2,



lim

∆θ→0

p

R ̇ 2 + R2 cos (π − Θw)



θ=π/2

=



lim

∆θ→0

R (θ) − R (θ + ∆θ)

∆θ



θ=π/2

≡ −R ̇

π

2



(54)

Therefore, using the fact that cos(π − x) = − cos(x), we have

R ̇



π

2



q

R ̇ 2



π

2



+ R2



π

2



= cos Θw (55)

10
Comparing to Eq. (48) shows that

cos Θw = − (1 + α) Boundary Condition 2 (56)

Substituting for α from Table 2.1 gives:

γlv cos Θw = γlv + λl (57)

Searcy gives this equation in the same paper in which he introduced the partial free energy of

bonding, λl[1], and it now has been shown how it arises naturally as a boundary condition on the

solution for the free energy-minimizing drop shape.

Remark 1 If we make the assignment

λl ≡ (γls − γsv) − γlv (58)

then Eq. (57) is identical to the Young-Dupr ́e result in Eq. (1).

Remark 2 Eq. (57) is independent of the influence of gravity, even though gravitational contributions to the thermodynamic potential have been included in the problem.

2.2.3 Evaluation of the Euler-Lagrange Equation

Substituting for f in Eq. (44) gives, after a considerable amount of algebraic manipulation,







2R ̇ 2 + R2 − RR ̈



R ̇ 2 + R2

3/2

+

R − R ̇ cot θ

R

p

R ̇ 2 + R2





 + β R cos θ = η (59)

Recall that η, the Lagrange multiplier imposed by the volume constraint on the drop, must be

constant. Therefore the grouping of terms on the left side of Eq. (59) must have fundamental

significance. In fact, R(θ) cos θ is the dimensionless elevation (h/L) above the substrate surface at

a point on the drop surface specified by θ. And because in Table 2.1 β is the dimensionless form

of ρg, the last term on the left side of Eq. (59) is the dimensionless contribution of the hydrostatic

pressure at any height within the drop due to gravity.

Principles of differential geometry [6] can be used to show that the term within square brackets

on the left side of Eq. (59) is the local (dimensionless) mean curvature of the liquid-vapor surface1

,K = LK, at any point2specified by θ (see Appendix A). Therefore, this term represents the

contribution to the pressure within the drop due to the presence of a curved surface as originally

derived by Laplace [7].

In summary, then,the form of the Euler-Lagrange equation is equivalent to

1We use the convention, common in the theory of capillary phenomena, of denoting by the term “mean curvature”

a property of surfaces that is actually two times the mean curvature as it is defined in differential geometry.

2For θ → 0, cot θ in Eq. (59) is singular. Appendix A deals with this in greater detail.

11
η

V 1/3

= −

ρgh(θ)

γlv

+ K(θ) = constant (60)

At any point on the liquid-vapor surface, then, the mean curvature must decrease as the height

increases in such a way as to keep the pressure within the drop constant everywhere. Thus, η is seen

to be the thermodynamic force conjugate to the volume of the liquid, which is constant everywhere

inside the drop at mechanical equilibrium.

Eq. (60) was derived by minimizing Searcy’s expression for the free energy (Eq. (7). However,

because Eq. (7) and Eq. (2) have exactly the same functional dependence on the areas of the various

interfaces, Eq. (60) does not depend on which of these two expressions for the thermodynamic

potential is used.

3 Numerical Results

The remaining task is to solve the Euler-Lagrange equation (59) subject to the boundary conditions

given in Eqs. (45) and (50). Eq. (59) is a second-order nonlinear ordinary differential equation

(ODE), and there are no analytical techniques available to solve it. Therefore, we appeal to a

numerical solution. To numerically solve an nth-order ODE, it is first convenient to re-express it

as n coupled first-order ODEs, and there are standard techniques to accomplish this. Furthermore,

Eq. (59) has a singularity at θ = 0. The singularity can be handled by reparameterizing the

problem in terms of arc length of the liquid-vapor surface and by transforming from spherical polar

coordinates to cylindrical coordinates. When this is done, the singularity is still present but can be

removed by expanding about the singular point. The details of the reparameterization, coordinate

transformation, and removal of the singularity are given in Appendix B.

The numerical approach used here involves two algorithms from the NIST core mathematical

library,3(“CMLIB”): SDRIV3 and SNSQE. These use a shooting method [8] from the apex of the

drop down to the contact line with the substrate. For a given drop volume, the solution is iterated

by adjusting the value of the Lagrange multiplier until the desired drop volume is obtained.

In the results that follow, we relate the equilibrium shapes of sessile drops to physical quantities

instead of the dimensionless parameters used in the previous section. To that end, we assume values

of γlv = 72.2 mJ/m2 and ρ = 1000 kg/m3, which are consistent with the value reported for pure

water in contact with air. We also assume g = 9.8 m/s2.

Figure 3 shows calculated cross-section profiles of a sessile liquid drop with volume of 1 mL. The

profile is shown for three different values of the thermodynamic contact angle Θw. The main effect

of increasing Θw is to pull the drop into a more equiaxed shape having greater mean curvature.

The influence of gravity on the equilibrium shape becomes increasingly pronounced as the

drop volume increases. Fig. 4 shows the calculated equilibrium profiles for different drop volumes

when Θw = 90◦. Larger drops adopt more oblate shapes relative to smaller drops because the

gravitational contribution to the free energy is reduced by having the liquid distributed at lower

elevations above the substrate.

It is interesting to compare the calculated equilibrium profiles of drops to experimental observations of sessile drops. Recently, Barberis and Beruto performed measurements of water drops

resting on a flat surface of polytetrafluoroethylene (PTFE) [2]. Both the water and PTFE were

carefully treated to avoid surface contamination. Measurements made of the water-air surface

3The NIST core mathematical library can be found at the NIST Guide to Available Mathematical Software,

http://gams.nist.gov/.

12
Figure 3: Cross-sectional profiles of equilibrium shapes of liquid sessile drops with volume of 1.0 mL.

Profiles are shown for Θw = 45◦, 90◦, and 135◦.

Figure 4: Cross-sectional profiles of equilibrium shapes of liquid sessile drops with Θw = 90◦.

Profiles are shown drop volumes of 0.1 mL, 0.5 mL, 1.0 mL, and 2.0 mL.

13
Figure 5: Photograph of a 2 mL drop of pure water resting on a clean, flat surface of PTFE. The red

dashed line shows the profile calculated using the same physical parameters as in the experiment.

Photograph courtesy of F. Barberis [2].

energy gave γlv = 72 mJ/m2, with very little change as a function of time. Figure 5 shows a

photograph of a side view of a 2-mL drop from that study. The height of the drop was measured 10

times, giving an average value of 4.61 mm and an uncertainty of 0.2 mm, expressed as one standard

deviation [2]. The macroscopic contact angle was measured with image analysis software and was

reported to be 104◦ ± 0.5

◦

, where the uncertainty is expresssed in terms of the precision of the

measuring device [2].

Superimposed on the photograph in Fig. 5 is a calculated equilibrium profile. For the calculations, the conditions of the experiment were reproduced as closely as possible (γlv = 72 mJ/m2

,

ρ = 1000 kg/m3, and g = 9.8 m/s2). Different values of Θw were calculated, and the best fit to

the experimentally observed profile was obtained for Θw = 110◦. This is a modestly greater value

than that reported experimentally [2]. However, it should be remembered that only the apparent,

macroscopic contact angle was measured. The distortion of the drop due to the influence of gravity

can cause the apparent contact angle to deviate somewhat from the actual thermodynamic contact

angle, the latter of which can reliably be measured only at much higher magnifications. Taking

these considerations into account, the agreement between theory and experiment is remarkably

good.

Eq. (60) indicates that the pressure within a liquid drop should be a sharply decreasing function

of the drop volume. Fig. 6 shows this relation for water (γlv = 72.2 mJ/m2) when Θw = 45◦, Θw =

90◦, and Θw = 135◦. As expected, the pressure difference increases with increasing contact angle

because the drop assumes an equilibrium shape with increasing mean curvature as Θw increases,

as shown in Fig. 3.

4 Discussion and Conclusions

Initially, the analysis presented in this report was undertaken with the hope of comparing Eq. (7),

recently proposed by Searcy [1], to the more commonly accepted form in Eq. (2). It quickly became

evident that the only difference between these two thermodynamic potentials, from a macroscopic

point of view, is the slightly different expression that is derived for the thermodynamic contact

angle given in Eq. (57). In fact, as already described, the simple substitution λl = (γls − γsv) − γlv

makes that equation identical to the equation of Young [9]. Other than this, the consequences of the

Searcy formulation and the Gibbs formulation are mathematically equivalent. That is, there is no

way to distinguish between the two thermodynamic potentials using only macroscopic observations

14
Figure 6: Pressure difference across the liquid-vapor surface as a function of volume of sessile liquid

drops when γlv = 72.2 mJ/m2, ρ = 1000 kg/m3, and g = 9.8 m/s2.

15
of drop shapes. Adamson [10] has pointed out that there are no apparent experimental means to

measure γls and γsv independently. In fact,if γlv is known from a separate experiment, measurement

of the contact angle in a wetting experiment can provide only the difference between γls and γsv.

Searcy’s λl parameter provides another way to look at the wetting phenomenon that basically avoids

the issue of the individual values of γls and γsv. In principle, the thermal free energy contribution

to λl could be independently measured by immersion calorimetry, but there are also likely to

be configurational entropic contributions that likely would make an independent measurement of

λl difficult. Nevertheless, despite the conceptual differences between Searcy’s approach and the

traditional approach to wetting, this report demonstrates that they both make precisely the same

predictions about the equilibrium shapes of sessile drops on rigid solid substrates.

The common approach, employed by Young [9], for relating the contact angles is to view γsv,

γls, and γlv as force vectors that must balance at mechanical equilibrium. The Young-Dupr ́e

equation follows by summing the horizontal components of these three forces to zero. In contrast,

the present analysis shows how Eq. (1) arises as a necessary boundary condition imposed by the

minimization of the system’s thermodynamic potential. The same conclusions have been reached

using variational methods on other systems [11, 12]. In the usual force-balance approach for deriving

the Young equation, the balance of forces in the vertical direction is neglected. It seems to be

implicitly assumed that the solid substrate is capable of balancing the vertical component of the

liquid-vapor surface tension. In the energy-minimization approach used in the present study, the

analogous assumption is that the substrate is rigid, i.e. it stores no elastic energy because it cannot

be deformed. In reality, all solids are deformable, and at equilibrium the substrate surface will

adopt a shape that minimizes a modified thermodynamic potential in which the elastic energy of

the substrate is included along with all the other contributions that have been considered in this

analysis. In the limiting case for which the substrate has zero shear modulus (i.e. a liquid) it has

been shown [12] that the condition at the three-phase junction is that of a von Neumann triangle

when gravity is neglected, regardless of whether the solution is obtained by minimization of the

thermodynamic potential or by a force balance approach. Similar results also have been obtained

for liquid drops on a thin elastic membrane [13].

In the present case, it has been shown that gravity has no influence on the thermodynamic

contact angle, although clearly gravity can influence the macroscopic, or apparent, contact angle

observed at lower magnifications. More generally, no body forces can exert an influence on the

thermodynamic contact angle. Previously, it also was shown that the thermodynamic contact

angle is not influenced by elastic mismatch stresses between a wetting solid particle and its solid

substrate [14, 15]. Therefore, the contact angle is a robust property of a given system, and depends

only on the energy densities of the participating interfaces.

Acknowledgments

Alan Searcy is gratefully acknowledged for many hours of discussion about the testable differences

between the two theories. G.B. McFadden, of the Mathermatical and Computational Sciences

Division at NIST, provided invaluable help with the numerical computations and also provided a

number of suggestions for improving the manuscript.

A Appendix

In this section we derive the form for the mean curvature of the surface shown in Fig. 1. In spherical

polar coordinates, the position vector ~x specifying the surface is parameterized by φ and θ. We

16
will restrict attention to surfaces of revolution like that considered in this study, for which

~x (φ, θ) =





R(θ) sin θ cos φ

R(θ) sin θ sin φ

R(θ) cos θ





(61)

A.1 First Fundamental Form

The first fundamental form for the arc length, s, is then given by [6]

ds

2 = gφφ (dφ)

2 + 2gφθdφdθ + gθθ (dθ)2

(62)

where the components of the tensor g are

gφφ =

∂~x

∂φ ·

∂~x

∂φ

= R2sin2θ (63)

gφθ =

∂~x

∂φ ·

∂~x

∂θ

= 0 (64)

gθθ =

∂~x

∂θ ·

∂~x

∂θ

= R2(θ) + R ̇ 2(θ) (65)

A.2 Surface Normal Vector

The unit vector ˆn(φ, θ) normal to the surface at any point is given by

nˆ ≡

∂~x

∂θ ×

∂~x

∂φ







∂~x

∂θ ×

∂~x

∂φ







(66)

Evaluating the derivatives in Eq. (66) and simplifying gives

nˆ =

1



R ̇ 2 + R2

1/2





R ̇ cos θ cos φ − R sin θ cos φ

R ̇ cos θ sin φ − R sin θ sin φ

R ̇ sin θ + R cos θ





(67)

17
A.3 Second Fundamental Form

The second fundamental form, II, relates to the local deviation of the surface away from its tangent

plane at the point of tangency [6]:

II = bφφ (dφ)

2 + 2bφθdφ dθ + bθθ (dθ)2

(68)

where the coefficients of the tensor b are

bφφ =

∂

2~x

∂φ2

· nˆ

=

R2sin2θ − RR ̇ sin θ cos θ



R ̇ 2 + R2

1/2

(69)

bφθ =

∂

2~x

∂φ∂θ · nˆ

= 0 (70)

bθθ =

∂

2~x

∂φ2

· nˆ

=

R2 + 2R ̇ 2 − RR ̈



R ̇ 2 + R2

1/2

(71)

A.4 Mean Curvature

The mean curvature of any smooth surface is given in terms of the tensors g and b according to [6]

K =

gφφbθθ − 2gφθbφθ + gθθbφφ

gφφgθθ − g

2

φθ

(72)

For a surface of revolution, this can be simplified considerably using the results of the previous two

sections:

K (θ) = R − R ̇ cot θ

R



R ̇ 2 + R2

1/2

+

R2 + 2R ̇ 2 − RR ̈



R ̇ 2 + R2

3/2

(73)

The boundary condition on the surface at θ = 0, given in Eq. (46), is that R ̇ (0) = 0. Therefore,

the term containing cot θ in the previous equation will be well-behaved in the neighborhood about

θ = 0 as long as R ̇ approaches zero faster than tan θ approaches zero. We can examine two limiting

18
cases. First, for small enough drops, experiments show that the surface is shaped like a section of

a sphere with radius rs. In this limiting case,

R ̇ (θ) −→ rs (1 − cos Θw) sec θ

h

1 − (1 − cos Θw)

2

sec2θ

i

tan θ

The tan θ in this expression cancels cot θ, so the first term in Eq. (73) remains finite as θ approaches

zero.

Another limiting case is that of very large drops. Experimentally obtained cross-sectional

profiles on large drops indicate that they are approximately planar near the center of symmetry (see,

for example, Fig. 5. If the drop is planar in a neighborhood near θ = 0, then in this neighborhood,

R(θ) ≈ R(0) sec θ

The the derivative with respect to θ is

R ̇ (θ) ≈ R(0) sec θ tan θ

Again, the tan θ in this expression and cot θ terms cancel, so the first term in Eq. (73) remains

finite when the drop is nearly flat at its apex.

At the point θ = 0, the assumption of rotational symmetry dictates that the two principal

curvatures must be equal. Therefore, for any drop size or shape, the mean curvature at θ = 0 must

be given by

K(0) = 2gθθ

bθθ

=

2R2 + 4R ̇ 2 − 2RR ̈



R ̇ 2 + R2

3/2

(74)

B Appendix

Because cot θ in Eq. (59) is singular at θ = 0, a change both in parameterization of the liquid-vapor

surface and a transformation of the coordinate system is used, as shown in Fig. 7.

The relations between quantities in the old and new coordinate system are:

R sin φ −→ R(t) (75)

R cos φ − R(0) −→ −Z(t) (76)

R

π

2



−→ R(tf ) (77)

R(0) −→ −Z(tf ) (78)

where t = s/l is the dimensionless arc length, tf is the total arc length of the drop profile, and

where the new functions R(s) = r(s)/l and R(s) − z(s)/l are understood to be dimensionless.

It is readily shown [16] that, in the new coordinate system, the principle curvatures are

19
Figure 7: Cylindrical polar coordinate system used in the numerical solution of the Euler-Lagrange

equation. The liquid-vapor surface is parameterized by its arc length s. ψ(s) is the tangent angle

to the drop (ψ ≤ 0).

1

R1

= −

dψ

ds

(79)

1

R2

= −

sin ψ

dR

(80)

(81)

Because the mean curvature is invariant under coordinate transformations, we may substitute this

result into Eq. (59). Thus the Euler-Lagrange equation (59) is transformed to

η ̃ = β Z (t) −

∂ψ

∂t −

sin ψ

R

(82)

In this equation, we use  ̃η = η + β Z (tf ). Furthermore, from Fig. 7, the tangent angle ψ is related

to Z and R according to

dR

dt

= cos ψ (83)

dZ

dt

= sin ψ (84)

(85)

Eqs. (82)–(84) thus are a system of three coupled first-order ordinary differential equations that

must be solved simultaneously for Z(t), R(t), and ψ(t). Furthermore, using Eqs. (75)–(78), the

boundary conditions given in Eq. (46) and (56) are transformed to:

20
ψ (tf ) = −Θw (86)

R(0) = 0 (87)

Z(0) = 0 (88)

(89)

In Eq. (82), the term sin ψ/R has a singularity at t = 0, which can be removed by expanding R,

Z, and ψ about that point:

R(t) = r1t +

r3

3!t

3 +

r5

5!t

5 +

r7

7!t

7 + O(t9

) (90)

Z(t) = z2

2!t

2 +

z4

4!t

4 +

r6

6!t

6 + O(t8

) (91)

ψ(t) = ψ1t +

ψ3

3! t

3 +

ψ5

5! t

5 +

ψ7

7! t

7 + O(t9

) (92)

(93)

which assumes that R(t) and ψ(t) both are odd functions of t, and that Z(t) is even. The coefficients

in these expansions are determined by substituting the expansions into Eqs. (82), (83), and (84).

From Eq. (83), we find

r1 = 1; r3 = −ψ

2

1

;

r5 = −4ψ3ψ1 + ψ

4

1

; r7 = −6ψ5ψ1 − 10ψ

2

3 + 20ψ3ψ

3

1 − ψ

6

1

From Eq. (84),

z2 = ψ1; z4 = ψ3 − ψ

2

1

; z6 = ψ5 − 10ψ3ψ

2

1 + ψ

5

1

From Eq. (82),

ψ1 =

−η ̃

2

; ψ3 =

−3β η ̃

8

;

ψ5 =

β η ̃

3

16

−

5β

2η ̃

16

; ψ7 =

−11β η ̃

5

256

+

33β

2η ̃3

64

−

35β

3η ̃

128

Substituting these coefficients, the singular term sin ψ/R becomes

sin ψ

R

=

−η ̃

2

−

β ηt ̃

2

16

+



β η ̃

3 − β2η ̃

384 

t

4 +



β η ̃

5

61440

+

31β

2η ̃3

46080

−

β

3η ̃

18432

t

6 + O(t8

) (94)

This expression is substituted for sin ψ/R over a small interval 0 ≤ t ≤ t0 in defining the ODE for

numerical solution. If t0 = 0.01, the expansion should give 14 or 15 digits of accuracy.

21
References

[1] A. W. Searcy, The dependence of particle shapes on partial free energies of bonding to inert

substrates, Scripta Mater. 40 (8) (1999) 979–982.

[2] F. Barberis, D. Beruto, Personal communication, university of Genoa, Italy (2003).

[3] J. W. Gibbs, The Collected Works of J. Willard Gibbs, Volume 1, Thermodynamics, Yale

University Press, New Haven, CT, 1928.

[4] C. Herring, Surface tension as a motivation for sintering, in: R. Gomer, C. S. Smith (Eds.),

Structure and Properties of Solid Surfaces, University of Chicago Press, Chicago, IL, 1953, pp.

3–63.

[5] L. A. Segel, Mathematics Applied to Continuum Mechanics, Dover Publications, New York,

1987.

[6] J. J. Stoker, Differential Geometry, Wiley-Interscience, New York, 1969.

[7] P. S. Laplace, Mecanique Celeste, Volume 4, Supplements 1 (1806) and 2 (1807), Little and

Brown, Boston, MA, 1839.

[8] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes in C, 2nd

Edition, Cambridge University Press, London, 1997.

[9] T. Young, in: G. Peacock (Ed.), Miscellaneous Works of the Late Thomas Young, Vol. 1, John

Murray, London, England, 1855.

[10] A. W. Adamson, A. P. Gast, Physical Chemistry of Surfaces, 6th Edition, Wiley-Interscience,

New York, 1997.

[11] R. E. Collins, J. C. E. Cooke, Fundamental basis for the contact angle and capillary pressure,

Trans. Faraday Soc. 55 (1959) 1602–1606.

[12] J. W. Bullard, Fundamental studies of sintering, PhD dissertation, University of California,

Berkeley, CA (1993).

[13] M. E. R. Shanahan, The equilibrium of a sessile drop on a thin solid, J. Chim. Phys. 84 (3)

(1987) 459–464.

[14] J. W. Bullard, M. Menon, Equilibrium shapes of solid particles on elastically mismatched

substrates, J. Colloid Interface Sci. 219 (2) (2000) 320–326.

[15] D. J. Srolovitz, S. H. Davis, Do stresses modify wetting angles?, Acta Mater. 49 (6) (2001)

1005–1007.

[16] H. M. Princen, in: E. Matijevic (Ed.), Surface and Colloid Science, Vol. 2, Wiley-Interscience,

New York, 1969.

22