
==== Front
Langmuir
Langmuir
la
langd5
Langmuir
0743-7463
1520-5827
American Chemical Society

37708252
10.1021/acs.langmuir.3c01051
Article
Harvesting Energy from Changes in Relative Humidity Using Nanoscale Water Capillary Bridges
https://orcid.org/0000-0001-5958-8521
Tang Binze *†
Buldyrev Sergey V. *‡
Xu Limei *†§∥
https://orcid.org/0000-0003-1149-0693
Giovambattista Nicolas *⊥#
† International Center for Quantum Materials, School of Physics, Peking University, Beijing 100871, China
‡ Department of Physics, Yeshiva University, 500 West 185th Street, New York, New York 10033, United States
§ Collaborative Innovation Center of Quantum Matter, Beijing, 100190, China
∥ Interdisciplinary Institute of Light-Element Quantum Materials and Research Center for Light-Element Advanced Materials, Peking University, Beijing 100871, China
⊥ Department of Physics, Brooklyn College of the City University of New York, Brooklyn, New York 11210, United States
# Ph.D. Programs in Chemistry and Physics, The Graduate Center of the City University of New York, New York, New York 10016, United States
* E-mail: bztang@pku.edu.cn.
* E-mail: buldyrev@yu.edu.
* E-mail: limei.xu@pku.edu.cn.
* E-mail: ngiovambattista@brooklyn.cuny.edu.
14 09 2023
26 09 2023
14 09 2024
39 38 1344913458
19 04 2023
16 08 2023
© 2023 The Authors. Published by American Chemical Society
2023
The Authors
https://creativecommons.org/licenses/by-nc-nd/4.0/ Permits non-commercial access and re-use, provided that author attribution and integrity are maintained; but does not permit creation of adaptations or other derivative works (https://creativecommons.org/licenses/by-nc-nd/4.0/).

We show that nanoscale water capillary bridges (WCB) formed between patchy surfaces can extract energy from the environment when subjected to changes in relative humidity (RH). Our results are based on molecular dynamics simulations combined with a modified version of the Laplace–Kelvin equation, which is validated using the nanoscale WCB. The calculated energy density harvested by the nanoscale WCB is relevant, ≈1700 kJ/m3, and is comparable to the energy densities harvested using available water-responsive materials that expand and contract due to changes in RH.

Office of Advanced Cyberinfrastructure 10.13039/100000105 HRD-1547380 National Key Research and Development Program of China 10.13039/501100012166 2021YFA1400501 National Natural Science Foundation of China 10.13039/501100001809 12204039 National Natural Science Foundation of China 10.13039/501100001809 11935002 Division of Chemistry 10.13039/100000165 2223461 Division of Chemistry 10.13039/100000165 1856704 document-id-old-9la3c01051
document-id-new-14la3c01051
ccc-price
==== Body
pmcIntroduction

The capillary condensation of water is a common phenomenon in Nature, with important scientific and technological applications.1 For example, the formation of water capillary bridges (WCB) in porous rocks can affect the flow of supercritical carbon dioxide, which is important for the development of carbon capture/storage technologies.2−5 WCB can also induce strong capillary forces. This is particularly relevant in the production of micro- and nano-electro-mechanical systems (MEMS and NEMS); in such devices, capillary forces can lead to damages during the fabrication process.6 WCB can also form spontaneously in atomic force microscopy (AFM) experiments, between the AFM tip and the substrate of interest, affecting the corresponding measurements.7,8 Adhesive forces induced by WCB among colloidal9 and nanoscale particles10 can induce self-assembly process or trap the system in specific structures.

The phenomenon of water capillary condensation is usually described well by the Laplace–Kelvin (LK) equation. The LK equation relates the surface curvature of the WCB to the temperature and relative humidity (RH) of the environment. The LK equation is based on classical thermodynamics, where the liquid is modeled as a continuum and is built on the concept of surface tension. As such, the LK equation is expected to break down at the nanoscale, at which point the system dimensions become comparable to the size of the water molecules. Surprisingly, numerous studies indicate that the LK equation remains valid down to 10 nm,11−17 and recent studies indicate that it may hold at even smaller scales, <2 nm.18−22 In this work, we perform molecular dynamics (MD) simulations of WCB confined by two parallel, silica-based walls at different values of RH. The first aim of this work is to test whether the LK equation holds at length scales of nm. The LK equation is built on the assumption that the vapor obeys the ideal gas equation of state (EOS) and that the liquid is incompressible (with a density much larger than the density of the vapor). However, our simulations indicate that the ideal gas EOS does not hold for the case of water. Here, we show that a modified LK equation can be used to characterize the capillary condensation of water on the nanoscale.

The second aim of this work, and our main motivation, is to study the energy exchanges between a WCB and its environment (vapor) upon changes in the RH, a topic that has been overlooked in the past. This is motivated by a previous work where we show that, surprisingly, nanoscale WCB can be used to store energy at a given temperature, with non-negligible energy densities.23,24 If the energy exchanges between the vapor and the nanoscale WCB due to changes in RH are indeed relevant, then it would be possible, in principle, to create a water responsive (WR) material that would host a large number of nanoscale WCB per unit volume, which could harvest energy from variations in the RH of the surroundings. For example, such a material could be composed of stacked monolayers decorated with hydrophobic and hydrophilic nanoscale domains. WR materials that can harvest energy from changes in RH are already available.25−30 In most cases, these materials respond to RH changes due to complex molecular-level mechanisms that are not always clear, and their efficiency can vary considerably from one system to another. Interestingly, WR materials that can expand and contract upon variations of RH in the environment have been used to generate mechanical work, for the development of muscle-like actuators for miniature and biomimetic robotics, and to produce electricity.26,31

Computer Simulations Details

The molecular dynamics (MD) simulations and numerical procedures employed in this work are described in detail in refs (23 and 32). Briefly, we consider two identical hydrophobic silica-based walls of side length L = 14 nm that are decorated with a hydrophilic stripe-shaped patch; see Figure 1a. The hydrophilic patch (hydroxilated silica) extends across the wall and is 5 nm wide. The two walls are placed parallel to each other, with the two hydrophilic patches mirroring one another. MD simulations are performed with wall–wall separations h ranging from h = 3.0 to 7.0 nm, at 0.5 nm intervals. We consider a cubic system with side length L so that the walls expand across the system parallel to the x–y plane; periodic boundary condition apply along the x-, y-, and z-direction. At each fixed h, an orthorhombic volume of liquid water is placed between the walls, fully covering the patches and expanding from one hydrophilic patch to the other. Soon after the MD simulation starts, a water capillary bridge (WCB), periodic along the y-direction, forms between the walls. Snapshots of the system at h = 3 and 7 nm are included in Figure 1b,c; a schematic diagram of a WCB is included in Figure 1d. Depending on h and the RH, the WCB is composed of N water molecules with 5500 ≤ N ≤ 17000.

Figure 1 Water capillary bridge, confining surfaces, and thermodynamic path considered. (a) Top view of the walls employed in this work. The walls are hydrophobic, with a silica-based structure, and are hydroxilated over a 5 nm wide, stripe-shaped patch that extends along the wall; gray, red, and white spheres represent, respectively, the Si, O, and H atoms. The walls expand across the simulation box with periodic boundary conditions along the directions parallel to the walls. (b, c) Snapshots of a water capillary bridge expanding between two of the walls shown in (a) for wall separations: (b) h = 3 nm and (c) h = 7 nm. The WCB is in contact with the walls hydrophilic patches and is translationally symmetric along the direction of the patch (i.e., it has an effective infinite length along the y-axis). The boxes indicated in light blue are the volumes used to estimate the density of the vapor. (d) Schematic diagram of a capillary bridge showing all relevant variables used in this work. (e) Thermodynamic closed path considered to calculate the energy harvested by the nanoscale WCB. The relative humidity varies in the range RHmin ≤ RH ≤ RHmax; the wall separation varies within the range hmin ≤ h ≤ hmax. Path 1 → 2 (3 → 4) corresponds to a compression (expansion) of the WCB at constant low-RH (high-RH). Along path 2 → 3 (4 → 1) the RH increases (decreases) while the height of the WCB remains unchanged.

All MD simulations are performed using the LAMMPS software package.33 The equations of motion are integrated using the velocity Verlet algorithm with a time step of dt = 1 fs. Each simulation was performed at constant N, V, and T; a Nosé–Hoover thermostat is used to control the temperature of the system, T = 400 K. Simulations are performed for 5 ns. The first 2 ns are for equilibration, and the remaining 3 ns are used for data analysis. Configurations are saved every 1 ps for data analysis.

Water molecules are represented by using the SPC/E water model. The structure of the silica walls as well as their interactions with water are described in detail in refs (32 and 34). Briefly, the silica walls are composed of Si and O atoms arranged in a β-cristobalite structure. The O and Si atoms interact with the water O atoms via Lennard–Jones (LJ) interactions; there are no LJ interactions with the water H atoms. In the case of the hydrophobic regions of the walls, the O and Si atoms have no partial charges. The hydrophilic patch is obtained by hydroxilating the corresponding wall surface area. In this case, the surface silanol groups have no net charge, but the corresponding O, Si, and H atoms do have partial charges, and hence, they interact with water molecules via electrostatic interactions. The partial charges of the Si, O, and H atoms of the surface silanol groups are qO = −0.71 e, qSi = 0.31 e, and qH = 0.4 e. We note that the values of qO and qH are comparable to the partial charge of the O and H atoms in the SPC/E water model, qO,W = −0.8476 e and qH,W = 0.4238 e. As shown in Figure 1a, the silanol groups are distributed homogeneously over the hydrophilic patch of the walls with a surface density of x = 5 silanol groups per nm2. The (SPC/E) water contact angle with our fully hydrophobic wall is θ = 108° and θ ≈ 0° for the case of a fully hydroxilated wall.32,35

Results

The main goal of this work is to study the energy exchanged by a WCB and its environment (vapor) in the thermodynamic loop shown in Figure 1e. To do so, we first search for an expression that relates the curvature of the WCB and the RH of the surrounding vapor. If the vapor obeyed the ideal gas equation of state, such a relationship would be provided by the LK equation (see Supporting Information),1

where γ is the liquid–vapor surface tension, νl is the number density of liquid water at the working conditions, and r is the radius of the capillary bridge (r > 0 for convex capillary bridges, r < 0 for concave capillary bridges; see Figure 1d);P and P0 are the pressure and saturation pressure of the vapor. However, we find that the SPC/E vapor does not behave as an ideal gas. As shown in Figure 2, the equation of state (EOS) of the vapor at T = 400 K can be approximated by using different expressions, including (i) the van der Waals EOS, (where V′ and A are constants; Nv and V are the number of molecules and volume of the vapor), (ii) the rescaled ideal gas EOS, (where z is a constant), (iii) the truncated Virial EOS, (where B2 is a constant), and (iv) the shifted-volume EOS proposed in ref (36)2

Although all the EOS (i)–(iv) could be used to model the pressure of the vapor at the conditions studied, for simplicity, we will only consider eq 2. From Figure 2, we obtain B = −2.526 nm3. Since B < 0, one cannot interpret B as an excluded volume correction due to the effective volume occupied by the molecules in the gas. However, it is possible to interpret B as an effective, temperature-dependent volume correction that considers the nonideal nature of the gas caused by molecular attractions.

Figure 2 Pressure equation of state of SPC/E water in the vapor state at T = 400 K. Lines correspond to fits using different equations of state (see text).

In the Supporting Information, we obtain the modified Laplace–Kelvin (mLK) equation that follows when the vapor EOS is given by eq 2. When applied to the WCB shown in Figure 1, the mLK equation reads (see Supporting Information)3

where f is an auxiliary function defined as4

which depends on the chemical potential of the system, μ; f0 is the auxiliary function f evaluated at the saturation pressure P0 (at which the equilibrium vapor and liquid are separated by a flat interface) and b′ = B/N. Although eq 3 does not directly involve the RH, the RH can be easily obtained from eq 4,5

or, using eq 3,6

In the next section, we first validate the mLK equation, eq 3, for the nanoscale WCB used in this study. This equation is then used to estimate the energy harvested by the WCB that is subjected to the thermodynamic loop shown in Figure 1e.

Modified Laplace–Kelvin Equation

In order to test eq 3, we perform more than 63 independent MD simulations of WCB composed of different number of water molecules 5500 ≤ N ≤ 17000 and at various wall separations 3 nm ≤ h ≤ 7 nm (see Figure S1 in Supporting Information; the WCB is unstable and breaks for h ≥ 7.5 nm and N ≤ 17000). From the MD simulations, we evaluate independently both sides of eq 3. The auxiliary function f(μ) for the vapor is evaluated using eq 4 with the vapor pressure P given by eq 2. Since eq 2 depends on Nv and V, we consider two volumes within the system that exclude the walls and the WCB; see, e.g., the blue boxes in Figure 1b,c. The values of Nv and V are then calculated directly from the MD simulations by averaging the number of molecules found in these blue boxes and dividing them by the corresponding volume.

The right-hand side of eq 3 depends on γ, νl, and r. To calculate γ, we use the Kirkwood method,37 as explained in our previous work.23 At T = 400 K, we find that γ = 0.0383 N/m, which is consistent with values reported in the literature for SPC/E water (e.g., in ref (37), it is reported that γ = 0.0376 at T = 400 K). To obtain νl, we perform an independent MD simulation of a water slab (>3 nm thick) expanding across the system box and in coexistence with its vapor (T = 400 K). In equilibrium, the system reaches its saturation pressure (P0 = 0.13–0.14 MPa) and the number density of water within the slab defines νl; we obtain νl = 30.28 nm–3 (T = 400 K; the corresponding value reported in ref (37) is νl = 30.44 nm–3). To calculate the radius of curvature of the WCB, r, we note that capillarity theory (CT) predicts that the WCB profile is a circle of radius r centered at a point (xc, zc = 0); see Figure 1d. The procedure to obtain the WCB profile from the MD simulations is detailed in refs (23, 32, and 38), where it is also shown that the CT works remarkably well down to approximately h > 3 nm. Accordingly, we obtain r by fitting the average profile of the WCB by a circle of radius r. As an example, we include in Figure 3a various WCB profiles obtained from MD simulations (symbols) with the corresponding CT prediction (lines) for the case h = 5.5 nm. (In the WCB profiles shown in Figure 1(a), the point closest to the upper and lower walls can be misleading.23 Accordingly, when calculating the theoretical profile using CT, we omit these two points during the fitting procedure.) Note that (i) all WCB cover completely the hydrophilic patch and (ii) the contact angle of water varies with N (and h); see ref (23).

Figure 3 (a) Profiles of selected WCB with different number of water molecules N; the wall separation is h = 5.5 nm in all cases. Symbols represent the WCB profiles obtained from MD simulations. In all cases, the prediction from CT (lines) fits remarkably well the MD results. (b) Force exerted on the walls by the WCB shown in (a) as a function of the corresponding water contact angle θ. Black circles are from the MD simulations. The blue squares are the CT predictions using eqs 8 and 9, with b obtained from (a). The solid red line is the best fit to the MD simulation data using eq 8 and considering b as a fitting parameter (b = beff = 6.32 nm). (c) Auxiliary function f(μ) for the vapor surrounding the WCB (eq 4) as a function of 1/r, where r is the radius of curvature of the WCB profile. Data points are obtained from independent MD simulations at varied wall separations h and number of water molecules N. The values of f(μ) obtained from the MD simulations are noisy (see text). The red line is the best fit of the MD data; the blue line is the prediction of the modified Laplace–Kelvin (eq 3). (d) Relative humidity as a function of 1/r obtained from (c) using eq 5.

Figure 3c shows the auxiliary function f(μ) for the vapor as a function of 1/r (symbols) collected from all our MD simulations performed at different N and h (see Figure S1 in the Supporting Information). The red line in Figure 3c is the best linear fit to the data and has a slope of k′ = 0.305 nm. The blue line is the prediction for the mLK equation; its slope is nm. While k′ is somewhat larger than k, we note that the MD simulation data points in Figure 3c have large error bars, and hence, the value of k′ has a large uncertainty. For a better comparison between the mLK and MD simulations, we implement the following analysis. (i) To each data point in Figure 3c, located at (1/ri, ln(fi/f0)) and with standard deviation δi (i = 1, 2, ..., n, where n is the total number of points), we associate a Gaussian distribution Gi(xi, σi). This Gaussian distribution is centered at xi = ln(fi/f0) and has a standard deviation of σi = δi. (ii) For each point i, we choose a random value yi that obeys the corresponding Gaussian probability distribution Gi(xi, δi). The set of values (y1, y2, ..., yn) so produced represent a potential distribution of data points that could be obtained from our MD simulations. (iii) This set of data points is also fitted by a straight line and the corresponding slope is calculated. (iv) By repeating steps (ii) and (iii) for a total of m = 10000 times, we obtain m different values of k′. By averaging all the values of k′ so obtained, we calculate the average slope k̅ = 0.355 nm and the corresponding standard deviation σ = 0.083 nm. While k̅ is still larger than the theoretical slope k, the deviations from MD simulations and the mLK prediction are within 1.5σ. Hence, we conclude that while the agreement between the MD simulations and the mLK equation is not perfect, the results from the MD simulations are consistent with the mLK. Accordingly, by using eq 6, one can estimate the RH of the vapor from the curvature of the WCB; see Figure 3d. In the next section, we will use this relationship to estimate the ability of the WCB to harvest energy from changes in RH.

Energy Harvested from Changes in Relative Humidity

The energy harvested by the WCB along the thermodynamic closed path shown in Figure 1e, Ehvd, is equal to the corresponding work done by the WCB on the walls along cycle 1 → 2 → 3 → 4. However, work on the walls is done only along the paths 1 → 2 and 3 → 4. Briefly, Ehvd = Wby WCBcycle = W1→2 + W3→4, where7

Here, F(h) is the force that the WCB exerts on the walls (F > 0 if the walls are being pushed away from each other). On path 3 → 4, the work is negative because the force acting on the wall is always pointing toward the WCB and the walls displacement is opposite to the force direction. (In our study, the plate separation is always h ≥ 3 nm. At these separations, the force on the plates is a negative and monotonic decaying function of h.23,39 However, for approximately h < 2 nm, AFM experiments40 and MD simulations41−43 have observed an oscillatory behavior in the water mediated interactions between confining surfaces.) Instead, the work produced on the walls is positive along path 1 → 2. Note that the curvature of the bridge is larger along the 1 → 2 path than along the 3 → 4 path. As we will show below, the force induced by the WCB with a curved interface (low RH, path 1 → 2) is larger than the force induced by the WCB with a flat interface (high RH, path 3 → 4). Therefore, the energy harvested by the WCB is positive.

As explained in refs (23 and 32), F(h) can be calculated directly from the MD simulations. In addition, MD simulations23 show that, for approximately h ≥ 3 nm, F(h) is in excellent agreement with the force predicted by CT,8

where w = L = 14 nm is the WCB length; θ = θ(h) is the contact angle of the WCB in contact with the patchy walls and is given by19

In eq 8, b is the thickness of the WCB at the walls (b = 2 x(±h/2)). As an example, we include in Figure 3b the values of F(h) obtained from MD simulations with h = 5.5 nm and for the WCB composed of different numbers of water molecules (circles). The blue squares are the prediction from CT using eqs 8 and 9, with b and r obtained from the corresponding WCB profile (solid lines in Figure 3a); depending on N, the WCB base thickness vary slightly within the range 6.0–7.4 nm.23 These values of b are slightly larger than the hydrophilic patch width (5 nm) and it is due to the ability of the water molecules of the WCB to wet the hydrophobic boundaries of the hydrophilic patch.23 Also included in Figure 3b is the best fit to the forces obtained from the MD simulations, using eq 8 and considering b as an effective, single fitting parameter, independent of N (red line). We find that for the effective WCB base thickness beff = 6.32 nm, the red line fits slightly better the force values calculated from the MD simulations.

To calculate Ehvd, all we need is F(h) at RH = RHmin and RHmax for all wall separations, hmin ≤ h ≤ hmax. Alternatively, since the mLK equation allows one to express the RH as a function of the WCB curvature, r = r(RH), to calculate Ehvd, we only need F(h) at constant WCB radii rmin(RHmin) and rmax(RHmax) and for hmin ≤ h ≤ hmax. In our MD simulations, we can control only N and h, but not r. Therefore, we proceed as follows. (i′) Perform p = 63 MD simulations at different values of h and N (see Figure S1 in the Supporting Information). (ii′) For each state point (N, h), the force F(h, N) is calculated from the MD simulations and the corresponding r(h, N) (as shown in Figure 3a). The variables are substituted to get F(h, r). (iii′) From steps (i′) and (ii′), one obtains p triplets of values [r, h, F(h, r)]. These sets of points are then fitted using eqs 8 and 9 with the effective WCB base thickness beff = 6.32 nm. This procedure results in a smooth function for the force, Fr(h), that is then used to calculate W1→2 and W3→4 via eq 7.

Figure 4a shows the forces obtained from the MD simulations for all values of N and h considered (symbols). The solid lines are the fitted forces Fr(h), for the constant r. It follows from Figure 4a that the force induced by the WCB on the walls is larger for small RH than for large RH, implying that the energy extracted by the WCB in the thermodynamic cycle of Figure 1e is indeed positive (see discussion above).

Figure 4 (a) Magnitude of the force exerted by the WCB on the walls as a function of the wall separation h. Circles are results form MD simulations; lines correspond to the forces predicted by CT at constant r, obtained from eqs 8 and 9; the value of r is indicated by the color bar on the right-hand side of the panel. Using the mLK, the values of r are also converted to values of RH. (b) Magnitude of the force exerted on the walls by the WCB as a function of h and RH. See also Figure S2 of the Supporting Information. Symbols indicate the magnitude of the force obtained from the MD simulations (see color bar); The contour map represents the forces calculated from eqs 8 and 9. The red closed path corresponds to the thermodynamic path for which the harvested energy density ρE is maximum (see also Figure 1e). The values of 1/r (upper axis) are obtained from the RH values (lower axis) using eq 6.

In order to calculate Ehvd, we consider the values of hmin, hmax, RHmin, and Rmax that maximize Wcycle for which the WCB is stable. The largest energy harvested by the WCB is Ehvd = 1.863 × 10–18 J, which corresponds to hmin = 3.0 nm, hmax = 6.5 nm, RHmin = 0.921, and RHmax = 0.995. For comparison with available water-responsive materials, we also calculate the corresponding energy density ρE = Ehvd/V0, where V0 is the volume of the system. Considering that, conservatively, V0 = L2 × hmax, one obtains ρE = 1462 kJ/m3. We can also optimize the closed path of Figure 1 to maximize ρE, instead of Ehvt. We find that the best value for the energy density harvested by the WCB is ρE = 1723 kJ/m3 and corresponds to hmin = 3 nm, hmax = 5 nm, RHmin = 0.865, and Rmax = 0.995 (For these values of (hmin, hmax, Rmin, Rmax), the number of water molecules at the states 1, 2, 3, and 4 indicated in the thermodynamic loop of Figure 4b are, respectively, 8028, 6500, 7500, and 12000. Accordingly, 1528 water molecules evaporate during the compression process 1 → 2, and 4500 water molecules condense during the decompression process 3 → 4.) This energy density is not negligible. Indeed, it is comparable to the energy densities found in water-responsive materials based on carbon nanotubes, ≈1800 kJ/m3 (subjected to a much larger change in RH, 0.1 ≤ RH ≤ 0.9), which have potential applications in the development of artificial muscles and actuators.44 See Table 1 for a comparison of the energy densities of the WCB and other water-responsive materials that can be used to harvest energy from changes in RH.

Table 1 Energy Densities ρE of Water-Responsive Materials That Can Be Used to Harvest Energy from Changes in RH

system	ρE (kJ/m3)	
spider silk29	≈500	
bacterial spore25	≈21300	
carbon nanotubes hybrid yarn44	≈1800	
regenerated silk film47	≈1626	
bioinspired polymer (PEE–PPy)48	≈114	
graphene fiber49	≈520	

For comparison, we also indicate in Figure 4a the RH value corresponding to the WCB radius r (using the mLK equation). The force acting on the walls as a function of the RH and walls separation is shown Figure 4b. Also included is the thermodynamic closed path for which the energy density extracted by the WCB is maximum (red lines).

Role of Temperature and Wall Patch Thickness on Ehvd and ρE

The findings discussed so far indicate that the results from MD simulations based on nanoscale WCB (h ≥ 3 nm) are consistent with (i) the (macroscopic) CT and (ii) the mLK equation. Accordingly, without performing further MD simulations, one can use eqs 2, 6, 8, and 9 to estimate the energy and energy density that a WCB can harvest from changes in RH (a) for different patch thicknesses b and (b) at different temperatures.

a b-Effects

Figure 5a shows the energy harvested by the WCB for selected values of hydrophilic patch width b at T = 400 K. Ehvd is calculated for the thermodynamic closed path of Figure 1e as a function of the minimum RH of the cycle. In these calculations, we set hmin = 3.0 nm since this is approximately the smallest wall separations at which CT holds.23,45 We also set RHmax = 1 since this value of RH minimizes the force produced by the WCB on the walls along path 3 → 4, which in turn maximizes the Ehvd. The value of hmax depends on RHmin, and it is set to and for b = 2r(RHmin), where r(RHmin) is the negative radius of curvature found from eq 6. For b > −2r(RHmin), this condition assures that at the maximum wall separations hmax = −2r(RHmin), the WCB does not detach from the walls. For b ≤ −2r(RHmin), this condition assures that, at the maximum walls separation , the two liquid–vapor interfaces of the WCB do not touch one another. Using eqs 2, 6, 7, 8, and 9, one can show that the analytical expression for the energy harvested by the WCB in the thermodynamic closed path of Figure 1e is given by10

where and .

Figure 5 Estimated (a) energy Ehvd(RHmin) and (b) energy density ρE(RHmin) harvested by the WCB for the thermodynamic closed path showed in Figure 1e at T = 400 K (hmin = 3 nm, RHmax = 1.0; hmax depends on RHmin, see text). The color bar indicates the width b of the wall hydrophilic patch; see Figure 1. Lines correspond to (bottom to top) b = 6, 10, 20, 40, 80, 150, 250, 350, 450, and 500 nm. The maximum harvested energy occurs for minimum relative humidity RHmin ≈ 0.93–0.99, depending on b. And the maximum harvested energy density occurs for minimum relative humidity RHmin ≈ 0.87–0.90, depending on b. (c) Maximum harvested energy Ehvd(b) for a wall hydrophilic patch of with b [defined by the maxima in (a)]. (d) Minimum RH and (e) maximum wall separations h for which Ehvd(b) reaches its maximum value in (a). (f) Maximum harvested density energy ρE*(b) for a hydrophilic wall patch of with b resulting from (b). (g) Minimum RH and (h) corresponding maximum wall separations for which ρE(b) reaches its maximum value in (b).

The maxima in Figure 5a correspond to the maximum energy that can be harvested by the WCB at a given b, Ehvd*(b); see Figure 5c. The corresponding minimum RH at which Ehvd(b) reaches its maximum, RHmin*, and the resulting maximum plate separation, hmax*, are shown in Figure 5d,e as a function of b. Interesting, RHmin* ≈ 0.93–0.99 for all values of b, which is rather large; hmax* < 83 nm, even for large b ≈ 400 nm.

As shown in Figure 5c, Ehvd*(b) ∝ b. To understand the linear behavior of Ehvd*(b), we note that the force produced by the WCB can be expressed as F(h) = FP(h) + Fγ(h), where (see eq 8)11

is the force due to the pressure of the liquid within the WCB, and12

is the force produced by the liquid–vapor interface of the WCB. Along paths 1 → 2 and 3 → 4 in Figure 1e, RH is constant and, hence, the mLK equation implies that r is constant as well. Therefore, for a given h and RH, θ is also constant (eq 9). Accordingly, for a given h and RH, FP(h) ∝ b, while Fγ(h) remains practically constant. This implies that the work of FP in the closed path of Figure 1e increases linearly with b, while the corresponding work done by Fγ does not change with b; see Figure 6. It follows that for large b, Ehvd* = Wby WCBcycle ∝ b; i.e., Ehvd*(b) increases linearly with increasing b because the dominant force produced by the WCB is FP. We note that the work produced by Fγ, while small, is negative while the work produced by FP is positive. This means that the energy harvested by the WCB is solely due to the work performed by FP.

Figure 5b shows the energy density harvested by the WCB for selected values of the hydrophilic patch width b at T = 400 K. Specifically, ρE(b) = Ehvd(b)/V(b), where V(b) = Lhmax(b + 2.0) × 1.2. This arbitrary volume is sufficiently large to contain the WCB at all values of h and RH studied. The maxima in Figure 5f correspond to the maximum energy density that can be harvested by the WCB at a given b, ρE*(b); see Figure 5. The corresponding minimum RH at which ρE(b) reaches its maximum, RHmin*(b), and the resulting maximum plate separation, hmax*(b), are included in Figure 5f–h. The values RHmin* ≈ 0.87–0.90 are still large, while hmax* = 5.5–6.0 nm is rather small. It follows from Figure 5f–h that all quantities, ρE*(b), RHmin*(b), and hmax*(b), saturate for large b. The fact that ρE*(b) saturates for large values of b may not be surprising since V(b) ∝ b for large b. Interestingly, the harvested energy density is no larger than ρE*(b) ≈ 5200 kJ/m3. This value is ≈3 times the energy harvested by our nanoscale WCB with b = 5 nm.

b T-Effects

Varying the temperature of the system is expected to alter key parameters in eqs 2, 6, and 8, which are used to estimate Ehvd(b) and ρE(b). Specifically, the fitting parameter B in eq 2, the saturation pressure P0, νl, and γ in eq 6, and the effective wall patch beff in eq 8 are, in principle, all T-dependent. Safely, one can assume that effective b is not affected by T. To calculate B, P0 and γ at different temperatures, we performed additional MD simulations; the corresponding values for these quantities are reported in Table 2.

Figure 6 Work harvested by the WCB forces as function of the wall hydrophilic patch width b. FP and Fγ are the forces due to the WCB inner pressure and liquid–vapor surface tension, respectively. At a given b, the work of FP and Fγ are evaluated for the thermodynamic closed path shown in Figure 1e and for which the harvested energy is maximum (see Figure 5c). The work produced by Fγ is negative and roughly constant for large b, meaning that Fγ does not contribute to the energy harvested by the WCB during the thermodynamic cycle considered. The WCB can harvest energy due to FP since its net work is positive; the work produced by FP increases with increasing the wall patch width, FP ∝ b. Inset: Magnification of the main panel for small values of b.

Table 2 Temperature Dependence of Various Parameters Used in This Worka

T (K)	300	350	400	450	500	
ρbulkliq (g/cm3)	0.987	0.952	0.906	0.847	0.769	
γ (N/m)	0.0551	0.0472	0.0383	0.0295	0.0205	
P0 (Pa)	2500.7	61004.4	133120.5	552769.8	1899326.2	
νv,0 (nm–3)	0.0006	0.0137	0.0257	0.1017	0.3581	
B (nm–3)	20.032	–6.133	–2.526	–1.411	–8.422	
a γ and ρbulkliq are the surface tension and density of bulk water estimated from MD simulations of a 3 nm thick film of liquid water coexisting with its vapor. P0 and νv,0 are the corresponding saturation pressure and vapor number density. B is the fitting parameter in the shifted-volume EOS of SPC/E water (eq 2).

Figure 7 Effects of temperature on (a) Ehvd*(b) and (b) ρE*(b) for the case of a WCB in contact with a hydrophilic patch of width b = 5 nm. Decreasing the temperature enhances the energy and energy density harvested by the WCB.

Using the parameters given in Table 2, we obtain the Ehvd*(b) and ρE*(b) at T = 300–500 K; see Figure 7. It follows that one can improve the harvested energy and energy density by working just at lower temperatures. As shown in Table 2, the considerable T-dependence of Ehvd*(b) and ρE*(b) is due to the large increase in the surface tension of water upon cooling; γ increases by ≈2.7 times as the temperature decreases from T = 500 K to T = 300 K.

Summary and Discussion

In this work, we studied the energy exchanged between a WCB and its surrounding vapor along the closed thermodynamic path, shown in Figure 1e. By following the close path in Figure 1e, the WCB is subjected to changes in RH as well as expansions and contractions. Similar closed thermodynamic paths have been used to study the efficiency of WR materials to convert changes in RH into work.26,28 In our simulations, the WCB expands between two hydrophilic patches located on opposite parallel surfaces (Figure 1a). The hydrophilic patches pin the WCB so it cannot diffuse along the directions parallel to the surfaces; the water contact angle, however, is not constrained, and it varies as the WCB expands and contracts (Figure 3a).

The main finding of this work is that a nanoscale WCB, subjected to the closed thermodynamic path of Figure 1e, can harvest a relevant amount of energy density (energy per unit volume), ρE ≈ 1700 kJ/m3. At T = 400 K, this energy is harvested while the wall separation varies in the range 3–5 nm, and the RH changes in the range ≈0.86–0.99%. For comparison, the ρE harvested by water-responsive materials based on carbon nanotubes is ≈1800 kJ/m3 (subjected to a much larger change in RH, 0.1 ≤ RH ≤ 0.9). Our results indicate that, at least theoretically, it may be possible to construct WR materials that host a large number of WCB in their interior. For example, such a material could be composed of stacked monolayers, each decorated with alternating hydrophobic/hydrophilic patches of nanoscale dimensions. The hydrophilic patches could pin a large number of WCB which could harvest energy as the RH of the environment varies. Interestingly, our MD simulations also show that the main source of work produced by the WCB along the closed loop of Figure 1e is not due to the liquid–vapor surface interface, but by the internal pressure in the WCB (Figure 6). In other words, the large value of ρE harvested by the nanoscale WCB is not due to a nanoscale effect, i.e., due to the large surface-to-volume ratio of the WCB, but because the WCB is just very small. A detailed analysis of the effects of WCB size on the energy harvested was also included.

From a practical point of view, our MD simulations show that macroscopic thermodynamics (capillarity theory) holds for nanoscale WCB formed between patchy walls separated by only h ≥ 3 nm.23,24 In particular, we find that a simple modified version of the macroscopic Laplace–Kelvin equation is consistent with the MD simulations. This implies, for example, that one could use macroscopic thermodynamics to explore capillarity phenomena at (1) nm-scales, and at different RH. For example, our theoretical calculation indicate that the energy harvested by a nanoscale WCB is expected to increase considerably with decreasing temperatures. The possibility of using thermal machines based on interfaces that convert heat into work have been noticed in the past.46

Supporting Information Available

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.langmuir.3c01051.Derivation of the modified Laplace–Kelvin equation; Two complementary figures showing (i) the state points simulated and (ii) the force produced by the WCB on the walls at different RH and for selected values of h (PDF)

Supplementary Material

la3c01051_si_001.pdf

The authors declare no competing financial interest.

Acknowledgments

L.X. and B.T. acknowledge the financial support by the National Natural Science Foundation of China (NSFC Grant No. 11935002 and No. 12204039) and the National Key R&D Program (Grant No. 2021YFA1400501). S.V.B. is grateful to NSF for financial support (Grant No. 1856704). N.G. is thankful for support from the NSF (Grant No. CHE-2223461) and to the NSF-CREST “Center for Interface Design and Engineered Assembly of Low Dimensional systems (IDEALS)” (Grant No. HRD-1547380). We are also grateful for computational resources provided by the supercomputer TianHe-1A in Tianjin and the High Performance Computing Platform of Peking University, China, and by a grant of computer time from the City University of New York High Performance Computing Center under NSF Grants CNS-0855217, CNS-0958379, and ACI-1126113.
==== Refs
References

de Gennes P. ; Brochard-Wyart F. ; Quéré D. Capillarity and Wetting Phenomena: Drops, Bubbles, Pearls, Waves; Springer-Verlag: New York, 2004.
Li B. ; Benson S. M. Influence of small-scale heterogeneity on upward CO2 plume migration in storage aquifers. Adv. Water Resour. 2015, 83 , 389–404. 10.1016/j.advwatres.2015.07.010.
Zhao B. ; MacMinn C. W. ; Juanes R. Wettability control on multiphase flow in patterned microfluidics. Proc. Natl. Acad. Sci. U. S. A. 2016, 113 , 10251–10256. 10.1073/pnas.1603387113.27559089
Arif M. ; Abu-Khamsin S. A. ; Iglauer S. Wettability of rock/CO2/brine and rock/oil/CO2-enriched-brine systems: Critical parametric analysis and future outlook. Adv. Colloid Interface Sci. 2019, 268 , 91–113. 10.1016/j.cis.2019.03.009.30999164
Sun E. W.-H. ; Bourg I. C. Molecular dynamics simulations of mineral surface wettability by water versus CO2: Thin films, contact angles, and capillary pressure in a silica nanopore. J. Phys. Chem. C 2020, 124 , 25382–25395. 10.1021/acs.jpcc.0c07948.
Maboudian R. ; Ashurst W. R. ; Carraro C. Tribological challenges in micromechanical systems. Tribol. Lett. 2002, 12 , 95–100. 10.1023/A:1014044207344.
Yang L. ; Tu Y. S. ; Tan H. L. Influence of atomic force microscope (AFM) probe shape on adhesion force measured in humidity environment. Appl. Math. Mech. (Engl. Ed.) 2014, 35 , 567–574. 10.1007/s10483-014-1813-7.
Weeks B. L. ; Vaughn M. W. ; DeYoreo J. J. Direct imaging of meniscus formation in atomic force microscopy using environmental scanning electron microscopy. Langmuir 2005, 21 , 8096–8098. 10.1021/la0512087.16114907
Kralchevsky P. A. ; Denkov N. D. Capillary forces and structuring in layers of colloid particles. Curr. Opin. Colloid Interface Sci. 2001, 6 , 383–401. 10.1016/S1359-0294(01)00105-4.
Laube J. ; Dormann M. ; Schmid H.-J. ; Madler L. ; Colombi Ciacchi L. Dependencies of the adhesion forces between TiO2 nanoparticles on size and ambient humidity. J. Phys. Chem. C 2017, 121 , 15294–15303. 10.1021/acs.jpcc.7b05655.
Fisher L. ; Israelachvili J. Direct experimental verification of the Kelvin equation for capillary condensation. Nature 1979, 277 , 548–549. 10.1038/277548a0.
Fisher L. R. ; Israelachvili J. N. Experimental studies on the applicability of the Kelvin equation to highly curved concave menisci. J. Colloid Interface Sci. 1981, 80 , 528–541. 10.1016/0021-9797(81)90212-5.
Fisher L. ; Gamble R. ; Middlehurst J. The Kelvin equation and the capillary condensation of water. Nature 1981, 290 , 575–576. 10.1038/290575a0.
Kohonen M. M. ; Christenson H. K. Capillary condensation of water between rinsed mica surfaces. Langmuir 2000, 16 , 7285–7288. 10.1021/la991404b.
Mitropoulos A. C. The kelvin equation. J. Colloid Interface Sci. 2008, 317 , 643–648. 10.1016/j.jcis.2007.10.001.17950748
Zhong J. ; Riordon J. ; Zandavi S. H. ; Xu Y. ; Persad A. H. ; Mostowfi F. ; Sinton D. Capillary condensation in 8 nm deep channels. J. Phys. Chem. Lett. 2018, 9 , 497–503. 10.1021/acs.jpclett.7b03003.29323911
Yang G. ; Chai D. ; Fan Z. ; Li X. Capillary condensation of single- and multicomponent fluids in nanopores. Ind. Eng. Chem. Res. 2019, 58 , 19302–19315. 10.1021/acs.iecr.9b04116.
Yang Q. ; Bi R. ; Banerjee D. ; Nasrabadi H. Direct Observation of the Vapor–Liquid Phase Transition and Hysteresis in 2 nm Nanochannels. Langmuir 2022, 38 , 9790–9798. 10.1021/acs.langmuir.2c00590.35913840
Yang Q. ; Sun P. ; Fumagalli L. ; Stebunov Y. ; Haigh S. ; Zhou Z. ; Grigorieva I. ; Wang F. ; Geim A. Capillary condensation under atomic-scale confinement. Nature 2020, 588 , 250–253. 10.1038/s41586-020-2978-1.33299189
Kim S. ; Kim D. ; Kim J. ; An S. ; Jhe W. Direct evidence for curvature-dependent surface tension in capillary condensation: Kelvin equation at molecular scale. Phys. Rev. X 2018, 8 , 041046 10.1103/PhysRevX.8.041046.
Zandavi S. ; Ward C. Contact angles and surface properties of nanoporous materials. J. Colloid Interface Sci. 2013, 407 , 255–264. 10.1016/j.jcis.2013.06.062.23890594
Factorovich M. H. ; Molinero V. ; Scherlis D. A. Vapor pressure of water nanodroplets. J. Am. Chem. Soc. 2014, 136 , 4508–4514. 10.1021/ja405408n.24580542
Tang B. ; Buldyrev S. V. ; Xu L. ; Giovambattista N. Energy Stored in Nanoscale Water Capillary Bridges between Patchy Surfaces. Langmuir 2020, 36 , 7246–7251. 10.1021/acs.langmuir.0c00549.32460499
Tang B.-Z. ; Yu X.-J. ; Buldyrev S. V. ; Giovambattista N. ; Xu L.-M. Energy stored in nanoscale water capillary bridges formed between chemically heterogeneous surfaces with circular patches. Chin. Phys. B 2020, 29 , 114703 10.1088/1674-1056/abb664.
Chen X. ; Mahadevan L. ; Driks A. ; Sahin O. Bacillus spores as building blocks for stimuli-responsive materials and nanogenerators. Nat. Nanotechnol. 2014, 9 , 137–141. 10.1038/nnano.2013.290.24463362
Chen X. ; Goodnight D. ; Gao Z. ; Cavusoglu A. H. ; Sabharwal N. ; DeLay M. ; Driks A. ; Sahin O. Scaling up nanoscale water-driven energy conversion into evaporation-driven engines and generators. Nat. Commun. 2015, 6 , 7346 10.1038/ncomms8346.26079632
Cavusoglu A.-H. ; Chen X. ; Gentine P. ; Sahin O. Potential for natural evaporation as a reliable renewable energy resource. Nat. Commun. 2017, 8 , 617 10.1038/s41467-017-00581-w.28951541
Park Y. ; Chen X. Water-responsive materials for sustainable energy applications. J. Mater. Chem. A 2020, 8 , 15227–15244. 10.1039/D0TA02896G.
Agnarsson I. ; Dhinojwala A. ; Sahni V. ; Blackledge T. A. Spider silk as a novel high performance biomimetic muscle driven by humidity. J. Exp. Biol. 2009, 212 , 1990–1994. 10.1242/jeb.028282.19525423
Boudot M. ; Elettro H. ; Grosso D. Converting water adsorption and capillary condensation in usable forces with simple porous inorganic thin films. ACS Nano 2016, 10 , 10031–10040. 10.1021/acsnano.6b04648.27792305
Wang H. ; Liu Z.-L. ; Lao J. ; Zhang S. ; Abzalimov R. ; Wang T. ; Chen X. High Energy and Power Density Peptidoglycan Muscles through Super-Viscous Nanoconfined Water. Adv. Sci. 2022, 9 , 2104697 10.1002/advs.202104697.
Giovambattista N. ; Almeida A. B. ; Alencar A. M. ; Buldyrev S. V. Validation of capillarity theory at the nanometer scale by atomistic computer simulations of water droplets and bridges in contact with hydrophobic and hydrophilic surfaces. J. Phys. Chem. C 2016, 120 , 1597–1608. 10.1021/acs.jpcc.5b10377.
Plimpton S. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 1995, 117 , 1–19. 10.1006/jcph.1995.1039.
Giovambattista N. ; Rossky P. J. ; Debenedetti P. G. Effect of pressure on the phase behavior and structure of water confined between nanoscale hydrophobic and hydrophilic plates. Phys. Rev. E 2006, 73 , 041604 10.1103/PhysRevE.73.041604.
Giovambattista N. ; Debenedetti P. G. ; Rossky P. J. Effect of surface polarity on water contact angle and interfacial hydration structure. J. Phys. Chem. B 2007, 111 , 9581–9587. 10.1021/jp071957s.17658789
Elliott J. R. ; Lira C. T. ; Lira C. T. Introductory chemical engineering thermodynamics; Prentice Hall: Upper Saddle River, NJ, 2012; Vol. 668 .
Ismail A. E. ; Grest G. S. ; Stevens M. J. Capillary waves at the liquid-vapor interface and the surface tension of water. J. Chem. Phys. 2006, 125 , 014702 10.1063/1.2209240.16863319
Giovambattista N. ; Debenedetti P. G. ; Rossky P. J. Hydration behavior under confinement by nanoscale surfaces with patterned hydrophobicity and hydrophilicity. J. Phys. Chem. C 2007, 111 , 1323–1332. 10.1021/jp065419b.
Almeida A. B. ; Giovambattista N. ; Buldyrev S. V. ; Alencar A. M. Validation of capillarity theory at the nanometer scale. II: Stability and rupture of water capillary bridges in contact with hydrophobic and hydrophilic surfaces. J. Phys. Chem. C 2018, 122 , 1556–1569. 10.1021/acs.jpcc.7b08577.
Uhlig M. R. ; Benaglia S. ; Thakkar R. ; Comer J. ; Garcia R. Atomically resolved interfacial water structures on crystalline hydrophilic and hydrophobic surfaces. Nanoscale 2021, 13 , 5275–5283. 10.1039/D1NR00351H.33624666
Leoni F. ; Calero C. ; Franzese G. Nanoconfined fluids: Uniqueness of water compared to other liquids. ACS Nano 2021, 15 , 19864–19876. 10.1021/acsnano.1c07381.34807577
Calero C. ; Franzese G. Water under extreme confinement in graphene: Oscillatory dynamics, structure, and hydration pressure explained as a function of the confinement width. J. Mol. Liq. 2020, 317 , 114027 10.1016/j.molliq.2020.114027.
Engstler J. ; Giovambattista N. Different temperature-and pressure-effects on the water-mediated interactions between hydrophobic, hydrophilic, and hydrophobic–hydrophilic nanoscale surfaces. J. Chem. Phys. 2022, 157 , 064701 10.1063/5.0097908.35963716
Kim S. H. ; Kwon C. H. ; Park K. ; Mun T. J. ; Lepró X. ; Baughman R. H. ; Spinks G. M. ; Kim S. J. Bio-inspired, moisture-powered hybrid carbon nanotube yarn muscles. Sci. Rep. 2016, 6 , 1–7. 10.1038/srep23016.28442746
Almeida A. B. ; Buldyrev S. V. ; Alencar A. M. ; Giovambattista N. How Small Is Too Small for the Capillarity Theory?. J. Phys. Chem. C 2021, 125 , 5335–5348. 10.1021/acs.jpcc.0c11140.
Laouir A. ; Luo L. ; Tondeur D. ; Cachot T. ; Le Goff P. Thermal machines based on surface energy of wetting: Thermodynamic analysis. AIChE J. 2003, 49 , 764–781. 10.1002/aic.690490320.
Park Y. ; Jung Y. ; Li T.-D. ; Lao J. ; Tu R. S. ; Chen X. β-Sheet Nanocrystals Dictate Water Responsiveness of Bombyx Mori Silk. Macromol. Rapid Commun. 2020, 41 , 1900612 10.1002/marc.201900612.
Ma M. ; Guo L. ; Anderson D. G. ; Langer R. Bio-inspired polymer composite actuator and generator driven by water gradients. Science 2013, 339 , 186–189. 10.1126/science.1230262.23307738
Cheng H. ; Hu Y. ; Zhao F. ; Dong Z. ; Wang Y. ; Chen N. ; Zhang Z. ; Qu L. Moisture-activated torsional graphene-fiber motor. Adv. Mater. 2014, 26 , 2909–2913. 10.1002/adma.201305708.24449196
