**1. Introduction**

Optimum design of heat exchanger using nanotechnology is a burning field to reduce the energy consumption. Recently, application of nanosolids, nanofluids,

and nanogases is the promising nanolevel research area of interest for energy savings in heat exchanger. Investigation of the nanolevel heat transfer using molecular dynamics (MD) simulation is only a new pioneer concept for the last few years [1–5]. Here investigations are based on atomic movement within a nanosystem during MD simulation. Experimental study on the heat transfer of nano-size devices or equipments are very time-consuming and expensive for the existing testing capabilities [6–9]. Many studies have been done for enhancing the energy consumption of heat exchanger by improving thermal properties. This thermal behavior of nanocomposite can bring a huge transformation and innovation in the heat transfer. The ultrathin thermal polyurethane heat transfer material can be applied to a number of different fabrics in heat exchanger. The use of graphene (Gr) in thermoplastic polyurethane (TPU), that is, Gr/TPU nanocomposite in place of traditional material in heat exchanger, increases the heat transfer rate in a significant manner. Hussein studied to calculate thermal properties of metals and nonmetals at room temperature for applications in heat exchanger and found that metallic materials are preferably suitable for heat transfer application [10]. But there is some limitation for application of metals due to corrosion. To overcome the situation, implementation of polymer heat exchanger technology for the past decades is a pioneering innovation for heat exchanger materials [11]. The major limitation of polymer for application in heat exchanger is very low thermal conductivity. To improve that property graphene-reinforced polymer nanocomposite is a suitable candidate material in evaporation and condensation applications within heat exchanger [12]. Thermal performance of polymer nanocomposite heat exchangers mainly deals with on shell and tube heat exchangers, plate heat exchangers, finned tube heat exchangers, immersed tube heat exchangers, and hollow fiber heat exchangers [13]. Currently, thermoplastic elastomer is used in heat exchanger applications. Thermoplastics elastomer can be repeatedly softened by heating and solidified by cooling as long as the material is not thermally damaged by overheating [14]. The thermal expansion of thermoplastic polymer can be beneficial with regard to fouling because repeated expansion and contraction of the polymer channels can lead to scale detachment [15]. It is seen from previous studies that there are less number of simulation-based studies like MD simulation on enhanced heat transfer in heat exchanger materials. Detail results of MD simulation are obtained by solving Newton's equation of motion of every atom within nanoscopic system. The basic dynamics parameters of all atoms, that is, position, velocity, and interaction force, play a vital rule during MD simulation. Nanoheat transfer problems of nanocomposites are related to thermo-mechanical properties of nanomaterials. For the design of nanodevices like nanoheat exchanger (NHE), the concepts of the nanothermal properties with temperature variation and dimension of the nanodevice is very much promising ideas. Determination of the thermal properties of nanocomposites by MD simulations is very time-consuming and a challenging task. The objective of the chapter is to characterize thermal properties like thermal conductivity and thermal expansion coefficient of graphene-reinforced polyurethane nanocomposite using molecular dynamics (MD) modeling for nanoheat exchanger material application.

*vi*ð Þ¼ *t* þ *Δt vi*ð Þþ*t*

*DOI: http://dx.doi.org/10.5772/intechopen.86527*

*Molecular Dynamics Simulation-Based Study on Enhancing Thermal Properties…*

the gradient of the total potential energy (*E*) on atom *i*

step.

following equation:

b is the actual bond length.

and θ is the actual bond angle.

ity, and φ<sup>0</sup> is the reference torsional angle.

below:

**91**

*ai*ðÞþ*t ai*ð Þ *t* þ *tΔt* 2

where (*a*i) is the acceleration of atom *i*, *t* is the current time, and Δ*t* is the time

*ai* <sup>¼</sup> *Fi mi*

where *mi* is mass of atom *i*th and *F*<sup>i</sup> is the force vector of atom *i*th obtain from

*Fi* <sup>¼</sup> *dE dri*

In MD simulation, a model system is built at the atomistic level with prescribed potentials (also known as the force field) acting between the atoms. The potential energy is dependent on the force field that is applied to the system. The potential energies of the system are determined from both bonded and nonbonded energies. The total potential energy combines all energetic contributions shown in the

<sup>E</sup> <sup>¼</sup> Ebond <sup>þ</sup> Eangle <sup>þ</sup> Etorsion <sup>þ</sup> Einversion <sup>þ</sup> Ecross‐term <sup>þ</sup> ð Þ EvdW <sup>þ</sup> Ecoulomb (7)

The different types of bond and nonbond energy equations have been shown

*kb*ð Þ *b* � *b*<sup>0</sup>

*kθ*ð Þ *θ* � *θ*<sup>0</sup>

where kb is the stretching force constant, b0 is the equilibrium bond length, and

where k<sup>θ</sup> is the angle-bending force constant, θ<sup>0</sup> is the equilibrium bond angle,

where kφ is the torsional barrier, φ is the actual torsion angle, n is the periodic-

where K<sup>ω</sup> is the force constant and ω is the angle between the axis and the plane.

*r*12 *ij*

*kω*ð Þ *ω* � *ω*<sup>0</sup>

� *Bij r*6 *ij*

2

*Evdw* <sup>¼</sup> <sup>∑</sup> *Aij*

*Ebond* <sup>¼</sup> <sup>1</sup> 2

*Ebend* <sup>¼</sup> <sup>1</sup> 2

2

*Einversion* <sup>¼</sup> <sup>1</sup>

*Etorsion* <sup>¼</sup> <sup>1</sup>

<sup>E</sup> <sup>¼</sup> Ebonded <sup>þ</sup> Enon‐bonded (5)

<sup>2</sup> (8)

<sup>2</sup> (9)

<sup>2</sup> (11)

(12)

*k<sup>ϕ</sup>* 1 þ cos *nϕ* � *ϕ*<sup>0</sup> ð Þ ð Þ (10)

<sup>E</sup> <sup>¼</sup> <sup>ð</sup>Evalance <sup>þ</sup> Ecross‐termÞ þ Enon‐bonded (6)

*Δt* (2)

(3)

(4)

MD simulation consists of many parts which are mainly (a) molecular interactions, (b) molecular minimization, (c) algorithms, (d) ensemble, (e) boundary conditions, (f) atomistic stress calculation, etc. MD simulation helps to determine the position (ri) and velocity (vi) vectors of atom *i* with the time integration method (e.g., the velocity Verlet algorithm):

$$r\_i(t + \Delta t) = r\_i(t) + \nu\_i(t)\Delta t + \frac{1}{2}a\_i(t)\Delta t^2 \tag{1}$$

*Molecular Dynamics Simulation-Based Study on Enhancing Thermal Properties… DOI: http://dx.doi.org/10.5772/intechopen.86527*

$$v\_i(t + \Delta t) = v\_i(t) + \frac{a\_i(t) + a\_i(t + t\Delta t)}{2}\Delta t \tag{2}$$

where (*a*i) is the acceleration of atom *i*, *t* is the current time, and Δ*t* is the time step.

$$a\_i = \frac{F\_i}{m\_i} \tag{3}$$

where *mi* is mass of atom *i*th and *F*<sup>i</sup> is the force vector of atom *i*th obtain from the gradient of the total potential energy (*E*) on atom *i*

$$F\_i = \frac{dE}{dr\_i} \tag{4}$$

In MD simulation, a model system is built at the atomistic level with prescribed potentials (also known as the force field) acting between the atoms. The potential energy is dependent on the force field that is applied to the system. The potential energies of the system are determined from both bonded and nonbonded energies. The total potential energy combines all energetic contributions shown in the following equation:

$$\mathbf{E} = \mathbf{E}\_{\text{bonded}} + \mathbf{E}\_{\text{non-bonded}} \tag{5}$$

$$\mathbf{E} = (\mathbf{E}\_{\text{valence}} + \mathbf{E}\_{\text{cross-term}}) + \mathbf{E}\_{\text{non-bonded}} \tag{6}$$

$$\mathbf{E} = \left(\mathbf{E}\_{\text{bond}} + \mathbf{E}\_{\text{angle}} + \mathbf{E}\_{\text{torsion}} + \mathbf{E}\_{\text{inversion}}\right) + \mathbf{E}\_{\text{cross-term}} + \left(\mathbf{E}\_{\text{vdW}} + \mathbf{E}\_{\text{coulomb}}\right) \tag{7}$$

The different types of bond and nonbond energy equations have been shown below:

$$E\_{bond} = \frac{1}{2}k\_b(b - b\_0)^2\tag{8}$$

where kb is the stretching force constant, b0 is the equilibrium bond length, and b is the actual bond length.

$$E\_{bend} = \frac{1}{2}k\_{\theta}(\theta - \theta\_0)^2\tag{9}$$

where k<sup>θ</sup> is the angle-bending force constant, θ<sup>0</sup> is the equilibrium bond angle, and θ is the actual bond angle.

$$E\_{\text{torsion}} = \frac{1}{2}k\_{\phi}(\mathbf{1} + \cos\left(n\phi - \phi\_0\right))\tag{10}$$

where kφ is the torsional barrier, φ is the actual torsion angle, n is the periodicity, and φ<sup>0</sup> is the reference torsional angle.

$$E\_{inversion} = \frac{1}{2}k\_w(\rho - \rho\_0)^2\tag{11}$$

where K<sup>ω</sup> is the force constant and ω is the angle between the axis and the plane.

$$E\_{vdw} = \sum \frac{A\_{\vec{ij}}}{r\_{\vec{ij}}^{12}} - \frac{B\_{\vec{ij}}}{r\_{\vec{ij}}^6} \tag{12}$$

and nanogases is the promising nanolevel research area of interest for energy savings in heat exchanger. Investigation of the nanolevel heat transfer using molecular dynamics (MD) simulation is only a new pioneer concept for the last few years [1–5]. Here investigations are based on atomic movement within a nanosystem during MD simulation. Experimental study on the heat transfer of nano-size devices or equipments are very time-consuming and expensive for the existing testing capabilities [6–9]. Many studies have been done for enhancing the energy consumption of heat exchanger by improving thermal properties. This thermal behavior of nanocomposite can bring a huge transformation and innovation in the heat transfer. The ultrathin thermal polyurethane heat transfer material can be applied to a number of different fabrics in heat exchanger. The use of graphene (Gr) in thermoplastic polyurethane (TPU), that is, Gr/TPU nanocomposite in place of traditional material in heat exchanger, increases the heat transfer rate in a significant manner. Hussein studied to calculate thermal properties of metals and nonmetals at room temperature for applications in heat exchanger and found that metallic materials are preferably suitable for heat transfer application [10]. But there is some limitation for application of metals due to corrosion. To overcome the situation, implementation of polymer heat exchanger technology for the past decades is a pioneering innovation for heat exchanger materials [11]. The major limitation of polymer for application in heat exchanger is very low thermal conductivity. To improve that property graphene-reinforced polymer nanocomposite is a suitable candidate material in evaporation and condensation applications within heat exchanger [12]. Thermal performance of polymer nanocomposite heat exchangers mainly deals with on shell and tube heat exchangers, plate heat exchangers, finned tube heat exchangers, immersed tube heat exchangers, and hollow fiber heat exchangers [13]. Currently, thermoplastic elastomer is used in heat exchanger applications. Thermoplastics elastomer can be repeatedly softened by heating and solidified by cooling as long as the material is not thermally damaged by overheating [14]. The thermal expansion of thermoplastic polymer can be beneficial with regard to fouling because repeated expansion and contraction of the polymer channels can lead to scale detachment [15]. It is seen from previous studies that there are less number of simulation-based studies like MD simulation on enhanced heat transfer in heat exchanger materials. Detail results of MD simulation

*Inverse Heat Conduction and Heat Exchangers*

are obtained by solving Newton's equation of motion of every atom within nanoscopic system. The basic dynamics parameters of all atoms, that is, position, velocity, and interaction force, play a vital rule during MD simulation. Nanoheat transfer problems of nanocomposites are related to thermo-mechanical properties of nanomaterials. For the design of nanodevices like nanoheat exchanger (NHE), the concepts of the nanothermal properties with temperature variation and dimension of the nanodevice is very much promising ideas. Determination of the thermal properties of nanocomposites by MD simulations is very time-consuming and a challenging task. The objective of the chapter is to characterize thermal properties like thermal conductivity and thermal expansion coefficient of graphene-reinforced

polyurethane nanocomposite using molecular dynamics (MD) modeling for

*ri*ð Þ¼ *t* þ *Δt ri*ð Þþ*t vi*ð Þ*t Δt* þ

MD simulation consists of many parts which are mainly (a) molecular interactions, (b) molecular minimization, (c) algorithms, (d) ensemble, (e) boundary conditions, (f) atomistic stress calculation, etc. MD simulation helps to determine the position (ri) and velocity (vi) vectors of atom *i* with the time integration method

> 1 2 *ai*ð Þ*t Δt*

<sup>2</sup> (1)

nanoheat exchanger material application.

(e.g., the velocity Verlet algorithm):

**90**

where *Aij* and *Bij* are the repulsive and attractive term coefficients, respectively, and *rij* is the distance between the two atoms.

$$E\_{coulomb} = \frac{1}{\varepsilon} \frac{q\_1 q\_2}{r\_{ij}} \tag{13}$$

graphene-reinforced nanocomposites is considered in this study. The constructed nanocomposite lattice parameters with dimensions a = 22 Å, b = 22 Å, and c = 95 Å are constrained in such a way that after dynamic equilibrium process density of Gr/ TPU, nanocomposite is 1.38 g/cc (nearly experimental value). MD simulation run is divided into two parts, namely equilibration run and production run. Equilibration run helps to develop the molecular structure under the condition of the desired thermodynamic state, while the production run helps to calculate different thermal parameters, namely specific heat, thermal expansion, and thermal conductivity. In equilibration run, two important conditions have to be fulfilled. One is the minimum energy stabilized condition at a prescribed temperature, and another is initial stress-free structure within periodic boundary condition. So, first steps nanocomposite models are moved through energy minimization, canonical

*Molecular Dynamics Simulation-Based Study on Enhancing Thermal Properties…*

ensemble (constant number of atoms, volume, and temperature) (NVT)) dynamic simulations, and temperature annealing cycle, respectively. The duration of the dynamic run is considered 200 ps with an integration time step of 1 fs (femtosecond). This process is followed by graphene as a rigid structure so that the lattice dimension (c) in the z-direction will be changed. Temperature annealing cycle involves temperature up and down from 300 to 600 K to get the minimization of energy in the structure. The annealing time is set for 500 ps, during which the temperature is raised from 300 to 600 K with a rate of 6 K/ps and cool down to 300 K with the same rate. In the second step, the non-constrained parts (TPU)

nanocomposite will be 1.38 g/cc (nearly experimental value) after using isothermalisobaric (NPT) ensembles. The lowest energy structure models are fully relaxed under an isothermal-isobaric NPT ensemble (i.e., constant numbers of atoms, pressure, and temperature) at 300 K and 1 atm for 500 ps. The isothermal-isobaric (NPT) ensembles help to relax the lattices parameters and angles in order to obtain a final reasonable equilibrated structure. These steps generate various curves of various parameters such as energy, pressure, volume, and temperature versus simulation run time. These curves are very important to study thermal properties of the nanocomposite. **Figure 1** shows the developed 1% Gr/TPU nanocomposite model by

After dynamic equilibrium process density of 1% Gr/TPU, the nanocomposite is

Heat capacity is one of the important thermal properties for the nanocomposite system. In this work, MD simulation is applied to calculate the isobaric heat capacity (Cp), and the value of Cp can be determined according to the following equation:

within lattices are compressed in such a way so that the final density of

MD simulation.

**93**

**Figure 1.**

*Developed 1% Gr/TPU nanocomposite model.*

*DOI: http://dx.doi.org/10.5772/intechopen.86527*

1.38 g/cc which is shown in **Figure 2**.

where q1 and q2 are the charges on the interacting atoms, ε is the dielectric constant, and rij is the interatomic distance.

The Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies (COMPASS) [16] which is incorporated in both amorphous and forcite plus atomistic simulation modules in the Material Studio is used for this present study. COMPASS functional form has 11 valence terms (including diagonal and offdiagonal cross coupling terms) and 2 nonbond interaction terms (the Coulombic and Lennard-Jones functions for electrostatic and van der Waals (vdW) interactions, respectively). During all the simulations, the temperature and pressure are maintained by Andersen and Berendsen method. The calculation of nonbonded interactions is simulated by applying a cutoff distance of 12.5 Å. The spline and buffer widths are 1 and 0.5 Å, respectively.

Experimental methods for prediction of the enhanced thermal properties of graphene-reinforced nanocomposites are limited because nanometer scale measurements are difficult and costly. Thus, MD simulation techniques are only an economical path to characterize nanomaterial and nanocomposites for heat exchanger material within small length and small time.

### **2. MD simulation models and methods for thermal property calculation**

Heat can transfer through electrons and phonons, by excitations and by scattering in the nanocomposites [17]. This different ways heat transfer methods help to understand the mechanism of thermal properties enhancement in graphene-based nanocomposites. In order to study the enhanced thermal properties of graphenereinforced thermoplastic polyurethane nanocomposites, the wide range thermal parameters of the nanocomposite material (like thermal conductivity, coefficient of the thermal expansion, glass transition temperature, etc.) have to be calculated. To calculate these parameters, different simulation models are constructed using molecular modeling software package by Materials Studio 2017. In MD simulation study, there are three types of models which are developed, namely concentrated model (CM), layer model (LM), and interfacial model (IM). Different models have been used to study different properties of the nanocomposite with respect to different parameters. In this study, we have focused mainly on layer models to characterize enhanced thermal properties. Both graphene and polyurethane models are simulated separately before constructing the graphene-reinforced thermoplastic polyurethane (Gr/TPU) nanocomposites. The condensed-phase optimized molecular potential for atomistic simulation studies (COMPASS) force field has been selected to describe the atomistic behavior for the simulation models. In MD simulation, each atom is modeled as a point mass and interacts with other atoms through force field. The position and momentum of atoms are updated based on Newton's equation of motion.

After the construction all of the graphene and TPU model using amorphous cell module within the Material Studio, build layer option is used to construct the layer models of the Gr/TPU nanocomposites. It is seen that with high weight fraction condition of graphene, nanocomposites behave likely more brittle than polymer matrix due to growth of void and chain disentanglement. So, 1% weight fraction of *Molecular Dynamics Simulation-Based Study on Enhancing Thermal Properties… DOI: http://dx.doi.org/10.5772/intechopen.86527*

**Figure 1.**

*Developed 1% Gr/TPU nanocomposite model.*

where *Aij* and *Bij* are the repulsive and attractive term coefficients, respectively,

*ε q*1*q*<sup>2</sup> *rij*

The Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies (COMPASS) [16] which is incorporated in both amorphous and forcite plus atomistic simulation modules in the Material Studio is used for this present study. COMPASS functional form has 11 valence terms (including diagonal and offdiagonal cross coupling terms) and 2 nonbond interaction terms (the Coulombic and Lennard-Jones functions for electrostatic and van der Waals (vdW) interactions, respectively). During all the simulations, the temperature and pressure are maintained by Andersen and Berendsen method. The calculation of nonbonded interactions is simulated by applying a cutoff distance of 12.5 Å. The spline and

(13)

*Ecoulomb* <sup>¼</sup> <sup>1</sup>

where q1 and q2 are the charges on the interacting atoms, ε is the dielectric

Experimental methods for prediction of the enhanced thermal properties of graphene-reinforced nanocomposites are limited because nanometer scale measurements are difficult and costly. Thus, MD simulation techniques are only an economical path to characterize nanomaterial and nanocomposites for heat

**2. MD simulation models and methods for thermal property calculation**

Heat can transfer through electrons and phonons, by excitations and by scattering in the nanocomposites [17]. This different ways heat transfer methods help to understand the mechanism of thermal properties enhancement in graphene-based nanocomposites. In order to study the enhanced thermal properties of graphenereinforced thermoplastic polyurethane nanocomposites, the wide range thermal parameters of the nanocomposite material (like thermal conductivity, coefficient of the thermal expansion, glass transition temperature, etc.) have to be calculated. To calculate these parameters, different simulation models are constructed using molecular modeling software package by Materials Studio 2017. In MD simulation study, there are three types of models which are developed, namely concentrated model (CM), layer model (LM), and interfacial model (IM). Different models have been used to study different properties of the nanocomposite with respect to different parameters. In this study, we have focused mainly on layer models to characterize enhanced thermal properties. Both graphene and polyurethane models are simulated separately before constructing the graphene-reinforced thermoplastic polyurethane (Gr/TPU) nanocomposites. The condensed-phase optimized molecular potential for atomistic simulation studies (COMPASS) force field has been selected to describe the atomistic behavior for the simulation models. In MD simulation, each atom is modeled as a point mass and interacts with other atoms through force field. The position and momentum of atoms are updated based on Newton's

After the construction all of the graphene and TPU model using amorphous cell module within the Material Studio, build layer option is used to construct the layer models of the Gr/TPU nanocomposites. It is seen that with high weight fraction condition of graphene, nanocomposites behave likely more brittle than polymer matrix due to growth of void and chain disentanglement. So, 1% weight fraction of

and *rij* is the distance between the two atoms.

*Inverse Heat Conduction and Heat Exchangers*

constant, and rij is the interatomic distance.

buffer widths are 1 and 0.5 Å, respectively.

equation of motion.

**92**

exchanger material within small length and small time.

graphene-reinforced nanocomposites is considered in this study. The constructed nanocomposite lattice parameters with dimensions a = 22 Å, b = 22 Å, and c = 95 Å are constrained in such a way that after dynamic equilibrium process density of Gr/ TPU, nanocomposite is 1.38 g/cc (nearly experimental value). MD simulation run is divided into two parts, namely equilibration run and production run. Equilibration run helps to develop the molecular structure under the condition of the desired thermodynamic state, while the production run helps to calculate different thermal parameters, namely specific heat, thermal expansion, and thermal conductivity. In equilibration run, two important conditions have to be fulfilled. One is the minimum energy stabilized condition at a prescribed temperature, and another is initial stress-free structure within periodic boundary condition. So, first steps nanocomposite models are moved through energy minimization, canonical ensemble (constant number of atoms, volume, and temperature) (NVT)) dynamic simulations, and temperature annealing cycle, respectively. The duration of the dynamic run is considered 200 ps with an integration time step of 1 fs (femtosecond). This process is followed by graphene as a rigid structure so that the lattice dimension (c) in the z-direction will be changed. Temperature annealing cycle involves temperature up and down from 300 to 600 K to get the minimization of energy in the structure. The annealing time is set for 500 ps, during which the temperature is raised from 300 to 600 K with a rate of 6 K/ps and cool down to 300 K with the same rate. In the second step, the non-constrained parts (TPU) within lattices are compressed in such a way so that the final density of nanocomposite will be 1.38 g/cc (nearly experimental value) after using isothermalisobaric (NPT) ensembles. The lowest energy structure models are fully relaxed under an isothermal-isobaric NPT ensemble (i.e., constant numbers of atoms, pressure, and temperature) at 300 K and 1 atm for 500 ps. The isothermal-isobaric (NPT) ensembles help to relax the lattices parameters and angles in order to obtain a final reasonable equilibrated structure. These steps generate various curves of various parameters such as energy, pressure, volume, and temperature versus simulation run time. These curves are very important to study thermal properties of the nanocomposite. **Figure 1** shows the developed 1% Gr/TPU nanocomposite model by MD simulation.

After dynamic equilibrium process density of 1% Gr/TPU, the nanocomposite is 1.38 g/cc which is shown in **Figure 2**.

Heat capacity is one of the important thermal properties for the nanocomposite system. In this work, MD simulation is applied to calculate the isobaric heat capacity (Cp), and the value of Cp can be determined according to the following equation:

$$\text{Cp} = \frac{(\text{KE} + \text{PE} + \text{PV})^2}{\text{K}\_\text{B} \text{ T}} \tag{14}$$

where V0 is the equilibrated system volume before the cooling simulation starts. The fractional change of volume with respect to temperature (ΔV/ΔT) is obtainable from volume versus temperature relationship curve. It is seen that change of volume with respect to temperature (ΔV/ΔT) is a different value above glass transition

nanocomposite. So, volumetric coefficient of thermal expansion (VCTE) has two values due to glass transition temperature (Tg) [20]. The glass transition temperature and volumetric coefficient of thermal expansion (VCTE) of nanocomposite can also be obtained by calculating the free volume as a function of temperature, since the free volume undergoes an abrupt change when the material goes through glass transition. By probing the lattice cell with a spherical probe, using the "atom volume and surfaces" tool of the Materials Studio (MS), the free volume in the nanocomposite is calculated as a function of temperature during the annealing

Free volume is the volume that is not occupied by either the graphene or the TPU chains. The free volume fraction (FVF) can be obtained by the following

*FVF* <sup>¼</sup> *Vf*

where Vf is the free volume and V0 is the occupied volume of the polymer

contribution. Therefore, total thermal conductivity (K)

*Free volume calculation using atom volume and surfaces tool in MS.*

Thermal conductivity is the sum of the phonon contribution and the electronic

where Kp and Ke are the phonon contribution and the electronic contribution to

the thermal conductivity, respectively. Though, electrons contributed thermal conductivity is neglected in most graphene-reinforced nanocomposite. Thermal conductivity is an important thermal property relevant to thermal management applications. Thermal conductivity is generally calculated using equilibrium or nonequilibrium MD approaches. Equilibrium molecular dynamic (EMD) facilitates thermal conductivity prediction in all directions using one simulation, whereas nonequilibrium molecular dynamic (NEMD) requires the use of a thermal

gradient and therefore only enables the calculation of thermal conductivity in one

*Vf* þ *Vo*

*K* ¼ *Ke* þ *Kp* (18)

(17)

temperature (Tg) for graphene-reinforced thermoplastic polyurethane

*Molecular Dynamics Simulation-Based Study on Enhancing Thermal Properties…*

process as shown in **Figure 3**.

*DOI: http://dx.doi.org/10.5772/intechopen.86527*

equation [21]:

chains.

**Figure 3.**

**95**

where KE is the kinetic energy, PE is the potential energy, P is the pressure, V is the volume, KB is the Boltzmann constant, and T is the temperature. The specific heat at constant volume (Cv) is obtained from the following equation:

$$\mathbf{C}\mathbf{v} = \frac{\left(\delta E^2\right)\_{\mathrm{NVT}}}{\mathbf{k}\_{\mathrm{B}}\mathbf{T}^2} \tag{15}$$

where δE is the fluctuation of the energy, kB and T are Boltzmann constant, and absolute temperature, respectively.

In order to study the glass transition temperature and coefficient of thermal expansion (CTE), a high-temperature annealing protocol is followed. At each temperature, the system is equilibrated by isothermal-isobaric (NPT) ensemble in MD simulation at atmospheric pressure for 500 ps. The temperature is raised up to 600 K and equilibrated for 500 ps using an NPT ensemble under atmospheric pressure and then dropped by 20 K each time until it reached 300 K. The cooling down method is applied after the heating up method by decreasing the temperature with the same settings and simulation time. Since each temperature drop is only 20 K, the structure is re-equilibrated very quickly every time its temperature is decreased. For each temperature, the volume of the simulation box V is examined over the duration time of the MD simulation, and the average value is calculated. From the volume versus temperature relationship curve, it is seen that there is a discontinuity in the volume versus temperature slope, which gives the glass transition temperature (Tg) of nanocomposite. The volume versus temperature (V-T) results are important to know two factors; first, this provides a means of determining the quality of the force field used in the simulations, and second, prediction of the volumetric glass transition temperature (Tg) [18]. The simulation result is in good agreement with experiment and demonstrates the accuracy of COMPASS force field. The volumetric coefficient of thermal expansion (VCTE) is defined by (α) [19]:

$$a = \frac{1}{V\_0} \frac{\Delta V}{\Delta T} \tag{16}$$

**Figure 2.** *Density (g/cc) versus time (ps) in dynamic equilibrium run.*

*Molecular Dynamics Simulation-Based Study on Enhancing Thermal Properties… DOI: http://dx.doi.org/10.5772/intechopen.86527*

where V0 is the equilibrated system volume before the cooling simulation starts. The fractional change of volume with respect to temperature (ΔV/ΔT) is obtainable from volume versus temperature relationship curve. It is seen that change of volume with respect to temperature (ΔV/ΔT) is a different value above glass transition temperature (Tg) for graphene-reinforced thermoplastic polyurethane nanocomposite. So, volumetric coefficient of thermal expansion (VCTE) has two values due to glass transition temperature (Tg) [20]. The glass transition temperature and volumetric coefficient of thermal expansion (VCTE) of nanocomposite can also be obtained by calculating the free volume as a function of temperature, since the free volume undergoes an abrupt change when the material goes through glass transition. By probing the lattice cell with a spherical probe, using the "atom volume and surfaces" tool of the Materials Studio (MS), the free volume in the nanocomposite is calculated as a function of temperature during the annealing process as shown in **Figure 3**.

Free volume is the volume that is not occupied by either the graphene or the TPU chains. The free volume fraction (FVF) can be obtained by the following equation [21]:

$$FVF = \frac{V\_f}{V\_f + V\_o} \tag{17}$$

where Vf is the free volume and V0 is the occupied volume of the polymer chains.

Thermal conductivity is the sum of the phonon contribution and the electronic contribution. Therefore, total thermal conductivity (K)

$$K = K\_{\varepsilon} + K\_{p} \tag{18}$$

where Kp and Ke are the phonon contribution and the electronic contribution to the thermal conductivity, respectively. Though, electrons contributed thermal conductivity is neglected in most graphene-reinforced nanocomposite. Thermal conductivity is an important thermal property relevant to thermal management applications. Thermal conductivity is generally calculated using equilibrium or nonequilibrium MD approaches. Equilibrium molecular dynamic (EMD) facilitates thermal conductivity prediction in all directions using one simulation, whereas nonequilibrium molecular dynamic (NEMD) requires the use of a thermal gradient and therefore only enables the calculation of thermal conductivity in one

**Figure 3.** *Free volume calculation using atom volume and surfaces tool in MS.*

Cp <sup>¼</sup> ð Þ KE <sup>þ</sup> PE <sup>þ</sup> PV <sup>2</sup>

Cv <sup>¼</sup> *<sup>δ</sup>E*<sup>2</sup>

heat at constant volume (Cv) is obtained from the following equation:

absolute temperature, respectively.

*Inverse Heat Conduction and Heat Exchangers*

(α) [19]:

**Figure 2.**

**94**

*Density (g/cc) versus time (ps) in dynamic equilibrium run.*

where KE is the kinetic energy, PE is the potential energy, P is the pressure, V is the volume, KB is the Boltzmann constant, and T is the temperature. The specific

where δE is the fluctuation of the energy, kB and T are Boltzmann constant, and

In order to study the glass transition temperature and coefficient of thermal expansion (CTE), a high-temperature annealing protocol is followed. At each temperature, the system is equilibrated by isothermal-isobaric (NPT) ensemble in MD simulation at atmospheric pressure for 500 ps. The temperature is raised up to 600 K and equilibrated for 500 ps using an NPT ensemble under atmospheric pressure and then dropped by 20 K each time until it reached 300 K. The cooling down method is applied after the heating up method by decreasing the temperature with the same settings and simulation time. Since each temperature drop is only 20 K, the structure is re-equilibrated very quickly every time its temperature is decreased. For each temperature, the volume of the simulation box V is examined over the duration time of the MD simulation, and the average value is calculated. From the volume versus temperature relationship curve, it is seen that there is a discontinuity in the volume versus temperature slope, which gives the glass transition temperature (Tg) of nanocomposite. The volume versus temperature (V-T) results are important to know two factors; first, this provides a means of determining the quality of the force field used in the simulations, and second, prediction of the volumetric glass transition temperature (Tg) [18]. The simulation result is in good agreement with experiment and demonstrates the accuracy of COMPASS force field. The volumetric coefficient of thermal expansion (VCTE) is defined by

> *<sup>α</sup>* <sup>¼</sup> <sup>1</sup> *V*<sup>0</sup>

*ΔV*

NVT

KB <sup>T</sup> (14)

kBT<sup>2</sup> (15)

*<sup>Δ</sup><sup>T</sup>* (16)

direction. The indirect method is an equilibrium molecular dynamic (EMD) method which is derived from Green-Kubo approach [22, 23], where current fluctuations are used to compute the thermal conductivity via the fluctuation-dissipation theorem [24]:

$$K\_{a\beta} = \frac{V}{k\_B T^2} \int\_0^\infty \langle J\_a(t) J\_\beta(\mathbf{0}) \rangle dt \tag{19}$$

*Jq* <sup>¼</sup> *<sup>E</sup>*

*Molecular Dynamics Simulation-Based Study on Enhancing Thermal Properties…*

where E is the subtracted energy from the heat sink. Ac is cross-sectional area

Interface thermal conductance (ITC) is developed between graphene and polymer matrix within nanocomposite due to their weak interactions. However, the thermal conductivity (TC) of graphene-reinforced TPU nanocomposites are far below the thermal conductivity (TC) of graphene. Such a low filler efficiency is most likely attributed to the low interface thermal conductance (ITC) between graphene and polymer. The simulating method (pump-probe transient thermoreflectance method) is used to calculate the interface thermal conductance (ITC) between the filler (graphene) and matrix (TPU). The graphene-TPU interfacial

> *<sup>R</sup>* <sup>¼</sup> *<sup>Δ</sup><sup>T</sup> Jq*

In the equilibration stage within MD simulation, NPT dynamic run is carried out for 500 ps at room temperature and pressure to generate curves of energy, density, pressure, and temperature versus time. These curves are used to decide the cutoff between equilibration and production runs. It is observed that all the models are equilibrated at 50 ps, that is, no fluctuations after 50 ps. At the end of the equilibrations, the density of the nanocomposites stabilized at an average density of 1.3 g/cc with a standard deviation of 0.02 g/cc. The reason behind for the density differs from experimental value because MD simulation deals with material

defect-free and impurities. It is seen in the volume versus temperature (V-T) curve during NPT dynamic run that free volume change affects the thermal properties of the nanocomposites. Further it is also observed that free volume within lattice increases according to the strain application. The glass transition mainly depends on two factors: (a) free volume and (b) the mobility of chain segments. Initially there is no graphene in TPU chains for that free volume is zero in the simulated cell and polymers are free to move within this cell volume. As a result, high values for the entropy and low values for the bulk, shear, and Young's moduli. Further, the incorporation of the graphene into the simulated cells increases free volume and decrease the entropy and increase the values of the bulk, shear, and Young's moduli. Since the entropy of the system is related to the free volume and the Connolly surface of the TPU nanocomposites, the prediction of these parameters may give a concept on the enhancement of the nanocomposite thermal properties. The simulation results show that the TPU matrix reinforced with graphene have a tendency to increase glass transition temperature (Tg) for stronger interlocking between the graphene and TPU molecules. The Connolly surface to volume ratio is small values for the neat TPU and increase with the increase in the graphene loading. The free volume is defined as the volume on the side of the Connolly surface without atoms. This simulation study reveals that Tg of TPU is in the range of 285–305°C and for the Gr/TPU nanocomposite is 350 K (experimental 223–323 K). From the V-T curve, volumetric coefficient of thermal expansion (α) is evaluated from the slope

above and below Tg, which is between 2.6 � <sup>10</sup>�<sup>5</sup> and 2.4 � <sup>10</sup>�<sup>4</sup>

). The free volume change rapidly when the material goes

through glass transition, according to Fox and Flory's theory of glass transition [20].

thermal resistance (R) is then calculated by following equation:

and Δt are the time step.

*DOI: http://dx.doi.org/10.5772/intechopen.86527*

**3. Results and discussion**

tal 3.15 � <sup>10</sup>�<sup>4</sup>

**97**

/K�<sup>1</sup>

<sup>2</sup>*AcΔ<sup>t</sup>* (22)

(23)

/K�<sup>1</sup> (experimen-

where kB is the Boltzmann constant, V and T denote the volume and temperature of the system, J<sup>α</sup> is the heat flux in the *α* direction, and the angular brackets denote the ensemble average. The heat flux vector can be written as

$$\mathbf{J(t)} = \frac{\mathbf{1}}{V} \frac{d}{dt} \sum\_{i=1}^{N} r\_i E\_i \tag{20}$$

where ri and Ei are the position and total energy of the ith atom, respectively. EMD is particularly useful for geometries where periodic boundary conditions can be applied. EMD is often computationally more expensive, and the results are more sensitive to the simulation parameters. In EMD, the system is set to the desired temperature, and then a constant energy scheme is used with the wellknown Green-Kubo relations to calculate the thermal conductivity tensor. The direct method is a NEMD method in which a temperature difference is introduced into the simulation domain and the thermal conductivity is computed according to Fourier's law as K = �J/ΔT, where J and ΔT are heat flux and temperature gradient across the system, respectively. For nonequilibrium MD (NEMD) methods, a long slab of polymer nanocomposite is constructed, and a difference in temperature is established between a heat source and a sink at the ends of the slab, and the flux is calculated. Equilibrium systems are simulated by nonequilibrium MD (NEMD) based on our homemade PERL script to compute their thermal conductivities. The nonequilibrium state can be established either by applying two thermostats at different temperatures to maintain a constant temperature at the two ends of the system or by artificially swapping atom velocities in different regions to impose a constant heat flux also known as the reverse nonequilibrium MD (RNEMD) method based on Muller-Plathe's approach [25]. In the reverse nonequilibrium MD (RNEMD) method, the energy exchange is carried out by exchanging the kinetic energy of two particles: the hottest particle in the cold layer and the coldest particle in the hot layer. The energy E is therefore variable and needs averaging over many exchanges. In the RNEMD method, the simulation box was divided into a number of slabs in the heat flux (z) direction with the same thickness. The heat flux was generated by exchanging the kinetic energy between the highest kinetic (the hottest) atom in the heat sink and the lowest kinetic energy (the coldest) atom in the heat source. The larger momentum exchange rate in RNEMD method suggests higher energy exchange frequency between heat source and heat sink. The thermal conductivity was calculated using Fourier's law:

$$K = -\frac{J\_q}{\nabla T} \tag{21}$$

where ∇T denoted the time-integrated temperature gradient from least squares approximation of the discrete temperatures according to the heat flow direction. The temperature of each slab was calculated using the virial theorem, and Jq is the heat flux given as

*Molecular Dynamics Simulation-Based Study on Enhancing Thermal Properties… DOI: http://dx.doi.org/10.5772/intechopen.86527*

$$J\_q = \frac{E}{2A\_c\Delta t} \tag{22}$$

where E is the subtracted energy from the heat sink. Ac is cross-sectional area and Δt are the time step.

Interface thermal conductance (ITC) is developed between graphene and polymer matrix within nanocomposite due to their weak interactions. However, the thermal conductivity (TC) of graphene-reinforced TPU nanocomposites are far below the thermal conductivity (TC) of graphene. Such a low filler efficiency is most likely attributed to the low interface thermal conductance (ITC) between graphene and polymer. The simulating method (pump-probe transient thermoreflectance method) is used to calculate the interface thermal conductance (ITC) between the filler (graphene) and matrix (TPU). The graphene-TPU interfacial thermal resistance (R) is then calculated by following equation:

$$R = \frac{\Delta T}{J\_q} \tag{23}$$

#### **3. Results and discussion**

direction. The indirect method is an equilibrium molecular dynamic (EMD) method which is derived from Green-Kubo approach [22, 23], where current fluctuations are used to compute the thermal conductivity via the fluctuation-dissipation

> ð<sup>∞</sup> 0

where kB is the Boltzmann constant, V and T denote the volume and temperature of the system, J<sup>α</sup> is the heat flux in the *α* direction, and the angular brackets

where ri and Ei are the position and total energy of the ith atom, respectively. EMD is particularly useful for geometries where periodic boundary conditions can be applied. EMD is often computationally more expensive, and the results are more sensitive to the simulation parameters. In EMD, the system is set to the desired temperature, and then a constant energy scheme is used with the wellknown Green-Kubo relations to calculate the thermal conductivity tensor. The direct method is a NEMD method in which a temperature difference is introduced into the simulation domain and the thermal conductivity is computed according to Fourier's law as K = �J/ΔT, where J and ΔT are heat flux and temperature gradient across the system, respectively. For nonequilibrium MD (NEMD) methods, a long slab of polymer nanocomposite is constructed, and a difference in temperature is established between a heat source and a sink at the ends of the slab, and the flux is calculated. Equilibrium systems are simulated by nonequilibrium MD (NEMD) based on our homemade PERL script to compute their thermal conductivities. The nonequilibrium state can be established either by applying two thermostats at different temperatures to maintain a constant temperature at the two ends of the system or by artificially swapping atom velocities in different regions to impose a constant heat flux also known as the reverse nonequilibrium MD (RNEMD) method based on Muller-Plathe's approach [25]. In the reverse nonequilibrium MD (RNEMD) method, the energy exchange is carried out by exchanging the kinetic energy of two particles: the hottest particle in the cold layer and the coldest particle in the hot layer. The energy E is therefore variable and needs averaging over many exchanges. In the RNEMD method, the simulation box was divided into a number of slabs in the heat flux (z) direction with the same thickness. The heat flux was generated by exchanging the kinetic energy between the highest kinetic (the hottest) atom in the heat sink and the lowest kinetic energy (the coldest) atom in the heat source. The larger momentum exchange rate in RNEMD method suggests higher energy exchange frequency between heat source and heat sink. The thermal conductivity was calculated using

*<sup>K</sup>* ¼ � *Jq*

where ∇T denoted the time-integrated temperature gradient from least squares approximation of the discrete temperatures according to the heat flow direction. The temperature of each slab was calculated using the virial theorem, and Jq is the

*<sup>J</sup>α*ð Þ*<sup>t</sup> <sup>J</sup>β*ð Þi <sup>0</sup> *dt* � (19)

*riEi* (20)

<sup>∇</sup>*<sup>T</sup>* (21)

*<sup>K</sup>αβ* <sup>¼</sup> *<sup>V</sup> kBT*<sup>2</sup>

denote the ensemble average. The heat flux vector can be written as

J tðÞ¼ <sup>1</sup> *V d dt* <sup>∑</sup> *N i*¼1

theorem [24]:

*Inverse Heat Conduction and Heat Exchangers*

Fourier's law:

heat flux given as

**96**

In the equilibration stage within MD simulation, NPT dynamic run is carried out for 500 ps at room temperature and pressure to generate curves of energy, density, pressure, and temperature versus time. These curves are used to decide the cutoff between equilibration and production runs. It is observed that all the models are equilibrated at 50 ps, that is, no fluctuations after 50 ps. At the end of the equilibrations, the density of the nanocomposites stabilized at an average density of 1.3 g/cc with a standard deviation of 0.02 g/cc. The reason behind for the density differs from experimental value because MD simulation deals with material defect-free and impurities. It is seen in the volume versus temperature (V-T) curve during NPT dynamic run that free volume change affects the thermal properties of the nanocomposites. Further it is also observed that free volume within lattice increases according to the strain application. The glass transition mainly depends on two factors: (a) free volume and (b) the mobility of chain segments. Initially there is no graphene in TPU chains for that free volume is zero in the simulated cell and polymers are free to move within this cell volume. As a result, high values for the entropy and low values for the bulk, shear, and Young's moduli. Further, the incorporation of the graphene into the simulated cells increases free volume and decrease the entropy and increase the values of the bulk, shear, and Young's moduli. Since the entropy of the system is related to the free volume and the Connolly surface of the TPU nanocomposites, the prediction of these parameters may give a concept on the enhancement of the nanocomposite thermal properties. The simulation results show that the TPU matrix reinforced with graphene have a tendency to increase glass transition temperature (Tg) for stronger interlocking between the graphene and TPU molecules. The Connolly surface to volume ratio is small values for the neat TPU and increase with the increase in the graphene loading. The free volume is defined as the volume on the side of the Connolly surface without atoms. This simulation study reveals that Tg of TPU is in the range of 285–305°C and for the Gr/TPU nanocomposite is 350 K (experimental 223–323 K). From the V-T curve, volumetric coefficient of thermal expansion (α) is evaluated from the slope above and below Tg, which is between 2.6 � <sup>10</sup>�<sup>5</sup> and 2.4 � <sup>10</sup>�<sup>4</sup> /K�<sup>1</sup> (experimental 3.15 � <sup>10</sup>�<sup>4</sup> /K�<sup>1</sup> ). The free volume change rapidly when the material goes through glass transition, according to Fox and Flory's theory of glass transition [20]. Further simulation study reveals that graphene/TPU nanocomposite thermal conductivity is 1.5 W/mK, whereas TPU thermal conductivity is 0.2 W/mK. NEMD simulations are used to calculate the thermal conductivity either by imposing a thermal gradient into the system of particles or by introducing a heat flux flow in the reverse nonequilibrium MD (RNEMD) method. In the present study, heat flux flow method is used to calculate thermal conductivity in the longitudinal direction. The total number of layers is 40. Two types of exchange method are used in the present study, namely VARIABLE and FIXED. About 1 kcal/mol energy exchange is taken in FIXED method. The number of exchanges is taken as 500 for equilibrium stage under NVT and 1000 for production stage under NVE. The time steps in between two exchanges are fixed at 100. Due to the presence of the graphene-TPU interface, there exists a temperature jump ΔT at the interface. The present study obtained values are in good agreement with previous values obtained from simulations and experimental measurements in the literatures [26–31]. The present study contributes some novel procedures during MD simulation work which will not be done in previous researchers.

Cp specific heat of nanoparticle bulk material

*DOI: http://dx.doi.org/10.5772/intechopen.86527*

*Molecular Dynamics Simulation-Based Study on Enhancing Thermal Properties…*

α the volumetric coefficient of thermal expansion

k thermal conductivity kB Boltzmann constant Kp the phonon contribution Ke the electronic contribution

R interfacial thermal resistance

ρ density of nanomaterial ri the position of the ith atom Ei total energy of the ith atom

Jq the heat flux

**Author details**

**99**

Animesh Talapatra<sup>1</sup>

\* and Debasis Datta<sup>2</sup>

\*Address all correspondence to: animesh\_talapatra@yahoo.co.in

© 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium,

1 MCKVIE, Howrah, West Bengal, India

provided the original work is properly cited.

2 IIEST, Howrah, West Bengal, India
