Influence of Confinement on Flow and Lubrication Properties of a Salt Model Ionic Liquid Investigated with Molecular Dynamics
IInfluence of Confinement on Flow and Lubrication Properties of a Salt Model IonicLiquid Investigated with Molecular Dynamics
Miljan Daˇsi´c and Igor Stankovi´c ∗ Scientific Computing Laboratory, Center for the Study of Complex Systems,Institute of Physics Belgrade, University of Belgrade, Pregrevica 118, 11080 Belgrade, Serbia
Konstantinos Gkagkas
Advanced Technology Division, Toyota Motor Europe NV/SA,Technical Center, Hoge Wei 33B, 1930 Zaventem, Belgium (Dated: October 2, 2018)We present a molecular dynamics study of the effects of confinement on the lubrication and flowproperties of ionic liquids. We use a coarse grained salt model description of ionic liquid as a lubricantconfined between finite solid plates and subjected to two dynamic regimes: shear and cyclic loading.The impact of confinement on the ion arrangement and mechanical response of the system has beenstudied in detail and compared to static and bulk properties. The results have revealed that thewall slip has a profound influence on the force built–up as a response to mechanical deformationand that at the same time in the dynamic regime interaction with the walls represents a principaldriving force governing the behaviour of ionic liquid in the gap. We also observe a transition froma dense liquid to an ordered and potentially solidified state of the ionic liquid taking place undervariable normal loads and under shear.
PACS numbers: 68.35.Af, 47.85.mf,47.27.N,83.10.-y,83.50.Jf
INTRODUCTION
In this work, the lubricating ability and flow proper-ties of salt model ionic liquids (ILs) containing salt–likespherical cations and anions are studied. ILs are moltensalts typically consisting of large–size organic cationsand anions. The thermochemical stability, negligible va-por pressure, viscosity, wetting performance and otherphysicochemical properties of ILs are important factorscontributing to the interest in their research for lubri-cant applications [1, 2]. In addition, their propertiescan be modified by an applied voltage using confiningcharged surfaces in order to build–up an electric fieldacross the nanoscale film. The applied potential canaffect the structure of IL layers and lead to externallycontrollable lubricating properties [3–5]. There is also asignificant flexibility in tuning the physical and chemi-cal properties of ILs which can affect lubrication such asviscosity, polarity and surface reactivity, by varying theiratomic composition as well as the cation–anion combina-tion. The thermal stability and negligible vapor pressureof ILs enable their usage at a high temperature.Regarding the ability of ionic liquids to dynamicallypenetrate between surfaces, i.e. wetting, sometimes it isconsidered that a low contact angle of the lubricant indi-cates the affinity between the liquid and the surface, sincethe liquid is more likely to stay in the area in which itwas initially placed. It is also expected that a lubricant isgoing to penetrate into small–gap components. However,the effect of wettability of the ionic liquids is not under-stood well. The wetting of plate surfaces such as mica isknown to be partial by at least some ILs [6, 7]. Lubri-cation necessarily involves intimate molecular features of the liquid–solid plate interface, related with those mecha-nisms determining the ionic liquid’s wetting of the plate.When ILs are used as lubricants, their ions are orderedinto layers and adsorbed onto surfaces [8]. These adsorp-tion layers can reduce friction and wear, particularly inthe case of boundary lubrication [8].An important observation is that ILs confined be-tween surfaces feature alternating positive and negativeionic layers, with an interlayer separation correspond-ing to the ion pair size [9, 10]. However, determiningthe structure of ILs during flow and the mechanism ofnanoscopic friction with ILs as lubricants, poses a greatscientific challenge, and so far a few studies in this di-rection have been performed [3]. ILs involve long–rangeCoulombic interactions inducing long–range order on fargreater scales than the IL itself [11–13]. Recent studies ofIL lubricants [3–5, 10] have shown that if the moleculesinteract via non–bonded potentials (Lennard–Jones andCoulombic potential), this can capture all main physicalattributes of the IL–lubricated nanotribological system.Therefore, molecular–scale simulations can provide im-portant insights which are necessary for understandingthe differences in flow behaviour between bulk and con-fined liquid lubricants and the mechanisms behind, suchas boundary layers formation in case of shearing and/orapplied normal load.For this study, we utilize our previously developedcoarse grain molecular dynamics (MD) simulation setupconsisting of two solid plates, and an ionic liquid lubri-cant placed between them [10]. The motivation for thechosen values of relevant model parameters (i.e. veloc-ities, pressures, temperatures, time duration of simula-tions) comes from potential applications of ILs as lubri-Typeset by REVTEX a r X i v : . [ c ond - m a t . s o f t ] S e p cants in automotive industry. Under typical operationof internal combustion engines, the conditions inside thecombustion chamber vary significantly. Temperature canrange from 300 K to the values higher than 2000 K, whilepressure ranges from atmospheric to the values higherthan 10 MPa [14]. The piston reciprocates with a sinu-soidal velocity variation with speeds varying from zeroto over 20 m/s, with a typical speed being around 1m/s.The time required for one revolution of the engine is ofthe order of 10 − s, while the total distance travelled bythe piston over this period is of the order of 0 . MODEL
The model used in this work is a coarse–grained modelof IL which has already been exploited in previous stud-ies [3–5, 10] and it is known as SM model (salt–likemodel). It is a charged Lennard–Jones system consistingof cations and anions. There are two types of interatomicinteractions in our system and both of them are non–bonded: Lennard–Jones (LJ) potential and Coulombicelectrostatic potential. In the current work we are com-paring bulk and confined IL properties. Therefore, thereare three different atom types taken into consideration:( i ) cations, ( ii ) anions and ( iii ) solid plate atoms. Be-tween all types of atoms we apply full LJ 12-6 potential,with the addition of Coulombic electrostatic potential forthe interactions between ions. In our system the cationsand the anions are charged particles, while the solid plate atoms are electroneutral. Accordingly, we have imple-mented a LJ 12-6 potential combined with Coulombicelectrostatic potential: V ( r ij ) = 4 (cid:15) ij (cid:34)(cid:18) σ ij r ij (cid:19) − (cid:18) σ ij r ij (cid:19) (cid:35) + 14 π(cid:15) (cid:15) r q i q j r ij (1)Parameters { (cid:15) ij , σ ij } define the LJ potential betweendifferent types of particles: i, j = A , C , P which refer toanions, cations and solid plate atoms, respectively. Thediameter of cations and anions is set to σ CC = 5 ˚A and σ AA = 10 ˚A, respectively. The mass of cations and anionsis m C = 130 g/mol and m A = 290 g/mol, respectively.The asymmetry of ion sizes is typical in many experimen-tally explored systems and the parameters have alreadybeen explored in literature, cf. Ref. [5, 10]. The atoms ofthe solid plates have a diameter of σ PP = 3 ˚A. The massof the solid plate atoms is m P = 65 g/mol. The LJ po-tential has a short–range impact, since it vanishes rapidlyas the distance increases ∝ r − , while the Coulombic po-tential has a long–range impact, ∝ /r . To handle long–range interactions, we have used a multi–level summationmethod (MSM) [15], since it scales well with the numberof ions and allows the use of mixed periodic (in x and y di-rections) and non-periodic (in z direction) boundary con-ditions, which are present in our simulation setup withconfined IL. On the other hand, in our simulation setupwith bulk IL, periodic boundary conditions are appliedin all three directions ( { x, y, z } ). Ions are modelled ascoarse grain particles, the charge of which is set equal toelementary: q C = + e and q A = − e , i.e., e = 1.6 · -19 C.The dielectric constant is set to (cid:15) r = 2 to account for thedielectric screening, as in Refs. [4, 5, 10].In this study, modelling the elasticity of metallic platesplays a secondary role (central role belongs to the lat-eral and normal forces created by the lubricant). There-fore, we have selected a simplified model in which plateatoms interact strongly with each other if they belongto the same plate, i.e., (cid:15) PP = 120 kCal/mol, as op-posed, to a very weak interaction between the differentplates (cid:15) top / bottom = 0 .
03 kCal/mol. The parameter (cid:15) PP is so strong in order to ensure that the initial configura-tion of the solid bodies will basically remain unchanged(apart from high frequency oscillations). Furthermore,even though typical engineering systems are often metal-lic, our initial coarse grained MD studies of liquid be-haviour according to the applied conditions justified theimplementation of a simpler solid system which does notfeature substrate polarization, cf. Ref. [10]. Finally, itis possible that the actual surfaces might feature carboncoatings or depositions, in which case the surface polar-ization can be of secondary importance.In Table I we are presenting the values of { (cid:15) ij , σ ij } parameters used in our model. Arithmetic mixing rulesfor the LJ parameters are applied: (cid:15) ij = √ (cid:15) i · (cid:15) j , σ ij = pair ij (cid:15) ij [kCal/mol] σ ij [˚A]CC 0.03 5AA 0.03 10CA 0.03 7.5PC 0.3 4PA 0.3 6.5PP 120 3TABLE I: List of LJ parameters used in simulations. ( σ i + σ j ) / BULK IONIC LIQUID
All MD simulations in this study were performed usingthe LAMMPS software [16]. The bulk ionic liquid is im-plemented by randomly placing a chosen number of ions( N C = N A = 1000) into a 3D simulation box that is peri-odic in all three directions. In order to make the bulk ILcomparable with its confined counterpart, we have cho-sen a simulation box volume which enables the pressureexperienced by the confined IL. More specifically, for thepresent system, the pressure is p ≈ T = 330 K. The system was relaxedfor t tot = 3 · fs until the internal energy had convergedand the pressure had approached the desired value. Thesimulation timestep was dt = 0 . (cid:104) p (cid:105) = 1 . L = 99 ˚A. The energy re-laxed to a value of (cid:104) E int (cid:105) = 0 . ρ n = 3400 mol / m and ρ m = 719 kg / m respectively. FIG. 1: Configuration snapshot ( yz cross–section) of a bulk ILat the end of relaxation simulation. Cations are representedas smaller blue spheres and anions as larger red spheres. We have calculated the viscosity in two ways: using theGreen–Kubo relation since the viscosity of a system can
FIG. 2: Dependence of Green–Kubo (GK) viscosity coefficient η GK on simulation time t s in case of bulk ionic liquid. Thetime needed to obtain the viscosity coefficient is around t rel =5 ns.FIG. 3: Average stress tensor component τ in function of theshear rate ˙ γ of a bulk SM ionic liquid. We have conductedshearing simulations on four orders of magnitude of the shearrate ˙ γ , therefore with three orders of magnitude span, whichis followed by three orders of magnitude span of τ . Points areobtained via shearing simulations and solid line is obtainedaccording to: τ = η GK · ˙ γ , where η GK is obtained via Green–Kubo relation. be represented as an integral of the autocorrelation func-tion [17], and using non–equilibrium molecular dynamicssimulations with different shear strains.In the non–equlibrium shearing simulations, the bulkIL is placed into a triclinic (non–orthogonal) simula-tion box with periodic boundary conditions applied inall three directions. Due to the deformation of the sim-ulation box, every point in the box has an additionalvelocity component (a so called streaming velocity ). Inorder to prevent the streaming velocity from affecting thethermal kinetic energy, we use the so-called SLLOD ther-mostat [18, 19]. The
SLLOD thermostat accounts for thestreaming velocity which depends on an atom’s positionwithin the simulation box and it needs to be accountedfor controlling the temperature.
FIG. 4: (a) Schematic of the simulation setup shown as yz cross–section. There are two solid plates at the top and bottomof the system (i.e. Top and Bottom plate, names are chosen according to their position along the z axis), consisting of tworegions: at the outermost ones the particles are moving as a single entity (Top/Bottom Action) and at the innermost ones theparticles are thermalized at a controlled temperature (Top/Bottom Thermo). The thermalized layers are in direct contact withthe lubricant while the action layers are used to impose external velocity and/or force to the solid plates. (b) Front ( yz ) viewwith respect to the shear direction. A typical simulation configuration and key dimensions of the geometry are given. Thesolid plates are made up of FCC (111) atomic layers. The ionic liquid is composed of an equal number of cations and anions(cations: smaller blue spheres; anions: larger red spheres). Controlled shearing of the simulation box resultsin a stress acting on IL, which is quantified via thestress tensor. The relation between the stress tensor τ ij components and the shear rate ˙ γ ij of correspond-ing shear strain (cid:15) ij , with coefficient of viscosity η ij asa proportionality constant is: τ ij = η ij · ˙ γ ij where ij = { xy, xz, yz } . We have applied three indepen-dent shear strains ( (cid:15) xy , (cid:15) xz , (cid:15) yz ). For each of them wehave calculated its corresponding stress tensor compo-nent ( τ xy , τ xz , τ yz ). All shear strains were the same: (cid:15) xy = (cid:15) xz = (cid:15) yz = (cid:15) = 1 leading to the shear rate of˙ γ = (cid:15) · t tot = t tot , where t tot is the total simulation timeof the shearing simulations. We have performed simula-tions at four orders of magnitude of the total simulationtime: t tot = { . , , , } ns, and thus at four orders ofmagnitude of the corresponding shear rate. In this waywe wanted to check the quality of our relaxation proce-dure and if there are shear rate dependence changes in thesystem. We have iterated the shearing simulations (at ashearing velocity of 1 m/s) using the output of the pre-vious run as the input of the next run, obtaining higherstrains (up to a strain of 5). We did not observe a straindependence in the response of the system, meaning thatthe result is unaffected if the strain is further increased.In Figure 2, we show the time relaxation of the Green–Kubo viscosity coefficient, which stabilizes around η GK =0 . · s. The configuration snapshot of the bulkIL at the end of the simulation (cf. Figure 1) shows thatthe ions remain randomly positioned, like they were atthe start of simulation, which indicates the liquid stateof the bulk ionic liquid. The simulations for all threeshear strains give similar values of stress components,and resulting values are shown in Figure 3. The points { ˙ γ, τ } were obtained via shearing simulations and thesolid line was obtained according to τ = η GK · ˙ γ , where η GK was obtained via Green–Kubo relation. Hence, weconclude that the results of shearing simulations are inagreement with the results of relaxation simulation andtherefore there are no changes in the bulk system whichare shear rate dependent. CONFINED IONIC LIQUID
In order to study the properties of our ionic liquid un-der confinement, we use a setup consisting of two solidplates (so called Top and Bottom plate) and ionic liquidlubricant placed between them. Such simulation setuphas been introduced and described in detail in our pre-vious paper [10], hence at this point we will describe itbriefly. The geometry is shown as a schematic in Fig-ure 4(a) together with the number of the coarse grainparticles used. In Figure 4(b) we show a configurationsnapshot of our system in yz cross–section. By imple-menting such a geometry we have attempted to achievea realistic particle squeeze–out behaviour with the forma-tion of two lateral lubricant regions in a similar mannerto the simulations of Capozza et al. [5]. For the descrip-tion of the solid surfaces we have combined rigid layers ofparticles moving as a single entity on which the externalforce or motion is imposed, denoted by ”Top Action” and”Bottom Action” in Figure 4(a), with thermalized lay-ers, denoted by ”Top Thermo” and ”Bottom Thermo” inwhich particles vibrate thermally at T = 330 K. The par-ticles in the Top and Bottom action layers are moved asrigid bodies and particles in the thermo layers are allowedto move thermally. In this way, we prevent a progressivedeformation of the plates during the cyclic movement.The thermo layers only vibrate thermally since a strongLJ interaction holds them together. The ionic liquid isneutral in total, so the total number of cations and an-ions is the same: N C = N A = N IL /
2. In the presentsimulations the total number of IL atoms is N IL = 2500.The plates are driven along the x direction at a con-stant velocity V x , as shown in Figure 4(a). The solidplates are made up of nine atomic layers at a FCC (111)lattice arrangement. Periodic boundary conditions areapplied in the x and y directions, while the simulationbox is kept fixed in the z direction. The Bottom platecan therefore be considered to be infinite, while the Topplate is surrounded by the lateral reservoirs, in which thelubricant can freely expand. The lateral reservoirs areimplemented as a mechanistic way for allowing the lubri-cant to be dynamically squeezed in or out as an exter-nal load or velocity is applied. The number of lubricantmolecules effectively confined inside the gap can dynam-ically change depending on the loading conditions. Thisis important for exploring the possible states of a me-chanical system of particles that is being maintained inthermodynamic equilibrium (thermal and chemical) witha lubricant reservoir (i.e., void spaces in tribological sys-tem). The nano–tribological system is open in the sensethat it can exchange energy and particles, realizing aneffectively grand–canonical situation [20, 21].We have shown that our bulk IL is a Newtonian fluid:the validity of τ = η GK · ˙ γ relation over the whole range ofshearing rate ˙ γ supports that fact. Our model does notassume the nature of viscous response of IL. Only basedon simulation results we conclude that bulk salt model(SM) IL behaves as a Newtonian fluid. For a differentchoice of parameters one could obtain power law or solidlike behaviour. On the other hand, confinement has aprofound influence on the structure of ILs in thin films [9,10, 20, 22], therefore when the same IL is confined it doesnot behave as a Newtonian fluid, as we will show in therest of the paper.The confining surfaces can induce ordering of the parti-cles in their vicinity. We have used simulations to obtainthe static force–distance characteristic [10]. In order todetermine a reliable static force–distance characteristic,at each calculation point we have to ensure the system isin equilibrium. Concerning the realization of those sim-ulations the inter–plate gap is modified in the followingmanner: the gap is increased or decreased (i.e., the Top–Bottom plate distance is changed) with a constant veloc-ity V z = 5 m/s for a move period of time t move = 20 ps;thereafter, we apply conjugated gradient minimizationon the ions in order to minimize their internal energyand relax them after the move period. As the energyminimization is performed, the ions take positions whichensure their minimal internal energy and the Top platestays fixed for a stay period of time t stay = 50 ps, dur- ing which period the average value of the normal force iscalculated; that value is presented as a simulation pointin F z ( d z ) static characteristic, cf. Figure 5. In orderto avoid systematic errors due to the initial position ordirection, the plate movement is performed in differentdirections and from different initial configurations, hencethe Figure 5 shows the averaged values of several realiza-tions. FIG. 5: Dependence of normal force F z acting on the Topplate on plate-to-plate distance d z . Five characteristic points { I , I , I , II , II } with corresponding interplate distances d z ≈ { , , , , } ˚A are marked on the F z ( d z ) curve.Also, the two characteristic intervals of d z are labeled, wherethe interval I is bounded by the points I and I , while theinterval II is bounded by the points II and II . The hori-zontal solid line denotes F z = 0 pN. The dashed line connectsthe points obtained from the simulation and serves as a visualguide. Equilibrium Behaviour of Confined Ionic Liquid
A non–monotonous behaviour of the normal force F z acting on the Top plate can be observed in Figure 5 asthe plate-to-plate distance is changing from one point atwhich system is equilibrated to another, using the previ-ously described procedure. The points ( d z , F z ) have beenobtained through our simulations, while the dashed lineserves as a visual guide. It can be seen that the normalforce strongly depends on the inter–plate distance. Thepresence of negative values of normal force F z can be un-derstood as the IL trying to reduce the plate-to-platedistance due to adhesion phenomena. These changesof the normal force are correlated with the extractionand inclusion of IL layers into the gap, as already ob-served experimentally, cf. Ref. [9]. During the performedstationary state simulations, the cationic–anionic layerswere squeezed out in pairs, in order to keep the sys- FIG. 6: Configuration snapshots ( yz cross–section) accompanied with ionic density distribution along the z direction in threerepresentative points of the interval I : { I , I , I } . Left panels correspond to the static case of Top plate’s movement, whileright panels correspond to the dynamic case of Top plate’s movement. tem locally neutral, as observed in experimental stud-ies [8, 9, 20, 22, 23].There is a strongly decreasing trend of the maximalnormal force which can be sustained by the system asthe number of ionic layers confined between the platesincreases, i.e., for the two ionic layers the maximal force F I z,max = 3 pN, while for the three ionic layers it is F II z,max = 0 .
25 pN. In our model, the Lennard-Jones in-teraction between the plates and the ions is ten timesstronger than between the ions themselves. The ioniclayers closest to the plates are therefore more stable thanthe layers in the midpoint of the gap (interval II ). As aresult, the three-layer system becomes less dense and canbuild up a lower normal force compared to the two-layersystem. We have selected two intervals of interest forthe interplate distance which capture the presence of lo-cal maxima and subsequent minima of the normal force F z accompanied with the compression of IL. This cor-responds to the expulsion of a cation–anion layer pairfrom the gap. The intervals are: d Iz = [14 . ,
20] ˚A, d IIz = [22 ,
27] ˚A, and they are labeled as I and II respec- tively. In order to understand the changes of the systemconfigurations and to correlate them with the changes ofthe inter–plate distance, snapshots of the system fromthe MD simulations corresponding to several character-istic points of the intervals I and II have been selectedand studied in more detail: I , , II , which correspondto the limits of the intervals, and the local maximum ofthe interval I , labeled as I .The left vertical panel of Figure 6 shows the systemconfiguration in the yz cross-section and the ionic densitydistribution along the z–direction obtained in the equilib-rium force–distance simulations for the three character-istic points of the interval I : { I , I , I } , correspond-ing to plate-to-plate distances d z = { . , , . } ˚A,respectively. In Figure 7 the left vertical panels showanalogous results for the two characteristic points of theinterval II : { II , II } , corresponding to plate-to-platedistances d z = { , } ˚A, respectively. The ions are de-picted smaller than their LJ radii in order to allow adirect observation of the layering. The positions of theatomic centres of the innermost atomic layers of the Topand Bottom plate are labeled in Figures 6 and 7 as z T FIG. 7: Configuration snapshots ( yz cross–section) accompanied with ionic density distribution along the z direction in tworepresentative points of the interval II : { II , II } . Left panels correspond to the static case of Top plate’s movement, whileright panels correspond to the dynamic case of Top plate’s movement. and z B respectively. As the Bottom plate is kept fixedduring the whole simulation, z B remains constant while z T changes with the Top plate displacement. A generalfeature observed under all conditions was the formationof cationic layer close to the plates. The reason is thesmaller size of the cations ( σ CC = 5˚A) compared to theanions ( σ AA = 10˚A). Following, the second layer getsinduced by the first one (due to Coulombic interaction)and it is populated by anions. The distance between thefirst and the second layer from the bottom is in the rangeof 1 − . / √ yz cross–sectionconfiguration snapshots together with the ionic densitydistribution along the z axis, shown in the left panelsof Figures 6 and 7 for the cases of intervals I and II ,respectively, we have prepared the xy cross–section con-figuration snapshots, shown in the left panels of Figures 8and 9. Cyclic Compression of Confined Ionic Liquid
We have investigated the dynamic behaviour of theIL during a periodic linear movement of the Top plate along the z axis, between the two limiting points of theintervals I and II . The space between the solid plateswas in this way periodically expanded and compressed.Periodic movements of the Top plate were performed atthree constant velocities V z = { . , , } m/s but novelocity dependent differences in the system behaviourwere observed. We have performed ten cycles in orderto determine how much the cycles differ and to deter-mine a statistically reliable average cycle. The confinedionic liquid lubricant responds to the cyclic movementwith a hysteresis in the normal force F z ( d z ) as shownin Figure 10. We present both raw data of all cycles(thin solid lines) and a smooth average cycle (thick solidline). In the case of interval I there are three points ofinterest { I , I , I } , corresponding to the points noted inFigure 5. Points I and I are the starting and endingpoint respectively and the point I corresponds to themaximum of the normal force F z in the smooth averagecycle. We observe that between each two of those pointsthere are clear tendencies in the average cycle of the nor-mal force F z ( d z ). First, in the segment I → I , i.e. inthe extension half of the cycle, there is a continuous in-crease of normal force F z from negative values up to thevalue around zero in point I . In point I there is oneanionic layer confined in the gap and normal force F z hasa negative value. With the dynamic increase of the gapions are pulled–in from lateral reservoirs into the gap. Inpoint I an additional cationic–anionic layer pair is fullyformed in the gap, hence increasing the number of con-fined anionic layers to two. Next, there is the segment FIG. 8: Configuration snapshots ( xy cross–section) in threerepresentative points of the interval I : { I , I , I } . Left pan-els correspond to the static case of Top plate’s movement,while right panels correspond to the dynamic case of Topplate’s movement. We have highlighted the confined regionwith dashed lines (Top plate’s width along the y axis is ahalf of the total system’s width) and also we have sketchedcrystallization patterns with solid lines. Periodic boundaryconditions are applied in the x and y directions, while simu-lation box, which is cubic, is kept fixed in the z direction. I → I where the ions are compressed within the gap,which is consistent with the continuous increase of nor-mal force F z . In this segment, the normal force F z takespositive values meaning that the ionic liquid shows resis-tance to the compression but does not flow out. Afterthat, in segment I → I there is a sharp decrease of nor-mal force F z which is correlated with the squeezing–outof the additional cationic-anionic layer taken in from thelateral reservoirs during the extension half–cycle. Duringthe compression half–cycle there is a return to the initialstate I , where the gap contains one compact anioniclayer.We should note that the distributions of cations andanions in the dynamic case for interval I bear close re-semblance. Let us now discuss the changes in the num-ber of confined ionic layers as a function of the inter–plate distance and correlate them with the changes inthe normal force F z acting on the Top plate: in therange d z = [11 , .
2] ˚A the normal force F z acting onthe Top plate has a steep decrease, reaching the mini-mum at point I . For the point I at d z = 14 . FIG. 9: Configuration snapshots ( xy cross–section) in tworepresentative points of the interval II : { II , II } . Left pan-els correspond to the static case of Top plate’s movement,while right panels correspond to the dynamic case of Topplate’s movement. We have highlighted the confined regionwith dashed lines (Top plate’s width along the y axis is ahalf of the total system’s width) and also we have sketchedcrystallization patterns with solid lines. plete cationic layers inside the gap. The value of normalforce F z is negative and in point I it has the deepestminimum when considering the whole F z ( d z ) character-istic. With increasing plate-to-plate distance d z the nor-mal force F z is increasing, with a sign change of normalforce F z around d z = 15 . d z = 17 . F z > F z reaches a local maximum in the point I at d z = 17 . d z beyond the point I there is a continuous de-crease of the normal force up to the distance d z = 14 . i ) the sign of the normal force in point I and ( ii ) themagnitude of the normal force at local maximum I . Inthe case of cyclic (dynamic) movement of the plates, thenormal force is positive F z >
0, i.e. the IL keeps pullingapart the plates at point I and the maximum of thenormal force in the point I ( F dynz = 1 pN) is lower thanin the equilibrium case ( F eqlz = 3 pN). Both observationsindicate that the plate’s motion is preventing the ionicliquid to fully fill the void space of the gap. Also, thereis a substantial slip during the ejection of IL from thegap, which results in a lower normal force. Otherwise,if no slip would be present the maximum normal force FIG. 10: This figure presents the results of dynamicextension–compression cycles in the intervals I and II . Inpanel (a) we present dynamic F z ( d z ) characteristic in theinterval I : thin lines represent the hystereses of ten dynamiccycles, solid line on top of them is the smooth average hys-teresis. There is also a solid horizontal line which correspondsto F z = 0. Panel (b) is analogous to the panel (a), just itpresents the results in the interval II . at velocity V z = 1 m/s should be about two orders ofmagnitude higher based on the bulk viscosity coefficientcalculated in the previous section.Partial filling of the gap due to the motion of the wallsis even better observable in the results for the interval II . While the equilibrium characteristic has a local max-imum, cf. Figure 5, in the dynamic case there are onlytwo characteristic points (starting and ending point { II , II } and monotonously increasing normal force betweenthem. At point II at d z = 22 ˚A in the static case, wenotice that at the mid point between the plates there isa broad maximum of the cation density distribution, seeFigure 7. In the static case we notice that, similar tothe transition from one to two anionic layers within theinterval I , there is a transition from two to three anioniclayers within the interval II , which happens in proximityof the point d z = 24 ˚A. At point II we notice two sharpanionic layers in the proximity of the plates and the thirdanionic layer which is broader, less sharp and positionedin the middle of the inter–plate gap, cf. Figure 7. In thedynamic case the number of layers remains the same in the interval II , they just get separated during the exten-sion; and a creation of additional ionic layers by the ionsflowing from the lateral reservoirs into the gap does nottake place, cf. Figure 7.We can conclude that in a confined system with stronginteraction between the walls and the IL, the major driv-ing force that pulls IL into the gap between the plates isthe interaction with the wall atoms rather than the inter–IL interactions. In order to visualize what happens inthe vicinity of the plates, we are presenting snapshots of xy cross–section configurations in the intervals I and II ,check the Figures 8 and 9, respectively. Even on a cursorylook, one sees that the phase behaviour of the confined ILis complex: in Figure 8, we observe a salt–like orderingtaking place in all representative points { I , I , I } of theequilibrium configurations. In the dynamic case the ILexhibits some level of ordering for a small gap ( I ) and itis amorphous in the other two points. On the other hand,in Figure 9 there was no movement of the IL in and outof the gap and the IL formed a two–dimensional squarecrystal { II , II } on both surfaces during the dynamiccase. In the equilibrium configurations, there are prob-ably enough ions in the gap that allow the IL to obtainits liquid–like character.At this point, we would like to quantify how could theprocesses described above contribute to the energy losses.If two macroscopically smooth surfaces come into con-tact, initially they only touch at a few of these asperitypoints. A motion of two bodies in contact lubricated byan ionic liquid would involve the generation of new con-tacts and the separation of existing ones. Ionic liquids arecharacterized by strong Columbic interactions betweenthe particles. By calculating the area covered within theaverage cycle of the F z ( d z ) curves in Figure 10, we cal-culate the amount of work invested per average dynamiccycle, i.e., the hysteretic energy losses. There is a bigdifference in the amount of invested work in the two in-tervals: 3 . · ˚A for the interval I compared to0 . · ˚A for the interval II , where the vertical dis-placement of the Top plate in the two intervals is roughlythe same ∆ d z ≈ F I z,max = 3 pN,while for the three ionic layers it is F II z,max = 0 .
25 pN, cor-responding to the two maxima of the equilibrium force–distance characteristic in Figure 5.
Shear Behaviour of Confined Ionic Liquid
In order to study the behaviour of our confined ILunder shearing, we apply a relative motion between theplates along the x direction. The Bottom plate is keptfixed and a constant velocity V x is imposed on the Top0plate. We are interested in establishing how does the lat-eral (frictional) force F x depend on the confinement gap d z = { , , , , , } ˚A. FIG. 11: Dependence of the frictional force divided by thecontact area of the Top plate with IL lubricant (cid:104) F x (cid:105) /S xy onthe interplate distance d z . The three representative points { P , P , P } are marked. Points obtained in simulations areshown as circle markers, accompanied with errors along the y axis. Linear fit through those points is shown as a solidline. In the inset dependence of specific friction (cid:104) F x (cid:105) / (cid:104) F z (cid:105) on the interplate distance d z is shown, with y axis in logscale. Simulation points are shown as circle markers, whilethe dashed line serves as a visual guide.FIG. 12: Dependence of the frictional force divided by thecontact area of the Top plate with IL lubricant (cid:104) F x (cid:105) /S xy onthe Top plate’s lateral velocity V x = 0 . −
10 m/s. The errorbars represent the standard deviation of the average valuesobtained from the simulation data. The lines showing thefriction trends are obtained by linear regression.
In Figure 11 we are showing the dependence of the timeaveraged frictional force divided by the contact area ofthe Top plate and the IL lubricant, i.e. (cid:104) F x (cid:105) /S xy on theinterplate distance d z . We observe a linear increase of the frictional force per contact area with the increase ofthe interplate distance, with a slope of 4 nN /µ m . Inthe inset of Figure 11, we are showing the dependenceof specific friction defined as the ratio of the time av-eraged frictional and normal force (cid:104) F x (cid:105) / (cid:104) F z (cid:105) on the in-terplate distance d z . By comparing the Figure 11 withthe results for the bulk liquid in Figure 3 we observethat there is no correlation with the lubricant viscos-ity (i.e., otherwise frictional force would be three ordersof magnitude higher). This leads us to the assumptionthat our pressurized systems, whether they form a crys-talline lattice or not, do not lie in a typical hydrody-namic regime and operate under full slip conditions inwhich the ionic liquid moves together with one of thewalls. As there is no solid–solid contact between thetwo surfaces, but lubrication through very thin, highlyviscous films which are solid–like, mixed or dry lubri-cation are the two potential regimes that can describethe observed conditions. A parametric study on differ-ent shearing velocities V x = 0 . −
10 m/s at two wallseparations d z = 17 ,
27 ˚A provides additional informa-tion for the characterization of the tribological regime ofour system. In Figure 12 one can observe a logarithmic(weak) dependence of the frictional force per contact areaon lateral velocity of the Top plate’s movement which isconsistent with the observations of previous studies of ILlubrication, cf. Refs. [11, 13].From Figure 11 we have selected three representativepoints with d z = { , , } ˚A labeled as { P , P , P } respectively. We provide an overview of the yz config-uration cross–sections together with ionic density distri-butions along the z axis (cf. Figure 13) at the simulationonset t = 0 and after t = 3 ns. In the panels of Figure 14we have highlighted the confined region with dashed lines(the Top plate’s width along the y axis is half of the totalsystem’s width) and we have also sketched crystallizationpatterns with solid lines. In Figures 13 and 14 we showinitial configurations at the input of friction simulations,together with the final configurations obtained after thefriction simulations. We observe that any initial crys-tallization is not lost due to the lateral motion of theTop plate, but only slightly modified due to the motion,which suggests that the lateral movement does not al-ter the ordering. This is a significant finding since thelongitudinal movement does alter the local order (it de-stroys the crystal structure for small gaps and induces itin larger ones). CONCLUSIONS
In the current work we have used a molecular dynam-ics simulation setup in order to study the response ofa model ionic liquid to imposed mechanical deformation.The properties of bulk and confined ionic liquid have beeninvestigated under both static and dynamic conditions.1
FIG. 13: Configuration snapshots ( yz cross–section) accompanied with ionic density distribution along the z direction in threerepresentative points { P , P , P } . Left panels correspond to the start of friction simulations t = 0, while right panels correspondto the end of friction simulations t = 3 ns. Top plate’s lateral velocity is set to V x = 2 m/s, total simulation time is t tot = 3 ns,hence all friction simulations have run until the Top plate had covered a distance of d x = V x · t tot = 60 ˚A along the x direction. First, we have shown that the Green–Kubo viscosity coef-ficient fits the shearing simulation results of our bulk saltmodel ionic liquid, indicating its liquid state. Our simula-tion results have shown the significant impact of the con-finement and interaction with the walls on the ionic liquidresponse to mechanical deformation. The force–distancehysteresis surface under cyclic loading is smaller than onewould expect considering only the viscosity value of theliquid. The simulations have also shown the transitionfrom a liquid to a highly dense and ordered, potentiallysolidified state of the IL taking place under variable nor-mal load and under shear. The wall slip has a profoundinfluence on all the forces which arise as a response tothe mechanical deformation. We also observe that theinteraction of the IL with the walls represents a principaldriving force for all processes observed in the dynamicregime for a range of studied velocities. If sufficient timeis allowed for the system to reach equilibrium, inter–ionic interactions pull more ionic liquid inside the confinementgap.Ionic liquids feature strong long–ranged Coulombicforces and their models require significant computationaleffort. Coarse grained models, such as the salt modelimplemented in the current study, are useful for bridg-ing the gap between the molecular processes that controlthe lubrication phenomena and the macroscopic perfor-mance in engineering applications. The implementationof simplified models that describe fundamental physico-chemical phenomena at a reduced computational cost canprovide deep insights which shed light onto the mecha-nisms and processes that can render ILs as potentiallyinteresting lubricant candidates.2
FIG. 14: Configuration snapshots ( xy cross–section) in threerepresentative points { P , P , P } . Left panels correspond tothe start of friction simulations t = 0, while right panels cor-respond to the end of friction simulations t = 3 ns. We havehighlighted the confined region with dashed lines (Top plate’swidth along the y axis is a half of the total system’s width)and also we have sketched crystallization patterns with solidlines. Top plate’s lateral velocity is set to V x = 2 m/s, to-tal simulation time is t tot = 3 ns, hence all friction simula-tions have run until the Top plate had covered a distance of d x = V x · t tot = 60 ˚A along the x direction. ACKNOWLEDGEMENTS
The work of M.D. and I.S. was supported in part bythe Serbian Ministry of Education, Science and Techno-logical Development under Project No. OI171017 and byCOST Action MP1305. Numerical simulations were runon the PARADOX supercomputing facility at the Scien-tific Computing Laboratory of the Institute of PhysicsBelgrade.
Author Contributions
K.G. and I.S. designed the study. M.D. and I.S. per-formed the simulations. They also wrote the paper withinputs from K.G..
References ∗ Electronic address: [email protected] [1] F. Zhou, Y. Liang, and W. Liu, Chemical Society Re-views , 2590 (2009).[2] R. Hayes, G. G. Warr, and R. Atkin, Physical ChemistryChemical Physics , 1709 (2010).[3] O. Y. Fajardo, F. Bresme, A. A. Kornyshev, and M. Ur-bakh, The journal of physical chemistry letters , 3998(2015).[4] O. Fajardo, F. Bresme, A. Kornyshev, and M. Urbakh,Scientific reports (2015).[5] R. Capozza, A. Vanossi, A. Benassi, and E. Tosatti, TheJournal of chemical physics , 064707 (2015).[6] Z. Wang and C. Priest, Langmuir , 11344 (2013).[7] D. A. Beattie, R. M. Espinosa-Marzal, T. T. Ho, M. N.Popescu, J. Ralston, C. J. Richard, P. M. Sellapperu-mage, and M. Krasowska, The Journal of Physical Chem-istry C , 23676 (2013).[8] A. M. Smith, K. R. Lovelock, N. N. Gosvami, T. Welton,and S. Perkin, Physical Chemistry Chemical Physics ,15317 (2013).[9] A. E. Somers, P. C. Howlett, D. R. MacFarlane, andM. Forsyth, Lubricants , 3 (2013), ISSN 2075-4442.[10] K. Gkagkas, V. Ponnuchamy, M. Dasic, and I. Stankovic,Tribology International , 83 (2017).[11] A. C. F. Mendon¸ca, A. A. H. P´adua, and P. Malfreyt,Journal of Chemical Theory and Computation , 1600(2013).[12] N. Voeltzel, A. Giuliani, N. Fillot, P. Vergne, and L. Joly,Phys. Chem. Chem. Phys. , 23226 (2015).[13] F. Federici Canova, H. Matsubara, M. Mizukami,K. Kurihara, and A. L. Shluger, Phys. Chem. Chem.Phys. , 8247 (2014).[14] K. Holmberg, P. Andersson, and A. Erdemir, TribologyInternational , 221 (2012), ISSN 0301-679X.[15] D. J. Hardy, J. E. Stone, and K. Schulten, Parallel Com-puting , 164 (2009).[16] S. Plimpton, Journal of computational physics , 1(1995).[17] S. Viscardy, J. Servantie, and P. Gaspard, The Journalof chemical physics , 184512 (2007).[18] D. J. Evans and G. Morriss, Physical Review A , 1528(1984).[19] P. J. Daivis and B. Todd, The Journal of chemical physics , 194103 (2006).[20] J. Gao, W. D. Luedtke, and U. Landman, Phys. Rev.Lett. , 705 (1997).[21] J. Gao, W. D. Luedtke, D. Gourdon, M. Ruths, J. N.Israelachvili, and U. Landman, The Journal of PhysicalChemistry B , 3410 (2004).[22] S. Perkin, Phys. Chem. Chem. Phys. , 5052 (2012).[23] R. Hayes, N. Borisenko, M. K. Tam, P. C. Howlett, F. En-dres, and R. Atkin, The Journal of Physical ChemistryC , 6855 (2011).[24] B. Bhushan, J. N. Israelachvili, and U. Landman, Nature374