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

Superhydrophobic surfaces (SHPos), i.e., hydrophobic surfaces that are topographically rough, may support a droplet or a liquid column in an unwetted (Cassie) state [1]. Then, in the case of a droplet, the apparent contact angle exceeds the Young angle and can approach 180 deg. It is also well known that the use of SHPos in pressure-driven flows in microchannels reduces frictional drag [2]. However, it is challenging to design and fabricate SHPo surfaces that can permanently maintain a stable Cassie state. Indeed, due to phenomenon such as evaporation [3], condensation [4], pressurization [1], impact [5], and gas diffusion [6], the Cassie state may be gradually compromised resulting in a transition to a wetted (Wenzel) state over a time period known as longevity [6,7].

Limited longevity presents a key challenge in the industrial application of SHPo surfaces. In the context of microgap cooling, structured surfaces that have high longevity at high operating pressure are desirable because they enable high channel pressure drops which result in high volumetric flow rate and, thus, low caloric resistance in comparison to planar Poiseuille flow [8]. Similarly, in underwater applications, structured surfaces that have high longevity under high hydrostatic pressure may enable vehicles to operate at high speeds with low energy consumption [9,10]. Among the methods of enhancing longevity suggested in the literature are gas injection and electrolysis of water to produce inert gas [7,10,11]. To implement such techniques, there is a need to model the physics of Cassie to Wenzel state transition as a result of gas diffusion.

Theoretical models and experimental studies show that the shape of the meniscus changes in two primary stages. Initially, the meniscus curvature increases and, thus, the contact angle (θ) increases. A stable Cassie state may be reached in which the contact angle reaches a terminal value. Otherwise, it becomes equal to the advancing contact angle (θ_A) at the critical time (t_cr). Then, the triple contact line depins and moves down into the cavity with the contact angle fixed at θ_A [12,13]. The meniscus touches the bottom of the cavity at the final time (t_f) in a symmetric or asymmetric fashion, where the latter is thought to be caused by the presence of impurities in the cavity walls [14]. A broad literature review on gas diffusion-induced Cassie to Wenzel state transition was done by Xue et al. for structures whose gas cavities are closed to the atmosphere [15]. To validate the aforementioned theory, researchers have experimentally tracked the menisci using techniques such as photography [11], light reflection [6], confocal microscopy [14], and acoustic tracking [16].

Neglecting inertial and viscous effects in the liquid and applying the ideal gas law in the gas cavity, Emami et al. predict the longevity of a submerged ridge-type SHPo surface using Eq. (1), an equilibrium force balance on the menisci and the ideal gas law [17]. Their results predict that, for a given geometry, longevity decreases with increasing hydrostatic pressure because the solubility of gases increases with the gas pressure in the cavity [17]. Prior experimental studies by Samaha et al. [6] and Lv et al. [18] show that longevity decreases exponentially with hydrostatic pressure in superhydrophobic fibrous coatings and micropore arrays, respectively. The use of Eq. (1) by Emami et al. [13,17] and others [11,18] in experimental studies implies a convective as opposed to purely diffusive mass transfer process in the liquid phase. Rahn and Paganelli show that the value of $ξ¯$ is equivalent to the ratio of the diffusion coefficient and an assigned effective boundary layer thickness [19]. Therefore, the accuracy of the empirical model depends on the careful choice of $ξ¯$⁠, which can be arbitrary [11].

For applications dominated by diffusion, longevity can be predicted by solving the transient diffusion equation, with no need for computing $ξ¯$⁠, as per the Epstein–Plesset model for gas diffusion in bubbles surrounded by liquid [20]. Sogaard et al. successfully correlated experimental results of diffusion-induced Cassie-to-Wenzel state transition on a microcavity array with a purely numerical solution of the 1D, transient gas-diffusion equation [21]. Results in their study show that one can predict, with reasonable accuracy, the longevity of a SHPo surface based on geometry, thermophysical properties, and initial and boundary conditions [21]. In this paper, we analytically solve the 1D transient diffusion equation for the dissolved gas concentration field in the liquid domain during a Cassie-to-Wenzel state transition to predict the longevity of a parallel ridge-type SHPo surface. We, and not Sogaard et al., also consider the cases where (a) gas diffusion into the liquid is driven by hydrostatic pressure rather than an external pressure due to a piston or high-pressure gas pressurizing the liquid phase and (b) the initial dissolved gas concentration in the liquid is arbitrary.

where p_g(t) is the pressure of the gas in the cavity, ρ_ℓ is the density of the liquid, g is the acceleration due to gravity, h(t) is the height of the liquid domain, σ is the surface tension along the meniscus, and R(t) is the radius of curvature of the meniscus.

We consider gas diffusion in quiescent liquid initially in the Cassie state above a ridge-type structured surface as depicted in Fig. 1. The liquid is supported above one period of the ridges (2d), the gas cavity height is s₀, and the width of the meniscus is 2a. The meniscus has a contact angle of θ(t) measured from the vertical and the corresponding radius of curvature is R(t). The thin dashed lines represent the shape of the meniscus at critical time t_cr, intermediate time t where t_cr < t < t_f, and final time t_f. The triple contact line is initially pinned at y = 0, but it is displaced by s(t) for t_cr < t < t_f. We analyze two representative studies where the boundary condition at y = h and the role of hydrostatic pressure are the main differences. In case study I, a massless, impermeable piston applies a constant pressure at y = h (for t > 0) which far exceeds the hydrostatic pressure of the liquid, i.e., p_∞ ≫ ρ_ℓgh as per Fig. 1(a); therefore, the gas pressure in the cavity is initially equal to the liquid pressure. In case study II, the top of the liquid is exposed to the atmosphere and hydrostatic pressure of the liquid is important as per Fig. 1(b); therefore, there is an initial pressure difference across the meniscus. It follows from the Young–Laplace equation that the meniscus is initially flat for case study I and has a finite radius of curvature R₀ for case study II as shown in the figure.

The gas concentration field in the gas phase is always assumed to be uniform because the molecular diffusivity in the gas phase is much larger than that in the liquid phase.

where the − sign corresponds to a downward protruding meniscus and vice versa. Only the downward protruding case is considered herein, so only the − sign is needed in subsequent equations.

where R(t) can be expressed in terms of p_g (t) as per Eq. (2).

where p_g is initially the same as the liquid pressure, i.e., p_g (0) = p_∞. It follows that the meniscus is initially flat, i.e., R(0) = ∞ and θ(0) = π/2.

where $ρg*(0,τ−1)=0$ to account for the (discontinuous) change in the partial density of the gas along the meniscus from ρ_g,0 to p_g (0)H at τ₀ = 0⁺. Overall, the time-varying dissolved gas concentration on the liquid side of the meniscus is a continuous function, ρ_g (0, t), as depicted by the thick solid line in Fig. 2. This function has been approximated by discrete values evaluated at times τ₁, τ₂,…, τ_j−1, τ_j where j denotes the time-step.

where p_g (0⁺) = p_∞ and the last term on the right-hand side, i.e., corresponding to i = j = 1 or t₁ − τ₁ = 0, is zero as per the initial condition.

where p_g (τ₁) = p_g (t₁).

where $mg,ℓ′=m0′−mg,c′+a h ρg,0$ is the mass of gas in the liquid domain and N denotes one plus the number of terms in the sum over n.

where p_g (t_j−1) is the gas pressure in the cavity for the interval t_j−1 < t < t_j and R(t_j−1) can be expressed in terms of p_g (t_j−1) as per the Young–Laplace equation. The variable p_g (t_j−1) is then the only unknown and it is computed numerically and then substituted into Eq. (34) to compute the dissolved gas concentration field, ρ_g,bc (y, t_j). If, at the end of time-step j, the total mass of gas left in the cavity is such that p_g (t_j) > p_cr, then ρ_g,bc (y, t_j+1) is computed in a similar manner. Conversely, if p_g (t_j) < p_cr, then the critical time is located in the interval t_j < t < t_j+1 in which case the dissolved gas concentration field is computed using the “after critical time” model.

Critical time depends on the geometry of the structures, height of the liquid domain, and initial dissolved gas concentration in the liquid. If the liquid solution is saturated with gas before critical time, then there is no diffusion-induced loss of gas in the cavity. It follows that the gas pressure in the cavity and meniscus deflection are constant as per the ideal gas law and Young–Laplace equations, respectively.

and solving for t_f.

The pressure of the compressed gas can be computed numerically by substituting Eq. (2) into Eq. (56) and substituting the resulting expression into Eq. (55). The temperature in the gas cavity may rise due to the moving boundary work done on the gas, therefore, we considered both adiabatic and isothermal compression.

with $mg,c′(tj−1)$ given by Eq. (76), but evaluated at t_j−1. The pressure of gas, p_g (t_j−1), is computed numerically, then ρ_g (y, t_j) is computed from Eqs. (71), (72), and (74).

The before and after critical time solutions described herein were implemented using MATLAB^® for both case studies. The gas in the cavity and above the liquid was assumed to be 100% nitrogen, M_g = 0.0280 kg/mol [27], and the liquid to be water. Thermophysical properties, evaluated at an ambient temperature of 25 °C and at pressure p_∞ = 1 atm, were set to the following values: $H=1.82×10−7 kg/(m3Pa)$ [27], $Dgl=2.01×10−9 m2/s$ [28], and σ = 0.073 N/m [29].

In order to demonstrate the functionality of the model presented in case study I, a sample set of dimensions and initial condition of $s0=10 μm, a=5 μm, h=800 μm$ and ρ_g,0 = 0 was considered. The step sizes before and after critical time were set to 0.001 s and 1 s, respectively. Figure 3 shows that the gas pressure in the cavity decreases exponentially from a maximum pressure of p_∞ to p_cr at t_cr = 1.82 s. The longevity is 176.7 s; therefore, t_cr is 1% of t_f and thus insignificant. Nevertheless, it is important to use the “before critical time” model for applications where the meniscus is always stable, for example, when ρ_g,0 = 0.018 kg/m³ as shown in Fig. 3.

The computed dissolved gas concentration fields at time t_cr and t_f are presented in Fig. 4 where two cases are compared: (a) a simplified case where the gas pressure is held constant at p_g (t) = p_∞ and (b) a case where the gas pressure is allowed to decrease from p_∞ to p_cr during the interval 0 < t < t_cr as depicted in the plot of p_g (t) in Fig. 3. It follows that, in the latter case, within the interval 0 < t < t_cr, the surface boundary condition decreases from p_∞H to p_crH as shown in the inset of Fig. 5. In the vicinity of the meniscus, the dissolved gas concentration gradient is initially steep and surface concentration is at its maximum (p_∞H). Gradually, the dissolved gas concentration gradient flattens as the gas diffuses into the bulk of the liquid domain and as the surface concentration decreases. It is for that reason that the plot of ρ_g (y, 0.001 s) crosses that of ρ_g (y, 1.82 s) in Fig. 5.

The mass of gas in the cavity, $mg,c′$⁠, is presented in Fig. 6. It decreases with time until complete wetting is achieved. Figure 7 presents the longevity for the case at hand as a function of h. The effect of h on t_f is shown to be more important when the height is close to the minimum height $hmin$⁠; longevity approaches infinity as h approaches $hmin$ from the right. Consequently, it is possible to achieve infinite longevity, i.e., stability, when $h<hmin$⁠.

Case study II accounts for hydrostatic pressure and gas diffusion at the exposed top surface. The dimensions s₀ = 85 μm, a = 73.5 μm, and h = 16.5 cm were chosen based on an experimental study by Xu et al. [11]. The temperature rise in the gas cavity as a result of adiabatic compression evaluates to only 1.3 K and the difference in initial gas pressure in the cavity is 0.01% as per Eq. (55) where γ = 1.4; therefore, the process was assumed to be isothermal. For samples of the same geometry, longevity varies with the initial dissolved gas concentration as shown in Fig. 8. As the initial gas concentration approaches 0.8 p_∞H, critical time and longevity approach 25 min and 48 h, respectively. For initial gas concentration greater than 0.8p_∞H, the rate of gas diffusion across the meniscus is so slow that the meniscus can be maintained in the Cassie state for indefinite periods. Xu et al. measured critical time and longevity, on a single trench submerged under a water column of the same dimensions and ρ_g,0 = p_∞H, to be about 2 h and 10 h, respectively [11]. The reason for the differences is that the model in case study II is applicable only to SHPo surfaces with low solid fraction, whereas the experimental setup in Xu et al. implies large solid fraction.

Figure 9 shows that t_cr decreases with the increasing cavity depth while t_f increases. The cavity depth is limited to a minimum of 15 μm to prevent the sagging meniscus from contacting the substrate before critical time. For s₀ = 15 μm, the triple contact line depins at t_cr = 110 s and transition to Wenzel state starts almost immediately, at t_f = 170 s. This highlights another extreme case where it is necessary to consider the time-dependent boundary condition for 0 < t < t_cr. However, t_cr ≪ t_f for a deep cavity, i.e., s₀ = 85 μm, where t_cr = 89 s and t_f = 2.7 h. In a deeper cavity, the initial mass of gas in the cavity is larger, hence it takes less time to reach critical time, but more time to deplete the gas in the cavity through gas diffusion.

The computed results before and after critical time were validated using the partial differential equation Toolbox in MATLAB^®. In this case, Eq. (21) subjected to the boundary conditions (24) and (25) was solved iteratively using a finite element method and the solution process was as follows. At the time instance t_n (t₀ = 0 s), the algorithm computes the pressure of the gas p_g (t_n) from Eq. (38) using the current mass of the gas in the cavity $mg,c′(tn)$⁠, where $mg,c′(t0)=m0′$⁠. If the computed p_g (t_n) < p_cr, p_g (t_n) is set equal to p_cr for all subsequent time-steps. Next, the algorithm updates the boundary condition (24) and Eq. (43) is solved over one time interval Δt, i.e., until t_n+1 = t_n + Δt, to compute the new dissolved gas concentration field in the liquid ρ_g (y, t_n+1). Consequently, the algorithm integrates ρ_g (y, t_n+1) over the domain to compute the total mass of the gas that has diffused into the liquid from t₀ until t_n+1 denoted as $mg,ℓ′(tn+1)$⁠. Thence, the net change of this mass of the diffused gas over the current time-step $dmg,ℓ′(tn)$ is computed by subtracting $mg,ℓ′(tn)$ from $mg,ℓ′(tn+1)$⁠. It is noted that $mg,ℓ′(t0)=0$⁠. Next, $mg,c′(tn+1)$ is computed by subtracting $dmg,ℓ′(tn)$ from $mg,c′(tn)$⁠. Then, $mg,c′(tn+1)$ is set equal to $mg,c′(tn)$ and the solution process is repeated until the maximum specified time instance. The domain was discretized with 1920 elements. The time-step was set equal to 5 × 10⁻⁴ s and 5 × 10⁻² s for t_n < 2 s and t_n ≥ 2 s, respectively, i.e., equal to the half of the corresponding time-steps that were used in the semi-analytical analysis. For case study I, the results are in excellent agreement, with the discrepancies for t_cr and t_f to be less than 0.09% and 0.14%, respectively, and the maximum discrepancy for $mg,ℓ′$ found to be less than 0.49%. For case study II, as represented in Fig. 9, the discrepancies for t_cr and t_f were less than 1.15% and 0.9%, respectively.

We present a semi-analytical, 1D, transient model for gas diffusion-induced Cassie to Wenzel state transition on ridge-type structured surfaces with low solid fractions. Our model accounts for meniscus curvature insofar as surface tension forces and gas cavity volume. For water on fluoropolymer-coated surfaces, it is sufficient to assume 1D gas diffusion because the meniscus is almost flat given that the advancing contact angle of water on a fluoropolymer-coated silicon surface is close to 110 deg. In addition, our model captures the effects of time-dependent surface concentration on the meniscus, hydrostatic pressure, and initial dissolved gas concentration on the longevity of SHPo surfaces. In a future publication, we will extend the analysis to 2D liquid domains in order to fully capture the effects of meniscus curvature and solid fraction.

Based on results from selected case studies, we demonstrate that, with appropriate dimensions and initial conditions, it is possible to stabilize the meniscus in a pinned or depinned Cassie state (metastable) for long periods if the liquid domain reaches saturation. We then conclude that precharged liquid domains combined with SHPo surfaces that have deep cavities tend to last longer in the Cassie state. We see potential applications of this 1D solution in the design of water-based microgap cooling systems where the primary concern is to optimize channel pressure drop, longevity, geometry, and materials for superior cooling performance.

The work of Toby Kirk was supported by an EPSRC doctoral scholarship. The computations in this paper were executed on the Tufts High-Performance Computing Research Cluster at Tufts University.

• Google (Google Faculty Research Award to Professor Marc Hodes)

• National Science Foundation (1402783)

a =

width of the gas cavity (m)

A =

surface area (m²)

B_n =

nth coefficient of the series for ρ_g,ac (kg/m³), case I

C_n =

nth coefficient of the series for ρ_g,ac (kg/m³), case II

d =

pitch of parallel ridges (m)

D_gl =

molecular diffusivity of gas in liquid (m²/s)

g =

acceleration due to gravity (m/s²)

h =

height of liquid domain above triple contact line (m)

H =

Henry's constant (kg m⁻³Pa⁻¹)

M =

molar mass (kg/mol)

$m′$ =

mass per unit depth (kg/m)

n =

positive integer

N =

positive integer

p =

pressure (Pa)

p_∞ =

atmospheric pressure (Pa)

R =

meniscus radius of curvature (m)

R_u =

universal gas constant (J mol⁻¹ K⁻¹)

s =

distance moved by the triple contact line (m)

t =

time (s)

T =

absolute temperature (K)

$V′$ =

volume per unit depth (m³/m)

Greek Symbols

β =

eigenvalue

γ =

exponent for adiabatic or isothermal compression

θ =

contact angle (rad)

λ =

eigenvalue

$ξ¯$ =

average invasion coefficient (m² s/kg)

ρ_ℓ =

density of liquid phase (kg/m³)

ρ_g =

partial density of gas in liquid phase (kg/m³)

ρ_sat =

equilibrium concentration of gas (kg/m³)

$ρg*$ =

instantaneous partial density of gas minus the initial partial density of gas in the liquid phase (kg/m³)

σ =

surface tension (N/m)

τ =

time variable (s)

Φ =

auxiliary solution