Two-Dimensional Numerical Analysis of Gas Diffusion-Induced Cassie to Wenzel State Transition

Mayer, Michael D.; Kadoko, Jonah; Hodes, Marc

doi:10.1115/1.4051320

Abstract

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.

1 Introduction

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 [4]. In liquid flow in the Cassie state over SHs, this reduction in the solid–liquid interfacial area reduces frictional drag [5]; 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 $p_{g, c} (0^{+})$ [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 (⁠ $θ_{a}$ ⁠) [8,11]. We refer to this time as a critical time $(t_{cr})$ ⁠. 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 [8]. 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 $t_{f}$ [10,12]².

The upward component of the surface tension force during the metastable Cassie state depends on the shape of the cavity sidewalls [13]. 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 [13]. 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 [15]. Lv et al. [8] 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. [3] 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. [11]. 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. [14] 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 [16]. 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. [3] 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. [3] 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. [3]. However, we note that Xu et al. [3] 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. [3] 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

We consider a 2D quiescent liquid column that is initially in the Cassie state above a ridge-type structured surface as per the two configurations shown in Fig. 1. Our focus is restricted to half of a single cell due to periodicity, as discussed in detail later. According to the saturation table for water and the Kelvin Equation, at 25 deg, in the limiting case of a meniscus with a contact angle of 110 deg, the vapor pressure is about 4% of standard atmospheric pressure. 110 deg is the advancing contact angle of water on silicon treated with a flouroploymer such as PTFE, [17], the scenario considered here. Due to the relatively small effect of vapor pressure, we can ignore the vapor component of the mixture which we refer to as the “gas” phase. Additionally, to simplify the model we consider a single gas species in the cavity, in our case nitrogen. In Fig. 1,

h (t)

is the height of the portion of the liquid domain, Ω, between the top of a ridge and the (assumed flat) top surface of the liquid, s₀ is the ridge height (or cavity depth), 2a is the width of the cavity and 2d is the ridge period such that half of the ridge thickness is (d–a). We additionally define the solid fraction of the ridge as

ϕ = (d - a) / d

⁠. Assuming the pressure on the liquid side of the meniscus is equal to or larger than that on the gas side of it, a stress balance along the meniscus is given by

p_{g, c} (t) = p_{\infty} + ρ_{l} g [h (t) - y_{m} (x, t)] - \frac{σ}{R (x, t)}

(1)

where

p_{g, c}

is the gas pressure in the cavity,

p_{\infty}

is atmospheric pressure,

ρ_{l}

is the density of the liquid, g is the acceleration due to gravity,

y_{m} (x, t)

is the y-coordinate of the meniscus, σ is the liquid–gas surface tension and R is the radius of curvature of the meniscus.

y_{m} (x, t)

is a moving boundary and its location must be tracked during our solution. We assume that

h (t) ≫ | y_{m} (x, t) | \to h (t) - y_{m} (x, t) \approx h (t)

⁠; therefore, changes in hydrostatic pressure along the meniscus are negligible and R is independent of x, i.e., the meniscus is a circular arc. Consequently, Eq. simplifies to

p_{g, c} (t) = p_{\infty} + ρ_{l} g h_{0} - \frac{σ}{R (t)}

(2)

Fig. 1

View large Download slide

Schematics with liquid initially in the Cassie state over one pitch of ridges: (a) hydrostatic pressure neglected and top surface of liquid is impermeable, (b) hydrostatic pressure considered and top surface of liquid surface exposed to noncondensible gas. This is a single cell of a periodic surface. As the lines of periodicity are through the middle of the parallel ridges we refer to the solid thickness as a “half-ridge.”

where

h_{0} = h (t = 0)

³. Also, it follows from the geometry of the problem that

y_{m} (x, t) = \pm [\sqrt{R {(t)}^{2} - x^{2}} - \sqrt{R {(t)}^{2} - a^{2}}] - s (t)

(3)

where s(t) is the distance (a positive quantity) from the top of a ridge to the triple contact line and the ‘+’ and ‘–’ signs preceding the bracketed term correspond to menisci protruding upward and downward, respectively. Henceforth, we restrict our attention to horizontal and downward protruding menisci. As per Fig. 1, the contact angle, i.e., the (positive) angle between the vertical and the tangent to the meniscus at the triple contact line, is between

π / 2

and the advancing contact angle. Before critical time, i.e, when

s (t) = 0

⁠, the meniscus is pinned and R decreases with time. After the critical time it assumes a constant value corresponding to the advancing contact angle. The geometric relationship between R and θ is

R (t) = {\begin{matrix} \frac{a}{cos [π - θ (t)]} & for 0 < t \leq t_{cr} \\ \frac{a}{cos (π - θ_{a})} & for t_{cr} \leq t \leq t_{f} \end{matrix}

(4)

As mentioned before, ridge-type SHs are periodic structures so we can focus our solution on a single half-cell of the surface as shown in Fig. 1. The cavities are closed to the environment and thus the only mechanism for gas transport is diffusion into the liquid. Furthermore, we analyze the liquid domain on the interval

0 \leq x \leq d

because it's symmetric about x = 0. The volume of gas per unit depth in one half of a cavity,

V_{g, c}^{'} (t)

⁠, is given by

V_{g, c}^{'} (t = 0) = a s_{0}

plus the integral of

y_{m} (x, t)

on

0 \leq x \leq a

⁠. It follows that

V_{g, c}^{'} (t) = a s_{0} - a s (t) - \frac{1}{2} {R {(t)}^{2} \arcsin [\frac{a}{R (t)}] - a \sqrt{R {(t)}^{2} - a^{2}}}

(5)

where

R (t)

can be represented in terms of gas pressure as per Eq. . The gas in the cavity obeys the ideal gas law

p_{g, c} (t) V_{g, c}^{'} (t) = \frac{m_{g, c}^{'} (t)}{M_{g}} R_{u} T

(6)

where $m_{g, c}^{'}$ is the residual mass of gas in the cavity per unit depth, $M_{g}$ is the molecular weight of the gas, $R_{u}$ is the universal gas constant and T is the absolute temperature of the gas.

We consider a scenario where a liquid is slowly, as to not allow any impingement of the surface, poured onto an SH, and then immediately gas begins to diffuse into the liquid. Our initial time is the time immediately after the liquid has been poured. In Fig. 1(a), we depict a case study I, where an impermeable membrane separates the liquid from the pressure of the surrounding atmosphere (⁠

p_{\infty}

⁠), assumed to be much greater than the hydrostatic pressure of the liquid column, i.e.,

p_{\infty} ≫ ρ_{L} g h (t)

⁠. In this case, at t = 0, gas pressure in the cavity,

p_{g, c} (0)

⁠, is equal to

p_{\infty}

and thus, the meniscus is flat. Figure 1(b) shows a schematic of case study

II

⁠, where we consider the effects of hydrostatic pressure and gas diffusion at the top of the liquid domain. Because of hydrostatic pressure, the initial pressure in the cavity will be higher than that of the atmosphere and the meniscus has an initial radius of curvature of

R (0^{+})

[11]. To solve for the initial volume, we assume that the gas pressure in the cavity is initially (i.e., at t = 0)

p_{\infty}

⁠. However, we assume that the gas is instantly compressed to a pressure of

p_{g, c} (0^{+})

and a volume per unit depth of

V_{g, c}^{'}

due to hydrostatic pressure. Assuming isothermal compression,

p_{g, c} (0^{+})

is given by

p_{\infty} a s_{0} = p_{g, c} (0^{+}) V_{g, c}^{'} (0^{+})

(7)

The meniscus has been assumed flat $p_{\infty}$ 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., $s (0^{+}) \neq 0$ and $t_{cr} = 0$ ⁠, if $V_{g, c}^{'} (0^{+}) < V_{g, c}^{'} (θ = θ_{a})$ ⁠.

We neglect the advective mass transport of the dissolved gas species, a valid assumption in analogous mass transfer problems [18]. Thus, the two-dimensional, transient diffusion equation governs the dissolved gas concentration field as per

\frac{\partial ρ_{g}}{\partial t} = D_{gl} (\frac{\partial^{2} ρ_{g}}{\partial x^{2}} + \frac{\partial^{2} ρ_{g}}{\partial y^{2}})

(8)

where

ρ_{g} (x, y, t)

is the partial density of the dissolved gas in the liquid and

D_{gl}

is the binary molecular diffusivity of the dissolved gas in the liquid. We prescribe symmetry boundary conditions at x = 0 and x = d as per Eqs. and . The cavity is sealed on all sides; therefore, the gas leaves it only through diffusion across the meniscus. Dissolved gas concentration on the liquid side of the meniscus is determined from Henry's Law as per Eq. , where H is Henry's constant. Case study I and case study II only differ in the boundary condition at the top of the domain and share the following boundary conditions.

\frac{\partial}{\partial x} ρ_{g} (0, y_{m} (x, t) < y < h (t), t > 0) = 0

(9)

\frac{\partial}{\partial x} ρ_{g} (d, 0 < y < h (t), t > 0) = 0

(10)

ρ_{g} (0 < x < a, y_{m} (x, t), t) = p_{g, c} (t) H

(11)

\frac{\partial}{\partial y} ρ_{g} (a < x < d, 0, t > 0) = 0

(12)

\frac{\partial}{\partial x} ρ_{g} (a, - s (t) < y < 0, t_{cr} < t < t_{f}) = 0

(13)

For case study

I

⁠, we assume that the top of the domain is impermeable to gas as per

\frac{\partial}{\partial y} ρ_{g} (0 < x < d, h (t), t > 0) = 0, case study I

(14)

For case study

II

⁠, we assume that the top of the domain is open to the environment such that the dissolved gas concentration on the liquid side of the top of the domain given by Henry's Law as per

ρ_{g} (0 < x < d, h (t), t > 0) = p_{\infty} H, case study II

(15)

The initial condition is uniformly dissolved gas concentration as per

ρ_{g} (x, y, 0) = ρ_{g, 0}

(16)

Recall that, for case studies I and II, the initial pressures of gas in the cavity are

p_{\infty}

and

p_{g, c} (0^{+})

⁠, respectively. In the liquid domain, the total mass of the solvent species is preserved. Moreover, since the partial molar volume of the dissolved gas is much smaller than that of the liquid, the total change in volume of the liquid domain is negligible. It follows from the geometry of the problem that the height of the liquid above the ridge is

h (t) = h_{0} - \frac{a}{d} s (t) - \frac{1}{2 d} {R {(t)}^{2} \arcsin [\frac{a}{R (t)}] - a \sqrt{R {(t)}^{2} - a^{2}}}

(17)

We bound the gas in the cavity by a deforming or moving control volume (of unit depth),

{cv}_{g} (t)

⁠. The only flux of gas through the corresponding deforming control surface,

{cs}_{g} (t)

⁠, is due to the diffusion of it along with the meniscus. Applying conservation of mass to the gas in the cavity yields

0 = \frac{d}{d t} \int_{c v_{g} (t)} ρ_{g, c} d A - D_{gl} \int_{c s_{g} (t)} \nabla ρ_{g} \cdot \hat{n} d s

(18)

where

ρ_{g, c} = m_{g, c}^{'} / V_{g, c}^{'}

is the density of gas in the cavity as per Eq. , assumed uniform because of the high molecular diffusivity of gases,

d A

is a differential area in the domain,

\hat{n}

is the unit vector normal to the meniscus and pointing into the liquid as per

\hat{n} = \frac{- \partial y_{m} / \partial x i + j}{\sqrt{1 + {(\partial y_{m} / \partial x)}^{2}}}

(19)

and s is the coordinate along with the meniscus. Upon substituting Eq. into Eq. and transforming the surface integral along

d s

in Eq. to one along

d x

and integrating over time we obtain

0 = m_{g, c}^{'} (t) - m_{g, c}^{'} (0) - D_{gl} \int_{0}^{t} \int_{0}^{a} ({\frac{\partial ρ_{g}}{\partial y} |}_{y = y_{m}} - \frac{\partial y_{m}}{\partial x} {\frac{\partial ρ_{g}}{\partial x} |}_{y = y_{m}}) d x d τ

(20)

Here,

m_{g, c}^{'} (t)

and

m_{g, c}^{'} (t = 0)

may be expressed as functions of the pressure and volume of gas in the cavity as per Eq. , the ideal gas law. Moreover, the volume of gas in the cavity,

V_{g, c}^{'} (t)

⁠, is a function of its pressure, from which the radius of curvature of the meniscus follows from Eq. , and the displacement of the triple contact line. It follows that Eq. is equivalent to

0 = \frac{M_{g}}{R_{u} T} {p_{g, c} (t) V_{g}^{'} [p_{g, c} (t), s (t)] - p_{\infty} s_{0} a} - D_{gl} \int_{0}^{t} \int_{0}^{a} {[\frac{\partial ρ_{g}}{\partial y} - \frac{\partial y_{m}}{\partial x} \frac{\partial ρ_{g}}{\partial x}]}_{y = y_{m}} d x d τ

(21)

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

3.1 Introduction.

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 $q = 1, 2, \dots, Q, Q + 1, \dots, 2 Q$ for $t \leq t_{cr}$ and $q = 1, 2, \dots, Q, Q + 1, \dots, 2 Q, 2 Q + 1$ for $t > t_{cr}$ as shown in Fig. 2. The columns of subdomains to the left and right of the triple contact point at $x = a, y = 0$ for $t \leq t_{cr}$ and $x = a, y = - s (t)$ for $t > t_{cr}$ are separated by a (vertical) line emanating from it. There are Q rows and 2 columns of subdomains up until $t_{cr}$ ⁠, when subdomain $2 Q + 1$ 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 y_q, where $q = 0, 1, 2, \dots, Q - 1$ ⁠, except at the top of the domain, where the value of y decreases with time and corresponds to $h (t)$ ⁠. 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, $Ω_{1}$ ⁠, is bounded by a deforming boundary, $y = y_{m} (x, t)$ ⁠, and a fixed one, y = y₁. At critical time, the meniscus depins and the triple contact line moves down the impermeable ridge; therefore, we separate $Ω_{1}$ into two subdomains as per Fig. 2(b). Then, subdomain $Ω_{1}$ becomes rectangular and bounded from below by y = 0 and subdomain $Ω_{2 Q + 1}$ contains the rest of the liquid previously in subdomain $Ω_{1}$ ⁠. 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.

Fig. 2

Labeling of (a) 2Q subdomains in case study I (where subdomain 1 is rectangular at t = 0) when 0≤t≤tcr and II when 0+≤t≤tcr and (b) 2Q+1 ones when tcr<t≤tf. At t=tcr, the portion of subdomain 1 below y = 0 becomes subdomain 2Q+1.

View large Download slide

Labeling of (a) 2Q subdomains in case study I (where subdomain 1 is rectangular at t = 0) when $0 \leq t \leq t_{cr}$ and II when $0^{+} \leq t \leq t_{cr}$ and (b) $2 Q + 1$ ones when $t_{cr} < t \leq t_{f}$ ⁠. At $t = t_{cr}$ ⁠, the portion of subdomain 1 below y = 0 becomes subdomain $2 Q + 1$ ⁠.

We map each subdomain to the canonical

[- 1, 1] \times [- 1, 1]

(square) domain suitable for Chebyshev collocation and obtain machinery to effectively accommodate meniscus curvature. The transformed variables ξ_q and η_q are given by the maps

ξ_{q} (x) = 2 \frac{x - x_{q, -}}{x_{q, +} - x_{q, -}} - 1

(22)

η_{q} (x, y, t) = 2 \frac{y - y_{q, -} (x, t)}{y_{q, +} (t) - y_{q, -} (x, t)} - 1

(23)

where $x_{q, -}, x_{q, +}, y_{q, -} (x, t)$ ⁠, and $y_{q, +} (t)$ are the positions of the left, right, lower and upper boundaries, respectively, on subdomain q [19]. Here $y_{q, -}$ only depends upon x and t for subdomain 1 when $t \leq t_{cr}$ and subdomain $2 Q + 1$ when $t > t_{cr}$ and $y_{q, +}$ only depends on t for subdomains Q and 2Q.

Next, we transform the dependent variable from

ρ (x, y, t)

to

ρ [ξ (x), η (x, y, t), t]

⁠. The first-order partial derivatives of

ρ (x, y, t)

with respect to x and y in terms of (ordinary and partial) derivatives of the maps and those of

ρ (ξ, η, t)

are

ρ_{x} = ξ_{x} ρ_{ξ} + η_{x} ρ_{η}

(24)

ρ_{y} = η_{y} ρ_{η}

(25)

where the various derivatives are denoted by subscripts. Noting that ξ_xx, ξ_y, ξ_yy and η_yy equal 0 in all subdomains, the second-order partial derivatives are

ρ_{x x} = ξ_{x}^{2} ρ_{ξ ξ} + 2 ξ_{x} η_{x} ρ_{ξ η} + η_{x}^{2} ρ_{η η} + η_{x x} ρ_{η}

(26)

ρ_{y y} = η_{y}^{2} ρ_{η η}

(27)

Too, η_xx is only nonzero on subdomain 1 before critical time and subdomain

2 Q + 1

after it. Furthermore,

{ρ_{t} |}_{ξ, η} = {ρ_{t} |}_{x, y} + ρ_{η} η_{y} y_{t}

(28)

where the second term on the right side accounts for the deformation rates in y of physical subdomains 1, Q and 2Q for

t \leq t_{cr}

and Q, 2Q, and

2 Q + 1

for

t > t_{cr}

when ξ and η are held constant [20,21]. It follows that the transient diffusion equation, Eq. , in the (Cartesian)

ξ η

plane is

{ρ_{t} |}_{ξ, η} - ρ_{η} η_{y} y_{t} = D_{gl} [ξ_{x}^{2} ρ_{ξ ξ} + 2 ξ_{x} η_{x} ρ_{ξ η} + (η_{x}^{2} + η_{y}^{2}) ρ_{η η} + η_{x x} ρ_{η}]

(29)

The dissolved gas concentration profile in each subdomain is approximated as a Lagrangian polynomial as

ρ_{q} (ξ_{q}, η_{q}) = \sum_{ℓ = 0}^{L_{q}} \sum_{n = 0}^{N_{q}} l_{ℓ}^{(L_{q})} (ξ_{q}) l_{n}^{(N_{q})} (η_{q}) ρ_{q, ℓ, n}

(30)

where we have added the subscript q to the dependent variable ρ to emphasize our splitting of our domain into subdomains. Here,

l_{ℓ}^{(L_{q})} (ξ_{q})

and

l_{n}^{(N_{q})} (η_{q})

are the

ℓ

-th and n-th terms of L_q-th and N_q-th order Lagrange-basis polynomials in ξ_q and η_q, respectively, as per

l_{ℓ}^{(L_{q})} (ξ_{q}) = \prod_{k = 0, k \neq ℓ}^{L_{q}} \frac{ξ_{q} - ξ_{q . k}}{ξ_{q, ℓ} - ξ_{q, k}} for k = 0, 1, 2, \dots, L_{q}

(31)

l_{n}^{(N_{q})} (η) = \prod_{k = 0, k \neq n}^{N_{q}} \frac{η_{q} - η_{q, k}}{η_{q . n} - η_{q, k}} for k = 0, 1, 2, \dots, N_{q}

(32)

and

ρ_{q . l, n}

are unknown values of the dissolved gas concentration at predetermined interpolation points,

ξ_{q, ℓ}

and

η_{q . n}

⁠. These are chosen as the Chebyshev points as per

ξ_{q, ℓ} = - cos (\frac{π}{L} ℓ), for ℓ = 0, 1, 2, \dots, L_{q}

(33)

η_{q, n} = - cos (\frac{π}{N} n), for n = 0, 1, 2, \dots, N_{q}

(34)

such that they vary from –1 to 1 and they are collocated with points in the physical domain as per the inverses of the maps given by Eqs. and . The uneven spacings of the Chebyshev points reduce Runge's phenomenon at the edges of the computational domain [22]. We note that

l_{ℓ}^{(L_{q})} (ξ_{q})

and

l_{n}^{(N_{q})} (η_{q})

equal 1 where

ξ_{q} = ξ_{q, ℓ}

and

η_{q} = η_{q, n}

⁠, respectively, and 0 on all other Chebyshev points. The dissolved gas concentrations at the Chebyshev nodes in each subdomain are expressed in the form of a

(L_{q} + 1) (N_{q} + 1)

column vector as per

ρ_{q} = {[ρ_{0, 0}, \dots, ρ_{L_{q}, 0}, ρ_{0, 1}, \dots, ρ_{L_{q}, 1}, \dots, ρ_{0, N_{q}}, \dots, ρ_{L_{q}, N_{q}}]}^{T}

(35)

Given Eq. , the column vectors of first-order partial derivatives of the dissolved gas concentration field are [23]

\frac{\partial ρ_{q}}{\partial x} = [I_{y} \otimes D_{x}] ρ_{q}

(36)

\frac{\partial ρ_{q}}{\partial y} = [D_{y} \otimes I_{x}] ρ_{q}

(37)

Here, I_y and I_x are the

(N + 1) \times (N + 1)

and

(L + 1) \times (L + 1)

identity matrices, respectively, D_x is an

(L + 1) \times (L + 1)

Chebyshev differentiation matrix as per

D_{i, j} = \frac{{\bar{c}}_{i}}{{\bar{c}}_{j}} \frac{{(- 1)}^{i + j}}{x_{i} - x_{j}}, i \neq j

(38)

D_{i, i} = - \frac{x_{i}}{2 (1 - x_{i}^{2})}, 1 \leq i \leq L - 1

(39)

D_{0, 0} = - \frac{2 L^{2} + 1}{6}

(40)

D_{L, L} = \frac{2 L^{2} + 1}{6}

(41)

where the coefficients

{\bar{c}}_{i}

and

{\bar{c}}_{j}

are

{\bar{c}}_{ℓ} = {\begin{matrix} 2 & ℓ = 0 or ℓ = L \\ 1 & 1 \leq ℓ \leq L - 1 \end{matrix}

(42)

and D_y is the analogous

(N + 1) \times (N + 1)

Chebyshev differentiation matrix.⁴ The symbol ⊗ denotes the tensor product, which produces matrices of the requisite

(L + 1) (N + 1) \times (L + 1) (N + 1)

dimensions. The second-order partial derivatives are

\frac{\partial^{2} ρ_{q}}{\partial x^{2}} = {[I_{y} \otimes D_{x}]}^{2} ρ_{q}

(43)

\frac{\partial^{2} ρ_{q}}{\partial y^{2}} = {[D_{y} \otimes I_{x}]}^{2} ρ_{q}

(44)

3.3 Time Discretization.

We discretized the rate of change of dissolved gas concentration at the Chebyshev nodes using a 3-point formula. Then, we compute the dissolved gas concentration using an implicit scheme. On a given k-th time-step, we march forward with a constant time-step size, δt. For constant δt,

{(ρ_{t})}_{ξ, η}

in Eq. is given by

{[\frac{\partial ρ^{(k + 1)}}{\partial t}]}_{ξ, η} = \frac{3 ρ^{(k + 1)} - 4 ρ^{(k)} + ρ^{(k - 1)}}{2 δ t} + O (δ t^{2})

(45)

where the superscripts

(k + 1)

⁠, (k), and

(k - 1)

denote the number of time steps elapsed since t = 0. The same formula is used to evaluate

d R / d t

in subdomains 1 and 2Q + 1 and

d h / d t

in subdomains Q and 2Q. The analogous two-point formula is used for the first time-step. As the contact angle approaches the advancing contact angle before critical time, we reduce the time-step significantly for better precision. This reduces the three-point formula to the form,

\begin{matrix} \frac{\partial ρ^{(k + 1)}}{\partial t} = {[δ t (\frac{δ t}{δ t^{'}} + 1)]}^{- 1} {(2 + \frac{δ t}{δ t^{'}}) \frac{δ t}{δ t^{'}} ρ^{(k + 1)} \\ - {(1 + \frac{δ t}{δ t^{'}})}^{2} ρ^{(k)} + ρ^{(k - 1)}} + O [{(δ t)}^{2}] \end{matrix}

(46)

where $δ t^{'}$ is the new time-step size and $d R / d t$ and $d h / d t$ are computed in a similar manner. We note that Eq. (46) simplifies to Eq. (45) for $δ t' = δ t$ ⁠. The first-time-step after critical time was treated with a two-point time discretization before reverting to a 3-point formula.

3.4 Transformed Diffusion Equation.

The transformed diffusion equation is expressed in matrix form for each subdomain. Changing the x's and y's to ξ's and η's, respectively, in Eqs. , , , and , provides

ρ_{ξ}

⁠,

ρ_{η}, ρ_{ξ ξ}

⁠, and

ρ_{η η}

⁠, respectively. It then follows from Eqs. – that

\frac{\partial ρ^{(k + 1)}}{\partial x} = [ξ_{x} (I_{η} \otimes D_{ξ}) + η_{x} (D_{η} \otimes I_{ξ})] ρ^{(k + 1)}

(47)

\frac{\partial ρ^{(k + 1)}}{\partial y} = [η_{y} (D_{η} \otimes I_{ξ})] ρ^{(k + 1)}

(48)

\frac{\partial^{2} ρ^{(k + 1)}}{\partial x^{2}} = {{[ξ_{x} (I_{n} \otimes D_{ξ}) + η_{x} (D_{η} \otimes I_{ξ})]}^{2} + η_{x x} (D_{η} \otimes I_{ξ})} ρ^{(k + 1)}

(49)

\frac{\partial^{2} ρ^{(k + 1)}}{\partial y^{2}} = {[η_{y} (D_{η} \otimes I_{ξ})]}^{2} ρ^{(k + 1)}

(50)

where the superscript

(k + 1)

denotes the dissolved gas concentration in the liquid at time

t = (k + 1) δ t

and the various partial derivatives of η and ξ are

(L + 1) (N + 1) \times (L + 1) (N + 1)

diagonal matrices. Expressing the diffusion equation in the form

A_{q} ρ_{q}^{(k + 1)} = F_{q}

(51)

where

A_{q}

and

F_{q}

are an

(L + 1) (N + 1) \times (L + 1) (N + 1)

matrix and an

(L + 1) (N + 1) \times 1

vector, respectively, it becomes

\begin{matrix} (D_{gl} {{[ξ_{x} (I_{η} \otimes D_{ξ}) + η_{x} (D_{η} \otimes I_{ξ})]}^{2} + η_{x x} (D_{η} \otimes I_{ξ}) + {[η_{y} (D_{η} \otimes I_{ξ})]}^{2}} + \\ η_{y} y_{t} (D_{η} \otimes I_{ξ}) - \frac{3}{2 δ t} I_{q}) ρ^{(k + 1)} = \frac{ρ^{(k - 1)} - 4 ρ^{(k)}}{2 δ t} \end{matrix}

(52)

for the case of three-point time discretization. Here, I_q is an $(L + 1) (N + 1) \times (L + 1) (N + 1)$ identity matrix and y_t has the same dimensions. We note that, in certain subdomains, $η_{x}$ ⁠, η_xx, $η_{y}$ ⁠, and $y_{t}$ 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.

Rows in matrices

A_{1} - A_{2 Q + 1}

and the corresponding rows in

F_{1} - F_{2 Q + 1}

corresponding to internal nodes are assembled according to Eq. . We apply the necessary boundary conditions by modifying the

2 (L_{q} + N_{q})

rows in them corresponding to external boundaries or internal ones where subdomains intersect and, denoting the revised matrices by

{\tilde{A}}_{1} - {\tilde{A}}_{2 Q + 1}

⁠, and by generating the matrices

{\tilde{A}}_{i, j}

and

{\tilde{A}}_{j, i}

⁠. The latter are sparse matrices necessary to impose the continuity of dissolved gas concentration or flux at the Chebyshev nodes located on an Ω_i-Ω_j interface. (If Ω_i and Ω_j do not share an interface, then

{\tilde{A}}_{i, j} = 0

⁠.) The resulting linear system of equations is

[\begin{matrix} {\tilde{A}}_{1} & \dots & {\tilde{A}}_{1, q} & \dots & {\tilde{A}}_{1, 2 Q + 1} \\ ⋮ & ⋱ & ⋮ \\ {\tilde{A}}_{q, 1} & {\tilde{A}}_{q} & {\tilde{A}}_{q, 2 Q + 1} \\ ⋮ & ⋱ & ⋮ \\ {\tilde{A}}_{2 Q + 1, 1} & \dots & {\tilde{A}}_{2 Q + 1, q} & \dots & {\tilde{A}}_{2 Q + 1} \end{matrix}] [\begin{matrix} ρ_{1}^{(k + 1)} \\ ⋮ \\ ρ_{q}^{(k + 1)} \\ ⋮ \\ ρ_{2 Q + 1}^{(k + 1)} \end{matrix}] = [\begin{matrix} {\tilde{F}}_{1} \\ ⋮ \\ {\tilde{F}}_{q} \\ ⋮ \\ {\tilde{F}}_{2 Q + 1} \end{matrix}]

(53)

In discussing this system of equations we refer to the matrix operator as $A_{global}$ ⁠, the full solution vector as $ρ_{global}$ ⁠, and the forcing vector as $F_{global}$ ⁠. Examples of how boundary nodes are accommodated are as follows. First, to prescribe an (external) homogeneous Neumann condition in x on the $i^{th}$ (noncorner) Chebyshev node in the column vector $ρ_{q}$ ⁠, we set the $i^{th}$ row of $A_{q}$ equal to the right-side of Eq. (47) and the $i^{th}$ row of $F_{q}$ 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 $A_{1}, A_{2}, F_{1}$ ⁠, and $F_{2}$ and utilize the relevant rows in ${\tilde{A}}_{1, 2}$ and ${\tilde{A}}_{2, 1}$ 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. [22] and Trefethen [24].

3.6 Sumdomain Types.

The subdomains fall into 7 groups as follows.

The left, bottom subdomain when $0 < t \leq t_{cr}$ in case study I and $0^{+} < t \leq t_{cr}$ in case study II bounded by $0 < x < a$ and $y_{m} < y < y_{1}$ ⁠, i.e., $Ω_{1}$ as per Fig. 2(a). After critical time, it subdivides into $Ω_{1}$ of group (B) and $Ω_{2 Q + 1}$ of group (G) as per Fig. 2(b).
Middle subdomains on $0 < x < a$ ⁠, i.e., from $Ω_{2}$ to $Ω_{Q - 1}$ and, additionally, $Ω_{1}$ for $t > t_{cr}$ ⁠.
The top subdomain on $0 < x < a$ ⁠, i.e., Ω_Q.
The bottom subdomain on a < x < d, i.e., $Ω_{Q + 1}$ ⁠.
Middle subdomains on a < x < d, from $Ω_{Q + 2}$ to $Ω_{2 Q - 1}$ ⁠.
The top subdomain on a < x < d, i.e., $Ω_{2 Q}$ ⁠.
The bottom subdomain on $0 < x < a$ after critical time, $Ω_{2 Q + 1}$ ⁠. Before critical time, the portion of the solution vector corresponding to $Ω_{2 Q + 1}$ 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, $ξ (x)$ and $η (x, y, t)$ ⁠, 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 Solution

3.7.1 Solution Up to Critical Time.

We iteratively solve for the pressure of the gas in the cavity and the dissolved gas concentration in the liquid. At a given

k^{th}

time-step we numerically solve for an improved estimate of the pressure of gas in the cavity,

p_{g, c} (t_{k, i + 1})

⁠, as per

0 = f_{1} [p_{g, c} (t_{k, i + 1})] + f_{2} [p_{g, c} (t_{k, i}), ρ_{1} (t_{k, i})]

(54)

where f₁ is the first term in Eq. (evaluated at

s (t) = s_{0}

⁠) and f₂ is the second one. An inner loop, denoted by i is used to update the maps

η_{x}, η_{x}

⁠, and

y_{t}

⁠. The first iteration of the inner loop is the solution of the linear system of equations where we assume no motion at the meniscus since the last time-step and is denoted by

ρ_{1} (t_{k, 1})

⁠. Subsequently, the solution

ρ_{1} (t_{k, i = 1})

is then substituted into Eq. to update the pressure of gas in the cavity, as denoted by

p_{g, c} (t_{k = 1, i = 2})

⁠. We compute the integral in f₂ in Eq. from

t_{k - 1}

to t_k using the Clenshaw–Curtis quadrature method as it preserves accuracy of the Chebyshev spectral method [22]. With this improved estimate of

p_{g, c} (t_{k, i + 1})

⁠, we update the maps and the corresponding terms

η_{x}, η_{y}

⁠, and

y_{t}

in Eq. for use in the next iteration. The boundary condition on the meniscus is updated to

p_{g, c} (t_{k, i = 2}) H

⁠. In the next iteration, i = 2,

p_{g, c} (t_{k, i = 2})

is the initial condition to compute a more accurate dissolved gas concentration field in the liquid,

ρ_{q} (t_{k, i = 3})

⁠. Finally, we iterate until a convergence criteria is met

| \frac{p_{g, c} (t_{k, i + 1}) - p_{g, c} (t_{k, i})}{p_{g, c} (t_{k, i + 1})} | \leq ϵ

(55)

where ϵ is user-defined.

The algorithm is:

Determine the ${\tilde{A}}_{q}$ and ${\tilde{F}}_{q}$ matrices using solution of the previous iteration as an initial condition.
Numerically solve for $ρ_{num} (t_{k, i})$ and $ρ_{sing}$ as per Eq. (B9).
Compute $p_{g, c} (t_{k, i + 1})$ as per Eqs. (54).
Update the geometric parameters, $η_{x}, η_{y}$ ⁠, and $y_{t}$ ⁠. 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 $(k + 1)$ ⁠. If the convergence criterion is not met, proceed to the next iteration (i + 1) for the current time step (k).

At

t = 0^{+}

⁠, the shape of the domain and dissolved gas concentration in the liquid are given by the initial conditions. At the end of the first time-step and first (inner) iteration,

t = δ t

and i = 1, we approximate the shape of the liquid domain and the boundary condition at the meniscus to be the same as those at

t = 0^{+}

p_{g, c} (t_{k = 1, i = 1}) \approx {\begin{matrix} p_{\infty}, & for case study I \\ p_{g, c} (0^{+}), & for case study II \end{matrix}

(56)

When critical time has been reached, we fix the gas pressure and allow the meniscus to depin.

3.7.2 After Critical Time.

After critical time, the gas pressure and radius of curvature of the meniscus are constant; therefore,

p_{g, c} (t \geq t_{cr}) = p_{cr}, θ (t \geq t_{cr}) = θ_{A}

and

R (t \geq t_{cr}) = R_{cr}

⁠. However, the meniscus depins such that

s (t)

is no longer constant and drives the variation of

y_{m} (x, t)

in the

y_{2 Q + 1}

to

η_{2 Q + 1}

map. At the first time-step after critical time, t_k

(for t_{cr} < t_{k} \leq t_{cr} + δ t)

⁠, we use an initial guess for the finite displacement of the triple contact line

s (t_{k, i = 1})

⁠. We perform a mass balance on the gas cavity in a similar manner to Eq. to obtain

0 = f_{1} [s (t_{k, i + 1})] + f_{2} [s (t_{k, i}), ρ_{1} (t_{k, i})]

(57)

For subsequent time steps, we employ the same approach as described in the previous section, except we compute $s (t_{k, i + 1})$ instead of $p_{g, c} (t_{k, i + 1})$ ⁠.

3.8 Longevity.

We assume that the transition to the Wenzel state occurs when the meniscus reaches the bottom of the gas cavity as depicted in Fig. 1, i.e.,

y_{m} (0, t_{f}) = - s_{0}

⁠, where

t_{f}

is the final time or longevity. In a limiting case of

ϕ_{s} \to 0

and

θ_{A} \to 90 deg

⁠, the solution to

t_{f}

approaches that of 1D model described in Kadoko et al. [16]. Therefore, in the 1D model, mass conservation of the gas in the cavity yields

\begin{matrix} 0 = \frac{p_{\infty} s_{0} M_{g}}{R_{u} T} - ρ_{L} g H \cdot D_{gl} t_{f} \\ + \frac{2 h}{π^{2}} \sum_{n = 1}^{\infty} (\frac{1}{n^{2}} {ρ_{L} g h H + (p_{\infty} H - ρ_{g, 0}) [1 - cos (π n)]} \\ [exp (- β_{n}^{2} D_{gl} t_{f}) - 1]) \end{matrix}

(58)

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 $H = 1.82 \times 10^{- 7} kg / (m^{3} Pa)$ [25], $D_{gl} = 2.01 \times 10^{- 9} m^{2} / s$ [26], $θ_{A} = 110 deg$ [17] and $σ = 0.073 N / m$ [27] at ambient conditions of $T = 25 ° C$ and $p_{\infty} = 1 atm$ ⁠. 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, $110 deg$ is the advancing contact angle of water on silicon treated with a flouroploymer [17]. Unless noted, we set the geometric parameters to be $s_{0} = 10 μ m$ and $a = 5 μ m$ with the liquid initially degassed. To generate the computational domain for case study $I$ ⁠, we partitioned the liquid domain into 4 subdomains (before critical time), i.e., Q = 2, where the interfaces y_q are located on ${0, 3 a, h}$ ⁠. 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 y_q are located on ${0, 3 a, 3 a + 3 (h_{0} - 3 a) / 4, h}$ ⁠. 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 $δ t = 10^{- 3} s$ was initially used before critical time. However, as the angle of the deforming meniscus neared the advancing contact angle it was changed to $δ t = 10^{- 4} s$ to more accurately capture $t_{cr}$ ⁠. A time-step of $δ t = 0.01 s$ after critical time was needed to accurately capture $t_{f}$ ⁠. A mesh dependence study was used to calculate the necessary spatial grid size. The total number of nodes was doubled until $t_{cr}$ 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 $II$ ⁠. For both cases $t_{f}$ varied by less than 0.25 percent between the final two iterations. The difference in grid construction between case study I and case study $II$ was in order to better capture the sharp gradients at the top of the domain that occur in case study $II$ ⁠.

4.1 Case Study I.

To examine case study $I$ we ran a series of computations that spanned the parameter space of solid fraction (i.e., $ϕ = (d - a) / d$ varied between 0 and 1) over a large range of values of the initial height of the liquid column (h₀). Figure 3 plots critical time versus solid fraction when h₀, 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, h₀, for these specific parameters because the time it takes for gas to diffuse through the entire domain is significantly larger than $t_{cr}$ ⁠. Also plotted is critical time computed from the quasi-one-dimensional model in our earlier publication [16]. This model is exact in the limit as $ϕ \to 0$ and $θ_{a} \to 90 deg$ 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 $110 deg$ leads to an underestimation of the critical time by 6% as $ϕ \to 0$ ⁠.

Fig. 3

$Plot of critical time, tcr, versus solid fraction, ϕ, for case study I. The total height of the liquid domain is h0=800 μm, the gas cavity half width is a=5 μm and the cavity depth is s0=10 μm.$

View large Download slide

Plot of critical time, $t_{cr}$ ⁠, versus solid fraction, $ϕ$ ⁠, for case study $I$ ⁠. The total height of the liquid domain is $h_{0} = 800 μ m$ ⁠, the gas cavity half width is $a = 5 μ m$ and the cavity depth is $s_{0} = 10 μ m$ ⁠.

Fig. 4

$Plot of longevity versus initial height of the liquid table, h0, for case study I for various solid fractions, ϕ. The gas cavity half width is a=5 μm and the cavity depth is s0=10 μm. The vertical dashed lines indicate the minimum height in which full Cassie to Wensel state transition can occur.$

View large Download slide

Plot of longevity versus initial height of the liquid table, h₀, for case study $I$ for various solid fractions, $ϕ$ ⁠. The gas cavity half width is $a = 5 μ m$ and the cavity depth is $s_{0} = 10 μ m$ ⁠. The vertical dashed lines indicate the minimum height in which full Cassie to Wensel state transition can occur.

Figure 4 plots longevity (⁠

t_{f}

⁠) versus h₀ for solid fractions ranging from 0.1 to 0.9. Again, it is evident that increasing solid fraction speeds up the diffusive process, this time exhibited by the decrease in

t_{f}

⁠. This plot quantifies the increase in

t_{f}

with decreasing h₀. The (vertical) dashed lines correspond to the limit at which a full Cassie to Wenzel state transition can occur. They represent the minimum liquid column height that has the capacity to hold

m'_{0} - m'_{f}

mass of gas in a saturated solution, where

m'_{f}

is the mass in the cavity at time

t_{f}

⁠, i.e.

h_{min} = \frac{(1 - ϕ) (m'_{0} - m'_{f})}{a (p_{\infty} H - ρ_{g, 0})}

(59)

which is a modification to a result found in our earlier publication, i.e., Kadoko et al. [16]. The computed data is again compared to that of Kadoko et al. [16] 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. $V_{min}$ is the volume of gas in the cavity when the lowest part of the meniscus makes contact with the cavity bottom.

Fig. 5

$Plot of volume of gas in the cavity versus time for various solid fractions as per case study I. The total height of the liquid domain is h=700 μm, the gas cavity half width is a=5 μm and the cavity depth is s0=10 μm. The inset focuses on the time before critical time.$

View large Download slide

Plot of volume of gas in the cavity versus time for various solid fractions as per case study $I$ ⁠. The total height of the liquid domain is $h = 700 μ m$ ⁠, the gas cavity half width is $a = 5 μ m$ and the cavity depth is $s_{0} = 10 μ m$ ⁠. The inset focuses on the time before critical time.

Figure 6 shows contour plots of the dissolved gas concentration field in the lower region of the liquid domain, $- 10 μ m < y < 15 μ m$ ⁠, with $ϕ = 0.5$ ⁠. Figure 6(a) shows ρ at $t = 0.05 s$ ⁠, 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 $t = 25 s$ and $t = t_{f}$ ⁠, respectively. At this point the concentration gradient has decreased significantly and diffusion rates are low. This helps illuminate why $t_{f}$ is much larger than $t_{c}$ ⁠.

Fig. 6

View large Download slide

Contour plot of the dissolved gas concentration in the liquid for case study $I$ at (a) $t = 0.05 s$ ⁠, (b) $t = 0.539 s$ (⁠ $t_{cr}$ ⁠), (c) $25 s$ ⁠, and (d) $t = 38.34 s$ (⁠ $t_{f}$ ⁠) showing a portion of the liquid domain on interval $- 10 μ m < y < 50 μ m$ for initially degassed water in case study $I$ ⁠. The total height of the liquid domain is $h = 800 μ m$ ⁠, the gas cavity half width is $a = 5 μ m, ϕ = 0.5$ ⁠, and the cavity depth is $s_{0} = 10 μ m$ ⁠. The sharp decrease in the magnitude of concentration gradient between a) and d) show that diffusion has slowed greatly by the time $t_{f}$ is reached.

4.2 Case Study II.

Figure 7 plots longevity as a function of cavity depth for a characteristic initial liquid height, $h_{0} = 3000 μ m$ and cavity half-width of $a = 5 μ m$ ⁠. The first thing to note is that a larger initial column height of $h_{0} = 3000 μ m$ in case study $II$ leads to less longevity, $t_{f}$ ⁠, than in case study $I$ for $h_{0} = 1500 μ m$ and the same cavity depth, $s_{0} = 10 μ m$ ⁠, as per Fig. 4. This is because in case study $II$ 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 $II$ from Kadoko et al. [16]. In this case, it fulfills the solution requirement for $ϕ = 0$ ⁠, demonstrating much the same shape as our results.

Fig. 7

View large Download slide

Plot of cavity depth versus longevity as per case study $II$ for $a = 5 μ m$ and $h_{0} = 3000 μ m$ ⁠. Inset: Plot of cavity depth versus critical time as per case study $II$ for $a = 5 μ m$ and $h_{0} = 3000 μ m$ ⁠.

Figure 8 plots longevity for varying initial gas concentration for a characteristic initial liquid height, $h_{0} = 3000 μ m$ ⁠. 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. [16] for $thet a_{a} = 90 deg$ are included.

Fig. 8

View large Download slide

Plot of initial dissolved gas concentration in the liquid versus longevity as per case study II for $a = 5 μ m$ and $h_{0} = 3000 μ m$ ⁠. The cavity depth is $s_{0} = 10 μ m$ ⁠.

5 Conclusion

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. [3]. 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.

Acknowledgment

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.

Funding Data

Google Award in Faculty Research.
National Science Foundation under (Grant No. 1402783; Funder ID: 10.13039/100000001).

Footnotes

2

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) [9]. A complete Wenzel state is achieved when the nano-structures are wetted.

3

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.

4

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).

The impermeability boundary conditions and Henry's Law, i.e., Eqs. and , respectively, apply on exterior faces and on the

Ω_{1}

-

Ω_{Q + 1}

and

Ω_{1}

-

Ω_{2}

interfaces, we match fluxes. Thus

\frac{\partial}{\partial x} ρ_{1} (0, y, t) = 0

(A1)

ρ_{1} [x, y_{m} (x, t)] = p_{g, c} (t) H

(A2)

\frac{\partial}{\partial x} ρ_{1} (a, y, t) = \frac{\partial}{\partial x} ρ_{Q + 1} (a, y, t)

(A3)

\frac{\partial}{\partial y} ρ_{1} (x, y_{1}, t) = \frac{\partial}{\partial y} ρ_{2} (x, y_{1}, t)

(A4)

We transform these conditions to be in terms of ξ and η as per Eqs. and using the maps

ξ_{q} (x) = 2 \frac{x}{a} - 1, for 0 < x < a

(A5)

η_{1} (x, y, t) = 2 \frac{y - y_{m} (x, t)}{y_{1} - y_{m} (x, t)} - 1, for t < t_{cr}

(A6)

where

ξ_{q} (x)

is valid for groups (A), (B), (C), and (G) and

η_{1} (x, y, t)

is only valid for the interval

0 < t \leq t_{cr}

because subdomain

Ω_{1}

splits into an

a \times y_{1}

rectangle (also denoted by

Ω_{1}

and belonging to group (B)) and

Ω_{2 Q + 1}

after it. The derivative of

ξ_{q} (x)

is given by

\frac{d ξ_{q}}{d x} = \frac{2}{a}, for 0 < x < a

(A7)

The partial derivatives of

η_{1} (x, y, t)

with respect to x and y are

\frac{\partial η_{1}}{\partial x} = \frac{η - 1}{y_{1} - y_{m} (x, t)} \frac{x}{\sqrt{R {(t)}^{2} - x^{2}}}

(A8)

\frac{\partial η_{1}}{\partial y} = \frac{2}{y_{1} - y_{m} (x, t)}

(A9)

Also, required, for groups (A) and (G), is the partial derivative of y with respect to time, which follows from Eqs. and as

{(\frac{\partial y_{(A, G)}}{\partial t})}_{ξ, η} = \frac{1 - η}{2} {[\frac{R (t)}{\sqrt{R {(t)}^{2} - a^{2}}} - \frac{R (t)}{\sqrt{R {(t)}^{2} - x^{2}}}] \frac{d R}{d t} - \frac{d s}{d t}}, for 0 < t \leq t_{cr}

(A10)

A.2 Subdomain Type (B).

The boundary conditions for subdomains in group (B) consist of symmetry on x = 0 (for nodes solely in one subdomain) and continuity of dissolved gas concentration and flux of it on the interfaces share by two subdomains. Thus

\frac{\partial}{\partial x} ρ_{q} (0, y, t) = 0

(A11)

\frac{\partial}{\partial x} ρ_{q} (a, y, t) = \frac{\partial}{\partial x} ρ_{q + Q} (a, y, t)

(A12)

ρ_{q} (x, y_{q - 1}, t) = ρ_{q - 1} (x, y_{q - 1}, t)

(A13)

\frac{\partial}{\partial y} ρ_{q} (x, y_{q}, t) = \frac{\partial}{\partial y} ρ_{q + 1} (x, y_{q}, t)

(A14)

where

q = {\begin{matrix} 2, 3, \dots, Q - 1, & for 0 < t \leq t_{cr} \\ 1, 2, \dots, Q - 1, & for t_{cr} < t \leq t_{f} \end{matrix}

(A15)

The vertical map and its derivative with respect to y are

η_{q} (y) = 2 \frac{y - y_{q - 1}}{y_{q} - y_{q - 1}} - 1

(A16)

\frac{\partial η_{q}}{\partial y} = \frac{2}{y_{q} - y_{q - 1}}

(A17)

As for the Q – 1 corner points that are shared by four subdomains each, we prescribe three concentration and one flux continuity condition as per

ρ_{q} (a, y_{q}, t) = ρ_{q + 1} (a, y_{q}, t)

(A18)

ρ_{q + 1} (a, y_{q}, t) = ρ_{q + Q + 1} (a, y_{q}, t)

(A19)

ρ_{q + Q + 1} (a, y_{q}, t) = ρ_{q + Q} (a, y_{q}, t)

(A20)

\frac{\partial}{\partial x} ρ_{q + Q} (a, y_{q}, t) = \frac{\partial}{\partial x} ρ_{q} (a, y_{q}, t)

(A21)

where $q = 1, 2, \dots, Q - 1$ ⁠.

A.3 Subdomain Type (C).

Subdomain type (C), i.e., Ω_Q, is the topmost one on

0 < x < a

⁠. At the top side of it we apply impermeability and constant dissolved gas concentration boundary conditions for case studies I and II, respectively. On the bottom and right sides, we apply dissolved gas concentration and flux continuity, respectively. At the left side we apply symmetry conditions. Thus

\frac{\partial}{\partial x} ρ_{Q} (0, y, t) = 0

(A22)

\frac{\partial}{\partial x} ρ_{Q} (a, y, t) = \frac{\partial}{\partial x} ρ_{q + Q} (a, y, t)

(A23)

ρ_{Q} (x, y_{Q - 1}, t) = ρ_{Q - 1} (x, y_{Q - 1}, t)

(A24)

\frac{\partial}{\partial y} ρ_{Q} (x, h (t), t) = 0, for case study I

(A25)

ρ_{Q} (x, h (t), t) = p_{\infty} H, for case study II

(A26)

The vertical map and its derivative are

η_{Q} (y, t) = 2 \frac{y - y_{Q - 1}}{h (t) - y_{Q - 1}} - 1

(A27)

\frac{\partial η_{Q}}{\partial y} = \frac{2}{h (t) - y_{Q - 1}}

(A28)

For groups (C) and (F), the rate of change of y with respect to t is

{(\frac{\partial y_{(C, F)}}{\partial t})}_{ξ, η} = \frac{1}{2} (1 + η) \frac{d h}{d t}

(A29)

A.4 Subdomain Type (D).

Subdomain type (D), i.e.,

Ω_{Q + 1}

⁠, is directly above the ridge. The symmetry and impermeability boundary conditions apply to the exterior. We apply dissolved gas concentration and flux continuity on the left and top interfaces. Thus

ρ_{Q + 1} (a, y, t) = ρ_{1} (a, y, t)

(A30)

\frac{\partial}{\partial x} ρ_{Q + 1} (d, y, t) = 0

(A31)

\frac{\partial}{\partial y} ρ_{Q + 1} (x, 0, t) = 0

(A32)

\frac{\partial}{\partial y} ρ_{Q + 1} (x, y_{1}, t) = \frac{\partial}{\partial y} ρ_{Q + 2} (x, y_{1}, t)

(A33)

The horizontal and vertical maps are

ξ_{q} (x) = 2 \frac{x - a}{d - a} - 1, for a < x < d

(A34)

η_{Q + 1} (y, t) = 2 \frac{y}{y_{1}} - 1

(A35)

where Eq. is valid for subdomain groups (D), (E) and (F). The derivative of

ξ_{q}

with respect to x is

\frac{\partial ξ_{q}}{\partial x} = \frac{2}{d - a}

(A36)

and the derivative of the vertical map with respect to y is

\frac{\partial η_{Q + 1}}{\partial y} = \frac{2}{y_{1}}

(A37)

A.5 Subdomain Type (E).

Similar to subdomain group (B), boundary conditions for group (E) consist of a symmetry condition on one exterior side (x = d) and continuity on three internal interfaces. Thus

ρ_{q} (a, y, t) = ρ_{q - Q} (a, y, t)

(A38)

\frac{\partial}{\partial x} ρ_{q} (d, y, t) = 0

(A39)

ρ_{q} (x, y_{q - Q - 1}, t) = ρ_{q - 1} (x, y_{q - Q - 1}, t)

(A40)

\frac{\partial}{\partial y} ρ_{q} (x, y_{q - Q}, t) = \frac{\partial}{\partial y} ρ_{q + 1} (x, y_{q - Q}, t)

(A41)

where

q = Q + 2, Q + 3, \dots, 2 Q - 1

(A42)

The map used for subdomains of type (E) is

η_{q} (y) = 2 \frac{y - y_{q - Q - 1}}{y_{q - Q} - y_{q - Q - 1}} - 1

(A43)

where

\frac{\partial η_{q} (y)}{\partial y} = \frac{2}{y_{q - Q} - y_{q - Q - 1}}

(A44)

A.6 Subdomain Type (F).

Subdomain group (F) is the topmost portion of the liquid domain on a < x < d. The boundary condition at the top of

Ω_{2 Q}

is the same as that on Ω_Q as per Eqs. and . On the left and bottom sides, we prescribe dissolved gas concentration and flux continuity as per

ρ_{2 Q} (a, y, t) = ρ_{Q} (a, y, t)

(A45)

\frac{\partial}{\partial x} ρ_{2 Q} (d, y, t) = 0

(A46)

\frac{\partial}{\partial y} ρ_{2 Q} (x, y_{Q - 1}, t) = \frac{\partial}{\partial y} ρ_{2 Q - 1} (x, y_{Q - 1}, t)

(A47)

\frac{\partial}{\partial y} ρ_{2 Q} (x, h (t), t) = 0, for case study I

(A48)

ρ_{2 Q} (x, h (t), t) = p_{\infty} H, for case study II

(A49)

The map in the vertical direction and its derivatives are the same as those given in Eqs. (A27)–(A28).

A.7 Subdomain Type (G).

Subdomain type (G), i.e.,

Ω_{2 Q + 1}

exists after critical time as shown in Fig. 2(b). The boundary conditions are

\frac{\partial}{\partial x} ρ_{2 Q + 1} (0, y, t) = 0

(A50)

\frac{\partial}{\partial x} ρ_{2 Q + 1} (a, y, t) = 0

(A51)

ρ_{2 Q + 1} (x, y_{m} (x, t), t) = p_{cr} H

(A52)

\frac{\partial}{\partial y} ρ_{2 Q + 1} (x, 0, t) = \frac{\partial}{\partial y} ρ_{1} (x, 0, t)

(A53)

The vertical mapping function is

η_{2 Q + 1} (x, y, t) = 1 - 2 \frac{y}{y_{m} (x, t)}

(A54)

The first-order partial derivatives of

η_{2 Q + 1} (x, y, t)

are

\frac{\partial η_{2 Q + 1}}{\partial x} = \frac{1 - η}{y_{m} (x, t)} \frac{x}{\sqrt{R {(t)}^{2} - x^{2}}}

(A55)

\frac{\partial η_{2 Q + 1}}{\partial y} = - \frac{2}{y_{m} (x, t)}

(A56)

The partial derivative in y with respect to t is given Eq. (A10), where $d R / d t = 0$ ⁠, 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. [28] 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 ${\tilde{A}}_{global}$ 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 ${\tilde{A}}_{global}$ 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.

The local problem in the vicinity of the triple contact point, cast in polar (⁠

r, ϕ

⁠) coordinates to facilitate its solution, is shown in Fig. 9. Only the Laplacian terms in the diffusion equation are singular; therefore, we need only resolve Laplace's equation in our local analysis. We decompose our solution into a numerically-resolved component,

ρ_{num}

⁠, and an analytically-resolved singular one as per

ρ_{global} = ρ_{num} + ρ_{sing}

(B1)

Fig. 9

View large Download slide

Local problem before critical time

We impose the in homogemenous boundary condition along the meniscus as per

ρ_{num} = H p_{g, c}

⁠. Then, both boundary conditions on the local problem are homogeneous and it reads

\nabla^{2} ρ_{sing} = 0

(B2)

{\frac{\partial ρ_{sing}}{\partial ϕ} |}_{ϕ = 0} = 0

(B3)

{ρ_{sing} |}_{ϕ = θ + π / 2} = 0.

(B4)

The solution is

B r^{λ} cos (λ ϕ)

(B5)

where B is the unknown singularity strength, and the eigencondition reads

cos [λ (θ + π / 2)] = 0.

(B6)

We need only consider the strongest singularity such that

ρ_{sing} = B r^{λ_{0}} cos (λ_{0} ϕ)

(B7)

where

λ_{0} = π / (2 θ + π)

⁠. In Cartesian coordinates this result is

ρ_{sing} = B {(x^{2} + y^{2})}^{1 / (2 λ_{0})} cos [λ_{0} arctan (y / x)]

(B8)

We subtract

ρ_{sing}

from the numerical problem by recasting the overall problem as

[\begin{matrix} A_{global} & v_{sing} \\ v_{diff} & 0 \end{matrix}] [\begin{matrix} ρ_{num} \\ B \end{matrix}] = [\begin{matrix} F_{global} \\ 0 \end{matrix}]

(B9)

Here, on an interior point,

v_{sing, i}

takes the form

v_{sing, i} = y_{t} \times \frac{1}{B} {\frac{\partial ρ_{sing}}{\partial y} |}_{x = x (i), y = y (i)} - α \frac{1}{B} {ρ_{sing} |}_{x = x (i), y = y (i)}

(B10)

where α is a coefficient that depends on the finite difference scheme used. It assumes the values

α = {\frac{3}{2 δ t}, \frac{1}{δ t}, [δ t (\frac{δ t}{δ t'} + 1)] (2 + \frac{δ t}{δ t'}) \frac{δ t}{δ t'}}

(B11)

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 $ρ_{sing}$ 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.

When matching flux or concentration on adjacent subdomains singularity subtraction is not necessary and thus the appropriate element of

v_{sing}

is set to 0. Finally, for points on external boundaries,

v_{sing}

is populated with the appropriate contribution of

ρ_{sing}

⁠. For example, when the

i^{th}

row of

A_{global}

corresponds to a node operated by the symmetry boundary condition at x = 0, then

v_{sing, i} = \frac{1}{B} {\frac{\partial ρ_{sing}}{\partial x} |}_{x = 0, y = y (i)}

where i indicates the i-th index and y(i) is the y location of the i-th node. Thus, the solution, $ρ_{global} = ρ_{num} + ρ_{sing}$ satisfies the external boundary conditions.

We remove the singularity from the numerical solution, $ρ_{num}$ via the a row vector $v_{diff}$ in Eq. (B9). It applies the differential operator $\partial / \partial x$ at the (singular) triple contact point in subdomain Q + 1. The end result is a new system of $(2 Q + 1) (L + 1) (N + 1) + 1$ equations. The element 0 in the right-side vector restricts $\partial ρ_{num} / \partial x = 0$ and therefore removes the singularity from $ρ_{num}$ ⁠. 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 $x = a, y = 0$ and at the triple contact point, respectively, again cast in polar (⁠ $r, ϕ$ ⁠) 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.

Fig. 10

View large Download slide

Local problems at (a) x = a, y = 0 and (b) the triple contact point after critical time

A separation of variables solution shows that the form of the solution at x = a, y = 0 most singular in the partial derivative with respect to r is

ρ_{sing, 1} = C r^{2 / 3} cos (ν_{1} ϕ)

(B12)

where C is the unknown singularity strength. The singular problem at the triple contact point is

\nabla^{2} ρ_{sing, 2} = 0

(B13)

Subject to

{\frac{\partial ρ_{sing, 2}}{\partial ϕ} |}_{ϕ_{2} = π / 2} = 0

(B14)

{ρ_{sing, 2} |}_{ϕ_{2} = θ_{a} + π / 2} = 0.

(B15)

It follows from the separation of variables and trigonometric identities that the form of the solution at the triple contact line that is most singular in the partial derivative with respect to r is

ρ_{sing, 2} = r_{2}^{π / (2 θ_{a})} [D sin (\frac{π ϕ_{2}}{2 θ_{a}}) + E \cot (\frac{π}{2 θ_{a}}) cos (\frac{π ϕ_{2}}{2 θ_{a}})]

(B16)

where D is the unknown singularity strength.

Subtracting

ρ_{sing, 1}

and

ρ_{sing, 2}

from the numerical problem changes the overall one to

[\begin{matrix} A_{global} & v_{1, sing} & v_{2, sing} \\ v_{1, diff} & 0 & 0 \\ v_{2, diff} & 0 & 0 \end{matrix}] [\begin{matrix} ρ_{num} \\ C \\ D \end{matrix}] = [\begin{matrix} F_{global} \\ 0 \\ 0. \end{matrix}]

(B17)

The role of the column vectors $v_{1, diff}$ and $v_{2, diff}$ are analogous to that of $v_{diff}$ in Eq. (B9) before critical time. Moreover, the row vector $v_{1, diff}$ causes the operator $\partial / \partial x$ to act on the element of $ρ_{num}$ corresponding to the bottom left corner of subdomain Q + 1 and the row vector $v_{2, diff}$ causes the operator $\partial / \partial y$ to act on the element of $ρ_{num}$ corresponding to the triple contact spot in subdomain $2 Q + 1$ ⁠. Both these partial derivatives are set equal to 0. The result is a system of $(2 Q + 1) (L + 1) (N + 1) + 2$ equations. The solution process is analogous to that before critical time.

Nomenclature

Greek Symbols

Α =: coefficient for singularity removal
$ρ$ =: vector of dissolved gas concentration
δt =: time step
ϵ =: convergence criteria
η =: map to vertical of unit square
Ω =: domain
$ϕ$ =: solid fraction, $(d - a) / d$
ρ =: density (⁠ ${mol / m}^{2}$ ⁠)
σ =: surface tension along meniscus (⁠ $N / m$ ⁠)
θ =: contact angle of penetration into cavity, measured from the horizontal
ξ =: map to horizontal of unit square

Roman

$A$ =: matrix operator
${\tilde{A}}_{i, j}$ =: matrix for coupling subdomains i and j
a =: half-width of the cavity (m)
$B, C$ =: singularity strengths
D =: binary molecular diffusion coefficient (⁠ $m^{2} / s$ ⁠) or differentiation matrix
d =: half of the ridge period (m)
$F$ =: forcing vector
g =: gravitational constant (⁠ $m^{2} / s$ ⁠)
H =: Henry's constant (⁠ ${mol / m}^{3} \cdot Pa$ ⁠)
h =: transient height of the liquid domain (m)
h₀ =: initial height of the liquid domain (m)
$m^{'}$ =: mass per unit depth of gas in cavity (⁠ ${kg / m}^{2}$ ⁠)
$M_{g}$ =: molecular weight of gas (⁠ $kg / mol$ ⁠)
p =: pressure (⁠ $Pa$ ⁠)
Q =: number of subdomains in vertical direction
R =: radius of curvature of meniscus (m)
$R_{u}$ =: universal gas constant (⁠ $J / mol \cdot K$ ⁠)
s =: coordinate along arc spanning meniscus (m)
s =: transient location (in y) of the triple contact point (m)
s₀ =: depth of the cavity (m)
T =: temperature (K)
t =: time (s)
v =: contributions from singularity removal to global matrix
$V^{'}$ =: volume per unit depth of gas in cavity (⁠ $m^{2}$ ⁠)
y_m =: meniscus location (m)

Subscripts

$\infty$ =: far-field (atmosphere)
$a$ =: advancing (contact angle)
$cr$ =: critical
$diff$ =: differential operator on the triple contact point
$f$ =: final
$g, c$ =: gas in the cavity
$global$ =: referring to global matrix and vectors
$gl$ =: gas–liquid
$g$ =: gas
$l$ =: liquid
$num$ =: referring to numerical nonsingular solution
$q$ =: subdomain number
$sing$ =: referring to quasi-analytical singular solution

References

1.

Lee

,

C.

, and

Kim

,

C. J.

,

2009

, “

Maximizing the Giant Liquid Slip on Superhydrophobic Microstructures by Nanostructuring Their Sidewalls

,”

Langmuir

,

25

(

21

), pp.

12812

–

12818

2.

Bolognesi

,

G.

,

Cottin-Bizonne

,

C.

, and

Pirat

,

C.

,

2014

, “

Evidence of Slippage Breakdown for a Superhydrophobic Microchannel

,”

Phys. Fluids

,

26

(

8

), p.

082004

.10.1063/1.4892082

Google Scholar

Crossref

3.

Xu

,

M. C.

,

Sun

,

G. U.

, and

Kim

,

C. J.

,

2014

, “

Infinite Lifetime of Underwater Superhydrophobic States

,”

Phys. Rev. Lett.

,

113

(

13

), p.

136103

.10.1103/PhysRevLett.113.136103

Google Scholar

Crossref

PubMed

4.

Lafuma

,

A.

, and

Quere

,

D.

,

2003

, “

Superhydrophobic States

,”

Nat. Mater.

,

2

(

7

), pp.

457

–

60

5.

Rothstein

,

J. P.

,

2010

, “

Slip on Superhydrophobic Surfaces

,”

Annu. Rev. Fluid Mech.

,

42

(

1

), pp.

89

–

109

.10.1146/annurev-fluid-121108-145558

Google Scholar

Crossref

6.

Yahui

,

X.

,

Pengyu

,

L.

,

Hao

,

L.

, and

Huiling

,

D.

,

2016

, “

Underwater Superhydrophobicity: Stability, Design and Regulation, and Applications

,”

ASME Appl. Mech. Rev.

,

68

(

3

), p.

030803

.10.1115/1.4033706

7.

Rahn

,

H.

, and

Paganelli

,

C. V.

,

1968

, “

Gas Exchange in Gas Gills of Diving Insects

,”

Respiration Physiol.

,

5

(

1

), pp.

145

–

164

.10.1016/0034-5687(68)90083-2

Google Scholar

Crossref

8.

Lv

,

P.

,

Xue

,

Y.

,

Shi

,

Y.

,

Lin

,

H.

, and

Duan

,

H.

,

2014

, “

Metastable States and Wetting Transition of Submerged Superhydrophobic Structures

,”

Phys. Rev. Lett.

,

112

(

19

), p.

196101

.

Google Scholar

Crossref

PubMed

9.

Lee

,

C.

, and

Kim

,

C.-J.

,

2011

, “

Underwater Restoration and Retention of Gases on Superhydrophobic Surfaces for Drag Reduction

,”

Phys. Rev. Lett.

,

106

(

1

).

10.

Dilip

,

D.

,

Jha

,

N. K.

,

Govardhan

,

R. N.

, and

Bobji

,

M. S.

,

2014

, “

Controlling Air Solubility to Maintain “Cassie” State for Sustained Drag Reduction

,”

Colloids Surf. A: Physicochem. Eng. Aspects

,

459

, pp.

217

–

224

.10.1016/j.colsurfa.2014.07.006

Google Scholar

Crossref

11.

Emami

,

B.

,

Hemeda

,

A. A.

,

Amrei

,

M. M.

,

Luzar

,

A.

,

Gad-el Hak

,

M.

, and

Vahedi Tafreshi

,

H.

,

2013

, “

Predicting Longevity of Submerged Superhydrophobic Surfaces With Parallel Grooves

,”

Phys. Fluids

,

25

(

6

), p.

062108

.10.1063/1.4811830

Google Scholar

Crossref

12.

Lv

,

P.

,

Xue

,

Y.

,

Liu

,

H.

,

Shi

,

Y.

,

Xi

,

P.

,

Lin

,

H.

, and

Duan

,

H.

,

2015

, “

Symmetric and Asymmetric Meniscus Collapse in Wetting Transition on Submerged Structured Surfaces

,”

Langmuir

,

31

(

4

), pp.

1248

–

1254

13.

Hensel

,

R.

,

Finn

,

A.

,

Helbig

,

R.

,

Killge

,

S.

,

Braun

,

H.-G.

, and

Werner

,

C.

,

2014

, “

In Situ Experiments to Reveal the Role of Surface Feature Sidewalls in the Cassie–Wenzel Transition

,”

Langmuir

,

30

(

50

), pp.

15162

–

15170

14.

Søgaard

,

E.

,

Andersen

,

N. K.

,

Smistrup

,

K.

,

Larsen

,

S. T.

,

Sun

,

L.

, and

Taboryski

,

R.

,

2014

, “

Study of Transitions Between Wetting States on Microcavity Arrays by Optical Transmission Microscopy

,”

Langmuir

,

30

(

43

), pp.

12960

–

12968

15.

Dufour

,

R.

,

Saad

,

N.

,

Carlier

,

J.

,

Campistron

,

P.

,

Nassar

,

G.

,

Toubal

,

M.

,

Boukherroub

,

R.

,

Senez

,

V.

,

Nongaillard

,

B.

, and

Thomy

,

V.

,

2013

, “

Acoustic Tracking of Cassie to Wenzel Wetting Transitions

,”

Langmuir

,

29

(

43

), pp.

13129

–

13134

16.

Kadoko

,

J.

,

Karamanis

,

G.

,

Kirk

,

T.

, and

Hodes

,

M.

,

2017

, “

One-Dimensional Analysis of Gas Diffusion-Induced Cassie to Wenzel State Transition

,”

ASME J. Heat Transfer-Trans. ASME

,

139

(

12

), p.

122006

.10.1115/1.4036600

Google Scholar

Crossref

17.

Goswami

,

S.

,

Klaus

,

S.

, and

Benziger

,

J.

,

2008

, “

Wetting and Absorption of Water Drops on Nafion Films

,”

Langmuir

,

24

(

16

), pp.

8627

–

8633

.http://www.agcce.com/cytop-technical-information/

Google Scholar

Crossref

PubMed

18.

Epstein

,

P. S.

, and

Plesset

,

M. S.

,

1950

, “

On the Stability of Gas Bubbles in Liquid-Gas Solutions

,”

J. Chem. Phys.

,

18

(

11

), pp.

1505

–

1509

.10.1063/1.1747520

Google Scholar

Crossref

19.

Oberkampf

,

W. L.

,

1976

, “

Domain Mappings for the Numerical Solution of Partial Differential Equations

,”

Int. J. Numer. Methods Eng.

,

10

(

1

), pp.

211

–

223

.10.1002/nme.1620100116

Google Scholar

Crossref

20.

Furzeland

,

R. M.

,

1977

, “

The Numerical Solution of Two-Dimensional Moving Boundary Problems Using Curvilinear co-Ordinate Transformations

,” Brunel University Mathematics Technical Papers Collection.

21.

John

,

C.

,

1987

,

Free and Moving Boundary Problems

,

Clarendon Press

,

New York/Oxford, UK

22.

Trefethen

,

L. N.

,

2000

,

Spectral Methods in MATLAB

, Vol.

10

,

SIAM, Society for Industrial and Applied Mathematics

,

Philadelphia, PA

.

Google Scholar

Crossref

23.

Weidong

,

G.

,

Gérard

,

L.

, and

Ranga

,

N.

,

2012

,

The Application of the Chebyshev-Spectral Method in Transport Phenomena

,

68

,

Springer

,

New York; Berlin

.

24.

Canuto

,

C.

,

Hussaini

,

M. Y.

,

Quarteroni

,

A.

, and

Zang

,

T. A.

,

1988

,

Spectral Methods in Fluid Dynamics

,

Springler-Verlag

,

New-York

.

Google Scholar

Crossref

25.

Sander

,

R.

,

2016

,

Henry's Law Constants

,

National Institute of Standards and Technology

,

Gaithersburg, MD

, accessed June 29, 2016, http://webbook.nist.gov

26.

Ferrell

,

R. T.

, and

Himmelblau

,

D. M.

,

1967

, “

Diffusion Coefficients of Nitrogen and Oxygen in Water

,”

J. Chem. Eng. Data

,

12

(

1

), pp.

111

–

115

.10.1021/je60032a036

Google Scholar

Crossref

27.

Lemmon

,

E. W.

,

McLinden

,

M. O.

, and

Friend

,

D. G.

,

Thermophysical Properties of Fluid Systems

,

National Institute of Standards and Technology

,

Gaithersburg, MD

, accessed June 29, 2016, http://webbook.nist.gov

28.

Game

,

S.

,

Hodes

,

M.

,

Keaveny

,

E.

, and

Papageorgiou

,

D.

,

2017

, “

Physical Mechanisms Relevant to Flow Resistance in Textured Microchannels

,”

Phys. Rev. Fluids

,

2

(

9

), p.

094102

.10.1103/PhysRevFluids.2.094102

Google Scholar

Crossref

Two-Dimensional Numerical Analysis of Gas Diffusion-Induced Cassie to Wenzel State Transition

Abstract

1 Introduction

2 Mathematical Model

3 Numerical Method

3.1 Introduction.

3.2 Domain Decomposition and Mapping.

3.3 Time Discretization.

3.4 Transformed Diffusion Equation.

3.5 Assembly of Linear System of Equations.

3.6 Sumdomain Types.

3.7 Solution

3.7.1 Solution Up to Critical Time.

3.7.2 After Critical Time.

3.8 Longevity.

4 Results and Discussion

4.1 Case Study I.

4.2 Case Study II.

5 Conclusion

Acknowledgment

Funding Data

Footnotes

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).

Appendix B: Singularity Removal

B.1 Singularity Removal.

B.1.1 Before Critical Time.

B.1.2 After Critical Time.

Nomenclature

Greek Symbols

Roman

Subscripts

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

Journals

Conference Proceedings

eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only