Molecular Dynamics Simulation of Effects of Stretching and Compressing on Thermal Conductivity of Aligned Silicon Oxygen Chains

Wen-Xue Xu(徐文雪), Xin-Gang Liang(梁新刚)$^{\hyperlink{s}{**}}$

doi:10.1088/0256-307X/37/4/046601

⇦Chinese Physics Letters, 2020, Vol. 37, No. 4, Article code 046601⇨ Molecular Dynamics Simulation of Effects of Stretching and Compressing on Thermal Conductivity of Aligned Silicon Oxygen Chains * Wen-Xue Xu (徐文雪), Xin-Gang Liang (梁新刚)** Affiliations School of Aerospace Engineering, Key Laboratory for Thermal Science and Power Engineering of Ministry of Education, Tsinghua University, Beijing 100084 Received 18 December 2019, online 24 March 2020 ^*Supported by the National Science Fund for Creative Research Groups (Grant No. 51621062).
^**Corresponding author. Email: liangxg@tsinghua.edu.cn Citation Text: Xu W X and Liang X G 2020 Chin. Phys. Lett. 37 046601 Abstract The effects of stretching and compressing on the thermal conductivity (TC) of silicon oxygen chain are studied by means of non-equilibrium molecular dynamics simulation. It is found that stretching can improve TC, and compressing may reduce the TC and can also increase the TC. This mechanism is explained based on the variation of phonon group velocity and the specific heat per volume with stretching and compressing. The distributions of bond angle and bond length under different normalized chain lengths are given. It is found that the bond length and bond angle in the skeleton chain would deviate from their original position. In addition, the phonon density of states (PDOSs) of silicon and oxygen atoms in the chains under different normalized chain lengths are analyzed. The overall trend is that the TC increases and the peaks of PDOSs move towards higher frequency with increasing stretch strain. DOI:10.1088/0256-307X/37/4/046601 PACS:66.10.cd, 66.30.hk, 66.70.-f © 2020 Chinese Physics Society Article Text It is very common and usual to find high heat flux in electronic devices with the rapid development of semiconductor technology. Devices with high heat flux have increasing application scenarios, which also pose severe challenges to heat dissipation. To reduce thermal resistance, the gap between a device and its substrate of cooling is usually filled with thermal interfacial materials with high TC.[1] Products of thermal interfacial materials on the market include thermal conductive paste, thermal conductive elastic cloth, thermal conductive gel, and phase-changing thermal conductive adhesive, etc.[2] Silicone rubber is the most widely used thermal interfacial material because it is of excellent performance such as high temperature resistance, chemical inertia, elasticity, high TC, etc. To improve the TC of silicone rubber, a lot of studies have been conducted,[1] which are mainly concerned with the TC of materials from the aspects of particle additives,[3–6] interfacial interaction[7] and process treatment.[5,8,9] However, there are few studies on the matrix structure influence of the silicone rubber on the TC. The study of polyethylene has shown that amorphous polyethylene has a low TC while the polyethylene after directional alignment could improve 2–3 orders of magnitude.[10–12] Gao and Meng[13] used the non-equilibrium molecular dynamics simulation (NEMD) method to study the stretching effect on the TC of three kinds of nanometer materials, i.e., C, BN and SiC, and found that their TC would increases first and then decreases as the axial stretching strain increases. Lin et al.[14] calculated the TC of polyethylene chains with different stretching strains along the direction of skeleton chain by the NEMD method. They found that the TC increases with the stretching strain, and discussed the thermal transport mechanism of polymer chains by analyzing the vibration density of states. Wu et al.[15] applied molecular dynamics method to investigate the influence of stretching strain on the TC of single layer black phosphorus nanoribbon at room temperature in different directions. The results showed that stretching strain affects the TC of nanoribbon by changing PDOSs and mean free path (MFP). Dong et al.[16] studied the TC of germanium with different stretching strains with the equilibrium molecular dynamics (EMD) method. Their results showed that the TC either in the direction of strain or perpendicular to the strain would decrease with the increase of stretching strain. Zhang et al.[17] simulated the effect of stretching strain on the TC of germanium thin films using the NEMD method, and found that TC of germanium decreases with the increase of stretching strain, while increases with the increase of compressing strain. Tang[18] studied the stretching strain of graphene and found that the TC of graphene decreases with the increase of stretching strain, and the higher the peak frequency value of PDOSs are, the higher the TC is. The available references that investigated the effect of stretching strain on TC are listed in Table 1.[13–22] From Table 1, it seems that the TC for crystals (especially for metal) decreases monotonically as the stretching strain increases,[16–21] with three exceptions[13,15,22] in which the graphene material shows two opposite conclusions.[18,22] This is because the material in Ref. [22] is serpentine. There are few references concerned with amorphous materials and only in Ref. [14] the authors found that the TC increases as the stretching strain increases. To the our best knowledge there has been no study on effects of stretching and compressing strain on the TC of silicone rubber with aligned silicon oxygen chains up to date, which is the investigation purpose of the present study. Therefore, whether stretching and compression would affect the TC of silicone rubber and how they would affect the TC should be understood.

**Table 1.** The variation trend of TC with the increase of stretching strain.
References	Variation tendency with enhancing stretching strain	Materials	Structure
[13]	Increase first & Then decrease	C, BN, SiC	Nonmetal, crystal
[14]	Monotonically increase	PE	Nonmetal, amorphous
[15]	Monotonically increase	Black phosphorus	Nonmetal, crystal
[16]	Monotonically decrease	Germanium	Metal, crystal
[17]	Monotonically decrease	Germanium thin films	Metal, crystal
[18]	Monotonically decrease	Graphene	Nonmetal, crystal
[19]	Monotonically decrease	Graphyne nanotubes	Nonmetal, crystal
[20]	Monotonically decrease	Graphyne nanotubes	Nonmetal, crystal
[21]	Monotonically decrease	Monolayer WSe$_2$	Nonmetal, crystal
[22]	Monotonically increase	Graphene nanoribbons	Nonmetal, crystal

There are two common molecular dynamics methods,[23] namely EMD and NEMD. The NEMD method has the advantage of fast convergence. It is essential to apply temperature bias to make the system in non-equilibrium state. After the system is stable, the heat flux and temperature gradient are counted to calculate the TC by Fourier's law. A common NEMD method[24] is to establish the temperature distribution and calculate heat flux, such as the dual-temperature method.[25] In this Letter, the dual-thermostat NEMD method[26,27] is adopted to calculate the TC of silicone rubber. The boundary condition is set to be thermal bath boundary, and the two regions are named as hot and cold regions, respectively. The energy entering the hot region and the energy leaving the cold region are calculated to obtain an average heat flux. The system is divided into 40–50 segments along the heat flux direction, depending on the chain length. The temperature of each segment is calculated and plotted to obtain temperature gradient. Then Fourier's law is adopted to get the TC: $$\begin{align} &J=-\lambda \nabla {\rm T},~~ \tag {1} \end{align} $$ $$\begin{align} &J=\frac{\Delta E}{A\Delta t},~~ \tag {2} \end{align} $$ $$\begin{align} & \lambda =-\frac{\left\langle {J(t)} \right\rangle }{\left\langle {{\partial T} / {\partial z}} \right\rangle }.~~ \tag {3} \end{align} $$ Equation (1) is the relationship between TC $\lambda$ and heat flux $J$. The heat flux $J$ is the energy per unit time through per unit cross-sectional area that is perpendicular to the direction of the heat flux. It can be calculated by Eq. (2) for a given time duration $\Delta t$, where $\Delta E$ is the total heat flow through area $A$. $\nabla T$ is the temperature gradient along the direction of the heat flux, the angular brackets mean the ensemble average. The main components of silicone rubber and silicone oil are silicon oxygen chains. The skeleton chains of silicon oxygen chains are the silicon atoms and oxygen atoms, and two methyl groups are connected to the silicon atoms. Firstly, a monomer is established by the parameters of bond lengths and bond angles. Then, the monomer is annealed and its energy is minimized in the Materials Studio (MS) software to determine the length of the monomer in chain direction, which is 5.0 Å. Then, the Visualizer Module of the MS software is used to duplicate the monomer and set methyl as the end group. Single chains with lengths of 10, 25, 50 and 80 nm are obtained, and the crystal structure as shown in Fig. 1 is established with a density of 0.963 g/cm$^3$.[28] The focus of the present study is to explore the stretch and compression on the influence of TC, the influence of point defect, environment temperature etc. on the TC[29] is not considered here. The force field is set in COMPASS.[30–35] Then the Minimizer Model is used to optimize its structure by annealing operation until a stable structure is obtained. For the same relaxed structure, relaxation time steps of 0.05 fs, 0.1 fs, 0.25 fs and 0.4 fs are tried to observe energy conservation of the NVE (micro-canonical) ensemble. It is found that discontinuity of energy integration in the NVE ensemble during the simulation happens if time step is too large, which may be caused by fast vibrating hydrogen atoms in the silicon-oxygen chain. Therefore, the time step is set to be 0.05 fs. Next, the calculation of TC of the system is carried out. Steps of molecular dynamics simulation from two million to four million under NVT ensemble are applied at 300 K to fully relax the structure of the simulation system again. Steps of 20 to 30 million are run to calculate the TC in the NVE ensemble. In the cold and hot regions, the Langevin method[36,37] is used to control their temperature, respectively. After the heat flux and temperature distribution are stable, the heat flux and temperature gradient are used to calculate the TC. Both the non-bond cutoff radius and the Coulomb force truncation radius are 9.5 Å in the calculation. The TC of silicone rubber at 300 K calculated by the above process is 0.219 W/(m$\cdot$K), which is consistent with the measured value of 0.25 W/(m$\cdot$K) at 298.15 K,[38] indicating the accuracy of the method.

Fig. 1. Single chain of silicon oxygen and cross-sectional view and side view of the simulation box. The yellow atoms are silicon, the red ones are oxygen, the gray ones are carbon, and the white ones are hydrogen, the direction of the heat flow is shown by the arrow.

The crystal structures with a total length of 10, 25, 50 and 80 nm in the $Z$ direction are stretched and compressed to different degrees. The degree of the stretching or compressing is expressed by the normalized length with respect to the original length $L$. The TC of the silicone rubber with the normalized length of 0.90, 0.95, 1.00, 1.01, 1.05, 1.10 and 1.15 (from compressing to stretching) are calculated, respectively. The cross-sectional area of the system is 21.5 Å $\times 21.5$ Å, and the average temperature is set to be 300 K, which is fixed in this study. A typical temperature distribution along the chain is depicted in Fig. 2. Many researchers have shown that properties depend on temperature.[39,40] However, the present research focuses on the effects of stretching and compressing of the silicone rubber. The temperature dependence is not considered here. The temperature is linear between the thermal baths, indicating a stabilized state. The variation of TC with the normalized length is shown in Fig. 3. As can be seen, the TC increases with the increase of stretching length. On the other hand, the variation of the TC with the compressing is not consistent. For longer chains (50 nm, 80 nm), their TC first decreases and then increases while for shorter chains (10 nm, 25 nm) it increases with the increase of the compressing ratio.

Fig. 2. A typical temperature profile.

Fig. 3. Effects of stretching/compressing ratio on TC.

The commonly used expression $\lambda =cvl/3$[14,15,41] is used to explain the variation of the TC, where $c$ is specific heat of the phonons per volume, $v$ is the group velocity, $l$ is the phonon MFP. In the case of stretching, the tension of the molecular chain is larger. According to the phonon theory of one-dimensional atomic chain, it can be known that the phonon group velocity $v$ increases with the increase of the tension. In addition, when the chain is stretched, the transverse vibration of atoms is depressed and there will be less disturbance on the longitudinal vibration and the MPF for longitudinal phonon will increase. The decrease of $c$ is a negative factor but with less effect on TC in terms of the value of TC. The increase of $v$ and $l$ in stretching case increases the TC. When the structure is compressed, the tension between atoms decreases and hence the group velocity decreases. On the other hand, the number of atoms per unit volume, i.e., the density, increases and the phonon modes per volume increases too, which increases $c$ (specific heat per volume increases positively with phonon modes) in turn. The combination of these two effects makes the TC increase or decrease, or first decrease then increase with the compression. In addition, these two competing effects make compression have less effect on TC than stretching.

Fig. 4. Distribution of bond angles at different normalized lengths.

To explore the change of bond length and bond angle under stretching or compressing and its influence on the TC, the probabilities of bond and angle distributions at different stretching strains is calculated by analyzing trajectory document in MS. An NVE ensemble calculation is performed and a configuration file every certain time step is output. In each configuration file, there is information about atoms' position, so that the corresponding bond angle and bond length information through geometric calculations can be obtained. By analyzing the information of each configuration, the distribution probabilities of the bond angle or bond length are obtained, as shown in Fig. 4. It can be seen that the distributions of bond angles under different stretching strains are slightly different. The angles of O–Si–O and Si–O–Si become generally larger than the original ones when stretching strain is applied. When compressing strain applies to the system, these two angles become smaller than the originals. In short, stretching strain and compressing strain will deviate bond angles. On the other hand, the bond angles off the skeleton chain are less affected by stretching and compressing. The distribution of the Si–C–H angle remains almost unchanged. However, the distribution of the O–Si–C angle is complicated. The section of angle concentration does not change much, with one exception. At the normalized length of 0.90, the O–Si–C angle even has three peaks, two of them deviate from their original positions. Figure 5 is the distribution of bond length with the normalized chain length. As can be seen, stretching and compressing have little effect on bond length of C–H that is not on the skeleton chain, the distributions of bond lengths of Si–O and Si–C are affected to some degree. As can be seen from Fig. 5, at the normalized length of 0.90, bond Si–O distribution is almost unchanged, but the bond length of Si–C becomes smaller. The bond length distribution of Si–C is almost unchanged at the normalized length of 0.95, but that of Si–O becomes larger. In summary, the bond length distribution of Si–O or bond Si–C has certain deviation in a same stretching ratio or compressing ratio.

Fig. 5. Bond length distribution under different stretching ratios.

The variation of bond length could help to explain the variation of TC with normalized chain length. From the theory of solid physics, the group velocity $v=a(g/m)^{0.5}$ for one-dimensional atomic chain where $a$ is the lattice constant, $g$ is the force constant between adjacent atoms and $m$ is the mass of the atoms. When the chain is stretched, $a$ becomes larger, which increases the group velocity. When the chain is compressed and the compression is not very much, the reduction of $a$ is apparently smaller than the increase of $a$ in the stretched case, and the bond angle changes more. The distance between the O–Si bond length (along the main force direction) in Fig. 5 and the figure of the O–Si–O angle (along the main force direction) in Fig. 4 indicate these changes. Therefore, there is an apparent increase in TC with increasing chain length. Kong[42,43] developed a method that is based on the theory of vibration scattering phonon theory to get the information of phonon. The method uses the displacement information of molecular dynamics to obtain dynamic matrix, and PDOSs are obtained by calculating the eigenvalues of the dynamic matrix. This method can accurately generate phonon distribution while taking into account the non-harmonic effects of phonon.

Fig. 6. Phonon DOSs of silicon atoms in the stretching/compressive siloxane structure and a partial enlargement of it.

Fig. 7. Peak frequency and TC in Fig. 6.

Fig. 8. Phonon DOSs of oxygen atoms in the stretching /compressive siloxane structure and a partial enlargement of it.

Figure 6 is the PDOS of silicon atoms in a silicon oxygen chain length of 10 nm at the different normalized chain lengths. The results at other lengths are almost similar, but not listed here. The longer the normalized length is, the higher the peak value of PDOS is. In addition, the corresponding frequencies of peaks for different normalized lengths are different. The higher peak frequency corresponds to the higher TC, as seen in Fig. 7. The reason behind this is the change in the MFP of phonon in the chain direction. The peaks will move to higher frequencies when the chain is stretched. The interaction force in the chain direction is enhanced and the transverse interaction becomes relatively weaker. The MFP in the chain direction will become larger under reduced transverse interaction, which is consistent with the change of MFP in positive Poisson's ratio materials.[15] Figure 8 is the PDOS of oxygen atoms. The results are similar to those of silicon atoms, which are consistent with the previous studies on polyethylene chains.[18] There is also an apparent size effect in TC as shown in Fig. 3. The silicone rubber with longer chain length has larger TC. The size effect on TC has been discussed in many reports[29,44] and will not be discussed here. In summary, we have applied the NEMD method to study the effects of stretching and compressing on the TC of crystalline silicone rubber. It is found that the TC gradually increase with the increase of the stretching length, which is consistent with the previous research results.[14] With the increase of compressing, the TC also tends to increase. The bond length and the bond angle on the skeleton show different variations with that off the skeleton. In addition, the distributions of bond length and bond angle under stretching and compressing are also analyzed. The PDOS shows that higher peak value generally corresponds to higher TC. The introduction of stretching strain in chain direction strengthens the force between atoms in this direction, moves the frequency peaks to higher value, and increases the MFP in this direction by reducing the transverse disturbance. References Synergetic effect of thermal conductive properties of epoxy composites containing functionalized multi-walled carbon nanotubes and aluminum nitrideThermal conductivity of silicone rubber filled with ZnOExperimental and theoretical studies of effective thermal conductivity of composites made of silicone rubber and Al2O3 particlesPolymer Nanowire Arrays With High Thermal Conductivity and Superhydrophobicity Fabricated by a Nano-Molding TechniqueHigh Thermal Conductivity of Single Polyethylene Chains Using Molecular Dynamics SimulationsDouble-slit photoelectron interference in strong-field ionization of the neon dimerMOLECULAR DYNAMICS SIMULATION ON THERMAL CONDUCTIVITY OF ONE DIMENISON NANOMATERIALSStrain Effect Analysis on Thermal Conductivity of Ge Thin FilmsEffects of tensile strain and finite size on thermal conductivity in monolayer WSe ₂Unusual thermal conductivity behavior of serpentine graphene nanoribbons under tensile strainMolecular Dynamics Calculations of the Thermal Conductivity of Molecular Liquids, Polymers, and Carbon NanotubesNonequilibrium molecular dynamics methods for computing the thermal conductivity: Application to amorphous polymersNon-equilibrium molecular dynamics calculation of heat conduction in liquid and through liquid-gas interfaceAn enhanced version of the heat exchange algorithm with excellent energy conservation propertiesPolymer Data Handbook, 2nd ed. Polymer Data Handbook, 2nd ed . Edited by James E. Mark (University of Cincinnati, OH). Oxford University Press: New York. 2009 . xii + 1250 pp. 195.00. ISBN 978-0-19-518101-2 .A local resonance mechanism for thermal rectification in pristine/branched graphene nanoribbon junctionsMolecular dynamics simulation of thermal energy transport in polydimethylsiloxaneAb initio calculations and force field development for computer simulation of polysilanesPolysiloxanes: ab initio force field and structural, conformational and thermophysical propertiesCOMPASS II: extended coverage for polymer and drug-like molecule databasesApplication of the G-JF discrete-time thermostat for fast and accurate molecular simulationsMolecular-dynamics study of a three-dimensional one-component model for distortive phase transitionsNumerical and experimental studies on the measurement of thermal conductivity of a silicone rubber plate as a reference materialPure spin current generated in thermally driven molecular magnetic junctions: a promising mechanism for thermoelectric conversionMolecular Dynamics Simulations of Thermal Properties of Solid Uranium DioxideMean free path spectra as a tool to understand thermal conductivity in bulk and nanostructuresPhonon dispersion measured directly from molecular dynamics simulationsImplementation of Green's function molecular dynamics: An extension to LAMMPSThermal conductivity of silicon nanowire by nonequilibrium molecular dynamics simulations

[1]	Mao L, Zhao D, Zou X, Wang J H and Shi L Y 2017 Insul. Mater. 8 9 (in Chinese)
[2]	Yang B C, Chen W Y, Zeng L and Hu Y D 2006 Annual Electronic Components Conference (Xining, China 15–19 August 2006) p 111 (in Chinese)
[3]	Teng C C, Ma C C, Chiou K C and Lee T M 2012 Compos. Part B 43 265
[4]	Wang J B, Bao Y B, Li Q Y and Wu C F 2012 Acta Mater. Composit. Sin. 29 6 (in Chinese)
[5]	Lei H J, Wang M L and Gong W F 2006 Henan Chem. Industry 23 20 (in Chinese)
[6]	Mu Q H, Feng S Y and Diao G Z 2007 Polym. Compos. 28 125
[7]	Gao B Z, Xu J Z, Peng J J, Kang F Y, Du H D, Li J, Chiang S W, Xu C J, Hu N and Ning X S 2015 Thermochim. Acta 614 1
[8]	Mou Q H, Feng W Y and Sun W Y 2008 China Patent No. 101284925A
[9]	Chen F 2008 China Patent No. 101168620A
[10]	Cao B Y, Kong J, Xu Y, Yung K L and Cai A 2013 Heat Transfer Eng. 34 131
[11]	Henry A and Chen G 2008 Phys. Rev. Lett. 101 235502
[12]	Xu Y F, Kraemer D, Song B, Jiang Z, Zhou J W, Loomis J, Wang J J, Li M D, Ghasemi H, Huang X P, Li X B and Chen G 2019 Nat. Commun. 10 1
[13]	Gao Y F and Meng Q Y 2010 Acta Metall. Sin. 46 1244
[14]	Lin Y P, Zhang M Y, Gao Y F, Mei L Y, Fu Y Z and Liu Y Q 2014 Acta Polym. Sin. 6 789
[15]	Wu J W, Tao Y, Chen C, Chen Y W and Chen Y F 2018 Southeast Univ. (Engl. Ed.) 34 43
[16]	Dong H K, Zhao D, Xu L, Li M B and Qi Y H 2018 J. Bohai Univ. 39 41 (in Chinese)
[17]	Zhang X L, Gong C Z and Wu G Q 2017 Rare Met. Mater. Eng. 46 370
[18]	Tang K 2015 MS thesis (Wuhan: Huazhong University of Science and Technology) (in Chinese)
[19]	Wen J H 2014 MS thesis (Xiangtan: Xiangtan University) (in Chinese)
[20]	Gong T 2013 MS thesis (Xiangtan: Xiangtan University) (in Chinese)
[21]	Yuan K P, Zhang X L, Li L and Tang D W 2019 Phys. Chem. Chem. Phys. 21 468
[22]	Gao Y, Yang W and Xu B 2016 Carbon 96 513
[23]	Feng X L, Li Z X and Guo Z Y 2001 J. Eng. Thermophys. 22 195
[24]	Elena A A and Florian M P 2012 Soft Mater. 10 42
[25]	Terao T, Lussetti E and Florian M P 2007 Phys. Rev. E 75 057701
[26]	Ikeshoji T and Hafskjold B 1994 Mol. Phys. 81 251
[27]	Wirnsberger P, Frenkel D and Dellago C 2015 J. Chem. Phys. 143 124104
[28]	Mark J E 1999 Polymer Data Handbook (New York: Oxford University Press) pp 417–424
[29]	Chen X K, Liu J, Xie Z X, Zhang Y, Deng Y X and Chen K Q 2018 Appl. Phys. Lett. 113 121906
[30]	Mao J H 2014 MS thesis (Chongqing: Chongqing university) (in Chinese)
[31]	Zhang M Y, Wang R H, Lin Y P, Li S M, Fu Y Z and Liu Y Q 2015 Polym. Mater. Sci. Eng. 31 68 (in Chinese)
[32]	Luo T F, Esfarjani K, Shiomi J, Henry A and Chen G 2011 J. Appl. Phys. 109 074321
[33]	Sun H 1995 Macromolecules 28 701
[34]	Sun H and Rigby D 1997 Spectrochim. Acta Part A 53 1301
[35]	Sun H, Jin Z, Yang C W, Akkermans R L C, Robertson S H, Spenley N A, Miller S and Todd S M 2016 J. Mol. Model. 22 47
[36]	Grønbechjensen N, Hayre N R and Farago O 2014 Comput. Phys. Commun. 185 524
[37]	Schneider T and Stoll E 1978 Phys. Rev. B 17 1302
[38]	Fujino J I, Honda T and Yamashita H 1997 Heat. Transf-Jap Res. 26 435
[39]	Wu D, Cao X H, Chen S Z, Tang L M, Feng Y X, Chen K Q and Zhou W X 2019 J. Mater. Chem. A 7 19037
[40]	Li J K and Tian X F 2010 Chin. Phys. Lett. 27 036501
[41]	Fan Y and Dames C 2013 Phys. Rev. B 87 035437
[42]	Kong L T 2011 Comput. Phys. Commun. 182 2201
[43]	Kong L T, Bartels G, Campañá C, Denniston C and Müser M H 2009 Comput. Phys. Commun. 180 1004
[44]	Wang S C, Liang X G, Xu X H and Ohara T 2009 J. Appl. Phys. 105 014316