Letter pubs.acs.org/Langmuir
Hybrid Energy-Minimization Simulation of Equilibrium Droplet Shapes on Hydrophilic/Hydrophobic Patterned Surfaces Hiroyuki Matsui,* Yuki Noda, and Tatsuo Hasegawa National Institute of Advanced Industrial Science and Technology (AIST), Tsukuba 305-8562, Japan S Supporting Information *
ABSTRACT: We have developed an efficient algorithm for simulating the equilibrium shape of a microdroplet placed on a flat substrate that has a fine, discontinuous, and arbitrarily shaped hydrophilic/hydrophobic patterned surface. The method uses a hybrid energy-minimization technique that combines the direct search method to determine the droplet shape around solid/liquid contact lines with the gradient descent method for the other parts of the droplet surface. The method provides high-convergence at a low computational cost with sufficient mesh resolution, providing a useful tool for the optimal design of printed electronic devices.
■
INTRODUCTION Wetting phenomena of microliquids on solid surfaces has attracted considerable attention in recent times, because understanding these phenomena forms a fundamental basis for the emerging field of “printed electronics.”1−6 Printing technology is used for the patterned deposition of functional inks that contain dissolved or dispersed semiconducting or metallic materials, which enables the construction of active or passive electronic components. Particularly, the technology is expected to realize the low-cost production of large area and flexible electronic products under ambient conditions. One of the crucial challenges in printed electronics is to achieve printed thin-film patterning on flat substrate surfaces with sufficient resolution, although this is not feasible with conventional processes such as inkjet printing. It has been suggested that a predefined hydrophilic/hydrophobic patterned surface is promising for finer thin-film patterning because the hydrophilic area should primarily be wetted by the inks.7 However, it is also true that the shape of a deposited droplet is not always coincident with the predefined pattern, depending on the wettability conditions as well as the fluidic nature of the microdroplet. Therefore, it is necessary to simulate and predict the shape of droplets on substrate surfaces for optimal circuit design of printed electronic devices. The most popular method of droplet shape simulation is to solve the Navier−Stokes (NS) equation and continuity equation numerically by one method from among the volume of fluid (VOF),8 level set,9 and phase field methods.10 Another recently developed technique is the lattice Boltzmann method, which solves the discrete Boltzmann equation instead of the NS equation.11−13 Although these approaches can predict a temporal variation in the droplet shape, there is a strict limitation in applying them to realistic three-dimensional (3D) problems in terms of accuracy and calculation cost. Particularly, approximations of the liquid−vapor interface in these methods have intrinsic risks of yielding incorrect results. The calculation © 2012 American Chemical Society
based on the methods also needs at least several hours, even for quite rough calculations, and needs much more time for obtaining reliable results by using software such as COMSOL.14 In contrast, the energy-minimization (EM) approach should be less time-consuming and afford reliable results in cases wherein the focus is only on the equilibrium shape of the droplet on the patterned substrate.15−25 Some of the work on the EM calculation for free surfaces has been carried out using a software such as Surface Evolver.15,16 An elementary step in EM is to evolve the droplet shape so that the volume remains constant while the total surface energy at the liquid−vapor, liquid−solid, and solid−vapor interfaces is reduced. After the iteration of the small evolutions, we obtain the minimal-energy state at equilibrium. Since a droplet surface is a system that has a large degree of freedom, the calculation speed and accuracy strongly depend on the algorithm that defines how to direct the evolution of the droplet shape. Although the gradient descent and conjugate gradient methods are simple and well-established algorithms as utilized in Surface Evolver,15,16 they are not applicable to systems with discontinuous energy functions such as a hydrophilic/hydrophobic boundary. This problem has been avoided, thus far, by smoothing the wettability patterns19 or by adding the restriction that the contact lines cannot intercept discontinuous boundaries.24 Under the mentioned restriction, the initial state of the droplet must be sufficiently close to the unknown final state, which significantly limits the scope and range for applying droplet shape simulations. In the present paper, we report a hybrid EM algorithm, which provides high convergence with a short calculation time for simulating the equilibrium shape of microdroplets placed on fine and discontinuous hydrophilic/hydrophobic patterned surfaces. We apply the direct search method to the droplet Received: September 15, 2012 Revised: October 10, 2012 Published: October 23, 2012 15450
dx.doi.org/10.1021/la303717n | Langmuir 2012, 28, 15450−15453
Langmuir
Letter
the total number of vertices), which is effective in minimizing the calculation cost of the direct search method. When the droplet is significantly deformed during the EM procedure, the triangles may become elongated and the density of mesh becomes inhomogeneous. Hence, we performed a mesh reconstruction process for every 10 updates. The reconstruction consists of two steps. The first step is to merge two neighboring vertices if the distance between them is less than a specific length, L. In this step, one or two triangles are removed. The second step is to place a new vertex at the midpoint between two vertices if the distance is larger than 2L. In this step, one or two triangles are added. After the mesh reconstruction process, the sides of all triangles have similar lengths, between L and 2L. To improve calculation efficiency further, the initial droplet surface is represented by a relatively small number of triangles, typically a few hundred, while the number of vertices is about half that number. When the first minimization converges, we increase the density of the triangles by dividing each triangle into four small triangles (this process is the same as “refinement” in Surface Evolver15,16). These refinements and subsequent EM updates are conducted several times, and the final number of triangles is typically as large as 104−105. The program was coded using Java SE 6 and executed on a commercial laptop [CPU: Intel Core i7−2640M (2.80 GHz); Memory: 8 GB; OS: Windows 7 (64-bit)].
shape around the solid/liquid contact lines, and we apply the gradient descent method to the other parts of the droplet surface, whose variation causes only a smooth change in the surface energy. One of the advantages of this method is that the initial shape of the droplets can be defined independently of the substrate patterns because the contact lines are allowed to cross the boundaries of the discontinuous patterns, which is in sharp contrast to the conventional algorithm that places several restrictions on the initial droplet shape as well as on the evolution of droplet on such problems. The entire calculation is completed within a few seconds in many cases by simply using a laptop with a standard configuration.
■
SIMULATION METHODS
A liquid−air interface is implemented as a mesh composed of triangular elements.15−25 Each vertex is usually shared by several triangles. The initial shape of a droplet is simply assumed as a partial sphere that contacts the substrate surface at an arbitrary contact angle. The substrate is assumed to be flat, and the area is divided into two or more domains (specified by k), each of which has a different wettability with contact angles, θk. By summing the liquid−vapor, liquid−solid, and solid−vapor interface energies, the total energy of the droplet is calculated as
■
RESULTS AND DISCUSSION Figure 1a shows the decrease in the total energy in a 50 pL droplet placed on a homogeneous substrate with an equilibrium
M
E = γLV(ALV −
∑ ASL , k cos θk)
(1)
k=1
where γLV is the surface tension of the liquid; ALV is the area of the liquid−vapor interface, and ASL,k is the area of the wetted region in the kth domain of the substrate surface.19,21,22,24,25 We can neglect the effect of gravitation because the volume of our droplet is less than 1 μL. In order to identify the equilibrium shape of the droplet, we attempt to create small changes in the droplet shape by shifting vertices to decrease the total energy. In this process, we shift the vertices only along the normal vector n⃗i of the droplet surface because the movements of vertices along the surface do not change the shape effectively. A shift of the ith vertex is represented as Δr⃗i = λin⃗i, where λi is the length of the shift. The process is iterated until the change of energy becomes sufficiently small, at which point the shape is regarded as the equilibrium shape of the droplet. The essence of the EM algorithm is to determine a good set of {λi} in each step. Since the energy function E(λ1,···, λN) is not smooth when some of the vertices are on the boundary between the surfaces with different wettabilities, the gradient methods would not result in a good convergence. To avoid this issue, we use a combination of the gradient descent and direct search methods. We divide all the vertices into two groups: those which form solid/liquid contact lines (Group A) and those which do not contact the substrate (Group B), and we update the two groups alternately. First, the vertices in Group B are updated by the gradient descent method, while the vertices in Group A are fixed. In order to maintain a constant volume of the droplet, the gradient vector (∂E/∂λi) is projected onto the subspace perpendicular to the gradient vector of the volume (∂V/∂λi):
⎧ ⎡ ⎪ ∂E ∂V ⎢ ⎛⎜ ∂E ∂V ⎞⎟ λi = C ⎨ − ∑ ∂λi ⎢⎣ j ⎜⎝ ∂λj ∂λj ⎟⎠ ⎪ ∂λi ⎩
⎛ ∂V ⎞2 ⎤⎫ ⎪ ∑ ⎜⎜ ⎟⎟ ⎥⎥⎬ ∂ λ j ⎠ ⎦⎪ j ⎝ ⎭
Figure 1. Variation of total surface energy of a droplet (a) on a homogeneous substrate and (b) on a predefined hydrophilic surface (50 × 200 μm; contact angle of 10° inside and 90° outside) as a function of the number of steps for the energy-minimization calculation. Solid and dashed lines represent the results obtained by the hybrid and gradient descent methods, respectively. Initial droplet shape and the shapes after minimization (500 steps) are shown in the inset.
(2)
Here, the summation is taken for the vertices in Group B. The scaling constant, C, is determined such that the decrease in energy is as large as possible. Next, the vertices in Group A are updated in a one-by-one direct search method. For each vertex in Group A, we select a counter vertex in Group B at random and move the two vertices together so that the volume does not change. The direction of movements is again normal to the surface, and the distance is determined so that the energy decreases as much as possible. The direct search updates are performed for all the vertices in Group A. It should be noted that the number of vertices in Group A is only on the order of √N (here, N is
contact angle of 45°. The initial shape of the droplet was assumed as a hemisphere having a contact angle of 90°. During EM calculations, the energy decreased exponentially and converged to the same value in both the gradient descent and the hybrid methods. The energy residual after the 500 steps was estimated at
