Molecular dynamics simulation study on the mechanisms of liquid-phase permeation in nanopores

Separation and Purification Technology 220 (2019) 259–267

Contents lists available at ScienceDirect

Separation and Purification Technology journal homepage: www.elsevier.com/locate/seppur

Molecular dynamics simulation study on the mechanisms of liquid-phase permeation in nanopores

T

⁎

Tomohisa Yoshiokaa, , Rina Kunimorib, Ikumi Hisaokab, Hiroki Nagasawab, Masakoto Kanezashib, Toshinori Tsurub a b

Center for Membrane and Film Technology, Graduate School of Science, Technology and Innovation, Kobe University, 1-1, Rokkodai, Nada, Kobe 657-8501, Japan Department of Chemical Engineering, Hiroshima University, 1-4-1 Kagami-yama, Higashi-Hiroshima 739-8527, Japan

A R T I C LE I N FO

A B S T R A C T

Keywords: Ceramic membrane Molecular dynamics Nanofiltration Permeation mechanism

Liquid-phase permeation was simulated in nano-scale pores via non-equilibrium molecular dynamics (NEMD). Two virtual cristobalite membranes were prepared with pore diameters of 1.7 and 2.4 nm. NEMD simulation system was employed as an ideal experimental system to calculate the affinity between liquid argon molecules and membrane materials during permeation. When argon-membrane interactions decreased, permeation flux increased. With a smaller interaction the permeation flux exceeded the value posited by the Hagen-Poiseuille theorem, while a lower-than-expected level of permeation flux was observed when the interactions with the pore surface became greater. We focused on the viscosity change of liquid in a nano-scale pore due to attractive or repulsive interactions with the pore surface, and a mathematical model for describing the liquid permeation flux in a nanopore was proposed by solving the Navier-Stokes equation by considering the viscosity distribution of a liquid confined in a pore. The local viscosity of a liquid confined in a pore was calculated from the total potential distribution in the pore based on the Andrade equation. The predicted level of permeation flux, the velocity profiles of different pore sizes, and the interactions of the pore models all showed good agreement with the NEMD simulation results.

1. Introduction Due to recent industry developments and to a rapid increase in the populations of emerging nations, maintaining water resources has become an important global issue. To meet the demand for the development of highly efficient water purification processes, commercial membrane water treatment and separation technologies for the desalination of seawater and filtration of inland water is receiving everincreasing amounts of attention with respect to features that save energy and reduce cost. Nano-filtration membranes with nanometer pore sizes that range in the single digits are expected to be applied to solvent separation and ion separation processes in addition to high-grade water treatment. Researchers have found, however, that the permeation flux of liquids in microporous nanofiltration membranes would deviate from the values predicted by the ordinal Hagen-Poiseuille theorem, and the rapid decline in water flux that occurs with decreases in pore size is a serious problem for nanoporous ceramic membranes [1,2]. On the contrary, permeation within carbon nanotubes has been reported with water flux that is more than 10,000-fold higher than that in nanoporous ceramic membranes [3]. If the deviation of a permeation flux from the ⁎

estimated value by the Hagen-Poiseuille flow of an ordinal continuous bulk liquid can be predicted quantitatively based on the interactions of a permeating liquid with the surface and structures of nanopores, that would be useful and indispensable in the design and development of high-performance microporous nanofiltration membranes. Extremely high levels of water flux in carbon tubes has been energetically studied under both experimental examination and molecular dynamics simulation [4–10]. The most popular mathematical model for describing enhanced water flux is the slip length model [11]. This model can explain high levels of water flux by introducing the concept of slip length in the vicinity of a pore wall. The slip length of the model Lz in the direction of flow, z, is defined by Eq. (1), and the velocity profile, uz(r), as a function of location, r, in the radial direction in a cylinder with radius, R, with pressure gradient, dp/dz, is given by Eq. (2), which is derived by solving the Navier-Stokes equation using Lz, where μ is the viscosity of a permeating liquid.

Ls =

uz (r ) duz /dr

Corresponding author. E-mail address: [email protected] (T. Yoshioka).

https://doi.org/10.1016/j.seppur.2019.03.014 Received 31 July 2018; Received in revised form 19 February 2019; Accepted 5 March 2019 Available online 07 March 2019 1383-5866/ © 2019 Elsevier B.V. All rights reserved.

r=R

(1)

Separation and Purification Technology 220 (2019) 259–267

T. Yoshioka, et al.

ideal experimental system.

uz (m/s)

2. Theory

Hagen-Poiseuille The momentum of the transport in the z- to r-direction in a cylinder radius, R, under a pressure gradient, dp/dz, in the z-direction is given by Eq. (3) using Newton’s law of viscosity for the Navier-Stokes equation with elimination of the unnecessary terms for simplification.

uz = 0

r (nm) r = -R

∂p ∂ ⎛ ∂uz ⎞ rμ = r ∂z ∂r ⎝ ∂r ⎠

LS r=R

Eq. (3) can be easily solved by applying the following boundary conditions with an assumption of constant viscosity, μ, in a cylindrical pore, which leads to a parabolic-shaped velocity distribution in the zdirection as a function of r (Eq. (5)).

Fig. 1. Conceptual diagram of slip boundary conditions.

uz (r ) =

R2 ⎡ r 2 2Ls ⎤ ∂p 1−⎛ ⎞ + 4μ ⎢ R R ⎥ ⎝ ⎠ ⎦ ∂z ⎣

(3)

∂uz (r = 0) = 0uz (r = R) = 0 ∂r

(2)

uz (r ) =

As shown in Fig. 1, Eq. (2) can be drawn as a parabolic curve with a finite value of uz on the pore wall (at r = R) that expresses the enhancement of permeation flux with regards to the predicted flux according to the Hagen-Poiseuille theorem. However, the physical meaning of slip length is as yet unclear, in addition, it remains difficult to explain the unexpected lower levels of water flux in the nanopores of ceramic membranes. To solve the above problem, in 2010 Myers proposed a mathematical model by solving the Navier-Stokes equation using two different viscosities of a liquid: the first was for a cylindrical pore in the main bulk region around the pore center (1), viscosity μ1; and, the second was for a limited annulus ring region next to the pore surface (2), viscosity μ2 [12]. Although an artificial assumption of phenomenologically curious discontinuous change in the viscosity of an identical liquid persisted in this model, the physical meaning of the model was clear and was consistent with the deviation of density and diffusivity coefficients from the bulk liquid properties calculated by MD simulations [13,14]. Bonthuis et al. proposed a model for the behaviors of ions and water molecules in the vicinity of solid surface by taking into account of the effect of electrostatic and hydrodynamic forces [15,16]. They introduced an increased fluid viscosity for the case of hydrophilic solid surface and a decreased viscosity for hydrophobic one to the hydrodynamic term in the model. In their system, since the viscosity drastically changed with approaching to the solid surface, the original continuous viscosity profile was approximately treated as two values, one was within a thin layer whose thickness of water molecular-scale on the solid surface and the other was for the boundary beyond. The thickness of the boundary layer was 0.4 nm for hydrophilic solid and 0.15 nm for hydrophobic one. The model could successfully describe the apparent excess surface conductivity. In our system, however, there was no electrostatic interaction, so the distinction of boundary layer on the pore solid surface and other region was not employed. The feature of our model was introduction of a continuous viscosity distribution along with radial direction according to the interaction from the pore wall for estimation of liquid permeability in a nano-scale pore. The present work presents a mathematical transport model for a liquid confined in a nanopore with a continuous viscosity profile as a function of location in the radial direction in a cylindrical pore based on the potential distribution in the pore. Liquid-phase permeation simulations were also conducted using non-equilibrium molecular dynamics (NEMD) simulation to obtain molecular scale permeation characteristics and macroscopic permeation flux for two different cylindrical pore sizes (1.7 and 2.4 nm), and for several repulsive and attractive interactive conditions with a pore surface. The validity of the permeation model was established by comparing the predicted velocity profiles in nanopores and the permeation flux for different conditions with MD simulation results as an

R2 ⎡ r 2 ∂p 1−⎛ ⎞⎤ ⎢ 4μconst. ⎣ ⎝R⎠ ⎥ ⎦ ∂z

(4)

(5)

The volumetric flow rate, QH.P. can be described as Eq. (6) by integrating Eq. (5). Eq. (6) is well known as the Hagen-Poiseuille equation [17].

QH.P. =

π ΔpR 4 8μconst.l

(6)

In general, when a permeating fluid can be treated as a continuous and homogeneous fluid, transport momentum can be expressed using Eq. (3). This is often true for liquid permeation in relatively large pores in the micro and sub-micro meters scale, and a permeating liquid flow rate can be correctly estimated using the Hagen-Poiseuille equation [18]. Concerning the liquid state of the flow rate through microporous membranes with single nanometer-scale pores, the deviations of the flow rate from the value predicted by Eq. (6) were experimentally observed [1,3]. We assumed that the deviation could be attributed to a change in the physical properties of a fluid in the vicinity of a pore wall where a confined liquid would have a particular orientation, density, and diffusivity that would differ from that of a bulk liquid due to interactions with the molecules that compose the nanopore [13,14]. The liquid volume ratio of a portion under the influence of an interaction with the pore wall to the total fluid volume in the nanopore is relatively high to the extent of being significant. Ortiz-Young et al. reported that the viscous shear forces of confined water in the nanospace near hydrophilic surfaces were greater than those of bulk water [19]. Therefore, to study the liquid permeation mechanisms of viscous flow in nanopores, inhomogeneous viscosity profiles of nano-confined fluid that are caused by interactions with the pore wall should be considered. When viscosity, μ, is not constant, and can change as a function of the distance from the pore center, r, Eq. (3) is modified to Eq. (7) using μ(r). The local velocity at r can be easily derived by integrating Eq. (7) with the boundary conditions in Eq. (4) to give Eq. (8).

∂uz ⎞ ∂p ∂ ⎛ rμ (r ) r = ∂z ∂r ⎝ ∂r ⎠ uz (r ) =

pu − pd 2l

∫0

d /2

(7)

r dr μ (r )

(8)

3. MD simulation models and methods In order to investigate fine-scale phenomena such as high-density fluid permeation in nanopores, molecular dynamics simulation is very useful for both molecular-scale direct observation and quantitative minute analysis. We adopted a non-equilibrium molecular dynamics permeation simulation of liquid argon as an ideal experimental system 260

Separation and Purification Technology 220 (2019) 259–267

T. Yoshioka, et al.

of nanofiltration in order to verify the above liquid permeation model in nanopores. A virtual membrane was modeled as a cristobalite silica structure with the dimensions of 4.3 × 4.3 × 8.6 nm3 with one penetrating cylindrical nanopore [20]. An unit cube composed of eight tetrahedral structure of SiO4/2 whose SieO bond length was 0.164 nm, OeSieO angle was 109.47°, and SieOeSi angle was 180.0° was repeated by 6 times in the x- and y-direction and 12 times in the z-direction, respectively, to obtain a cuboid-shaped β- cristobalite SiO2 structure of 4.3 × 4.3 × 8.6 nm3. Si and O atoms whose center were located in a cylinder from the center of (x, y) = (0, 0) to radius, r = 0.85 nm and 1.2 nm were removed to form an artificial cylindrical pore models. The diameters were d = 1.7 nm and 2.4 nm, respectively. Due to this preparation method of cylindrical pores, the pore surfaces were not smooth but had unevenness in the oxygen atomic size. Since a real ceramic membrane has pore size distribution and the pore shape is different from an ideal cylinder, actual liquid permeation phenomena could not be simulated by the presented membrane model. The ideal model was adopted in this work to study simply the effect of interaction between permeating liquid molecules and the pore surface on liquid permeation properties. The length of the model pore was 8 nm, and this straight channel part would be a model of a narrow inter-particle neck with the highest resistance for permeation. Fig. 2 shows two membrane models with different sizes of nanopores. The smaller pore diameter was 1.7 nm, and the greater was 2.4 nm. Argon was adopted as the permeating molecule, and the temperature was set at 110 K to maintain liquid-state permeation. The Lennard-Jones (LJ) potential given by Eq. (9) was used to represent the interactions between permeating molecules and molecules that comprise a pore wall. 12

⎡ σij φ (rij ) = 4εij ⎢ ⎜⎛ ⎟⎞ r ⎣ ⎝ ij ⎠

Table 1 Lennard-Jones potential parameters. ε/(kT) [-]

ε/k [K]

σ [nm]

Ar Atoms in membranes Atoms in fluctuating walls

1.13 0.113–4.52 1.13

124.0 0.1–4.0 εAr/k 124.0

0.3418 0.2700 0.3418

liquid state argon. In comparison with the polar molecule of water, since a long-range electrostatic interaction doesn’t work for LJ particles, the viscosity is lower. For example, the viscosity of water at room temperature and pressure is around 1 mPa s, while that of argon under the simulation conditions in this work was around 0.1 mPa s, 1/10 of water viscosity. Therefore, the liquid permeability of argon is expected to be relatively higher than normal water. In addition, no strong electrostatic interaction between permeating argon and the pore surface might prevent a specific liquid state in the vicinity of the pore wall but might cause a more continuous change of the interaction field towards pore center, which possibly affects permeation properties of liquid argon. The energy parameters of the membranes, εmem, determined the interaction between permeating molecules and the pore walls. The εmem and the energy parameters of argon, εAr, established a ratio, εmem/εAr, that ranged from 0.1–4.0 and determined the liquid-permeation characteristics in nanopores. A smaller value for εmem/εAr would virtually represent a system where the permeating molecule-membrane attractive interactions were smaller such as water in a hydrophobic carbon nanotube, while a larger value would correspond to water in a hydrophilic SiO2-ZrO2 ceramic material. A fluctuating-wall MD method was adopted for our liquid state permeation simulations [21]. A schematic image of the non-equilibrium MD simulation cell is shown in Fig. 3. The pressure gradient was kept constant by pressurizing the outside of two plates (fluctuating walls) settled at each end of the unit cell. The upstream side spaces and downstream side spaces were filled with liquid-phase molecules. An equation representing the motion of a fluctuating wall was suggested by the balance of internal and external pressures, namely the motion of the fluctuating wall that was dominated by the balance of applied forces from outside the cell, Apext, where A was an area of the wall, and the force, Fwall, was the summation of forces from fluid molecules. Liquidphase guest molecules were introduced into the pore and permeated along with the movements of the two fluctuating walls. This MD simulation technique enabled us to simulate liquid-phase permeation during a quasi-steady state under given pressures at each end of the cell with a pressure gradient: upstream pressure, pu = 30 MPa; and, downstream pressure, pd = 26 MPa. The temperature was kept constant at 110 K by scaling the kinetic energy of the permeating molecules. At the initial state of simulation (t = 0), argon molecules were settled outside the membrane (upstream side: 3332 molecules, downstream side: 1960 molecules), and the pressure was set as pu = pd = 30 MPa to introduce molecules into a pore by diffusion from both entrances of the membrane (Fig. 4(a)). Within 3 ns, a sufficient number of molecules was

6

σij ⎤ − ⎜⎛ ⎟⎞ ⎥ ⎝ rij ⎠ ⎦

Molecule

(9)

In Eq. (9), rij is the intermolecular distance between molecules i and j, σij and εij are the size and interaction parameters, respectively. The values of the LJ parameters used in this work are summarized in Table 1. In order to focus on the effect of interaction between permeating liquid molecules and the pore surface on liquid permeation properties, argon atom modeled as the Lennard-Jones particle was employed as a permeating molecule with the simple van der Waals interaction. The thermodynamic properties of argon were as follows. The triple point, Tt = 83 K, boiling point, Tb = 87.3 K, critical temperature, Tc = 150.9 K, and critical pressure, pc = 4.90 MPa. In order to simulate liquid state argon and steady state permeation through a nanoscale pore within a finite calculation time, the operating temperature, T was set to be 110 K, a temperature below Tc, and a high enough pressure of pu = 30 MPa was adopted. The liquid density of argon calculated by using the potential parameters listed in Table 1 was 1378 kg/ m3 at 100 K and 15 MPa. It was very close to the actual density of argon, 1349 kg/m3 at 100 K and 10 MPa, which indicated that the employed LJ potential and its potential parameters were suitable for simulation of

Membrane structure Cristobalite silica

d nm

d = 1.7 nm

y

y z x w = 4.3 nm

l = 8.6 nm

x

Fig. 2. Membrane models. 261

d = 2.4 nm

Separation and Purification Technology 220 (2019) 259–267

T. Yoshioka, et al.

Permeation of fluid

Fluctuating wall (pressure control)

Apext

Fwall

pext1

pext2

Upstream side

membrane

Downstream side

Fig. 3. Schematic image of a fluctuating wall MD simulation cell.

Nperm [s−1], could be estimated from the slope of the straight approximated line, and the volumetric permeation rate, Qsim [m3/s], was calculated using Eq. (10).

Qsim. =

Nperm. ρ

(10) −3

In Eq. (10), ρ [m ] indicates the molecular number density of the permeating fluid that was obtained from a bulk MD simulation of the molecule under conditions of T = 110 K and p = 30 MPa. The details of the bulk simulation method are described in Section 4.2.2. 4. Results and discussion 4.1. Permeation simulation Fig. 6 shows the permeation flow rates of argon fluid for NEMD permeation simulations at a steady state plotted against a normalized interaction parameter, εmem/εAr. According to the simple Hagen-Poiseuille equation, the permeation flow rate of a liquid depends on the geometrical factors of a membrane such as pore size and membrane thickness, on the conditions under which it is operated such as temperature and pressure gradient as a driving force, and on the physical properties as viscosity of a permeating fluid, but is independent of an affinity of a permeating fluid with membrane material. The two broken constant lines in Fig. 6 indicate the flow rates of liquid argon, as predicted by Eq. (6), for two different sizes of pores: 1.7 and 2.4 nm. To calculate the theoretical flow rates, the μconst used in Eq. (6) was established via MD simulation of bulk argon molecules. On the contrary, the simulated permeating flow rate showed a clear dependence on the given interaction parameters between a permeating molecule and a membrane molecule, εmem. When the flow rate was decreased, the interaction was increased. With less interaction, where εmem./εAr < 0.5, the observed permeation flow rate exceeded the value predicted by Eq. (6), while an increased interaction caused the flow rates to fall short of

Fig. 4. Schematic image of liquid permeation simulation.

4000 upstream

N [-]

3000 2000

downstream

1000 0

in-pore

0

10

20

Time [ns]

10-17

Fig. 5. Example of the number of fluid molecules as a function of time (d = 2.4 nm, εmem./εAr = 0.1).

d = 2.4 nm Q [m3/s]

transported into the pore to create a liquid-like high density phase in the nanopore, and we confirmed the number of molecules in each region upstream-outside/downstream-outside, and inside the pore. After the system attained a state of equilibrium under a uniform pressure, the downstream pressure was changed to pd = 26 MPa, and a liquid-state permeation simulation was begun under a pressure gradient, Δp = 4 MPa (Fig. 4(b)). The time step for integration was Δt = 10 fs and the interval of velocity scaling was 100 steps. Fig. 5 shows the time course for the number of molecules in each region upstream-outside/ downstream-outside, and inside the pore for permeation through the d = 2.4 nm membrane model when εmem/εAr = 0.1. When a pressure gradient was applied after 3 ns, the number of molecules outside the membrane either increased or decreased at an identical rate, which meant the system was under a steady-state flow. The permeation rate,

10-18

d = 1.7 nm

10-19

0

1

2

3

4

εmem./εAr[-] Fig. 6. Relationship between permeability and reduced interaction (Broken: Hagen-Poiseuille, Plot: Simulation, Solid curve: Model). 262

Separation and Purification Technology 220 (2019) 259–267

T. Yoshioka, et al.

3

pore wall, but for a larger interaction (εmem./εAr = 3.0), the slope became smaller in the vicinity of the pore wall. In both cases, the shape of uz(r) obviously did not coincide with an ideal parabolic curve. These results suggest that an ordinary Hagen-Poiseuille flow model that uses constant viscosity may fail to describe the liquid-state flow in a nanopore. There may be limitations to the application of a bulk homogeneous fluid model to a nano-confined liquid.

εmem./εAr = 0.1

2

uz(r) [m/s]

εmem./εAr = 1.0 4.2. Calculation by the proposed model As a first step in modifying the normal Hagen-Poiseuille flow model, the inhomogeneous viscosity of a permeating fluid confined in a nanopore was introduced in the present work. We assumed that the fluid viscosity could change in the vicinity of a pore wall due to the pore wall potential acting on the fluid molecules as a function of location in the nanopore, r, or as a function of the distance from the pore wall. This assumption seems plausible because, as a physical property, viscosity largely depends on primitive intermolecular interactions. In a nanopore, a considerable portion of a confined fluid must be under the influence of interactions with the pore wall, which may be an additional interaction factor for permeating molecules that could cause meaningful change in the local viscosity of nano-confined liquid. If the viscosity of a nano-confined liquid can be estimated as a function of r, a velocity distribution in the direction of flow, uz(r), can be given by Eq. (8) based on the basic concept of Newton’s law of viscosity (Eq. (7)). The viscosity distribution in a nanopore, μ(r), was estimated using the following procedures. In general, the viscosity of a fluid can be calculated from a self-diffusion coefficient via the Stokes-Einstein equation [22]. However, it is difficult to directly calculate the local diffusivity of permeating molecules with sufficient accuracy in a region with thickness, Δr, on almost a sub-nanometer scale in a nanopore. Therefore, the local viscosity, μ(r), was estimated by combining the local potential in a nanopore, E(r), that included both the inter-permeating molecular interactions and the external potential acting on the permeating molecules from the pore wall, along with the viscosity of Lennard-Jones fluid, μ(E), as a function of potential, E(ε). For model LJ molecules with different inter-molecular energy parameters, ε, in Eq. (9), total potential energy, E, and viscosity, μ, were calculated from MD simulations in bulk.

1

εmem./εAr = 3.0 0

-1

0

1

r-direction [nm]

Fig. 7. Observed velocity profiles by MD simulation (d = 2.4 nm); the data have been simply reflected to negative values of r for ease of visualization (the same in Figs. 8, 13, and 14).

the expectations for both nanopore models. Decreases in the simulated liquid flow rate caused by an increase in the affinity of a permeating molecule for a membrane was qualitatively similar to experimental observations of alcohol permeation through the nanopores of SiO2-ZrO2 nanofiltration membranes [1,2]. If the permeation state deviated from a normal Hagen-Poiseuille flow under an assumption of homogeneous viscosity, the velocity distribution in the pore could be shifted or distorted from an ideal parabolic form. In this MD simulation, the local velocity in the direction of flow, uz(r), was easily calculated by averaging the velocity in z, uz(r, t, i), for a sufficient number of molecules located between r − Δr/2 and r + Δr/2 for a certain simulation period, where i denotes the i-th molecule. Fig. 7 shows the velocity distributions in the direction of flow as a function of the location in the pore, r, for three different interactions in the d = 2.4 nm pore. The data have been simply reflected (symmetrically) to negative values of r for ease of visualization (the same in Figs. 8, 13, 14). The broken curve indicates a parabolic profile calculated using Eq. (5). For a smaller interaction (εmem./εAr = 0.1), the slope of the velocity distribution suddenly increased with approaches to the

4.2.1. Potential profiles, E(r) Fig. 8 shows the potential distributions from differing interactions of liquid state argon molecules as they fill nanopores. The potential

0

0

-10

E(r) [kJ/mol]

E(r) [kJ/mol]

εmem./εAr = 0.1 εmem./εAr = 1.0

-20

εmem./εAr = 3.0

-10

-20

(a) -30

-1

0

(b) -30

1

r-direction [nm]

-1

0

r-direction [nm]

Fig. 8. Potential profile in nano-pores: (a) d = 2.4 nm, (b) d = 1.7 nm. 263

1

Separation and Purification Technology 220 (2019) 259–267

T. Yoshioka, et al.

60

minimum in the vicinity of the pore wall was not seen in weaker permeating molecules with a membrane pore surface interaction of εmem/ εAr = 0.1, and the potential energy increased exponentially from the center of the pore to the pore wall. On the contrary, for pores showing εmem/εAr > 1.0, the potential minimum was obvious, and the depth increased as increases in the interaction between permeating molecules and the membrane pore surface. When εmem/εAr = 1.0, where interaction corresponded to Ar molecular interaction, the potential minimum was observed due to higher molecular density on the pore surface than that of liquid Ar. The Lennard-Jones dispersion interaction was known as a short-range force, and 3.5σAr was adopted as a cut-off distance [23] in this work. Therefore, at the center of a pore with a diameter of 2.4 nm, the potential energy on the fluid was almost constant and was independent of the interaction between the permeating molecule and the pore surface (εmem/εAr), which indicated that the permeating molecules had no energetic effects caused by the pore surface (Fig. 8(a)). However, for a pore diameter of 1.7 nm, the interaction from molecules constituting the pore wall reached the center of the pore, and the observed potential energy depended on the value of εmem/εAr (Fig. 8(b)). The potential distributions in nanopores, E(r), were estimated as a function of the distance from the pore center, r, by numerical adjustment of the simulated data to an empirical polynomial expression.

MSD [nm2]

50

kB T 1 3πα D

MSD =

t →∞

N

∑ 〈 (ri (t ) − ri (0))2〉 i=1

1 1 〈 (r (t ) − r (0))2〉 = lim 〈MSD〉 t →∞ 6t 6t

ln μ = A +

(13) (14)

Ea λE =A+ RT RT

(15)

The Andrade equation was theoretically derived by Andrade in 1934 [26]. This equation is widely used since it describes the experimentally observed temperature dependency of liquid viscosity data. The value of A in Eq. (15) is a constant that depends on the liquid material, and Ea is the apparent activation energy for fluidization. This apparent activation energy can be reasonably assumed to correlate with the potential energy. Therefore, we have introduced a proportional constant, λ, to the Andrade equation and described liquid viscosity, μ, as a function of potential energy, E, at a constant temperature, T. In Fig. 12, the estimated viscosity in this work (110 K) was compared with experimental (100–145 K) and other MD (100–145 K) data [25]. The viscosity of liquid argon (110 K, 30 MPa) based on the StokesEinstein equation was close to the experimental values at densities around 1300–1400 kg/m3. In addition, the estimated values of viscosity in this work were also showed agreement with the direct calculation results from microscopic stress tensor data obtained by MD simulation

RDF [-] 1

10

3

In order to evaluate the viscosity dependency of the Ar-like fluid caused by the potential energy of the surrounding molecules, several virtual Ar fluids with different levels of interaction or potential energy were simulated by varying the value of the ε parameter with that of the fixed σ parameter shown in Table 1. Fig. 11(a) shows the calculated viscosity plotted against the potential energy of the fluid corresponding to the varied ε parameters. A lower potential energy meant a greater attractive interaction between molecules, which corresponded to a larger ε parameter. Therefore, the observed viscosity of the virtual fluid drastically increased with an increase in the inter-molecular attractive interaction. Fig. 11(b) shows the relationship between viscosity and density of liquid argon. As the figure demonstrates, there was a strong positive correlation between the viscosity and density. Since larger inter-atomic interaction caused increase in density and higher density would lead to frequent exchange of momentum among argon atoms to increase viscosity, this seems to be a reasonable result. The exponential increase in viscosity was fitted by the Andrade equation given by Eq. (15).

2

5

1 N

D = lim

(12)

α

2

and n(r) is the number of molecules in the region between r and r + Δr. The mean squared displacement (MSD) is given by Eq. (13) where N is the number of molecules and ri(t) − ri(0) is the displacement of i molecules for a certain period, t. The self-diffusivity was calculated from the slope of the MSD curve in Fig. 10 using Einstein’s equation (Eq. (14)) [25].

3

0

1

Fig. 10. Mean squared displacement of liquid Ar.

For Eq. (12), ρ is the average density of the molecules in the system

0

0

Time [ns]

(11)

1 n (r ) ρ 4πr 2Δr

20

0

In Eq. (11), k is the Boltzmann constant, α is the effective molecular size in the fluid, and D is the self-diffusivity. In order to obtain α and D of the simulated fluid, MD simulations of bulk fluid were conducted for a constant number of molecules (N = 864). An arbitrary initial density was set and NPT (N, p, T: constant) MD simulation was carried out under the same conditions as for permeation simulation (T = 110 K, p = 30 MPa) to stabilize the system. After reaching a steady state with an almost constant density, the simulation conditions were switched to NVT (N, V, T: constant) MD simulation mode, and the location of the molecules was recorded along with time in order to calculate the radialdistribution function (RDF) and mean squared displacement (MSD). The radial distribution function (RDF) is defined by Eq. (12) [24], and the effective molecular size, α was estimated from the location of the first RDF peak, as shown in Fig. 9. As shown in Fig. 9, the value of α was 0.371 nm in this case, which was close to the size parameter of argon in Table 1, 0.3418 nm.

PDF(r ) =

30

10

4.2.2. Fluid viscosity, μ(E) The viscosity of fluid was estimated using the Stokes-Einstein equation (Eq. (11)).

μ=

40

15

r [Å] Fig. 9. Radial distribution function of liquid Ar. 264

Separation and Purification Technology 220 (2019) 259–267

T. Yoshioka, et al.

(a)

(b)

10-4

μ [Pa s]

10-4

10-5

-10

-5

0

10-5 0

500

E [kJ/mol]

1500

1000

Density, ρ

[kg/m3]

Fig. 11. Relationship between (a) potential energy and viscosity and (b) density and viscosity (T = 110 K, p = 30 MPa).

10-3

Exp. (100-145 K) MD (100-145 K) This work (110 K)

μ(r) [Pa s]

μ [Pa s]

2

1

εmem./εAr = 3.0 εmem./εAr = 1.0

10-4

εmem./εAr = 0.1 10-5

(a)

3

εmem./εAr = 0.1 0 500

1000

Density, ρ [kg/m3]

1500

(b)

2

εmem./εAr = 1.0 uz(r) [m/s]

Fig. 12. Comparison of estimated viscosity in this work (110 K) with experimental (100–145 K) and other MD (100–145 K) data [27].

of LJ argon particles. These results suggested that practical values of liquid argon could be estimated by our method. However, Eq. (11) used here is a relationship between diffusivity coefficient and viscosity for the case of no slip at the surface of the liquid particles modeled as hard sphere. It should be noted that this calculation method doesn’t always give an accurate viscosity, since Eq. (11) is an approximated model for evaluation of viscosity from diffusivity coefficient.

1

εmem./εAr = 3.0 0

4.2.3. Viscosity profiles μ(r) Fig. 12(a) shows the calculated viscosity distribution, μ(r), in a pore of the 2.4 nm model, which was obtained by combining the potential distribution, E(r), and the fluid viscosity, μ(E), as a function of intermolecular potential, E. The viscosity distribution depended not only on the inter Ar molecular interaction but also Ar-membrane pore surface interaction. A greater level of interaction between the molecules of the fluid and the pore wall tended to deepen the overall potential depth in the vicinity of the wall, which increased the viscosity of the fluid. Smaller interactions between fluid molecules and the pore wall, on the other hand, lowered the viscosity of the fluid around the pore wall due to either shallower or rising potential. Numerical computation suggested that the fluid viscosity in the vicinity of the pore wall could be several times higher or lower than the viscosity, which was very close to the normal bulk fluid viscosity at the center of the pore.

-1

0

r-direction [nm]

1

Fig. 13. Calculated viscosity profiles and velocity profiles (Plot: Simulation, Solid line: Calculated).

4.2.4. Velocity profiles uz(r) and flux The velocity distributions in the pores of the 2.4 nm model were calculated by integrating the results from Eq. (8) with the predicted viscosity distributions of the fluid in the pores, μ(r), as shown in Fig. 12(b) as solid curves. The calculated velocity distributions deviated from the parabolic curve derived from the Hagen-Poiseuille flow model, and showed good agreement with the MD simulation results. The fluid viscosity was lowered when the attractive interaction was small due to the high level of potential energy in the vicinity of the pore wall. Therefore, a large velocity gradient occurred there, which led to an increase in the overall volumetric flow rate of the fluid. On the 265

Separation and Purification Technology 220 (2019) 259–267

T. Yoshioka, et al.

40

the volumetric liquid flow rate of nanopores in this study deviated from that posited by the Hagen-Poiseuille theorem, as reported in other experimental studies. In this study, the flow rate was higher for smaller attractive interactions between permeating molecules and the pore wall, while it became smaller with increases in the interactions. We decided that the inhomogeneity of the liquid viscosity in a nano-pore should be taken into account to explain these phenomena, and numerically predicted a velocity profile in nano-pores by solving the Navier-Stokes equation with viscosity distribution as a function of the location in the radial direction in a pore. The viscosity distribution in a nano-pore was estimated quantitatively from the potential distribution in the pore. The calculated velocity distribution formed a slightly different shape from that of a normal parabolic curve, as predicted by the Hagen-Poiseuille theorem, but showed good agreement with the observed velocity distribution in the MD simulation results. The volumetric flow rate as an integrated value for velocity distribution also well reproduced the simulated flow rate. Deviation of the flow rate from the Hagen-Poiseuille theorem can be explained by using a viscosity distribution that is based on the potential distribution in a nano-pore.

εmem./εAr = 3.0

30

ρN [/nm3]

εmem./εAr = 1.0 20

εmem./εAr = 0.1

10

0

-1

0

r-direction [nm]

References

1

[1] T. Tsuru, S. Wada, S. Izumi, M. Asaeda, Silica-zirconia membranes for nanofiltration, J. Membr. Sci. 149 (1998) 127–135, https://doi.org/10.1016/S0376-7388(98) 00163-X. [2] T. Tsuru, T. Sudou, S. Kawahara, T. Yoshioka, M. Asaeda, Permeation of liquids through inorganic nanofiltration membranes, J. Colloid Interface Sci. 228 (2000) 292–296, https://doi.org/10.1006/jcis.2000.6955. [3] F. Du, L. Qu, Z. Xia, L. Feng, L. Dai, Membranes of vertically aligned superlong carbon nanotubes, Langmuir 27 (2011) 8437–8443, https://doi.org/10.1021/ la200995r. [4] M. Majumder, A. Stinchcomb, B.J. Hinds, Towards mimicking natural protein channels with aligned carbon nanotube membranes for active drug delivery, Life Sci. 86 (2010) 563–568, https://doi.org/10.1016/j.lfs.2009.04.006. [5] D. Mattia, F. Calabro, Explaining high flow rate of water in carbon nanotubes via solid-liquid molecular interactions, Microfluid. Nanofluid. 13 (2012) 125–130, https://doi.org/10.1007/s10404-012-0949-z. [6] J.A. Thomas, A.J.H. McGaughey, Reassessing fast water transport through carbon nanotubes, Nano Lett. 8 (2008) 2788–2793, https://doi.org/10.1021/nl8013617. [7] J.A. Thomas, A.J.H. McGaughey, Water flow in carbon nanotubes: transition to subcontinuum transport, Phys. Rev. Lett. 102 (2009) 184502, https://doi.org/10. 1103/PhysRevLett. 102.184502. [8] J.A. Thomas, A.J.H. McGaughey, O.K. Arnebeck, Pressure-driven water flow through carbon nanotubes: insights from molecular dynamics simulation, Inter. J. Therm. Sci. 49 (2010) 281–289, https://doi.org/10.1016/j.ijthermalsci.2009.07. 008. [9] X. Qin, Q. Yuan, Y. Zhao, S. Xie, Z. Liu, Measurement of the rate of water translocation through carbon nanotubes, Nano Lett. 11 (2011) 2173–2177, https://doi. org/10.1021/nl200843g. [10] S.K. Kannam, B.D. Todd, J.S. Hansen, P.J. Daivis, How fast does water flow in carbon nanotubes? J. Chem. Phys 138 (2013) 094701, https://doi.org/10.1063/1. 4793396. [11] C. Neto, D. Evans, E. Bonaccurso, H. Butt, V. Craig, Boundary slip in Newtonian liquids: a review of experimental studies, Rep. Progr. Phys. 68 (2005) 2859, https:// doi.org/10.1088/0034-4885/68/12/R05. [12] T.G. Myers, Why are slip lengths so large in carbon nanotubes? Microfluid Nanofluid 10 (2011) 1141–1145, https://doi.org/10.1007/s10404-010-0752-7. [13] J. Thomas, A. McGaughey, Density, distribution, and orientation of water molecules inside and outside carbon nanotubes, J. Chem. Phys. 128 (2008) 084715, https:// doi.org/10.1063/1.2837297. [14] A. Botan, B. Rotenberg, V. Marry, P. Turq, B. Noetinger, Hydrodynamics in clay nanopores, J. Phys. Chem. C 115 (2011) 16109–16115, https://doi.org/10.1021/ jp204772c. [15] D.J. Bonthuis, R.R. Netz, Unraveling the combined effects of dielectric and viscosity profiles on surface capacitance electro-osmotic mobility, and electric surface conductivity, Langmuir 28 (2012) 16049–16059, https://doi.org/10.1021/la3020089. [16] D.J. Bonthuis, R.R. Netz, Beyond the continuum: how molecular solvent structure affects electrostatics and hydrodynamics at solid-electrolyte interfaces, J. Phys. Chem. B 117 (2013) 11397–11413, https://doi.org/10.1021/jp402482q. [17] S.P. Sutera, R. Skalak, The history of Poiseuille’s law, Ann. Rev. Fluid Mech. 25 (1993) 1–20, https://doi.org/10.1146/annurev.fl.25.010193.000245. [18] M. Cheryan, Ultrafiltration and microfiltration handbook second edition Chapter 4, CRC Press, New York pp. 116–120, ISBN: 1-56676-598-6; 1998. [19] D.O. Young, H.C. Chiu, S. Kim, K. Voitchovsky, E. Riedo, The interplay between apparent viscosity and wettability in nanoconfined water, Nat. Commun. 4 (2013) 2482, https://doi.org/10.1038/ncomms3482. [20] T. Yoshioka, M. Asaeda, T. Tsuru, A molecular dynamics simulation of pressuredriven gas permeation in a micropore potential field on silica membranes, J. Membr. Sci. 293 (2007) 81–93, https://doi.org/10.1016/j.memsci.2007.01.039.

Fig. 14. Number density profiles of molecules in a nano-pore.

contrary, when the interaction with the wall was large, a deep potential minimum in the vicinity of the pore wall resulted in a high level of viscosity. Since high viscosity prevents the formation of a high velocity gradient, the overall flow rate is lowered. According to the literature, the possibility for inhomogeneity in fluid viscosity is much different in the vicinity of a nano-pore wall than it is in the center of the pore [5,12,13]. That theory of fluid viscosity was supported by the density distributions in the radial direction of the pores in this study, as shown in Fig. 13. Calculating the velocity distribution by considering viscosity as a function of r capably predicted the MD-simulated velocity distribution, which means the proposed model accurately described the liquid-state, fluid-flow phenomena in a nano-pore. Finally, overall volumetric flow rates were calculated by integrating the predicted velocity distribution, uz(r), as given by Eq. (16), for two nano-pore models with variable fluid molecule-wall interaction parameters, εmem/εAr.

Qmodel =

∫ 2πruz (r ) dr

(16)

The calculated flow rate showed both qualitative and quantitative agreement with that obtained by the MD simulations for different εmem/ εAr parameters, as shown in Fig. 6. The deviation of the flow rate from the Hagen-Poiseuille theorem could be explained by basing the local viscosity of fluid on the potential distribution within a nano-pore. The overall trend of the volumetric flow rates and velocity distributions could be well explained by considering the local viscosity as a function of location in the pore, r. On the other hand, the finer oscillatory structure seen in the velocity profiles obtained from the simulations was not perfectly recovered using the velocity profiles calculated by the proposed model. This might suggest the limitation of this continuous fluid approximation model for fluid behaviors in nanopores. In order to estimate velocity profile and flow rate more precisely, consideration of some discrete model with, for example, molecular excluded volume might be required. 5. Conclusions In this work, two cylindrical pore models with pore diameters of 1.7 and 2.4 nm, respectively, were prepared, and the liquid-phase permeation phenomena in nano-pores were studied using a non-equilibrium MD (NEMD) simulation method. The MD simulation results for 266

Separation and Purification Technology 220 (2019) 259–267

T. Yoshioka, et al. [21] H. Takaba, Y. Onumata, S. Nakao, Molecular simulation of pressure-driven fluid flow in nanoporous membranes, J. Chem. Phys. 127 (2007) 054703, https://doi. org/10.1063/1.2749236. [22] I. Avramov, Relationship between diffusion, self-diffusion and viscosity, J. NonCrystal. Solids 355 (2009) 745–747, https://doi.org/10.1016/j.jnoncrysol.2009.02. 009. [23] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford Science Publications, Oxford pp. 7–12, ISBN: 0-19-855645-4; 1989. [24] C.Y. Gao, D.M. Nicholson, D.J. Keffer, B.J. Edwards, A multiscale modeling demonstration based on the pair correlation function, J. Non-Newton. Fluid Mech.

152 (2008) 140–147, https://doi.org/10.1016/j.jnnfm.2007.05.003. [25] D. Coglitore, S.P. Edwardson, P. Macko, E.A. Patterson, M. Whelan, Transition from fractional to classical Stokes-Einstein behaviour in simple fluids, R. Soc. Open Sci. 4 (2017) 170507, https://doi.org/10.1098/rsos.170507. [26] E.N.C. Andrade, A theory of the viscosity of liquids-part I, Phil. Mag. 17 (1934) 497–511, https://doi.org/10.1080/14786443409462409. [27] G.A. Fernández, J. Vrabec, H. Hasse, A molecular simulation study of shear and bulk viscosity and thermal conductivity of simple real fluids, Fluid Phase Equilibria 221 (2004) 157–163, https://doi.org/10.1016/j.fluid.2004.05.011.

267