Phase-field model for quantitative analysis of droplet wetting

Abstract

A general phase-field formalism is presented, aimed at capturing the excess surface energies and geometric profiles of three-dimensional, asymmetric liquid droplets on solid surfaces. To ensure strict volume conservation of the droplet, a nonlinear definition is employed for the internal volume. To facilitate modeling in systems where the interface widths approach zero, an extrapolation method is proposed to interpret data obtained at finite interface widths. The numerically tractable algorithm yields accurate predictions for the excess surface energies and state diagrams of known wetting configurations on both the external and internal surfaces of cylinders, demonstrating quantitative agreement with established results from other models and computational techniques.

---

I. Introduction

Take the classical problem of a liquid droplet partially wetting a solid surface.¹ The calculation of the surface free energy is directly related to the curved liquid–gas and liquid–solid surface areas in three dimensions. This becomes a nontrivial computational task when the liquid droplet takes on an arbitrary, nonaxisymmetric shape. Theoretically, the equilibrium droplet shape is associated with the surface free energy minimum with respect to shape variations. In this work, it is demonstrated that this minimization task can be quantitatively addressed by proposing a phase-field model and implementing a volume-conserved algorithm.

Assume that the three-dimensional (3D) shape of the droplet can be represented by a single scalar droplet-shape function, a phase field, denoted by ϕ(r), where r represents the Cartesian coordinates of the spatial domain. By supplying an additional, fixed phase field ϕ_S(r) to describe the solid surface structure, the excess surface free energy ΔE can be expressed as a functional of ϕ(r). This functional is then minimized to determine the equilibrium configuration of the system. In this framework, all interfaces are treated as diffuse, where the value of ϕ(r) varies smoothly from +1 to 0 across the interface. Most phase-field theories are based on variations of the classical Cahn–Hilliard model.^2–4

This sound idea has been utilized in different versions of phase-field models for wetting.^5–10 Indeed, the adaptability and flexibility of these models have been demonstrated across a wide range of physical scenarios, encompassing both static and dynamic problems. The purpose of the current study is three-fold.

Firstly, a new model is introduced that possesses fundamentally correct properties in all limiting cases. This approach is based on a multicomponent phase-field theory originally proposed for the coexistence of three fluids, which exhibits all the desired characteristics for modeling coexistence of gas (G), liquid (L), and solid (S) phases.¹¹ When applied to wetting problems where the substrate ϕ_S(r) is fixed, the model simplifies to a functional representation of the excess surface energy ΔE as a functional of ϕ(r) alone. This formalism avoids complications present in some earlier models, such as the need to treat the liquid–gas and solid–liquid interfaces separately. A detailed description of the employed model in this study is provided in Section II A and its derivation is documented in Appendices A, B and C.

Secondly, the volume-constraint problem of the liquid droplet is addressed. One key advantage of the phase-field representation is the elimination of the need to explicitly parameterize the curved, two-dimensional (2D) surface of the droplet, as the interface is now represented by a smoothly varying 3D field. A given liquid droplet has a conserved volume V₀, but may adopt different shapes. A conventional approach to handling this volume constraint is to introduce a Lagrange multiplier, λ, and add the term λ(V[ϕ(r)] − V₀) to the excess energy ΔE, where V[ϕ(r)] is the field-based definition of volume, expressed as a linear function of ϕ. The Euler–Lagrange (EL) equation is then derived and solved simultaneously with λ.¹² While convenient, this approach leads to an undesirable coupling between the bulk properties in the 3D phase field and the interface. Appendix E outlines the issues when a simple expression is used for V[ϕ(r)]. An alternative strategy involves directly minimizing ΔE[ϕ(r)] without introducing λ, using the so-called “model-B dynamics,” which conserves the mass of the droplet. In this formulation, a fourth-order differential equation in ϕ(r) must be solved to describe the mass-conserving flow dynamics.^{3,6,8,13–19} The present paper addresses the equilibrium, static problem only.

Can the volume-constraint issue be resolved in a simpler manner? Here, a non-linear expression for V[ϕ(r)] is introduced in combination with an energy-penalty approach, without the Lagrange multiplier and the use of model-B dynamics. The fundamental idea of using V[ϕ(r)] is unchanged: in the sharp interface limit, it still precisely defines the volume. The non-linear form, however, ensures the decoupling of the bulk properties from the interface, an important feature to avoid the pitfalls associated with the linear form. The approach can be straightforwardly implemented within an energy-penalty framework described in Section II B.

Thirdly, the usefulness of the model is demonstrated through a comparison with known results. As inherited from the Cahn–Hilliard (CH) model, a finite interface width κ arises in most phase-field models. In the simplest case, an interfacial profile ϕ(z) has the form [1 − tanh(−2z/κ)]/2 where z is the distance from a flat interface. When applying the model to real-world problems, it becomes necessary to push the asymptotic limit κ → 0. This introduces a numerical-precision challenge: κ is always finite and competes with the mesh size of the computational grid. How to extrapolate useful information from finite κ-calculations is explored in Section II C. The extrapolated surface excess energy is benchmarked against known theoretical results using the simple example of a spherical cap wetting a flat wall surface, as described in Section II D.

The flexibility of the algorithm presented here lies in its minimal input requirements and general applicability to complex geometries. All one needs to do is to provide the function ϕ_S(r) to describe the solid surface and an initial, rough guess of the liquid shape ϕ(r). An example is provided in Section II E, demonstrating the convergence of the algorithm in modeling the reverse clamshell conformation. No specification of the wetting angles is needed, as these naturally emerge in the optimized final result. The 3D computation can be automated and carried out on a single-thread computer within a manageable amount of time. This makes the method particularly well-suited for high-throughput studies or exploratory investigations of wetting phenomena on arbitrarily shaped substrates.

The success of the algorithm presented in this paper is further validated through a comprehensive study of the wetting behavior of a liquid droplet on both the exterior and interior surfaces of a cylindrical solid. The subject has been extensively studied over the past few decades. Typically, analytical and/or numerical solutions are obtained using an axisymmetric model for the axisymmetric droplet conformation, while fully 3D numerical solutions can be performed using the Surface Evolver triangulation tool.^20–28 Section III covers a quantitative comparison between the state diagrams and phase stabilities calculated here and those reported in the literature. Some technical details are presented in Appendices F, G, and H. Generally, there is excellent quantitative agreement.

By fulfilling the three-fold task, I hope readers recognize the usefulness of this algorithm in addressing complex wetting phenomena with broader applicability. It can be developed into a powerful tool, offering a viable alternative to Surface Evolver for the quantitative calculation of 3D conformations of wetting liquid droplets.

II. Basic formalism

A. Wetting model

The equilibrium shape and total surface energy of a liquid droplet in contact with both solid and gas phases can be determined from minimization of the excess surface energy,¹

ΔE = γ_LGA_LG + (γ_SL − γ_GS)A_SL, (1)

with respect to the 3D liquid shape variation. Here, γ_LG, γ_GS, and γ_SL are the interfacial tensions between liquid–gas (LG), gas–solid (GS) and solid–liquid (SL) interfaces, respectively. The total interface areas between liquid–gas and solid–liquid are A_LG and A_SL on curved surfaces. Recall Young's relation between the liquid–solid wetting angle θ and surface tensionswhich is given for materials making up the system. To calculate a droplet shape, the energy

ΔE/γ_LG = A_LG − cos θA_SL (3)

is minimized with respect to shape deformation, with a fixed internal volume V₀ as the constraint.

In this work, a scalar phase field ϕ(r) is used to model the droplet profile, shown schematically in Fig. 1. It has the value +1 in the spatial regions where the liquid is present and transits to 0 otherwise. The shape of the droplet is described by the interface profile of the function ϕ(r), used to minimize the excess surface energy. The presence of a solid wall surface is by a phase field ϕ_S(r) that typically has the value 1 inside the wall; ϕ_S(r) is not varied in the calculation. The gaseous state is described by the third phase field, ϕ_G(r) = 1 − ϕ(r) − ϕ_S(r), which is written in terms of the first two. Smooth interfaces typically have a characteristic length scale κ.

Fig. 1 Basic concepts in this study. The 3D liquid-droplet shape is represented by ϕ_L(r) = ϕ(r), the primary focus of the calculation. The solid profile is fixed at a pre-specified ϕ_S(r) and the gas profile is deduced from the relation ϕ_G(r) = 1 − ϕ_L(r) − ϕ_S(r). Far from the interfaces, each of these fields approaches the idealized bulk values, as illustrated in the sketch. A three-component phase-field model is used for ϕ, ϕ_S and ϕ_G.

The form of the phase-field functional used here is adapted from the three-component theory,¹¹ reviewed in Appendix A. For the current wetting theory, from the derivation in Appendices B and C,

where

U_LG(ϕ) = 6ϕ²(1 − ϕ²) − U_SL (6)

U_SL(ϕ) = 6ϕ²ϕ_S[3(1 − ϕ_S) − 2ϕ]. (7)

In problems involving an adsorbing solid surface, ϕ_S(r) must be prespecified. Examples of ϕ_S(r) are provided in Appendix D.

It is obvious that the model reduces to the CH model for the LG system when ϕ_S = 0. It can also be shown that the model gives rise to the CH energy for the SL interface when ϕ_G = 0. In the far-field region, the conditions ∂U_LG/∂ϕ → 0 and ∂U_SL/∂ϕ → 0 are both satisfied independently for any value of ϕ_S. Physically, this ensures that the entire liquid droplet remains stabilized within the spatial region of interest. The above condition also reflects the requirement that ϕ_S(r) must satisfy its own EL equation at a pure GS interface. It is a delicate yet crucial point, discussed further in Appendix C, and is essential for constructing a physically consistent and properly functioning wetting phase-field model.

B. The EL equation

In this study the volume of the liquid droplet is fixed at V₀ through the constraintThe particular nonlinear-ϕ formula avoids the pitfalls commonly encountered when using the more conventional linear formula .

Following the formalism in the previous section, and incorporating an energy-penalty term to enforce the volume constraint above, a cost function is constructed,

which is to be minimized with respect to ϕ, for a large parameter Λ. The variation of the functional J with respect to ϕ, i.e., δJ/δϕ, gives rise to the EL equation,where the symbol ′ is used to denote derivatives of U_LG and U_SL with respect to ϕ.

Various computational methods can be employed to tackle the task of solving such a differential equation. Here, I adopt a straightforward approach, using the five-point central finite-difference scheme to approximate the Laplacian operator, and utilizing the limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) algorithm to solve the numerical system.²⁹ The goal is to demonstrate the many useful features of the model after incorporating the volume constraint. Naturally, this procedure leaves room for further improvement through more advanced numerical techniques.

What about the undesirable numerical artifacts discussed in Appendix E concerning the linear-ϕ form of the volume constraint? From the conceptional standpoint, the EL equation is a variational derivative, which encapsulates the same physical principles as the Young–Laplace equation, commonly expressed in terms of surface curvatures and, crucially, the Laplace pressure.¹ The Laplace pressure is an interfacial property and should not extend into the bulk region, far from the interface. However, when a linear-ϕ constraint is used, the resulting EL equation incorrectly propagates this interfacial information into the bulk. In contrast, the ϕ(1 − ϕ) factor in eqn (10) suppresses the last term in the bulk region, making it active only near the interface, where it effectively reproduces the Laplace pressure. The coefficient, i.e., the product of penalty Λ and the constraint term in brackets […] is proportional to the Laplace pressure (a finite quantity). When a large Λ is employed, the constraint term is driven to vanish numerically, yet the product remains finite, thereby preserving the physical meaning of the pressure while avoiding unphysical behavior in the bulk.

The fact that the linear-ϕ constraint leads to unphysical behavior in a phase-field wetting theory was already recognized in ref. 5. There, a ϕ(1 − ϕ) term was empirically introduced into the EL equation, which aligns with the nonlinear form suggested above. However, the coefficient of their term was introduced in a different manner, without incorporating the physical considerations outlined here. As a result, the method in ref. 5 exhibited significant variation in the droplet volume.

C. Numerical treatment

The droplet volume defines the length scalewhich is used as the basic unit to reduce length-related variables and parameters. In the remainder of this section, R₀ = 1 (equivalently, V₀ = 4π/3) is taken. In some later sections, V₀ (or R₀) is explicitly reintroduced for consistency. To enforce the volume constraint, the calculation was done in 4 periods, with the value of Λ at 10, 10², 10³, and an eventual Λ = 10⁴ to improve the convergence. The final computation maintains a volume to within a relative error of approximately 10⁻⁴, which is sufficiently precise for most practical applications.

The spatial domain is discretized using a 3D grid with equally spaced nodes, characterized by a mesh size h. A competing length scale is the interfacial width κ. Ideally, both h → 0 and κ → 0 are desirable to achieve an accurate calculation of the droplet shape. A common question arises: What values of h and κ should be chosen for the numerical computation?

The answer is not straightforward. For a given h, one can progressively reduce κ to compute a sequence of solutions that approach the thin-interface limit κ → 0. However, when the interfacial width becomes smaller than the mesh size h, the grid scheme introduces significant errors due to the sharp variations in ϕ(r) across the interfacial region, which cannot be adequately resolved using finite differences.

The mesh size h cannot be too small either. The number of nodes directly impacts both memory usage and computational speed. Although modern high-performance computing environments offer substantial memory resources, the limits are still there. Consequently, there are practical constraints on how small h can be in actual computations.

Rather than discarding valuable results obtained from finite κ and h, it is more informative to analyze and understand the associated errors. Typically, the leading terms are,

ΔE(κ, h)/γ_LG = f₀ + f₁κ + f₂h²/κ + f₃κ²… (12)

The third term stems from the contributions involving U_LG and U_SL in (4) and (5). Generally, the integral carries an error of order , which, when combined with 2/κ, results in the h²/κ behavior. The procedure then becomes clear: after obtaining the full set of ΔE(κ, h) for various combinations of κ and h, the numerical results can be fitted to (12), for the purpose of extrapolating f₀. This value is then considered as the final numerical result for an extremely thin interface, independent of h and κ. In the following presentation, most of these extrapolated values are depicted as yellow-filled circles in plots.

D. Validation

To validate the model, numerical treatment and interpolation method, the trivial example of a droplet wetting a flat surface is taken as the first test case. The wall profile is modeled by a tanh function, taken from Appendix D. A series runs was conducted for h = 0.005, 0.010, and 0.020. In each run with a given h, a range of κ was adopted, sequentially reduced from a moderate to a low value, and applied in stages. Then, the profile converged from a larger κ serves as the initial guess for the subsequent calculation with a smaller κ, facilitating smoother convergence.

For selected wetting angles θ, the fitted results [f₀ from (12)] are shown in Fig. 2(a) as yellow circles. In the background, the exact solution, derived from the analytic expressions for the surface area and volume of a spherical cap, is represented by the solid curve. From the raw data in Fig. 2(b), it is evident that there is a significant κ-dependence; therefore, directly adopting the calculation from any given κ would give rise to misleading conclusion. The excellent agreement between the interpolated data points in Fig. 2(a) and the exact analytical solution confirms the validity of the procedures used to obtain the final data. Further validation against other known results is provided in the next section.

Fig. 2 Benchmarking the final numerical data from the current algorithm against the exact solution, when a droplet of volume V₀ wets a flat surface. In plot (a), the circles represent the final computed data of the excess energy ΔE, which is reduced as ΔẼ = ΔE/(4πR₀²γ_LG) and plotted as a function of selected −cos θ where θ is the wetting angle and R₀ is defined in (11). The exact analytic result is shown as the solid curve for comparison. In plot (b), unfilled circles (h = 0.005), squares (0.01), and diamonds (0.02) represent the raw data, ΔẼ(κ, h), which are used to determine the asymptotic values at κ = 0, plotted on the left vertical axis as yellow-filled circles. The lines in (b) represent the full fitted functions from (12), with the f₂-term removed after the fitting.

The entire calculation, for a given h, was sequentially performed on a single core CPU (AMD, Zen 4). For the three different h cases, each including approximately 12 κ values, the clock times were roughly a few days, a few hours, and about a hour, respectively.

E. Emergence of the solution

The starting point of the calculation is an arbitrary proposed blob of liquid, represented by ϕ(r) within the confined space of a given surface. Nearly any shape for the blob, even with a volume far from V₀, can serve as the initial guess, as long as it has the correct morphology. In the example shown in Fig. 3, a phase-field profile for ϕ_S(r), representing a hollow tube of solid, is specified according to Appendix D. The initial guess for ϕ(r) is an off-centered ellipsoid, which is shown in the cross-sectional view (left) and the axis-cutting view (right) in Fig. 3(a). Notably, in this particular example, there is significant overlap between the initial liquid droplet and the interior of the wall. Additional examples are given in Appendix F.

Fig. 3 Example of the calculation process for determining the reverse-clamshell profile of a liquid droplet within a confined hollow tube. The data is taken from a single run of the case κ/R₀ = 0.15 and h/R₀ = 0.02 for cos θ = 0.5, where the tube radius R/R₀ = 1.1. Four snapshots, labeled (a) through (d), are taken at different optimization time steps, 0, 8, 64 and 361, respectively. Orange and navy-blue colors represent the droplet and the wall, respectively. In each snapshot, the left and right panels show views of the cross-section, cut perpendicular to the tube axis and through the droplet center of mass, and the plane containing the tube axis and the droplet center of mass, respectively. The reduced excess energy, ΔẼ = ΔE/(4πγ_LG), and the cost function, J̃ = J/4π, are plotted as red circles and white circles, respectively, in (e) as functions of the minimization step t.

The overlap is eliminated through the minimization of the free-energy model, as it leads to a large, undesirable interaction energy in the U_SL term. After a few steps of the minimization, the overlap is substantially removed from the wall interior. This is illustrated in in Fig. 3(b) where a clear reduction of ϕ(r) in the wall space is shown in step 8 of the LBFGS updates. Note that the overall liquid shape is not fully optimized in the early stages of the process.

The liquid shape undergoes significant improvement, particularly in the contact regions between the droplet and the wall, as shown in step 64 of the LBFGS minimization, Fig. 3(c). One notable feature is the considerable expansion of the droplet-wall contact area, indicating that the liquid shape is adapting more smoothly to the wall surface. At the same time, the profile near the axial region of the tube experiences much smaller distortions, as it has not been significantly influenced. In this early stage of minimization, an approximate, yet recognizable, liquid shape begins to emerge. While the overall shape is not yet fully optimized, this stage marks a meaningful step toward a more physically accurate configuration.

The remainder of the calculation progressively refines the phase-field profile ϕ(r) as the minimization continues. By step 360, the minimization reaches convergence, as determined by the LBFGS convergence criterion, signaling that the solution has stabilized. At this point, the calculation is terminated, and the final, fully optimized 3D liquid shape is obtained, which is presented in Fig. 3(d).

The minimization trajectory of the total cost function J and the reduced energy ΔE/γ_LG are both displayed in Fig. 3(e) represented by circles and red circles, respectively, as a function of the minimization step t. The curves illustrate the progress of the minimization process over time, offering insight into how the system evolves toward its optimal configuration. When entering the calculation, at the beginning, J and ΔE/γ_LG differ from each other, due to the large discrepancy between the initial guess and the desired volume. The gap is quickly reduced. During the later stages of the calculation, the focus shifts primarily to optimizing the liquid shape, with the energy terms becoming increasingly refined. Arrows on the plot indicate the specific steps at which the four snapshots in Fig. 3(a)–(d) are taken, marking key stages when the calculation is driven to convergence.

The algorithm is quite stable and can accommodate a wide range of initial conditions for the liquid droplet, provided that the initial morphology is reasonably chosen and remains within the boundaries of the computational domain. In particular, when appropriate values of h and Λ are used (for example, h/R₀ = 0.02 and Λ = 10²), the algorithm generates the droplet profile in interactive, real-time scales, allowing for rapid decision-making in subsequent steps. This versatility makes the method especially valuable for modeling complex systems where the initial state may be uncertain or not well defined.

III. Results on cylindrical geometries

In this section, the usefulness of the above algorithm is demonstrated by calculating the wetting surface energies when a droplet wets the exterior of a cylinder and the interior of a hollow cylinder of radius R. The former is sometimes referred to as a “fiber”, and the latter as a “tube”, “channel”, “pipette”, or “capillary tube”. A detailed comparison with previously known results, obtained through other theoretical and numerical approaches, is presented. This comparison includes the state diagram and state stability analyses. Five states are described below: barrel (B), clamshell (C), reverse barrel (RB), plug (P), and reverse clamshell (RC). They are schematically represented in Fig. 4.

Fig. 4 The five wetting states, stable in parameter regions discussed in Sections III B and III C, are illustrated for a liquid droplet on solid-cylinder and hollow-cylinder surfaces. These are presented through conceptual plots (upper panel) and corresponding 3D snapshots obtained from the solutions to the phase-field model (lower panel). To produce the snapshots in (a) and (b), R/R₀ = 1 and cos θ = 0.7 are used; in (c)–(e), R/R₀ = 1.2 and cos θ = 0.9 are used.

There is a large body of literature on the physical properties calculated from axisymmetry theory (initially analytic, then transitioning to numerical methods for the complete solution) and numerical computations (primarily using Surface Evolver²⁰).^{23,26–28,30–37} These studies have provided valuable insights, some of which have been verified by carefully designed experiments.^{22,26,27,32,34,36–39} In some references, the P state is referred to as a “slug”, the RB state as a “ring” or “annulus”, and the RC state as an “adhered drop” or “droplet”.

A. Surface energies

Examples of the five branches of the surface energy calculated from this study are shown in Fig. 5 by yellow symbols. Plot (a) contains the energies for B and C [Fig. 4(a) and (b)] at cos θ = 0.7, and plot (b) RB, P, and RC [Fig. 4(c)–(e)] for cos θ = 0.95, as functions of the cylinder radius R.

Fig. 5 Reduced surface free energy ΔẼ = ΔE/(4πR₀²γ_LG) for the five states in Fig. 4, for wetting on (a) solid-cylinder and (b) hollow-cylinder surfaces, respectively. The solid curves, representing the B state in (a) and the RB and P states in (b), correspond to exact solutions obtained from axisymmetric models, independent from the current approach. The symbols filled by yellow are data points that have been extrapolated from the κ- and h-dependent data. The dashed curves are energies of the unstable states. The cyan symbols indicate stability limits of these states. In all plots, the data error bars do not exceed the symbol size and R₀ = (3V₀/4π)^1/3.

For comparison, the black curve in Fig. 5(a), behind circles for B, is generated from an independent axisymmetry theory (see Appendix G), which is treated here as a precise numerical result. The numerical solutions, from both the phase-field theory and the axisymmetry theory, can extend to larger region of R, as shown in the figure by dashed symbols and dashed curves. The stability limit, depicted here by the cyan circle, is discussed in the next subsection.

The calculated energy for the C state, represented by squares in Fig. 5(a), is completely asymmetric in shape; hence, there is no simpler theory to compare it to. The green curve in the background is drawn based on the squares to guide the eye. The C branch terminates at a lower R, indicated by the cyan square, which is discussed in the next subsection.

For the hollow-cylinder case, the example shown in Fig. 5(b) is for cos θ = 0.95, which demonstrates three competing energy branches, RC, P, and RB. The axisymmetric P state has an exact analytic solution, based on the assumption of a cylinder droplet in the tube, capped by two negative spherical caps, as described in Appendix H and ref. 23. The analytic solution is shown by the red curve in the background, where the yellow diamonds represent data points calculated from the 3D phase-field theory.

The calculated RB branch in Fig. 5(b) is represented by the circles. In the background, a black solid curve shows the precise numerical results from an independent axisymmetry theory [Appendix G]. This energy branch has two stability limits, high and low, represented by the cyan circles. The droplet shape of the RC state is completely asymmetric. The squares in Fig. 5(b) are the data points produced from the current phase-field study. The green curve in the background is drawn to guide the eye. The RC state terminates at the cyan square. The discussion on the stability limits can be found in Section III C.

In both plots, only minor differences between the red/black curves and the yellow symbols can be observed. These differences can serve as a guideline for error estimates in further applications, such as using the phase-field data to predict state diagrams in the next section. The close agreement between the data calculated here and independent theories further validates the algorithm introduced in this study. This also supports the reliability of the data points for the asymmetric states, C and RC. The axisymmetry of the B, RB, and P states emerges naturally from the current numerical study; it is not an imposed condition on the phase-field theory.

B. States outside a solid cylinder

The last subsection described one example of the energy curve for cos θ = 0.7 in Fig. 5(a). A phase transition radius can be identified when the two energy branches intersect, indicated by an arrow in the plot. This data point enters into Fig. 6(a) as a yellow circle. Other phase-field data are also represented by yellow circles in the same plot. The dotted line in the background is drawn based on these data points.

Fig. 6 State diagram and stability limits of the B and C states. The radius of the cylinder R is reduced by R₀ in (11). In plot (a), the state boundary is determined by comparing the surface energies. The dotted curves in all plots represent this boundary, interpolated from known data points. In plot (b), the circles represent the B state stability limit determined using the off-center test (see text), and the solid line corresponds to the stability limit predicted by the inflection-point theory²¹ (see Appendix G). In plot (c), the circles indicate the stability limit of the C state, determined by locating the closure of the wetting points in the central cross-section (see text). In all plots, the crosses, diamonds and squares correspond to data reported in ref. 22, 25 and 26, respectively. See a summary in Table 1.

Table 1 Summary of data plotted in Fig. 7 (upper tablet) and Fig. 8 (lower tablet). The first column lists the references by the first-author basis. The third and fourth columns list the theoretical and numerical methods used, and the properties determined. Abbreviations used are: AS (analytic solution), AT (vanishing angle test, see Sections III B and III C), AX (axisymmetry numerical method, see Appendix G), BL (B state limit), CL (C state limit), CM (off-center-of-mass test), DP (deflection-point analysis), PF (phase field theory), PL (P stability limit), RBL (RB stability limit), RCL (RC stability limit), SB (state boundary), SE (Surface Evolver), and YL (axisymmetric Young–Laplace)

Ref.	Symbols in plots	Theory	Properties
McHale^22	× (a)	YL & SE	SB
	× (b)	YL & DP	BL
Eral^25	◊ (a)	YL & SE	SB
	◊ (b)	YL & DP, SE & CM	BL
	◊ (c)	SE & CM	CL
Chou^26	□ (a)	SE	SB
	□ (b)	SE & CM	BL
	□ (c)	SE & AT	CL
This work	○ (a)	PF	SB
	○ (b)	PF & CM	BL
	○ (c)	PF & AT	CL
	Solid (b)	AX & DP	BL
Collicott^23	+ (a) and (d)	AS, YL, & SE	SB
	+ (b)	AS	PL
	+ (c) and (f)	SE	CL
	+ (e)	YL	Upper RBL
Liang^27	□ (a) and (b)	SE	SB, PL
	□ (c) and (f)	SE & AT	RCL
Lv^28	◊ (a), (c), (d) and (f)	YL, SE	SB, RCL
	◊ (e)	YL	RBL
This work	○ (a), (b), and (d)	PF	SB, PL
	○ (c) and (f)	PF & AT	RCL
	○ (e)	AX & DP	RBL
	Solid (b)	AS	PL

---

The B–C state boundary has been determined by McHale and Newton.²² Subsequently, two independent studies revisited the problem and provided more detailed datasets.^25,26 In all these studies, the B state was computed using the axisymmetry theory²² and the C state was obtained using Surface Evolver to model the full 3D conformation. Their state boundaries are reproduced in Fig. 6(a), crosses corresponding to tabulated data from ref. 22, whereas diamonds and squares representing data extracted from figures in ref. 25 and 26, respectively. Notably, all four sets of symbols show strong agreement with each other.

A conjecture was proposed that, for a given contact angle θ, the B state reaches a high-R stability limit, which manifests as the inflection point in the axisymmetry profile of the liquid droplet.²¹ A more detailed discussion can be found in the supplementary materials of ref. 25. This stability limit is depicted in Fig. 3 and 5 of ref. 25 and 26, respectively. The same assessment is also reproduced in the present work, using the axisymmetry theory proposed in Appendix G. All results converge onto the same solid curve shown in Fig. 6(b).

An insightful numerical test for validating the inflection-point theory was proposed by Eral et al.²⁵ In this approach, numerical solutions for the B state, obtained in their case using Surface Evolver, are perturbed by displacing the droplet's center of mass slightly away from the central axis, thereby introducing an axisymmetry-breaking deformation. The surface energies of the perturbed and unperturbed (axisymmetric) configurations are then compared. If the energy difference becomes negative, the B state is deemed to have reached its stability limit. In the present work, the same testing procedure is employed, but using phase-field theory instead. An illustrative example is shown in Fig. 7(a) for cos θ = 0.7, where the arrow marks the identified stability limit; this value is also plotted as a yellow data point in Fig. 6(b). By aggregating all data points from Fig. 6(a) and (b), a consistent agreement is observed between theoretical predictions and numerical results across both diagrams, from all different approaches.

Fig. 7 (a) Off-center test for the B state and (b) closure limit for the C state. For a given θ, the droplet in a B state [see (c) for a central cross section example] is axisymmetric and coaxial with the cylinder axis, possessing a surface energy ΔE_B. The off-center test forces the droplet's center of mass to displace from the axis by 10%R [plot (d)] and in the meantime, the current algorithm produces a full 3D nonsymmetric profile with an associated surface energy ΔE_off. The point at which ΔE_off becomes lower than ΔE_B marks the B state stability limit, indicated by an arrow in the plot. In (b), the contact angle α² [see the central cross section view in (e)] is plotted as a function of R/R₀. The radius at which α = 0 is identified as the C stability limit, and is indicated in (b) by yellow circles. All energies are reduced by ΔẼ = ΔE/(4πR₀²γ_LG).

The C state reaches its stability limit at low values of R when the clamshell profile (see Fig. 7(e)) can no longer be maintained. In this study, the angle α, defined as the angle between the contact point (marked by a red circle) and the vertical mirror symmetry line in the central cross-section of the C profile, is measured as a function of R. By gradually reducing R in a stage-by-stage manner to ensure numerical convergence at each step, the value of α is calculated and analyzed. The closure of the wetting edges, indicated by the vanishing of α, is taken as the criterion for identifying the stability limit of the C state. Fig. 7(b) presents α² as a function of R, showing a linear behavior that extrapolates to the yellow circles at α² = 0 for given cos θ. In contrast, a direct plot of α (not shown) exhibits clear non-linear behavior. The extrapolated values from this numerical analysis are used to determine the stability limits of the C state, which are reported in Fig. 6(c).

A similar numerical procedure was employed in ref. 26, where the closure of the C state was identified by tracking the height of the contact point (indicated by the red dot in Fig. 7(e)). Their data, extracted from Fig. 5 of that study, are represented as squares in Fig. 6(c). In a different approach, the center-of-mass perturbation testing method, used to test the B state stability above, was applied in ref. 25 to assess the C state stability. The corresponding data, taken from their Fig. 3, are shown as diamonds in the same plot. While all three data sets exhibit a qualitatively similar trend, they differ quantitatively at finer scales.

C. States inside a hollow tube

The current theory predicts three wetting states of a liquid droplet confined within a hollow tube at cos θ = 0.95, in the vicinity of R/R₀ = 1: RC, P, and RB. As illustrated in Fig. 4, two of these are axisymmetric about the tube axis. These three states have competing surface energies ΔE in the range of R plotted in Fig. 5(b).

Region 1. Fig. 8 presents the state diagram and stability limits of relevant wetting states on the R–θ plane. For larger contact angles, corresponding to −cos θ ≳ −0.90, the RB energy branch lies above RC and may even become unstable. In this regime, only two energy branches, P and RC, intersect at the state boundary. The circles shown in Fig. 8(a) delineate the boundary between the P and RC states, as determined from the present study. These results can be compared with the same boundary indicated by crosses, reported in ref. 23, squares, reported in ref. 27, and diamonds, reported in ref. 28.

Fig. 8 State diagram of a droplet inside a hollow tube as a function of tube radius R and contact angle θ. Plots (a)–(c) provide an overview of the state diagram and stability limits of P and RC. The wetting regime is further detailed in (d) where three states occupy different parameter regimes. The stability limits of RB is presented in (e) and that of RC in (f). In all plots, yellow circles represent data obtained from the current study, squares are readoff from Fig. 9 of ref. 27, diamonds are readoff from Fig. 5 and 6 of ref. 28, plus symbols are readoff from Fig. 15 of ref. 23. The solid curves in (b) are derived from the exact solution in Appendix H, and those in (e) from high-accuracy numerical solutions of the axisymmetry model presented in Appendix G. The circle symbols in (b) and those near the upper RB limit in (e) are determined directly from the droplet profiles produced by the phase field theory; the circle symbols near the lower RB limit in (e) are obtained from the off-center test. The circles in (c) and (f) are generated using the contact-point closure extrapolation method described in Fig. 9(b). In all plots, the dotted lines indicate the interpolated state boundaries in (a) and (d). See a summary in Table 1.

The P stability limit can be traced to the physical picture of geometrically maintaining the conformation within a hollow tube.²³ The solid curve in Fig. 8(b) is obtained from the exact analysis presented in Appendix H. Beyond the stability threshold, the numerically computed phase-field profile undergoes a discontinuous transition to a metastable hollow-core RB profile. The yellow circles in this plot indicate the largest R for which the P profile remains stable. For comparison, the squares in the same plot represent the P stability limit adapted from ref. 27, showing reasonable agreement with the current assessment.

The RC stability limit, shown in Fig. 8(c) and (f) as yellow circles, is determined using a numerical procedure analogous to that employed for the C state. For each specified value of θ, the contact angle α in the central cross section, illustrated in Fig. 9(e), is measured across a sequence of decreasing R runs, each case allowing sufficient time to reach equilibrium. The resulting data for α² are then analyzed to extrapolate the point at which α² vanishes, corresponding to the intersection with the R/R₀ axis. These extrapolated limits are represented by the yellow circles in Fig. 9(b).

Fig. 9 (a) Off-center test for the RB state and (b) closure limit for the RC state. For a given θ, the droplet in a RB profile [see (c) for a central cross section example] is axisymmetric and coaxial with the cylinder axis, with a surface energy ΔE_RB. The off-center test forces the droplet's center of mass to displace from the axis by 10%R [plot (d)] and in the meantime, the current algorithm produces a 3D nonsymmetric profile and its surface energy ΔE_off. The arrow indicates the RB stability limit. In (b), the square contact angle α² [see the central cross section view in (e)] is plotted as a function of R. The radius corresponding to the closure of the contact points, when α² = 0, is considered as the RC stability limit. These values are shown in (b) by yellow circles.

Liang et al. also determined the RC stability limit using a similar criterion, based on the observation of the closing of the upper “edges” in Fig. 9(e), obtained through Surface Evolver calculations.²⁷ For comparison, their results are represented by square symbols in Fig. 8(c) and (f). Also using Surface Evolver as a tool, Collicott et al.,²³ as well as Lv and Hardt²⁸ determined the RC stability limit by observing the point where the computation is no longer sustainable to produce the correct morphology. The crosses are read-off from Fig. 16 of ref. 23, and diamonds from Fig. 6 of ref. 28. A reasonable agreement can be seen between these calculations and the present phase-field data indicated by the yellow circles.

Region 2. This region in the state diagram primarily lies below −cos θ ∼ −0.90, corresponding to small wetting angles θ. Shown in the state diagram in Fig. 8(d), the RB state emerges as the energetically favored conformation over a significantly broad range of R. This was also suggested in ref. 23; for comparison, their determination of the phase boundary is shown in Fig. 8(d) by crosses, which agree well with the yellow circles, given the computational difficulties encountered for small θ in both approaches. However, this stability region was not reported in ref. 27, where the P–RC boundary was extended to lower values of θ, as indicated by squares in the plot. The fact that the RB energy can be lower than those of P and RC is demonstrated in Fig. 5(b), where the intersections of the energy branches are marked by two arrows.

The RB state is characterized by a hollow-core profile. Solutions obtained from the current phase-field model, as well as those derived from the axisymmetry theory detailed in Appendix G, include large-R configurations. These are represented in Fig. 5(b) by the long-dashed curve and the dashed circles. A stability limit at large R truncates the RB energy branch, indicated by the cyan circle at high R/R₀.

Although the original inflection-point theory was developed for analyzing the limit of the B state,^21,25 here it is assumed also valid for the RB state. The blue curve in Fig. 8(e) shows the corresponding numerical result, obtained by analyzing the axisymmetric RB profiles, which is discussed in Appendix G. In addition, using the testing technique introduced in ref. 25, the stability limit can be evaluated by examining the surface energy of an off-centered RB configuration. Such a profile is generated by displacing the droplet's center of mass slightly from the central axis of the tube (see Fig. 9(d)). The resulting energy difference between the centered and off-centered configurations can then be analyzed; an example for cos θ = 0.8 is shown in Fig. 9(a), where the stability threshold is indicated by the arrow. Returning to Fig. 8(e), the results of the off-center test, represented by the yellow circles overlaid on the blue curve, exhibit good agreement with the inflection-point theory.

The RB state also exhibits a second stability limit at low R. As R decreases, the hollow core of the RB profile gradually shrinks and ultimately closes at this low-R stability threshold. Based on the axisymmetry theory presented in Appendix G, the closure limit can be identified through numerical solutions of the droplet profile, and is depicted by the solid curve in Fig. 8(e). Alternatively, this stability limit can also be determined directly from the phase-field model. The yellow circles in the same figure represent results obtained from phase-field calculations. The shaded RB region in Fig. 8(d) falls inside the stability region presented in Fig. 8(e).

Lv and Hardt have conducted a commendable and thorough analysis of the stability limit of the RB state,²⁸ exploring a wide parameter space in θ and R. They numerically solved the Young–Laplace equation, which in this case is a differential equation derived from the exact axisymmetric theory. For comparison purposes, the diamond symbols in Fig. 8 are plotted based on data adapted from their Fig. 5. In general, the stability limits of the RB state determined here are in good agreement with their calculations, although minor systematic deviations can be observed in plot (e). Their low RB limit was obtained by analyzing the merging point of the two branches of RB solution to the Young–Laplace equation. Note that, the three approaches, inflection condition, off-center-of-mass test, and merging of two solution branches, all agree well for the determination of the RB stability.

Other states. During the search for the RB state, starting from an initial condition with an azimuthal-angle dependence about the tube axis, three other asymmetric states were obtained as solutions of the Euler–Lagrange equation of the phase-field model. These are illustrated in Fig. 10, each containing holes directed along the tube's axial direction. The solutions display 2-, 3-, and 4-fold rotational symmetries. When angular perturbations in the azimuthal direction are introduced, they converge back to the same profiles.

Fig. 10 Snapshots of the three asymmetric states, in (a), (b), and (c), discovered for droplet wetting in a hollow cylinder of radius R/R₀ = 1.7 and cos θ = 0.5, from the 3D phase-field calculation.

The center of mass of the droplets in the figure remains on the tube axis. However, when the center of mass is slightly shifted to produce an off-axis solution, the surface energy ΔẼ is reduced. This indicates that these are unstable profiles prone to off-axis perturbations, in the same way that RB states behave beyond their lower boundary presented in Fig. 8(e). In the cos(θ) − R parameter space, no stable regions for these states have been found, so far.

IV. Summary

This paper focuses on developing a phase-field model and its associated computational algorithm, tailored for the quantitative calculation of excess surface energy and droplet profiles when a 3D, asymmetric droplet wets a solid surface. The challenge of maintaining a fixed droplet volume is effectively addressed through the introduction of a nonlinear formalism. The approach to the limit of vanishing interfacial width, required for modeling physical systems where the droplet scale is much larger than the interfacial thickness, is managed through a proposed data extrapolation technique. The validity of the algorithm is demonstrated by recalculating energy branches, state diagrams, and stability limits for several known states, for which exact solutions, precise numerical results, or three-dimensional computations using Surface Evolver are available. This paper offers an integrated presentation of a range of physical and numerical ideas from various origins, delivered within a cohesive framework.

Several technical improvements could further enhance the efficiency of the algorithm. For clarity in conveying the core ideas, a finite-difference method is employed in Section II A. However, more sophisticated numerical techniques, such as those developed for related interface problems,^{8,13,18,40–45} could certainly be adopted and implemented within the same theoretical framework. For simplicity, a typical computation performed above used a single-core CPU. The use of distributed computing can offer significant speedup in real time.

The fact that a neural network offers an effective way to represent a function can be further exploited in developing a new approach to solving the Euler–Lagrange equation in this study. For example, the phase-field model for the fluid membrane problem, which shares many fundamental features with the current surface problem, has recently been solved using a unsupervised machine-learning approach.⁴⁶ In such a framework, the function ϕ(r) can be directly represented by a neural network and the energy model, which is described in Section II A, remains the same. This avoids the need to discretize the computational domain into a grid, thereby eliminating the associated storage and resolution limitations. A recent advancement in solving the wetting problem inside a hollow tube is presented in ref. 47, where the 2D surface profile is embedded within a neural network. The suitability of phase-field models for neural network-based approaches represents a promising direction for future numerical implementations.

The algorithm presented in this paper offers an alternative to the Surface Evolver software. Section III explored droplet wetting on both solid and hollow cylindrical substrates. The current approach can be readily extended to other recently studied wetting systems, such as wetting on two parallel or perpendicular cylindrical surfaces,^48,49 a cone-shaped substrate,⁵⁰ and a junction between a cylinder and a plane.⁵¹ Gravitational effects can be easily incorporated by adding a potential energy term to U_LG in (6). Going further, since the model is based on a multi-component phase-field theory, it is potentially generalizable to wetting problems involving two immiscible fluid droplets on a substrate,²⁴ starting from using multiple ϕ_i in Appendix A instead of 3. The robustness of the theory, coupled with its ease of implementation, can facilitate the applications in a wide range wetting phenomena.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data contained in this paper are available upon reasonable request.

Appendices

A Three-component phase-field model

A fundamental problem is to examine the coexistence of three fluids, described by phase fields ϕ₁(r), ϕ₂(r), ϕ₃(r), respectively. The surface free-energy model has been previously introduced and well documented,¹¹in which the integral spans the entire space of interest. The functional g is designed to follow the Cahn–Hilliard form,
where κ is related to the interface width. The relative surface-tension coefficients Γ₁ = γ₁₂ + γ₃₁ − γ₂₃, Γ₂ = γ₂₃ + γ₁₂ − γ₃₁, and Γ₃ = γ₃₁ + γ₂₃ − γ₁₂. Additionally, the incompressibility condition
is imposed. Note that the functional is symmetric with respect to the rotation of indices 1, 2, and 3.

B Two-component coexistence

One can easily check, by eliminating ϕ₁, that the above model reduces to the standard Cahn–Hilliard functional for the coexistence of fluids 2 and 3,The EL equation of a freely standing interface (that is, no other constraints) followsA special solution to the above, for a flat interface, is the well-known shifted hyperbolic tangent (or sigmoid) profile between the two bulk phases, represented by (ϕ₂ = 1 − ϕ₃ = 1, ϕ₃ = 0) and (ϕ₂ = 1 − ϕ₃ = 0, ϕ₃ = 1),

Here, the flat interface is located at z₀ and z − z₀ is the distance from it. κ is the characteristic interface width.

C Phase-field wetting model

In the current wetting problem, the three-component model is adapted by setting

ϕ₁ = ϕ(r),

to describe the liquid profile,

ϕ₃ = ϕ_S(r),

a pre-defined solid surface, and

ϕ₂ = 1 − ϕ(r) − ϕ_S(r),

the corresponding gas density. Additionally, the surface energy (B1), now having ϕ₃ = ϕ_S(r), is used as the reference for a pure GS interface. The excess surface energy, after the insertion of the liquid droplet, can be finally written as

ΔE = E − E₂₃ = γ_LGA_LG + (γ_SL − γ_GS)A_SL. (C1)

From the three-component model, the LG and SL interface areas are expressed byTwo functions,andhave been introduced.

It is no surprise that the EL equation of the model reproduces the differential equation for LG coexistence when setting ϕ_S = 0. A more subtle property emerges in the spatial region where ϕ ≪ 1 and ϕ_S ≠ 0. In this limit, one can show that the EL equation reduces to the EL eqn (B2), governing the GS interface alone. In other words, to linear order in ϕ, each of the area expressions in (C2) also contains the EL equation for the standalone GS interface. Physically, this means that the formalism must correctly capture the GS energy even before the droplet is introduced, a fine but important detail often overlooked in many other phase-field models for wetting

To proceed further, ϕ_S(r) is assumed to follow the EL equation, (B2), so that the areas can be simplified to

where

U_LG(ϕ) = 6ϕϕ_G(ϕϕ_G + ϕϕ_S + ϕ_Gϕ_S − ϕ_S²) − 6ϕ(ϕ_S − 3ϕ_S² + 2ϕ_S³), (C7)

U_SL(ϕ) = 6ϕϕ_S(ϕϕ_S + ϕϕ_G + ϕ_Sϕ_G − ϕ_G²) + 6ϕ(ϕ_S − 3ϕ_S² + 2ϕ_S³). (C8)

Simplification yields (6) and (7). In the ϕ ≪ 1 regions, i.e., far away from the liquid droplet,
which is a necessary condition for correct description of the GS interface, discussed in the last paragraph. So far, no approximations have been made.

D Solid surface profiles

In arriving at the above, the prescribed profile ϕ_S(r) is assumed to follow the EL eqn (B2) itself. Below, the function forms for the surface profiles used in the text are listed. Among these, the flat-wall surface is an exact solution to the EL problem in (B2). The solid cylinder and hollow tube profiles are approximations for finite κ, which become exact solutions to (B2) only in the asymptotic limit κ → 0.

Flat surface. Assume that a flat solid surface is located at z = 0. The profile

is used in (C7) and (C8). The numerical results are presented in Section II D.

Hollow tube. Assume that the tube axis of the cylinder is located at y = y_tube and points in the z-direction (see Fig. 4) and the tube has a radius R. Let

r = [x² + (y − y_tube)²]^1/2 (D2)

be the distance from the tube axis. ThenNote that inside the tube, ϕ_S → 0, and outside, →1. The profile is used in Section III C of the text. It serves as an approximation for finite κ and becomes an exact solution as κ → 0.

Solid tube. Assume the same tube geometry described above and the same r in (D2). Instead of a hollow tube,

describes a solid tube, which has ϕ_S → 1 inside and 0 outside. The profile is used in Section III B of the text. It serves as an approximation for finite κ and becomes an exact solution as κ → 0.

E Volume constraint

This study is based on a volume constraint, as given in (8), which differs from the commonly usedThe limitations of using the above, in particular, its tendency to produce erroneous bulk limits, are discussed here. Two numerical approaches are typically employed to incorporate the constraint, Lagrange multiplier method and energy penalty method.

Lagrange multiplier. In this method, the interfacial energy from (3) is used in

which is optimized with respect to ϕ and λ. Even in the elementary case of LG two-fluid coexistence, a physical inconsistency arises. In the far-field region, where the gradient term can be neglected, the EL equation becomes,
asymptotically. For a finite-volume problem, in the interfacial region, λ is related to the non-zero Laplace pressure of the system. Instead of the expected ϕ = 0, a residual ϕ = −κλ/24 arises in the bulk, which, for a large calculation box, causes the integral in (E1) inaccurate.

Energy penalty. What happens if one instead adopts the energy-penalty method? The expression becomes

where Λ is a pre-specified, large parameter. In this case, the EL equation takes the form:
For a physically meaningful solution in the interfacial region, the entire last term corresponds to the non-zero Laplace pressure λ.

This leads to the same residual issue as described above. One might argue that with an asymptotically large Λ, the term would shrink such that λ/Λ → 0. However, for a given Laplace pressure λ, set by the system's physical conditions, the far-field value of ϕ remains finite. The consequence is that the entire ϕ(r) profile must be adjusted to satisfy this numerical condition, leading to non-physical outcomes.

The key issue lies in the fact that the linear constraint, (E1), enforces a constant λ (reduced Laplace pressure) across the entire spatial domain. On the other hand, the Laplace pressure is a concept intrinsically tied to the droplet region alone; extending it uniformly into the far-field introduces a conceptual inconsistency caused by the linear constraint. In contrast, the nonlinear constraint introduced in (8) naturally decouples the interfacial and far-field regions, thereby avoiding this problem.

F Initial guess

The algorithm presented here accommodates rough initial guesses for the different states discussed in the text. Fig. 11 illustrates the initial profiles for (a) a droplet on a flat surface wall and (b)–(f) the B, C, P, RC, and RB states on curved surfaces, used in the calculations of this study. There are significant overlaps between the droplets and the wall profiles, which are shown as wall-line crossings in these plots. The contact angles are initially set incorrectly too. In particular, the RB initial guess has pronounced sharp edges. The algorithm automatically corrects these initial errors during the minimization process.

Fig. 11 Central cross-sectional view (left) and axis-cutting view (right) of the initial profiles ϕ(r) used as the inputs for all calculations in this study, plots (a)-(f). The navy-blue regions, some of which are behind the orange-colored droplet regions, represent the profile defined by ϕ_S(r) in Appendix D.

G Axisymmetry theory

The calculation of the surface areas and volume of a droplet in an axisymmetric configuration is a much simpler task. Previous authors have followed an axisymmetry theory, written in terms of the Young–Laplace equation.^21,22,28,38 Here, a numerically oriented axisymmetry model is designed to take a different approach, intentionally. It shares the same principle as the Young–Laplace equation, where the thin interface is described by a 2D surface.

Consider the B state as an example. The tube axis is divided into a series of nodes labeled by i = 0, 1, 2, …, N, where the neighboring nodes are spaced by a distance h. The droplet's axisymmetric shape is approximated by a series of circles of radii r_i, each concentric with the corresponding nodes. The radii of the first and last circles are fixed at r₀ = r_N = R. Assume now that adjacent circles, of radii r_i and r_i+1, are connected by straight lines to form a hollow frustum. Then, the total surface area and volume can be expressed as functions of r₁, r₂, …, r_N−1 and h. This construction is illustrated in Fig. 12.

Fig. 12 Schematic representation of the geometry used in the axisymmetry theory. In plot (a), a series of hollow frustums are used, with adjustable base radii for the calculations of the surface-area and volume elements. These radii are subject to change in the optimization. Plots (b) and (c) illustrate the first three radius nodes near the cylinder exterior and interior surfaces, used for defining the second derivative of the profile numerically. A change of the sign marks the inflection condition.

For given R and θ, the surface energy, which is a function of multivariables r₁, r₂, …, r_N−1 and h, can be optimized, under the constraint that the droplet volume remains fixed at a prescribed value V₀. A multivariable optimization computational package can then be employed for this purpose. Note that if one takes the Lagrange multiplier method at this stage, the desired solution corresponds to a first-order saddle point of the total cost function J in the multivariable space.

The concept of using the inflection point to determine the stability of the B state was first suggested in ref. 22. The second derivative of the profile near the contact point is examined. As R increases for a given θ, near the contact point with the wall, the B profile transits from a convex to a concave shape, while the RB profile transits from concave to convex, sketched in Fig. 12(b) and (c). The value of R that separates these two distinct curvatures is referred to as the stability limit for B in Section III B and for RB in III C according to the axisymmetry theory. As a consistency check, the stability threshold of the B profile determined through this approach precisely matches the off-center-of-mass test, explained in Section III B. This matching was also verified in ref. 25.

H Exact theory for P

For given values of θ and R, the surface profile of the axisymmetric P state can be determined analytically, without any need for adjustable parameters.²³ This profile consists of a cylinder with radius R, capped by either negative spherical caps with a base radius R (for cos θ > 0) or positive spherical caps with a base radius R (for cos θ < 0). See illustrations in Fig. 13. The height H is then a function of V₀. The red solid energy curve in Fig. 5(b) is derived from this model.

Fig. 13 Illustrations of the axis-cutting view of the P profiles used in Appendix H.

The geometric limit of a profile can also be evaluated based on this model. For cos θ > 0, increasing R causes the two negative spherical caps to approach one another until they eventually touch, signaling the P state limit. On the other hand, for cos θ < 0, a limit is reached when H vanishes. These stability limits are represented by the solid curves in Fig. 8(b). The same procedure of touching negative spherical caps was used in ref. 27 for the cos θ > 0 case.

I Surface energies in a hollow tube

In ref. 28, Lv and Hardt presented the surface energies of the RB, RC, and P states for selected θ values in their Fig. 6. Plots (b)–(d) in Fig. 14 provide a comparison between the results from the current study and their calculations, for θ = 30°, 60°, and 90°, respectively. Plot (a) is reproduced from Fig. 5(b) above for completeness and the same figure legends there are used in all four plots. Added to plots (b)–(d), are gray symbols, which represent the data taken from Fig. 6 of ref. 28.

Fig. 14 Reduced surface free energies, ΔẼ = ΔE/(4πR₀²γ_LG), for the P, RB, and RC states in Fig. 4, for wetting on a hollow-cylinder surfaces, for selected θ values. The same legends used in Fig. 5 are followed here. The gray circles, gray squares, and gray diamonds are representative data for the P, RB, and RC states, read-off from Fig. 6 of ref. 28. Dashed lines and symbols represent unstable solutions, tested by the inflection condition (Appendix G) and by the off-center-of-mass numerical test (Section III C); they agree well with the gray dashed symbols from ref. 28, produced there by solving the Young–Laplace equation. The data from ref. 23 is plotted by plus and cross symbols. Plots (a)–(d) correspond to cos θ = 0.95, 0.886, 0.5, 0, respectively.

In ref. 23, Collicott et al. calculated the surface energies of all three states for θ in an increment of 10°. The plus and cross symbols in Fig. 14(b)–(d) are plotted according to the data read-off from their Fig. 8 and 12 for the selected θ values, of the RB and RC energy branches. Going from high-R to low-R, their data terminates at an earlier R/R₀ value on both branches, near the blue symbols. This explains the apparent deviation of their RC and RB stability limits from the current calculation, in Fig. 8(c), (e) and (f).

Acknowledgements

Financial support from the Natural Sciences and Engineering Research Council of Canada and computational resources from the Digital Research Alliance of Canada are gratefully acknowledged.

References

1. P.-G. de Gennes, F. Brochard-Wyart and D. Quéré, Capillarity and Wetting Phenomena, Springer, New York, NY, 2004.
2. J. W. Cahn and J. E. Hilliard, J. Chem. Phys., 1958, 28, 258–267.
3. N. Provatas and K. Elder, Phase-Field Methods in Materials Science and Engineering, Wiley-VCH, Weinheim, Germany, 2010.
4. J. Kim, S. Lee, Y. Choi, S.-M. Lee and D. Jeong, Math. Probl. Eng., 2016, 2016, 9532608.
5. Y. Di, R. Li and T. Tang, Commun. Comput. Phys., 2008, 3, 582–602.
6. R. Borcia, I. D. Borcia and M. Bestehorn, Eur. Phys. J.: Spec. Top., 2009, 166, 127–131.
7. L. Cueto-Felgueroso and R. Juanes, Phys. Rev. Lett., 2012, 108, 144502.
8. H. G. Lee and J. Kim, Phys. A, 2015, 423, 33–50.
9. S. Y. Yeh and C. W. Lan, Langmuir, 2015, 31, 9348–9355.
10. Y. Wu, F. Wang, W. Huang, M. Selzer and B. Nestler, Phys. Rev. Fluids, 2022, 7, 054004.
11. F. Boyer and C. Lapuerta, Math. Modell. Numer. Anal., 2006, 40, 653–687.
12. D. M. Anderson, G. B. McFadden and A. A. Wheeler, Annu. Rev. Fluid Mech., 1998, 30, 139–165.
13. C. Liu and J. Shen, Phys. D, 2003, 179, 211–228.
14. F. Campeloa and A. Hernandez-Machado, Eur. Phys. J. E: Soft Matter Biol. Phys., 2006, 20, 37–45.
15. R. Teshigawara and A. Onuki, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 021603.
16. A. Badillo, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 033005.
17. H. Wu, Electron. Res. Arch., 2022, 30, 2788–2832.
18. J. Yang, Y. Li and J. Kim, J. Comput. Phys., 2023, 491, 112345.
19. Z. Wang, W. Zhang, X. Mao, K.-S. Choi and S. Li, J. Comput. Phys., 2024, 516, 113297.
20. K. A. Brakke, J. Klinowski and A. L. Mackay, Philos. Trans. R. Soc. London, Ser. A, 1996, 354, 2143–2157.
21. G. McHale, M. I. Newton and B. J. Carroll, Oil Gas Sci. Technol., 2001, 56, 47–54.
22. G. McHale and M. I. Newton, Colloids Surf., A, 2002, 206, 79–86.
23. S. H. Collicott, W. G. Lindsley and D. G. Frazer, Phys. Fluids, 2006, 18, 087109.
24. X.-F. Wu, M. Yu, Z. Zhou, A. Bedarkar and Y. Zhao, Appl. Surf. Sci., 2014, 294, 49–57.
25. H. B. Eral, J. de Ruiter, R. de Ruiter, J. M. Oh, C. Semprebon, M. Brinkmann and F. Mugele, Soft Matter, 2011, 7, 5138–5143.
26. T.-H. Chou, S.-J. Hong, Y.-E. Liang, H.-K. Tsao and Y.-J. Sheng, Langmuir, 2011, 27, 3685–3692.
27. Y.-E. Liang, C.-J. Wu, H.-K. Tsao and Y.-J. Sheng, J. Phys. Chem. C, 2015, 119, 25880–25886.
28. C. Lv and S. Hardt, Soft Matter, 2021, 17, 1756–1772.
29. D. Liu and J. Nocedal, Math. Program. B, 1989, 45, 503–528.
30. S. P. Lin and W. C. Liu, AIChE J., 1975, 21, 775–782.
31. P. A. Gauglitz and C. J. Radke, Chem. Eng. Sci., 1988, 43, 1457–1465.
32. B. Song, A. Bismarck, R. Tahhan and J. Springer, J. Colloid Interface Sci., 1998, 197, 68–77.
33. L. A. Slobozhanin, J. I. D. Alexander and A. I. Fedoseyev, Phys. Fluids, 1999, 11, 3668–3677.
34. X. F. Wu and Y. A. Dzenis, Acta Mech., 2006, 185, 215–225.
35. Z. Lu, T. W. Ng and Y. Yu, Int. J. Heat Mass Transfer, 2016, 93, 1132–1136.
36. B. Zheng and R. F. Ismagilov, Angew. Chem., Int. Ed., 2005, 44, 2520–2523.
37. C. W. Extrand, Langmuir, 2007, 23, 1867–1871.
38. B. Carroll, J. Colloid Interface Sci., 1976, 57, 488–495.
39. H. D. Wagner, J. Appl. Phys., 1990, 67, 1352–1355.
40. O. Wodo and B. Ganapathysubramanian, J. Comput. Phys., 2011, 230, 6037–6060.
41. M. Fernandino and C. A. Dorao, Appl. Math. Modell., 2011, 35, 797–806.
42. M. Hintermüller, M. Hinze and M. H. Tber, Optim. Methods Software, 2011, 26, 777–811.
43. D. Lee, J.-Y. Huh, D. Jeong, J. Shin, A. Yun and J. Kim, Comput. Mater. Sci., 2014, 81, 216–225.
44. X. Yang, J. Zhao, Q. Wang and J. Shen, Math. Models Methods Appl. Sci., 2017, 27, 1993–2030.
45. Q. Cheng, C. Liu and J. Shen, J. Comput. Appl. Math., 2021, 394, 113532.
46. Y. Rohanizadegan, H. Li and J. Z. Y. Chen, Soft Matter, 2024, 20, 5359–5366.
47. C. Pan, S. Feng, S. Tao, H. Zhang, Y. Zheng and H. Ye, Phys. Fluids, 2024, 36, 022116.
48. X.-F. Wu, A. Bedarkar and K. A. Vaynberg, J. Colloid Interface Sci., 2010, 341, 326–332.
49. S. Protiere, C. Duprat and H. A. Stone, Soft Matter, 2013, 9, 271–276.
50. Y.-E. Liang, H.-K. Tsao and Y.-J. Sheng, Langmuir, 2015, 31, 1704–1710.
51. Y.-E. Liang, C.-C. Chang, H.-K. Tsao and Y.-J. Sheng, Soft Matter, 2013, 9, 9867–9875.

---