Evaluation of Effective Thermal Conductivity of Graphite Flake/Aluminum Composites by Two-Dimensional Image Simulation under a Correction Function

Abstract

This study aims to evaluate the effective thermal conductivity of graphite flake/aluminum composites using two-dimensional (2D) image simulations. However, the effective thermal conductivity calculated from the two-dimensional microstructure images may not be equivalent to that measured using the experimental methods. The reason is that the two-dimensional microstructure image cannot reveal depth information based on the observation surface, which leads to the orientation difference between the graphite flakes in the 2D microstructure image and the experimental sample. Here, the orientation of the graphite flakes relative to heat flow direction was characterized by the angle between the graphite flake basal plane and the heat flow direction. The relationship between the angles in the 2D cross-sections extracted from three-dimensional (3D) models, angles in the 3D models, and aspect ratios of graphite flakes displayed in the 2D cross-sections were studied by computer simulation. We found that the angle in the 2D cross-section was larger than that in the corresponding 3D model, and the difference between the angles can result in a thermal conductivity error of up to 840 W m⁻¹ K⁻¹. In addition, all the angles and aspect ratios were distributed on a curved surface, and the curved surface function could convert the angle in the cross-section into the corresponding angle in the 3D model. Finally, the effective thermal conductivity of the graphite flake/aluminum composite with 10 vol% graphite flakes was determined using a 2D image simulation, and the interfacial thermal conductance was calculated by the reversed method.

Fig. 5 Orientation difference between the GFs in the 2D cross-sections and in the 3D models (a); TC difference between the GFs in the 2D cross-sections and in the 3D models (b).

1. Introduction

Graphite flakes (GFs)-reinforced aluminum (Al) (GFs/Al) matrix composites are lightweight, have low cost, and exhibit superior thermal properties, such as high effective thermal conductivity (ETC), and low coefficient of heat expansion (CTE). GFs/Al composites have attracted considerable attention from researchers who study high thermal conductivity (TC) materials, and GFs/Al composites have been proposed as one of the most promising candidates for thermal management materials. In previous experimental studies on the ETC of GFs/Al composites, researchers focused mainly on the orientation of GFs and the interface reaction between Al and GFs.^1–5) The interface reaction was effectively resolved through low-temperature solid sintering.¹⁾

To evaluate the ETC of GFs/Al composites, common theoretical methods include mixing rules,^1,2,6) Fricke’s equation,⁷⁾ effective medium approach (EMA),⁸⁾ and some three-dimensional (3D) theoretical models.^9,10) However, the ETC of GFs/Al composites was affected by the microstructural characteristics, such as the size, distribution, and orientation of the GFs. These theoretical methods cannot consider the effects of these microstructural characteristics on the ETC of composites. An image-based simulation proved effective in evaluating the ETC of composites from microstructure images, such as 3D image simulation, which uses ultra-high-resolution X-ray computed tomography images to reconstruct the composite microstructure.¹¹⁾ Nevertheless, 3D numerical simulations are extremely time-consuming and expensive.

Thus, in this study, two-dimensional (2D) image simulation^12,13) was employed to evaluate the ETC of GFs/Al composites. However, the orientation of the GFs in the 2D microstructure image may not be equivalent to that in the corresponding experimental sample due to a 2D microstructure image cannot provide the information in the depth direction with respect to the viewing surface. The ETC calculated using the 2D image-based simulation was not reliable. Thus, through simulation, the relationship between the orientations of GFs in 2D cross-sections extracted from 3D models, orientations of the GFs in the corresponding 3D models, and aspect ratios of GFs displayed in 2D cross-sections were investigated. Moreover, the ETC of the 10 vol% GFs/Al composite was calculated using 2D image-based simulation, and the interfacial thermal conductance at the GFs-Al interface was calculated using the reversed method.

2. Calculation Procedure

2.1 Creating 3D and 2D models

Figure 1 shows a flowchart of the calculation procedure for comparing the orientations of the GFs in the 2D model and the corresponding 3D model. First, 3D models were created, as shown in Fig. 2(a). The 3D model’s dimensions were 110 × 100 × 100 elements along the X-, Y-, and Z-axis directions, respectively. One GF was embedded in the center of the 3D model. The shape of the GFs was an elliptical cylinder, in which the semi-major axis (D) was 40 elements, semi-minor axis (d) was 20 elements, and thickness (T) was 10 elements. As shown in Fig. 2(b), the GFs in the 3D model can rotate θ_x, θ_y, and θ_z around the X-, Y-, and Z-axis, respectively, according to the right-hand rotation rule. The rotated coordinates (x, y, z) were calculated using the following rotation matrix:   

\begin{equation} \begin{bmatrix} x\\ y\\ z \end{bmatrix} = \begin{bmatrix} \textit{cos}\theta_{z} & -\textit{sin}\theta_{z} & 0\\ \textit{sin}\theta_{z} & \textit{cos}\theta_{z} & 0\\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} \textit{cos}\theta_{y} & 0 & \textit{sin}\theta_{y}\\ 0 & 1 & 0\\ -\textit{sin}\theta_{y} & 0 & \textit{cos}\theta_{y} \end{bmatrix} \begin{bmatrix} 1 & 0 & 0\\ 0 & \textit{cos}\theta_{x} & -\textit{sin}\theta_{x}\\ 0 & \textit{sin}\theta_{x} & \textit{cos}\theta_{x} \end{bmatrix} \begin{bmatrix} x_{0}\\ y_{0}\\ z_{0} \end{bmatrix}, \end{equation} (1)

where (x₀, y₀, z₀) are the initial coordinates. The x and y were set to coincide with the semi-major axis (D) and semi-minor axis (d) of the GFs, respectively. In this study, the orientation of GF in 3D models was represented by the angles α, β, and γ, which is the difference between the global coordinate X-Y-Z and the local coordinate x-y-z. In other words, α, β, and γ were the angles between the axes (X-x, Y-y, and Z-z), and their values were calculated as follows:   

\begin{equation} \alpha = \arccos \left(\frac{\overrightarrow{Ox} \cdot \overrightarrow{OX}}{|\overrightarrow{Ox|}|\overrightarrow{OX}|}\right),\quad \beta = \arccos\left(\frac{\overrightarrow{Oy} \cdot \overrightarrow{OY}}{|{\overrightarrow{Oy|}}|\overrightarrow{OY}|} \right),\quad \gamma = \arccos \left(\frac{\overrightarrow{Oz} \cdot \overrightarrow{OZ}}{|{\overrightarrow{Oz|}}|\overrightarrow{OZ}|} \right) \end{equation} (2)

where $\overrightarrow{OX}$, $ \overrightarrow{OY}$ and $\overrightarrow{OZ}$ are the unit coordinate vectors, and their values are (1,0,0), (0,1,0), and (0,0,1), respectively. $\overrightarrow{Ox} $, $\overrightarrow{Oy}$, and $\overrightarrow{Oz}$ were calculated using eq. (1).

Fig. 1

Flowchart of the calculation procedure for comparing the orientations of GFs in two- and three-dimensional models.

Fig. 2

Three-dimensional model and orientation of GFs. (a) 3D model of GFs/Al composites; (b) Orientation of GFs in 3D model; (c) Orientation of GFs in the 2D cross-sectional image; (d) Describing the angle between of GFs basal plane and heat flow; α, β, γ and θ′ are the angles between the axes, x-X, y-Y, z-Z, and z-X, respectively; α, β, and γ are within [0,90°]; OM is the line perpendicular to OY and Oz, and θ_2D is the angle between OX and OM; Rectangle O′ABC is the minimum bounding rectangle of GFs, θ_2D is the angle between X-axis and AB; XP is perpendicular to the basal plane of GFs.

Second, as many 2D cross-sections were extracted from each 3D model as possible, the cross-sections were extracted parallel to the X-O-Z plane, as shown in Fig. 2(b). The shape of GFs may be irregular in the cross-sectional images; therefore, a minimum bounding rectangle (O′ABC, see Fig. 2(c)) of the GFs was introduced to measure the orientation of the GFs. The orientation of GFs in those cross-sectional images was defined by the angle (θ_2D) between GFs and the X-axis, as shown in Fig. 2(c). θ_2D was calculated from the 2D cross-sectional images, excluding GFs smaller than 5 × 5 elements. Moreover, the aspect ratio (R) of GFs in each extracted 2D cross-sectional image was also calculated.

2.2 TC calculation of GFs

Here, the heat flow direction was set in the 3D model and in the 2D cross-sectional image to be along the X-axis. The TC of GFs along the X-axis can be calculated as follows:¹⁴⁾   

\begin{equation} \lambda_{\mathrel{/\!/}} = \lambda_{a} \left[1 - \left(1 - \frac{\lambda_{c}}{{\lambda_{a}}}\right) \textit{sin}^{2}\theta\right] \end{equation} (3)

where λ_// is the TC of GFs parallel to the X-axis direction, λ_a is the TC of GFs in the basal-plane direction (880 W m⁻¹ K⁻¹), and λ_c is the TC of GFs in the out-of-plane direction (38 W m⁻¹ K⁻¹). θ in the 3D model refers to the orientation of GFs relative to heat flow direction, which is equal to the angle between the X-axis and basal plane of the GFs (Fig. 2(d)), and it is denoted as θ_3D. θ_3D can be obtained as follows:   

\begin{equation} \theta_{3D} = |90-\theta{'}|,\quad \theta{'} = \arccos \left(\frac{\overrightarrow{Oz} \cdot \overrightarrow{OX}}{|\overrightarrow{Oz|}|\overrightarrow{OX}|}\right) \end{equation} (4)

where θ′ is the angle between X-axis and z-axis (Fig. 2(b)). In other words, θ′ is the angle between the heat flow direction and the normal vector of the basal plane, and θ_3D is the angle between the heat flow direction and basal plane.

In the 2D cross-section images, θ refers to the angle between the cross-section of the GFs and the X-axis, i.e., θ = θ_2D. θ_2D can be calculated mathematically, as follows:   

\begin{equation} \theta_{2D} = \arccos \left(\frac{\overrightarrow{OM} \cdot \overrightarrow{OX}}{|\overrightarrow{OM}||\overrightarrow{OX}|}\right) \end{equation} (5)

where OM is the line perpendicular to OX and Oz (Fig. 2(b)), and θ_2D is in the range [0, 90).

2.3 Image-based simulation method

The finite-volume method was employed to calculate the 2D temperature distributions. The temperature of the elements can be calculated using the following equation:¹³⁾   

\begin{equation} T_{x,y}^{n+1} = T_{x,y}^{n} + \frac{\Delta t}{\rho c} \left(\frac{q_{x+1,y}^{n} - q_{x-1,y}^{n}}{\Delta x} + \frac{q_{x,y+1}^{n} - q_{x,y-1}^{n}}{\Delta y} \right) \end{equation} (6)

where $T_{x,y}^{n}$ is the temperature of the element at (x, y) coordinates, $T_{x,y}^{n + 1}$ is the temperature of the element at (x, y) coordinates after a time Δt, ρ is the density, c is the specific heat, and Δx and Δy are the sizes of the elements. The heat flow qⁿ can be calculated as follows:   

\begin{align*} &q_{x+1,y}^{n} = \lambda \left(\frac{T_{x+1,y}^{n} - T_{x,y}^{n}}{\Delta x} \right), \quad q_{x-1,y}^{n} = \lambda \left(\frac{T_{x,y}^{n} - T_{x-1,y}^{n}}{\Delta x} \right) \\ &q_{x,y+1}^{n} = \lambda \left(\frac{T_{x,y+1}^{n} - T_{x,y}^{n}}{\Delta y}\right), \quad q_{x,y-1}^{n} = \lambda \left(\frac{T_{x,y}^{n} - T_{x,y - 1}^{n}}{\Delta y} \right) \end{align*}

When the heat flow move at the GFs-Al interface:   

\begin{align} &q_{x+1,y}^{n} = h (T_{x+1,y}^{n} - T_{x,y}^{n}), \quad q_{x-1,y}^{n} = h(T_{x,y}^{n} - T_{x-1,y}^{n}) \\ &q_{x,y+1}^{n} = h(T_{x,y+1}^{n} - T_{x,y}^{n}),\quad q_{x,y-1}^{n} = h(T_{x,y}^{n} - T_{x,y-1}^{n}) \end{align} (7)

where h is the interfacial thermal conductance between different materials.

Figure 3 shows the schematics of the simulation model. The model is a sandwich structure consisting of two heat sources and a composite part. The composite part was obtained from the microstructure image of the GFs/Al composites. The size of the composite part (N_x × N_y) was 450 × 600 elements, and that of the heat source was 5 × 600 elements (N_L = N_R = 5 elements). The size of each element was 1.18 × 10⁻⁶ m. Both the upper and lower sides were the periodic boundary, while the left and right sides were the adiabatic boundary. The temperatures of the left and right edge elements were fixed at 301 K and 300 K, respectively. The initial temperature of the other elements was set at 300 K, and the temperature was iteratively updated until the temperature variation was lower than 10⁻¹³ K. The temperature distribution at this moment was in a steady state, and the ETC of the GFs/Al composite, λ_eff, was calculated as follows:   

\begin{equation} \lambda_{\textit{eff}} = \frac{\lambda_{\textit{Al}} \Delta T_{12}N_{x}}{\Delta T_{LR}- N_{L}\Delta T_{12} - N_{R}\Delta T_{12}} \end{equation} (8)

where λ_Al is the TC of the Al matrix, ΔT₁₂ is the average value of the temperature difference between the first and second columns, and N_L, N_R, N_x, and N_y are the number of elements along the corresponding direction.

Fig. 3

Simulation model for the ETC calculation.

3. Results and Discussions

3.1 Orientations of GFs in 2D and 3D models

The 2D orientations of GFs, i.e., θ_2D, were obtained from eq. (5) and the 2D-cross-sectional images. The θ_2D from the 2D-cross-sectional images extracted from the same 3D model showed the same values. Figure 4 shows a comparison of θ_2D obtained from the mathematical calculation and 2D-cross-sectional images. The solid line is the θ_2D from eq. (5), and the open circles are the θ_2D from the 2D-cross-sectional images. The open circles are on the solid line (x = y), which confirms that the θ_2D values obtained from the 2D cross-sectional images were reliable.

Fig. 4

Comparison of the θ_2D obtained from mathematical calculation and 2D cross-sectional image.

The angle difference between θ_2D and θ_3D (θ_2D-3D) was calculated. Figure 5(a) shows θ_2D-3D in a standard regular triangle. The black dots represent the calculated points. The black color indicates that the value of θ_2D-3D is zero. The other colors indicate that the value of θ_2D-3D is between 0° and 90°. A few black dots were in the black region, whereas the rest were distributed in other colored regions. In the triangle, most of the θ_2D-3D in the triangle was under 20°, but there were some θ_2D-3D that exceeded 40°. Accordingly, the difference between the TCs of GFs in the 2D and 3D models (i.e., TC_3D-2D) is shown in Fig. 5(b). TC_3D-2D reached a maximum of 840 W m⁻¹ K⁻¹. Additionally, based on a comparison of Figs. 5(a) and 5(b), it was concluded that the value of TC_3D-2D may be large, even if the value of θ_2D-3D is small in the triangle.

Fig. 5

Orientation difference between the GFs in the 2D cross-sections and in the 3D models (a); TC difference between the GFs in the 2D cross-sections and in the 3D models (b).

3.2 Relationship between aspect ratio of GFs and the orientation of GFs

The aspect ratio of the GFs (R) in each 2D cross-section extracted from the same 3D model showed different values, while the corresponding θ_2D was the same. To investigate the correlation between the aspect ratios and orientations, the arithmetic average value of R was calculated and marked as $\bar{R}$, i.e., $\bar{R} = \frac{\sum R }{\sum i }$ (i is the number of R). Figure 6 shows the relationship between θ_2D-3D, $\bar{R}$ and θ_2D. According to Fig. 6(a), all data were distributed on a curved surface, which means that any of the θ_2D-3D, $\bar{R}$, and θ_2D can be expressed using the curved surface function. Subsequently, the surface fitting method showed a similar curved surface, as shown in Fig. 6(b). The fitted curved surface consists of the following functions:   

\begin{align} \theta_{\text{fitted-3D}} &= \theta_{\text{2D}} - [\text{A}_{0} + \text{A}_{1} \times \theta_{\text{2D}} \times f - A_{4} \times \theta_{\text{2D}} \times {e^{(\text{A}_{5}\times \bar{R})}}] \\ f &= \frac{1}{1+{e}^{-(\text{A}_{2} + \text{A}_{3} \times \theta_{2\text{D}} - \bar{R})}} \\ &\quad (0^{\circ} < \theta_{\text{2D}} < 90^{\circ}, 1 < \bar{R} < 6.7) \end{align} (9)

The coefficients, A₀, A₁, A₂, A₃, A₄, and A₅, are listed in Table 1, and θ_fitted-3D represents the angle between the GFs basal plane and heat flow direction in the 3D model. This function can effectively convert θ_2D into θ_3D and improve the reliability of the ETC calculated using 2D image simulation. The deviation between θ_fitted-3D and θ_3D was less than 6°. Figure 7 shows the GFs TC calculated using eq. (3). The black solid line represents the TC calculated using eq. (3), where θ is equal to θ_3D, and the gray circles represent the TC calculated using eq. (3), where θ was equal to θ_fitted-3D. The circles were tightly bound around the black solid line, and the errors were negligible.

Fig. 6

Relationship between the θ_2D, θ_2D-3D, and $\bar{R}$. Scatter image (a); Nonlinear polynomial fitting surface (b).

Table 1 Parameters in eq. (9).

Fig. 7

TC of GFs calculated using eq. (3), θ was equal to the θ_3D for the black solid line and θ_fitted-3D for the gray circle.

3.3 ETC of GFs/Al composites obtained from the 2D image-based simulation

10 vol% GFs/Al composites were fabricated via spark plasma sintering. The diameter of the GFs was 137.02 µm (98% purity), and the average particle of Al powder was 30 µm (99.9% purity). The relative density of the samples was 99.4%. The ETC of the sample, measured using a steady-state thermal-conductivity-measuring device, was 238 W m⁻¹ K⁻¹. The microstructures of the samples were observed using an optical microscope. Figure 8(a) shows the microstructures of the 10 vol% GFs/Al composites. To evaluate the TC distributions of GFs in the microstructures, different colors were used to mark the GFs according to the orientations of the GFs, as shown in Fig. 8(b). For example, the red color indicates that the orientation of the GFs was parallel to the heat flow direction, and the GFs had high TC; the other colors indicate that the GFs tilt from the heat flow direction, decreasing TC. To calculate the ETC of 10 vol% GFs/Al composites using 2D image-based simulation, the orientations of GFs in the microstructures were corrected using eq. (9), and the TC of GFs was calculated using eq. (3). The TC of the Al matrix (λ_Al-matrix) was calculated using the following equation:¹⁵⁾   

\begin{align} \lambda_{\text{Al-matrix}} &= \frac{1}{4}[\lambda_{p} (3v_{p} - 1) + \lambda_{Al} (3v_{\textit{Al}} - 1) \\ &\quad + ([\lambda_{p}(3v_{p}-1) + \lambda_{\textit{Al}} (3v_{\textit{Al}} - 1)]^{2} + 8\lambda_{p}\lambda_{\textit{Al}})^{\frac{1}{2}}] \end{align} (10)

where the TCs of the raw material Al powder (λ_Al) was taken as 228 W m⁻¹ K⁻¹ as in our previous works, TC of the pores (λ_p) was 0.214 W m⁻¹ K⁻¹, volume fraction of Al (v_Al) was 0.994, and volume fraction of pores (v_p) was 0.006. The pores in this study were assumed to be embedded in the Al matrix. The λ_Al-matrix was 226 W m⁻¹ K⁻¹. Figure 9 shows the calculated ETC considering the interfacial heat conduction (h) from 1.1 × 10³ to 1.1 × 10⁹ W m⁻² K⁻¹. When the calculated ETC was equal to the experimental ETC, the value of h was 1.8 × 10⁷ W m⁻² K⁻¹. This result was not equal to the h given by Xue, Zhou and Schmidt (4.5 ∼ 5 × 10⁷ W m⁻² K⁻¹).^2,8,16) To view the Al-GFs interface (see Fig. 10), pores can be observed. Additionally, Li,⁴⁾ Zhou,⁵⁾ and Huang¹⁷⁾ reported that an amorphous layer interface could be readily formed at the Al-GFs interface. Therefore, the small difference between h calculated in this work and h reported by Xue, Zhou, and Schmidt can be attributed to the pores and amorphous layer at the interface.

Fig. 8

Microstructure images of the 10 vol% GFs/Al composite; (a) The OM images; (b) Distribution of the GFs TC.

Fig. 9

ETC of the 10 vol% GFs/Al composite considering the interfacial thermal conductance (h) from 1.1 × 10³ to 1.1 × 10⁹ W m⁻² K⁻¹.

Fig. 10

A SEM image of the microstructure at the Al-GFs interface.

4. Conclusions

This study aimed to calculate the ETC of GFs/Al composites using 2D image-based simulations. However, the ETC calculated from the microstructure images might not be equivalent to that measured using experimental methods. Because the 2D microstructure image cannot reveal the depth information related to the observation surface, the orientation of GFs in the 2D microstructure image differs from that of the experimental sample. We studied the relationships between the orientation of GFs in 2D cross-sectional images, orientation of GFs in 3D models, and aspect ratio of GFs in the 2D cross-sectional images. Based on the results, the following conclusions were drawn:

The orientation of the GFs was defined by the angle between the GFs basal plane and X-axis. The angle was marked as θ_2D in the 2D cross-section and θ_3D in the 3D model.

1. (1)    θ_2D was generally larger than θ_3D, and the difference between them could cause the TC error to be as high as 840 W m⁻¹ K⁻¹. No regularity was observed between the θ_2D and θ_3D differences.
2. (2)    The data for θ_2D, θ_3D and aspect ratio were on a curved surface, and the curved surface function can be used to convert θ_2D into θ_3D. The difference between the actual θ_3D and θ_3D values calculated using the curved surface function was less than 6°.
3. (3)    The ETC of the 10 vol% GFs/Al composite was evaluated using 2D image-based simulation, where the suggested converting θ_2D into θ_3D was employed.
4. (4)    The interfacial thermal conductance at the GFs-Al interface was calculated as 1.8 × 10⁷ W m⁻² K⁻¹ which has the same order of magnitude as reported in other papers (4.5 ∼ 5 × 10⁷ W m⁻² K⁻¹).

REFERENCES

• 1)   J.K.  Chen and  I.S.  Huang: Compos., Part B 44 (2013) 698–703. doi:10.1016/j.compositesb.2012.01.083
• 2)   C.  Xue,  H.  Bai,  P.F.  Tao,  J.W.  Wang,  N.  Jiang and  S.L.  Wang: Mater. Des. 108 (2016) 250–258. doi:10.1016/j.matdes.2016.06.122
• 3)   T.  Etter,  P.  Schulz,  M.  Weber,  J.  Metz,  M.  Wimmler,  J.F.  Löffler and  P.J.  Uggowitzer: Mater. Sci. Eng. A 448 (2007) 1–6. doi:10.1016/j.msea.2006.11.088
• 4)   W.  Li,  Y.  Liu and  G.  Wu: Carbon 95 (2015) 545–551. doi:10.1016/j.carbon.2015.08.063
• 5)   C.  Zhou,  G.  Ji,  Z.  Chen,  M.  Wang,  A.  Addad,  D.  Schryvers and  H.  Wang: Mater. Des. 63 (2014) 719–728. doi:10.1016/j.matdes.2014.07.009
• 6)   C.  Xue,  H.  Bai,  P.F.  Tao,  N.  Jiang and  S.L.  Wang: J. Mater. Eng. Perform. 26 (2017) 327–334. doi:10.1007/s11665-016-2447-z
• 7)   H.  Kurita,  T.  Miyazaki,  A.  Kawasaki,  Y.  Lu and  J.-F.  Silvain: Compos., Part A 73 (2015) 125–131. doi:10.1016/j.compositesa.2015.03.013
• 8)   C.  Zhou,  W.  Huang,  Z.  Chen,  G.  Ji,  M.L.  Wang,  D.  Chen and  H.W.  Wang: Compos., Part B 70 (2015) 256–262. doi:10.1016/j.compositesb.2014.11.018
• 9)   F.  Gori and  S.  Corasaniti: Int. J. Heat Mass Transfer 77 (2014) 653–661. doi:10.1016/j.ijheatmasstransfer.2014.05.047
• 10)   T.  Wejrzanowski,  M.  Grybczuk,  M.  Chmielewski,  K.  Pietrzak,  K.J.  Kurzydlowski and  A.  Strojny-Nedza: Mater. Des. 99 (2016) 163–173. doi:10.1016/j.matdes.2016.03.069
• 11)   M.W.  Czabaj,  M.L.  Riccio and  W.W.  Whitacre: Compos. Sci. Technol. 105 (2014) 174–182. doi:10.1016/j.compscitech.2014.10.017
• 12)  K. Sugio, K. Kono, Y.B. Choi and G. Sasaki: Materials Science Forum, (Trans Tech Publ, 2018) pp. 1939–1943.
• 13)   K.  Sugio,  Y.-B.  Choi and  G.  Sasaki: Mater. Trans. 57 (2016) 582–589. doi:10.2320/matertrans.MC201505
• 14)  B.T. Kelly: Physics of Graphite, (Applied Science Publishers, London, 1981).
• 15)   R.  Landauer: J. Appl. Phys. 23 (1952) 779–784. doi:10.1063/1.1702301
• 16)   A.J.  Schmidt,  K.C.  Collins,  A.J.  Minnich and  G.  Chen: J. Appl. Phys. 107 (2010) 104907. doi:10.1063/1.3428464
• 17)   Y.  Huang,  Q.  Ouyang,  Q.  Guo,  X.  Guo,  G.  Zhang and  D.  Zhang: Mater. Des. 90 (2016) 508–515. doi:10.1016/j.matdes.2015.10.146