Multi-scale modeling of thermal conductivity of SiC-reinforced aluminum metal matrix composite

Abstract

High thermal conductivity is one important factor in the selection or development of ceramics or composite materials. Predicting the thermal conductivity would be useful to the design and application of such materials. In this paper, a multi-scale model is developed to predict the effective thermal conductivity in SiC particle-reinforced aluminum metal matrix composite. A coupled two-temperature molecular dynamics model is used to calculate the thermal conductivity of the Al/SiC interface. The electronic effects on the interfacial thermal conductivity are studied. A homogenized finite element model with embedded thin interfacial elements is used to predict the properties of bulk materials, considering the microstructure. The effects of temperatures, SiC particle sizes, and volume fractions on the thermal conductivity are also studied. A good agreement is found between prediction results and experimental measurements. The successful prediction of thermal conductivity could help a better understanding and an improvement of thermal transport within composites and ceramics.

Keywords

Multi-scale modeling ab initio calculations molecular dynamics thermal conductivity composites

Introduction

Due to their excellent properties including high strength-to-weight ratio, good wear resistance, and high strength retained at elevated temperatures, ceramics-reinforced metal matrix composites are desired for many critical components in structures, engineering systems, and medical applications. Among these materials, silicon carbide (SiC) reinforcement is widely used due to its high thermal conductivity for various semiconductor devices, capsule materials for nuclear equipment, and gas seal rings in compressor pump.

Extensive experimental studies have been conducted to measure the thermal conductivity of SiC ceramics or SiC-reinforced metal matrix composites for different microstructure.^1–6 Considering the wide range of material compositions under various conditions, a systematic understanding of thermal transport is needed for SiC-reinforced metal matrix composites to reduce the cycle of trial-and-error methods. To this end, it is necessary to develop a predictive model to accurately calculate the thermal conductivity for various compositions at different temperatures. As the interface greatly influences the thermal conductivity, an appropriate way of calculating interfacial thermal conductivity is needed. One approach that has been widely used is the Hasselman and Johnson⁷ model. Hasselman and coworkers^6,7 considered the interfacial conductance between different constituents and proposed an analytical solution for the effective conductivity of particles-reinforced composites. Extensive studies have used the Hasselman–Johnson model to evaluate the effect of interfacial thermal conductance. With an assumed range of interfacial conductance, the effective thermal conductivity of SiC_f/SiC composites was studied.⁸ The Hasselman–Johnson model can further be used to calculate the interfacial conductance by fitting it to experimental results. By varying the size of the reinforcement particles,^6,9 the interfacial thermal conductance of particle-reinforced metal matrix composites can be calculated using the Hasselman–Johnson⁷ model. Using the derived inverse equation, Nan et al.¹⁰ calculated the interfacial thermal resistance from effective thermal conductivities, despite the implementation of calculation equations different from the Hasselman–Johnson model. Using finite element methods (FEM) modeling, Islam and Pramila¹¹ studied the effect of interfacial thermal conductance on the thermal conductivities of fiber-reinforced composites, where the interfacial thermal conductance was obtained by curve-fitting of experimental results.

Despite the convenience of the Hasselman–Johnson model in evaluating interfacial thermal conductance, it provides no understanding of interfacial thermal transport, such as temperature gradients, the size of interfacial region, and interfacial thermal conductivity. Additionally, the calculation accuracy largely depends on experimental measurements. Hence, a numerical method is needed to directly simulate the interfacial thermal transport and calculate the interfacial thermal conductivity.

Several recent studies have dealt with the interfacial thermal conductance by the molecular dynamics (MD) method.^12–16 Depending on whether the simulations are for conductivity under constant temperature or not, classical MD simulations can be classified into equilibrium molecular dynamics (EMD) and nonequilibrium molecular dynamics (NEMD) simulations. A well-known method of calculating the thermal conductivity of materials based on the atomic trajectories from EMD simulations is the Green–Kubo formula,^17,18 which correlates thermal conductivity to the heat current autocorrelation function. Different from EMD simulations, NEMD simulations are realized by imposing a heat flux or a temperature difference across a material system.^17–19 In NEMD simulations, the temperatures can be controlled by adding energy or rescaling temperatures of atoms. With the temperature gradient $\nabla T$ and heat flux q, the thermal conductivity k can be calculated according to the Fourier’s law. Comparing two methods discussed above, the thermal conductivity of structural defects, such as grain boundaries, can be relatively easily determined using the NEMD method by the suitable placement of heated and cooled regions. In controlling the system temperature of NEMD simulations, it is advantageous to use so-called reverse nonequilibrium molecular dynamics (RNEMD).⁷ In RNEMD simulations, the heat flux is imposed on the system, and the resulting temperature gradient is measured. The RNEMD could minimize the computation time since the calculation of heat fluxes takes much more time in the average manipulation.^20,21 Hence, the RNEMD method is used to compute the thermal conductivity in this study.

While MD simulations can model thermal transport, most of studies are focused on the phononic contribution due to the incapability of MD simulations in explicitly modeling the contribution of electrons. The two-temperature model (TTM) has been introduced to include electron–phonon coupling. The electron and phonon are described as two separate subsystems with an assigned temperature for each. To incorporate electronic effects into the modeling of radiation damage simulations, Duffy and Rutherford^22,23 introduced a coupled TTM-MD scheme to account for electron–phonon coupling, electronic stopping, and both the temporal and spatial evolution of the phononic and electronic subsystems. Based on Duffy and Rutherford’s scheme, Phillips and Crozier²⁴ added a communication process between the electronic and atomic subsystems and achieved energy-conserving TTM-MD simulations. Wang et al.¹³ extended the conventional NEMD simulation by coupling it with the TTM model to analyze the interfacial thermal transport across two different metal/nonmetal interfaces. A better agreement with experimental data was obtained compared with conventional NEMD simulations with only phonons considered. Using a coupled TTM-MD scheme, Jones et al.¹⁴ studied the size and electronic effects on the thermal conductance of Al/GaN interface. A linear relationship was found between the interfacial resistance and inverse system length while electron transport had little effect on the interfacial thermal conductance, which is also supported by the studies on several other metal/nonmetal interfaces.¹⁵ To reliably predict the thermal conductivity of SiC particle-reinforced aluminum matrix composites, the TTM-MD simulations were implemented in this study to evaluate the electronic contributions to the thermal transport of Al/SiC interface.

Despite the convenience of MD simulations in predicting thermal properties, the empirical interatomic potentials may not be available for MD simulations, which necessitates the need of ab initio calculations^25,26 to predict the interatomic potentials. Additionally, it is limited by the simulation system size and computation costs for bulk materials. Finite element (FE) modeling provides a convenient way to evaluate the effective thermal conductivity of bulk materials through homogenization of one representative volume element (RVE).^27–29 To fully take advantage of ab initio calculations, MD simulation, and FE analysis, it is necessary to develop a multi-scale model that can incorporate the benefits of calculations at each scale.

In this study, a multi-scale model is developed to predict the thermal conductivity of SiC particle-reinforced aluminum matrix composites by bridging ab initio calculations, MD simulations, and FEM. The interatomic potentials are derived by the ab initio method and subsequently used in MD simulations. The multi-scale model implements the temperature-dependent interfacial thermal conductivity calculated from TTM-MD simulations, where electronic effects are explicitly evaluated. The actual microstructures of composites are considered by varying SiC particle sizes and volume fractions in the FE analysis. The developed multi-scale model can not only implement the benefits of ab initio calculations and MD simulations in considering material microstructure at small scale and predicting interfacial thermal conductivity but also greatly reduces computation costs through FE analysis.

Multi-scale modeling of thermal conductivity

The derivation and selection of interatomic potentials

At different levels of temperatures, SiC particle-reinforced aluminum metal matrix composites of varying particle sizes and volume fractions are studied. To calculate material properties in MD simulations, it is very important to use accurate interatomic potential models. Without adding any additional element in this type of composites, only two constituents, i.e. SiC and Al, exist and hence, the interface between SiC particle and Al matrix is a mixture of SiC and Al.

Since little information is available for the measured potential data of typical elements in SiC-reinforced aluminum metal matrix composites, i.e. Al-Si/C in this study, interatomic potentials need to be determined using a predictive method. Ab initio methods provide a convenient way of obtaining interatomic potentials. A predefined interatomic potential based on density functional theory is usually used with curve-fitted parameters.^30,31 By fitting the total energy of interatomic potentials to the results by ab initio calculation, the parameters of the interatomic potential model can be determined.

To calculate the interatomic potentials, the first step is to get the adhesive energies of the established system as the source data. Then, an analytical formula can be derived by using the inversion method,³² which expresses the interatomic potential in terms of adhesive energies. It is assumed that the adhesive energy is equal to the sum of interactions across the interface, which is a first-order but practical approximation for the metal/semiconductor interface. Because the charge transfer at the interface is mainly confined to the first Al and SiC layers,³³ the three-body interaction j–i–k is limited in this region, which means the central atom i is located in the first SiC layer, and j, k are chosen from the neighboring Al and SiC layers. For the three-body interactions, the modified Stillinger–Weber (MSW) potential form has been widely used in covalent matters with zinc blende structures, which is the typical interface structure of SiC ceramics. The MSW form describing three-body interaction is shown as

Φ_{jik} (r_{ij}, r_{ik}, θ_{jik}) = λ_{jik} \exp (\frac{γ_{ij}}{r_{ij} - R_{ij}} + \frac{γ_{ik}}{r_{ik} - R_{ik}}) \times (\cos θ_{jik} - \cos θ^{0})^{2}

(1)

where

λ_{jik}, γ_{ij}, R_{ij}, γ_{ik}, R_{ik}

, and

θ^{0}

are the parameters that need to be determined. More details about potential derivation could be referred to the work done in the literature.^32,33

Based on the interatomic potential derivation, ab initio calculation of adhesive energy of the interface structure was carried out using a CASTEP program based on the Al/SiC interface structure. A CASTEP program with generalized gradient approximation was used to calculate the adhesive energy of the interface structure, where a plane-wave cutoff energy of 400 eV was used with a K-point mesh of 7 × 7 × 3. To validate the accuracy of the ab initio calculation, the interatomic potentials of Al-Si and Al-C pairs were calculated in order to compare with the ab initio calculation results of Al-SiC (100) interface in literature.³² As can be seen from the comparison in Figure 1, the calculated interatomic potential results agree very well with the data obtained in the literature.

Figure 1.

Interatomic potentials comparison between results by Zhao and Chen³² and ab initio validation results.

In the selection of interatomic potential form, it was found that the popular Morse potential could fully represent the ab initio data of the Al–Si and Al–C interfaces.

V = D_{0} [e^{- 2 α (r - r_{0})} - 2 e^{- α (r - r_{0})}]

(2)

where r is the distance between the atoms, r₀ is the equilibrium bond distance, D₀ is the well depth of the potential, α is the width of the potential, and V is the interatomic potential. The derived parameters by curve-fitting the function to ab initio data are summarized in Table 1. For three-body interatomic potentials, the parameters shown in Table 2 were determined for the MSW potential equation by curve-fitting the contribution of three-body potentials to the ab initio calculation results for the hex-site configuration.³³

Table 1.

Morse potential function parameters parameterized to the ab initio data.

Pair interaction	Parameters	Morse potential
Al-Si	$D_{0} (eV)$	0.320
	α(1/Å)	1.532
	r₀Å	2.820
Al-C	$D_{0} (eV)$	0.992
	α(1/Å)	1.738
	r₀(Å)	2.100

Table 2.

Stillinger–Weber potential parameters for the Al/SiC interface.

$λ_{Al - Si - C}$ (eV)	$γ_{Al - Si}$ (Å)	$R_{Al - Si}$ (Å)	$λ_{Al - Si - Al}$ (eV)	$θ_{Al - Si - Al}^{0}$ (°)
0.9203	0.1760	3.78	0.0362	181.0210
$λ_{Al - C - Si}$ (eV)	$γ_{Al - C}$ (Å)	$R_{Al - C}$ (Å)	$λ_{Al - C - Al}$ (eV)	$θ_{Al - C - Al}^{0}$ (°)
0.2196	0.0613	3.39	0.5099	2.2609

SiC: silicon carbide.

The embedded-atom method (EAM) as shown in equation (3) was used to model aluminum in the MD simulations through a summation of pairwise interaction and the energy required to embed an atom in the electron density produced by the local site.³⁴

E = \frac{1}{2} \sum_{i, j, i \neq j} φ_{ij} (r_{ij}) + \sum_{i} F_{i} (ρ_{i})

(3)

where

φ_{ij}

represents the pair energy between atoms i and j seperated by the distance

r_{ij}

. F_i denotes the embedding energy required to embed an atom i into a local site with electron density

ρ_{i}

. It should be noted that the EAM describes the dynamics of atoms without explicitly considering the constribution of the free electrons. However, as the parameters in the EAM model are fitted to properties of the real metal, part of the electronic constribution to the bulk thermodynmic properties is included in the simulations,³⁵ hence yielding a larger thermal conductivity compared to MD simulations without considering electrons. To fully consider the electronic contribution in the thermal transport, the TTM model can explicitly describe the electron–phonon coupling effects and hence is implemented here.

Tersoff potentials were implemented in the MD simulation to represent the atomic interactions within SiC. The necessary constants and corresponding values for aluminum and silicon carbide could be found from the literature^34,36,37 and are summarized in Tables 3 and 4, respectively.

Table 3.

EAM potential parameters for aluminum.

Constant	Value	Constant	Value
r_e	2.886166	$F_{n 0}$	−2.806783
f_e	1.392302	$F_{n 1}$	−0.276173
$ρ_{e}$	20.226537	$F_{n 2}$	0.893409
α	6.942419	$F_{n 3}$	−1.637201
β	3.702623	F ₀	−2.83
A	0.251519	F ₁	0
B	0.313394	F ₂	0.929508
κ	0.395132	F ₃	−0.68232
λ	0.790264	η	0.779208
		F_e	−2.829437

EAM: embedded-atom method.

Table 4.

Tersoff potential parameters for silicon carbide.

	C	$Si$	$C - Si$
A	1393.6	1830.8	1597.3
B	346.7	471.18	404.20
λ	3.4879	2.4799	2.9839
μ	2.2119	1.7322	1.9720
β	1.5724 × 10 ^− 7	1.1000 × 10 ^− 6	–
n	7.2751 × 10 ^− 1	7.8734 × 10 ^− 1	–
c	3.8049 × 10⁴	1.0039 × 10⁵	–
d	4.384	16.217	–
h	−057058	−0.59825	–
R	1.8	2.7	2.2
S	2.1	3.0	2.5
χ	1.0	1.0	9.9776

MD equations of thermal conductivity in solid materials

The RNEMD algorithm comprises two components: a means of moving energy from the cold to the hot slab and a way of calculating the temperature gradient. For the calculation of thermal conductivity, the numerical system of heat conduction in a solid material is considered as shown in Figure 2. In the heating and cooling layers, the atoms are heated and cooled, respectively, and the heat is conducted through the atoms between them, i.e. conduction layer.

Figure 2.

Configuration of heat conduction layers.³⁸

The heat flux is given by the sum of all energy transfer per time per area

q = \frac{1}{2 tA} \sum \frac{m}{2} (v_{heat}^{2} - v_{cool}^{2})

(4)

where t is the length of the simulation, A is the cross-sectional area of the simulation box perpendicular to flow direction,

v_{heat}

and

v_{cool}

denote the velocities of the heating and cooling layers, whose velocities are exchanged, and the factor 2 arises from the periodicity. In RNEMD, by rescaling particle velocities at each time step, heat is added into the heating layer and removed from the cooling layer. The result of such a process is that a constant energy current, i.e. heat flux q, is given.

To calculate the thermal conductivity, the second task is the calculation of the temperature gradient. The temperature T within each layer can be evaluated by

\frac{3 {Nk}_{B} T}{2} = \sum \frac{m}{2} v^{2}

(5)

where the number of atoms in the layer is N, and k_B is Boltzmann constant. The temperatures

T_{heat}

and

T_{cool}

in the heating and cooling layers could be determined as

T_{heat} = \frac{2}{3 N_{heat} k_{B}} \sum \frac{m}{2} v_{heat}^{2}

(6)

T_{cool} = \frac{2}{3 N_{cool} k_{B}} \sum \frac{m}{2} v_{cool}^{2}

(7)

A linear response in the material system could be assumed as can be seen for all cases^17,18,20,39 in the MD calculation of thermal conductivities. Hence, the temperature gradient in the heat conduction direction can be assumed as

\nabla T = \frac{T_{heat} - T_{cool}}{L}

(8)

Combining equations (4)–(8) with the Fourier’s law, an approximation of thermal conductivity can be obtained as

k = \frac{q}{\nabla T} = \frac{qL}{(1 - \frac{N_{heat}}{N_{cool}}) T_{heat} + \frac{4 tA}{3 N_{cool} k_{B}} q}

(9)

In equation (9), L and A are related to the dimensions of the material system, and t is determined by the simulation time. These parameters are pre-determined and remain the same once the simulation setup is fixed.

In order to induce heat conduction through the atom layers, there are two types of thermal systems for the heating and cooling layers: constant temperature and constant heat flux. In the system of constant temperature, the heating and cooling layers are kept at constant high and low temperatures, respectively, whereas the constant flux system maintains a constant heat flux at the heating layer and subtracts the heat flux of the same amount from the cooling layer. In the current study, the latter method is employed to minimize computation time since the calculation of heat fluxes takes much time in the averaging manipulation. As can be seen in equation (9), after the heat flux q is given, the thermal conductivity k is determined by $N_{heat}$ and $N_{cool}$ , which are determined by the composition of the material system, and the material temperature $T_{heat}$ . Hence, the thermal conductivity k is related to the material composition and temperature.

The initial velocities of the simulation system are determined from the desired temperature for thermal conductivity calculation. The periodic boundary conditions are applied to all directions of the simulation system. The heat flux is added and subtracted by rescaling the kinetic energy of the heating and cooling layers, respectively. A constant-pressure-and-temperature ensemble first is used to fully relax the simulation system. The simulation is subsequently run for 2 × 10⁶Δt to reach a steady heat flow across the simulations system in a constant-volume-and-energy ensemble.

Two-temperature nonequilibrium MD simulation

In this study, the TTM-MD simulation scheme is implemented to perform nonequilibrium MD simulation for phonon degree of freedom, which is coupled with finite difference (FD) calculation for the electron degree of freedom in the LAMMPS package.⁴⁰

TTM depicts the coupled electronic and phononic thermal transport through separated temperature fields and a shared coupling term.¹³ The temporal and spatial evolution of temperature fields can be described by two coupled heat diffusion equations.

ρ_{e} c_{e} \frac{\partial T_{e}}{\partial t} = \nabla (k_{e} \nabla T_{e}) - g_{ep} (T_{e} - T_{p}) + g_{s} T_{p}'

(10)

ρ_{p} c_{p} \frac{\partial T_{p}}{\partial t} = \nabla (k_{p} \nabla T_{p}) + g_{ep} (T_{e} - T_{p}) - g_{s} T_{p}'

(11)

where ρ, c, T, and k represent the density, specific heat, temperature, and thermal conductivity of electrons (subscript e) and phonons (subscript p), respectively.

g_{ep}

is the electron–phonon coupling parameter, while g_s represents the electronic stopping effect. The temperature term

T_{p}'

denotes the temperature of phonons with velocity greater than the threshold for the electron stopping interaction.

The TTM-MD simulation coupled electronic and phononic subsystems simultaneously as shown in Figure 3. Phononic heat diffusion (equation (11)) is still modeled by the MD technique. As MD simulations only need the inputs of the atomic structure and interatomic potentials, k_p with its dependence on temperature and composition is implicitly calculated by MD. The electronic subsystem is modeled by solving equation (10) iteratively with the FD method in the TTM-MD simulations. Phonons and electrons are coupled via the terms $g_{ep} (T_{e} - T_{p})$ and $g_{s} T_{p}'$ in equations (10) and (11).

Figure 3.

Illustration of the coupled electronic and phonic subsystems in TTM-MD simulations. Dots represent atoms in MD simulations, and lines show the FD grid in solving equation (10) using the FD method.

In the implementation of the TTM-MD simulation,⁴⁰ the equation of motion for an atom i in the MD simulation is in the form of a Langevin thermostat.¹³

m_{i} \frac{\partial v_{i}}{\partial t} = F_{i} (t) - γ_{i} v_{i} + {\tilde{F}}_{i} (t)

(12)

where m_i,

v_{i}

, and

F_{i}

are the atomic mass, velocity, and total force on atom i, respectively. The total force

F_{i}

is evaluated by the interatomic potentials. The friction term

γ_{i}

represents the energy loss by electron-ion interactions

γ_{p}

and electronic stopping

γ_{s}

, which is related to

g_{ep}

and g_s in the electronic subsystem as

γ_{p} = \frac{m_{i} g_{ep}}{3 n_{i} k_{B}}

(13)

γ_{s} = \frac{m_{i} g_{s}}{3 n_{i} k_{B}}

(14)

where n_i and k_B denote the atom number density and the Boltzmann constant, respectively.

{\tilde{F}}_{i} (t)

is a random force in Langevin thermostat. For the TTM-MD simulation in this study

{\tilde{F}}_{i} (t) = \sqrt{\frac{24 k_{B} T_{e} γ_{i}}{Δ t}} {\tilde{R}}_{i}

(15)

where

Δ t

denotes the temporal interval of the simulation, and

{\tilde{R}}_{i}

is a random vector. The coupled electronic subsystem and phononic subsystems are equilibrated to a shared temperature through the last two terms in equation (12).

In the TTM-MD simulation, the MD part only needs the inputs of atomic structure and interatomic potentials, which can be either measured in the experiments or derived through ab initio methods as discussed in Section “The derivation and selection of interatomic potentials.” On the contrary, the implementation of TTM requires several input parameters as shown in equation (10), i.e. c_e, k_e, $g_{ep}$ , and g_s. One way to determine c_e, $g_{ep}$ , and g_s is the derivation from ab initio methods. Lin et al.⁴¹ conducted extensive studies on their temperature dependence for various metals from ab initio calculations, and a linear relation was found for a temperature below 1000 K. The calculated results for aluminum were implemented here in the TTM simulations. The Wiedemann–Franz law can be used to approximate k_e based on the measured Lorenz number⁴² and electrical conductivity⁴³ through extensive experimental studies.

Multi-scale modeling of thermal conductivity in FE framework

A homogenized 2D FE model is developed to calculate the effective thermal conductivity. The 2D FE model is able to calculate the effective thermal conductivity as a function of temperatures and compositions. The bulk material is represented by continuum elements with embedded thin interfacial elements based on the microstructure characteristics. The thermal conductivity could be evaluated from the FE calculation performed on one RVE as shown in Figure 4. Such a homogenization procedure was used to estimate the effective thermal conductivity of bulk materials, which include large amounts of polycrystals with more details about the FE multi-scale model in the previous study.⁴⁴ The SiC particle-reinforced aluminum matrix composites were studied in this paper. More experimental details with respect to SiC particle size and temperature could be referred to the literature.⁶ The obtained SiC particle was angular in shape. To facilitate the simulations, spherical SiC particle was assumed as shown in Figure 4. The SiC particle size was represented by the diameter of spherical SiC particle in the multi-scale model. Relatively uniform SiC particle size was found in the experiments for the composites studied here.⁶ Our preliminary studies showed that even after considering the uneven distribution of SiC particles in the actual materials,⁶ the predicted thermal conductivities still lied within the uncertainties of simulation results shown in this study. Hence, a uniform diameter of SiC particles was modeled in accordance with the observation from the experiments.⁶

Figure 4.

Multi-scale modeling of SiC-reinforced aluminum matrix composites with varying particle size.

A uniform heat flux was applied at each point of the boundary in Abaqus. For RVE, calculations in the heat flow direction were performed to determine the thermal conductivity, and the thermal conductivity is calculated as $κ = q / (\partial T / \partial z) A$ , where q represents input heat flux, $\partial T / \partial z$ shows temperature gradient, and A denotes the cross-sectional area.

Results and discussion

System size effects in thermal conductivity calculation

One important factor in calculating thermal conductivity is the finite-size effect in MD simulations.^14,45 When the length of the simulation cell is not significantly longer than the phonon mean-free path, the actual thermal conductivity will be affected by the system size.⁴⁵ This is caused by scattering that occurs at the interfaces with the heat source and sink. To avoid such effect in calculating the interfacial thermal conductivities, studies on the effect of material system size were first performed.

A typical time-average temperature profile used to calculate the thermal conductivity is shown in Figure 5. Within the regions near the heat source or heat sink, nonlinear temperature profiles are observed, which are attributed to the strong scattering caused by the heat source or heat sink.⁴⁵ In the intermediate regions, the temperature profiles are fitted with a linear function to obtain the temperature gradients, which are further used to calculate the thermal conductivity by Fourier’s law. The differences in the computed gradients for the different regions are used as an error estimate for the value of thermal conductivity.

Figure 5.

Typical temperature profile for 16 × 16 × 2000 Å Al system at room temperature.

The dependence of k on the dimensions perpendicular to the heat flow, i.e. x/y direction, is first studied. It is expected that k will not depend on the dimensions perpendicular to the heat flow as sensitively as in the direction of the length L_z. Due to periodic boundary conditions on the simulation cell, phonons are free to travel across the simulation cell perpendicular to the heat flow direction without scattering from any boundaries. Hence, changing the dimension perpendicular to the heat flow does not change the scattering in any obvious way. Figure 6 shows a comparison of k obtained for systems of different sizes in the direction perpendicular to the heat flow. The predicted k is almost the same for all the cases. As no obvious dependence on x/y dimensions is observed, the smallest x/y dimensions shown here are used for the rest of studies, i.e. 16 Å.

Figure 6.

Comparison of k as for different cell sizes perpendicular to the direction of the heat flow. Each system was 400 Å long parallel to the heat flow in MD simulations. The calculated results at the room temperature for aluminum were shown for demonstration. Similar trend was observed for SiC system.

Simulations to evaluate the size dependency of the simulation box in the dimension parallel to the heat flow L_z on the results were conducted. The thermal conductivities of the Al were calculated for simulation systems with varying length L_z. The simulated system sizes and the corresponding k are summarized in Table 5. As the system size increases to 16 × 16 × 2000 Å, the finite-size effect of the material system decreases, and the simulated thermal conductivity of Al tends to converge to the experimental observation⁴⁶ of

237 W / mK

. The simulation of a larger system size, i.e. 16 × 16 × 2500 Å, also yields a value close to the experimental observation. To reduce computation costs, the smaller system size, 16 × 16 × 2000 Å, is selected for the rest of studies.

Table 5.

System size dependency of the thermal conductivity of Al.

Case	System size (Å)	Thermal conductivity ( $W / mK$ )
1	16 × 16 × 400	27.3 ± 2.0
2	16 × 16 × 1500	193.9 ± 2.7
3	16 × 16 × 2000	237.6 ± 7.5
4	16 × 16 × 2500	236.8 ± 3.7

Calculations of interfacial thermal conductivity across the Al/SiC interface

To verify the capabilities of the TTM-MD simulations in modeling coupled electronic and phononic thermal transport across the Al/SiC interface, the thermal conductivity of aluminum was first calculated and compared with experimental measurements. Both conventional MD and TTM-MD simulations were conducted. Thermal conductivity is obtained according to the Fourier’s law of heat conduction. The experimental value and simulation results of k are summarized in Figure 7. It is worth noting that as electrons will not contribute much to thermal transport within SiC and Si, MD simulations without electron effects are sufficient to predict relatively accurate thermal conductivities for SiC and Si.¹³ In contrast, the contribution of electrons to thermal transport cannot be ignored for Al, wherein many free electrons exist. k from MD simulations without considering electrons for Al is $13.5 W / mK$ , which only accounts for the phononic thermal transport. The TTM-MD model predicts k = 237.6 W/mK for Al, which is very close to the experimental value and suggesting the capability of the TTM-MD model in modeling the electronic contribution to k.

Figure 7.

Thermal conductivity comparison between experimental measurement^46–48 and prediction results at the room temperature.

The thermal conductivities of SiC and Si at room temperature were also calculated using the MD simulation method. The MD calculation results were compared with the experimental data in order to validate the proper selection of interatomic potentials. Figure 7 shows the comparison of the predicted thermal conductivities with the experimental data.^46,47 It is evident that MD simulations match very well with the experimental measurement with average prediction errors well below 3.0%, thus confirming the validity of interatomic potential selection.

The thermal conductivities of Al/SiC interface at different temperatures were also calculated using the TTM-MD simulation method discussed above. For the NEMD method, it is important to establish a steady-state heat flow. This amounts to obtaining a stationary temperature profile as a function of time, thus ensuring that only steady-state heat flux is flowing. The typical temperature profile of Al/SiC structure is shown in Figure 8. The interface region is defined as the region between Al and SiC based on linear fitting results.⁴⁹ The existence of interface contributed to the significant temperature drop between Al and SiC by acting as a thermal barrier for the heat transport between Al and SiC. From the temperature profiles, the temperature drop at the interface, ΔT, can be obtained, which can be subsequently used to calculate the interfacial thermal resistance

R_{I} = \frac{Δ T}{q}

(16)

Figure 8.

A typical temperature profiles of Al/SiC interface structure at room temperature.

Figure 9 shows a linear relationship between interfacial resistance and inverse system length L_z, as supported by the literature.¹⁴ At the system limit, $L_{z} = L_{Al} + L_{SiC} \to \infty$ , the interfacial resistance should be independent of system size and thus can be used for the interfacial resistance of bulk materials.^14,49 The predicted interfacial thermal resistance at room temperature is 1.01 $m^{2} K / GW$ , which yields a thermal conductance of $9.86 \times 10^{8} W / m^{2} K$ . The predicted thermal conductance is apparently larger than the measured value in the literature.^6,50 However, the discrepancy lies within the uncertainties typically observed between the measured value and the predicted thermal conductance by MD simulations in the literature.^13,27,51,52 It is worth noting that despite the relatively large discrepancy found in MD simulations, the predicted interfacial thermal properties still can help understand the thermal transport at the atomic scale and further calculate effective thermal conductivities of composites under various conditions. In contrast, it is usually relatively difficult to examine thermal transport and measure interfacial properties at all conditions in the experiments. As ideal Al/SiC interface was selected in the simulations of this study, the defects like dislocations,²⁷ misfit grain boundaries, and voids⁵³ may exist within the actual interface and increase phonon scattering, hence lowering the actual interface thermal conductivity. The electronic contribution to the thermal transport of Al/SiC interface can also be evaluated by comparing the predicted interfacial resistance by TTM-MD simulations with those by MD simulations. While MD simulations only considered the contributions of phonons to the interfacial thermal transport, TTM-MD added the electronic contribution in the simulations. As can be seen in Figure 9, the predicted values by TTM-MD methods were very close to the results predicted by MD simulations in all three cases, showing a small effect of electronic contribution to interfacial thermal transport. These findings were also observed in several other studies of Al/nonmetal interface structures.^14–16

Figure 9.

Al/SiC interfacial resistance as a function of inverse system length at room temperature.

It is worth noting that the interfacial thermal conductivity was also calculated from Fourier’s law as required by the developed multi-scale FE model in this study. The experimental values for Al and SiC are extracted from the measured results from the literature,⁶ which are expected to be lower than those of pure Al and SiC materials (as discussed above) due to the existence of impurity in the actual materials. Both the MD results for the interface and the experimental results for Al and SiC are plotted in Figure 10. The interfacial thermal conductivity decreases with increasing temperature as the effect of phonon scattering increases. Due to the increased phonon scattering caused by the interface, the thermal conductivity of interfacial phase is apparently much lower than those of Al and SiC, which is the reason for the much lower thermal conductivity of SiC-reinforced aluminum matrix composites compared with each constituent, i.e. Al matrix and SiC particle, due to the existence of interface between them.

Figure 10.

Calculated interfacial thermal conductivities from MD simulations and comparison with experimental measurements of SiC and Al.⁶ The inset is an enlargement of the predicted thermal conductivity for Al/SiC interface with respect to temperature.

Temperature effects on the thermal conductivities of SiC-reinforced aluminum matrix composites

After interfacial thermal conductivities were obtained, the effective thermal conductivities of SiC particle-reinforced composites at different temperatures were calculated using the developed multi-scale model in Section “Multi-scale modeling of thermal conductivity in FE framework.” The temperature-dependent thermal conductivities for both SiC particles and Al matrix⁶ were also considered in the simulations. The calculated results are shown in Figure 11, where it can be seen that the simulated results of the multi-scale model yield a maximum 6.21% deviation compared to experimental results. The predicted results are consistently higher than the experimental measurements, which can be explained by the defects and voids existing in the actual composites. Since only ideal composite microstructure was considered in the simulations, the existence of defects and voids can lead to phonon scattering, lowering the thermal conductivity of the actual material.

Figure 11.

Comparison between effective thermal conductivities and experimental measurements ⁶ for 40 vol% particulate-SiC-reinforced aluminum matrix composites. The SiC particle size is 10.2 µm in diameter.

Effect of SiC reinforcement particle size on thermal conductivities of aluminum matrix composites

The effects of SiC particle size on the effective thermal conductivity were studied at room temperature. As shown in Figure 12, calculations of effective thermal conductivities show an increase in effective thermal conductivities with increasing SiC particle size as found experimentally by Hasselman et al.⁶ While particle sizes increased, the SiC volume fraction remained the same as 40 vol% in the composite. The prediction results show a good agreement with the experimental measurements as shown in Figure 12.

Figure 12.

Comparison between effective thermal conductivities and experimental measurements⁶ for 40 vol% SiC-reinforced aluminum matrix composites with different SiC particle size at room temperature.

The increasing thermal conductivity with increasing particle size can be explained by the fact that with increasing SiC particle size and associated decrease of total interfacial area, the relative contribution of interfacial thermal barrier to the total conductivity of the composite decreases, thus yielding an increasing effective thermal conductivity. The maximum discrepancy between experimental measurements and prediction results is within 9.88%, thus validating the capability of the multi-scale modeling in predicting the effective thermal conductivity.

Effect of SiC volume fraction on thermal conductivities of aluminum matrix composites

Studies on the effect of SiC volume fraction were also carried out using the developed multi-scale model. A nominal particle size of 37 µm⁵⁴ was implemented in modeling SiC-reinforced aluminum matrix composites. SiC particles were uniformly distributed in the aluminum metal matrix composite with varying SiC volume fractions from 5% to 25%. The variations of effective thermal conductivities are plotted in Figure 13. It can be seen in Figure 10 that the thermal conductivity of SiC is relatively larger than that of Al. As the SiC volume fraction increases, its relative contribution to the total conductivity of the composites increases, hence enhancing the effective thermal conductivity of the composites. As seen in Figure 13, the prediction results agreed well with experimental measurements with differences being smaller than 7.48%.

Figure 13.

Effect of SiC volume fraction on effective thermal conductivities and comparison with experimental results.⁵⁴

Conclusion

A multi-scale model has been developed to predict the thermal conductivity of SiC-reinforced aluminum metal matrix composites. An ab initio method was used to calculate the interatomic potential for MD simulations. The interfacial thermal conductivity was calculated in the TTM-MD simulations, which was subsequently used in the proposed multi-scale model to predict the thermal conductivity of bulk materials. The model successfully predicted the decreasing thermal conductivity with increasing temperature. A small effect of electronic contribution to the thermal transport of the Al/SiC interface was observed. It has also been shown that the increase of SiC particle size and volume fraction can increase the thermal conductivity due to the reduced interfacial thermal barrier and the increased contribution of SiC reinforcement, respectively. The prediction results agreed well with the measurements for different cases, showing a maximum of 9.88% discrepancy. The successful predictions of thermal conductivity proved the efficacy of the developed multi-scale model, which can be used to design thermal properties of ceramics and composites with different compositions.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Kim

Y-W

Lim

K-Y

Seo

W-S

. Microstructure and thermal conductivity of silicon carbide with yttria and scandia. J Am Ceram Soc 2014; 97: 923–928.

Nakano

Watari

Kinemuchi

et al.

Microstructural characterization of high-thermal-conductivity SiC ceramics. J Eur Ceram Soc 2004; 24: 3685–3690.

Taguchi

Igawa

Yamada

et al.

Effect of thick SiC interphase layers on microstructure, mechanical and thermal properties of reaction-bonded SiC/SiC composites. J Phys Chem Solids 2005; 66: 576–580.

Watari

Nakano

Sato

et al.

Effect of grain boundaries on thermal conductivity of silicon carbide ceramic at 5 to 1300 K. J Am Ceram Soc 2003; 86: 1812–1814.

Zhou

Hirao

Watari

et al.

Thermal conductivity of silicon carbide densified with rare-earth oxide additives. J Eur Ceram Soc 2004; 24: 265–270.

Hasselman

DPH

Donaldson

Geiger

. Effect of reinforcement particle size on the thermal conductivity of a particulate-silicon carbide-reinforced aluminum matrix composite. J Am Ceram Soc 1992; 75: 3137–3140.

Hasselman

DPH

Johnson

. Effective thermal conductivity of composites with interfacial thermal barrier resistance. J Compos Mater 1987; 21: 508–515.

Youngblood

Senor

Jones

et al.

The transverse thermal conductivity of 2D-SiCf/SiC composites. Compos Sci Technol 2002; 62: 1127–1139.

Kida

Weber

Monachon

et al.

Thermal conductivity and interfacial conductance of AlN particle reinforced metal matrix composites. J Appl Phys 2011; 109: 064907–064907.

10.

Nan

C-W

X-P

Birringer

. Inverse problem for composites with imperfect interface: determination of interfacial thermal resistance, thermal conductivity of constituents, and microstructural parameters. J Am Ceram Soc 2000; 83: 848–854.

11.

Islam

Pramila

. Thermal conductivity of fiber reinforced composites by the FEM. J Compos Mater 1999; 33: 1699–1715.

12.

Bagri

Kim

S-P

Ruoff

et al.

Thermal transport across twin grain boundaries in polycrystalline graphene from nonequilibrium molecular dynamics simulations. Nano Lett 2011; 11: 3917–3921.

13.

Wang

Ruan

Roy

. Two-temperature nonequilibrium molecular dynamics simulation of thermal transport across metal-nonmetal interfaces. Phys Rev B 2012; 85: 205311–205311.

14.

Jones

Duda

Zhou

et al.

Investigation of size and electronic effects on Kapitza conductance with non-equilibrium molecular dynamics. Appl Phys Lett 2013; 102: 183119–183119.

15.

Singh

Seong

Sinha

. Detailed consideration of the electron-phonon thermal conductance at metal-dielectric interfaces. Appl Phys Lett 2013; 102: 181906–181906.

16.

Zhou

X-W

Jones

Hopkins

et al.

Thermal boundary conductance between Al films and GaN nanowires investigated with molecular dynamics. PCCP 2014; 16: 9403–9410.

17.

Wang

Shimono

et al.

Molecular dynamics simulation of thermal conductivity of silicon thin film. Mater Trans 2007; 48: 2419–2421.

18.

Jund

Jullien

. Molecular-dynamics calculation of the thermal conductivity of vitreous silica. Phys Rev B 1999; 59: 13707–13707.

19.

Maruyama

. Molecular dynamics method for microscale heat transfer, New York: Taylor & Francis, 2000, pp. 189–226.

20.

Müller-Plathe

. A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity. J Chem Phys 1997; 106: 6082–6085.

21.

Kotake

Wakuri

. Molecular dynamics study of heat conduction in solid materials. JSME Int J, Ser B 1994; 37: 103–108.

22.

Duffy

Rutherford

. Including the effects of electronic stopping and electron–ion interactions in radiation damage simulations. J Phys Condens Matter 2007; 19: 016207–016207.

23.

Rutherford

Duffy

. The effect of electron–ion interactions on radiation damage simulations. J Phys Condens Matter 2007; 19: 496201–496201.

24.

Phillips

Crozier

. An energy-conserving two-temperature model of radiation damage in single-component and binary Lennard-Jones crystals. J Chem Phys 2009; 131: 074701–074701.

25.

Dmitriev

Yoshikawa

Kohyama

et al.

Modeling interatomic interactions across Cu/α–Al2O3 interface. Comput Mater Sci 2006; 36: 281–291.

26.

Gao

Deng

Heinisch

et al.

A new Fe–He interatomic potential based on ab initio calculations in α-Fe. J Nucl Mater 2011; 418: 115–120.

27.

Crocombette

J-P

Gelebart

. Multiscale modeling of the thermal conductivity of polycrystalline silicon carbide. J Appl Phys 2009; 106: 083520–083520.

28.

Mortazavi

Hassouna

Laachachi

et al.

Experimental and multiscale modeling of thermal conductivity and elastic properties of A/expanded graphite polymer nanocomposites. Thermochim Acta 2013; 552: 106–113.

29.

Shin

Yang

Chang

et al.

Multiscale homogenization modeling for thermal transport properties of polymer nanocomposites with Kapitza thermal resistance. Polymer 2013; 54: 1543–1554.

30.

Dmitriev

Yoshikawa

Kohyama

et al.

Modeling interatomic interactions across Cu/α–Al 2 O 3 interface. Comput Mater Sci 2006; 36: 281–291.

31.

Gao

Deng

Heinisch

et al.

A new Fe–He interatomic potential based on ab initio calculations in α-Fe. J Nucl Mater 2011; 418: 115–120.

32.

Zhao

Chen

. An inverse adhesion problem for extracting interfacial pair potentials for the Al(0 0 1)/3C–SiC(0 0 1) interface. Inverse Prob 2008; 24: 035019–035019.

33.

Zhao

Chen

Long

. Interfacial potentials for Al/SiC(111). J Phys Condens Matter 2009; 21: 225002–225002.

34.

Zhou

Wadley

HNG

Johnson

et al.

Atomic scale structure of sputtered metal multilayers. Acta Mater 2001; 49: 4005–4015.

35.

Richardson

Clancy

. Contribution of thermal conductivity to the crystal-regrowth velocity of embedded-atom-method-modeled metals and metal alloys. Phys Rev B 1992; 45: 12260–12268.

36.

Tersoff

. Carbon defects and defect reactions in silicon. Phys Rev Lett 1990; 64: 1757–1760.

37.

Tersoff

. Modeling solid-state chemistry: interatomic potentials for multicomponent systems. Phys Rev B 1989; 39: 5566–5568.

38.

Kotake

Wakuri

. Molecular dynamics study of heat conduction in solid materials. JSME Int J Ser B 1994; 37: 103–108.

39.

Thomas

Turney

Iutzi

et al.

Predicting phonon dispersion relations and lifetimes from the spectral energy density. Phys Rev B 2010; 81: 081411–081411.

40.

Plimpton

. Fast parallel algorithms for short-range molecular dynamics. J Comput Phys 1995; 117: 1–19.

41.

Lin

Zhigilei

Celli

. Electron-phonon coupling and electron heat capacity of metals under conditions of strong electron-phonon nonequilibrium. Phys Rev B 2008; 77: 075133–075133.

42.

Olafsson

Sandstrom

Karlsson

. Comparison of experimental, calculated and observed values for electrical and thermal conductivity of aluminium alloys. J Mater Sci 1997; 32: 4383–4390.

43.

Desai

James

. Electrical resistivity of aluminum and manganese. J Phys Chem Ref Data 1984; 13: 1131–1172.

44.

Dong

Shin

. Coupled thermomechanical multiscale modeling of alumina ceramics to predict thermally induced fractures under laser heating. J Am Ceram Soc 2015; 98: 920–928.

45.

Schelling

Phillpot

Keblinski

. Comparison of atomic-level simulation methods for computing thermal conductivity. Phys Rev B 2002; 65: 144306–144306.

46.

GoodFellow. Aluminium. http://www.goodfellow.com/E/Aluminium.html (2015, accessed 22 July 2015).

47.

Saurav

. The current understanding on the diamond machining of silicon carbide. J Phys D: Appl Phys 2014; 47: 243001–243001.

48.

Slack

. Thermal conductivity of pure and impure silicon, silicon carbide, and diamond. J Appl Phys 1964; 35: 3460–3466.

49.

Yang

Luo

Esfarjani

et al.

Thermal interface conductance between aluminum and silicon by molecular dynamics simulations. J Comput Theor Nanosci 2015; 12: 168–174.

50.

Hsieh

W-P

Lyons

Pop

et al.

Pressure tuning of the thermal conductance of weak interfaces. Phys Rev B 2011; 84: 184107–184107.

51.

Diao

Srivastava

Menon

. Molecular dynamics simulations of carbon nanotube/silicon interfacial thermal conductance. J Chem Phys 2008; 128: 164708–164708.

52.

Ong

Z-Y

Pop

. Molecular dynamics simulation of thermal boundary conductance between carbon nanotubes and SiO2. Phys Rev B 2010; 81: 155408–155408.

53.

Stevens

Zhigilei

Norris

. Effects of temperature and disorder on thermal boundary conductance at solid–solid interfaces: nonequilibrium molecular dynamics simulations. Int J Heat Mass Transf 2007; 50: 3977–3989.

54.

Abdul-Ameer

. Studying some of mechanical and thermal properties of Al-SiC composites. J Al-Nahrain Univ 2013; 16: 5–5.