# Synchronized molecular-dynamics simulation for the thermal lubrication of a polymeric liquid between parallel plates

## Abstract

The Synchronized Molecular-Dynamics simulation which was recently proposed by authors [Phys. Rev. X 4, 041011 (2014)] is applied to the analysis of polymer lubrication between parallel plates. The rheological properties, conformational change of polymer chains, and temperature rise due to the viscous heating are investigated with changing the values of thermal conductivity of the polymeric liquid. It is found that at a small applied shear stress on the plate, the temperature of polymeric liquid only slightly increases in inverse proportion to the thermal conductivity and the apparent viscosity of polymeric liquid is not much affected by changing the thermal conductivity. However, at a large shear stress, the transitional behaviors of the polymeric liquid occur due to the interplay of the shear deformation and viscous heating by changing the thermal conductivity. This transition is characterized by the Nahme-Griffith number which is defined as the ratio of the viscous heating to the thermal conduction at a characteristic temperature. When the Nahme-Griffith number exceeds the unity, the temperature of polymeric liquid increases rapidly and the apparent viscosity also exponentially decreases as the thermal conductivity decreases. The conformation of polymer chains is stretched and aligned by the shear flow for , but the coherent structure becomes disturbed by the thermal motion of molecules for .

###### keywords:

Multiscale modeling, Polymeric liquid, Lubrication, Rheology^{1}

[cor1]Email: yasuda@sim.u-hyogo.ac.jp, Tel/Fax +81(0)783031990

## 1 Introduction

Molecular dynamics (MD) simulations are used to investigate the materials properties and internal structures of the complex liquids and soft matters. In MD simulations, the dynamics of all the model molecules which compose of the materials are calculated. The material properties and internal structures of the materials are analyzed from the results of the MD simulations. The MD simulation is applicable for a wide variety of complex liquids and soft matters consist of the complicated molecules such as the polymer molecules and the biological molecules. However, due to the enormous computational time to calculate the dynamics of all the molecules, the MD simulation is not yet applicable to problems which concern mesoscopic and macroscopic scales motions far beyond the molecular size, where the dynamics of molecules are strongly coupled to the macroscopic transports of the complex liquids. At microscopic scales, the macroscopic transports are usually ignored and the MD simulation is performed under a certain ideal environment. To develop a simulation technology to analyze the molecular dynamics coupled with the macroscopic transports for the complex liquids is a significant challenge from both a scientific and engineering point of view. Multiscale modeling is a promising candidate to address this challenge.

The multiscale simulation for the flow behaviors of complex fluids was first advanced in the CONNFFESSIT approach for polymeric liquids by Laso and Öttinger(4); (5); (6), where the local stress in the fluid solver is calculated using the microscopic simulation instead of using any constitutive relations. The GENERIC approach is also presented for the non-isothermal polymeric flows by making important corrections and clarifications to the CONNFFESSIT scheme.(7) The strategy exploited in the CONNFFESSIT approach is also introduced into the heterogeneous multiscale modeling (HMM), which was proposed by E and Enquist,(8) as a general methodology for the efficient numerical computation of problems with multiscale characteristics. HMM has been applied to various problems, such as the simple polymeric flow(9), coarsening of a binary polymer blend(10) and the channel flow of a simple Lennard-Jones liquid(11). The equation-free multiscale computation proposed by Kevrekids et al. is also based on a similar idea and has been applied to various problems.(12); (13) De et al. proposed the scale-bridging method, which can correctly reproduce the memory effect of a polymeric liquid, and demonstrated the non-linear viscoelastic behavior of a polymeric liquid in slab and cylindrical geometries.(14); (15) The multiscale simulation for polymeric flows with the advection of memory in two and three dimensions was developed by Murashima and Taniguchi.(16); (17); (18) We have also developed a multiscale simulation coupling of MD and transport equations. The multiscale method was first developed for simple fluids(19) and subsequently extended to polymeric liquids with the memory effect.(20); (21); (22); (23) Recently, we proposed the synchronized molecular dynamics (SMD) method in which the multiscale method was extended to treat the coupled heat and momentum transfer of complex liquids.(24)

In the previous study, the SMD method was applied to the polymer lubrication generating the viscous heating between parallel plates. The rheological properties, conformation of polymer chains, and temperature rise due to the viscous heating were investigated for various applied shear stresses on the plate. We found that an interesting transitional behavior of the conformation of polymer chains occurs with a rapid temperature rise for a large applied shear stress due to the coupling of the shear deformation and heat generation under the shear flow.

In this study, we investigate the effects of changing the thermal conductivity of the polymeric liquid on the behaviors of rheological properties, conformation of polymer chains, and heat generation in the polymer lubrication between parallel plates. In the following, we first describe the problem considered in the present paper. The simulation method is briefly explained after the presentation of the problem. The SMD simulation of polymer lubrication is performed, and the results are discussed; these results are mainly the rheological properties and the coupling of the conformations of the polymer chains with heat generation under shear flows. Finally, a short summary is given.

## 2 Problem and Method

We consider a polymeric liquid contained in a gap of width between parallel plates with a constant temperature (see Fig. 1(a)). The polymeric liquid is composed of short Kremer-Grest chains(25) of ten beads, in which all of the bead particles interact with a truncated Lennard-Jones potential defined by

(1) |

and consecutive beads on each chain are connected by an anharmonic spring potential,

(2) |

with =30 and =. The polymeric liquid is in a quiescent state with a uniform density and a uniform temperature before a time . Hereafter, the -axis is perpendicular to the parallel plates, and the boundaries between the upper and lower plates and the polymeric liquid are located at and 0, respectively. The upper plate starts to move in the direction with a constant shear stress at a time =0, while the lower plate is at rest.

The macroscopic behavior of the polymeric liquid is described by the following transport equations:

(3a) | ||||

(3b) |

where is the velocity, is the stress tensor, is the internal energy per unit mass, and is the shear rate, i.e., . Hereafter, the subscripts , , and represent the index in Cartesian coordinates, i.e., {,,}{,,}. Here, we assume that the macroscopic quantities are uniform in the and directions, ==0, and the density of the polymeric liquid is constant. Fourier’s law for a heat flux with a constant and uniform thermal conductivity is also considered in Eq. (3b). Note that the thermal conductivity of polymeric liquids is anisotropic under shear flows in general(26); (27); (28); (29), and some experimental studies have reported that the linear stress-thermal relation between the stress tensor and thermal conductivity holds.(30); (31); (32); (33) However, in the present study, we only consider the isotropic thermal conductivity as the first step because the effect of shear thinning of the viscosity is thought to be more crucial to viscous heating under strong shear flows than that of the anisotropy of the thermal conductivity. We also assume that the velocity and temperature of the polymeric liquid are the same as those of the plates at the boundaries, i.e., the non-slip and non-temperature-jump boundary conditions.

The effect of viscous heating is estimated using the ratio of the first and second terms in Eq. (3b) to be . Here, is the gross shear rate of the system, which is defined by the ratio of the velocity of the upper plate to the width of the gap , , and is a characteristic temperature rise for the polymeric liquid. In the present problem, we consider a characteristic temperature necessary to substantially change the viscosity of the polymeric liquid, i.e., , where is the characteristic viscosity of the polymeric liquid at a temperature of . Thus, the Nahme-Griffith number , defined as

(4) |

represents the effect of viscous heating on the changes in the rheological properties.(3); (34) Usually, in lubrication systems and in high-speed processing operations with polymeric liquids, the Nahme-Griffith number is not negligibly small because of the large viscosity and the small thermal conductivity of the polymeric liquid.(34) For example, when a lubrication oil in a gap with a width of 10 m is subjected to shear deformation with a strain rate of , the Nahme-Griffith number is estimated to be . Thus, the rheological properties of the lubricant in such micro devices must be significantly affected not only by the large velocity gradient but also by the temperature increase caused by local viscous heating. To predict the rheological properties of the polymeric liquid in these systems, one must consider the temperature variation in Eq. (3b) coupled with Eq. (3a).

The viscous heating is one of the most fundamental and important problems in the lubrication systems. However, it is difficult to analyze the behaviors of polymer chains under the coupling of strong shear flows and viscous heating with using the conventional simulation technologies. As is seen in Eq. 4, the effect of viscous heating on the temperature and behaviors of polymer chains is enhanced as the system size is large. Thus, the thermal lubrication is one of the typical problems that the macroscopic transports and microscopic molecular dynamics are strongly interacted each other. We solve this problem by using the Synchronized Molecular-Dynamics (SMD) simulation which was proposed recently by authors(24).

In this paper, we briefly explain the SMD method. The molecular dynamics cells are assigned to small fluid elements with a interval and calculate the local stresses, temperatures, and conformations of polymer chains in each MD cell by using the non-equilibrium molecular dynamics (NEMD) simulations according to the local shear rates . (see Fig. 1(b)). The NEMD simulations are performed independently in each MD cell for a certain time duration , but at each time interval , the NEMD simulations are synchronized to satisfy the macroscopic heat- and momentum-transport equations Eq. (3). Thus, the couplings of the shear flow, viscous heading, and conformational dynamics of polymer chains are autonomously calculated in each MD cell with satisfying the macroscopic heat- and momentum-transfers.

## 3 Results

In this study, we investigate the effect of changing the thermal conductivity of the polymeric liquid on the rheological property, conformation of polymer chains, and temperature rise. We carry out the SMD simulations with changing the value of the thermal conductivity . (In the previous study(24), a fixed value of the thermal conductivity is only considered while the applied shear stress to the upper plate is varied variously.) Hereafter, we measure the length, time, temperature and density in units of , , , and , respectively. Here, is the Boltzmann constant, and is the mass of the LJ particle. In the following simulations, the density of the polymeric liquid is fixed to be , and the temperature of the plates and the width of the gap between the plates are fixed to be and , respectively. At this density and this temperature , the conformation of the bead particles becomes severely jammed and results in complicated rheological properties.(22); (35) Two values of the applied shear stress on the upper plate are considered in this study, i.e., =0.01 and 0.05. The SMD simulations are performed with 50, 100, 150, 200, and 300 for 0.01 and 50, 75, 100, 125, 150, 200, and 500 for , respectively. In the following, we present quantities averaged for a long period of time at the steady state, in which the shear stress is spatially uniform and the time derivative of the local internal energy, i.e., the left-hand side of Eq. (3b), is negligible.

Figure 2 shows the dependency of the apparent viscosity , defined as , on the thermal conductivity and the spatial average of temperature . For a small applied shear stress, i.e., , the apparent viscosity only slightly decreases as the thermal conductivity decreases (the inverse of thermal conductivity increases). The average temperature increases linearly as the function of the inverse of thermal conductivity , but the variation is very small; the temperature rise is at most 2.5% of the wall temperature in the figure. The facts that the variations of apparent viscosity and average temperature are small and the average temperature is a linear function of the inverse of thermal conductivity can be explained by the energy transport equation, Eq. (3b), for small Nahme-Griffith numbers; when the Nahme-Griffith number is small the viscosity of the polymeric liquid is not much affected due to the viscous heating, so that the local viscosities and shear rates are uniform between the parallel plates and the temperature rise is inversely proportional to the thermal conductivity. In Fig. 2(a), the Nahme-Griffith number is at most 0.06.

For a large applied shear stress, i.e., , the behaviors of the apparent viscosity and average temperature are quite different from those for a small applied shear stress. For , the Nahme-Griffith number of the polymeric liquid can be larger than unity, , at small thermal conductivities. In Fig. 2(b), the points at which the Nahme-Griffith number equals unity are indicated by the vertical arrows and the Nahme-Griffith number is larger than unity on the right side of the vertical arrow. It is seen that the average temperature rapidly increases when the Nahme-Griffith number exceeds the unity, and the apparent viscosity also rapidly decreases as the function of the inverse of thermal conductivity. When the Nahme-Griffith number exceeds the unity, the viscosity substantially deceases due to the viscous heating, so that the shear velocity increases and the viscous heating is enhanced. This positive feedback causes the rapid changes of the apparent viscosity and average temperature for small thermal conductivities.

Figure 3 shows the conformational changes in polymer chains as a function of the inverse of thermal conductivity. Here, the bond-orientation tensor is defined as,

(5) |

where is the number of polymer chains in each MD cell, is the number of bead particles in a polymer chain, for is the bond vector between consecutive beads in the same chain, and is the distance at which the sum has a minimum and is calculated to be . As the polymer chains are stretched and aligned by the shear flows driven by the applied shear stress to the upper plate, the difference of the - and -components of the bond orientation tensor, , and the component of the bond orientation tensor becomes large. It is seen that, for a small applied shear stress (Fig. 3(a)), the conformations of polymer chains are slightly stretched and aligned by the shear flows. The variations of each component of the bond orientation to the thermal conductivity are very small; the relative variations are at most 0.3% for and 5% for .

For a large applied shear stress (Fig. 3(b)), the transitional behaviors occurs around the point at which the Nahme-Griffith number equals unity. For large thermal conductivities (on the left side of downward arrow), the conformations of polymer chains are highly stretched and aligned by the shear flows. However, for small thermal conductivities (on the right side of downward arrow), the alignment of the polymer chains is disturbed due to the thermal motion of molecules and the conformation of the polymer chains becomes more isotropic although the shear velocities are enhanced as the temperatures increase. This indicates that the transition between the regime dominated by the shear flow and dominated by the thermal motion of molecules occurs by changing the thermal conductivity.

## 4 Summary

We carry out the SMD simulations for the thermal lubrication of a model polymeric liquid composed of short chains between parallel plates. The SMD simulations for the same problem were also carried out in Ref. (24), where the behaviors of the polymeric liquid for various applied shear stresses are investigated. In the previous study we found that the transitional behaviors occur in the dynamics of polymer chains by changing the applied shear stress; when the Nahme-Griffith number is small, , the conformation of polymer chains is stretched and aligned by the shear flow, but when the Nahme-Griffith number exceeds unity, , the coherent structure becomes disturbed by the thermal motion of molecules.

In the present study, the effects of changing the thermal conductivity of the polymeric liquid on the rheological properties, conformation of polymer chains, and temperature rise due to the viscous heating are investigated. It is found that at a small applied shear stress on the plate, i.e., , the temperature of polymeric liquid only slightly increases in inverse proportion to the thermal conductivity and the apparent viscosity of polymeric liquid is not much affected by changing the thermal conductivity since the Nahme-Griffith number is very small. However, at a large shear stress, i.e., , the Nahme-Griffith number of the polymeric liquid can be larger than unity, , when the thermal conductivity is sufficiently small. When the Nahme-Griffith number exceeds unity, the temperature of polymeric liquid increases rapidly and the apparent viscosity exponentially decreases as the thermal conductivity decreases. The conformation of polymer chains is stretched and aligned by the shear flow for , but the coherent structure becomes disturbed by the thermal motion of molecules for . Thus, in this study we demonstrate that the transitional behaviors which were found in the previous study also occur by changing the thermal conductivity under a fixed applied shear stress.

## acknowledgement

This study was financially supported by JSPS KAKENHI Grant Numbers 26790080 and 26247069. The computations have been performed using the supercomputer of ACCMS, Kyoto University.

### Footnotes

- journal:

### References

- M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, (Oxford University Press, New York, 1989).
- D. J. Evans and G. Morris, Statistical mechanics of nonequilibrium liquids, (Cambridge university press, New York, 2008).
- R. B. Bird, R. C. Armstrong, and O. Hassager, Dynamics of polymeric liquids Vol. 1 (John Wiley and Sons, New York, 1987).
- M. Laso and H. C. Öttinger, “Calculation of viscoelastic flow using molecular models: the CONNFFESSIT approach”, J. Non-Newtonian Fluid Mech. 47, 1 (1993).
- K. Feigl, M. Laso, and H. C. Öttinger, “CONNFFESSIT approach for solving a two-dimensional viscoelastic fluid problem”, Macromolecules 28, 3261 (1995).
- M. Laso, M. Picasso, H. C. Öttinger, “2-D time-dependent viscoelastic flow calculations using CONNFFESSIT”, AIChE J. 43, 877 (1997).
- M. Dressler, B. J. Edwards, Öttinger, “Macroscopic thermodynamics of flowing polymeric liquids”, Rheol. Acta 38, 117 (1999).
- W. E and B. Engquist, “The heterogeneous multi-scale methods”, Comm. Math. Sci. 1, 87 (2003).
- W. Ren and W. E, “Heterogeneous multiscale method for the modeling of complex fluids and micro-fluidics”, J. Compt. Phys. 204, 1 (2005).
- M. M̈uller and K. C. Daoulas, “Speeding Up Intrinsically Slow Collective Processes in Particle Simulations by Concurrent Coupling to a Continuum Description”, Phys. Rev. Lett. 107, 227801 (2011).
- M. K. Borg, D. A. Lockerby, J. M. Reese, “A multiscale method for micro/nano flows of high aspect ratio”, J. Compt. Phys. bf 233, 400 (2013).
- I. G. Kevrekidis, C. W. Gear, J. M. Hyman, P. G. Kevrekidis, O. Runborg, and C. Theodoropoulos, “Equation-free, coarse-grained multiscale computation: enabling microscopic simulations to perform system-level analysis”, Comm. Math. Sci. 1, 715 (2003).
- I. G. Kevrekidis and G. Samaey, “Equation-free multiscale computation: algorithms and applications”, Annu. Rev. Phys. Chem. 60, 321 (2009).
- S. De, J. Fish, M. S. Shephard, P. Keblinski, and S. K. Kumar, “Multiscale modeling of polymer rheology”, Phys. Rev. E 74, 030801(R) (2006).
- S. De, “Computational study of the propagation of the longitudinal velocity in a polymer melt contained within a cylinder using a scale-bridging method”, Phys. Rev. E 88, 052311 (2013).
- T. Murashima and T. Taniguchi, “Multiscale Lagrangian Fluid Dynamics Simulation for Polymeric Fluid” J. Polym. Sci. B 48, 886 (2010).
- T. Murashima and T. Taniguchi, “Multiscale Simulation of History Dependent Flow in Polymer Melt”, Europhys. Lett. 96, 18002 (2011).
- T. Murashima and T. Taniguchi, “Multiscale Simulation of History Dependent Flow in Polymer Melt”, J. Phys. Soc. Jpn. 81, SA013 (2012).
- S. Yasuda and R. Yamamoto, “A model for hybrid simulation of molecular dynamics and computational fluid dynamics”, Phys. Fluids 20, 113101 (2008).
- S. Yasuda and R. Yamamoto, “Rheological properties of polymer melt between rapidly oscillating plates: an application of multiscale modeling”, Europhys. Lett. 86, 18002 (2009).
- S. Yasuda and R. Yamamoto, “Multiscale modeling and simulation for polymer melt flows between parallel plates”, Phys. Rev. E 81, 036308 (2010).
- S. Yasuda and R. Yamamoto, “Dynamic rheology of a supercooled polymer melt in nonuniform oscillating flows between rapidly oscillating plates”, Phys. Rev. E 84, 031501 (2011).
- T. Murashima, S. Yasuda, T. Taniguchi, and R. Yamamoto, “Multiscale modeling for polymeric flow: particle-fluid bridging scale methods”, J. Phys. Soc. Jpn. 82, 012001 (2013).
- S. Yasuda and R. Yamamoto, “Synchronized molecular-dynamics simulation via macroscopic heat and momentum transfer”, Phys. Rev. X 4, 041011 (2014).
- K. Kremer and G. S. Grest, “Dynamics of entabgled linear polymer melts: A molecular-dynamics simulation,” J. Chem. Phys. 92, 5057 (1990).
- B.H.A.A. van den Brule and S.B.G. O’Brien, “Anisotropic conduction of heat in a flowing polymeric material”, Rheol Acta 29, 580 (1990).
- H. C. Öttinger and F. Petrillo, “Kinetic theory and transport phenomena for a dumbbell model under nonisothermal conditions”, J. Rheol. 40, 857 (1996).
- R. B. Bird and C. F. Curtiss, “Nonisothermal polymeric fluids”, Rheol. Acta 35, 103 (1996).
- R. B. Bird, C. F. Curtiss, and K. J. Beers, “Polymer contribution to the thermal conductivity and viscosity in a dilute solution”, Rheol. Acta 36, 269 (1997).
- D. C. Venerus, J. D. Schieber, H. Iddir, J. Guzmán, and A. Broerman, “Anisotropic Thermal Diffusivity Measurements in Deforming Polymers and the Stress-Thermal Rule”, Int. J. Thermophysics 22, 1215 (2001).
- J. D. Schieber, D. C. Venerus, K. Bush, V. Balasubramanian, and S. Smoukov, “Measurement of anisotropic energy transport in flowing polymers by using a holographic technique”, PNAS 101, 13142 (2004).
- J. D. Schieber, D. C. Venerus, and S. Gupta, “Molecular origins of anisotropy in the thermal conductivity of deformed polymer melts: stress versus orientation contributions”. Soft Matter 8, 11781 (2012).
- S. Gupta, J. D. Schieber, and D. C. Venerus, “Anisotropic thermal conduction in polymer melts in uniaxial elongation flows”, J. Rheol. 57, 427 (2013).
- C. J. Pipe, T. S. Majmudar, and G. H. McKinley, “High shear rate viscometry”, Rheol Acta 47, 621 (2008).
- R. Yamamoto and A. Onuki, “Dynamics and rheology of a supercooled polymer melt in shear flow,” J. Chem. Phys. 117, 2359 (2002).