Nanofluidic Manifestations of Structure in Liquids: a Toy Model | The

Jun 18, 2019 - Interaction with a solid boundary can alter the local order in the liquid near the boundary- an effect that is of special importance in...
0 downloads 0 Views 6MB Size
Article Cite This: J. Phys. Chem. C 2019, 123, 16787−16795

pubs.acs.org/JPCC

Nanofluidic Manifestations of Structure in Liquids: A Toy Model A. Z. Patashinski,† M. A. Ratner,† R. Orlik,‡ and A. C. Mitus*,§ †

Northwestern University, Department of Chemistry, 2145 Sheridan Road, Evanston, Illinois 60208, United States Orlik Software, ul. Przestrzenna 35/1, 50-534 Wroclaw, Poland § Wroclaw University of Science and Technology, Department of Theoretical Physics, Wybrzeze Wyspianskiego 27, 50-370 Wroclaw, Poland

Downloaded via NOTTINGHAM TRENT UNIV on August 27, 2019 at 09:27:18 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.



ABSTRACT: Interaction with a solid boundary can alter the local order in the liquid near the boundaryan effect that is of special importance in geometries where these layers of altered structure occupy a significant part of the liquid. Currently, the local order is unambiguously confirmed and visualized only in two-dimensional (2D) liquids. A 2D liquid may play the role of a toy model for more complex liquids. We use molecular dynamic simulations to study instantaneous and time-averaged particle positions and velocities, quantitative characteristics of the local order, and drag forces in a 2D Lennard-Jones liquid occupying a narrow channel with moving walls. In layers of the liquid adjacent to the walls, long-time-averaged velocity significantly deviates from the predictions of conventional hydrodynamics. In the same layers, the time-averaged local order characteristics and the particle density show significant changes, with the amplitudes of these changes depending on the thermodynamic state of the liquid, the wall roughness, and commensurability with the wall. rheological22 local-structure-based effects. One suggests8,18 that these effects are determined by the commensurability of the local structures of the liquid with that of the wall. Numerous studies used computer simulations as well as theoretical models to provide a connection of the atomic-scale picture of the liquid−solid interface with the molecular theory of transport in inhomogeneous systems. The parameters of common interest are the density and velocity profiles and their relation to various types of boundary conditions, slips, and structure of the walls. An exhaustive discussion of an atomic interfacial friction and its relation to slip in nanoscopic Couette flows can be found in ref 23. Those topics were studied using nonequilibrium molecular dynamics simulations,18,23−25 the lattice version of Boltzmann equations,26 the kinetic Enskoglike model,27 or the generalized hydrodynamic model.28 The setups were similar to that used in our study. Those models offer, apart from a large body of detailed results, interesting physical pictures of nanoflows, based in part on the concept of global order in layers of fluid close to the walls, expressed in terms of static structure factors18,23 or spatial probability density.18 The critical difference between our study and those found in the literature is our choice of the locally ordered mosaic state of the liquid (see section 4 for details). It was found in earlier studies (see below) that 2D liquids have a recognizable local structure only in a narrow band of states in the thermodynamic (density−temperature) plane; this band is

1. NANOFLUIDICS AND LOCAL ORDER IN LIQUIDS: A PHYSICAL PICTURE In solid materials, channels having the width of a few nanometers (nanochannels) naturally appear as cracks and gaps. In recent decades, nanochannels are engineered for special nanofluidic applications.1 Experiments2−6 and computer simulations7,8 show that the predictions of the Navier− Stokes9 (NS) hydrodynamics give a valid approximation for channels of a width larger than (typically) a few tens of molecular diameters but fail for narrower geometries. One possible mechanism leading to these violations is the coupling of hydrodynamics and electrochemistry. This mechanism is characterized by length scales (typically) in the tens of nanometer range: the Debye, Dukhin, or slip lengths. For nanochannels, some specific nanofluidic effects (e.g., flow enhancement 10−12 ) are expected that violate the NS approximation. A critical discussion of these concepts can be found in the review paper.1 Here, we propose and prove a different mechanism of violations of the NS picture in nanometer-scale geometries: a change of the local liquid structure near the walls due to wall− liquid interaction. It is known from experiment and computer simulations that both the characteristics of the boundaries and those of the liquid become important in nanoscopic geometries: the characteristics of the boundary influence the flow of water in confinement;4,13 slippage on surfaces depends crucially on their wettability.2,6,8,14−16 We assume that interaction with the walls can change the local order in the liquid near the walls of a channel17,18 and thus increase the tangential pressures19,20 (see also ref 21), promoting non© 2019 American Chemical Society

Received: April 19, 2019 Revised: June 17, 2019 Published: June 18, 2019 16787

DOI: 10.1021/acs.jpcc.9b03697 J. Phys. Chem. C 2019, 123, 16787−16795

Article

The Journal of Physical Chemistry C

refs 33 and 34) due to the competition between close packed local structures: fcc, hcp, and icosahedron. The objective of the study is to prove the suggested physical picture using the 2D LJ liquid as a mostly schematic (toy) model for a locally ordered liquid. To this end, we present a quantitative analysis of the impact of the inherent and boundary-modified local order in a 2D LJ liquid on the nanofluidic 2D Couette flow. In the general case, ions and/or impurities may be present in a locally ordered liquid. Then, coupling of electrochemical and structure-related effects is expected to result in a more complex physical picture. For a first study, we simplify the task by choosing a liquid where local order-related effects are the only generic mechanisms that can violate the NS picture.

approximately parallel to the liquidus line and includes this line as well as the states close to it but at lower densities. The choice of temperature and average density of the liquid serves the main goal of our study: the manifestation of the local structure of the liquid. Another important new element is the choice of the value of the induced velocity gradients. In our study, it is orders of magnitude smaller than in those found in the literature. The small value of the gradient excluded disruption of the local structure near the wall, formation of layers substantially differing in physical properties, slip phenomena, and other strong nonlinear effects. In the linear regime, the effects are not as strong; to study these effects, it was necessary to control the relaxation to a steady flow state and, once this steady state was reached, to average the velocity and other characteristics over times sufficient to suppress fluctuations in this steady state. The analysis of local structural aspects of the nanoflow offers complementary to the literature studies information which, in turn, can support the aforementioned physical pictures. The physical picture of local-structure mediated nonNewtonian flow is as follows. The solid wall of the nanochannel changes the parameters of the liquid in a layer of the liquid adjacent to it. This includes a nonstructural effectchange of the densityand a structural onechange of the local (short-range) order of the liquid. The origin of the latter is the fact that the wall plays the role of an external field, with some (crystal) symmetry, which acts on and modifies liquid’s local structure, characterized by its own (crystal) symmetry. The strength of this effect decays as the distance from the walls increases. The range of the direct interaction between atoms is about the atom size, but correlations in the liquid can transfer the influence of the wall and change the local order in a layer of the liquid that has a larger width. This crossover spatial scale depends on the quality (degree of crystallinity) of local structure of the liquid, on the compatibility between wall-induced and inherent liquid’s symmetries, and on the state (defined by temperature and bulk density) of the liquid. Namely, a compatible wall promotes crystalline order and clusterization at the wall, while a noncompatible wall suppresses those effects. In a nanofluidic system, boundary layers of altered structure may occupy a substantial part of the channel, so the transport properties of the entire channel become substantially changed. The presented physical picture is a consequence of the statistical mechanics of locally ordered liquids29 and is expected to be generic for other locally structured liquids. One notes that the local structure of liquids is still a rather intuitive and nonquantitative concept, and the problem of liquid structure is still open. The suggestion of local order in some three-dimensional (3D) liquids is supported by lightscattering data30 and observations of liquid−liquid phase transitions,31,32 and also computer simulations.33−37 However, quantification of this order is still a challenge. Thus, it is desirable to find a simple but at least partially realistic case where the local order in the liquid can be unambiguously defined and quantified. As a possible candidate for this task, we study here the 2D Lennard-Jones (LJ) liquid. The local (triangular) order in this liquid is unique and well quantified.38−42 Another advantage of a 2D system is that a 2D configuration can be easily visualized and analyzed. In contrast, analysis of the local structures in even simple 3D liquids encounters serious methodological problems (see, e.g.,

2. THE SYSTEM: 2D COUETTE FLOW IN A NANOCHANNEL Constant temperature molecular dynamics43,44 was used to simulate a 2D system of Lennard-Jones atoms (gray circles in Figure 1 (top)) occupying the space between two parallel walls

Figure 1. Schematics of the system (top). Phase diagram of 2D LJ liquid close to T = 2 (bottom). L and S stand, respectively, for liquidus (red) and solidus (blue) lines.

made of wall particles (black circles in Figure 1 (top)). Wall particles excerpted LJ forces on the liquid atoms, but the relative positions of wall particles have been frozen in approximately crystalline relative positions (section 5.1). In our simulations, the wall is modeled on the atomic scale, and that eliminates the necessity to impose special boundary conditions. Our choice of the boundary corresponds to the wall atoms having much larger masses than the atoms of the liquid. Since the densities of the wall atoms and of the liquid were approximately equal (section 5.1), one expects18 that noslip boundary conditions are valid. For shear velocity gradients larger than those used in this study, the relation between friction and slip becomes more complex.23 The walls moved along the x-axis as solid bodies having velocities of +v for the upper wall and −v for the bottom one. A wall consisted of two layers of particlesfor densities of the wall particles used in the study, an increase of the number of layers in a wall did not 16788

DOI: 10.1021/acs.jpcc.9b03697 J. Phys. Chem. C 2019, 123, 16787−16795

Article

The Journal of Physical Chemistry C

(∇v)* ≡ v /Ly ·τ /σ 2 = 10−4 . For a system of 50 × 25 atoms, simulating under steady conditions for about 2.5 × 107 h brought the system close to a (statistical) steady state, for which velocity profiles vx(y) and other properties were analyzed. The instantaneous velocity vx(x, y, t) averaged over small parts of the liquid had large fluctuations. To obtain the profiles vx(y), we partitioned the liquid into narrow (of the width ∝ 2.5σ) strips separated by horizontal lines y = const. The velocity profile vx(y) was calculated by first averaging over the strip centered by the corresponding y-coordinatethis yields the still significantly fluctuating velocity vx(y, t)and then averaging over a long enough (see below) time θ: vx(y) = ⟨vx(y, t)⟩θ. As a practical definition, we identify a state as steady when the profiles vx(y) are steady within the accuracy of 2%, and also reproducible in repeating simulations. The smoothing time θ necessary for the steady-state fluctuations to fit into the 2% limit depended on temperature and density; for T = 2 and ρ = 0.92, θ ≥ 5 × 106 h. We assume that reproducibility of the smoothed profiles indicates that the time θ is sufficient for the system to sample a representative ensemble of configurations and local structures. Then, on these asymptotically large time scales, the liquid can be treated as a continuum medium. Navier−Stokes hydrodynamics9 predicts that for slow large scale motions vx(y) has a linear dependence on y:

lead to significant changes of the results. The LJ potential has the form ULJ(r ) = 4ϵ[(σ /r )12 − (σ /r )6 ]; large values of the microscopic wetting parameter αW ≫ 1 indicated strong wetting conditions for the system at the temperature and density of simulations.21 The relative motion of the walls induced an average flat nanoscopic shear (Couette-type) flow9 in the liquid. We use the vibration period τ = ϵ/(mσ 2)1/2 for a particle with mass m in the harmonic part of the LJ potential as the time unit; the integration step was h = τ/500. Systems of various sizes were studied (up to 4 × 104 atoms), but to obtain representative results, simulating much smaller systems of a few thousand atoms appeared sufficient. The initial configuration for liquid particles was that of an ideal (triangular) solid composed of Ly horizontal layers each of Lx atoms, with the distance between the layers equal to approximately 21/6σ 3 /2 ≈ σ ; melting of this crystalline configuration took about 102τ. Reaching a stationary flow (defined below) in the liquid took a time about an order of magnitude larger. Below, we report the results for Lx = 50, Ly = 25 and comment on the case Ly = 40. The shear flow induced by the moving walls releases heat in the liquid, so some form of a thermostat is necessary to keep the temperature in a steady state constant. To produce the NVT ensemble, we used the Nosé−Hoover thermostat;43,44 the applicability of this approach to liquids under shear-flow conditions was demonstrated in a number of studies (see, e.g., ref 45). This choice of thermostatting, from many thermostatting strategies,23,43,44 was motivated by its robustness noneqilibrium results agree with other methods of thermostatting studied in ref 45. To ensure that the system is in a steady state, the short- and long-time-averaged temperature profiles T(y, t) have been continuously controlled together with the velocity and density profiles. In the steady states studied, the long-time-averaged temperature profiles were constant in the whole system within the 2% fluctuation range. The reduced temperature of the liquid T* ≡ kBT/ϵ = 2 was substantially larger than the critical temperature Tc* ≈ 0.56, while the reduced number density of the liquid ρ* ≡ ρσ2 ranged from 0.8 to 0.96 (in what follows, we use the notations T and ρ for the reduced temperature and number density, respectively). The equilibrium phase diagram of the 2D Lennard-Jones system close to T = 2 is shown in Figure 1 (bottom). The nature of the melting transition in 2D systems remains unclear and includes the possibility of a discontinuous phase transition as well as two continuous phase transitions with an intermediate hexatic phase (see, e.g., refs 38, 39, and 46 and the literature cited therein), but in both cases, there exist two densities, ρL < ρS (ρL ≈ 0.94, ρS ≈ 0.98 for T = 2), which separate liquid and crystal phases. In the former case, they are referred to as liquidus and solidus lines, which constitute the boundaries of the coexistence region, while, in the latter case, they determine the boundaries of the hexatic case. For historical reasons, we use the former notions. The current study is oriented onto local structural aspects of a liquid in a small system, and thus, the aforementioned ambiguity is of no importance.

vx(y) = 2vy/Ly( −Ly /2 ≤ y ≤ Ly /2)

(1)

Figure 2 presents the asymptotic steady-state profiles vx(y) for ρ = 0.88, 0.92, and 0.94. For densities ρ = 0.88 and lower,

Figure 2. Velocity profiles vx(y) for a few densities shown in the legend. Solid red line: hydrodynamic prediction. Dashed black line: hyperbolic tangent fit. Black solid lines: analysis of data for ρ = 0.94; see text.

the profiles appeared close to the hydrodynamic predictions; as the density of the liquid increased, deviations from these predictions became apparent, and layers of decreased mobility appeared near the walls. Increase of the density led to widening of these layers. In a certain sense, these layers acted as extensions of the solid boundaries. In this case, the slopes of the profiles in the central part became higher than predicted by Navier−Stokes equations. The characteristic width of the boundary layer, where the flow is hindered by the walls, was estimated by fitting the profiles with a hyperbolic tangent function: vx(y) = −tanh(α(y − y0)) and then using the linear approximations in the middle part of the flow and at the walls and finding their intersection point. This procedure is illustrated for the case ρ = 0.94 in Figure 2. The total width (2) (1) (2) of both boundary layers Δv = Δ(1) v + Δv , where Δv and Δv represent, correspondingly, the widths of the upper and lower boundary layers, reads Δv = 15σ, which constitutes around

3. FLOW AT CHANNEL WALLS 3.1. Velocity Profiles: Deviations from the Navier− Stokes Picture. In this study, the average shear velocity was kept small. In reduced units, v* ≡ vτ/σ = 0.01 and 16789

DOI: 10.1021/acs.jpcc.9b03697 J. Phys. Chem. C 2019, 123, 16787−16795

Article

The Journal of Physical Chemistry C 60% of the channel’s width. The choice of total width Δv as a characteristic of boundary layers instead of separate widths (2) Δ(1) v and Δv results from the observation that they are not exactly equal due to fluctuations. Velocity profiles like those shown in Figure 2 are well-known in the literature; a detailed study of the influence of fluid density, channel width, fluid-wall potential, or temperature on the shape of the profiles in planar Couette flows was reported in ref 28 (see also ref 27). In our study, they serve as a starting point for a study of a manifestation of the local structure of the liquid in nanoscopic flows. 3.2. Drag Force. A direct measure of a channel’s resistance to a shear flow is the force F⃗ exerted by the wall on the liquid. The instantaneous value F⃗ (t) of this force strongly fluctuates. Averaging (smoothing) F⃗ (t) over a time θ = 106 h gives a weakly (

(4)

Figure 3 shows the plot of reduced viscosity η(ρ)/η(0.8) in the boundary layer against density. For low densities, this ratio

Figure 3. Plot of reduced viscosity η(ρ)/η(0.8) in the boundary layer. Solid line: rational fit; see text. 16790

DOI: 10.1021/acs.jpcc.9b03697 J. Phys. Chem. C 2019, 123, 16787−16795

Article

The Journal of Physical Chemistry C

Figure 4. Typical instantaneous equilibrium/flow configurations of the 2D LJ system at T = 2 for liquid: ρ = 0.8 (a), 0.92 (b), and liquidus line (ρ = 0.94 (c)). (d) Ideal crystalline wall (T = 2, ρ = 0.92). Red (green) circles represent SLA (LLA) atoms; see text.

0.555.47 This criterion is additionally validated by an observation that in computer simulations of a 2D LJ crystal only a small fraction of seven-atom clusters had Q6 < 0.555.39 The value Q6 = 0.555 corresponds to a crystal in close proximity to the solidus line. We refer to the central atom of a seven-atom cluster, recognized as crystalline (hexagonal), as a solid-like atom (SLA). Central atoms of clusters not classified as crystalline are referred to as liquid-like atoms (LLAs). Typical instantaneous configurations of a 2D Lennard-Jones liquid are shown in Figure 4, where SLA atoms are represented by red circles and LLA atoms by green circles. For ρ = 0.8, well below the liquidus density ρL, the fraction of SLA atoms is small and the system is a matrix of LLA material with few small islands of crystalline order (case a). On the contrary, close but below the liquidus line L (ρ = 0.92, case b) and on the liquidus line (ρ = 0.94, case c), SLA atoms aggregate and form connected clusters (crystallites) of crystalline-ordered matter. Thus, close to the liquidus line and in the coexistence region, a typical configuration of a 2D LJ liquid is a mosaic of crystallites of SLA atoms separated by amorphous matter of LLA atoms.39,40 The SLA component percolates very close or on the liquidus line. It is not immediately accompanied by the onset of long-range orientational order of the seven-atom SLA clusters; the latter takes place in a narrow interval of strip states : (ρS ≃ 0.96 for T = 2) and marks the onset of shear rigidity (solidity in terms of ref 48). In the narrow interval of strip states, the 2D LJ liquid becomes a complex system.40 The mosaic is a manifestation of diminished stability of the crystalline structure resulting from competition between attraction and repulsion forces.39 Crystallites and amorphous clusters do not constitute a thermodynamic phase but rather a dynamic, fluctuating structure (mosaic) that undergoes random but substantial changes on a time scale of τ0 ∼ τ;39 this time is the shortest in the hierarchy of structural relaxation times. The most basic parameter characterizing the local order in a 2D liquid is the equilibrium concentration c(b) 6 of SLA atoms in a large boundary-less (periodic) equilibrium liquid. The

system-average value of c(b) 6 is a function of state. The plot of c(b) 6 (T, ρ) for T = 2 is shown in Figure 5.

Figure 5. Concentration c(b) 6 as a function of density ρ for T = 2. Vertical lines denote liquidus line (L, red), solidus line (S, blue), and middle of strip states : (green).

5. STRUCTURE OF BOUNDARY LAYERS AT CHANNEL WALLS 5.1. The Profiles of Local Order and Density. To study the crossover of the symmetry of force field due to the walls, from perfect solid like (triangular) to a kind of disordered one, we freeze the wall atoms in their instantaneous positions corresponding to a fluctuating ideal triangular lattice with density equal to the current bulk density of the liquid. The fluctuations are Gaussian with rms equal to ξ (in units of σ).38 In the case ξ = 0, which corresponds to an ideal triangular crystal, the walls promote an increase of local crystal-like order close to themboth of concentration of SLA and the quality of local structures. For nonzero values of ξ, the amplitude of this constructive effect decreases, and for sufficiently large values of ξ, the destructive influence of the wall on the local structure may arise. Additionally, in both cases, the density can also be influenced by the walls. In fact, the stratification of the liquid normal to the wall, with strata separation of the order of molecular size, resulting from the interaction of the liquid with 16791

DOI: 10.1021/acs.jpcc.9b03697 J. Phys. Chem. C 2019, 123, 16787−16795

Article

The Journal of Physical Chemistry C

Figure 6. Density profile ρ̃(T, ρ; y) (a) and solid-like quality profile q6(T, ρ; y) (b) for ρ = 0.92. SLA concentration profiles c6(T, ρ; y) for ρ = 0.92 (c) and 0.80 (d). bulk, L, and strip mark the corresponding values for bulk, liquidus line, and strip states. T = 2, ideal wall (ξ = 0).

a planar wall for planar Couette flow is well-known.18,23,27,28 The important difference in our study is the very small (below the short-time fluctuation level) shear velocity gradient that excludes nonlinear effects leading to a possible strong stratification in temperature, density, and local structure. To quantify those predictions, we introduce the averaged over configurations/time profiles, calculated using the same partitioning (coarse-graining) of the liquid into horizontal layers and time-averaging as those described above for the velocity vx(y). The following coarse-grained profiles were studied: density ρ̃(T, ρ; y), concentration c6(T, ρ; y), and quality of local crystal order q6(T, ρ; y), the latter defined as the average of parameter Q6/Q6,max. We use 10 layers along the y direction which corresponds, on average, to Ly/10 = 2.5 atoms per layer. The coarse-graining of physical quantities requires some explanations. In particular, in the case of density, this procedure smooths out the aforementioned oscillations of the density (stratification) close to the wall (Figure 6a) and, in consequence, ignores some aspects of interfacial structure. The oscillations have the period of the distance between atoms, while coarse-graining accounts for the presence of the smallest element of local structure (SLA clusters). The size of an SLA cluster (approximately 2σ) introduces a new important local scale. Thus, we coarse-grain over a layer with a width approximately equal to 2.5σ. The higher resolution, with layer width smaller than 2σ, can lead to unnecessary artifacts. In particular, the SLA clusters near the wall are oriented, so the higher resolution density profiles show oscillations that can be, in principle, used to measure their orientation. This topic is beyond the scope of the current study. The coarse-grained profiles calculated for an ideal wall (ξ = 0) are shown in Figure 6 (we have found that the profiles in the static case and for the flow coincide within the statistical errorsit is a consequence of the relatively slow flow studied in this paper). Consider first the case ρ = 0.92, close to but below the liquidus line. The coarse-grained profiles a−c clearly mark, as

expected, the constructive effect of the wallsclose to them the density, the solid-like quality of local order, as well as the concentration of the SLA component are enhanced in comparison to their bulk values. The coarse-grained density profiles do not display the characteristic stratification effects, as discussed above. The quantitative analysis leads to an important conclusion. Namely, the density in the first layer (panel a) reads ρ̃1 ≈ 0.95. In the bulk system, this would result in an increase of the concentration to c6 = c(b) 6 (0.95) = 0.59 ± 0.01 (Figure 5). In the current case, the concentration in the first layer is noticeably higher: c6 = 0.85 ± 0.01 (panel c). This shows that the local concentration of the SLA component is modified by the walls in two ways. First, indirectly, as the result of the modification of the density ρ → ρ̃(T, ρ; y) (first term in eq 7). Second, directly, due to wall-induced external solid-like ordering potential (second term in eq 7): c6(T , ρ ; y) = c6(b)(T , ρ ̃(T , ρ ; y)) + δc6(T , ρ , ξ ; y)

(7)

An increase of c6 promotes the following qualitative effect. Close to the walls, the concentration of SLA atoms becomes higher than that on the liquidus line; in the first layer, it is higher than in the strip states (panel c). This shows that the walls promote the change of structure of sufficiently dense simple liquid close to the liquidus line to that of a complex one with enhanced tendency for creation of rigid local structures. The increase of the studied parameters close to the wall is accompanied by their decrease in the central area of the channel. In the case of density, it is a direct consequence of the conservation of the number of atoms; the density in the middle drops down to around 0.905 (panel a). The corresponding decrease of c6 in the bulk is estimated as follows. The linear approximation of the data from Figure 5 yields c(b) 6 (0.905) = 0.415 ± 0.01, in agreement with data from Figure 6c, where c6 = 0.42 ± 0.01. Thus, contrary to the previous case (proximity of the wall), the change of the concentration of the SLA component is driven by the change of the density, while the effect of the external ordering potential is negligible: c6(T, ρ; y) ≈ c(b) 6 (T, ρ̃(T, ρ; y)). 16792

DOI: 10.1021/acs.jpcc.9b03697 J. Phys. Chem. C 2019, 123, 16787−16795

Article

The Journal of Physical Chemistry C Δ ≈ Δv ≈ Δ c6

The general features of the profiles discussed above are representative for a range of densities around the liquidus line. For lower densities (ρ = 0.8), the wall-induced inhomogeneities become weak; see panel d. (2) (1) (2) The total width Δc6 = Δ(1) c6 + Δc6 of widths Δc6 and Δc6 of the upper and lower boundary layers with an enhanced local solid-like order is equal to the number of layers (for both walls) for which c6(T, ρ; y) > c(b) 6 (T, ρ). The characteristic uncertainty of this definition constitutes half of a layer. Figure 7 shows the plots of Δc6 and Δc6/Ly as a function of the density.

(8)

We have found that the widths Δ(T, ρ) for the 50 × 25 and 50 × 40 systems coincide within calculated uncertainties. Thus, the influence of the walls on the local structure of the liquid seems to be independent of the width of the nanochannel and the ratio Δ(T, ρ)/Ly becomes a universal parameter which characterizes the departure of the flow from hydrodynamic predictions: for Δ(T, ρ)/Ly ≪ 1, the flow is close to hydrodynamic predictions, while, for Δ(T, ρ)/Ly ≈ 1, the flow is strongly hindered in most of the channel’s width. 5.2. Influence of Density and Local Structure on Velocity Profiles. Having analyzed the flow and local structure of the liquid, we are ready to verify the main hypothesis which states that the local structure of the liquid, enhanced by the interactions with the walls, contributes, together with an inhomogeneous density profile, to the nonNewtonian flow in nanochannels. The interactions with the wallsthe origin of the non-Newtonian flowpromote an increase of both density and local crystallinity, the latter being only in part dependent on the increase of density; see eq 7. To provide the arguments in favor of this hypothesis, it would be necessary to study separately the influence of both factors on the flowa challenging topic because they are not fully independent. Instead, we produce a convincing plausibility argument. Namely, let us study the nonlinear flows in two cases: (i) c6 is not enhanced close to the wall and the density remains enhanced, and (ii) both c6 and the density are enhanced close to the wall. If the densities in both cases are approximately equal and the nonlinearity of the velocity profiles is weaker in the first case than in the second, then the local crystallinity (local structure) contributes to the nonlinearity. To verify this scenario, we study a specific case when the effect of the structure of the wall on the local structure of the liquid is missing, due to the incompatibility between local symmetry of the liquid and of the wall. To this end, we lower the quality of crystalline order of the wall by introducing the Gaussian fluctuations described above.

Figure 7. Total width Δ = Δc6 of both upper and lower boundary layers (in units of σ) as a function of density. The horizontal dashed line marks the width of the nanochannel. Inset: plot of Δ/Ly against density.

The boundary layer is approximately 5σ thick for low density and constitutes more than half of the nanochannel for ρ = 0.92. For still higher densities, it grows up rapidly, and for the density corresponding to the strip states (ρ = 0.96), it extends over the whole nanochannel. We find that the boundary layer widths estimated using the velocity profiles and local structure analysis have close values; e.g., for ρ = 0.94, Δv = 15σ and Δc6 = 17σ. We conclude that the non-Newtonian flow is closely related to enhanced local order in a boundary layer with width Δ:

Figure 8. Effect of Gaussian wall (ξ = 0.20) on the flow for T = 2, ρ = 0.92. (a) Instantaneous configuration, (b and c) profiles c6(T, ρ; y) and ρ̃(T, ρ; y), (d) velocity profiles for a few values of ξ shown in the figure. Red line: hydrodynamic prediction. 16793

DOI: 10.1021/acs.jpcc.9b03697 J. Phys. Chem. C 2019, 123, 16787−16795

Article

The Journal of Physical Chemistry C Figure 8 characterizes an influence of the Gaussian wall with ξ = 0.2 on the structure of the liquid and on its flow for ρ = 0.92. A comparison of instantaneous configuration (panel a) with its counterpart for ξ = 0 (Figure 4) clearly shows a destructive effect of the wall on the liquid’s structure close to it. Namely, the SLA atoms do not cluster on the fluctuated wall, contrary to the case of the nonfluctuating wall. This qualitative observation is supported by the profile c6(T, ρ; y) (panel b) which is constant (up to uncertainties): c6(T, ρ; y) = c(b) 6 (T, ρ). Contrary to the concentration profile, the density profile ρ̃(T, ρ; y) (panel c) remains inhomogeneous and reproduces the main features of its counterpart for ξ = 0 (Figure 6a). High density ρ̃ 1 in the first layer promotes an increase of concentration c6 in this layer; this effect is opposed by the destructive structural influence of the Gaussian wall (δc6 < 0 in eq 7), resulting in homogeneity of profile c6(T, ρ; y). Velocity profiles for this flow are shown in panel d. In the case of ξ = 0.2, the nonlinearity is still present but it is noticeably weaker than that for the crystalline wall (ξ = 0). This validates the aforementioned scenario. We conclude that the nonlinearity of the velocity profiles is partly enhanced by the local solid-like structure of the liquid, enhanced by the structured wall. All in all, this constitutes a manifestation of the local structure of liquids via directly measurable physical parameters.

and other strong nonlinear effects. These effects are interesting but need a separate study. The choice of temperature and density at no more than about 30−50% of SLA is very important. For this choice, at the length scale of the channel, the liquid behaves as a normal liquid. For higher densities, corresponding to a liquid between the mosaic band and solidus line, the liquid is, on the length scale of a narrow channel, totally crystalline and may be hexatically ordered; the atomic mechanisms of a shear flow in these states are substantially different from those in mosaic states. The study of those topics goes beyond the scope of the current paper. An important question is then, how does the physical picture in the 2D LJ liquid relate to real liquids? There is a large body of experimental facts that support the existence of local order in at least some 3D liquids, most importantly in the water at normal temperature and pressure. However, interactions in these liquids are much more complex than in the 2D LJ liquid, and recognizing and characterizing the local order in these liquids is still very challenging. Currently, to get insight into the local order in these liquids, one can only look for and study manifestations of this order in measurable properties. One assumes that, although interactions and thus the local order in a liquid are liquid-specific, generic manifestations of this order are similar for different liquids including the 2D LJ liquid. The possibility to simulate a 2D LJ liquid in states from significantly locally ordered to largely disordered, visualize the configurations, and easily recognize the ordered clusters makes this system a valid toy model to study those effects. Appearance of layers of altered structure near solid walls is one of these effects. Understanding the local structures in liquids is an ongoing task. An important part of this task is to define quantitative characteristics of these structures. We believe that studies of most simple but still locally ordered liquids can help in this task and hope that the framework developed in this study can be useful in studies of more complex liquids.

6. CONCLUSIONS To summarize, we have demonstrated an existence of a direct relation between the onset of the 2D non-Newtonian flow in nanochannels and the local structure of the liquid. This is, to the best of our belief, a first demonstration of the impact of the local structure of the liquid on directly measurable physical parameters of the liquid. The degree of commensurability of both structures, which can be tuned, e.g., by chemical engineering of the surfaces, influences the non-Newtonian flow. Full commensurability promotes an enhancement of the non-Newtonian flow in a boundary layer of the liquid with local structural properties different from those in the bulk of the liquidan increase of the concentration of solid-like atoms and the higher solid-like quality of small clusters is found. The width Δ(T, ρ) of the boundary layer constitutes a new spatial scale in nanofluidics. For sufficiently dense liquid, this scale can become comparable to the diameter of the nanochannel, resulting in strongly non-Newtonian flow. On the contrary, the walls which are incommensurable with the local structure of the liquid have a destructive effect on the width of the boundary layer and the velocity profiles approach the predictions of the macroscopic hydrodynamics. Moreover, while commensurable walls attract clusters of solid-like atoms, potentially resulting in jamming effects, incommensurable walls repel them and prevent the local-structure-driven clusterization on the walls. Those effects depend on the density. The non-Newtonian flow effects become weak for a low concentration of solid-like component of the liquid. For sufficiently dense liquid, its local viscosity in the boundary layer becomes very sensitive to small changes of the density. In our study, the long-time-averaged velocity gradient is orders of magnitude smaller than that found in the literature. We tested that larger gradients result in disruption of the local structure near the wall, with possible formation of layers substantially differing in physical properties, slip phenomena,



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

A. C. Mitus: 0000-0003-2316-780X Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We thank Dr. Creighton Thomas and Professor Monica Olvera de la Cruz for useful discussions and help in the early stages of this work.



REFERENCES

(1) Bocquet, L.; Charlaix, E. From Microfluidic Application to Nanofluidic Phenomena Issue. Chem. Soc. Rev. 2010, 39, 1073−1095. (2) Georges, J.-M.; Millot, S.; Loubet, J.-L.; Tonck, A. Drainage of Thin Liquid Films Between Relatively Smooth Surfaces. J. Chem. Phys. 1993, 98, 7345−7360. (3) Raviv, U.; Klein, J. Fluidity of Bound Hydration Layers. Science 2002, 297, 1540−1543. (4) Li, T.-D.; Gao, J.; Szoszkiewicz, R.; Landman, U.; Riedo, E. Structured and Viscous Water in Subnanometer Gaps. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 75, 115415. 16794

DOI: 10.1021/acs.jpcc.9b03697 J. Phys. Chem. C 2019, 123, 16787−16795

Article

The Journal of Physical Chemistry C (5) Maali, A.; Cohen-Bouhacina, T.; Couturier, G.; Aimé, J.-P. Oscillatory Dissipation of a Simple Confined Liquid. Phys. Rev. Lett. 2006, 96, 086105. (6) Becker, T.; Mugele, F. Nanofluidics: Viscous Dissipation in Layered Liquid Films. Phys. Rev. Lett. 2003, 91, 166104. (7) 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. (8) Bocquet, L.; Barrat, J.-L. Flow Boundary Conditions from Nanoto Micro-Scales. Soft Matter 2007, 3, 685−693. (9) Landau, L. D.; Lifshitz, E. M. Fluid Mechanics; Elsevier: Oxford, U.K., 2009. (10) Sui, H.; Han, B. G.; Lee, J. K.; Walian, P.; Jap, B. K. Structural Basis of Water-Specific Transport Through the AQP1 Water Channel. Nature 2001, 414, 872−878. (11) Thomas, J. A.; McGaughey, A. J. H. Reassessing Fast Water Transport Through Carbon Nanotubes. Nano Lett. 2008, 8, 2788− 2793. (12) Joseph, S.; Aluru, N. R. Why Are Carbon Nanotubes Fast Transporters of Water? Nano Lett. 2008, 8 (2), 452−458. (13) Goyon, J.; Colin, A.; Ovarlez, G.; Ajdari, A.; Bocquet, L. Spatial Cooperativity in Soft Glassy Flows. Nature 2008, 454, 84−87. (14) Pu, Q.; Yun, J.; Temkin, H.; Liu, S. Ion-Enrichment and IonDepletion Effect of Nanochannel Structures. Nano Lett. 2004, 4, 1099−1103. (15) Cottin-Bizonne, C.; Steinberger, A.; Cross, B.; Raccurt, O.; Charlaix, E. Nanohydrodynamics: The Intrinsic Flow Boundary Condition on Smooth Surfaces. Langmuir 2008, 24, 1165−1172. (16) Huang, D.; Sendner, C.; Horinek, D.; Netz, R.; Bocquet, L. Water Slippage Versus Contact Angle: A Quasiuniversal Relationship. Phys. Rev. Lett. 2008, 101, 226101. (17) Srivastava, D.; Santiso, E. E.; Gubbins, K. E. Pressure Enhancement in Confined Fluids: Effect of Molecular Shape and Fluid-Wall Interactions. Langmuir 2017, 33, 11231−11245. (18) Thompson, P. A.; Mark, O.; Robbins, M. O. Shear Flow Near Solids: Epitaxial Order and Flow Boundary Conditions. Phys. Rev. A 1990, 41, 6830−6837. (19) Long, Y.; Palmer, J. C.; Sliwinska-Bartkowiak, M.; Gubbins, K. E. Pressure Enhancement in Carbon Nanopores: a Major Confinement Effect. Phys. Chem. Chem. Phys. 2011, 13, 17163−17170. (20) Long, Y.; Palmer, J. C.; Coasne, B.; Sliwinska-Bartkowiak, M.; Jackson, G.; Mueller, E. A.; Gubbins, K. E. On the Molecular Origin of High-Pressure Effects in Nanoconfinement: The Role of Surface Chemistry and Roughness. J. Chem. Phys. 2013, 139, 144701. (21) Shi, K.; Gu, K.; Shen, Y.; Srivastava, D.; Santiso, E. E.; Gubbins, K. E. High-Density Equation of State for a Two-Dimensional Lennard-Jones Solid. J. Chem. Phys. 2018, 148, 174505. (22) Hansen, J. S.; Daivis, P. J.; Travis, K. P.; Todd, B. D. Parameterization of the Nonlocal Viscosity Kernel for an Atomic Fluid. Phys. Rev. E 2007, 76, 041121. (23) Yong, X.; Zhang, L. T. Slip in Nanoscale Shear Flow: Mechanisms of Interfacial Friction. Microfluid. Nanofluid. 2013, 14, 299−308. (24) Travis, K. P.; Todd, B. D.; Evans, D. J. Departure from NavierStokes Hydrodynamics in Confined Liquids. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1997, 55, 4288−4295. (25) Travis, K. P.; Gubbins, K. E. Poiseuille Flow of Lennard-Jones Fluids in Narrow Slit Pores. J. Chem. Phys. 2000, 112, 1984−1994. (26) Frydel, D.; Rice, S. A. Lattice-Boltzmann Study of the Transition from Quasi-Two-Dimensional to Three-Dimensional One Particle Hydrodynamics. Mol. Phys. 2006, 104, 1283−1297. (27) Guo, Z.; Zhao, T. S.; Shi, Y. Simple Kinetic Model for Fluid Flows in the Nanometer Scale. Phys. Rev. E 2005, 71, No. 035301(R). (28) Guo, Z.; Zhao, T. S.; Shi, Y. Generalized Hydrodynamic Model for Fluid Flows: From Nanoscale to Macroscale. Phys. Fluids 2006, 18, 067107.

(29) Patashinski, A. Z.; Mitus, A. C.; Ratner, M. A. Towards Understanding the Local Structure of Liquids. Phys. Rep. 1997, 288, 409−434. (30) Hansen, J.-P.; McDonald, I. R. Theory of Simple Liquids with Applications to Soft Matter; Academic Press: San Diego, CA, 2013. (31) Brazhkin, V. V.; Voloshin, R. N.; Popova, S. V. Metastable States in Simple Melts: the Kinetics of Transitions in Liquid Se and S Under High Pressure. Phys. Lett. A 1992, 166, 383−387. (32) Brazhkin, V. V.; Voloshin, R. N.; Popova, S. V.; Umnov, A. G. Nonmetal-Metal Transition in Sulphur Melt Under High Pressure. Phys. Lett. A 1991, 154, 413−415. (33) Mitus, A. C.; Smolej, F.; Hahn, H.; Patashinski, A. Z. Q446”Shape Spectroscopy” of Local f.c.c. Structures in Computer Simulations of Crystallization. Europhys. Lett. 1995, 32 (9), 777−782. (34) Mitus, A. C.; Smolej, F.; Hahn, H.; Patashinski, A. Z. Theory and Practice of ”Shape Spectroscopy” of Local FCC Structures in Computer Simulations of Nucleation and Crystallization. Phys. A 1996, 232, 662−685. (35) Steinhardt, P. J.; Nelson, D. R.; Ronchetti, M. BondOrientational Order in Liquids and Glasses. Phys. Rev. B 1983, 28, 784−805. (36) Klumov, B. A.; Khrapak, S. A.; Morfill, G. E. Structural Properties of Dense Hard Sphere Packings. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 184105. (37) Klumov, B. A. Structural Features of a Lennard-Jones System at Melting and Crystallization. JETP Lett. 2013, 97, 327−332. (38) Mitus, A. C.; Patashinski, A. Z.; Patrykiejew, A.; Sokolowski, S. Local Structure, Fluctuations, and Freezing in Two Dimensions. Phys. Rev. B: Condens. Matter Mater. Phys. 2002, 66, 184202. (39) Patashinski, A. Z.; Orlik, R.; Mitus, A. C.; Grzybowski, B. A.; Ratner, M. A. Melting in 2D Lennard-Jones Systems: What Type of Phase Transition? J. Phys. Chem. C 2010, 114, 20749−20755. (40) Patashinski, A. Z.; Ratner, M. A.; Grzybowski, B. A.; Orlik, R.; Mitus, A. C. Heterogeneous Structure, Heterogeneous Dynamics, and Complex Behavior in Two-Dimensional Liquids. J. Phys. Chem. Lett. 2012, 3, 2431−2435. (41) Patashinski, A. Z.; Orlik, R.; Mitus, A. C.; Ratner, M. A.; Grzybowski, B. A. Microphase Separation as the Cause of Structural Complexity in 2D liquids. Soft Matter 2013, 9, 10042−10047. (42) Xu, X.; Rice, S. A.; Dinner, A. R. Relation Between Ordering and Shear Thinning in Colloidal Suspensions. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 3771−3776. (43) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Oxford University Press: 1987. (44) Todd, B. D.; Daivis, P. J. Nonequilibrium Molecular Dynamics: Theory, Algorithms and Applications; Cambridge University Press: 2017. (45) Evans, D. J.; Holian, B. L. The Nose-Hoover Thermostat. J. Chem. Phys. 1985, 83, 4069−4074. (46) Strandburg, K. J. Two-Dimensional Melting. Rev. Mod. Phys. 1988, 60, 161−207. (47) Mitus, A. C.; Marx, D.; Sengupta, S.; Nielaba, P.; Patashinskii, A. Z.; Hahn, H. Locating Liquid-Solid Transitions in Computer Simulations Based on Local Structure Analysis. J. Phys.: Condens. Matter 1993, 5, 8509−8522. (48) de Souza, V. K.; Harrowell, P. Liquids and Structural Glasses Special Feature: Rigidity Percolation and the Spatial Heterogeneity of Soft Modes in Disordered Materials. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 15136−15141.

16795

DOI: 10.1021/acs.jpcc.9b03697 J. Phys. Chem. C 2019, 123, 16787−16795