Multifaceted design optimization for superomniphobic surfaces

See allHide authors and affiliations

Science Advances  21 Jun 2019:
Vol. 5, no. 6, eaav7328
DOI: 10.1126/sciadv.aav7328


Superomniphobic textures are at the frontier of surface design for vast arrays of applications. Despite recent substantial advances in fabrication methods for reentrant and doubly reentrant microstructures, design optimization remains a major challenge. We overcome this in two stages. First, we develop readily generalizable computational methods to systematically survey three key wetting properties: contact angle hysteresis, critical pressure, and minimum energy wetting barrier. For each, we uncover multiple competing mechanisms, leading to the development of quantitative models and correction of inaccurate assumptions in prevailing models. Second, we combine these analyses simultaneously, demonstrating the power of this strategy by optimizing structures that are designed to overcome challenges in two emerging applications: membrane distillation and digital microfluidics. As the wetting properties are antagonistically coupled, this multifaceted approach is essential for optimal design. When large surveys are impractical, we show that genetic algorithms enable efficient optimization, offering speedups of up to 10,000 times.


Superomniphobic surfaces show physical micro- or nanotexturing that enable even low–surface tension liquids to remain suspended atop a vapor-filled surface structure. This vapor-suspended state is prized for its liquid-shedding abilities, enabling high droplet mobility and low viscous drag (1). These surfaces have considerable potential to be transformative across a broad array of applications. These range from tackling current global-scale crises, via sustainable technologies for water purification (2, 3) and antimicrobial surfaces in biomedicine (4, 5), through everyday applications such as antifingerprint coatings (6) and packaging designed to reduce food waste (7), to digital microfluidics as a versatile biological and chemical technology (8).

Two promising textures aimed to enable these technologies are the reentrant and doubly reentrant geometries. Naturally occurring examples of these structures have been shown to imbue the cuticle of the springtail arthropod (Collembola) with superolephobic properties even for highly wetting, pressurized liquids, while further exhibiting abrasion-resistance and antimicrobial abilities (912). Recent breakthroughs in microfabrication techniques have also allowed these reentrant and doubly reentrant structures, as well as more complex textures, to be produced with micrometer-scale resolution, including the use of three-dimensional (3D) printing technology, fluidization of polymer micropillars, and lithographic methods (1316).

Despite these highly versatile techniques, a large obstacle still exists to widespread development: It is not known how to design the surface structures to enable optimal performance in real-world applications. Successful superomniphobic designs must exhibit three key wetting properties: (i) a low contact angle hysteresis (CAH) to maximize liquid mobility (17), (ii) a high critical pressure, the maximum sustainable pressure at which the superoloephobic state is stable (18), and (iii) a high energetic barrier to failure, in which liquid infiltrates the surface texture and the high liquid mobility is lost (19).

The complex surface designs mean that both computational and experimental studies are expensive and time consuming to perform, and so are largely restricted in scope to considering only single wetting properties—never all three. This is highly problematic as the structural parameters of the design couple each wetting property, often antagonistically. For example, two effective ways to increase the critical pressure are to decrease the pillar-pillar separation and decrease the system scale (17). However, decreasing the separation results in an increase in the CAH (20), whereas decreasing the system scale decreases the energetic barrier to the wetting transition due to the liquid-vapor interfacial area decreasing in squared proportion. Although this threefold perspective has been introduced and advocated before [see, for example, (17)], the true lowest-energy failure mechanisms have never been incorporated.

This work overcomes the optimization challenges for superomniphobic wetting property design. We begin by developing computational strategies to systematically survey the effect of the structural parameters on the CAH, critical pressure, and minimum energy barrier individually. These methods are highly general and can be applied to any conceivable surface design. This leads to the discovery of previously unidentified mechanisms for the receding contact line and failure of the suspended state, as well as the development and validation of quantitative analytical models. We correct a number of inaccurate assumptions in prevailing models. In particular, we highlight the development of a capillary bridge model to replace the grossly inaccurate prediction for the optimal texture height in the critical pressure study.

To illustrate the importance of multifaceted optimization, we then consider two relevant example applications: water purification via membrane distillation and droplet-based digital microfluidics (Results and Discussions). Membrane distillation shows considerable potential as a sustainable, low-energy water purification technology, capable of extracting potable drinking water from highly contaminated water sources [see (2, 3) for recent reviews]. A major challenge, however, is that oils readily foul the membranes, leading to a breakdown in device performance. Meanwhile, digital microfluidics is anticipated to enable reusable, reconfigurable, and material-efficient lab-on-chip devices (8). In this technology, the major challenge is that commonly used but low–surface tension solvents pin strongly to the surface, leading to drop immobilization and device failure (21). We find here that the doubly reentrant surface geometry is ideally situated to meet these challenges, as it is robust to pressure even for highly wetting or surfactant-contaminated liquids.

In such complex surface design featuring many antagonistically coupled wetting properties, we recognize that it is not always desirable to perform large-scale wetting property surveys. Thus, we develop a genetic algorithm to perform the simultaneous optimization with high efficiency, offering a speedup of up to 10,000 times. This versatile approach is highly complementary to recent innovations in complex surface microfabrication techniques (1316), such that, together, we anticipate a transformative approach to surface design.


Contact angle hysteresis

We begin by simulating the liquid-vapor interface advancing and receding along a single row of surface structures (setup detailed in the Supplementary Materials), to obtain the macroscopic advancing and receding contact angles θa and θr, respectively, and the CAH (CAH = θa − θr). These simulated structures are shown in Fig. 1, in which structures of variable dimensions are arranged in a square array. Throughout, all dimensions shown in Fig. 1 are reported relative to the system size B and indicated with a subscript “r”. For example, the reduced cap width is Wr = W/B. Unless otherwise stated, B = 60 lattice spacings, with the cap thickness tr and lip width lr remaining fixed at 0.05. For the reentrant geometries, the lip depth is Lr = 0. Throughout, the microscopic contact angle θ° = 60° is used as a representative contact angle for organic solvents wetting fluorinated surfaces [see, for example, (22)]. We also investigate the nonwetting case presented in the Supplementary Materials.

Fig. 1 Simulation surface configuration.

Illustration of the 3D simulation repeat unit (left), with 2D cross section showing labeled structural parameters (right).

We find that the CAH depends only on the area fraction Fr of the cap [Fr = (W/B)2], and the total cap height Dr, as shown in Fig. 2A (i). Separate advancing and receding plots are shown in the Supplementary Materials, alongside available comparison with previous experimental measurements. As the liquid-vapor interface never impinges under the cap, the hysteresis is identical for both reentrant and doubly reentrant geometries.

Fig. 2 Quantification and mechanisms leading to the CAH for reentrant and doubly reentrant geometries at zero applied pressure.

(A) (i) CAH dependence on both the area fraction Fr and total cap height Dr. Symbols indicate the depinning mechanism upon receding, with purple diamonds indicating a hybrid mechanism. (ii and iii) Comparison of the bridge-, edge-, and lip-depinning receding models (solid lines, color-coded) against the simulated θr (data points); examples shown with varying Fr at fixed Dr = 0.05 and 0.35. The ±1° error bars in the simulation data are too small to be seen. (B) 3D visualization of the advancing liquid-vapor interface (shown in blue); the advancing direction is indicated by a black arrow. Black and red lines indicate the center and edge 2D cross sections that are also presented (right). (C) (i to iv) Visualizations of the major four receding mechanisms. The receding direction is indicated by black arrows.

Across the simulated parameter range, the same advancing mechanism is observed, as illustrated in Fig. 2B. In this, the advancing occurs when the angle of the approximatively planar liquid-vapor interface results in the liquid contacting the cap of the neighboring structure.

In contrast, the receding mechanism exhibits considerable variation across the parameter range and is therefore primarily responsible for the large variation in CAH observed in Fig. 2A (i). Overall, four dominant receding mechanisms are observed: bridge, edge, lip depinning, and a fourth nondepinning mechanism. Characteristic examples of these are shown in Fig. 2C (i to iv), respectively. The operative regions of each mechanism are labeled in Fig. 2A (i), indicated by gray squares, red circles, blue triangles, and black pentagons, respectively. The hybrid depinning mechanisms (purple diamonds) indicate the regions in which the dominant mechanisms smoothly interpolate. We now describe and model each of these receding mechanisms in turn.

For the lowest area fractions, at the point of receding, the three-phase contact line is pinned to the outermost pillar around the top perimeter of the cap, as shown in Fig. 2C (i). A capillary bridge is strained between the cap and bulk liquid, such that receding occurs at the point of bridge depinning. This has been observed experimentally, (20, 23, 24), and we are now able to quantitatively test the receding model proposed by Butt et al. (20, 24). In this model, the receding liquid-vapor interface strains the capillary bridge parallel to the receding direction, thus tilting the bridge from the normal (detailed in the Supplementary Materials). The two suggested consequences of this are that the direction of the pinning force is tilted from the normal by π/2 − θr/2 and that this force depends on the average contact angle θ¯=θ°+π/2θr/2. However, to yield an accurate model for use with wetting liquids as in this work, we stress that the appropriate average contact angle to use is θ¯=θ°max+π/2θr/2, where θ°max=max(θ°,π/2). This is because for wetting liquids, the maximum pinning force is achieved when the contact angle reaches 90°. With this modification, at the point of receding this microscopic pinning force balances the macroscopic force required to move the contact line, yieldingtan(θr2)=24Wrα(1)

In the simplest model, it is assumed that at the point of depinning, the three-phase contact line closely follows the square cap perimeter of length 4Wr. However, to reflect the actual contact line morphology, the shape parameter α is introduced. α is equal to 1 if the contact line is perfectly square and π/4 for a circle. As the shape parameter cannot be predicted a priori, it is treated as a fitting parameter. An example of this is shown in Fig. 2A (ii, gray curve and square points) for Dr = 0.05. Here, α = 0.861, reflecting the contact line deviating from perfectly square by depinning at the cap corners. This yields an average agreement between the simulation and model of 0.4° (average agreement for all Dr tested is 0.6° and maximum is 2°).

At large Dr and Fr, a previously unidentified capillary bridge depinning mechanism is observed in which the bridge is strained between the cap edge and bulk liquid phase. This lip-depinning receding mechanism, as shown in Fig. 2C (iii), results in a substantial decrease in 𝜃r, demonstrated in Fig. 2A (iii, blue curve and triangular points). The model that we introduce is to approximate the receding interface as a capillary bridge pinned to the side of the cap, stretched parallel to the receding direction (detailed in the Supplementary Materials). As with the bridge-depinning model, however, we must account for the average contact angle around the contact line, yieldingcos(θr2)=2(Wr+Dr)α2(2)

The accuracy of this model is demonstrated in Fig. 2A (iii), at Dr = 0.35, for which α = 0.887. This yields an average deviation between the simulation and model of 1° (average deviation for all Dr tested is 2° and maximum is 6°).

This lip-depinning model also predicts the existence of systems in which a receding contact angle no longer exists. In these extreme cases of the lip-depinning mechanism, because depinning is not able to occur, a droplet caused to move across the surface would leave a trail of suspended liquid trapped between the caps. We observe this predicted nondepinning mechanism in simulations, indicated by black pentagons in Fig. 2A (i and iii). This nonreceding case is illustrated in Fig. 2C (iv).

At intermediate area fractions and lip depths, the depinning mechanism is no longer capillary bridge mediated. Instead, the edge-pinned receding mechanism is observed, as shown in Fig. 2C (ii). Here, the interface maintains approximately the same morphology as it depins laterally from the edge of the cap. Thus, we are able to analyze the energetic change of sliding the interface laterally by a small distance to obtain the angle at which receding becomes energetically favorable—the receding angle. This is derived in the Supplementary Materials. Taking into account the liquid receding from the cap top and edges and top surface of the microchannelcosθr=(Wr+2Dr)cosθ°+Wr1(3)

This represents a generalization of previous edge-depinning models (25, 26), in which by taking account of the liquid receding from the cap sides, we are now able to describe the edge-depinning mechanism accurately for wetting liquids. This is demonstrated in Fig. 2A (ii, red curve and circular points). Without any fitting parameters, the average deviation between the model and simulation results is 2° (average deviation for all Dr tested is 3° and maximum is 7°).

Critical pressure

Unlike with the CAH, the critical pressure is sensitive to whether the surface geometry is reentrant or doubly reentrant. Throughout, the critical pressure shown ΔPc is referenced with respect to the pressure γlv/B. Here, we find that ΔPc is only dependent on the area fraction Fr and the pillar height Hr. Although ΔPc is affected by the presence of a doubly reentrant lip, ΔPc does not depend on the precise lip depth Lr. The critical pressure dependencies on Fr and Hr are shown for the reentrant and doubly reentrant structures in Fig. 3A (i and ii), respectively. For both structural types, these dependencies change across the parameter space because of the presence of two different pressure-induced failure mechanisms: Base Failure and Cap Failure.

Fig. 3 Critical pressure analysis for reentrant and doubly reentrant geometries.

(A) Contour plots of ΔPc variation with Fr and Hr for reentrant (i) and doubly reentrant (ii) geometries. Data points mark the critical height at which the failure mechanism switches from Base Failure (BF) to Depinned Cap Failure (DCF) or Pinned Cap Failure (PCF), and error bars indicate the uncertainty in this height due to the diffuse interface width. Solid and dashed white lines show the critical height based on the capillary model and 2D model, respectively. (B) Model fits to ΔPc of the Cap Failure mechanisms at Hr = 0.25 for reentrant (i) and doubly reentrant (ii) geometries. (C to E) The three failure mechanisms shown in 3D, with associated diagonal cross sections. Critical pressure liquid morphologies are shown in blue, the vapor phase is shown in white, and the interface is indicated with a black solid line. Red regions show how the unstable meniscus evolves upon increasing ΔP above ΔPc. (D and E) Under-cap views, highlighting the shapes of the contact lines at the critical pressure. (F) Details of the 3D horizontal (3DD) and 3D diagonal (3DH) capillary bridge models used, showing the inner and outer circumferences (blue) against the system configuration. The 3D illustration compares the simulated liquid-vapor interface (light blue) to the horizontal capillary model (dark blue).

For a given area fraction Fr, at low pillar heights, the Base Failure mechanism is operative, as illustrated in Fig. 3C. In this, the suspended state fails because the sagging liquid-vapor interface touches the base of the system while the three-phase contact line remains pinned to the bottom of the cap. By increasing the pillar height, ΔPc is increased.

However, at Hr above the critical height Hc, the Cap Failure mechanism becomes operative. Here, at the critical pressure, the system can no longer simultaneously support the uniform mean curvature of the liquid-vapor interface and the contact line morphology. For the reentrant geometry, this results in the contact line depinning and sliding inward, as shown in Fig. 3D. For the doubly reentrant geometry with a thin lip width lr, this results in the liquid-vapor interface ballooning outward while the contact line remains pinned, as shown in Fig. 3E. In both of these cases, increasing the pillar height further now results in no change to ΔPc.

Therefore, at fixed Fr, the maximum ΔPc occurs for HrHc, in the Cap Failure region. However, it is detrimental for design performance if the height is increased above the critical height, as this mechanically weakens the structure without increasing ΔPc (27, 28). The optimum pillar height is therefore Hr = Hc, which defines the Base Failure–Cap Failure boundary. We therefore focus on discussing the critical pressure due to the Cap Failure mechanisms, before analyzing the critical height.

Depinned Cap Failure for reentrant geometries. To understand how ΔPc is influenced by the area fraction (or alternatively the cap width Wr) in the depinned Cap Failure mechanism, we begin by examining the rudimentary model proposed by Tuteja et al. (17, 29). In thisΔPc=4αsinθ°1WrαWr(4)

For convenience, we incorporate the shape parameter α, the same parameter as defined previously in the CAH section, to unify the critical pressure models on circular (α = π/4) and square (α = 1) geometries. To rationalize this model, at ΔPc, the pinning force of the contact line balances the force due to the pressure over the area between the pillars (17, 29). Two key assumptions are made. First, the contact line is supposed to follow the cap edge, while second, the contact angle around the contact line is presumed to be uniform and equal to θ°. We test this model in Fig. 3B (i). The square-cap model is observed to fit the simulation results very poorly, overestimating the critical pressure by between 26% and 95% in the tested range 0.016 ≤ Fr ≤ 0.8. If instead, a circular contact line model is used, we find that this agrees with the simulation data up to moderate area fractions (Fr < 0.6). Overall, by observing the contact line shape obtained through simulations, as shown in Fig. 3D, we find that the contact line varies in morphology, from circular at low Fr to approximately square at high Fr.

We now develop a more sophisticated model, capable of accurately describing the critical pressure for reentrant and doubly reentrant geometries, at all contact angles. Three modifications are introduced to Eq. 4. First, it cannot be assumed that the contact line follows the cap edge, leading to the introduction of Wr′, the corrected reduced width: Wr′ = Wra, where a is a parameter that describes the difference between the actual width that the contact line assumes and the width of the cap. Since a cannot be predicted a priori, we treat a as a second fitting parameter. For reentrant geometries and doubly reentrant geometries with θ° > 90°, we anticipate a ≈ 0, due to contact line pinning on the outer edge. For doubly reentrant geometries with θ° < 90°, we anticipate a ≈ 2Lr, due to contact line pinning on the inner edge. Second, we propose that the shape parameter α varies continuously as a function of Wr′ between the circular and square limits, such that α=π4+(1π4)(Wr)x. The exponent x describes the strength of this crossover and is a second fitting parameter. Third, if the contact angle on the hydrophilic reentrant geometry is increased, the pinning force of the contact line is maximized at θ° = 90°. For θ° > 90°, the pinning force remains at this maximal value (shown in the Supplementary Materials). Thus, to generally describe the critical pressure on all reentrant geometries, we replace θ° in Eq. 4 with the corrected contact angle θ°min=min[θ°,90°].

This crossover model is shown in Fig. 3B (i) to be in excellent agreement with the simulation data, yielding a = 0.023 and x = 6.7. As anticipated, a is small relative to the cap width and is of the order of the diffuse interface width (ϵ/B = 0.01). The large exponent x reflects the simultaneous change of both the perimeter and area of the contact line as the system crosses from a circular to square configuration.

Pinned Cap Failure for doubly reentrant geometries. Next, we consider the critical pressures of the doubly reentrant structure, as shown in Fig. 3B (ii). In the pinned Cap Failure region, the liquid-vapor interface is pinned to the inner cap edge, as shown in Fig. 3E. The doubly reentrant lip enforces an approximately square contact line across the entire range of Fr tested, such that excellent agreement between the simulations and the critical pressure model in Eq. 4 is achieved at x = 0, α = 1, and a = 0.080. This model is also successfully used for θ° = 110° in the Supplementary Materials. As anticipated, within the uncertainty introduced by the diffuse interface, a ≈ 2Lr (where 2Lr = 0.1). However, the data in Fig. 3B (ii) are best described by replacing θ° in Eq. 4 with 90°. In 2D, and for axisymmetric doubly reentrant wells, it is facile to show that the critical pressure occurs when the contact angle reaches 90° for thin lip widths lr (10). We conclude here that this remains true, even for 3D, nonaxisymmetric caps. This therefore verifies the proposition that the doubly reentrant lip maximizes the critical pressure for any surface wettability (16).

Critical heights. At a given Fr, the critical heights are obtained both analytically and through simulations by observing the maximum depth that the liquid-vapor interface sags under the pillar in the Cap Failure regions. If the pillar height Hr is equal to this sagging depth, the failure mechanism is simultaneously Base Failure and Cap Failure, thus defining the failure mode boundary and the critical height Hc. The salient observation based on the simulated critical heights (white data points in Fig. 3A, i and ii) is that regardless of cap area fraction, the optimal pillar height is unexpectedly short and should never exceed ≈ 0.2B. We now rationalize this observation.

In many critical pressure models, a 2D circular-arc model is used to estimate the sagging height of the liquid-vapor interfaceHc=Sr21cosθ°minsinθ°min(5)where the separation Sr = (1 − Wr) between horizontally adjacent pillars [see, for example, (17, 29)], or Sr=2(1Wr) for diagonally separated pillars [see, for example, (30)]. The latter model is shown in Fig. 3A (i and ii, dashed white lines), and exemplifies the conclusion that, except at very high Fr, a 2D estimation grossly overestimates the critical height. All currently manufactured, low-Fr structures relying on these 2D models are therefore substantially taller than necessary, which can be seen in (18) and (31), for example.

The actual, nonmonotonic critical height variation with Fr, as shown in Fig. 3A (i and ii), can be rationalized by considering that in 3D, two principal radii of curvature characterize the liquid-vapor interface at each point. At low Fr, the small contact line radius enforces a small, negative principal radius if curvature R1 on the liquid-vapor interface. Since the critical pressure is positive and proportional to 1/R1 + 1/R2, the second principal radius of curvature R2 must be smaller in magnitude than R1 and positive, resulting in a substantially reduced sag height compared to the 2D case. At large Fr, however, observable for the reentrant geometry in Fig. 3A (i), the 2D model is recovered as the principal radius of curvature approximated by the circular arc (R1) and is considerably smaller than the second principal radius of curvature of the interface (R2). In this case, the interface shape becomes well approximated by the single radius of curvature R1.

We therefore recognize that the liquid-vapor interface is able to be modeled as a capillary bridge, for which we define the inner radius to contact the cap at an angle θ°min and the outer radius to contact the simulation boundary with an angle equal to zero. Because the capillary bridge is axisymmetric, whereas the simulated system is square, there are two limiting cases of where the inner and outer radii contact the structure and simulation boundary, respectively, as shown in Fig. 3F (i). Either the inner and outer radii contact the central edges of the structure and simulation boundaries (the horizontal model, representing the minimum possible radii) or contact is made at the corners (the diagonal model, representing the maximum possible radii). Both potential capillary models are shown in Fig. 3A (i and ii, solid white lines). We reserve the detailed derivation to the Supplementary Materials.

In all cases, it is observed that the simulated critical heights are bounded by the horizontal and diagonal capillary bridge models. Furthermore, at low Fr, the horizontal capillary model accurately predicts the critical height. For structures where the interface is pinned to the outer cap edge, namely, reentrant geometries and doubly reentrant geometries with θ° > 90° (shown in the Supplementary Materials), the diagonal model is shown to closely predict the critical height at high Fr. An illustrated comparison of the horizontal and diagonal models is shown in Fig. 3F (ii). Through this, we are able to successfully capture the critical height suppression at low Fr and maximum Hc at intermediate Fr. We further validate the capillary bridge model in the Supplementary Materials by showing that this model is able to accurately reproduce an experimental interfacial profile from (32), whereas a circular arc model substantially overestimates the interfacial sagging height.

Minimum energy transition mechanisms

Transition states and pathways. To design surfaces that maintain a suspended state in challenging environments, it is not sufficient to only understand how susceptible a surface design is to pressure. Failure can be initiated through a broad range of additional perturbations, such as flow (33), vibration (34), evaporation (35), condensation (36), droplet impact (37), changes to electric (38) or magnetic (31) fields, or thermal fluctuations at the nanoscale (39). In a real application, several perturbations will be present simultaneously, meaning that failure is unlikely to be initiated by only a single perturbation but instead via their combination. In fabricating a texture resistant to failure, it is therefore vital to understand this combined failure in the worst-case scenario, the minimum-energy pathway (MEP) by which the suspended state collapses. This is a steepest-descent pathway between two metastable states, in which the maximum energy along the path occurs at a saddle point (the transition state). The minimum energy barrier is the difference in energies between the suspended and transition states. This places a lower bound on the collapse energy barrier. If this barrier cannot be overcome by the perturbations applied to a candidate surface design, then the suspended state can be guaranteed to remain stable.

Through utilizing the Doubly Nudged Elastic Band algorithm (40), three transition pathways are found: Base Contact (BC), Pillar Contact (PC), and Cap Contact (CC), visualized in movies S1 to S3, respectively. From these, the transition state morphologies are surveyed using the gradient-squared method, each of which is shown in Fig. 4 (A to C). In the large-scale structural surveys, as we only wish to obtain the transition states and minimum energy barriers rather than the full paths, the gradient squared method is considerably more efficient that the full-pathway algorithm, as only a single minimization is required, as opposed to evolving a string of multiple systems across the landscape. These transition search algorithms are detailed further in the Supplementary Materials. We begin by discussing the qualitative characteristics of each before quantifying the suspended-to-collapsed minimum energy barrier. All transition state searches are carried out at zero applied pressure.

Fig. 4 Reentrant and doubly reentrant transition state analysis.

(A to C) 3D visualizations of the transition states of each transition pathway (liquid-vapor interfaces shown in blue) with associated diagonal cross sections (liquid-vapor interface outlined in black). (C) Under-cap views show the three CC transition state morphologies. (D) 3D contour plots showing the energy barrier ΔEr of the lowest energy transition mechanism for the reentrant geometries. Each surface is a surface of constant ΔEr. The dividing surface between different transition mechanisms is shown in black. Impossible geometries with pillar widths Ar wider than the cap width Wr are shaded in dark gray. Geometries approaching this limit and requiring infeasibly large computational domains are shaded in light gray. (E) 3D contour plots showing the energy barrier ΔEr of the lowest energy transition mechanism for the doubly reentrant geometries.

The BC mechanism, as shown in Fig. 4A, is observed for both reentrant and doubly reentrant geometries. The mechanism proceeds via the liquid-vapor interface sagging toward the system base while pinned to the cap lip. The transition state is observed to occur after the liquid meniscus has contacted the base of the system. This mechanism is highly prevalent across a broad range of structural and chemical properties, such as on nonwetting geometries (see the Supplementary Materials) as well as pillars and nails (39, 4144).

The PC mechanism, illustrated in Fig. 4B, occurs only for the reentrant geometry. The mechanisms proceeds via the liquid-vapor interface impinging under the cap, with the transition state occurring before the interface detaches from the solid surface. Following this, the interface slides down the pillar to contact the base and complete the transition.

Last, the CC mechanism is observed only on the doubly reentrant geometries, as shown in Fig. 4C. Here, the expected PC mechanism is unstable with respect to condensation of liquid inside the cap structure. In this mechanism, the transition begins with the condensation of liquid in one corner of the cap underside, which subsequently grows to fill the cap. The precise location of the transition state can take one of three morphologies, as shown in Fig. 4C (i to iii). Which variant occurs is discussed in the Supplementary Materials. The transition then continues by filling the cap entirely, such that a new free energy minimum is obtained, a suspended state with a liquid-filled cap. An additional energy barrier is required to complete the wetting mechanism: the two separate liquid-vapor interfaces must coalesce by crossing the cap lip before the remainder of the mechanism proceeds exactly as with PC. However, this coalescence barrier is small relative to the condensation barrier presented and decreases further as the lip width is reduced. This mechanism is particularly important to understand as for many applications requiring the doubly reentrant geometry, the liquids are volatile, or else the suspended state needs to be maintained over long time scales. Two such applications are discussed further in subsection Simultaneous Optimization.

Of further note is that the transitions presented here are all MEPs regardless of liquid volatility. The pressure treatment in the free energy functional used in these simulations effectively contacts every point of the system with an external fluid reservoir at constant pressure. Thus, fluids may be exchanged anywhere within the system. For nonvolatile liquids, although the CC condensation mechanism is a minimum energy pathway on the doubly reentrant geometry, CC may not be realized on an experimental time scale. The transition will therefore occur via a noncondensing route, the minimum energy path of which is BC.

Minimum energy barriers. Overall, we find that each structure has, at most, two potential transition pathways. Throughout this work, the energy barrier ΔEr is expressed relative to the reference energy γlvB2. The barrier of each pathway is affected differently by the structural parameters, which we summarize in Table 1. In the style of traditional phase diagrams, we present the lowest energy barrier mechanism at each parameter value tested in Fig. 4 (D and E) and so are able to predict the dominant collapse mechanism.

Table 1 Geometrical parameters that affect ΔEr for each transition mechanism, indicated with filled circles.

View this table:

Beginning with the reentrant geometry, BC and PC compete for the lowest energy collapse mechanism, as shown in Fig. 4D. As shown in Table 1, for all BC transition states, the interface morphology depends on the pillar height Hr and cap width Wr. As both of these structural parameters increases, the liquid-vapor interface increases in area, leading to an increase in ΔEr. However, ΔEr cannot be increased indefinitely by increasing Hr; as at a critical pillar height, the height-independent PC mechanism becomes the lowest energy transition pathway. As the PC transition state is associated with the liquid wetting the reentrant cap underside, ΔEr is increased by expanding the liquid-vapor interfacial area required to do so. This requires Wr to be maximized and the pillar width Ar to be minimized. In designing a reentrant structure exhibiting the maximum energy barrier, the mechanistic switch upon increasing Hr is a key point to highlight as, assuming taller pillars are mechanically weaker than shorter pillars (27, 28), at a given Wr and Ar, the optimal structure height is that on the BC-PC boundary, reminiscent of the critical pressure case.

Figure 4E shows the lowest energy collapse mechanism for the doubly reentrant geometries, which is markedly different to the (singly) reentrant equivalent. Here, it is the BC and CC mechanisms that compete for the lowest energy path. In Table 1, we find that the minimum energy barrier depends on four parameters: the CC energy barrier depends on Wr, Ar, and the lip depth Lr, whereas the BC barrier depends on Wr and also Hr. However, ubiquitously, the lowest CC barriers are obtained by minimizing Ar, such that the barrier diagram shown in Fig. 4E is at constant Ar = 0.05.

Regarding the BC mechanism, the only effect on changing from a reentrant to doubly reentrant geometry is to increase the range of Wr and Hr for which the BC mechanism is operative. However, except for the smallest pillar heights, the CC mechanism has a considerably smaller ΔEr compared to the BC barrier. This is principally caused by the condensing critical nucleus having a relatively small, energetically unfavorable liquid-vapor interfacial area, compared to the large, energetically favorable solid-liquid interfacial area. Therefore, to maximize the CC barrier, the liquid-vapor interfacial area must be maximized, while minimizing the solid-liquid contact area. This is effectively realized by maximizing Wr and minimizing Lr and Ar. As the CC mechanism is independent of the height of the structure, in a similar manner to the reentrant structure, the optimal pillar height is located on the boundary between the two failure mechanisms.

Simultaneous optimization

Overall, six structural parameters influence the three key wetting properties: Ar (pillar width), Hr (pillar height), Lr (lip depth), tr (cap thickness), Wr (cap width), and the system scale B/Bref (where Bref = 1 μm). Having studied how these parameters affect each individual wetting property, the parameters that antagonistically couple the wetting properties become apparent. First, to reduce the CAH, Wr must be reduced; but this reduces ΔPc and ΔEr. Second, to increase ΔPc, the system scale must be reduced; but this reduces ΔEr. To overcome this unfavorable coupling, we simultaneously optimize the surface structures, which is demonstrated for two example applications: membranes for water purification via membrane distillation, and droplet-based digital microfluidics. To perform this simultaneous optimization, we begin by developing an application-specific scoring function that grades a candidate design against the desired wetting properties. We then optimize the scoring function using two methods. The first is to evaluate the scoring function over the entire parameter range tested, and from this, find the optimum structure. In the second method, we demonstrate that for designs where it is not practical to perform comprehensive wetting property surveys, because of increased surface complexity, for example, genetic algorithms can be used to efficiently perform the simultaneous optimization.

Conventionally in membrane distillation, purification is achieved by passing a heated-contaminated water source over a hydrophobic membrane, through which water vapor is able to pass to collect in a clean water reservoir. However, oils readily foul the membrane, leading to breakthrough of the contaminated liquid into the fresh water reservoir (2, 45). We overcome this using a doubly reentrant structure. To optimize the geometry, we construct a suitable scoring by recognizing that the first priority is for the membrane to be pressure resistant under typical operating conditions: water at 70°C (γlv = 64.4 mN m−1) under pressures of approximately 100 kPa (ΔPctarget) (46). To ensure that the suspended state remains stable, the minimum energy barrier must be of the order of 100 kBT (ΔErtarget). This is of particular importance for the doubly reentrant geometry to prevent failure via condensation within the texture (47). Last, to reduce viscous drag across the membrane, we minimize the CAH and impose the condition that the CAH should not exceed 90° (CAHcutoff). The critical pressure, energy barrier, and CAH conditions generate individual scoring functions SP, SE, and SC, respectivelySP=12[1+tanh(ΔPcΔPctargetΔPwidth)]SB=12[1+tanh(ΔErΔErtargetΔEwidth)]SC=max(CAHcutoffCAHCAHcutoff,0)(6)from which the total score is the geometric mean of these: Score=(SPSBSC)13. For SP and SB, a tanh profile is selected to appropriately localize critical pressures and energy barriers within a range of suitable operating conditions. This leads to the widths chosen here as ΔPwidth = 0.5 and ΔEwidth = 5 × 10−5. Meanwhile, the linear function for SC aims to ensure that low-CAH structures are always favored.

By maximizing this six-dimensional scoring function using either the results from the wetting property survey or a genetic algorithm, the optimal structure is obtained with a score of 0.794. The optimal parameters are (Ar, Hr, Lr, tr, Wr, B/Bref) = (0.05, ≥0.17, 0.17, 0.05, 0.27, 0.32). The optimum system scale of 320 nm is notably similar to that of springtail cuticles (9). Both the springtail cuticle and membrane design have been selected for pressure-resistant liquid shedding ability, while allowing the unimpeded movement of gasses through the surface. The membrane design proposed here may therefore reflect a natural optimum for robust gaseous diffusion. The optimum design yields the properties: ΔPc = 162 kPa, ΔE = 1.25 × 103 kBT, and CAH = 42° (θa = 165°, θr = 123°). This CAH is typical of currently manufactured reentrant microtextures [see, for example, (27, 29, 31)].

By studying the individual wetting properties, we can rationalize the optimal structural design. The optimal value of Ar represents the minimum pillar width tested, whose sole function is to maximize ΔEr. Hr reflects the observation that the maximum critical pressure is achieved at HrHc. Lr + tr optimizes the CAH, while the specific value of Lr maximizes ΔEr. Last, the small value of Wr reduces the CAH while retaining a high ΔPc due to the small system scale.

The scoring function at fixed optimal values of Ar, Hr, and tr is shown in Fig. 5A (i) as a 3D contour plot. A 2D cut through this is shown in Fig. 5A (ii) at the optimal Lr, to show that the optimal structure scale is bounded by the critical pressure criterion from above and the minimum energy barrier criterion from below. Also shown in Fig. 5A (ii) are projections of successive generations of the genetic algorithm. The optimal structure was located after 20 generations, requiring the sampling of only 0.01% of the 7.2 × 106 possible structures considered overall.

Fig. 5 Simultaneous optimization of the three wetting properties for membrane distillation and digital microfluidics applications.

(A) (i) 3D contour plot of the membrane distillation scoring function at fixed Hr = 0.3, Ar = 0.05, and tr = 0.05. Each surface is a surface of constant score. (ii) A 2D slice of the 3D contour plot at the optimal Lr = 0.17. Square data points show the initial (white), second (light gray), fifth (dark gray), and final (black) generations of the genetic algorithm, projected onto the 2D plane. (B) Scoring function for the digital microfluidics application, projected onto the Hr = 0.3 plane at fixed B = 100 μm, also showing the successive generations of the genetic algorithm population.

The four examples of manufactured doubly reentrant surfaces feature system sizes of 1 to 100 μm, as smaller scales are currently challenging to manufacture (1316). We now choose to optimize a structure whose manufacture has already been demonstrated, so we fix B = 100 μm in accordance with the texture designed in a recent work by Liu and Kim (16). A leading-edge application for this is surfaces designed for digital microfluidics (8, 21). Devices fail when droplets become immobilized by pinning to the surface, particularly problematic for low–surface tension solvents (21) or for reactive processes where the surface tension is variable and hard to predict.

Both of these problems are readily overcome using the doubly reentrant geometry. We demonstrate this by optimizing a surface structure for use within a particularly challenging scenario—digital microfluidics using microliter volume droplets of n-hexane (γlv = 27.4 mN m−1). The pressure within such a droplet (88 Pa) introduces the target pressure for use in the scoring scheme Eq. 6: ΔPctarget = 100 Pa, with the width ΔPwidth = 0.005. Furthermore, as the CAH should be minimized, we impose the condition CAHcutoff = 50°. At the imposed length scale used, the barrier score SB ≈ 1 (ΔE is of the order of 10 × 1010 kBT), meaning that we choose to optimize the score (SPSC)1/2.

Overall, the optimum structure, 0.508, is obtained using both the wetting property survey and the genetic algorithm. The optimum parameters are (Hr, Dr, and Wr) = (≥0.13, 0.22, and 0.16), yielding the properties: ΔPc = 105 Pa, CAH = 25°. The pillar width Ar and ratio of Lr to tr become free parameters to choose. A 2D contour plot of the scoring function at constant Hr = 0.3 is shown in Fig. 5B, in which projections of successive generations of the genetic algorithm are shown. Here, the algorithm converged after 14 generations, requiring 2.2% of the entire population to be sampled.

For both the membrane distillation and digital microfluidics applications, the sensitivity of the optimized structural dimensions and properties can be assessed relative to the choice of scoring function parameters. This is achieved through reoptimizing the geometries when each parameter in the scoring functions shown in Eq. 6 were varied individually by ±5%. It is found that the optimized membrane distillation geometry is insensitive to the parameter variation. This observation also applies to the optimized digital microfluidics geometry, except in the scenario where CAHcutoff is reduced by 5%. In this case, the optimal CAH is reduced by 13%, and the optimal critical pressure is reduced by 17%. However, this variation is due to the CAHcutoff reduction (2.5°) being relatively large compared to the low optimal CAH (25°) for this application.

The manufactured structure reported by Liu and Kim (16) is similar to the optimum geometry. However, the key difference is that the optimal geometry has a considerably shorter pillar height Hr than the manufactured geometry by a factor of 3.7 times. This is due to the unexpectedly short critical height required, as discussed in the critical pressure section.


Overall, to optimize the wetting properties of the reentrant and doubly reentrant surface texture for a vast variety of potential applications, we began by developing three computational strategies to comprehensively survey the key surface wetting properties: the CAH, critical pressure, and minimum energy barrier to the wetting transition. This was achieved for both wetting and nonwetting liquids (shown in the Supplementary Materials).

In the CAH study, we identified four major receding mechanisms, of which only two had previously been reported, and defined the structural dimensions where each is operative. For all receding mechanisms, we were able to develop and analyze quantitative models that were robustly validated against our simulation results.

In the critical pressure study, three failure mechanisms were observed and quantified as a function of the structural parameters. However, upon comparison with the simulation data, the prevailing and widely used critical pressure models were found to be considerably oversimplified. This leads to a particularly poor description of the liquid-vapor interface morphology, meaning that manufactured structures are many times taller (and mechanically weaker) than necessary. By developing a more sophisticated model, we were able to achieve both quantitative accuracy of the critical pressures and success at modeling the complex interface morphologies as capillary bridges.

In the minimum energy barrier study, we identified three failure mechanisms, quantified each barrier across the structural parameter space, and assessed which mechanism was most likely for a given geometry. Crucially, we not only showed how the doubly reentrant geometry is prone to condensation within the cap but also deduced effective designs to mitigate against this.

Last, we found the structural features that tend to maximize the critical pressure, minimize the energy barrier, and maximize the CAH. As it was not possible to optimize a surface geometry with respect to each individual wetting parameter, we performed the optimization by considering all three simultaneously. This was achieved in two ways for the optimum design of both membranes for water purification and surfaces for digital microfluidics. First, using the comprehensive wetting property surveys, we were able to evaluate and locate the maximum of a combined scoring function. However, we then demonstrated that a genetic algorithm was able to efficiently locate the optimum design in the six-dimensional parameter space. Although the designs tested here featured a relatively small number structural degrees of freedom, going forward, we highlight such optimization techniques as being powerful tools in designing more complex structures for special wettability applications. The computational techniques developed here are highly versatile and can be used for any mesoscopically structured surface in contact with multiple fluid phases.

Coupled with recent, substantial developments in fabrication techniques (including 3D printing, fluidization of polymer micropillars, and lithographic methods), we believe that the multifaceted optimization strategy presented here will be a powerful approach to designing real-world superomniphobic surfaces. In the future, an additional step will be to consider the mechanical reliability and scalability of manufactured designs, in which we have contributed to this discussion with our large reduction in the necessary pillar height.


Diffuse interface model

The simulations used to compute the CAH, critical pressure, and minimum energy barrier all used the same diffuse interface model and system discretization. Specific system setups for each wetting property are presented in the Supplementary Materials. Within the bifluidic diffuse interface model used, the order parameter ϕ(r) was chosen to represent the local composition at point r (ϕ = 1 in the pure liquid phase, and ϕ = −1 in the pure vapor phase at zero applied pressure). On the basis of a previous work, the free energy functional Ψ[ϕ] is composed of three terms (41)Ψ[ϕ]=Ψi+ΨsΔPVl(7)Ψi is the isotropic free energy term, expressed as an integral over the entire system volume VΨi=V(1ϵ(14ϕ412ϕ2+14)+ϵ2ϕ2)dV(8)enforcing the equilibrium values of ϕ = ±1 via the double well potential, and exacting an energetic penalty for forming an interface of width ϵ. This leads the liquid-vapor surface tension γlv=8/9.

Ψi is the fluid-solid interaction term, expressed as an integral over the surface areaΨs=Sh(16ϕs3+12ϕs+13)dS(9)where ϕs is the value of ϕ at the surface. h is the wetting parameter and is related to the microscopic contact angle θ° through h=2cosθ°. The cubic wetting potential negates spurious compositional changes close to the surface by ensuring that ϕs is equal to the bulk composition (48).

Within the external pressure term, the total liquid volume is calculated fromVl=Vϕ+12dV(10)

The simulation system is discretized into a cubic lattice of Nx × Ny × Nz nodes, in which each node is either located within the solid structure, on the solid surface, or in the fluid bulk.

Genetic algorithm

We began by randomly sampling the parameter space to generate an initial population of 40 surface structures. These were ranked on the basis of score, and the top 20 were retained for breeding. Candidate pairs for breeding were selected at random, and breeding occurred if the geometric mean of their scores was greater than a random number between 0 and 1. The offspring were equally likely to inherit each attribute from either parent. Each attribute was then mutated if a random number was less than the current mutation probability pi (set initially at 0.5). For the discrete structural variables, a mutation changed the attribute randomly by between −3 and 3 lattice spacings. For the continuous variable (system scale), the change in log(B/Bref) was selected randomly from the range −0.3 to 0.3. The mutation probability was reduced each generation, such that pi+1 = p0× (SD of scores in previous generation)1/2. Any offspring bred or mutated outside the testable parameter range was mutated back into the testable parameter range. The algorithm was deemed to have converged when the mutation rate decreased to zero.


Supplementary material for this article is available at

Section S1. Supplementary Materials and Methods

Section S2. Extended discussion for θ° = 60°

Section S3. Results and discussions for θ° = 110°

Fig. S1. CAH simulation setups.

Fig. S2. Decomposition of CAH into advancing and receding contact angles.

Fig. S3. Analytic models for the receding mechanisms.

Fig. S4. Illustrated critical pressure model parameters.

Fig. S5. Analysis and comparison of the circular arc and nodoid models of the sagging depth of the liquid-vapor interface at the critical pressure.

Fig. S6. Effects of the reentrant and doubly reentrant structural parameters on each transition type.

Fig. S7. Extended analysis of the lowest energy reentrant and doubly reentrant transition mechanisms.

Fig. S8. CAH mechanisms for θ° = 110°.

Fig. S9. Critical pressure analysis for θ° = 110°.

Fig. S10. Individual and lowest energy transition mechanisms for θ° = 110°.

Movie S1. BC mechanism for a doubly reentrant geometry at θ° = 60°.

Movie S2. PC mechanism for a reentrant geometry at θ° = 60°.

Movie S3. CC mechanism for a doubly reentrant geometry at θ° = 60°.

Movie S4. PC mechanism for a reentrant geometry at θ° = 110°.

References (4954)

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Acknowledgments: We would like to thank C. Semprebon, I. Liu, and E. Xi for useful discussions and C. M. Jones for updating the energy minimization software. Funding: We thank P&G and EPSRC (EP/P007139/1) for funding. Author contributions: H.K. and Y.G. conceived the idea and designed the research. J.R.P. performed the simulations and the data analysis with supervision from H.K. J.R.P., Y.G., and H.K. discussed the results and wrote the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Data are also available at doi:10.15128/r1m039k492r. The code used to generate the data may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article