Nonequilibrium molecular dynamics simulation of a dense confined nanofluid: Wall-nanoparticle interaction effects

Computational Materials Science 131 (2017) 11–20

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Nonequilibrium molecular dynamics simulation of a dense confined nanofluid: Wall-nanoparticle interaction effects C. Paredes-Arroyo a,c,⇑, R. Guzmán a,c,b a

Departamento de Ciencias Físicas, Universidad de La Frontera, Temuco, Casilla 54-D, Chile Center for Optics and Photonics, Universidad de Concepción, Concepción 4070386, Chile c Centro de Excelencia de Modelación y Computación Científica, CEMCC, Universidad de La Frontera, Temuco, Casilla 54-D, Chile b

a r t i c l e

i n f o

Article history: Received 11 August 2016 Received in revised form 6 November 2016 Accepted 10 January 2017

Keywords: Molecular dynamics Nano-colloidal suspensions Density profile Stiffness

a b s t r a c t In this article, we focus on the effects of wall-colloidal particle interaction intensity on both fluid and colloidal particles density profiles. Non-equilibrium molecular dynamics simulations were used to study a nanofluid confined between two explicit walls. Simulations were conducted keeping each wall at different temperatures by establishing a temperature gradient through a colloidal suspension. All interactions between walls, colloidal particles and fluid particles were explicitly considered in the simulations. We conclude that wall-nanoparticles interactions play a role on the features of both density and temperature profiles of nanofluids under confinement. Ó 2017 Elsevier B.V. All rights reserved.

1. Introduction Confined nanofluids or nano colloidal suspensions have been intensively studied during the last decade. Since they were proposed by Choi in 1995 [1], there have been numerous experimental and theoretical studies conducted to obtain insights of their exotic dynamical properties. Successful improvement in thermal conductivity, thermal diffusivity, and convective heat transfer of nanofluids, relative to base fluids, are constantly reported in different specialized journals. Potential applications on different technological developments include but are not limited to heat transfer intensification, mass transfer enhancements, energy storage and generation, friction reduction in mechanical devices, antibacterial activity and drug delivery [2–4]. However, until now, it is not a clear how the properties of nanofluids emerge. This question is strongly important in the case of thermodiffusive properties because they seem to display different behaviors under different preparation approaches [2]. Additionally, several inconsistencies between theoretical and computational models have been reported [3].

⇑ Corresponding author at: Departamento de Ciencias Físicas, Universidad de La Frontera, Temuco, Casilla 54-D, Chile. E-mail addresses: [email protected] (C. Paredes-Arroyo), robert. [email protected] (R. Guzmán). http://dx.doi.org/10.1016/j.commatsci.2017.01.018 0927-0256/Ó 2017 Elsevier B.V. All rights reserved.

Several theoretical approaches have been used for modeling the interactions present in nanofluids from an atomic point of view. In simple approaches, all nanoparticles interactions are modeled by using one of two interaction potentials: hard spheres [5] or shifted L-J potential expressions to the nanoparticle surface [6]. However, these approaches have resulted inadequate to properly reproduce some aspects of nanofluid dynamics. A more precise model considers each nanoparticle as a structure build by assembling a given number of atoms [7]. This approach has been suitable for constructing models of nanoparticles with different geometries and sizes. However, from a computational point of view, the size of the nanoparticle is limited by the maximum number of atoms that can be handled by a given amount of computational resources. In recent articles, this multi-atom composition model for nanoparticles has been used [8–13] but few nanoparticles were considered during these simulations. Whereas a system with one or two nanoparticles, surrounded of fluid, is useful as a model to understand some details of their interaction and dynamics. Here we propose to study a confined nanofluid with realistic values of nanoparticle volume fraction. Namely, we consider 80 nanoparticles and a volume fraction of 0:25. In order to reduce the computational cost, effective potentials were used for those interactions that involve nanoparticles [14]. Additionally, nanofluid were confined in between two walls which were kept at different temperature values. Under this conditions, we report some effects of wall-nanoparticle interaction on both density and temperature profiles of nanofluid.

12

C. Paredes-Arroyo, R. Guzmán / Computational Materials Science 131 (2017) 11–20

After a description of the computational model in Section 2, some results are presented in Section 3, and finally, conclusions are presented in Section 4. 2. Models and simulation method We consider a computational model consisting of a nanofluid confined between two parallel walls. The walls were defined perpendicular to the z axis, whereas periodic boundary conditions were imposed to x and y directions. This configuration allows the modeling of, for example, a nano-channel of infinite extent in the x and y directions. Fig. 1 shows a schematic diagram of a confined nanofluid. 2.1. Interaction potentials Each nanoparticle was considered as a uniform distribution of atoms interacting with each other by a Lennard-Jones potential


12
6  r r : UðrÞ ¼ 4  r r

ð1Þ

From this model, an effective interaction potential between nanoparticles was obtained. Of course, a nanoparticles-atom interaction naturally arises as a particular case. Hamaker in [14] reported a systematic deduction of this effective model for the case of spherical particles and applied it to study the attractive LondonVan der Waals interaction. By using this approach, it has been possible to model interfacial effects and accumulation of solvent atoms near nanoparticles boundaries with reduced computational cost [15]. Everaers et al. [16] proposed a generalization to ellipsoidal nanoparticles. Following [15] we can write effective potentials as: A U nn ðrÞ ¼ U nn ðrÞ þ U Rnn ðrÞ;

U ns ðrÞ ¼

2a3 r3ns Ans 3

9ða2  r 2 Þ " # ð5a6 þ 45a4 r 2 þ 63a2 r 4 þ 15r 6 Þr6ns  1 15ðr  aÞ6 ðr þ aÞ6

ð5Þ

where Ans ¼ 24pns qn r3ns ; rns ¼ ðrnn þ rss Þ=2. Here ss is the wall depth for solvent atoms Lenard-Jones potential, r is the separation between a given nanoparticle center and a solvent atom. The election nn ¼ ss ¼ ; Ann ¼ ; Ans ¼ 72 and rnn ¼ rss ¼ r was made in order to simulate an homogenous nanofluid under periodic boundary conditions [17]. Computational simulations were done using a classical molecular dynamics simulator, LAMMPS [18]. LAMMPS contains efficient algorithms to implement the effective potential showed in Eqs. (4) and (5). These algorithms were implemented by P. in ’t Veld and collaborators. LAMMPS implementation has been used to study properties in shear rheology of nanofluids [17] and to obtain an effective pair potential mediated by an explicit solvent [19]. Additionally, Verlet-Velocity algorithm was used to integrate equations of movement with a time step of Dt ¼ 0:005s. 2.2. Nanofluid preparation

where A U nn ðrÞ ¼ 

U Rnn ðrÞ ¼

ð2Þ

A Potentials U nn y U Rnn result after the integration of attractive [14] and repulsive [16] parts of Lenard-Jones (1) respectively. In these references, the integration was extended over all atomic components of the nanoparticle. Ann ¼ 4p2 nn q2n r6nn is called Hamaker constant and it depends on nanoparticles properties. The parameter r is the distance between 2 nanoparticles. In this article, nanoparticle radio was set at a ¼ 10rss where rss is the typical size of solvent atoms. Nanoparticle-solvent atom interaction is obtained in a similar way [15]:

Ann 6



 2  2a 2a r  4a2 þ 2 þ ln 2 2 2 r  4a r r 2

2

ð3Þ

Ann 37800 " 

#

r6nn r2  14ra þ 54a2 r2 þ 14ra þ 54a2 2ðr2  30a2 Þ r

ðr  2aÞ7

þ

ðr þ 2aÞ7



r7 ð4Þ

Fig. 1. A schematic representation of the nanofluid confined between two plates.

Following the prescription in [15] nanoparticles and solvent were prepared separately by using two identical simulation boxes and periodic boundary conditions in all directions. Then, solvent and nanoparticles were mixed to obtain the required nanofluid. In order to achieve this goal, a number of 80 nanoparticles (mass mn ¼ 83:33ms ) were prepared in a box whose volume fraction was set to /v ¼ 0:2492. Our preparation procedure departs from that reported in [15], because here, the potential

   pr ; UðrÞ ¼ A 1 þ cos rc

ð6Þ

was used as an interaction potential for nanoparticles, with A ¼ 500 and r c ¼ 11r. By using this potential in (6), the overlaps between nanoparticles were suppressed. As a result, nanoparticles are equilibrated at some fixed temperature T ¼ =kB . In a different box, a pure solvent was prepared with a number density qs ¼ 0:5441r3 . In this case, atom-atom interaction was modeled using a Lenard-Jones potential. At the end of this stage, solvent liquid phase was equilibrated at temperature T ¼ =kB . In a next stage, nanoparticles and solvent were merged in a single simulation box. As a consequence, some fluid atoms overlapped nanoparticles. In order to avoid the overlaps, all atoms separated by less than 5:00001r from a nanoparticle center were removed from the simulation box. Accordingly, a homogeneous liquid phase colloid of 80 nanoparticles suspended in a solvent with 68,700 atoms was obtained inside a simulation box with dimensions LX ¼ LY ¼ 55:073r and LZ ¼ 55:426r. This colloidal suspension was stabilized to P ¼ 0:1=r3 and T ¼ =kB by using a NosHoover T-P control [20,21]. In order to obtain the mean values for the nanofluid physical quantities, we integrated system

13

C. Paredes-Arroyo, R. Guzmán / Computational Materials Science 131 (2017) 11–20

2.3. Walls implementation In order to confine the nanofluid, it was necessary to implement a computer model for the walls. Atoms of mass mm ¼ 4:8834ms arranged in two fcc (1 1 1) planes were used to obtain each wall. A lattice density of qm ¼ 4r3 was enough to keep the solvent atoms confined. In order to obtain fcc (1 1 1) surfaces, axes orientations were defined at x ¼ ð1; 1; 0Þ; y ¼ ð1; 1; 2Þ y z ¼ ð1; 1; 1Þ, pffiffi pffiffiffi pffiffiffi with lattice periodicities lx ¼ 2r; ly ¼ 2 3 6 r; lz ¼ 3r. Fig. 2 shows the initial set up of the walls. During the simulation, atoms in the walls were restricted to harmonic oscillations around their lattice sites and all direct atom-atom interactions were neglected. Each atom in the lattice was supposed to be under a restoring force given as:

! ! ! F ðrÞ ¼ Kð r i  r 0 Þ

ð7Þ

The walls stiffness can be modified by changing values of K in (7). In this work, three different values of K were considered, namely K ¼ 500r2 ; K ¼ 850r2 and K ¼ 1200r2 . These values have been used in other simulations comprising walls [23]. It was considered that wall atoms interact with solvent atoms via LJ potential with sm ¼ 0:6 and rsm ¼ 0:75r [13]. In addition, the nanoparticle-wall interaction was modeled by using a similar approach to the case of the nanoparticle-solvent atom interactions (5). The intensities of wall-nanoparticle interaction were taken as: Anw ¼ ; 12; 24; 40 and 72. Because the results showed only marginal difference between Anw ¼ 40 and Anw ¼ 72, discussions in the follow will involve only to the three lower values of Anw . Temperature control of walls was performed using a velocityrescaling thermostat [24] and it was applied to each fcc (1 1 1) plane, independently [13]. The thermostat allowed us to fix the mean temperature of the walls and to correct the assumption of non-interacting lattice atoms. The temperature of each wall was set to a fixed value. Several simulations were performed with the temperature of the walls fixed at the same value. After that, a temperature difference DT was established between the walls. We explored three different values of DT, namely DT ¼ 0:04=kB ; DT ¼ 0:1=kB and DT ¼ 0:16=kB .

2.4. Confining the nanofluid The nanofluid was prepared under periodic boundary conditions and was permitted to evolve until its equilibrium state. The next stage was to confine the nanofluid between the two walls. An equivalent volume between the walls was filled with nanofluid components. This step was done with information about positions and velocities which were extracted from the last state in the previous stage of the simulation. In order to preserve the volume fraction of the nanofluid and to avoid the overlap between nanofluid and wall atoms, the nanofluid was compressed 2:5r along the z direction. This compression was done using two ideal walls. One specular wall was acting on solvent atoms, and a Hamaker type [14] wall, was acting on nanoparticles. After this stage, explicit walls were placed at the borders of the simulation domain. Then, ideal walls were removed and finally, the nanofluid was permitted to evolve until it adopted the new boundary condition. No thermostat was used for the nanofluid and the excess of temperature was removed by applying a temperature control to the walls. Cut-off distances for different potentials were chosen as follow: nanoparticle-nanoparticle 25r, nanoparticle-solvent atoms 9r, solvent-solvent 3r, solvent-wall 3r and nanoparticle-wall 9r. Evolution of the nanofluid under periodic boundary conditions was solved using 6  106 simulation steps. However, when the solution was under confinement, each configuration reported below, was allowed to evolve 5  106 simulation steps. In this case, averages

3 2.5 2

gnn(r)

movement equations coupled to a thermal bath with a damping constant C ¼ 0:01s1 [22].

1.5 1 0.5 0

10

15

20

25

30

35

40

r/σ

(a) 3.5 3

gns(r)

2.5 2 1.5 1 0.5 0

5

6

7

8

r/σ

9

10

11

12

(b) Fig. 2. Schematic representation of initial configuration of wall atoms.

Fig. 3. Pair correlation function for (a) nanoparticle-nanoparticle and (b) nanoparticle-solvent.

14

C. Paredes-Arroyo, R. Guzmán / Computational Materials Science 131 (2017) 11–20

of relevant quantities were obtained from the latest 2:5  106 steps of time. The nearest neighbor list was verified every 15 time step by using a neighbor skin of 2:4r. As comparison, a confined solvent was simulated following a very similar procedure to that describet above to confine the nanofluid. The confined nanofluid was simulated using 6 nodes Linux, each one equipped with 2 Intel Xeon E5-2660 10 cores processors linked with a Infiniband FDR 4X Mellanox 2 Ethernet 1 GB network. Simulations of pure confined solvent were executed on a Linux node equipped with two Intel Xeon E5-2660-10 core and a GPU Tesla k40c with 2880 cores and 12 GB RAM. For the simulations of the pure confined solvent, the Lamps GPU package was used [25–27].

3. Results and discussion 3.1. Nanofluid under periodic boundary conditions In an early stage of the simulation, we calculated the radial distribution function for each component of the nanofluid. This calculations were performed under periodic boundary conditions. Fig. 3 shows the nanoparticle-nanoparticle radial distribution function g nn ðrÞ (3a) and the nanoparticle-solvent atom radial distribution function g ns ðrÞ (3b). Our results are in correspondence with those reported before. Namely, the nanoparticle-nanoparticle RDF (Fig. 3a) is in agreement with Fig. 2 in Ref. [17]. However, little differences arise

Fig. 4. Schematic nanoparticles arrangement around the walls in steady state for DT ¼ 0:16=kB ; K ¼ 1200r2 . Here (a) Anw ¼  and (b) Anw ¼ 40. The size of solvent atoms were reduced in order to highlight the nanoparticles.

0.01

Δρn

ρn

0.008 0.006

K=850 ε/σ 2

K=1200 ε/σ 2

2.5

0.004

0.04

0.002

0.02

0

0.004

2 1.5

0 0

10

20

30

40

ρs

K=500 ε/σ 2

Δρs

0.012

1

50

0.5

0.002

0

0 0

10

20

30

40

50

z/σ Fig. 5. Density profiles calculated with walls at the same temperature and Anw ¼ . Continuous line for nanoparticles and dotted line for solvent.

Δρn

ρn

0.008 0.006

K=850 ε/σ2

K=1200 ε/σ2

2.5

0.004

0.04

0.002

0.02

0

0 0

0.004

10

20

30

40

2 1.5

ρs

K=500 ε/σ2

0.01

Δρs

0.012

1

50

0.5

0.002 0

0 0

10

20

z/σ

30

40

50

Fig. 6. Density profiles calculated with walls at the same temperature and Anw ¼ 40. Continuous line for nanoparticles and dotted line for solvent.

15

C. Paredes-Arroyo, R. Guzmán / Computational Materials Science 131 (2017) 11–20

ρn

Δρn

0.008

ΔT =0.10 ε/kB

ΔT =0.16 ε/kB

0.004

0.04

0.002

0.02

0.006

0

2 1.5

0 0

0.004

2.5

10

20

30

40

50

ρs

ΔT =0.04 ε/kB

0.01

Δρs

0.012

1 0.5

0.002 0

0 0

10

20

30

40

50

z/σ Fig. 7. Density profile for nanoparticles (continuous lines) and solvent (dotted lines) of the confined nanofluid. Three values of temperature difference between walls are considered. Anw was chosen as  and K ¼ 1200=r2 .

ρn

Δρn

0.008

ΔT=0.10 ε/kB

ΔT=0.16 ε/kB

2.5

0.004

0.04

0.002

0.02

0

0.006

2

0 0

10

20

30

40

1.5

50

0.004

1

0.002

0.5

0

ρs

ΔT=0.04 ε/kB

0.01

Δρs

0.012

0 0

10

20

30

40

50

z/σ

ρn

Fig. 8. Density profile for (a) nanoparticles and (b) solvent of the confined nanofluid. Three values of temperature difference between walls are considered. Anw was chosen as 40 and K ¼ 1200=r2 .

0.002

A=ε

ΔT=0.16 ε/kB

ΔT=0 ε/kB

0.001

0.08

ρn

0 0.002

A=12 ε

0.001

0.002

A=24 ε

ρn

ρn

0

0.001 0

0.07

0.004

0.06

0.003

0.05

0.002

0.04

0.001

0.03

ρn

0

0.002 0.001 0

A=40 ε

6.5

0.02

7

7.5

8

z/σ

0.01

15

20

25

30

35

40

z/σ

8.5 A=ε A = 12 ε A = 24 ε A = 40 ε

9

0 5.5

6

6.5

7

7.5

8

8.5

9

z/σ

Fig. 9. Density profiles of nanoparticles in middle of the nanofluid for different values of Anw (A in figure).

(a) 0.09

3.2. Nanofluid between two explicit walls Simulations of the confined nanofluid were performed following the procedure described in Section 2.

0.08 0.07 0.06

ρn

because they modeled 500 nanoparticles instead of 80 nanoparticles that we included in our calculations. Solvation effects appear in the systems because the particleparticle interactions and the nanoparticle-fluid interactions were explicitly considered for the simulations. In fact, nanoparticlesolvent radial distribution function shows that there are approximately three layers of solvent particles around each nanoparticle in the nanofluid. In this way, we verified that both nanofluid parameters and simulation procedures were properly implemented. After this verification, we proceeded to study the effects of the walls on the density and temperature profiles of the confined nanofluid.

0.05 0.04 0.03

0.005 0.004 0.003 0.002 0.001 0 46.5

47

47.5

48

48.5

49

49.5

0.02

A=ε A = 12 ε 0.01 A = 24 ε A = 40 ε 0 46.5 47

47.5

48

48.5

49

49.5

50

z/σ (b) Fig. 10. Density profiles of nanoparticles near to the coolest wall (a) and near to the hottest wall (b) of the nanofluid for different values of Anw (A in figure). Solid lines correspond to DT ¼ 0=kB and dotted lines correspond to DT ¼ 0:16=kB .

C. Paredes-Arroyo, R. Guzmán / Computational Materials Science 131 (2017) 11–20

Fig. 4 shows a schematic representation of a stationary state, obtained using LAMMPS tools. Once the stationary state was obtained, we calculated different density profiles along the z direction, perpendicular to the walls. The entire nanofluid volume was divided in 2000 slabs parallel to the xy plane. Each slab was defined as 0:0277r in width. Density profiles for nanofluid components represented by a density number in each slab (N=V). Each density number was obtained by averaging over time the total amount of nanoparticles and solvent atoms inside of each slab, and then dividing this numbers by the volume of the slab. 3.2.1. Wall stiffness Figs. 5 and 6 show density profiles of the nanofluid components. Continuous lines corresponds to the nanoparticles and dotted lines correspond to the solvent. In all cases, walls were kept at the same temperature during the simulation, and the results discussed here were obtained at an equilibrium state. Moreover, three different values of stiffness were explored for four different values of the Anw constant. Density profiles presented oscillations in the middle zone of the nanofluid for both solvent and nanoparticles. For each maximum in the nanoparticle density profile there was a minimum in the correspondent solvent density profile. Comparison between Figs. 5 and 6 showed that when Anw was increased, these maximums, were shifted symmetrically towards the walls. Consequently, there was an accumulation of nanoparticles near to the walls. On the other hand, from an analysis of Figs. 5 and 6 we concluded that the influence of the wall stiffness on density profiles was not appreciable for our calculations. In particular, differences between curves for distinct value of K were less than standard deviations (showed in the inset). Similar calculations were performed by setting each wall to different temperatures. In all these cases, there was not obtained appreciable influences of the stiffness on both density profiles and temperature profiles of the nanofluid. However, it is important to note that we worked with a special regime of fixed parameters. The influence of wall stiffness can be appreciable when another set of parameters is considered. Finally, we remark that the absence of influence of stiffness on density profiles was verified only for the walls and stiffness modeled as in Section 2. Detailed models of the wall stiffness could give account of a different behavior for this systems. Since stiffness presented little influence on the properties of the nanofluid explored here, we set the value of the stiffness to K ¼ 1200r2 in all results that are report from now on. Thus, we discuss how the nanofluid behaves when both external temperature differences and nanoparticle-wall interactions change.

3.2.2. Nanoparticle-wall interaction A global picture for the density profiles, within the entire volume of the nanofluid, is showed in Figs. 5–8. As we mentioned before, in Figs. 5 and 6 walls were kept the same temperature whereas the value of Anw was changed (and the stiffness). In Figs. 7 and 8 the values of Anw were fixed and the difference of temperature between walls were changed. In all cases, it was observed that, when Anw was raised, nanoparticles density profiles have higher peaks near to the walls. However, when a temperature difference were kept between the walls, nanoparticles were more likely to migrate towards hottest wall. This follow from the fact that the peaks of the density profiles near to the hottest wall were higher than those near to the coolest one. Density profiles near to the walls present a fine structure with until three maximums. For each minimum in the solvent density profile there were found a correspondent maximum in the density profile of the nanoparticles. When the walls where kept to different temperatures this maximums are highest near to the hottest wall an lower near to the coolest one. However, there were not substantial changes of the positions of this peaks along the z direction. For the largest value of Anw the fine structure of the density profiles of the solvent very near to the walls remains independent of the temperature difference DT between the walls, at least for the

1.6 Δ T=0.1ε/kB

Δ T=0.16ε/kB

1.2 1 0.8 0.6 0.4 0.2 0

0

10

20

30

z/σ

40

50

Fig. 12. Density profiles for pure confined solvent.

A=ε A = 12 ε A = 24 ε A = 40 ε

2.5

Δ T=0.04ε/kB

1.4

ρ

16

A=ε A = 12 ε A = 24 ε A = 40 ε

2.5

2

2

ρs

ρs

1.5 1.5

1

1

0.5

0.5

0 0.6

0.8

1

1.2

1.4

z/σ

(a)

1.6

1.8

2

2.2

0 53.2

53.4

53.6

53.8

54

54.2

54.4

54.6

54.8

z/σ

(b)

Fig. 11. Density profiles of solvente near to the coolest wall (a) and near to the hottest wall (b) of the nanofluid for different values of Anw (A in figure). Solid lines correspond to DT ¼ 0=kB and dotted lines correspond to DT ¼ 0:16=kB .

17

C. Paredes-Arroyo, R. Guzmán / Computational Materials Science 131 (2017) 11–20

1.08

0.35

1.06

0.25

vx

1.04

vy

1.02

T / (ε/kB)

0.3

0.2 vz

0.15 0.1

0.015

1

0.01

0.98 0.96

0.05

0.005

0.94

-1.5

-1

-0.5

0

v / (ε/ms)

0.5

1

1.5

2

0.92

1/2

0

10

20

0.4

1.08

P(vi)

0.35 0.3

T / (ε/kB)

vy

0.2 vz

0.15 0.1 0.05 -1.5

-1

-0.5

0

v / (ε/ms)

50

0

0.5

1

1.5

2

1/2

values explored with our simulations. On the other hand, the density profiles of the nanoparticles were reduced to a single very pronounced maximum near to the walls and the next maximum were found beyond of 15r from the wall. These considerations allows to consider that, for this regime of the nanoparticle-wall interaction energies, there will be a configuration with a layer of nanoparticles, almost fixed, near to the walls. The rest of nanoparticles will be likely to move in the middle of nanofluid around of certain fixed positions. This is consistent with the schematic picture presented in Fig. 4b. On the other hand, Figs. 5–8 shows that the height of the maximums of the density profiles for the solvent, near to the walls, decreased as the values of Anw raised from  to 40. This behavior indicates that when the nanoparticle-wall energy interaction was raised fluid particles are displaced from the walls. A possible explanation for this behavior is the formation of solvent layers around the nanoparticle and near to the walls [8]. When Anw was raised, walls attracted nanoparticles more intensively. Under this conditions, there will be interactions between solvent layers around the nanoparticle and the layer near the walls. As a consequence of this interaction some molecules of solvent were displaced away the walls and the nanoparticles approached to the walls. For comparative purposes, the density profiles for the pure solvent are shown in Fig. 12. The walls were kept at same temperature differences as in Fig. 8. Fig. 10 shows the density profiles, of the nanoparticles, for four different values of Anw . Continuous line corresponds to the walls kept at the same temperature (DT ¼ 0=kB ) and dotted lines corresponds to the walls kept at different temperature (DT ¼ 0:16=kB ). Fig. 11 shows a similar close up for the solvent density profile. Stiffness was set at K ¼ 1200r2 . Fig. 10a shows the behavior of

0.007

1.04

0.006

1.02

0.005

1

0.004

0.98

0.003

0.96

0.002

0.94

0.001

0.92

0

10

20

30

40

50

0

z/σ

(b) Fig. 13. Velocity distribution for slab adjacent to hotter wall with DT ¼ 0:16=kB y K ¼ 1200=r2 . (a) Nanofluid whit Anw ¼ 40 and (b) a pure confined solvent.

0.008

ΔT=0.04 ε/kB ΔT=0.1 ε/kB ΔT=0.16 ε/kB

1.06 vx

0.25

-2

40

z/σ (a)

(a)

0

30

ρn

0 -2

0.02

ΔT=0.04 ε/kB ΔT=0.1 ε/kB ΔT=0.16 ε/kB

P(vi)

ρn

0.4

(b) Fig. 14. Comparison between confined pure solvent () and nanofluid () temperature profiles for three different values of temperature difference between walls are showed. In (a) Anw ¼  and in (b) Anw ¼ 40. For all cases K ¼ 1200=r2 .

density profile near to coolest wall and Fig. 10b shows density profiles near to the hottest wall. Fig. 9 shows the behavior of density profiles of nanoparticles in the middle of the nanofluid. This figure complements the information showed in Fig. 10.

3.3. Temperature profiles In this section we report some effects of both walls stiffness and nanoparticle-wall interaction on temperature profile of the nanofluid. For this calculations the entire simulation box was divided in 24 slabs parallels to walls (xy plane) each one of width 2:29r. Under this conditions, we assume local equilibrium approximation to calculate temperature. In order to guarantee this local equilibrium assumption we calculate velocity distribution in three representative slabs, two adjacent to each wall and one in the middle of nanofluid. Fig. 13 shows velocity distribution (near to hottest wall) which were obtained by averaging 105 time steps. The velocity distribution present a Gaussian behavior and they are in close agreement with theoretical predictions in the middle zone of nanofluid. However, those near to the walls present slight deviations from normal behavior. This deviations increase for large values of applied temperature difference between walls, see for example the case DT ¼ 0:16=kB . Of course, this deviations are a consequence of interfacial effects which are amplified with larger values of applied temperature differences. Since this deviations are not significant for the regime of parameters explored here we assume local equilibrium in all cases.

18

C. Paredes-Arroyo, R. Guzmán / Computational Materials Science 131 (2017) 11–20

Anw = ε Anw = 40 ε solvent pure

T / (ε/kB)

1.06

0.014

1.04

0.012

1.02

0.01

1

0.008

0.98

ρn

1.08

0.006

0.96 0.004

0.94

0.002

0.92 0.9

0

10

20

30

40

50

0

60

z/σ

(a) Anw = ε Anw = 40 ε solvent pure

T / (ε/kB)

1.06

0.014

1.04

0.012

1.02

0.01

1

0.008

0.98

ρn

1.08

0.006

0.96 0.004

0.94

0.002

0.92 0.9

0

10

20

30

40

50

60

70

80

0

z/σ

(b)

On the other hand, when Anw  40, both nanofluid and solvent temperature profiles have the same behavior in regions near to the walls whereas in the middle regions nanofluid temperature profile keeps lower than solvent temperature profile. Moreover, nanofluid-wall temperature-jumps are the same as solvent-wall ones. A possible interpretation of this behavior in temperature profiles can be obtained from Figs. 5 and 6. As we mention before, density profile of nanoparticles shows that nanoparticles prefer to stay in certain specific positions through the nanofluid along the z direction. The separation between near nanoparticles is such that their interaction is weak [17]. Thus, this array of nanoparticles prevents to distribute thermal energy efficiently. However, when wall-nanoparticle interaction is less intense, density profiles shows an accumulation of nanoparticle near to hottest wall, so near this wall nanoparticles and solvent are very active because energy thermal energy is injected there. Furthermore, when Anw  40 Fig. 6 shows that nanoparticles are fastened to both wall so in this regions thermal energy distribution is dominate by solvent particles. Fig. 15 shows the temperature profiles for the nanofluid confined in two nano-channels. In Fig. 15a a nano-channel of width L ¼ 69r was considered, whereas, in Fig. 15b nano-channel width was set to L ¼ 82:5r. The temperature profiles present the same properties as that presented in Fig. 14. An interesting behavior is observed when we change the nanoparticles volume fraction /. As we can see from the Fig. 16, for low values of /, the temperature and density profiles were independent of Anw . As / is approaching to / ¼ 0:25 the behavior is similar to those discussed above.

Fig. 15. Temperature profiles of a nanofluid confined. Volume fraction of the nanoparticles were kept to / ¼ 0:25. Chanel width was set to L ¼ 69r in (a) and L ¼ 82:5r (b). Here DT ¼ 0:16=kB .

Anw=ε, ΔT=0 ε/kB Anw=40ε, ΔT=0 ε/kB

9000 8000

Anw=ε, ΔT=0.16 ε/kB Anw=40ε, ΔT=0.16 ε/kB

1.1

MSD / (σ2)

7000 0.03

nanofluid φV=0.25 nanofluid φV=0.14 nanofluid φV=0.03 solvent pure

0.025

1.05

0.015

1

ρn

T / (ε/k B)

0.02

6000 5000

9500

4000

9000

3000

8500

2000

8000

1000

7500

0 0.01

3.8 3.85 3.9 3.95

0

0.5

1

1.5

10

20

30

40

50

3

4

3.5

4

3.5

4

6

(a)

0.005

0

2.5

10 x Timestep

0.95

0.9

2

900

0

z/σ

800

Fig. 14 shows temperature profiles for both nanofluid and pure solvent (without nanoparticles) confined between the walls. First at all, when nanoparticles are added, temperature decrease in the middle zones of the nanofluid compared to solvent temperature. The difference between temperature profiles of nanofluid and those of solvent increase proportionally to the values of temperature difference between walls. When Anw is smallest than Ans temperature decrease near to coolest wall whereas it increase near hottest one. Moreover, for this values of Anw , nanofluid-wall interface temperature jumps decrease at both walls.

MSDz / (σ2)

700 Fig. 16. Temperature profiles for a the nanofluid at different volume fractions. Here DT ¼ 0:16=kB ; K ¼ 1200=r2 . Anw ¼  (continuous line) and Anw ¼ 40 (dotted line).

600 500 400 300 200 Anw=ε, ΔT=0 ε/kB Anw=40ε, ΔT=0 ε/kB

100 0

0

0.5

1

Anw=ε, ΔT=0.16 ε/kB Anw=40ε, ΔT=0.16 ε/kB

1.5

2

2.5

3

106 x Timestep

(b) Fig. 17. Mean square displacement for nanofluid. Total MSD (a) and MSD for z component (b).L ¼ 55:5r (solid), L ¼ 69r (dot-dot) and L ¼ 82:5r (dash-dash). For all K ¼ 1200=r2 .

C. Paredes-Arroyo, R. Guzmán / Computational Materials Science 131 (2017) 11–20 4.5

Anw=ε, ΔT=0 ε/kB Anw=40ε, ΔT=0 ε/kB Anw=ε, ΔT=0.16 ε/kB Anw=40ε, ΔT=0.16 ε/kB

4 3.5

g(r)

3 2.5 2 1.5 1 0.5 0 10

11

12

13

14

15

16

z/σ (a) 6 3.65 3.6 3.55 3.5 3.45 3.4 3.35 3.3

5

g(r)

4 3

5.82

5.84

5.86

5.88

5.9

5.92

2 1 Anw=ε, ΔT=0 ε/kB Anw=40ε, ΔT=0 ε/kB 6

solvent interaction, density profiles show that thermophoretic effects are present for all values of temperature between walls explored here. This thermophoretic effects cause an accumulation of nanoparticles at hottest wall and solvent particles at coolest one. Thermophoretic effects are attenuated, but no suppressed, when the nanoparticle-wall interaction is comparable to nanoparticlesolvent interaction. Temperature profiles inside of nanosuspension are affected by the presence of nanoparticle-wall interactions. If we compare temperature profiles for both nanofluid and pure solvent they departs at the middle zones. If nanoparticle-wall interaction is comparable to nanoparticle-solvent interaction nanofluid temperature profiles properties, near to the walls, are dominated by solvent. This results suggest that it is necessary a further exploration of thermal transport properties in the regime of parameter namely large nanoparticle numbers and confinement. Acknowledgments

5.8

0

19

6.5

Anw=ε, ΔT=0.16 ε/kB Anw=40ε, ΔT=0.16 ε/kB 7

7.5

8

8.5

C. Paredes-Arroyo dedicates this work to the Family ParedesArroyo and to Constanza Valenzuela. C. Paredes-Arroyo acknowledges partial financial support from the Master of Physics Program of Universidad de La Frontera. C. Paredes-Arroyo and R. Guzmán acknowledge CEMCC UFRO. R. Guzmán acknowledges partial financial support by the Center for Optics and Photonics (CEFOP) FB0824/2008. Powered@NLHPC: This research was partially supported by the supercomputing infrastructure of the NLHPC (ECM-02).

z/σ

(b) Fig. 18. Radial distribution function for the components of the nanofluid confined. Nanoparticle-nanoparticle (a) and nanoparticle-solvent (b). Here K ¼ 1200=r2 . L ¼ 55:5r (solid), L ¼ 69r (dot-dot) and L ¼ 82:5r (dash-dash).

Fig. 17 shows the mean square displacement for the nanofluid. The MSD for the z direction changed for different values of the nano-channel width. Fig. 18 shows the RDF calculated for the components of the nanofluid under confinement. It is interesting to note that the RDF for nanoparticle-solvent present a similar behavior as the RDF calculated for periodic boundary conditions. The RDF for nanoparticle-nanoparticle shows that the number of neighbors increase when the temperature difference between walls raise, the Anw decrease and the nano-channel width decrease. 4. Conclusions In this work we report an implementation of molecular simulation of a nanosuspension confined between two explicit walls. After replicate some results reported elsewhere we use this procedure to study effects of wall interactions on nanosuspension properties. We explore how wall stiffness and nanoparticle-wall interactions affect both particles density profiles and temperature profiles inside of suspension. From the previous section we can conclude that, in this regime of nanofluid parameters, there are a little influence of wall stiffness. However, it is important to note that it is expect that, for another set of nanofluid parameters, stiffness effects will be more important at nanofluid-wall interface. On the other hand, temperature and density profiles of confined nanofluid evidence a strong dependence on nanoparticle-wall interaction. If this interaction is much smaller than nanoparticle-

References [1] S. Choi, J. Eastman, Enhancing thermal conductivity of fluids with nanoparticles, Dev. Appl. Non-Newtonian Flows FED-231/MD-66 (1995) 99– 105. [2] W. Yu, H. Xie, A review on nanofluids: preparation, stability mechanisms, and applications, J. Nanomaterials 2012 (2012) 435873, http://dx.doi.org/10.1155/ 2012/435873, 17 pages. [3] M. Bahiraei, Particle migration in nanofluids: a critical review, Int. J. Therm. Sci. 109 (2016) 90–113, http://dx.doi.org/10.1016/j.ijthermalsci.2016.05.033. [4] S. Kumar, S. Choi, H. Patel, Heat transfer in nanofluids – a review, Heat Transf. Eng. 27 (2006) 3–19, http://dx.doi.org/10.1080/01457630600904593. [5] M. Dijkstra, R. van Roij, R. Evans, Phase diagram of highly asymmetric binary hard-sphere mixtures, Phys. Rev. E 59 (1999) 5744–5771, http://dx.doi.org/ 10.1103/PhysRevE.59.5744. [6] C. Powell, N. Fenwick, F. Bresme, N. Quirke, Wetting of nanoparticles and nanoparticle arrays, Colloids Surf. 206 (2002) 241–251, http://dx.doi.org/ 10.1016/S0927-7757(02)00079-1. [7] S. Knauert, J. Douglas, F. Starr, The effect of nanoparticle shape on polymernanocomposite rheology and tensile strength, J. Polm. Sci. B Polym. Phys. 45 (2007) 1882–1897, http://dx.doi.org/10.1002/polb.21176. [8] J. Lv, W. Cui, M. Bai, X. Li, Molecular dynamics simulation on flow behavior of nanofluids between flat plates under shear flow condition, Microfluid. Nanofluid. 78 (2011) 475–480, http://dx.doi.org/10.1007/s10404-010-0684-2. [9] C. Hu, M. Bai, J. Lv, P. Wang, X. Li, Molecular dynamics simulation on the friction properties of nanofluids confined by idealized surfaces, Tribol. Int. 78 (2014) 152–159, http://dx.doi.org/10.1016/j.triboint.2014.05.018. [10] C. Sun, W. Lu, B. Bai, J. Liu, Novel flow behaviors induced by a solid particle in nanochannels: Poiseuille and Couette, Chin. Sci. Bull. 59 (2014) 2478–2485, http://dx.doi.org/10.1007/s11434-014-0282-x. [11] H. Aminfar, M.A. Jafarizadeh, N. Razmara, Nanoparticles aggregation in nanofuid flow through nanochannels: insights from molecular dynamic study, Int. J. Mod. Phys. C 25 (2014) 1450066, http://dx.doi.org/10.1142/ S0129183114500661, 23 pages. [12] W. Cui, Z. Shen, J. Yang, S. Wu, Molecular dynamics simulation on flow behaviors of nanofluids confined in nanochannel, Case Stud. Therm. Eng. 5 (2015) 114–121, http://dx.doi.org/10.1016/j.csite.2015.03.007. [13] M. Frank, D. Drikakis, N. Asproulis, Thermal conductivity of nanofluid in nanochannels, Microfluid Nanofluid. http://dx.doi.org/10.1007/s10404-0151591-3. [14] H. Hamaker, The London–van der Waals attraction between spherical particles, Physica IV 10 (1937) 1058–1072, http://dx.doi.org/10.1016/S00318914(37)80203-7. [15] P. in ’t Veld, S. Plimpton, G. Grest, Accurate and efficient methods for modeling colloidal mixtures in an explicit solvent using molecular dynamics, Comput.

20

[16] [17]

[18] [19]

[20]

[21]

C. Paredes-Arroyo, R. Guzmán / Computational Materials Science 131 (2017) 11–20 Phys. Commun. 179 (2008) 320–329, http://dx.doi.org/10.1016/j. cpc.2008.03.005. R. Everaers, M. Ejtehadi, Interaction potentials for soft and hard ellipsoids, Phys. Rev. E 67 (2003) 1–8, http://dx.doi.org/10.1103/PhysRevE.67.041710. P. in ’t Veld, M. Petersen, G. Grest, Shear thinning of nanoparticle suspensions, Phys. Rev. E 79 (2009) 021401, http://dx.doi.org/10.1103/PhysRevE.79.021401, 4 pages. S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys. 117 (1995) 1–19, http://dx.doi.org/10.1006/jcph.1995.1039. G. Grest, Q. Wang, P. in ’t Veld, D. Keffer, Effective potentials between nanoparticles in suspension, J. Chem. Phys. 134 (2011) 144902, http://dx.doi. org/10.1063/1.3578181, 7 pages. S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys 81 (1984) 511–519, http://dx.doi.org/ 10.1063/1.447334. W. Hoover, Canonical dynamics: equilibrium phase-space distributions, Phys. Rev. A. 31 (1985) 1695–1697, http://dx.doi.org/10.1103/PhysRevA.31.1695.

[22] T. Schneider, E. Stoll, Molecular-dynamics study of a three-dimensional onecomponent model for distortive phase transitions, Phys. Rev. B 17 (1978) 1302–1322, http://dx.doi.org/10.1103/PhysRevB.17.1302. [23] N. Asproulis, D. Drikakis, Boundary slip dependency on surface stiffness, Phys. Rev. E. 81 (2010) 061503, http://dx.doi.org/10.1103/PhysRevE.81.061503, 5 pages. [24] L. Woodcock, Isothermal molecular dynamics calculations for liquid salts, Chem. Phys. Let. 10 (1971) 257–261, http://dx.doi.org/10.1016/0009-2614(71) 80281-6. [25] W.M. Brown, P. Wang, S.J. Plimpton, A.N. Tharrington, Implementing molecular dynamics on hybrid high performance computers – short range forces, Comp. Phys. Comm. 182 (2011) 898–911. [26] W.M. Brown, A. Kohlmeyer, S.J. Plimpton, A.N. Tharrington, Implementing molecular dynamics on hybrid high performance computers – particle-particle particle-mesh, Comp. Phys. Comm. 183 (2012) 449–459. [27] W.M. Brown, Y. Masako, Implementing molecular dynamics on hybrid high performance computers – three-body potentials, Comp. Phys. Comm. 184 (2013) 2785–2793.