Molecular Dynamics Simulation of Channel Size Dependence of the

Article pubs.acs.org/JPCC

Molecular Dynamics Simulation of Channel Size Dependence of the Friction Coefficient between a Water Droplet and a Nanochannel Wall Akinori Fukushima,*,† Toshiki Mima,‡ Ikuya Kinefuchi,‡ and Takashi Tokumasu§ †

Graduate School of Engineering, Tokyo University of Agriculture and Technology, 2-24-16 Nakamachi, Koganei, Tokyo 184-8588, Japan ‡ Department of Mechanical Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan § Institute of Fluid Science, Tohoku University, 2-1-1 Katahira, Aoba-ku, Sendai, Miyagi 980-8577, Japan ABSTRACT: The rapid spread of micro/nanoelectromechanical systems necessitates detailed understanding of fluidics within nanoscale structures. In this paper, the dynamics of a water droplet in nanochannels are analyzed using molecular dynamics simulations. As the channel size decreases, the shear stress between the droplet and the solid wall becomes much larger than predictions based on conventional slip boundary conditions. Our analysis shows that the Navier friction coefficient is quite sensitive to liquid pressure, which tends to be significantly large in hydrophobic nanochannels because of the Laplace pressure. We propose a modified version of the Young−Laplace equation that can accurately estimate the liquid pressure in nanochannels. By accounting for these nanochannel characteristics, we have successfully derived an expression that describes the channel size dependence of the shear stress between the droplet and the solid wall.

1. INTRODUCTION Recent advances in material science have enabled the development of new functional materials, which have subsequently stimulated numerous studies of devices constructed with these materials. Some of these materials have nanoscale structures. The transport properties of fluids through such small channels are different from those in macroscale channels,1 as can be seen in the unique flow fields within fuel cells. Because of high energy density and efficiency, fuel cells are candidate for next-generation energy systems. Polymer electrolyte fuel cells (PEFC) are particularly promising as energy sources in electric vehicles because of their low operating temperature and light weight.2−4 In this type of fuel cell, electrochemical reactions occur in a membrane electrode assembly (MEA), which consists of a gas diffusion layer, a microporous layer, a catalyst layer, and a polymer electrolyte membrane. Efficient proton transfer through the polymer electrolyte membrane requires adequate humidity inside the membrane. At low temperature and high humidity conditions, pores within the catalyst layer fill with condensed water, preventing oxygen gas diffusion. This phenomenon, called “flooding”, causes many serious problems such as a decline in the power generation efficiency and degradation of the membrane.5−12 Conversely, at high temperature and low humidity conditions, the deficiency of water reduces the ionic conductivity of the membrane and, as a result, decreases overall electric generation performance, a phenomenon known as “dry out”. Thus, water management that optimizes water distribution in MEAs is of critical importance for efficient operation of © 2015 American Chemical Society

PEFCs. The microporous layer is employed to prevent both of the phenomena and to maintain an efficient supply of oxygen gas.13−15 This layer has small pores with diameters between 10 and 100 nm. To control the amount of water in the membrane, not only vapor but also water droplets should be adequately exhausted outside the pores. To that end, it is important to understand droplet transport mechanisms, which remain unclear because of the large influence of surface tension in such small pores. Moreover, the catalyst layer has pores with diameters of the order of nanometers. Droplets in these pores are very small, comprising only a few thousand water molecules. Solid−liquid interfacial interactions almost entirely dominate the droplet characteristics. Additionally, the Laplace pressure is of the order of 10 MPa because of the significant curvature effect. These factors may effect on the water transport properties within nanochannels. However, details of these effects remain unknown. From both engineering and scientific aspects, insight into water transport phenomena in nanochannels is in high demand. Many studies on fluid transport through nanochannels have been published.16−19 Falk et al.16 evaluated the friction coefficient between water and carbon materials, such as graphene and carbon nanotubes (CNT), using molecular dynamics (MD) simulations, showing that the friction coefficient depends on the curvature of walls. Thomas et al.17 Received: August 15, 2015 Revised: November 11, 2015 Published: December 15, 2015 28396

DOI: 10.1021/acs.jpcc.5b07951 J. Phys. Chem. C 2015, 119, 28396−28404

Article

The Journal of Physical Chemistry C considered Poiseuille flow driven by pressure gradients within CNTs having diameters of several nanometers and compared the volumetric flow rate seen in the MD simulation with predictions of hydrodynamic equations. They found the volumetric flow rate to be larger than the predicted value because the slip length and viscosity depend on the tube diameter. Khademi et al.18 considered water flowing in a SiC nanotube under pressure gradient. They also showed that some water transport properties depend on the tube diameter, but this dependence differs from that in CNTs. The above-mentioned studies considered liquid dynamics in a channel whose entire region was filled with water. In order to deal with droplet dynamics, the effects of the liquid−vapor and vapor−solid interfaces should also be considered. The contact angle is one of the important characteristic quantities of these interfaces and has been intensively investigated.20−22 Thompson et al.20 considered two immiscible fluids confined between solid walls and evaluated the change of the contact angle due to the shear stress. They found that even if the channel size becomes a few nanometers, macroscopic equations can represent the shear stress dependence of the contact angle. Gentner et al.21,22 also evaluated the contact angle on both hydrophobic and hydrophilic walls under shear stress. They reported that the contact angle deviates from that predicted by hydrodynamic theory for large shear stress but that molecular kinetic theory reproduces the contact angle even in such conditions. In addition to the static properties of the contact line, dynamic properties of the contact line have also been examined and modeled.23−26 However, these studies focused only on contact line dynamics and did not consider overall droplet dynamics. Moreover, the effect of the channel size has not been clarified in detail. From such engineering and scientific demands, we have studied the channel size dependence of the drag force on a droplet that moves through a nanochannel by MD simulations.27−29 In the present paper, we analyze these results in detail and discuss the impact of the channel size on droplet transport properties in terms of the Navier friction coefficient. From these results, we obtain basic insight into the dynamics of droplets in nanochannels. The paper is organized as follows. Section 2 presents the macroscopic model of the dynamics of a droplet. In section 3, we describe the simulation systems and methods of this study. In section 4, we present our simulation results and discuss the channel size dependence of the Navier friction coefficient. We also present the modified Young− Laplace equation, which is applicable to microscopic channels. Section 5 summarizes our conclusions.

Figure 1. System considered in this study. Under a constant body force, f, a water droplet is moving rightward with a speed V between two parallel walls separated by a distance W. The droplet length, which is defined as the horizontal distance between the advancing and receding contact lines, is Wy. θL and θR are the receding and advancing contact angles, respectively.

velocity distribution is fully two-dimensional in the peripheral regions of the droplet. Hence, we divide the droplet into the Poiseuille flow region and the non-Poiseuille flow regions as shown in Figure 2.

Figure 2. Region classification of the droplet according to the velocity field. The velocity profile becomes parabolic in the central part of the droplet, called the Poiseuille flow region. The width of this region is Wyc. The velocity distribution becomes fully two-dimensional in the peripheral regions, called the non-Poiseuille flow regions. The width of these regions is Wyb. Fc is the shear force between the liquid and the walls in the Poiseuille flow region. FR and FL are the force between the liquid and the walls in the non-Poiseuille flow regions at the right and left sides of the droplet, respectively. The sum of these forces at the liquid−solid interface is balanced to the total body force F.

In the Poiseuille flow region, the velocity component in the y direction vy(z) is governed by

0=

d2vy dp +μ 2 +f dy dz

(1)

where μ is the viscosity, p is the pressure, and f is the body force on the droplet. (See the Appendix for a detailed derivation of this equation.) In this study, the velocity of the droplet is set at 50 m/s. When a droplet moves in nanometer order channels with such a high speed, Reynolds numbers are in 10−1 order. Thus, Peclet numbers are over 10, and the effect of diffusion can be ignored. The pressure gradient dp/dy remains constant within the Poiseuille flow region. The nonslip boundary condition is usually applied for flow analysis in large-scale channels. However, the relative velocity between the fluid and a solid wall significantly impacts the flow characteristics in nanoscale channels.30,31 We therefore employ a slip boundary condition32−35 on the top and bottom walls such that

2. DROPLET DYNAMICS MODEL BASED ON THE NAVIER−STOKES EQUATION Figure 1 shows a schematic illustration of the system considered in this study. A water droplet is located in a channel bounded by two parallel walls and is driven horizontally by a constant body force. After sufficient time, the system reaches a steady state with a constant droplet velocity. In the following, we derive a model equation that describes the shear force between the droplet and the walls at a steady state on the basis of the two-dimensional Navier−Stokes equation. Provided that the droplet is sufficiently long in the y direction (parallel to the walls), the fluid velocity distribution far from the droplet surfaces becomes independent of the y coordinate. In this region, the flow field is approximated as a Poiseuille flow with a slip velocity at the walls. In contrast, the

Ls 28397

∂vy ∂zz =±W /2

⎛ W⎞ = ∓vy ⎜ ± ⎟ ⎝ 2⎠

(2) DOI: 10.1021/acs.jpcc.5b07951 J. Phys. Chem. C 2015, 119, 28396−28404

Article

The Journal of Physical Chemistry C where Ls is a slip length and W is the channel size (the distance between the upper and the lower walls). The slip length Ls is usually constant in the analysis.31,32,35−37 Multiplying both sides of eq 2 by μ/L yields μ

∂vy ∂zz =±W /2

Ls =

⎛ W⎞ = ∓αvy ⎜ ± ⎟ ⎝ 2⎠

μ α

(3)

(4)

The constant α is called the Navier friction coefficient.30,31 Under a slip boundary condition, the shear stress on the fluid from the wall is given by the product of the Navier friction coefficient and the slip velocity. By solving eq 1 with the slip boundary condition and averaging the velocity over the z direction (see Appendix for the details), we can obtain the relationship between the velocity of the droplet motion V and the shear force on the droplet from the walls F as follows: F=

6αμ V (2WxWyc) + 2σlvWx(cos θL − cos θR ) − ΔPWxW 6μ + αW

(5)

where σlv is the surface tension coefficient, θL (θR) is the receding (advancing) contact angle, ΔP is the pressure loss due to friction forces on the walls in the non-Poiseuille flow regions (Figure 2), Wx is the channel length in the x direction, and Wyc is the length of the Poiseuille flow region. Given that 2WxWyc is the total area of the solid walls of the Poiseuille flow region, the coefficient of the first term of eq 5 represents the shear stress. Consequently, we obtain an expression that gives the shear stress as follows: τ=

6αμ V 6μ + αW

Figure 3. Molecular dynamics (MD) simulation systems. (a) Simulation system for the evaluation of the shear force between a droplet and the walls. (b) Schematic diagram of the structure of the solid wall. Solid atoms contained in two neighboring layers are shown. Filled and open circles represent solid atoms in the upper and lower layers, respectively. (c) Simulation system for the evaluation of the Navier friction coefficient. (d) Simulation system for the evaluation of the surface tension coefficient of the liquid−vapor interface.

(6)

⎧100 Å 3W < 100 Å Lz = ⎨ ⎩ 3W otherwise

This can be readily calculated using quantities evaluated by MD simulations. When two channels that have the different channel sizes and the same Navier’s friction coefficient are compared, the smaller channel has the larger slip velocity and the smaller velocity difference inside the fluid. However, by the effect of decreasing the channel size, the velocity gradient increases. Consequently, the shear stress increases with decreasing the channel size.

(7)

The periodic boundary condition is employed for all three directions. By changing the number of water molecules, five differentsized water droplets are represented in each channel. The average numbers of water molecules among five droplets are about 1056, 1760, 2484, and 3168 in channels with 20, 30, 40, and 50 Å, respectively. Details are shown in ref 27. In the following analysis, sampling for the Poiseuille flow region is performed within −5 Å < y < 5 Å, where the origin of the y coordinate is the center of the droplet. The SPC/E model39 is used to represent the water molecules. Coulombic interactions are evaluated using the smooth particle mesh Ewald method.40 The interaction between the solid wall and a water molecule is represented by the summation of the Lennard-Jones interactions between the wall atoms and the oxygen atom of the water molecule. The wall atoms interact mutually via harmonic oscillator potentials.

3. SIMULATION SYSTEM AND METHOD Figure 3a shows a snapshot of the simulation system. The channel is bounded by two flat slabs that have a similar structure to α-graphite (Figure 3b). These slabs are composed of five honeycomb layers each separated by 3.4 Å. The wettability of the solid wall is controlled by changing the bond length between nearest wall atoms, r1, rather than by modifying the pairwise interaction between the wall atoms and the water molecules. This method enables flexible control of the wettability while keeping the friction force between the droplet and the walls at a reasonable value. In the present study, the parameter r1 is set to 3 Å. This makes the contact angle about 150°, which corresponds to the typical value of carbon materials used in microporous layers of PEFCs.38 The channel size, which is defined as the distance between the surface layers of two slabs, is W = 20, 30, 40, or 50 Å. The volume of the simulation cell is 52 Å × 135 Å × Lz. Here, Lz is given by

E(rij) =

kn (rij − rn)2 2

(8)

where rij is the distance between the ith and jth atoms, and rn and kn (n = 1, 2, 3) are the equilibrium bond length and the force constant of the spring, respectively. As shown in Figure 28398

DOI: 10.1021/acs.jpcc.5b07951 J. Phys. Chem. C 2015, 119, 28396−28404

Article

The Journal of Physical Chemistry C 3b, r1 and k1 are the potential parameters for the covalent bond between nearest atoms in the same layer. Parameters with n = 2 or 3 represent interactions between atoms belonging to different layers. The mass of a wall atom is 120 g/mol. The initial geometry is built by randomly placing water molecules around the center of the channel. To make the water droplet move at a constant velocity relative to the wall, wall atoms in the outermost layers are moved along the y direction at the constant velocity Vw = −50 m/s, while the external force to balance the friction force is applied on the droplet. The external force is equally distributed to all of the water molecules. To keep the channel size constant during the simulation, we keep velocities in the z direction of atoms in the most outside layers zero. As a result, we can restrict the movement of the most outside layers and keep the channel size constant. To obtain a steady state, we carry out the simulation for 100 ps, controlling the temperature at 300 K by velocity scaling all water molecules and wall atoms except in the outermost layers. In this preparation run, the cutoff lengths of the Lennard-Jones interactions and the short-range component of the Coulomb interactions are 10 and 15 Å, respectively. Numerical integration of the equation of motion is carried out using the multi-time scale method.43,44 The time interval of the numerical integration is 1 fs except for the long-range component of Coulomb interactions, where a time interval of 10 fs is adopted to reduce the computational cost. After the system reaches a steady sate, the sampling simulation is carried out. During the sampling, we applied velocity scaling only for wall atoms in the second outermost layers of the slabs to control the temperature of the layers. To decrease errors in the numerical integration, the cutoff lengths for both LennardJones and shot-range Coulomb interactions are set at 15 Å. The time interval of the numerical integration for the long-range component of the Coulomb interactions is 5 fs, while that for other components is 1 fs. To evaluate the Navier friction coefficient between water and the solid wall, we perform additional simulations of the system without liquid−vapor interface (Figure 3c). The relative velocity between the liquid and the solid walls is the same as in the droplet case. After the preparative simulation for 100 ps to obtain the steady state, production runs are performed for 3 ns in total. The shear stress is evaluated using the method described in ref 27. The surface tension coefficient of the liquid−vapor interface, σlv, is calculated as the integrated value of the pressure distribution along the direction perpendicular to the flat liquid−vapor interface shown in Figure 3d.39 From this calculation, σlv is 59.3 × 10−3 N/m. The contact angle dependence on the bond length, r1, is shown in Figure 4. In this simulation, we use the simulation

model shown in Figure 3a. The procedure is shown in ref 38. A droplet (cylindrical filament) consisting of 4000 water molecules contacts on the lower wall of the channel whose size is set to 100 Å. The four different bond lengths (1.5, 2, 2.5, and 3) are employed. Figure 4 indicates the contact angle as a function of the bond length. The contact angle increases with increasing the bond length. When the bond length is 3 Å, the static contact angle is 149°. From Young’s equation (eq 9), we can obtain the difference between the surface tension coefficient of the liquid−solid interface, σls, and the surface tension coefficient of the liquid−solid interface, σsv. cos θ = −

σls − σsv σlv

(9)

The difference is 50.9 × 10−3 N/m.

4. RESULTS AND DISCUSSION 4.1. Navier Friction Coefficient. The Navier friction coefficient between water and the wall, α, is evaluated using the model shown in Figure 3c. The number of water molecules is chosen so that the normal stress on the wall becomes almost zero. The shear stress is 2.6 MPa when the average velocity of the liquid is 50 m/s relative to the wall. The velocity

Figure 5. Velocity distribution in the case that the pressure is highest. Red dots show MD results, and the black line shows the fitting curve.

distribution shown in Figure 5 is well approximated by a Poiseuille flow with a slip velocity: vy(z) = az 2 + b

(10)

The parameters a and b are determined to fit vy(z) of the MD result except in regions within 5 Å from the walls, where the liquid exhibits a layered structure. The slip velocity is defined as the difference between the velocity of water on the wall and the velocity of the wall and is given by ⎛ W⎞ Vslip = vy⎜ − ⎟ − Vw ⎝ 2⎠

(11)

Since the velocity field is symmetric about the centerline of the channel, the slip velocity (45 m/s) is evaluated only at the lower wall. From eq 3, the Navier friction coefficient, α, is evaluated to be 5.9 × 104 Pa/(m s). The coefficient of our hydrophobic surface is approximately 4 times larger than that of graphite, on which the friction force is known to be significantly small.37 4.2. Comparison of Predicted and Measured Shear Stresses. We examine the applicability of the droplet dynamics model based on the Navier−Stokes equation to droplets in nanochannels by comparing the shear stress predicted from eq 6 and that measured from MD simulations. The Navier friction

Figure 4. Contact angle of the water on the honeycomb wall as a function of the bond length r1. 28399

DOI: 10.1021/acs.jpcc.5b07951 J. Phys. Chem. C 2015, 119, 28396−28404

Article

The Journal of Physical Chemistry C coefficient obtained in section 4.1, and the viscosity of SPC/E water at 300 K (6.33 × 10−4 Pa·s),45 are substituted for α and μ in eq 6. Figure 6 shows the channel size dependence of the

Figure 7. Pressure dependence of the Navier friction coefficient between water and the solid wall. Filled circles represent the results of MD simulations. The dashed line is the least-squares line (eq 12). Figure 6. Channel size dependence of the shear stress between a droplet and solid walls. Filled circles represent the shear stress measured from MD simulations.24 The dot-and-dash curve is the predicted value from the droplet dynamics model based on the Navier−Stokes equation (eq 6), the solid curve is the predicted value from eq 18, and the dashed curve is the predicted value from eq 18 where the viscosity coefficient value is set at half of the bulk value. (a) Comparison of the expected values from eqs 6 and 18. (b) Comparison of MD results with expected values from eqs 6 and 18.

nanochannel becomes much higher than that of the gas phase due to the Laplace pressure, which is given by σ Pliq − Pgas = lv (13) r where Pliq and Pgas are the pressures of the liquid and gas phases, and r is the radius of curvature of the liquid−vapor interface. The difference between the advancing and receding contact angles is only about 4°; thus, the radii of curvature at both sides of the droplet are regarded as being identical in the following analysis. Because the shear force at the walls is almost balanced to the body force, the pressure inside the droplet is nearly uniform and can be estimated by eq 12. From a geometrical consideration, the radius of curvature is expressed by W r=− (14) 2 cos θ

shear stress. The chain line in Figure 6 indicates the value from eq 6. Since the channel size in this study is in nanometers order, the term αW is in 10−5 order. On the other hand, 6μ is in 10−3 order. As a result, the shear stress is almost independent of the channel size. Figure 6 also shows the comparison between MD results and the expected values from eq 6. This indicates that the model significantly underestimates the shear stress obtained by MD simulations. Such underestimation may be attributed to a channel size dependence of at least one of the Navier friction coefficients or the viscosity. First, we examine the channel size dependence of the Navier friction coefficient in detail. We then discuss the influence of the viscosity variation. 4.3. Pressure Dependence of the Navier Friction Coefficient. When a water droplet is formed in a nanochannel with a hydrophobic surface, the pressure inside the droplet can increase up to approximately 10 MPa. This is because of the surface tension effects, which play an insignificant role in transport properties within millimeter-scale channels.25 Such an increase in pressure would markedly affect the Navier friction coefficient, but its pressure dependence has not been well clarified in previous studies.25 Thus, we examined the pressure dependence of the friction coefficient using five simulation models with different pressures, as shown in Figure 3c. Namely, in this evaluation, we use the system that does not have a liquid−vapor interface. The obtained pressure is different from the pressure in the droplet case. The pressure is controlled by changing the number of water molecules in the channels with the same width (W = 50 Å). The Navier friction coefficient increases as the pressure increases (Figure 7), showing a good linear correlation. α(P) = αsP + αi

Combining eqs 13 and 14 yields 2 Pliq − Pgas = (σls − σsv) (15) W The left-hand side can be approximated as Pliq since Pliq ≫ Pgas in our system. As reported in ref 27, when the channel size is less than 50 Å, eq 15 underestimates the actual pressure inside the droplet. In the following discussion, we clarify the origin of this discrepancy and present a modified version of eq 13, which produces better pressure estimations. To check the validity of eq 13, we evaluated the average pressure in the Poiseuille flow region and the curvature radius of the droplet surface in systems with different channel sizes. The curvature radius is determined from the mass density distribution of water molecules averaged over the x direction. For each z coordinate, we detect the position of the liquid− vapor interface as the inflection point of a hyperbolic tangenttype function fitted to the mass density distribution along the y direction. We then obtained the curvature radius of the droplet surface by fitting a circle to these inflection points. As shown in Figure 8, the dependence of the pressure on the curvature radius is well reproduced by eq 13, where the surface tension coefficient is assumed to be the same as for the flat interface. In eq 14, the cosine of the contact angle cos θ is calculated from the curvature radius r and the channel size W. As shown in Figure 9, eq 14 gives a value of cos θ that is smaller than −1 when the channel size is 40 Å or less. This unreasonable result is attributable to the excluded thickness δ at the solid−liquid interface, where water molecules are not able to enter because of the strong repulsive force between water molecules and wall atoms. The channel size W in eq 14 should be replaced with the

(12) −3

where αs and αi are determined to be 1.5 × 10 s/m and 5.4 × 104 Pa/(m s) by the least-squares method. The Navier friction coefficient can be estimated correctly using this correlation, provided that we know the channel size dependence of the pressure inside the droplet, which is discussed in the following section. 4.4. Channel Size Dependence of Pressure inside a Droplet. The pressure inside a water droplet in a hydrophobic 28400

DOI: 10.1021/acs.jpcc.5b07951 J. Phys. Chem. C 2015, 119, 28396−28404

Article

The Journal of Physical Chemistry C

shows the excluded thickness as a function of the channel size. The excluded thickness slightly decreases as the channel size decreases. This is because the increase of the Laplace pressure in the narrower channel enables water molecules to move much closer to the walls. However, the difference in the excluded thicknesses among the channels considered here (0.2 Å at most) is so small that we approximate the excluded width as constant (δ = 5.1 Å). We evaluated the cosine of the contact angle again by replacing W in eq 14 with W − δ. As shown in Figure 8, the correction of the channel size yields a value of cos θ larger than −1. The modified version of eq 15 becomes 2 2 1 Pliq − Pgas = (σls − σsv) = (σls − σsv) δ W−δ W 1−

Figure 8. Interior droplet pressure plotted against the curvature radius of the liquid−vapor interface. Red markers show the result of MD simulations. The solid curve shows the pressure estimated from eq 13.

W

(16)

When the channel size is large enough, this equation agrees with eq 15. Figure 11 compares the average pressure in the

Figure 9. Channel size dependence of the cosine of the contact angle cos θ. Red markers are the estimated values of cos θ using eq 14. Green markers are also the estimated values using eq 14, but the channel size W is replaced with the effective channel size W − δ. Figure 11. Pressure inside a water droplet in different-sized channels. Red markers are the MD simulation results.24 The dashed and solid curves are the estimated pressures using eqs 15 and 16, respectively.

effective channel size W  δ to calculate cos θ correctly in nanochannels. Figure 10a shows the mass density distribution of water molecules along the z direction in the Poiseuille flow region when the channel size W is 50 Å. The positions of the liquid−solid interfaces are defined as the inflection points of hyperbolic tangent type functions fitted to the density distribution. The effective channel size is calculated as the distance between the upper and lower interfaces. Figure 10b

Poiseuille flow region obtained from the MD simulations26 and the estimated pressures by eqs 15 and 16. Equation 16 successfully reproduces the pressure obtained from the MD simulations, indicating that if we take into account the excluded thickness, the Young−Laplace equation still holds in nanochannels. 4.5. Channel Size Dependence of the Shear Stress Including the Pressure Effect. Using eqs 12 and 16, the Navier friction coefficient is expressed as a function of the channel size: α (W ) =

αi(W − δ) + 2αs(σls − σsv) W−δ

(17)

Substituting this result to eq 6 yields τ (W ) =

6μ[αi(W − δ) + 2αs(σls − σsv)]V 6μ(W − δ) + W [αi(W − δ) + 2αs(σls − σsv)] (18)

which gives the channel size dependence of the shear stress between the wall and the liquid. As shown in Figure 6, the modified expression of the shear stress (eq 18) predicts the shear stress evaluated from the MD simulations much better than the original expression (eq 6). The shear stress increases as the channel size decreases because the Laplace pressure (approximately 107 Pa) in the nanochannels noticeably increases the Navier friction coefficient. The difference between the predicted shear stresses by eqs 6 and 18 is about 10%, even when the channel size is 200 Å. The difference becomes much more pronounced for narrower channels. In such cases, eq 18, which accounts for the contribution of the excluded width in

Figure 10. (a) Density profile of the liquid along the direction normal to the wall. Red markers show the results of the MD simulations. The solid curve is a hyperbolic tangent type function fitted to the mass density distribution of the MD result. (b) Exclusion width between the liquid and the solid wall as a function of the channel size W. 28401

DOI: 10.1021/acs.jpcc.5b07951 J. Phys. Chem. C 2015, 119, 28396−28404

Article

The Journal of Physical Chemistry C

in the nanoscale channel or not. Thus, the assumptions used in the macroscopic system are employed. At the steady state, the following relations hold in the Poiseuille flow region:

the Laplace pressure and the pressure dependence of the Navier friction coefficient, gives a quantitative estimation of the shear stress between the liquid and the solid wall.

rn [Å] kn [N/m2]

2

3

3.00 1000

3.40 100

4.53 10

∂ν2 ∂ν = 2 = 0, ∂x1 ∂x 2

∂p ∂p = =0 ∂x1 ∂x3

n 1

∂ν2 = 0, ∂t

ν1 = ν3 = 0,

Table 1. Parameters of the Harmonic Oscillator Potential Model between Wall Atoms

(A2)

Since the surface tension does not appear in the Poiseuille flow region, we have

δs = 0

4.6. Viscosity Dependence of the Shear Stress. Finally, we focus on the viscosity dependence of the shear stress in eq 18 because the effective viscosity in nanochannels can be different from the bulk value. For example, the viscosity of water in a CNT with a diameter of several nanometers is about half of the bulk viscosity.17 To examine how the change of the viscosity impacts on the shear stress, the shear stress is evaluated, assuming the viscosity becomes half of the bulk viscosity (the dashed curve in Figure 6). The result indicates that the shear stress has relatively little dependence on the viscosity.

(A3)

The body force is parallel to the y direction and thus satisfies

f1 = f3 = 0,

f2 = f

(A4)

Using eqs A2−A4, the Navier−Stokes equations (eq A1) reduce to eq 1. The pressure gradient is evaluated as

dp P − PR ≈ L dy Wyc

(A5)

where PL and PR are the pressures at the left and right edges of the Poiseuille flow region. From eqs 13 and 14, these pressures are given by

5. CONCLUSIONS MD simulations are conducted to analyze the transport properties of water droplets in nanoscale channels in terms of the shear stress between the droplet and the solid walls. We find that the macroscopic equation is insufficient to estimate the friction force between the droplet and solid walls nanoscale channels. When the channel size is on the order of nanometers, the Navier friction coefficient, which is usually assumed constant, increases because of the large Laplace pressure. As a result, the shear stress between the droplet and the wall becomes much higher than that predicted by conventional slip boundary conditions. When the channel size is less than 40 Å, the exclusion width, which reduces the effective channel size, must be taken into account to estimate the curvature radius of the liquid−vapor interface. The Young−Laplace equation is valid even if the channel size is as small as 20 Å. An accurate expression, that reproduces the MD results well, for the shear stress between a droplet and the wall in hydrophobic nanochannels has been derived by considering these factors.

PL − Pgas = −

2σlv cos θL − ΔPL W

(A6)

PR − Pgas = −

2σlv cos θR + ΔPR W

(A7)

where ΔPL and ΔPR are the pressure drops within the nonPoiseuille flow regions at the left and right sides of the droplet. Combining eqs A5−A7 yields dp 2σlv ΔP ≈ (cos θL − cos θR ) − dy WWyc Wyc

(A8)

where ΔP = ΔPR + ΔPL. Substituting eq A8 to eq 1 and integrating eq 1 with the boundary condition eq 3, we obtain the velocity distribution: vy (z) = −

■

2σlv ΔP ⎤⎥ 1 ⎡⎢ f− (cos θL − cos θR ) + 2μ ⎢⎣ WWyc Wyc ⎥⎦

⎛ μ⎞ W2 × ⎜z 2 − −W ⎟ α⎠ 4 ⎝

APPENDIX The Navier−Stokes equations that include the surface tension term are written as

(A9)

Given that the mean velocity of the droplet is V

⎡ ⎛ ∂v ⎞⎤ ⎡ ∂v ⎤ j ∂vi ⎟⎟⎥ + 2σlvniδs + f ρ⎢ i + (vj·∇j )vi ⎥ = −∇i p + ∇j ⎢μ⎜⎜ i ⎣ ∂t ⎦ ⎣⎢ ⎝ ∂xi ∂xj ⎠⎥⎦

1 W

(A1)

W /2

∫−W /2 vy(z) dz = V

(A10)

From eqs A9 and A10, the body force is derived as follows:

where the subscripts i and j refer to the direction, vi is the velocity, ρ is the density, p is the pressure, σlv is the surface tension coefficient, and δs and ni are the δ-function and the normal vector on the liquid−vapor interface, respectively. As shown in Figure 2, the forces between the droplet and the solid wall are divided into those acting in the non-Poiseuille flow region and those acting in Poiseuille flow regions. The former is the sum of the force due to surface tension and the shear force generated by the two-dimensional flow field and is independent of the droplet length Wy. In this study, we focus on whether we can use the equation based on the macroscopic viewpoint also

f=

α W 2

+

α W2 μ 12

V+

2σlv ΔP (cos θL − cos θR ) − WWyc Wyc (A11)

The total force acting on the Poiseuille flow region is therefore Fc = fWxWycW =

6αμ V ·2WxWyc 6μ + αW

+ 2σlv × Wx(cos θL − cos θR ) − ΔPWxW 28402

(A12)

DOI: 10.1021/acs.jpcc.5b07951 J. Phys. Chem. C 2015, 119, 28396−28404

Article

The Journal of Physical Chemistry C

■

(18) Khademi, M.; Sahimi, M. Molecular Dynamics Simulation of Pressure-Driven Water Flow in Silicon-Carbide Nanotubes. J. Chem. Phys. 2011, 135, 204509. (19) Holt, J. K.; Park, H. G.; Wang, Y.; Stadermann, M.; Artyukhin, A. B.; Grigoropoulos, C. P.; Noy, A.; Bakajin, O. Fast Mass Transport through Sub-2-Nanometer Carbon Nanotubes. Science 2006, 312, 1034−1037. (20) Thompson, P. A.; Robbins, M. O. Simulations of Contact-Line Motion: Slip and the Dynamic Contact Angle. Phys. Rev. Lett. 1989, 63, 766. (21) Gentner, F.; Ogonowski, G.; De Coninck, J. Forced Wetting Dynamics: A Molecular Dynamics Study. Langmuir 2003, 19, 3996− 4003. (22) Werder, T.; Walther, J. H.; Jaffe, R. L.; Halicioglu, T.; Noca, F.; Koumoutsakos, P. Molecular Dynamics Simulation of Contact Angles of Water Droplets in Carbon Nanotubes. Nano Lett. 2001, 1, 697− 702. (23) Raj, R.; Maroo, S. C.; Wang, E. N. Wettability of Graphene. Nano Lett. 2013, 13, 1509−1515. (24) Ren, W.; Weinan, E. Boundary Conditions for the Moving Contact Line Problem. Phys. Fluids 2007, 19, 022101. (25) Ren, W.; Hu, D.; Weinan, E. Continuum Models for the Contact Line Problem. Phys. Fluids 2010, 22, 102103. (26) Bertrand, E.; Blake, T. D.; De Coninck, J. Influence of Solid− Liquid Interactions on Dynamic Wetting: A Molecular Dynamics Study. J. Phys.: Condens. Matter 2009, 21, 464124. (27) Fukushima, A.; Mima, T.; Kinefuchi, I.; Tokumasu, T. Molecular Dynamics Study for Channel Size Dependence of Shear Stress between Droplet and Wall. J. Nanosci. Nanotechnol. 2015, 15, 3224− 3228. (28) Fukushima, A.; Mima, T.; Kinefuchi, I.; Tokumasu, T. Molecular Dynamics Study for the Friction Coefficient between a Water Droplet and a Solid Wall. ECS Trans. 2013, 58, 87−93. (29) Fukushima, A.; Mima, T.; Kinefuchi, I.; Tokumasu, T. Molecular Dynamics Study of Water Transport Property in Micro Hydrophobic Pore. ECS Trans. 2013, 50, 197−206. (30) Joseph, P.; Tabeling, P. Direct Measurement of the Apparent Slip Length. Phys. Rev. E 2005, 71, 035303. (31) Hansen, J. S.; Todd, B.; Daivis, P. J. Prediction of Fluid Velocity Slip at Solid Surfaces. Phys. Rev. E 2011, 84, 016313. (32) Schnell, E. Slippage of Water over Nonwettable Surfaces. J. Appl. Phys. 1956, 27, 1149−1152. (33) Snoeijer, J. H.; Andreotti, B. Moving Contact Lines: Scales, Regimes, and Dynamical Transitions. Annu. Rev. Fluid Mech. 2013, 45, 269−292. (34) Thompson, P. A.; Troian, S. M. A General Boundary Condition for Liquid Flow at Solid Surfaces. Nature 1997, 389, 360−362. (35) Priezjev, N. V.; Darhuber, A. A.; Troian, S. M. Slip Behavior in Liquid Films on Surfaces of Patterned Wettability: Comparison between Continuum and Molecular Dynamics Simulations. Phys. Rev. E 2005, 71, 041608. (36) Werder, T.; Walther, J. H.; Jaffe, R.; Halicioglu, T.; Koumoutsakos, P. On the Water-Carbon Interaction for Use in Molecular Dynamics Simulations of Graphite and Carbon Nanotubes. J. Phys. Chem. B 2003, 107, 1345−1352. (37) Kannam, S. K.; Todd, B.; Hansen, J. S.; Daivis, P. J. Slip Length of Water on Graphene: Limitations of Non-Equilibrium Molecular Dynamics Simulations. J. Chem. Phys. 2012, 136, 024705. (38) Wang, X.; Zhang, H.; Zhang, J.; Xu, H.; Tian, Z.; Chen, J.; Zhong, H.; Liang, Y.; Yi, B. Micro-Porous Layer with Composite Carbon Black for Pem Fuel Cells. Electrochim. Acta 2006, 51, 4909− 4915. (39) Alejandre, J.; Tildesley, D. J.; Chapela, G. A. Molecular Dynamics Simulation of the Orthobaric Densities and Surface Tension of Water. J. Chem. Phys. 1995, 102, 4574−4583. (40) Scocchi, G.; Sergi, D.; D’Angelo, C.; Ortona, A. Wetting and Contact-Line Effects for Spherical and Cylindrical Droplets on Graphene Layers: A Comparative Molecular-Dynamics Investigation. Phys. Rev. E 2011, 84, 061602.

AUTHOR INFORMATION

Corresponding Author

*E-mail [email protected] (A.F.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by New Energy and Industrial Technology Development Organization and Japan Society for the Promotion of Science KAKENHI Grants H2248214 and 15K17968. All the simulations were carried out using the supercomputer system at the Advanced Fluid Information Research Center, Institute of Fluid Science, Tohoku University.

■

REFERENCES

(1) Karniadakis, G.; Beskok, A.; Aluru, N. Microflows and Nanoflows: Fundamentals and Simulation; Springer Science & Business Media: 2006; Vol. 29. (2) Lemons, R. A. Fuel Cells for Transportation. J. Power Sources 1990, 29, 251−264. (3) Shan, Y.; Choe, S.-Y. A High Dynamic Pem Fuel Cell Model with Temperature Effects. J. Power Sources 2005, 145, 30−39. (4) Wang, Y.; Chen, K. S.; Mishler, J.; Cho, S. C.; Adroher, X. C. A Review of Polymer Electrolyte Membrane Fuel Cells: Technology, Applications, and Needs on Fundamental Research. Appl. Energy 2011, 88, 981−1007. (5) Yang, X. G.; Zhang, F. Y.; Lubawy, A. L.; Wang, C. Y. Visualization of Liquid Water Transport in a PEFC. Electrochem. SolidState Lett. 2004, 7, A408. (6) Pasaogullari, U.; Wang, C.-Y. Two-Phase Modeling and Flooding Prediction of Polymer Electrolyte Fuel Cells. J. Electrochem. Soc. 2005, 152, A380. (7) Ge, S.; Wang, C.-Y. Liquid Water Formation and Transport in the Pefc Anode. J. Electrochem. Soc. 2007, 154, B998. (8) Sinha, P. K.; Wang, C.-Y. Pore-Network Modeling of Liquid Water Transport in Gas Diffusion Layer of a Polymer Electrolyte Fuel Cell. Electrochim. Acta 2007, 52, 7936−7945. (9) Wang, X.-D.; Duan, Y.-Y.; Yan, W.-M. Novel Serpentine-Baffle Flow Field Design for Proton Exchange Membrane Fuel Cells. J. Power Sources 2007, 173, 210−221. (10) Wang, X.-D.; Duan, Y.-Y.; Yan, W.-M.; Peng, X.-F. Local Transport Phenomena and Cell Performance of Pem Fuel Cells with Various Serpentine Flow Field Designs. J. Power Sources 2008, 175, 397−407. (11) Li, H.; Tang, Y.; Wang, Z.; Shi, Z.; Wu, S.; Song, D.; Zhang, J.; Fatih, K.; Zhang, J.; Wang, H.; Liu, Z.; Abouatallah, R.; Mazza, A. A Review of Water Flooding Issues in the Proton Exchange Membrane Fuel Cell. J. Power Sources 2008, 178, 103−117. (12) Nam, J. H.; Lee, K. J.; Hwang, G. S.; Kim, C. J.; Kaviany, M. Microporous Layer for Water Morphology Control in PEMFC. Int. J. Heat Mass Transfer 2009, 52, 2779−2791. (13) Wang, Z.; Wang, C.; Chen, K. Two-Phase Flow and Transport in the Air Cathode of Proton Exchange Membrane Fuel Cells. J. Power Sources 2001, 94, 40−50. (14) Weber, A. Z.; Newman, J. Effects of Microporous Layers in Polymer Electrolyte Fuel Cells. J. Electrochem. Soc. 2005, 152, A677− A688. (15) Gostick, J. T.; Ioannidis, M. A.; Fowler, M. W.; Pritzker, M. D. On the Role of the Microporous Layer in PEMFC Operation. Electrochem. Commun. 2009, 11, 576−579. (16) Falk, K.; Sedlmeier, F.; Joly, L.; Netz, R. R.; Bocquet, L. Molecular Origin of Fast Water Transport in Carbon Nanotube Membranes: Superlubricity Versus Curvature Dependent Friction. Nano Lett. 2010, 10, 4067−4073. (17) Thomas, J. A.; McGaughey, A. J.; Kuter-Arnebeck, O. PressureDriven Water Flow through Carbon Nanotubes: Insights from Molecular Dynamics Simulation. Int. J. Therm. Sci. 2010, 49, 281−289. 28403

DOI: 10.1021/acs.jpcc.5b07951 J. Phys. Chem. C 2015, 119, 28396−28404

Article

The Journal of Physical Chemistry C (41) Berendsen, H.; Grigera, J.; Straatsma, T. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. (42) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (43) Dullweber, A.; Leimkuhler, B.; McLachlan, R. Symplectic Splitting Methods for Rigid Body Molecular Dynamics. J. Chem. Phys. 1997, 107, 5840−5851. (44) Omelyan, I. P. Advanced Gradientlike Methods for Rigid-Body Molecular Dynamics. J. Chem. Phys. 2007, 127, 044102. (45) Medina, J. S.; Prosmiti, R.; Villarreal, P.; Delgado-Barrio, G.; Winter, G.; González, B.; Alemán, J. V.; Collado, C. Molecular Dynamics Simulations of Rigid and Flexible Water Models: Temperature Dependence of Viscosity. Chem. Phys. 2011, 388, 9−18.

28404

DOI: 10.1021/acs.jpcc.5b07951 J. Phys. Chem. C 2015, 119, 28396−28404