Molecular Dynamics Simulation of Heat Transfer ... - ACS Publications

Article pubs.acs.org/JPCC

Molecular Dynamics Simulation of Heat Transfer from a Gold Nanoparticle to a Water Pool Xiaoling Chen,† Antonio Munjiza,† Kai Zhang,§ and Dongsheng Wen*,†,‡ †

School of Engineering and Materials Science, Queen Mary University of London, London E1 4NS, U.K. School of Process, Environmental and Materials Engineering, University of Leeds, Leeds LS2 9JT, U.K. § Beijing Key Laboratory of Emission Surveillance & Control for Thermal Power Generation, North China Electric Power University, Beijing, China ‡

ABSTRACT: Heat transfer from nanoparticles to a surrounding liquid pool is an essential process for many applications. This work investigates the heat transfer process upon a continuous wave laser heating to a gold nanoparticle (GNP) in a water pool through molecular dynamics (MD) simulations based on realistic potentials. The interactions among gold atoms are described by the embedded atom method (EAM) potentials, both the rigid and flexible TIP3P water models are examined, and a modified Lennard-Jones potential is used for the interaction between the GNP and water. The results show that the interfacial thermal conductance is influenced by the selection of different water models and the interfacial wettability. The interfacial thermal conductance is smaller but increases slightly with the increase of heat flux for the flexible water model, while it keeps almost constant for the rigid model. Increasing the wettability between the particle and the fluid reduces the interfacial resistance. The rise of the interfacial water temperature is limited, due to the constraint of the boundary conditions, and no phase change in the water near the GNP is observed in the current simulation. for the formation of nanometer-sized vapor bubble was ∼85% of the critical temperature of water. Recently, molecular dynamics (MD) simulation has become a powerful tool to investigate the nanoscale heat transfer phenomena. MD simulations of gold nanoparticle immersed in liquid have been performed by Merabia et al.,20,21 who showed that a spherical nanoparticle can be heated well above its melting temperature without phase change occurring in the adjacent liquid. The simulation was conducted by applying the simple Lennard-Jones (L-J) interaction potential on all atoms but with additional FENE springs constraint for atoms inside the solid particle. It showed that the interfacial heat flux of a spherical nanoparticle increased linearly and then became saturated with further increasing particle temperature, which was different to a planar interface where the heat flux dropped abruptly due to the formation of a vapor layer. Similar observation was obtained in their studies of GNPs heating in octane and water fluid using a modified L-J potential with the consideration of liquid−solid interactions.21 The absence of the vapor layer was ascribed to be the curvature effect, which requires extreme high pressure close to the spherical nanoparticle. However, for metallic particles, heat is transferred mainly by electrons and phonons, and a more realistic potential, other than the L-J type, should be considered to simulate gold

1. INTRODUCTION Heat transfer from nanoparticles to a surrounding liquid pool is a fundamental process in many medical applications, i.e., controlled drug release, enhanced imaging and diagnosis, noninvasive hyperthermia treatment, and photothermal therapies, to name but a few,1−4 where an optimized heating effect is critical. Among all the particles, gold nanoparticles (GNPs) are increasingly used due to their excellent thermal, optical, chemical, and biological properties.5 Both experimental and theoretical methods have been used to study nanoparticlerelated heat transfer phenomena.6−9 For laser heating GNPs suspended in an aqueous medium, the optical energy is rapidly converted to thermal energy to raise the particle temperature through the surface plasmon resonance (SPR), and then the heat is transferred to the surrounding medium. During this possess, the GNPs and the surrounding medium experience a number of possible different thermo-physical phenomena, such as thermal expansion of particle, particle melting,10−12 and vapor bubble formation,9,13−19 depending on the interaction of the laser and GNPs. These phenomena such as the phase change (i.e., forming vapor bubble) may affect the application favorably or adversely. Recent experiments on laser heating nanoparticles in aqueous media9,13,14 demonstrated the vapor bubble formation phenomenon. For instance, a laser heating experiment based on the pump/X-ray probe method undertaken by Kotaidis et al.9 showed that water around the particles evaporated explosively at high laser fluence, and the threshold © 2013 American Chemical Society

Received: October 9, 2013 Revised: December 12, 2013 Published: December 16, 2013 1285

dx.doi.org/10.1021/jp410054j | J. Phys. Chem. C 2014, 118, 1285−1293

The Journal of Physical Chemistry C

Article

Figure 1. Simulation system. (a) Schematics of whole simulation system (about 80 Å in each dimension); (b) TIP3P water model with bond length b = 0.9572 and bond angle θ = 104.52, and the charge of O and H atoms are −0.830 and +0.415 separately, where σ is the molecular diameter.

where rij is the separation distance between atoms, Øαβ is a pairwise potential function, α and β are the element types of atoms i and j, ρβ is the contribution to the electron charge density from atom j of type β at the location of atom i, and Fα is an embedding function that represents the energy required to place atom i of type α into the electron cloud. Water model used in this work is a 3 site model TIP3P,26 Figure 1b, because of its wide applicability.27 The bonded and nonbonded potentials used to describe the interaction between water atoms are shown below:

atom interactions. In addition, there are many interaction potentials for water molecules, which should be examined more carefully. Consequently, it should be cautious to interpret the simulation results by Merabia et al.,20 especially regarding the comparison to the experimentally observed bubble phenomena around laser excited metal nanoparticles.9 This work will simulate for the first time the heat transfer process from one gold nanoparticle to a pool of water based on realistic potentials. Here, a many-body potential, embedded atom method (EAM), which captures better the physical reality of metallic bonding and is more preferable for FCC metallic structures such as gold, is used. The interaction among water molecules will be constructed by the TIP3P potentials, and a modified L-J will be used for the interaction between the water molecules and gold atoms. Nonequilibrium molecular dynamics simulation (NEMD) will be employed to investigate numerically the heat transfer phenomenon from a gold nanoparticle upon a continuous laser heating to a surrounding pool of water under ambient conditions (i.e., 300 K and 1 atm). Both steadystate density and temperature distribution profiles will be obtained, and the interfacial resistance will be calculated. The influence of the water model (rigid or flexible) and the surface wettability will be investigated.

i≠j

1 2

Ubend = K θ(θ − θ0)2

(3)

Uelectrostatic =

⎛⎛ ⎞12 ⎛ ⎞6 ⎞ σ σ = 4ε⎜⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎟ ⎜⎝ r ⎠ ⎝ rij ⎠ ⎟⎠ ⎝ ij

(4)

qiqj 4πε0rij

(5)

where Ustretch and Ubend are bonded potential terms that correspond to two types of atom movement, bond stretching and angle bending; UVan−der−Waals and Uelectrostatic are two nonbonded terms, the van der Waals interaction energy and the electrostatic interaction energy, respectively; UVan−der−Waals is described by a common expression, i.e., the L-J 12−6 function. Both the L-J and Columbic interactions are intermolecular interactions, which are calculated between nonbonded atoms. The interactions between H2O molecules and Au atoms are calculated only by considering the interactions between oxygen atoms and Au atoms using the L-J 12−6 potential. In order to reduce the computational amount, a cutoff radius of 10.0 Å is used for the force calculation. Table 1 lists the parameters used for TIP3P with a long-range Columbic solver28 and interaction functions. Both the rigid and flexible water models are used to describe water molecular interactions in this study. The flexible model allows the bond stretching and angle bending, while the rigid water model constrains the bond length and angle to specified values using the SHAKE algorithm. All simulations were carried out under a constant pressure of 1 atm using the Nose−Hoover barostat. First, the whole system was conducted under NPT ensemble (300 K, 1 atm) for a

∑ Øαβ (rij) i≠j

(2)

UVan − der − Waals

2. METHODOLOGY Similar to our previous studies,22−24 the simulations were conducted by the large-scale atomic/molecular massively parallel simulator (LAMMPS). Figure 1a presents the diagram of the whole simulation system, which consists of a GNP with diameter of 3 nm containing 820 atoms and 16 736 water molecules. The GNP was cut from the FCC lattice with a density of 19.32 g/cm3 and placed in the center of the cubic box with the edge length of about 80 Å. Periodic boundary condition was imposed on each dimension. For the interactions among gold atoms, different from the L-J combining FENE springs potential used in the simulations of Merabia et al.,20,21 the embedded atom method (EAM) potentials were used in this study, which represents the total energy of the system as two additive terms, a pairwise sum of interactions between atoms, and a term representing the electron density of each atomic site25 Ui(rij) = Fα(∑ ρβ (rij)) +

Ustretch = Kb(b − b0)2

(1) 1286



Article

Table 1. Parameters for TIP3P Water Molecular Model and Interaction Function

a

parameters

values

ε00 (kJ/mol) σ00 (Å) εHH (kJ/mol) σHH (Å) εOH (kJ/mol) σOH (Å) Kb (kJ/(mol Å2)) b0 (Å) Kθ(kJ/(mol rad2)) θ0 (deg) εAuO (kJ/mol)a σAuO (Å)a

0.4270 3.188 0.0 0.0 0.0 0.0 1883.8983 0.9572 230.2542 104.52 0.59 3.6

Figure 2. Transient evolution of the average temperature of a gold nanoparticle (1300 nW). After a short time period, the temperature fluctuates around a value, which indicates that the system reaches steady-state.

Reference 21.

period of 100 000 time-steps to reach the thermal equilibrium state. During this period, while a global Nose−Hover barostat was applied to control the pressure of the whole system, Nose− Hoover thermostat was used for gold and water separately. After the equilibration stage, the global thermostat was removed, and a heat source was applied. The particle was then heated by rescaling the velocities of gold atoms at every time step with a constant heat intensity, while the water at distances >22 Å away from the nanoparticle surface was thermostat at 300 K, acting as a heat sink. The GNP was not allowed to freely diffuse by neglecting the effect of the Brownian motion of particle. After a transient period, a steady state of the whole system was achieved, and the simulation data was collected over a period of 1 ns for the time averaging of required properties.

To obtain the temperature profile, the simulation domain was divided into a number of spherical shells along the radial direction with the origin residing at the center of the GNP. The radial distribution of the concerned properties, temperature and density, was calculated and averaged every 1000 time steps. The temperature definition and calculation are significant for NEMD simulation of the heat transfer problem, especially for the rigid water mode. The temperature was calculated separately for the solid and liquid in spherical shells with 1 Å interval. The thickness of the shell used for temperature calculation was 3 Å, which ensures that the fluctuation of the local temperature is not too large while still keeping the statistical accuracy. Examples of the transient temperature profile for both solid and liquid during the initial heating period are shown in Figure 3. It clearly shows that as the heating time

3. RESULTS ANALYSIS It is well-known that the melting point of a nanoparticle is size dependent. While the bulk melting temperature of gold is 1337 K, early melting would occur at the low dimensions due to its increased specific surface area. The melting depression phenomenon is an established phenomenon for nanoparticles, and a number of models have been proposed to predict the melting temperature at different particle dimensions. To avoid the melting of GNPs due to large heat flux and also as a proof of the reliability of the MD simulation, a separated MD study of the melting of individual gold nanoparticles was conducted. For the 3 nm GNP containing 820 atoms used in this study, the melting point was found to be about 870 K. The comparison of various melting-depression models and experiments shows that the global melting temperature of a 3 nm GNP particle would occur in the range of 800−900 K,29−31 which is consistent with our MD simulations. In order to avoid the melting of GNP, a laser power of 2100 nW was set as the up limit, and the gold particle was heated under a heat intensity varied from 700 to 2100 nW. During the heating, the average temperature of GNP increased rapidly to a steady state (Figure 2), typically in less than 100 ps. After that, the gold particle temperature oscillated around a constant value, indicating that the heat added to gold particle was equal to that leaving the system through the sink region. The steady-state average GNP temperature was shown to be proportional to the heat intensity, and the highest average particle temperature was about 800 K in this work, which is below the melting point.

Figure 3. Illustration of transient temperature distributions under a heat power of 1300 nW (rigid water model).

increases, the nonuniformity of temperature inside the particle increases, and the temperature gradient inside the fluid increases. The steady-state temperature profile, which was obtained by averaging the temperature of each spherical shell over a period of 1 ns after the steady state was reached, is used to calculate the interfacial thermal conductance. 3.1. Rigid Water Model. The steady-state temperature profile, based on the rigid water model, upon different heat fluxes are shown in Figure 4. It is noted that the temperature is not uniform inside the solid phase under all the heat fluxes, which shall be caused by the reduced thermal conductivity of the particle at the nanoscale. The temperature gradient inside 1287



Article

to the interfacial thermal conductance) is comparable to the characteristic length. Understanding the thermal resistance at the interface is very important for the applications containing high density of interfaces. The interfacial thermal conductance is defined as shown below: Q=

ΔT = G K ∇T RK

(6)

lK =

λ GK

(7)

where Q is the heat flux, ΔT is the temperature drop, R is the thermal boundary resistance, and GK is the interfacial thermal conductance, which is the inverse of RK . At the steady state, the interfacial heat flux flowing through the interface can be calculated by dividing the input heat intensity to the particle surface area. The interfacial thermal resistance is highly dependent on the accurate determination of the interface temperature jump. In the literature, the temperature jump is usually obtained by either extrapolating the temperature profile of the bulk region to the interface or by direct calculation over a finite thickness of the interface based on the simulation results. In this study, the temperature jump is calculated through the first approach by fitting the bulk fluid temperature with a theoretical equation for continuous laser heating of gold nanoparticle,8 as shown below:

Figure 4. Steady-state temperature profiles under different heat intensities for simulations (rigid water model). Both MD simulation results and analytical results for laser heating of gold (eqs 8 and 9) are presented for each case. The analytical results are plotted using colored lines.

gold particle would be almost uniform compared with that of water if assuming the bulk gold thermal conductivity. However, similar to the melting temperature described earlier, the thermal conductivity of gold particle is also strongly sizedependent.32−34 It has been well-known that the thermal conductivity of a metal is the sum of the electron and lattice thermal conductivity, and the electron contribution is dominant at the bulk scale. As the particle size becomes comparable or smaller than the intrinsic electron−phonon coupling length scale over which phonons and electrons can equilibrate near an interface, the phonon contribution becomes important. However, the EAM model can only capture the phonon contribution of the thermal conductivity, i.e., the simulated phonon thermal conductivity by EAM in this study is around 1.7 W/(m K), which is consistent with many literature values35,36 and should be the lower bound value of the thermal conductivity. We did not perform detailed electron−phonon coupling simulation of the thermal conductivity of gold, which itself alone is still an unsolved problem, but instead use some established phenomenological models to estimate the rough value. For instance when L≪ λe,b, where λe,b is the mean free path of electrons of bulk material and L is the characteristic size (the diameter of the particle size), the thermal conductivity of a metallic particle can be simply expressed as37 kp/kb = λe,p/λe,b = 1/(1 + Kn), where Kn is the Knudsen number, Kn = λe,b/L, λe,b and kb are the mean free path of electrons and thermal conductivity of the bulk material, respectively, and λe,p is the mean free path of electrons of particle. At 298.15 K, the λe,b value of gold is 36.14 nm, and kb is 315 W/(m K). The estimated electron thermal conductivity of a 3 nm GNP is about 24.14 W/(m K) at 298.15 K, which is much smaller than that of the bulk value. Clearly for 3 nm GNP, the much reduced thermal conductivity value should be responsible for the relatively large temperature gradient inside the gold particle. It is known that when heat flows through an interface, there is a resistance, so-called Kapitza resistance,38 resulting in a temperature jump at the interface due to different thermal properties in the two phases. The interfacial thermal resistance is caused by the property mismatch between the materials. Acoustic mismatch model (AMM) and/or diffusion mismatch model (DMM) are the two widely used methods to predict the interfacial thermal conductance. At the nanoscale, the interfacial resistance may dominate the heat transfer process since the Kapitza length (the ratio of thermal conductivity of one phase

T (r ) =

λ ⎛ λ ⎞ P ⎛ r2 ⎞ ⎜1 + w ⎜1 − 2 ⎟ + w ⎟ + B λAu ⎝ 4πλ w r ⎝ R ⎠ GK R ⎠ (8)

for r < R T (r ) =

P +B 4πλwr

for r > R

(9)

where P is the heat intensity, λw and λAu are the thermal conductivities of water and gold, respectively, R is the particle radius, and B is a constant parameter depending on the temperature boundary used. When the temperature boundary is infinite, the B value set for the boundary condition herein is 300 K. As shown in Figure 4, the water temperature does not change greatly, which approaches to the heat sink temperature in ∼1 nm distance from the interface. So, a constant thermal conductivity of water is a good assumption to fit the analytical equation to the MD results using eq 9. By assuming an average thermal conductivity of water of 0.78 W/(m K), the analytical temperature profile in water agrees with MD results quite well, as shown in Figure 4. At the interface, the gold temperature near the interface deviates slightly from the analytical value, which might be due to the thermal expansion at high heating temperature. The interfacial temperature jump ΔT is then determined by the difference between the solid temperature at the interface and the water temperature extrapolated from the analytical equation. The interfacial thermal conductance is calculated according to eq 6. As the heat intensity increases, the average GNP temperature at the steady state is proportional to the heat intensity. The interfacial thermal conductance is plotted as a function of the heat intensity in Figure 5. It shows clearly that the interfacial thermal conductance is almost a constant value around 212.05 MW/(m2 K), which corresponds to a Kapitza length of 3.68 nm in water according to λw = 0.78 W/(m K). This interfacial resistance value is within the broad 1288



Article

Figure 5. Interfacial thermal conductance GK (MW/(m2 K)) between a 3 nm gold and water (the rigid water model); the cyan filled triangle plot presents the results of GK for a larger simulation system, which will be stated in the discussion section.

Figure 7. Steady-state temperature profiles under different heat intensities by using the TIP3P flexible water model. Both MD simulation results and analytical results for laser heating of gold (eqs 8 and 9) are presented for each case. The analytical results are plotted using colored lines.

range of the previous experimental results, ranging from 100 to 300 MW/(m2 K).39,40 The details about the temperature dependence will be discussed in the sections below. In addition, the water density profile around gold nanoparticle is presented in Figure 6. There are two peaks and

Figure 8. Interfacial thermal conductance GK (MW/(m2 K)) between a 3 nm gold and water (flexible water model).

Figure 6. Number density profiles of water (based on the rigid model) under different heat intensities.

valleys observed for both hydrogen and oxygen atoms, implying the formation of two shell-like water molecules near the gold nanoparticle, due to the adsorption effect (i.e., van der Waals interaction). With the increase of the distance, the interaction between gold atoms and water molecules becomes weak. At a distance greater than 10 Å from the GNP surface, the water molecules exhibit the bulk property. With increasing heat intensity, the heights of the two peaks of hydrogen and oxygen atom number density profile decrease slightly. However, for the current cases, the local densities of oxygen and hydrogen are still larger than that of bulk, which indicates no vapor formed or phase change happened. Such a result is consistent with the temperature plot showing in Figure 4, which clearly shows that the maximum water temperature near the GNP surface is below 400 K. 3.2. Flexible Water Model. To study the effect of the water potential model on the interfacial thermal conductance, MD simulations using the flexible TIP3P water model, which allows the bond stretching and bond bending, were conducted under several heat intensities. The results are presented in Figures 7−9. Figure 7 shows the steady-state temperature distribution. Here, the water thermal conductivity used to fit the analytical temperature equation, eq 9, to the MD results is assumed to be 1.06 W/(m K), which is larger than that of the

Figure 9. Number density profiles of water (flexible water model) under different heat intensities.

rigid water model. Such a larger thermal conductivity is due to the nature of the water model. Similar values have been reported in other publications based on the Müller-Plathe (MP) method, a nonequilibrium molecular dynamics method where a known heat flux is imposed to create temperature gradient and the thermal conductivity is obtained from the temperature profile.41 Using the same TIP3P water model as in the current study, the determined thermal conductivity from the M-P method for rigid and flexible TIP3P water model was 0.814 and 1.063 W/(m K), respectively.41 These values are quite close to the thermal conductivity values used in this study. Different from the results of the rigid water model, there is a general trend that the interfacial thermal conductance increases 1289



Article

with the heat flux in the range of 160−200 MW/(m2 K) (Figure 8). A higher heat flux results in a higher interface temperature, consequently the interfacial thermal conductance increases with the interface temperature, which is consistent with a few other studies.42,43 The interfacial thermal conductance is slightly smaller than that obtained from the rigid water cases but is larger than those in ref 21 (150−170 MW/(m2 K)). This is due to the increased property mismatch between the solid/liquid for the flexible water model. The density of states (DOS) of the GNP and bulk water were calculated, Figure 10, to provide an additional explanation for

effect of wettability (i.e., the interactive strength between the particle and the medium), stronger Au−H2O interaction parameters were chosen from a previous MD study:44 εAuO = 8.3 kJ/mol and σAuO = 2.95 Å, which represents a wetting case comparable to the experimental adsorption energy data of water on gold. The originally used parameters are εAuO = 2.47 kJ/mol and σAuO = 3.6 Å,21 corresponding to a poor wetting case. The steady-state temperature distribution using the two sets of parameter is shown in Figure 12. It clearly shows that the NP

Figure 10. Comparison of the density of states (DoS) of GNP using the rigid and flexible water models.

Figure 12. Comparison of temperature profiles obtained from simulation with different Au−H2O interaction parameters under a heat intensity of 1300 nW. The colored lines are the analytical temperature distributions, which are fitted to simulation results.

the different interfacial thermal conductance values. The DOS of the flexible water model has two peaks at high frequency corresponding to the bond stretching and torsion, but at the low frequency, the magnitude of DOS for the rigid water is larger than that of the flexible water model. Since GNPs are predominantly in the low frequency region, the coupling between the GNP and water is stronger for the rigid water model, which results in a larger interfacial thermal conductance. This explains clearly the lower steady-state solid temperature by using the rigid water model, Figure 11.

temperature can be significantly affected by different solid/ liquid interaction parameters, but the water temperature profiles are similar quantitatively for both cases. As a result, the interfacial thermal conductance for the wetting case is much higher than that of the weaker interaction strength. The interfacial thermal conductance of the wetting interface is determined as 621.48 and 396.73 MW/(m2 K) for the case of the rigid water and flexible water models, respectively. Such a result is consistent with the theoretical description that the steady-state solid temperature is dependent on the interfacial thermal conductance, but the fluid temperature is mainly dependent on the heat intensity rather than the interfacial thermal conductance. For a given heat flux, the stronger interaction provides a fast cooling and results in a smaller increase of the steady-state NP temperature and hence a higher interfacial thermal conductance. The number density shows that as the GNP surface becomes wetter, the interfacial water layering becomes denser, and the peak is closer to GNP surface with an increased magnitude, Figure 13, which also implies a faster heat transfer from the nanoparticle to its surrounding medium.

4. DISCUSSION 4.1. Interfacial Thermal Conductance. Among the literature, the influence of temperature or heat flux on the interfacial thermal conductance is not consistent. While some studies reported an increase of the interfacial thermal conductance with the increase of particle temperature,42,43,45 others reported a reversed trend46,47 and a nondependent relationship.48 The results are likely influenced by the variation of the liquid layering structure near the interface upon temperature. For instance, Goicocha et al.48 reported constant Kapitza conductance values for water and self-assembly monolayer interfaces based on the rigid water model. The

Figure 11. Comparison of steady-state temperature profiles at 700 and 900 nW for the rigid and flexible water models; the colored lines are the analytical temperature distributions, which are fitted to simulation results.

3.3. Effect of the Wettability. Not only the water potentials but also the interaction potential between water and GNP can influence the wetting property and affect largely the interfacial conductance value. The Kapitza conductance/ resistance at the interface can be modified by the property mismatch (such as mass) and the interaction strength of the two materials forming the interface. In order to explore the 1290



Article

of nanobubbles or microbubbles around small particles9,51,52. For instance, by applying X-ray scattering to the laser excited particle solution, Plech and Kotaidis et al.9,53 showed that gold particles of 9 nm diameter could be excited by a nanosecond laser to produce large vapor bubbles that decayed within less than 1 nanosecond. Upon the generation of vapor bubbles, it insulated the particle and surrounding media that increased particle temperature. Kotaidis et al.9 showed that the water around the particles evaporated explosively at high laser fluence, and the threshold for the formation of nanometer-sized vapor bubble lied in ∼85% of the critical temperature of water. The existence of nanobubbles has also been analytically predicted. For instance, Wen54 showed that the formation of nanobubbles has an explosive nature, which can reach a few micrometers in around 10 ns, depending on the particle size and laser pulse selection. The rapid formation and contraction of bubbles around heated nanoparticles, associated with the propagation of pressure waves, were effective at a dimension much larger than that of a nanoparticle. In this study, no bubble formation was detected on top of a strongly heated 3 nm gold nanoparticle in water. The interfacial water temperature is far away from the critical point, and the current result is similar to that of Merabia et al.21 As a finite thermal boundary condition is typically used in MD simulations to reduce the computational cost, it would limit the rise of the interfacial temperature layer. To examine the boundary condition effect, MD simulation was conducted under two heat intensities in a large dimension. Herein, the thermal boundary is set 32 Å away from the GNP surface, rather than 22 Å as in the above simulations. The comparison shows that the interfacial thermal conductance value is similar to that obtained for the smaller system, Figure 5. Comparing the temperature distributions, however, the interfacial temperature becomes larger for the larger thermal boundary condition case, Figure 14. This is due to the extension of the heat sink position, removing the previous constraint of T = 300 K at 22 Å. The temperature distribution obtained by resolving the heat equation with infinite temperature boundary condition is also plotted in Figure 14, which shows that very high interface temperature could be reached. It suggests that some improvement in the simulation technique should be used in the future

Figure 13. Number density distribution (1300 nW) for the two different Au−H2O interactions.

temperature dependence of the interfacial thermal conductance is also affected by the surface wettability.45,49,50 At a hydrophobic interface, the density profile is more sensitive to the temperature. Consequently, the effect of the temperature on the interfacial thermal conductance for a less hydrophilic interface should be stronger than that for hydrophilic interface. It is found that the Kapitza length increases with the wall temperature47,48 when a solid wall is kept at the prescribed temperature. It is believed that such inconsistent results in the literature are due to the implementation of the simulations and the choice of the interaction potentials. From the simulation consideration, different results would be obtained if (i) inducing heat flow by adding and abstracting constant energy in given regions or (ii) by keeping the temperature at a given value using the thermostat methods. From the potential consideration, this study shows clearly that the choice of water model and the interaction potential between water and GNP can influence the result significantly. The interfacial thermal conductance becomes independent of the temperature, if the rigid water model was chosen (Figure 5), or increases with temperature if the flexible water was used, Figure 8. The temperature dependence shall be the results of the balance between the effect of phonon distribution variation and the effect of the variation in the interfacial water layer due to the change of the temperature. Shown when comparing the number density distributions of the two water models (Figures 6 and 9), it is clear that the density varies stronger at elevated surface temperatures for the flexible model. Not only the trend but also the choice of different fluid models affects significantly the value of the interfacial thermal conductance. For instance, the flexible water model produced overall smaller values, which is due to the higher overlap of density of states for the rigid water model. It should be noted that for the flexible model, many of the vibration modes of the water molecule are not thermally excited under low temperature, and their contribution to heat transfer is not evidenced in this work. Conclusively, the selection of different fluid models can influence the measurement of interfacial thermal conductance, the magnitude, and the temperature dependence. So it needs to be careful in analyzing the results regarding the interfacial thermal transport. 4.2. Nanobubble Phenomenon. As briefed in the introduction, there are still some debates about whether vapor bubbles could be formed on top of heated nanoparticles using MD simulation. There have been recently some experimental studies demonstrated that the physical existence

Figure 14. Temperature profiles at 700 and 1300 nW. For the larger system, heat sink is set at 32 Å away from the GNP surface; and for the infinite boundary condition, the heat sink is very far away. The temperature profile is computed through eqs 8 and 9 with B = 300 K. 1291



Article

(9) Kotaidis, V.; Dahmen, C.; Plessen, G.; Springer, F.; Plech, A. Excitation of Nanoscale Vapor Bubbles at the Surface of Gold Nanoparticles in Water. J. Chem. Phys. 2006, 124, 184702. (10) Link, S.; Burda, C.; Nikoobakht, B.; El-Sayed, M. A. LaserInduced Shape Changes of Colloidal Gold Nanorods Using Femtosecond and Nanosecond Laser Pulses. J. Phys. Chem. B 2000, 104, 6152−6163. (11) Lukianova-Hleb, E. Y.; Anderson, L. J. E.; Lee, S.; Hafner, J. H.; Lapotko, D. O. Hot Plasmonic Interactions: A New Look at the Photothermal Efficacy of Gold Nanoparticles. Phys. Chem. Chem. Phys. 2010, 12, 12237−12244. (12) Plech, A.; Kotaidis, V.; Grésillon, S.; Dahmen, C.; von Plessen, G. Laser-Induced Heating and Melting of Gold Nanoparticles Studied by Time-Resolved X-ray Scattering. Phys. Rev. B 2004, 70, 195423. (13) Dou, Y.; Zhigilei, L.; Winograd, N.; Garrison, B. J. Explosive Boiling of Water Films Adjacent to Heated Surfaces: A Microscopic Description. J. Phys. Chem. A 2001, 105, 2748−2755. (14) Garrison, B. J.; Itina, T. E.; Zhigilei, L. V. Limit of Overheating and the Threshold Behavior in Laser Ablation. Phys. Rev. E 2003, 68, 041501. (15) Neumann, J.; Brinkmann, R. Boiling Nucleation on Melanosomes and Microbeads Transiently Heated by Nanosecond and Microsecond Laser Pulses. J. Biomed. Opt. 2005, 10, 24001. (16) Vogel, A.; Linz, N.; Freidank, S.; Paltauf, G. Femtosecond-LaserInduced Nanocavitation in Water: Implications for Optical Breakdown Threshold and Cell Surgery. Phys. Rev. Lett. 2008, 100, 038102. (17) Lapotko, D. Optical Excitation and Detection of Vapor Bubbles Around Plasmonic Nanoparticles. Opt. Express 2009, 17, 2538−2556. (18) Taylor, R. A.; Phelan, P. E.; Otanicar, T.; Adrian, R. J.; Prasher, P. S. Vapor Generation in a Nanoparticle Liquid Suspension Using a Focused, Continuous Laser. Appl. Phys. Lett. 2009, 95, 161907. (19) Lukianova-Hleb, E.; Hu, Y.; Latterini, L.; Tarpani, L.; Lee, S.; Drezek, R. A.; Hafner, J. H.; Lapotko, D. O. Plasmonic Nanobubbles as Transient Vapor Nanobubbles Generated Around Plasmonic Nanoparticles. ACS Nano 2010, 4, 2109−2123. (20) Merabia, S.; Keblinski, P.; Joly, L.; Lewis, L. J.; Barrat, J. L. Critical Heat Flux Around Strongly Heated Nanoparticles. Phys. Rev. E 2009, 79, 021404. (21) Merabia, S.; Shenogin, S.; Joly, L.; Keblinski, P.; Barrat, J. L. Heat Transfer From Nanoparticles: A Corresponding State Analysis. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 15113−15118. (22) Song, P.; Wen, D. S. A Reactive Molecular Dynamic Simulation of Oxidation of a Silicon Nanocluster. J. Nanopart. Res. 2013, 15 (1), 1309. (23) Song, P.; Wen, D. S. Molecular Dynamic Simulation of Heating and Cooling of a Core-Shell Structured Metallic Nanoparticle. J. Phys. Chem. C 2010, 114 (19), 8688−8696. (24) Song, P.; Wen, D. S. Molecular Dynamics Simulation of the Sintering of Metallic Nanoparticles. J. Nanopart. Res. 2010, 12 (3), 823−829. (25) Daw, M. S.; Baskes, M. I. Embedded-Atom Method: Derivation and Application to Impurities, Surfaces, and Other Defects in Metals. Phys. Rev. B 1984, 29, 6443−6453. (26) Mahoney, M. W.; Jorgensen, W. L. A Five-Site Model for Liquid Water and the Reproduction of the Density Anomaly by Rigid, Nonpolarizable Potential Functions. J. Chem. Phys. 2000, 112, 8910− 8922. (27) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. 1998, 102, 3586−3616. (28) Price, D. J.; Brooks, C. L., III. A Modified TIP3P Water Potential for Simulation With Ewald Summation. J. Chem. Phys. 2004, 121, 10096. (29) Buffat, P.; Borel, J. P. S. Size Effect on the Melting Temperature of Gold Particles. Phys. Rev. A 1976, 13, 2287−2298. (30) Lewis, L. J.; Jensen, P.; Barrat, J. L. Melting, Freezing, and Coalescence of Gold Nanoclusters. Phys. Rev. B 1997, 56, 2248−2257.

to examine the heating phenomenon, including more realistic boundary conditions, such as the pressure transmission boundary condition, or though the multiscale approach covering a sufficiently large dimension.

5. CONCLUSIONS In this study, the heat transfer phenomenon from a 3 nm gold nanoparticle to a surrounding water pool is simulated by molecular dynamics employing EAM and TIP3P potentials. It shows that the simulation result is significantly affected by the choice of water models, the interfacial wetting, and the boundary conditions. It shall be cautious in comparing interfacial thermal conductance from different sources. The interfacial thermal conductance is found to be smaller but increases with the increase of heat flux for the flexible water model, while it keeps almost constant for the rigid model. The increased wettability between the particle and the fluid increases the interfacial thermal conductance significantly. The extension of the heat sink position to a larger distance increases the interfacial water temperature. However, the simulation does not show the phase change phenomenon as reported by some experimental studies, likely due to the limitation of the boundary conditions.

■

AUTHOR INFORMATION

Corresponding Author

*(D.W.) Tel: 0044 113 3431299. Fax: 0044 113 3435001. Email: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We wish to express our gratitude to the financial support from the China Scholarship Council, the National Science Foundation of China through its grant No. 51228601 and the 111 Project (B12034).

■

REFERENCES

(1) El-Sayed, I. H.; Huang, X.; El-Sayed, M. A. Surface Plasmon Resonance Scattering and Absorption of Anti-EGFR Antibody Conjugated Gold Nanoparticles in Cancer Diagnostics: Applications in Oral Cancer. Nano Lett. 2005, 5, 829−834. (2) Han, G.; Ghosh, P.; De, M.; Rotello, V. M. Drug and Gene Delivery Using Gold Nanoparticles. Nanobiotechnology 2007, 3, 40− 45. (3) Chakravarty, P.; Qian, W.; El-Sayed, M. A.; Prausnitz, M. R. Delivery of Molecules Into Cells Using Carbon Nanoparticles Activated by Femtosecond Laser Pulses. Nat. Nanotechnol. 2010, 5, 607−611. (4) Chen, H.; Wen, D. S. Electromagnetic heating of gold nanoparticle dispersions at 200 kHz. Nanomedicine 2013, 8 (2), 215−222. (5) Qin, Z.; Bischof, J. C. Thermophysical and Biological Responses of Gold Nanoparticle Laser Heating. Chem. Soc. Rev. 2012, 41, 1191− 1217. (6) Govorov, A. O.; Zhang, W.; Skeini, T.; Richardson, H.; Lee, J.; Kotov, N. A. Gold Nanoparticle Ensembles as Heaters and Actuators: Melting and Collective Plasmon Resonances. Nanoscale Res. Lett. 2006, 1, 84−90. (7) Del Fatti, N.; Arbouet, A.; Vallée, F. Femtosecond Optical Investigation of Electron−Lattice Interactions in an Ensemble and a Single Metal Nanoparticle. Appl. Phys. B: Laser Opt. 2006, 84, 175− 181. (8) Baffou, G.; Rigneault, H. Femtosecond-Pulsed Optical Heating of Gold Nanoparticles. Phys. Rev. B 2011, 84, 035415. 1292



Article

(31) Arcidiacono, S.; Bieri, N. R.; Poulikakos, D.; Grigoropoulos, C. P. On the Coalescence of Gold Nanoparticles. Int. J. Multiphase Flow 2004, 30, 979−994. (32) Fang, K. C.; Weng, C. I.; Ju, S. P. An Iinvestigation Into the Structural Features and Thermal Conductivity of Silicon Nanoparticles Using Molecular Dynamics Simulations. Nanotechnology 2006, 17, 3909−3914. (33) Teja, A. S.; Beck, M. P.; Yuan, Y.; Warrier, P. The Limiting Behavior of the Thermal Conductivity of Nanoparticles and Nanofluids. J. Appl. Phys. 2010, 107, 114319. (34) Yuan, S. P.; Jiang, P. X. Thermal Conductivity of Small Nickel Particles. Int. J. Thermophys. 2006, 27, 581−595. (35) Richardson, C. F.; Clancy, P. Contribution of Thermal Conductivity to The Crystal-regrowth Velocity of Embedded-atommethod-modeled Metals and Metal Alloys. Phys. Rev. B 1992, 45, 12260−12268. (36) Kuang, S.; Gezelter, J. D. A Gentler Approach to RNEMD: Nonisotropic Velocity Scaling for Computing Thermal Conductivity and Shear Viscosity. J Chem. Phys. 2010, 133, 164101. (37) Warrier, P.; Teja, A. Effect of Particle Size on the Thermal Conductivity of Nanofluids Containing Metallic Nanoparticles. Nanoscale Res. Lett. 2011, 6, 247:1−6. (38) Barrat, J. L.; Chiaruttini, F. Kapitza Resistance at the LiquidSolid Interface. Mol. Phys. 2003, 101, 1605−1610. (39) Ge, Z.; Gahill, D. G.; Braun, P. V. AuPd Metal Nanoparticles as Probes of Nanoscale Thermal Transport in Aqueous Solution. J. Phys. Chem. B 2004, 108, 18870−18875. (40) Wilson, O. M.; Hu, X.; Cahill, D. G.; Braun, P. V. Colloidal Metal Particles as Probes of Nanoscale Thermal Transport in Fluids. Phys. Rev. B 2002, 66, 224301. (41) Sirk, T. W.; Moore, S.; Brown, E. F. Characteristics of Thermal Conductivity in Classical Water Models. J. Chem. Phys. 2013, 138, 064505. (42) Murad, S.; Puri, I. K. Thermal Transport Across Nanoscale Solid-Fluid Interfaces. Appl. Phys. Lett. 2008, 92, 133105. (43) Murad, S.; Puri, I. K. Molecular Simulation of Thermal Transport Across Hydrophilic Interfaces. Chem. Phys. Lett. 2008, 467, 110−113. (44) Tay, K. A.; Bresme, F. Wetting Properties of Passivated Metal Nanocrystals at Liquid−Vapor Interfaces: A Computer Simulation Study. J. Am. Chem. Soc. 2006, 128, 14166−14175. (45) Ge, S.; Chen, M. Temperature Dependence of Thermal Resistance at a Solid/Liquid Interface. Mol. Phys. 2013, 111, 903−908. (46) Kim, B. H.; Beskok, A.; Cagin, T. Molecular Dynamics Simulations of Thermal Resistance at the Liquid-Solid Interface. J. Chem. Phys. 2008, 129, 174701. (47) Shi, Z.; Barisik, M.; Beskok, A. Molecular Dynamics Modeling of Thermal Resistance at Argon−Graphite and Argon−Silver Interfaces. Int. J. Therm. Sci. 2012, 59, 29−37. (48) Goicochea, J. V.; Hu, M.; Miechel, B.; Poulikakos, D. Surface Functionalization Mechanisms of Enhancing Heat Transfer at SolidLiquid Interfaces. J. Heat Transfer 2011, 133, 082401. (49) Ge, Z.; Cahill, D. G.; Braun, P. V. Thermal Conductance of Hydrophilic and Hydrophobic Interfaces. Phys. Rev. Lett. 2006, 96, 186101. (50) Xue, L.; Keblinski, P.; Phillpot, S. R.; Choi, S. U.-S.; Eastman, J. A. Two Regimes of Thermal Resistance at a Liquid−Solid Interface. J. Chem. Phys. 2003, 118, 337−339. (51) Kudryashov, S. I.; Allen, S. D. Submicrosecond Dynamics of Water Explosive Boiling and Lift-Off From Laser-Heated Silicon Surfaces. J. Appl. Phys. 2006, 100, 104908. (52) Neumann, J.; Brinkmann, R. Self-Limited Growth of LaserInduced Vapor Bubbles Around Single Microabsorbers. Appl. Phys. Lett. 2008, 93, 033901. (53) Plech, A.; Kotaidis, V.; Lorenc, M.; Wulff, M. Thermal Dynamics in Laser Excited Metal Nanoparticles. Chem. Phys. Lett. 2005, 401, 565−569. (54) Wen, D. Intracellular Hyperthermia: Nanobubbles and Their Biomedical Applications. Int. J. Hyperthermia 2009, 25, 533−541. 1293


Molecular Dynamics Simulation of Heat Transfer ... - ACS Publications

Recommend Documents