Meng Liao

Quy-Dong To
email: quy-dong.to@u-pem.fr

Céline Léonard

Vincent Monchiet

Van-Hoang Vo

Strain-induced friction anisotropy between graphene and molecular liquids

Keywords:

In this paper, we study the friction behavior of molecular liquids with anisotropically strained graphene. Due to the changes of lattice and the potential energy surface, the friction is orientation dependent and can be computed by tensorial Green-Kubo formula. Simple quantitative estimations are also proposed for the zero-time response and agree reasonably well with the Molecular Dynamics results. From simulations, we can obtain the information of structures, dynamics of the system and study the influence of strain, molecular shapes on the anisotropy degree. It is found that unilateral strain can increase friction in all directions but the strain direction is privileged. Numerical evidences also show that nonspherical molecules are more sensitive to strain and resulting more pronounced anisotropy effects.

I. INTRODUCTION

In micro-nanofluidic systems, the fluid-solid interface contributes a significant part in the overall behavior of the system [1-3]. Generally, the whole system can be modelled by the conventional macroscopic hydrodynamic equations combined with Navier boundary equations for the slip velocity v s at the wall. For a Newtonian fluid of viscosity η, the latter can be written in two equivalent forms

σ = λv s , v s = L s ∂v ∂n (1)
where σ = η∂v/∂n is the viscous shear stress, ∂v/∂n the normal derivative of fluid velocity v at the wall. The interface constants, λ the friction coefficient, η the viscosity and L s the slip length, are related by

L s = η λ (2)
The slip effects which depend on the nature of the fluid solid interaction, are enhanced when the hydrophobicity is involved. Both experiments with advanced techniques and computer simulations have provided supporting evidences of the phenomenon [4][5][6]. Atomistic simulations are generally based on either two techniques: Non Equilibrium Molecular Dynamics (NEMD) and Equilibrium Molecular Dynamics (EMD). The former concerns simulations of fluid flow and compute the slip length by comparing with Navier Stokes solution [START_REF] Karniadakis | Microflows and nanoflows: Fundamentals and simulation[END_REF][START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF][START_REF] Allen | Computer Simulation of Liquids[END_REF][START_REF] Frenkel | Understanding Molecular Simulation: From Algorithms to Applications[END_REF].

The latter is founded on the linear response theory and the friction coefficient is determined via Green-Kubo formula [START_REF] Chandler | Introduction to Modern Statistical Mechanics[END_REF][START_REF] Bocquet | [END_REF][START_REF] Evans | Statistical mechanics of nonequilibrium liquids[END_REF].

Regarding microfluidics systems involving water and carbon nanotubes (CNTs) or graphene based nano channels, experiments revealed that the water flow rates are many times higher than prediction using no slip boundary conditions [START_REF] Majumder | [END_REF][14][15]. Despite some scattering results in literature, Molecular Dynamics (MD) simulations agree with the exceptionally small friction found on CNTs and graphene [5,16,17]. Most of those results are obtained for pristine graphene or CNTs in laboratory while real graphene can have defects and be affected by other environmental mechanical conditions. Graphene can be subject to mechanical strain due to different reasons, for example lattice mismatch or thermal expansion in graphene/substrate systems. Graphite, that has similar structure as graphene, can be found in profound coal bed, under large compressive stress states. Atomistic simulations have shown that the friction can increase significantly with strain. For example, by MD simulations on graphene/water system, Xiong et al. [18] have observed the variation of slip effects by applying isotropic strain on the graphene sheet. Since graphene and graphite have many important applications, it is crucial to understand how the strain affects the slip effect.

The present work investigates graphene (or one graphite layer) and the interaction with different liquids namely water (H 2 O), carbon dioxide (CO 2 ) and methane (CH 4 ). Those fluids have different molecular shapes and important applications in technology (membrane and nanofluidics systems), energy and environment (carbon dioxide sequestration and methane recovery). Graphene or graphite can be naturally or artificially strained, which result anisotropic friction, i.e the friction along the strain direction is different from the other direction. This type of behavior can be described by a tensorial Green-Kubo expression which will be examined in this paper. The structure and the dynamics of those molecular liquids are also studied to find how they contribute to the friction. Simulation results show that strain has increased the friction but more dominantly along the strain direction. Non spherical molecules are more sensitive to those changes and relaxation time contributes a large part in the friction. Details of these findings will be presented in the following sections.

II. FRICTION AND SLIP TENSORS FOR ANISOTROPIC SURFACES

A. Green-Kubo expression for the friction tensor

Let us consider the situation where the fluid flowing past a solid surface of area A where the slip velocity v s is linearly proportional to the tangential friction force

F F = AΛv s . (3)
Based on linear response theory, Boquet and Barrat [START_REF] Bocquet | [END_REF]19] proposed to use the Green-Kubo formalism [20] to compute the friction coefficient. In the general situations, the expression for the friction tensor Λ is the following [21]

Λ = 1 Ak B T ∞ 0 F (0) ⊗ F (t) dt (4)
In ( 4), the ensemble average, notation .. , of the force correlation function is taken in the equilibrium state and integrated with time. Since the correlation is decaying with time, tensor Λ can be further decomposed into a static term and time decorrelation tensor τ t , for example

Λ = F (0) ⊗ F (0) Ak B T τ t , τ t = F (0) ⊗ F (0) -1 ∞ 0 F (0) ⊗ F (t) dt (5)
For isotropic-like surfaces, Λ can be replaced with a single friction coefficient λ [START_REF] Bocquet | [END_REF]19] and

τ t by the decorrelation time τ t λ = F (0) 2 Ak B T τ t , τ t = ∞ 0 F (0)F (t) dt F (0) 2 (6)
For surfaces with two axes of symmetry, say Ox and Oy, tensor Λ is diagonal in the Oxy system with principal values λ xx and λ yy . In this paper, we shall investigate the behavior of the friction tensor in the modelling of anisotropic surfaces.

Relations between friction coefficient and the fluid structure have been studied in the past. In particular, for FCC (100) crystal surfaces [22] or honeycomb lattice like graphene [16], simple working definitions of the friction coefficient at the solid-liquid interface can be derived. We note that the same theory can be extended to any Bravais lattice structures.

Given a surface lying on the plane xOy, the components λ αβ (α, β = x, y) of the friction tensor can be rewritten in the following form

λ αβ = 1 Ak B T ∞ 0 dt dzdz dxdx F α (x , z )F β (x, z) ρ(x, z, 0)ρ(x , z , t) (7)
with ρ(x, z, t) being the microscopic fluid density at planar coordinate x and distance z from the wall [START_REF] Hansen | Theory of simple liquids[END_REF]. In special cases where x and y are already axes of symmetry of the surface, the off-diagonal friction coefficients vanish, λ xy = 0, and it is sufficient to determine the principal values, λ xx and λ yy , of the friction tensor. The theoretical investigation of those friction coefficients will be detailed in the following section, with application to strained graphene system.

B. Approximation of fluid-wall potential for strained graphene

Let us assume that by some reason, the graphene sheet is uniformly strained in direction

x (armchair) and/or y (zigzag) [18] but still periodic. There are several ways of defining the unit cell and the pair of primitive reciprocal lattice vectors, for example (see Fig. 1).

FIG. 1. Lattice vectors and reciprocal vectors

(a 1 , a 2 ) → (q 2 , q 1 ), (a 1 , a 3 ) → (q 3 , q 1 ), (a 2 , a 3 ) → (q 3 , q 2 ),

with a 1 , a 2 , a 3 being the primitive lattice vectors and q 1 , q 2 , q 3 associated primitive reciprocal lattice vectors a 1 = l 0 (cos ϕe x + sin ϕe y ), a 2 = l 0 (cos ϕe x -sin ϕe y ), a 3 = 2l 0 sin ϕe y q 1 = q 0 e x 2 cos ϕ -e y 2 sin ϕ , q 2 = q 0 e x 2 cos ϕ + e y 2 sin ϕ , q 3 = q 0 e x cos ϕ .

The quantity ϕ is half angle between a 1 and a 2 , l 0 is the distance between hexagon centers as shown in Fig. 1 and q 0 = 2π/l 0 . In all cases, the following properties must hold

q i • a i = 0, q i • a j = ±2π (10)
We consider first the case of monatomic fluid and note that results for molecules can be obtained approximately by superpositions as suggested in Ref. [START_REF] Falk | [END_REF]. The analytical expression of friction force F (x, z) of a fluid atom with the wall is derived from the potential V (x, z).

To the first order Fourier series approximation (see Appendix B), we can write

V (x, z) = V 0 (z) -V 1 (z)[cos(q 1 • x) + cos(q 2 • x)] -α(z)V 1 (z) cos(q 3 • x). (11)
Depending on the choice of the unit cell, the original Fourier series contains only two of three terms cos(q 1 • x), cos(q 2 • x) or cos(q 3 • x). Here we adopt an expression that can account for the periodicity of the potential along all directions a 1 , a 2 and a 3 and the anisotropy effect via the parameter α(z). Comparing with the expression of Ref. [16]:

V (x, z) = V 0 (z) -V 1 (z)[cos(q 1 • x) + cos(q 2 • x)], (12)
the present expression is able to recover the symmetry of the graphene structure, which is important for the isotropy/anisotropy analysis (see Fig. 2). For given z the values of V 0 (z), V 1 (z) and α(z) are determined by fitting the exact results at 3 representative points, for example

V (0, z) = V (a 1 , z) = V (a 2 , z) = V (a 3 , z) = V 0 (z) -(α(z) + 2)V 1 (z), V (a 1 /2, z) = V 0 (z) + α(z)V 1 (z), V (a 3 /2, z) = V 0 (z) -(α(z) -2)V 1 (z). ( 13
)
C. Analytical model for fluid-wall interaction

After constructing the analytical expression for the potential, we can differentiate the latter with respect to x and y to derive the tangential force components

F x = f x (z)[sin(q 1 • x) + sin(q 2 • x) + 2α(z) sin(q 3 • x)], F y = f y (z)[sin(q 1 • x) -sin(q 2 • x)] f x (z) = V 1 (z)q 0 2 cos ϕ , f y (z) = -V 1 (z)q 0 2 sin ϕ (14)
Here V 1 (z)(1 + α(z))/2 and V 1 (z) are the energy barriers along direction a 1 (or a 2 ) and a 3 .

As a result, we can rewrite [START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF] in the form

λ xx = 1 Ak B T ∞ 0 dt dzdz f x (z)f x (z ) dxdx [sin(q 1 • x ) + sin(q 2 • x ) + 2α sin(q 3 • x )] [sin(q 1 • x) + sin(q 2 • x) + 2α sin(q 3 • x)] ρ(x, z, 0)ρ(x , z , t) (15)
and a similar expression for λ yy . Following the approach presented in Ref. [22], we make use of the Fourier transform and obtain the simple expression for λ xx

λ xx = q 2 0 4 cos 2 ϕAk B T ∞ 0 dt dzdz V 1 (z)V 1 (z ) × × { ρ q 1 (z)(0)ρ -q 1 (z )(t) + 2α(z)α(z )ρ q 3 (z)(0)ρ -q 3 (z )(t) } (16)
and for λ yy

λ yy = q 2 0 4 sin 2 ϕAk B T ∞ 0 dt dzdz V 1 (z)V 1 (z ) { ρ q 1 (z)(0)ρ -q 1 (z )(t) } (17)
In deriving (16,17), we assume that the fluid is homogeneous in the plane xOy and the symmetry between q 1 and q 2 . The notation stands for the real part and ρ q is given by the expression

ρ q (z, t) = dxe -iq.x ρ(x, z, t) (18)
It is clear that for perfect graphene surface where ϕ = π/6, α(z) = 1 and the roles of q 1 , q 2 , q 3 can be interchanged, one can show that λ xx = λ yy . This property can not be obtained using potential [START_REF] Evans | Statistical mechanics of nonequilibrium liquids[END_REF], as done in previous works. When only the first fluid layer at coordinate z 0 is considered, one can obtain the simple relations

F x (0)F x (t) q 2 0 V 2 1 (z 0 ) 4 cos 2 ϕ N (z 0 )[F(q 1 , z 0 , t) + 2α 2 (z 0 )F(q 3 , z 0 , t)] F y (0)F y (t) q 2 0 V 2 1 (z 0 ) 4 sin 2 ϕ N (z 0 )F(q 1 , z 0 , t) (19)
where F(q, z 0 , t) is the planar intermediate scattering function of N (z 0 ) molecules

F(q, z 0 , t) = 1 N (z 0 ) ρ q (z 0 )(0)ρ -q (z 0 )(t) (20)
Substituting t = 0 in the above equation yields the zero-time behavior and the relation with the planar structure factor S(q k , z 0 ) and the layer density

ρ s (z 0 ) = N (z 0 )/A F 2 x (0) /A q 2 0 V 2 1 (z 0 ) 4 cos 2 ϕ ρ s (z 0 )[S(q 1 , z 0 ) + 2α 2 (z 0 )S(q 3 , z 0 )] F 2 y (0) /A q 2 0 V 2 1 (z 0 ) 4 sin 2 ϕ ρ s (z 0 )S(q 1 , z 0 ) (21)
The 2D structure factor S(q, z 0 ) of the layer of the first liquid layer can be calculated by using the following expression [16]:

S(q, z 0 ) = 1 N (z 0 ) N j=1 cos(q • r j ) 2 + N j=1 sin(q • r j ) 2 (22)
In (22), N is the number of fluid molecules in the first fluid layer, r j is the coordinates of the jth fluid atoms. Although using equation ( 21) can reproduce correctly some important phenomena and tendencies, it considerably underestimates F 2

x and F 2 y by an order of magnitude. In reality, the first liquid layer is not truly localized in a plane, as considered by the theory, but rather spreads over a finite thickness. Using V (z 0 ) to represent the whole depth and the presence of many liquid layers can be responsible for those differences. To improve these issues, a possible treatment is to use the assumption of independent multi liquid layer at z 0 , z 1 , .., which results the equations

F 2 x (0) q 2 0 4 cos 2 ϕ i N (z i )V 2 1 (z i )[S(q 1 , z i ) + 2α 2 (z i )S(q 3 , z i )],
F 2 y (0)

q 2 0 4 sin 2 ϕ i N (z i )V 2 1 (z i )S(q 1 , z i ), (23)
in which calculation of structure factor must be done at different levels. In the present paper, we propose using a simpler approximate formula

F 2 x (0) q 2 0 4 cos 2 ϕ [S(q 1 , z 0 ) + 2α 2 (z 0 )S(q 3 , z 0 )] dzρ(z)V 2 1 (z), F 2 y (0) q 2 0 4 sin 2 ϕ S(q 1 , z 0 ) dzρ(z)V 2 1 (z). ( 24
)
This expression is also computationally simple while keeping all the important ingredients as before. As noted previously, the above analysis is applied to monatomic fluid model. For fluids composed of molecules, as example H 2 O or CO 2 , the full expression F 2 α (0) is the following

F 2 α (0) = i (F (i) α (0)) 2 + i =j F (i) α (0)F (j) α (0) ( 25
)
where the summations are done over the different atoms composing the molecule of the fluid.

In a more recent work [START_REF] Falk | [END_REF], numerical evidences have shown that we can neglect the cross correlation between the friction of different atom species. As a result, we can compute separately the friction of each species and then simply superpose,

F 2 α (0) i (F (i) α (0)) 2 (26)
The accuracy of those simplifications will be examined in comparison with the exact results from MD simulation presented in the next section.

III. MOLECULAR DYNAMICS SIMULATION A. Choice of systems and potentials

To study the friction theory presented previously and the influence of molecular sizes and shapes, and interactions, we considered different liquids confined between two strained graphene sheets. Those liquids are water (H 2 O, triangular shape), carbon dioxide (CO 2 , rod like shape) and methane (CH 4 , spherical particle), which are present abundantly in nature and involved in many technological, energetic and environmental problems. For example, the water-graphene interaction arises in fast transport nanofluidic systems with applications in desalinated membrane industry. Carbon dioxide sequestration and methane recovery are directly related to exploitation of natural gas reserve in coal bed while reducing green house effect. The presence of underground water is also an issue that must be taken into account.

In summary, the study of those systems from the atomic scale is of both theoretical and technological importances.

The choice of the different pairwise and many-body interaction potentials is crucial with a strong impact on the MD results. The fluid-wall interactions are modeled by Lennard-Jones (LJ) potentials whose parameters are given in Table I.

For the graphene sheet, denoted as G, we have used the adaptive intermolecular reactive bond order (AIREBO) potential [25] for the interaction between the carbon atoms.

In addition to the classical reactive empirical bond-order (REBO) functional form, LJ and 4-body torsional interactions are taken into account.

The model TIP4P-2005 [26] has been used for water, constituting of four rigid sites, three fixed point charges and one LJ center. The H-O-H bond angle 104.52 • and the O-H bond length 0.9572 Å are maintained using SHAKE algorithm [27]. It is worthy emphasizing that unlike many previous works on water/graphene systems, we take into account the realistic interaction between G-H based on Ref. [28] with LJ parameters provided in Table I. Numerical results in the later section also show that the friction coefficients change significantly and become closer to the experimental values in Ref. [15] with the use of G-H interaction.

Fluid-Fluid Fluid-Wall The CH 4 molecules are modeled as united spherical particles with the TraPPE force field [31] and the fluid-wall potential is also determined by Lorentz-Berthelot rules with the same σ G and G as for the CO 2 -Graphene interaction. All LJ parameters in CH 4 -Graphene model are presented in Table I.

σ f f [ Å] f f [meV] σ f w [ Å] f w [meV] H 2 O O-O 3
Since interaction potentials between the fluid molecules and graphene sheets are important for the derivation of meaningful results, we carry out some verifications. The global interaction potential between a fluid molecule and the graphene sheet have been recomputed from the pairwise LJ using parameters of Table I. The corresponding potential well-depths and equilibrium distances have then been compared with bibliographic entries. For some liquids, experimental data only exist for the molecule-graphite interaction. In that cases, we have also computed the global interaction potentials with graphite. The graphene slab is composed of 256 atoms and the graphite surface contains 8000 atoms divided in 5 layers.

Periodic boundary conditions in the x and y directions are applied in both cases. The C-C distance has been fixed to 2.46 Å for both structures. The geometries of H 2 O, CH 4 and CO 2 are fixed to those of the MD simulations. The corresponding results are displayed in Table II.

For H 2 O-Graphene, the potential parameters have been derived by Hugues et al. [28] from extensive plane-wave DFT calculations using the revPBE-vdW-DF functional in order to produce an efficiently implemented polarisable force-field (GRAPPA). The potential well depth and equilibrium geometry, obtained for the adsorption conformation where both H water atoms pointed towards a carbon atom of graphene, compare well with the results of coupled-cluster (DFT/CC) calculations [33].

V eq [meV] Ref.

Z For CO 2 and CH 4 , only their interactions with graphite have been previously reported to the best of our knowledge. The best estimate potentials have been reported from averaged interaction potential as deduced from analyses of experimental data and calculations by Vidali et al. [32]. The corresponding potential well depth and equilibrium geometry of the CO 2 -graphite system agree with the present values obtained for the parallel bridge conformation, that has been confirmed to be the lowest energy conformation by Xu et al. [34] with DFT calculations for the CO 2 -(44) six-ring aromatic surface system.

The present LJ potential for CH 4 -G leads to a hollow lowest energy conformation, i.e.

the CH 4 molecule preferentially adsorbs towards the center of a C ring, which is one of the lowest energy adsorption site found also by Xu et al. [34] for the CH 4 -(44) six-ring aromatic surface system. The recomposed global interaction potential of CH 4 with graphite also presents potential well depth and equilibrium geometry in good agreement with the Vidali et al. [32] best estimations.

From the present analysis, we can conclude that the LJ parameters given in Table I allow to reproduce the global interactions of the different molecules with graphene and the corresponding pairwise potentials can be safely used in MD simulations devoted to interface effects.

B. System setup

All our molecular dynamic simulations and friction calculations are done in equilibrium state using LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) package [35]. The systems, periodic in x, y, are first relaxed for 5 × 10 6 time step at zero lateral external pressure in order to reach the natural reference state. Next, the simulation box is deformed in x or y direction in accordance with the desired strain for graphene. In the present work, we keep the direction y (zigzag) undeformed while the strain in direction x (armchair) is varied from 0% to 10%. The systems undergo equilibration process again and only after 5 × 10 6 time step, we start to compute the auto correlation function based on data of the next 10 × 10 6 time step. During the whole simulation, an external uniform pressure P ext is applied to graphene sheet to equilibrate the liquid pressure inside the channel. Additionally, to avoid drifting, one atom of the lower graphene sheet is fixed. The simulation time step is 1 fs. The temperature of the graphene is maintained by Nose Hoover thermostat and the equations of motions of liquid atoms are integrated using Verlet algorithm. In order to obtain the liquid density for fluids, we list different temperature and different external pressure for the different systems in Tab III. To improve the reliability, the results are averaged over at least 10 independent runs of the same system. Since the immobile graphene nanochannel may lead to non-physical results, flexible graphene sheets have been used in all simulations. The periodic boundary conditions are adopted in

Liquids N ρ [nm -3 ] T [K] P ext [
x and y. The dimension of the graphene surface A is 42.6 Å × 24.6 Å and the height of nano-channels has been adjusted by external pressure and temperature of the given fluid.

A snapshot of our systems can be seen in Fig. 3.

In water-graphene simulations, the system is composed of 1001 water (H 2 O) molecules between two graphene sheets. For simulations involving methane (CH 4 ) in graphene chan- nel, 720 liquid molecules are considered. Regarding carbon dioxide-graphene system, there are 352 CO 2 molecules in the simulation box.

C. Results and discussion

The variation of the liquid density with respect to distance from the graphene surface can be seen in Fig. 4. We find that for spherical particles like CH 4 , the density profile is almost insensitive to the applying strain whereas for the most aspherical molecules, i.e CO 2 which adopts a rod shape, the influence is more visible. Specifically, both the first and the second peaks corresponding to the first and second liquid layers decrease. The effect is more pronounced for CO 2 simulations for which the second peak shifts slightly towards the bulk.

For all three different fluids, the positions of the first peak are however in agreement with the fluid-wall LJ distance, σ f w .

From the density profile, intuitively, one can suppose that less liquid molecules are present near the wall smaller will be friction of strained graphene. However, this prediction is disproved by detailed results presented in the following. The anisotropy effect, which can not be shown from the density profile, must also be studied.

The exact value of friction tensor can be computed using the time correlation presented in Fig. 5. To obtain the friction coefficients by the Green Kubo expression, the correlations F α (0)F β (t) are accumulated continuously until convergence. For given t, at least 800000 samples are collected for the ensemble average. Finally, results are averaged again over 10 independent simulations. The difference in behavior between strained and natural graphene can be observed, especially along direction x. All curves show an exponentially decaying behavior but they are quantitatively different. The zero-time value F 2 (0) at t = 0, representing the static friction intensity, is significantly higher for strained graphene and the correlation decays more slowly in that case. We will see later that the correlation time increases with the strain value. The same behavior can also be observed for friction along y but the change is less pronounced since the graphene has zero strain along y. The anisotropy effect due to unequal strain can also be seen from the x, y correlation behaviors which are identical for natural graphene and very different for strained graphene. Fig. 6 and Tab. IV shows the friction coefficient results for different fluid types at different strains. The isotropic behavior can be clearly seen for all types of liquid at zero strain and the anisotropy effect starts to increase as the strain increases. The friction increases for both directions but the increase is more important for direction x. It is interesting to note that the relation between λ xx and λ yy , and the applying strain is quasi-linear up to strain as large as 10%. The difference between λ xx and λ yy is more significant for non-spherical molecules like H 2 O and CO 2 , and is reduced for spherical molecules like CH 4 . Due to the direct connection between slip length and friction, the slip effect is minimum, maximum along x, y direction, respectively. The sensitiveness of non-spherical molecules suggests that it is possible to enhance or reduce the transport performance using engineering strain in one or two directions. However, it is less effective to use this method for spherical molecules. The viscosity for the TIP4P/2005 water model is 0.855 mPa.s at 298 K and 1 bar [36]. In order to compare our results with other researches, we converted all friction coefficients λ of water to the slip lengths L s using relation [START_REF] Karniadakis | Microflows and nanoflows: Fundamentals and simulation[END_REF]. At the strain-free state, the present L s value is 30.5 nm which is consistent with the value of about 30 nm obtained by Thomas et al. [37] and which is also close to the theoretical value of Myers, i.e. 35 nm [38]. In the recent experiment of Secchi et al. [15], measurements on very large radius (33 nm to 50 nm) CNTs show that L s varies from 45 nm to 17 nm, which is comparable with the present results. It can be concluded that accounting for the G-H interaction and the flexibility of the graphene sheet are very important. The present results show that the typical L s values are two times smaller than those resuting from simulations neglecting these aspects. More importantly, they are closer to experimental results. For example in Fig. 7, at all strain values, the present L s along direction x is significantly smaller than the isotropic results (same strain along x and y directions) of Xiong et al. [18]. Since the friction increase due to strain can be approximately decomposed by Eqs (5,[START_REF] Falk | [END_REF], we shall look at the influence of different terms in the overall behavior, including the static terms F x (0) , F y (0) depending on the structure factor S(q) and the integral involving potential strength and the decorrelation time τ x , τ y . Fig. 8(a) shows that the quasi-isotropic structure is observed at the first fluid layer at the maximal strain 10%. The changes of structure factor due to strain can be considered negligible Fig. 8(b), less than 2%. This remarks suggests that the contribution of S(q) on the anisotropy is not significant. The norm of reciprocal vectors q 1 , q 2 and q 3 at different strain states are presented in Table V.

0 2 4 6 8 10
Strain (%) 2 3.2 4.4 5.6 6.8 8 λ
(10 4 N • s/m 3 ) H 2 O (a)
Fig. 9 represents the dependence of coefficients α and V 1 in terms of the distance z of a methane molecule from the graphene sheet. As can be seen, both strain and z can affect α and V 1 , but the strain mainly contributes to the linear change of α while the potential ener- The structure factor of oxygen and hydrogen of water with different q. (c) The structure factor of oxygen (H 2 O), carbon (CO 2 ) and methane (CH 4 ) at strain-free state.

gies coefficient V 1 is very sensitive to z. The variation of V 1 (z) in the range, z ∈ [2.2, 3.6]( Å),
q 1 [ Å-1 ] q 3 [ Å-1 ] 0 2.
= q 2 .
better agrees with exponential form V 1 (z) = 3.031 × 10 4 exp{-5.797z}(eV) for strain free case and V 1 (z) = 2.104 × 10 4 exp{-5.67z}(eV) for 10% strain. It suggests that using the representative value at z 0 may not yield the best estimation. In this paper we consider the potential form of V 1 when calculating the static term in [START_REF] Falk | [END_REF].

In Fig. 10, we can find that the static forces F 2

x /A, F 2 y /A, and the decorrelation time depend on the strain. It is noted that the spherical molecules (CH 4 ) are less sensitive to the unilateral strain, which agrees with the friction results shown in Fig. 6. The variation of those quantities for H 2 O and CO 2 is much more significant. The theoretical results from [START_REF] Falk | [END_REF] have well predicted the variation trend and the value of F 2

x /A with strain. For F 2 y /A, there are still difference between theoretical and simulation results. The simulation results show slight changes of F 2 y /A while it is more visible according to the theoretical prediction. The discrepancies could be clearly understood from the analysis of the assumptions: i) the use of a time independent potential energy surface obtained from static graphene ii) the superposition hypothesis for molecules in Eq. ( 26) iii) the fact that we only account for the structure factor in the first liquid layer. The first assumption, adopted in most theoretical works [16,22,[START_REF] Falk | [END_REF], corresponds to situations where flexibility of the graphene sheet is negligible, e.g at small temperature or at large strain. Physically, stretching membrane reduces the out-of-plane vibration amplitude, which agrees with the reduced difference between MD solution and the theory at large strain for both F 2

x /A and F 2 y /A. The second superposition hypothesis have been tested for large alkanes [START_REF] Falk | [END_REF] and proved to yield good estimates.

Finally, the assumption that allows reducing the static part of ( 16) to ( 24) is relatively strong but necessary to obtain important connections to the molecular distribution and the interaction strength, appeared as two separated contributions. In [START_REF] Falk | [END_REF], although we does consider the variation of potential with depth for the whole slab, but only the structure factor of the first liquid layer is used for calculation. As a result, potential improvements to the model can be done from the analysis of second and third liquid layers based on the same procedure, in order to reduce further the gap between simulation and theory results.

Regarding the decorrelation time, numerical evidences have shown that it plays important role in the increase of the friction due to strain. Indeed, while the static term F 2

x /A is responsible for 10% increase, the decorrelation time can be doubled for H 2 O. It suggests that the motion of molecules has been affected by the changes of the environment. From atomistic viewpoint, the changes of lattice distances induce changes in the potential energy landscape and the diffusion mechanism of molecules. Spherical molecules like CH 4 are less affected, but molecules like H 2 O and CO 2 tend to have less mobility along x than along y. Strain (%) Strain (%) In the present paper, we have considered the friction between different liquids and an anisotropic surface, e.g graphene subject to anisotropic engineering strain. Due to the changes of the lattice structure, the potential field of strained graphene has lost its six fold symmetry and is responsible for anisotropic friction behavior between the liquid and the surface. Depending on the molecular shape and interaction strength, the anisotropy degree may vary from one fluid species to another. Using LAMMPS software, one can compute the friction tensor via time correlation integral and access to structure via post process routine.

Simple predictive estimation is proposed from a constructed surface potential for strained graphene.

The authors study three different liquids with distinct molecular shape, namely water, carbon dioxide and methane. Numerical evidences show that the strain induced anisotropy effect is significant, especially for non spherical molecules. When strain increases, the friction and the anisotropy degree also increase. For example, the friction coefficient can rise 100%

for water and 50% for carbon dioxide at 10% strain along x. The friction ratio between two directions increases from 1 to 1.6 for water and from 1 to 1.3 for carbon dioxide. However, for spherical molecules like methane, the variation of friction is insignificant. Investigation on the structure and dynamics of those liquids has revealed that the relaxation time increases considerably with strain and contributes an important part in the increase of the friction.

The contribution of this work helps better understanding and modelling the friction between liquids and anisotropic surfaces. Those surfaces can exist naturally, for example in the form of orthorombic crystal systems or cubic systems subject to misfit strain etc.. The study also contributes practical aspects to answer energetic and environmental challenges, specifically it is closely related to carbon dioxide sequestration process via Enhanced coal bed methane recovery method.

FIG. 2 .

 2 FIG. 2. Comparison of methane-graphene potentials V (x, z O ) at the first fluid layer, about one molecule diameter from graphene, produced by exact calculation and two analytical approximations (11) and (12) (a) Exact distribution of potential energy V (x, z O ) in the first layer. (b) Distribution of potential energy calculated by expression[START_REF] Bocquet | [END_REF]. (c) Distribution of potential energy calculated by expression[START_REF] Evans | Statistical mechanics of nonequilibrium liquids[END_REF]. It is easy to notice that (11) can recover the symmetry of the exact potential while[START_REF] Evans | Statistical mechanics of nonequilibrium liquids[END_REF] can not.

FIG. 3 .

 3 FIG. 3. Snapshots of the three systems (H 2 O, CO 2 and CH 4 ) in graphene nano-channel under consideration. The shape of three molecules are drawn according to molecules' bond lengths, bond angles and σ f w diameters.

FIG. 4 .

 4 FIG. 4. Fluid molecular density for different systems: (a) for H 2 O, (b) for CO 2 and (c) for CH 4 . (d) variation of liquid density around the first peaks (first liquid layer).

FIG. 5 .

 5 FIG. 5. Time correlation of friction forces along x and y directions for H 2 O and graphene.

2 OFIG. 7 .

 27 FIG. 7.Influence of the fluid models on the slip length L s of water-graphene system.

FIG. 8 .

 8 FIG. 8. (a) 2D structure factor of oxygen (H 2 O) with 10% unilateral strain on graphene sheet. (b)The structure factor of oxygen and hydrogen of water with different q. (c) The structure factor of oxygen (H 2 O), carbon (CO 2 ) and methane (CH 4 ) at strain-free state.

FIG. 9 .

 9 FIG. 9. Coefficient α (two upper figures) and potential energies V 1 (two lower figures) as a function of the strain and the distance z of methane molecule from the graphene sheet.

FIG. 10 .

 10 FIG.10. Factorization of friction coefficient λ to F 2 /A and τ t with strain on graphite ranging from 0% to 10%

TABLE I .

 I Interaction LJ parameters for H 2 O, CO 2 and CH 4 interacting with graphene.For CO 2 -CO 2 interaction, we have used EPM2 model[29] which consists of 12-6 LJ sites in conjunction with partial charges centered on each of sites. The O-C-O bond angle 180 • and the C-O bond length 1.16 Å are maintained constant by rigid body algorithm. For CO 2 -Graphene interaction, we have used LJ potential computed using Lorentz-Berthelot rules: σ AB = (σ A + σ B )/2 and AB = √ A B with LJ parameters of graphene from Ref.[30]: σ G = 3.36 Å, G = 2.413 meV (see TableI).

  Reactive force field[28]. b Best estimate[32].

				eq [ Å]	Ref.
	H 2 O-Graphene	130.37	139.92 a	3.28	3.19 a
	CO 2 -Graphene	158.70		3.13
	CO 2 -Graphite	177.21	178.2 b	3.11	3.2±0.1 b
	CH 4 -Graphene	102.68		3.49
	CH 4 -Graphite	118.59	130±0.1 b	3.47	3.45 b
	a

TABLE II .

 II Potential well depths and equilibrium geometries for the global interaction between one molecule of H 2 O, CO 2 , or CH 4 interacting with graphene or graphite.

TABLE III .

 III Simulation conditions for different liquid systems including N number of liquid molecules, ρ the average liquid density confined between the two graphene sheets, T the temperature and P ext the external pressure.

TABLE V .

 V The norm of reciprocal vectors for different strains. Note that q 1

ACKNOWLEDGMENTS

The authors Quy-Dong To and Van-Hoang Vo would like to thank program Hoa Sen Lotus 30569SK and the organization CampusFrance for the travel support. We also thank the Multi-Scale Modelling & Experimentation of Materials for Sustainable Construction (LabEx MMCD) for their support in CPU-time. We thank Dr. Laurent Joly for his helpful discussion with us.

[1] R. B. Schoch, J. Han, and P. Renaud, Rev. Mod. Phys. 80, 839 (2008).