We develop a two-dimensional model for the transient diffusion of gas from the cavities in ridge-type structured surfaces to a quiescent liquid suspended above them in the Cassie state to predict the location of the liquid vapor-interface (meniscus) as a function of time. The transient diffusion equation is numerically solved by a Chebyshev collocation (spectral) method coupled to the Young–Laplace equation and the ideal gas law. We capture the effects of variable meniscus curvature and, subsequently, when applicable, movement of triple contact lines. Results are presented for the evolution of the dissolved gas concentration field in the liquid and, when applicable, the time it takes for a meniscus to depin and that for longevity, i.e., the onset of the Cassie to Wenzel state transition. Two configurations are examined; viz., one where an impermeable membrane pressurizes the liquid above the ridges and one where hydrostatic pressure is considered and the top of the liquid is exposed to noncondensible gas.
Superhydrophobic surfaces (SHs) are structured surfaces that can support liquid in the unwetted (Cassie) state. They consist of a hydrophobic substrate on which liquid is supported between protruding structure elements by surface tension and (non-condensible) gas pressure, resulting in both gas and solid contact with the liquid (the vapor pressure of the liquid may also provide support). The distance between the protruding structures can vary but is often of the order of tens to hundreds of microns (see, Refs. [1–3]). A liquid is said to be in the Cassie state when the cavities between the surface elements are unwetted and in the Wenzel state when the liquid fills these cavities. For liquid droplets in the Cassie state, the reduction in total solid–liquid interfacial area results in apparent contact angles approaching 180 deg . In liquid flow in the Cassie state over SHs, this reduction in the solid–liquid interfacial area reduces frictional drag ; however, our focus in this paper is on a stationary liquid. Engineering and biological applications of SHs with stationary liquids include antifouling and respiration of underwater insects [6,7], respectively.
For a liquid column that is supported in the Cassie state, the balance of forces along with liquid–vapor interfaces (menisci), assuming they are stationary or move very slowly, is given by the Young–Laplace equation. However, applied pressure and gas diffusion can cause a liquid column that is initially in the Cassie state to transition to the wetted (Wenzel) one [3,4,8–10]. The commonly studied ridge-type SHs with closed cavities are our focus here. Upon submersion in a column of water, the gas in the cavity(s) of SHs with closed cavities is rapidly compressed, in what may be assumed an isothermal or adiabatic process, to a pressure we denote by [11,12]. The corresponding local increase in the dissolved gas concentration on the liquid side of the menisci in accordance with Henry's Law causes gas molecules to diffuse into the liquid. The diffusion of gas out of the cavity results in a decrease of the gas pressure as per the ideal gas law. Meniscus deformation can be separated into two regimes. In the first regime, the meniscus must satisfy a force balance as per the Young–Laplace equation and the contact angle (θ) increases. As gas continues to diffuse, the contact angle continues to increase and may reach its maximum value, i.e., the advancing contact angle () [8,11]. We refer to this time as a critical time . The second regime of meniscus deformation occurs after critical time. In this regime, when additional gas molecules diffuse into the liquid, the volume of gas in a cavity further decreases, and thus menisci depin and slide down the cavity walls . Finally, the menisci may reach the bottom surface (Wenzel state). The duration of the Cassie to Wenzel state transition is known as longevity and is denoted [10,12]2.
The upward component of the surface tension force during the metastable Cassie state depends on the shape of the cavity sidewalls . For example, for ridge-type structures that have straight side walls, it remains fixed during movement of the triple contact line. Conversely, for other sidewall profiles, such as cones or reentrant structures, it is dependent on the location of the triple contact line . In some experimental studies, the onset of the Wenzel state has been observed before the menisci fronts are theoretically expected to contact the bottom surface, possibly as a result of imperfections in the cavity walls or the presence of another favorable minimum energy state [8,12,14].
The mechanisms for Cassie to Wenzel state transition can be observed and validated experimentally using, for example, optical [10,14] and acoustic techniques . Lv et al.  reported the menisci shape during a pressure- and diffusion-induced transitions on micropore type structures using a confocal microscope. They measured the contact angle and displacement of a triple contact line as a function of time. Their results demonstrate the two regimes of meniscus deformation, initially showing an increase in contact angle followed by sliding into micropores upon achievement of the advancing contact angle. Similarly, Xu et al.  used a high-resolution camera to capture the gas diffusion-induced transition on a microscale trench sample. Their method is discussed in detail later. Additionally, several of the experimental studies showed that increasing the degree of gas saturation in the liquid can slow the diffusion process and thus enhance the longevity of SHs [10,12].
Early work on modeling gas diffusion-induced Cassie to Wenzel state transition on ridge-type structures was done by Emami et al. . In their model, the diffusion of dissolved gas into the liquid is proportional to the product of gas pressure in the cavity and an experimentally determined average mass transfer coefficient for air. This model has been used by others to correlate experimental data on gas diffusion-induced Cassie to Wenzel state transition by prescribing an appropriate value for the aforementioned mass transfer coefficient [3,10]. In contrast, for liquid over cylindrical microcavities, Søgaard et al.  numerically solved a one-dimensional transient gas diffusion model and successfully correlated their experimental data without the need for experimentally-determined fitting parameters. In their model, the radius of curvature of the meniscus was held constant at the advancing contact angle, ignoring meniscus deformation before critical time. The Young–Laplace equation was used to determine pressure within the cavity and this pressure was coupled to a one-dimensional diffusion equation via Henry's constant.
Previously, we developed a semi-analytical, one-dimensional solution for the evolving dissolved gas concentration field in the liquid, the critical time, and longevity on ridge-type SHs . We captured the effect of a time-dependent boundary condition on the menisci by applying Duhamel's theorem. Meniscus curvature was accounted for in the Young–Laplace equation and ideal gas law. However, in solving the transient diffusion equation, we assumed that the solid fraction of the ridges, defined as the meniscus length divided by ridge pitch, was vanishingly small and the meniscus was flat. In this numerical study, we relax the foregoing assumptions and thus refine our predictions for times required for menisci to depin and transition to the Wenzel state.
Xu et al.  experimentally capture the effects we model in this paper. They submerged a single superhydrophobic cavity whose length was much larger than its width in water and imaged the meniscus location over time. Initially, the meniscus was flat but immediately upon submersion the contact angle started to increase as gas diffused into the liquid. In the middle of the long cavity, far from the end walls, the deforming meniscus was unaffected by edge effects and the deformation appeared two-dimensional when viewed from the side. Upon reaching the advancing contact angle, the meniscus depinned from the cavity edges and proceeded to travel into the cavity with a constant radius of curvature. A complimentary model was developed using an empirical parameter to capture meniscus curvature effects. We build a more accurate model than the one implemented in Xu et al.  by fully resolving the two-dimensional gas distribution in the liquid, explicitly capturing meniscus curvature and transient, two-dimensional diffusion. Physically our model can be viewed as a two-dimensional slice of a periodic superhydrophobic surface in which the slice is taken in the middle of a long channel such as the one seen in the experiments by Xu et al. . However, we note that Xu et al.  focus on a single channel and our model is the first to our knowledge that resolves a fully periodic system of grooves and thus there is no direct experimental data available against which to compare our model. In contrast to our schematic configuration, Xu et al.  submerge their single grooves in a cylindrical container. Too, periodic systems of grooves are more common than single grooves in practice due to larger surface coverage.
2 Mathematical Model
where is the residual mass of gas in the cavity per unit depth, is the molecular weight of the gas, is the universal gas constant and T is the absolute temperature of the gas.
The meniscus has been assumed flat at the beginning of the isothermal compression process. These are reasonable assumptions when pouring water onto a uniformly-covered textured surface as air will be trapped in the trench at atmospheric pressure and initially the meniscus supports a negligible pressure difference. We note that the meniscus immediately depins if the hydrostatic pressure reduces volume in the cavity below the volume needed to reach the advancing contact angle, i.e., and , if .
A dimensionless analysis is forgone due to the appearance of numerous dimensionless numbers that appear when scaling the various constraints in this problem.
3 Numerical Method
We apply a spectral (Chebyshev collocation) method to numerically solve the transient, two-dimensional diffusion equation, Eq. (8), for the partial density of the dissolved gas in the liquid. This preserves discretization accuracy near curved boundaries. However, once transformed, the diffusion equation and boundary conditions are more complex.
3.2 Domain Decomposition and Mapping.
As per Fig. 2, we solve for the dissolved gas concentration field in domain, Ω by decomposing it into an arbitrary number of subdomains, Ωq, for for and for as shown in Fig. 2. The columns of subdomains to the left and right of the triple contact point at for and for are separated by a (vertical) line emanating from it. There are Q rows and 2 columns of subdomains up until , when subdomain is added to the bottom of the left column. Straight boundaries of subdomains parallel to the x axis correspond to fixed values of y denoted by yq, where , except at the top of the domain, where the value of y decreases with time and corresponds to . Care must be taken in choosing the height of the topmost subdomains to ensure the moving boundary there can be captured in that subdomain for the entirety of the solution. Heights of the subdomains are chosen so that those closest to the meniscus have aspect ratios near unity. Nodes are aligned on subdomain interfaces. Before critical time, the left, bottom subdomain, , is bounded by a deforming boundary, , and a fixed one, y = y1. At critical time, the meniscus depins and the triple contact line moves down the impermeable ridge; therefore, we separate into two subdomains as per Fig. 2(b). Then, subdomain becomes rectangular and bounded from below by y = 0 and subdomain contains the rest of the liquid previously in subdomain . In general, the subdomains fall into 7 categories based on their applicable boundary conditions and geometry and this is discussed in detail in Appendix A.
where , and are the positions of the left, right, lower and upper boundaries, respectively, on subdomain q . Here only depends upon x and t for subdomain 1 when and subdomain when and only depends on t for subdomains Q and 2Q.
3.3 Time Discretization.
3.4 Transformed Diffusion Equation.
for the case of three-point time discretization. Here, Iq is an identity matrix and yt has the same dimensions. We note that, in certain subdomains, , ηxx, , and change with time; therefore, at each time-step, we update them using the subdomain-level equations as described in Sec. 3.6.
3.5 Assembly of Linear System of Equations.
In discussing this system of equations we refer to the matrix operator as , the full solution vector as , and the forcing vector as . Examples of how boundary nodes are accommodated are as follows. First, to prescribe an (external) homogeneous Neumann condition in x on the (noncorner) Chebyshev node in the column vector , we set the row of equal to the right-side of Eq. (47) and the row of equal to 0. Secondly, for the node shared by subdomains 1 and 2 along the vertical line of symmetry on the left-side of the domain (see Fig. 2), we modify rows in , and and utilize the relevant rows in and to impose continuity of temperature and diffusive flux in the y direction. (The choice of the interfacial boundary condition, i.e., flux or concentration matching, assigned to a given subdomain is arbitrary.) We do not impose the homogeneous Neumann condition in x at this node. Finally, for a corner node shared by 4 subdomains, we match 3 concentrations and one diffusive flux. Additional information on the subject can be found in the texts by Canuto et al.  and Trefethen .
3.6 Sumdomain Types.
The subdomains fall into 7 groups as follows.
The left, bottom subdomain when in case study I and in case study II bounded by and , i.e., as per Fig. 2(a). After critical time, it subdivides into of group (B) and of group (G) as per Fig. 2(b).
Middle subdomains on , i.e., from to and, additionally, for .
The top subdomain on , i.e., ΩQ.
The bottom subdomain on a < x < d, i.e., .
Middle subdomains on a < x < d, from to .
The top subdomain on a < x < d, i.e., .
The bottom subdomain on after critical time, . Before critical time, the portion of the solution vector corresponding to is set to zero because it's only defined after the triple contact line depins.
In order to formulate the necessary system of equations, we describe boundary conditions and geometric maps, and , for each subdomain group. These are discussed in detail in Appendix A. Additionally, as part of the numerical method, we performed singularity removals on the singular points in the problem. The triple contact point is singular before and after critical time and the point at which the meniscus intersects the cavity wall is singular. The singularity removals were done in order to ensure correct capture of the singularity behavior, as well as increase the speed and accuracy of the solution by removing the need for large mesh densities in the vicinity of the singular points. The singularity removals are discussed in detail in Appendix B.
3.7.1 Solution Up to Critical Time.
where ϵ is user-defined.
The algorithm is:
Determine the and matrices using solution of the previous iteration as an initial condition.
Numerically solve for and as per Eq. (B9).
Compute as per Eqs. (54).
Update the geometric parameters, , and . Compute time derivatives using a 3-point formula of the form of Eq. (45).
If the convergence criterion given by Eq. (55) is met, then repeat steps 1 through 4 for the next time step . If the convergence criterion is not met, proceed to the next iteration (i + 1) for the current time step (k).
When critical time has been reached, we fix the gas pressure and allow the meniscus to depin.
3.7.2 After Critical Time.
For subsequent time steps, we employ the same approach as described in the previous section, except we compute instead of .
4 Results and Discussion
We implemented the model in MATLAB®. Assuming that the liquid and gas are pure water and nitrogen, respectively, and that liquid rests on a silicon substrate coated with a flouropolymer coating, we set the thermophysical properties to , ,  and  at ambient conditions of and . We note in practice that air would likely be the lgas in the cavity; however, since air is primarily nitrogen we assume the gas is pure nitrogen to restrict attention to a single species, simplifying the model. Additionally, is the advancing contact angle of water on silicon treated with a flouroploymer . Unless noted, we set the geometric parameters to be and with the liquid initially degassed. To generate the computational domain for case study , we partitioned the liquid domain into 4 subdomains (before critical time), i.e., Q = 2, where the interfaces yq are located on . This allocated most of the nodes to the region near the triple point where gradients are largest. For case study II, we partitioned the liquid domain into 6 subdomains (before critical time), i.e., Q = 3, where the interfaces yq are located on . This allocated nodes close to both the triple point and to the top-most meniscus. Before critical time the time stepping was handled as follows. A time-step of was initially used before critical time. However, as the angle of the deforming meniscus neared the advancing contact angle it was changed to to more accurately capture . A time-step of after critical time was needed to accurately capture . A mesh dependence study was used to calculate the necessary spatial grid size. The total number of nodes was doubled until varied by less than 0.1 percent. This resulted in L = 24 and N = 20, a total of 1920 nodes, for case study I and L = 36 and N = 18, a total 3888 nodes, for case study . For both cases varied by less than 0.25 percent between the final two iterations. The difference in grid construction between case study I and case study was in order to better capture the sharp gradients at the top of the domain that occur in case study .
4.1 Case Study I.
To examine case study we ran a series of computations that spanned the parameter space of solid fraction (i.e., varied between 0 and 1) over a large range of values of the initial height of the liquid column (h0). Figure 3 plots critical time versus solid fraction when h0, defined as the initial liquid column height, is fixed at 800 μm and the half-width of the gas cavity is a = 5 μm. Clearly, increasing solid fraction dramatically increases gas diffusion out of the cavity because of the increase in the width of the domain (2d) leading to lateral diffusion. Of note, critical time is almost entirely independent of chamber height, h0, for these specific parameters because the time it takes for gas to diffuse through the entire domain is significantly larger than . Also plotted is critical time computed from the quasi-one-dimensional model in our earlier publication . This model is exact in the limit as and and is far easier to solve, but begins to veer from the exact solution as these parameters change in value. As can be seen in Fig. 4 treating the meniscus as flat as opposed to having an advancing contact angle of leads to an underestimation of the critical time by 6% as .
which is a modification to a result found in our earlier publication, i.e., Kadoko et al. . The computed data is again compared to that of Kadoko et al.  and again shows underestimation with final time solutions veering from the 2D solution by about 9%.
Figure 5 plots the volume of gas in the cavity versus time for different solid fractions. The inset shows the results for the critical time and slightly thereafter. Since the cavity width is the same, the effect of different solid fractions can easily be explored. As gas diffuses from the cavity into the liquid, the meniscus deforms and eventually depins, reducing the volume of gas in the cavity. This happens significantly faster for domains with the same cavity width, but larger solid fraction because the total volume for diffusion is much larger, approaching that of a half-space for small . The discontinuity visible in the inset highlights the transition between the two types of meniscus movement: (1) meniscus deformation and (2) sliding of the contact line. It exists because after critical time the pressure in the cavity remains constant and s changes, whereas before critical time the pressure varies to accommodate the deforming meniscus. is the volume of gas in the cavity when the lowest part of the meniscus makes contact with the cavity bottom.
Figure 6 shows contour plots of the dissolved gas concentration field in the lower region of the liquid domain, , with . Figure 6(a) shows ρ at , soon after the diffusion process has started, such that the concentration gradient is relatively large and the meniscus has only begun to deform. Figure 6(b) depicts ρ at critical time, right before the meniscus begins to slide. Figures 6(c) and 6(d) show ρ at and , respectively. At this point the concentration gradient has decreased significantly and diffusion rates are low. This helps illuminate why is much larger than .
4.2 Case Study II.
Figure 7 plots longevity as a function of cavity depth for a characteristic initial liquid height, and cavity half-width of . The first thing to note is that a larger initial column height of in case study leads to less longevity, , than in case study for and the same cavity depth, , as per Fig. 4. This is because in case study hydrostatic pressure is taken into account. If the effect of hydrostatic pressure is significant the pressure in the cavity will be higher for the same atmospheric pressure, increasing the diffusive driving force. Secondly, Fig. 7 also shows that increasing cavity depth yields larger critical times and longevity because more gas must diffuse out of the cavity to reach the Wenzel state. This is potential valuable insight for design of these surfaces. Plotted alongside the results is the quasi-one-dimensional solution for case study from Kadoko et al. . In this case, it fulfills the solution requirement for , demonstrating much the same shape as our results.
Figure 8 plots longevity for varying initial gas concentration for a characteristic initial liquid height, . Increasing the initial gas concentration reduces the driving force for diffusion, increasing both critical time and longevity. Longevity grows rapidly with increasing initial concentration. This highlights that initial gas concentration, along with cavity depth and total chamber height, can be tuned to achieve the desired performance of a surface. Again to span the whole parameter space, the results from Kadoko et al.  for are included.
We present a two-dimensional numerical analysis of the gas diffusion-induced Cassie to Wenzel state transition in a quiescent liquid resting on parallel-ridge type SHs. This analysis is representative of the meniscus behavior in the middle of a long three-dimensional meniscus such as experimentally studied by Xu et al. . We demonstrate that the algorithm captures the deforming and moving boundaries of the liquid domain seen in experiments and the time-dependent boundary condition on the meniscus in parallel-ridge type SHs. Too, it computes the dissolved gas concentration field in the liquid. The algorithm can be extended to hierarchical reentrant SHPo surfaces where sidewalls of the structures have nonstraight profiles. Future work needed is resolution of the fully three dimensional problem, capturing wall effects near ends of the long channels, and treating the gas in the cavity as air is such that two species equations need to be resolved. Finally, the lack of experimental data to which to compare our results shows a need for experimental work capturing diffusion effects on periodic surfaces.
This material is based upon work supported in part by the National Science Foundation under Grant No. 1402783 and supported in part by a Google Award in Faculty Research to Professor Marc Hodes. We acknowledge useful discussions with Dr. Simon Game and Dr. Georgios Karamanis.
Google Award in Faculty Research.
National Science Foundation under (Grant No. 1402783; Funder ID: 10.13039/100000001).
For hierarchical SHs that are composed of micro-structures fabricated on top of a micro-structured substrate, only the nano-structures are wetted when the menisci reach the bottom surface (partial Wenzel state) . A complete Wenzel state is achieved when the nano-structures are wetted.
Although we ignore the effects of changes in meniscus shape and position on hydrostatic pressure, we do account for them in the rest of our formulation as per below as they are essential to, e.g., calculate the volume of the gas phase.
We have dropped the subscript q from L and N as they are the same in all subdomains.
Appendix A: Subdomain Types
A.1 Subdomain Type (A).
A.2 Subdomain Type (B).
A.3 Subdomain Type (C).
A.4 Subdomain Type (D).
A.5 Subdomain Type (E).
A.6 Subdomain Type (F).
A.7 Subdomain Type (G).
The partial derivative in y with respect to t is given Eq. (A10), where , because radius of the meniscus is constant after critical time.
Appendix B: Singularity Removal
B.1 Singularity Removal.
A singularity in diffusive flux exists at the triple contact point (x = a, y = 0) before critical time. After it, as the triple contact point advances down the ridge, one remains at this location and a second one follows the triple contact point. We remove them from our numerical solution (and, subsequently, reincorporate them into our final one) in the manner of Game et al.  to reduce the required mesh density as per the following steps. First, local analyses in the vicinity of the singular points are performed to resolve the forms of the singularities up to undetermined constants coefficients, the singularity strengths. Next, a column is added to that contains the contribution(s) to the diffusion equation and (external) boundary conditions from the solution(s) to the local (singular) problem. This requires that the unknown singularity strength(s) be appended to the solution vector. Finally, a row(s) is added to to bound the derivative of the numerical solution at the singularity(s). The resulting matrix system is one column (and one or two rows) larger than prior to singularity removal, a negligible increase in computational load. Below we discuss the specific implementation of each singularity.
B.1.1 Before Critical Time.
for the usual three-point formula, the two-point formula and the three-point formula with nonuniform time-step size, respectively. This is necessary, because resolves only the Laplacian terms in the diffusion equation, Eq. (52), and the penultimate and ultimate terms on its left side must thus be added.
where i indicates the i-th index and y(i) is the y location of the i-th node. Thus, the solution, satisfies the external boundary conditions.
We remove the singularity from the numerical solution, via the a row vector in Eq. (B9). It applies the differential operator at the (singular) triple contact point in subdomain Q + 1. The end result is a new system of equations. The element 0 in the right-side vector restricts and therefore removes the singularity from . The full solution is obtained by solving this new system of equations, constructing both the singular solution and the numerical solution from the solution vector (which includes the singularity strength) and adding them together. Additionally, all integrals and derivatives done in the analysis are split into a singular and nonsingular parts and then summed to ensure accuracy.
B.1.2 After Critical Time.
The local problems after critical time, denoted by the subscripts 1 and 2 at and at the triple contact point, respectively, again cast in polar () coordinates to facilitate their solution, are shown in Fig. 10. Again, locally, we need only resolve Laplace's equation and decompose our solution into numerically resolved and (singular) analytical components.
where D is the unknown singularity strength.
The role of the column vectors and are analogous to that of in Eq. (B9) before critical time. Moreover, the row vector causes the operator to act on the element of corresponding to the bottom left corner of subdomain Q + 1 and the row vector causes the operator to act on the element of corresponding to the triple contact spot in subdomain . Both these partial derivatives are set equal to 0. The result is a system of equations. The solution process is analogous to that before critical time.
- Α =
coefficient for singularity removal
vector of dissolved gas concentration
- δt =
- ϵ =
- η =
map to vertical of unit square
- Ω =
- ρ =
- σ =
surface tension along meniscus ()
- θ =
contact angle of penetration into cavity, measured from the horizontal
- ξ =
map to horizontal of unit square
matrix for coupling subdomains i and j
- a =
half-width of the cavity (m)
- D =
binary molecular diffusion coefficient () or differentiation matrix
- d =
half of the ridge period (m)
- g =
gravitational constant ()
- H =
Henry's constant ()
- h =
transient height of the liquid domain (m)
- h0 =
initial height of the liquid domain (m)
mass per unit depth of gas in cavity ()
molecular weight of gas ()
- p =
- Q =
number of subdomains in vertical direction
- R =
radius of curvature of meniscus (m)
universal gas constant ()
- s =
coordinate along arc spanning meniscus (m)
- s =
transient location (in y) of the triple contact point (m)
- s0 =
depth of the cavity (m)
- T =
- t =
- v =
contributions from singularity removal to global matrix
volume per unit depth of gas in cavity ()
- ym =
meniscus location (m)
advancing (contact angle)
differential operator on the triple contact point
gas in the cavity
referring to global matrix and vectors
referring to numerical nonsingular solution
referring to quasi-analytical singular solution