- Email: [email protected]
Effect of diamond shapes and associated thermal boundary resistance on thermal conductivity of diamond-based composites
Available online at www.sciencedirect.com
Computational Materials Science 41 (2007) 156–163 www.elsevier.com/locate/commatsci
Effect of diamond shapes and associated thermal boundary resistance on thermal conductivity of diamond-based composites J. Flaquer
a,b,*
, A. Rı´os b, A. Martı´n-Meizoso
a,b
, S. Nogales c, H. Bo¨hm
c
a
c
TECNUN, University of Navarra, Paseo de Manuel Lardiza´bal, 15, 20018 San Sebastia´n, Spain b CEIT, Paseo de Manuel Lardiza´bal, 13, 20018 San Sebastia´n, Spain Technical University of Vienna, Institute of Lightweight Design and Structural Biomechanics, Gusshausstrasse 27-29, A-1040 Vienna, Austria Received 17 January 2007; received in revised form 14 March 2007; accepted 15 March 2007 Available online 31 May 2007
Abstract This work presents techniques to generate densely packed diamond-like particles, for example, as present in a metal matrix composite of Al, Cu, or Ag reinforced with single crystals of synthetic diamonds; as intended to be used for heat sinks or very efficient thermal conductors, spreaders, etc. Two different ways are described to predict the thermal conductivity of the composite: line tracing and surface-averaged Hasselman and Johnson. The effects of diamond volume fraction, size and, in particular, shapes on the composite thermal conductivity are discussed. 2007 Elsevier B.V. All rights reserved. Keywords: Packing of polyhedra; Metal matrix composites; Particle-reinforced composites; Diamonds; Thermal conductivity; Simulated annealing
1. Introduction Materials with a huge thermal conductivity are required for many technical applications: heat sinks, rocket nozzles, heat spreaders for electric and microelectronic applications. . . Diamond is the isotropic material with the highest thermal conductivity, but diamonds are extremely expensive if a very large size is required, it is indeed hard to machine, etc. On the other hand, silver, copper and aluminium are the metals with the largest thermal conductivities. Within the ExtreMat [1] project a large effort has been made to develop new materials with improved thermal conductivities for different temperature ranges (room temperature and high temperature: above 500 C). One of the most promising materials is copper, silver or aluminium reinforced with diamonds. Synthetic monocrystalline diamonds, with particle diameter of up to about * Corresponding author. Address: TECNUN, University of Navarra, Paseo de Manuel Lardiza´bal, 15, 20018 San Sebastia´n, Spain. Tel.: +34 943219877; fax: +34 943213076. E-mail address: jfl[email protected] (J. Flaquer).
0927-0256/$ - see front matter 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2007.03.016
600 lm are produced (at very large temperature and pressure, by hot isostatic pressing) and are actually twice cheaper per kilogram than silver (see, for example, SDG121 to 126 grades from Indus-Global [2], with a grit size from 30 to 60; or MBD-4 grade diamonds, with a mesh 70/80 from Henan Famous Industrial Diamond Co., Ltd. [3]). If an aluminium matrix is used, and the diamonds surface is not specifically prepared, the aluminium matrix bonds to square {0 0 1} surfaces of the diamonds, but does not stick to the hexagonal diamond faces (with {1 1 1} orientations) [4,5]. Therefore, the thermal resistance at both interfaces between the metal matrix and either the square or the hexagonal diamond faces are expected to be very different. There are models to estimate the thermal conductivity of composite materials reinforced with particles. For example, Maxwell’s approximation [6], but only for dilute spherical particles, which does not match with a high volume fraction of cubo-octahedron particles and, even worse, a perfect thermal connection is assumed between the matrix and particles. Maxwell’s approximation is, in our case, an
J. Flaquer et al. / Computational Materials Science 41 (2007) 156–163
157
upper bound to the actual thermal conductivity [7]. A more accurate model is that of Hasselman and Johnson [8], where an interface resistance is introduced, but spherical particles are still assumed. 2. Model Line tracing is used to compute the composite thermal conductivity in our model. In this case, diamonds are generated as truncated octahedrons and truncated cubes; see Figs. 1 and 2, respectively. 2.1. Particle generation A first diamond is located at the origin and oriented along the axis. Afterwards different diamonds are generated at random (not colliding with the previous ones) and with a random spatial orientation. The generation of random orientations, uniformly distributed in the space, can be found elsewhere [9,10] hi ¼ cos1 ð1 randÞ; wi ¼ 2p rand; p /i ¼ rand 2
Fig. 2. Enumeration of vertices, edges and faces for a truncated cube. The numbers without background are for vertices; bold for those at the front and regular for those at the back. The numbers on a square background stand for the faces; black on the front and pale grey for those at the rear. Edge directions are indicated by arrows.
ð1Þ z
where rand stands for random numbers uniformly distributed in the range 0–1 (most computer programming languages provide subroutines for the generation of such pseudo-random numbers). To specify a crystal orientation, three (Euler) angles should be provided, as shown in Fig. 3 and Eq. (1). Because of the cubic symmetry of the consid-
θ
n3
n1
y
n2
ϕ ψ
x Fig. 3. Schematic of the axis and (Euler) angles used to specify a crystal orientation. n1, n2 and n3 are normal to {0 0 1} faces.
ered polyhedra, wi can be reduced to only the one-fourth of the range p ð2Þ wi ¼ rand: 2 The size of each truncated octahedron, Ri, is generated with a normal distribution about a mean size, Rmean, with a given standard deviation, rR: Ri ¼ randn rR þ Rmean ; Fig. 1. Enumeration of vertices, edges and faces for a truncated octahedron. The numbers without background are for vertices; bold for those at the front and regular for those at the back. The numbers on a square background stand for the faces; black on the front and pale grey for those at the rear. Edge directions are indicated by arrows.
ð3Þ
where randn stands for a random number generated from a typified normal distribution (with a mean equal to 0, and a standard deviation of 1). Most computer languages provide for subroutines for such a purpose. Anyway, computer codes for the generation of quasi-random numbers with
158
J. Flaquer et al. / Computational Materials Science 41 (2007) 156–163
normal deviates can be found in [11]. This size refers to the radius of the sphere that circumscribes the original octahedron (it will be now be truncated). Different diamond-looking shapes are generated by truncation of the octahedra. The parameter a decides the height of the cut in comparison with Ri. a = 0 results in the original octahedron; the Kelvin Cell (with six little squares and eight large hexagons) is obtained for a = 1/3; a second tetrakaidecahedron (cuboctahedron) is obtained, formed by six squares and eight equilateral triangles, for a = 1/2; and cubes are obtained for a P 2/3. Note that according to Kelvin’s conjecture [12] his tetrakaidecahedron fills (tessellates) the space completely (a solution with slightly smaller surface was found by Weaire and Phelan [13] in 1993). But these identical tetrakaidecahedra should be very carefully located to fill completely the volume. One single Kelvin’s tetrakaidecahedron only fills one half of the smallest in which it is inscribed. pffifficube ffi The side of that cube is 2 2 l, where l is the length of any of the edges pffiffiffi of 3the tetrakaidecahedron; the cube volume is then pffiffiffi16 3 2 l and the volume of the tetrakaidecahedron is 8 2 l . So, it will be shown how difficult is to get large volume fractions of such particles when they are located at random. The computer code generates shapes at random, with a normal distribution about any given mean shape, amean, with any value for its standard deviation, ra: ai ¼ randn ra þ amean :
ð4Þ
This covers most of actual diamond shapes. In particular, all these polyhedra belong to the uniform polyhedra category: they have all their vertices identical, and – as a result of that – all of them located on a second sphere of smaller radius qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ri ¼ Ri 2a2i 2ai þ 1 () 0 6 ai 6 1=2; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð5Þ Ri ¼ Ri 6a2i 8ai þ 3 () 1=2 6 ai 6 2=3; p ffiffi ffi Ri ¼ Ri ð1 ai Þ 3 () 2=3 P ai P 1 as shown in Fig. 4. That feature helps to decide about the collision possibility in a very efficient way. 2.2. Collision detection Collision detection among polyhedra is done in four steps. 2.2.1. Spheres interference The first and fast check is to compare distance between every possible diamond couple, dij, and to compare with the sum of their radii: Ri is the radius of the sphere which circumscribes the polyhedron, ith (truncated octahedron). If the circumscribed spheres do not touch each other, there is no possibility of collision and Eq. (6) is satisfied: Ri þ RJ < d ij :
ð6Þ
Fig. 4. Maximum diameter of the polyhedron after truncation of a regular octahedron. Those polyhedra with regular polygonal faces are indicated with arrows.
If Eq. (6) is not satisfied, there is a possible collision and the analysis should go much deeper. 2.2.2. Vertex within another polyhedron The second check is whether a vertex of polyhedron i is within the volume of the polyhedron j. If all the faces are well oriented, it is to say: their normals point outwards from the polyhedron, the vertex is within the polyhedron when it is to the left of every face, which means that when its coordinates are substituted into the equation of the plane containing a face, the result is always negative. If no vertex of polyhedron i is within the polyhedron j, next we should check the other way: Does a vertex of j lie within polyhedron i? If no vertex is found within the other polyhedron we might feel that there is no collision, and this is the case in a 2D world (if we were working with polygons), but two polyhedra may also collide with each other when one edge, belonging to one of them, cuts through any face of the second one. Therefore, a collision is only rejected when additionally no edge is found to cut any face. 2.2.3. One shortcut: split plane Of course the computation of edge/face collision checks consumes a lot of time (most of the time) and any trick to save this verification deserves the effort. The intermediate step between the vertices detection and edge to face computation is to find a split plane between both polyhedra. If it is possible to find a plane that keeps all the vertices of the first polyhedron to one side and all of the second polyhedron to the other, they certainly do not collide. If we work with convex polyhedra (as it is the case) the candidates for
J. Flaquer et al. / Computational Materials Science 41 (2007) 156–163
split planes may be reduced to only their faces. Therefore, this intermediate check is also fast and very efficient, saving a lot of computational time.
E¼
n1 X n X
159
maxðjxi xj j; jy i y j j;jzi zj jÞ þ maxðLx ;Ly ;Lz Þ;
j¼1 i¼jþ1
ð7Þ 2.2.4. Edge/face intersection Every edge (line) cuts a face (plane) at one point and it should be checked if this intersection point is within the polygon (actual face) or outside. A simple code in Matlab [14] is used to generate the computer programme. For every edge of the first polyhedron (i) the intersection with the plane defined by every face of polyhedron j is calculated. Now it is decided if the intersection point is located within the actual edge (between their vertices) or not; and afterwards if the intersection point is within the face or not. From the intersection point, two vectors are built towards two consecutive vertices of that face. Vertices are chosen in a clockwise sequence. The cross product of both vectors is positive if the intersection point is to the right of this face side, and negative otherwise. The point is only inside the face if all the cross products, from all the consecutives couples of vectors to all the edges of this face, are negative. 2.2.5. Note In theory, it is not necessary to check for the position of the vertices inside the other polyhedron, but only to check for the edge of one polyhedron intersecting a face of the second one. In fact it is not required to check the other way: the edges of the second intersecting the faces of the first. That edge/face check does all the required collision verifications. Problems arise when the intersection point is too close to some edge or vertex: the cross product tends to zero. So, the second step is to compute the distance of the intersection point to the vertices of the second polyhedron and its edges. Thus, the above described checking for vertices is most useful to avoid singular cases when doing the edge/face comparisons, for example, an edge lying on a face. (The closer the particle packing, the more frequently these border line situations are encountered.) So it is a good idea to clean up the edge/surface collision detection as much as possible of the difficult cases by means of those previous checks. 2.3. Packing As it is well known, a generalized movement is formed by a spatial movement (3 degrees of freedom) plus three additional twists in the space. Particle densification is performed by simulated annealing, a well-known technique [15] that will be explained in the next section. The main idea is to achieve a diamond distribution that minimizes the objective function, E, given in Eq. (7), which can be considered as the ‘‘energy’’ of the whole system.
where xi, yi and zi stand for the coordinates of the centre of the ith polyhedron; Lx, Ly and Lz represent the lengths of the minimum size box, containing all the particle’s centres. A generalized random movement is generated for a randomly selected particle; each moduli being normally distributed. The movement is always accepted when the value for the expression (7) decreases. The minimization of Eq. (7) clusters the initial distribution of diamonds into a cubic shape. Of course, other energy functions can be used, for example E¼
n X n 1X d3 2 j¼1 i¼1 ij
ð8Þ
or E¼
p 6
n X
ðx2i þ y 2i þ z2i Þ3=2 ;
i¼1
where dij represents the distance between the centres of particles i and j. In the second case, the particle distance to the absolute centre is penalized according to the volume of the sphere of this size. But both functions tend to form a sphere of particles and, if a cubic box is pursued, the corners will be empty spaces and the volume fraction not as large as with previous Eq. (7). Note that none of the energies proposed in Eqs. (7) or (8) changes when a particle is rotated, maintaining its centre; so, twists do not affect the system energy (but they may help to pack particles in a later displacement; so they are permitted and so are all movement without energy change). 2.3.1. Simulated thermal annealing Simulated annealing uses, in our case, generalized movements. The optimization procedure tries to replicate the minimization of the target function, as the energy of a system that cools down slowly and follows Boltzmann’s law. Initially the system is at a ‘‘temperature’’ T0 and the system’s total energy is E0. We, arbitrarily, take an initial value for T0, for example, T0 = 1001 K and the system ‘‘energy’’ is computed according to Eq. (7), as explained before. The algorithm proceeds as follows. A temperature drop DT is defined, which, also arbitrarily will be taken as 1 K. The actual ‘‘temperature’’ is set to its initial value: T = T0. A particle is chosen at random. A generalized displacement (displacements plus rotations: Dx, Dy, Dz, D/, Dh, Du) is chosen, also at random. This particle movement is checked for collision with the other particles. If the movement is not physically possible (there is a collision with another particle), the particle is returned to its original position and another particle is chosen for movement. If the proposed movement is physically possible (there is no collision) we
160
J. Flaquer et al. / Computational Materials Science 41 (2007) 156–163
Kelvin’s truncated octahedron (tetrakaidecahedron). Other size distributions, bimodal distributions and shape distribution are possible, but the results are being computed and the results will be shown in a future paper. 3.1. Particle volume fraction
Fig. 5. Cluster of 26 diamond-like polyhedra with different sizes and shapes.
shall compute the ‘‘energy’’ modification produced within the system, DE. If the energy of the system is reduced, DE < 0, the change is always accepted; the temperature is reduced by a further step: T = T DT; and another particle is chosen for movement. If the system energy is not reduced, DE P 0, a random probability (uniformly distributed between 0 and 1), p, is generated and compared with kT 0 DE q ¼ exp ; ð9Þ T
Volume fraction is computed within a parallelepiped or box, parallel to the coordinate axes, as shown in Fig. 6. This box contains all the particle centres and it is the smallest one. Its maximum dimension appears in the expression that calculates the ‘‘energy’’ of the system. The volume fraction is obtained by dividing the volume of the diamonds contained within this box by the volume of the box. Observe that for a small system, like the one shown in Fig. 6, most of diamonds intersect the counting box. Among 27 particles, there is usually only one – at the centre of the box – not intersecting the counting box, and the other 26 intersect the box. Even for a much larger system, for example 125 particles, only 27 will not cut the counting box; that is only a 22% of the particles are completely embedded in the box. 3.2. Thermal conductivity Thermal conductivity is computed, for a given diamond packing by line tracing. Diamonds and the aluminium metal matrix have different thermal conductivities: kd = 1800 W/m K [5,16], and km = 245 W/m K [17]. If a diamond volume fraction Vd = 0.64 is introduced, Maxwell’s equation [6]
where k is an arbitrary constant (we usually use k = 10). If p < q the new situation is accepted – even when there is no reduction in the energy; the temperature is updated to T = T DT, and another particle is chosen for movement. The algorithm ends when, after a significant consecutive number of trials, nothing changes (no movement is accepted). In deciding a figure for this number of trials it helps to observe the evolution of the system energy that usually, after this number of trials, flattens down into a horizontal line. The smallest temperature is kept to be 1 K and never drops to zero. Also, to avoid too many changes with energy increments, at the beginning of the process, the value of q is limited to 0.3. Thermal agitation allows the system to escape from local minima (traps) and, if appropriately performed, to get actual system minima. It is shown to allow for the densification of the particles in a very efficient way and to obtaining, in this way, very large particle volume fractions, see Fig. 5. 3. Results In the example a uniform size distribution is used for the diamonds (R = 100 lm), and their shape is a perfect
Fig. 6. Cluster of 27 diamond-like polyhedra with identical sizes and shapes. Volume fraction of diamonds is 0.64 within the indicated parallelepiped.
J. Flaquer et al. / Computational Materials Science 41 (2007) 156–163
2 kkmd 1 V d þ kkmd þ 2 kc ¼ km 1 kkmd V d þ kkmd þ 2
ð10Þ
predicts a thermal conductivity for the composite of kc = 810 W/m K, which is clearly an upper bound for the actual case. Maxwell assumes perfectly bonded spherical particles diluted in the matrix. According to Torquato [7], Maxwell’s approximation gives good estimates for non-dilute conditions (provided that the spheres are well separated from each other), and also for different types of spheres. Additionally two interfacial thermal conductances should be assumed for the contacts between the {0 0 1} diamond faces (very well bonded to the Al matrix) and the contacts between the {1 1 1} diamond faces (with a very poor bond) and the metal matrix. The following values are used: h001 = 1 · 108 W/m2 K and h111 = 1 · 107 W/ m2 K, so, about 10 times worse for the interfacially poorly bonded {1 1 1} diamond surfaces. These values have been estimated from typical values found in the literature for other material systems [18,19]. The model of Hasselman and Johnson [8] introduces the effect of a thermal barrier resistance at the interface between particles and matrix 2 kkmd 1 krhd V d þ kkmd þ 2 þ 2 krhd kc ¼ km ; ð11Þ 1 kkmd þ krhd V d þ kkmd þ 2 þ 2 krhd where r is the radius of the particles, which are assumed to be spherical and of identical size, and h the interface conductance. Of course, Eq. (11) was not derived for cubooctahedral-shaped particles, nor for different thermal barrier resistances at different crystal surface orientations. Actual measurements of the composite conductivity, by Ruch et al. [5] for the same volume fraction of diamonds (Vd = 0.62) provide values ranging from kc = 130 to 670 W/m K, depending on the production route: mechanically assisted infiltration (squeeze casting) or gas pressure infiltration, respectively. The reason for such large difference is that in the squeeze cast samples, no chemical bonding (no interfacial carbides) is observed at all, neither on {1 1 1} nor {1 0 0} faces, whereas there is some chemical bonding between aluminium and {1 0 0} diamond faces in the gas pressure infiltrated sample. Line tracing draws lines, parallel to the preferred axis, and computes the different lengths on the metal matrix and within diamonds, Lm and Ld, respectively. In addition, it counts how many square faces, N4, are cut and how many hexagonal, N6, are intersected by the line. The thermal resistance (in units of K/W) is estimated for this line as the serial addition: 1 Lm Ld N 4 N6 Ri ¼ þ þ þ ; ð12Þ si km kd h001 h111 where si represents the cross-sectional area assigned to ith line.
161
The volume fraction of diamonds within line i can be expressed as V d;i ¼
Ld Ld ¼ ; Lm þ Ld Li
ð13Þ
where Li = Lm + Ld is the total length of the ith line. In analogy the fraction of {0 0 1} faces over the total particle surface, S001, in line i can be written as S 001;i ¼
N4 N4 ¼ : N4 þ N6 Ni
ð14Þ
Accordingly, Eq. (12) can be modified as 1 1 V d;i V d;i S 001;i 1 S 001;i Ri ¼ þ þ Li þ Li ni ; si km kd h001 h111 ð15Þ where ni = Ni/li is the density of intersected interfaces per unit length. The thermal resistance of the whole composite is estimated as the parallel resistance for all the box, R, lines, in the present case a regular set of 20 · 20 = 400 lines. X 1 1 P ¼ : ð16Þ Ri s i R si i
i
A lower bound for the thermal conductivity of the composite, kc is then estimated as 400 1 1 X 1 ¼ Li R Li si R i si i¼1 1 400 X 1 V d;i V d;i S 001 1 S 001 ¼ þ þ ; þ ni km kd h001 h111 i¼1
kc ¼
ð17Þ where the Li are assumed to be equal for all lines. Eq. (17) can be evaluated with relatively low computational effort for given phase arrangements. For the one shown in Fig. 6, the obtained composite thermal conductivity is kc = 255 W/m K. The percentage of {0 0 1} faces over the total particle surface, S001 = N4/(N4 + N6), depends on the truncation ratio 3a2 pffiffiffi () 0 6 a 6 1=2; 3a2 þ 3ð1 3a2 Þ pffiffiffi 2 4 3ð2 3aÞ ¼ pffiffiffi 4 3ð2 3aÞ2 þ 6½4ð1 aÞ2 2ð2 3aÞ2
S 001 ¼ S 001
ð18Þ
() 1=2 6 a 6 2=3; S 001 ¼ 1 () 2=3 P a P 1; as it is shown in Fig. 7. Alternatively, a mean interface thermal conductance can be estimated as h ¼ h001 S 001 þ h111 ð1 S 001 Þ:
ð19Þ
This value, which depends on a, can be introduced into Hasselman and Johnson’s Eq. (11) to estimate the effect of the diamond shapes. The results are shown in Fig. 8.
162
J. Flaquer et al. / Computational Materials Science 41 (2007) 156–163
Fig. 7. Percentage of {0 0 1} faces against the truncation ratio of the octahedron.
Observe that line tracing ‘‘extrudes’’ the composite material in spaghetti-like lines and then joins their ends, with each ‘‘column’’ being perfectly isolated with respect to its neighbours. For sufficiently long lines through statistically homogeneous composites the first term on the righthand side of Eq. (17) can be recognized as the classical Wiener [20] lower bound for the thermal conductivity of composites with perfect interfaces. The second term accounts for the interfacial resistances in a consistent way. The Wiener bounds are one-point bounds [7] that are rather slack for composites with ideal interfaces, and Eq. (17) is very conservative for particle reinforced composites with interfacial resistances because heat flow around obstacles is not accounted for. This obviously can lead to very low predictions in the case of a conductive matrix containing particles that are weakly conducting or have interfaces of low conductance. In the limit of high defect densities, ni, or low conductances h001 and/or h111, the predictions of Eq. (17) are dominated by the second term on the right-hand side. Accordingly, there is a marked effect of the particle size. The predicted overall conductivity tends to vanish for very small particles with finite interfacial conductances despite the presence of a contiguous conducting matrix, i.e., the trivial lower bounds are attained. The Hasselman–Johnson model gives identical results to Mori-Tanaka estimates [21] and can, accordingly be viewed as being related to two-point (Hashin–Shtrikman) bounds. For radii exceeding a critical radius of rc ¼
Fig. 8. Effect of the diamond sizes and shapes (truncation ratio) on the thermal conductivity of the diamond/Al composite, for a constant diamond volume fraction Vd = 0.64 with a thermal conductivity, kd = 1800 W/m K. The short-dashed horizontal line represents the thermal conductivity of the aluminium matrix, km = 245 W/m K. The black star stands for the line tracing computation with interface conductances of h001 = 1 · 108 W/m2 K and h111 = 1 · 107 W/m2 K, a particle radius of r = 50 lm and truncation ratio a = 1/3. The solid and long-dashed lines represent predictions of the Hasselman–Johnson model for r = 50, 100, and 200 lm, respectively; evaluated with averaged conductances following Eqs. (18) and (19). The blue star pertains to the same conditions as the line tracing model. The vertical line indicates the range of experimental measurements [5]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
4. Discussion Different procedures are proposed to compute or bound the thermal conductivity of metal matrix composites. These procedures have different pros and cons, as follows.
1 kd km h kd km
ð20Þ
[22], its results are lower estimates and for r < rc upper estimates. When combined with Eq. (20) it can provide estimates at low computational cost, but cannot account for details of the actual particle shape and uses a coarse approximation for the effects of the distribution of the interfacial conductances on the particle surfaces. Details of the microgeometry may obviously be studied by combining the finite element method with an appropriate micromechanical approach [23]. Among these, periodic homogenization with unit cells requires periodic particle arrangements, making windowing methods and embedding models more appropriate choices for geometries as shown in Figs. 5 and 6. However, as the particle volume fraction increases, the average distance between particles decreases, which requires meshes with high number of small elements. This leads to models with high numbers of degrees of freedom that pose major demands in computer capacity. 5. Conclusions A fast and efficient way of generating microstructures composed of monocrystalline diamonds was described. A set of cubo-octahedral geometries has been used. The geometrical advantages of using uniform polyhedra have been exploited, and the way of computing possible collisions has been detailed. It may be used to model different properties
J. Flaquer et al. / Computational Materials Science 41 (2007) 156–163
of MMCs reinforced with single crystal synthetic diamonds, as used within ExtreMat project for improved thermal conductivity materials. A shape selection of diamonds, increasing the amount of {0 0 1} faces, improves the composite thermal conductivity from kc = 525 W/m K for Kelvin tetrakaidecahedra to 700 W/m K for cubical diamonds. A 33% improvement in the composite thermal conductivity may be expected in this way. The same improvement is obtained using larger diamonds (and reducing the interface area per unit of volume) to r = 180 lm, as follows from Fig. 8. Acknowledgements This work has been carried out within the framework of the European Integrated Project ‘‘ExtreMat’’ (New Materials for Extreme Environments). Thanks are given to the European Commission for funding this Integrated Project under FP6 (Contract No. NMP3-CT-2004-500253), as well as by the Spanish Ministry for Science and Technology (FTN2003-08228-C03-02). It only reflects the view of the authors and nor the European Community nor the Spanish Ministry for Science and Technology are liable for any use of the information contained therein. References [1] ExtreMat, Integrated European Project: New Materials for Extreme Environments, funded under FP6 (NMP3-CT-2004-500253), 2004. [2] Indus Global Superabrasives, India,.
163
[4] S. Kleiner, F.A. Khalid, P.W. Ruch, S. Meier, O. Beffort, Scripta Mater. 55 (2006) 291–294. [5] P.W. Ruch, O. Beffort, S. Kleiner, L. Weber, P.J. Uggowitzer, Compos. Sci. Technol. 66 (2006) 2677–2685. [6] J.C. Maxwell, A Treatise on Electricity and Magnetism, 3rd ed., Oxford University Press, Oxford, 1904. [7] S. Torquato, Random Heterogeneous Materials Microstructure and Macroscopic Properties, Springer, New York, 2002. [8] D.P.H. Hasselman, L.F. Johnson, J. Compos. Mater. 21 (1987) 508– 515. [9] R. Rodrı´guez-Martı´n, I. Ocan˜a, A. Martı´n-Meizoso, Mater. Sci. Forum 482 (2005) 171–174. [10] D. Zwillinger, Standard Mathematical Tables and Formulae, 30th ed., CRC Press, Boca Rato´n, FL, 1996. [11] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes The Art of Scientific Computing, Cambridge University Press, Cambridge, MA, 1986. [12] W. Thomson, Acta Math. 11 (1887) 121–134. [13] D. Weaire, R. Phelan, Phil. Mag. Lett. 69 (1994) 107–110. [14] Matlab 7.0, The MathWorks Inc., 1984. [15] P.S.M. van Laarhoven, E.H.L. Aarts, Simulated Annealing Theory and Applications, Kluwer Academic Publishers, Dordrecht, 1992. [16] J.E. Field, The Properties of Natural and Synthetic Diamond, 3rd ed., Academic Press, London, 1992. [17] CES Selector Version 4.5, Cambridge Engineering Selector, Granta Design Ltd., Cambridge, UK, 2004. [18] R.J. Stoner, H.J. Maris, Phys. Rev. B 48 (1993) 167373. [19] H.-K. Lyeo, D.G. Cahill, Phys. Rev. Lett. 96 (2006) 186101. [20] O. Wiener, Abh. Math.-Natur. Sa¨chs. Akad. Wiss. Leipzig 32 (1912) 509. [21] Y. Benveniste, J. Appl. Phys. 61 (1987) 2840. [22] D. Duschlbauer, Computational Simulation of the Thermal Conductivity of MMCs under Consideration of the Inclusion-Matrix Interface, VDI-Verlag, Du¨sseldorf, 2004. [23] H.J. Bo¨hm, A short introduction to continuum micromechanics, in: H.J. Bo¨hm (Ed.), Mechanics of Microstructured Materials, vol. 1, Springer, Vienna, 2004.
Computational Materials Science 41 (2007) 156–163 www.elsevier.com/locate/commatsci
Effect of diamond shapes and associated thermal boundary resistance on thermal conductivity of diamond-based composites J. Flaquer
a,b,*
, A. Rı´os b, A. Martı´n-Meizoso
a,b
, S. Nogales c, H. Bo¨hm
c
a
c
TECNUN, University of Navarra, Paseo de Manuel Lardiza´bal, 15, 20018 San Sebastia´n, Spain b CEIT, Paseo de Manuel Lardiza´bal, 13, 20018 San Sebastia´n, Spain Technical University of Vienna, Institute of Lightweight Design and Structural Biomechanics, Gusshausstrasse 27-29, A-1040 Vienna, Austria Received 17 January 2007; received in revised form 14 March 2007; accepted 15 March 2007 Available online 31 May 2007
Abstract This work presents techniques to generate densely packed diamond-like particles, for example, as present in a metal matrix composite of Al, Cu, or Ag reinforced with single crystals of synthetic diamonds; as intended to be used for heat sinks or very efficient thermal conductors, spreaders, etc. Two different ways are described to predict the thermal conductivity of the composite: line tracing and surface-averaged Hasselman and Johnson. The effects of diamond volume fraction, size and, in particular, shapes on the composite thermal conductivity are discussed. 2007 Elsevier B.V. All rights reserved. Keywords: Packing of polyhedra; Metal matrix composites; Particle-reinforced composites; Diamonds; Thermal conductivity; Simulated annealing
1. Introduction Materials with a huge thermal conductivity are required for many technical applications: heat sinks, rocket nozzles, heat spreaders for electric and microelectronic applications. . . Diamond is the isotropic material with the highest thermal conductivity, but diamonds are extremely expensive if a very large size is required, it is indeed hard to machine, etc. On the other hand, silver, copper and aluminium are the metals with the largest thermal conductivities. Within the ExtreMat [1] project a large effort has been made to develop new materials with improved thermal conductivities for different temperature ranges (room temperature and high temperature: above 500 C). One of the most promising materials is copper, silver or aluminium reinforced with diamonds. Synthetic monocrystalline diamonds, with particle diameter of up to about * Corresponding author. Address: TECNUN, University of Navarra, Paseo de Manuel Lardiza´bal, 15, 20018 San Sebastia´n, Spain. Tel.: +34 943219877; fax: +34 943213076. E-mail address: jfl[email protected] (J. Flaquer).
0927-0256/$ - see front matter 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2007.03.016
600 lm are produced (at very large temperature and pressure, by hot isostatic pressing) and are actually twice cheaper per kilogram than silver (see, for example, SDG121 to 126 grades from Indus-Global [2], with a grit size from 30 to 60; or MBD-4 grade diamonds, with a mesh 70/80 from Henan Famous Industrial Diamond Co., Ltd. [3]). If an aluminium matrix is used, and the diamonds surface is not specifically prepared, the aluminium matrix bonds to square {0 0 1} surfaces of the diamonds, but does not stick to the hexagonal diamond faces (with {1 1 1} orientations) [4,5]. Therefore, the thermal resistance at both interfaces between the metal matrix and either the square or the hexagonal diamond faces are expected to be very different. There are models to estimate the thermal conductivity of composite materials reinforced with particles. For example, Maxwell’s approximation [6], but only for dilute spherical particles, which does not match with a high volume fraction of cubo-octahedron particles and, even worse, a perfect thermal connection is assumed between the matrix and particles. Maxwell’s approximation is, in our case, an
J. Flaquer et al. / Computational Materials Science 41 (2007) 156–163
157
upper bound to the actual thermal conductivity [7]. A more accurate model is that of Hasselman and Johnson [8], where an interface resistance is introduced, but spherical particles are still assumed. 2. Model Line tracing is used to compute the composite thermal conductivity in our model. In this case, diamonds are generated as truncated octahedrons and truncated cubes; see Figs. 1 and 2, respectively. 2.1. Particle generation A first diamond is located at the origin and oriented along the axis. Afterwards different diamonds are generated at random (not colliding with the previous ones) and with a random spatial orientation. The generation of random orientations, uniformly distributed in the space, can be found elsewhere [9,10] hi ¼ cos1 ð1 randÞ; wi ¼ 2p rand; p /i ¼ rand 2
Fig. 2. Enumeration of vertices, edges and faces for a truncated cube. The numbers without background are for vertices; bold for those at the front and regular for those at the back. The numbers on a square background stand for the faces; black on the front and pale grey for those at the rear. Edge directions are indicated by arrows.
ð1Þ z
where rand stands for random numbers uniformly distributed in the range 0–1 (most computer programming languages provide subroutines for the generation of such pseudo-random numbers). To specify a crystal orientation, three (Euler) angles should be provided, as shown in Fig. 3 and Eq. (1). Because of the cubic symmetry of the consid-
θ
n3
n1
y
n2
ϕ ψ
x Fig. 3. Schematic of the axis and (Euler) angles used to specify a crystal orientation. n1, n2 and n3 are normal to {0 0 1} faces.
ered polyhedra, wi can be reduced to only the one-fourth of the range p ð2Þ wi ¼ rand: 2 The size of each truncated octahedron, Ri, is generated with a normal distribution about a mean size, Rmean, with a given standard deviation, rR: Ri ¼ randn rR þ Rmean ; Fig. 1. Enumeration of vertices, edges and faces for a truncated octahedron. The numbers without background are for vertices; bold for those at the front and regular for those at the back. The numbers on a square background stand for the faces; black on the front and pale grey for those at the rear. Edge directions are indicated by arrows.
ð3Þ
where randn stands for a random number generated from a typified normal distribution (with a mean equal to 0, and a standard deviation of 1). Most computer languages provide for subroutines for such a purpose. Anyway, computer codes for the generation of quasi-random numbers with
158
J. Flaquer et al. / Computational Materials Science 41 (2007) 156–163
normal deviates can be found in [11]. This size refers to the radius of the sphere that circumscribes the original octahedron (it will be now be truncated). Different diamond-looking shapes are generated by truncation of the octahedra. The parameter a decides the height of the cut in comparison with Ri. a = 0 results in the original octahedron; the Kelvin Cell (with six little squares and eight large hexagons) is obtained for a = 1/3; a second tetrakaidecahedron (cuboctahedron) is obtained, formed by six squares and eight equilateral triangles, for a = 1/2; and cubes are obtained for a P 2/3. Note that according to Kelvin’s conjecture [12] his tetrakaidecahedron fills (tessellates) the space completely (a solution with slightly smaller surface was found by Weaire and Phelan [13] in 1993). But these identical tetrakaidecahedra should be very carefully located to fill completely the volume. One single Kelvin’s tetrakaidecahedron only fills one half of the smallest in which it is inscribed. pffifficube ffi The side of that cube is 2 2 l, where l is the length of any of the edges pffiffiffi of 3the tetrakaidecahedron; the cube volume is then pffiffiffi16 3 2 l and the volume of the tetrakaidecahedron is 8 2 l . So, it will be shown how difficult is to get large volume fractions of such particles when they are located at random. The computer code generates shapes at random, with a normal distribution about any given mean shape, amean, with any value for its standard deviation, ra: ai ¼ randn ra þ amean :
ð4Þ
This covers most of actual diamond shapes. In particular, all these polyhedra belong to the uniform polyhedra category: they have all their vertices identical, and – as a result of that – all of them located on a second sphere of smaller radius qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ri ¼ Ri 2a2i 2ai þ 1 () 0 6 ai 6 1=2; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð5Þ Ri ¼ Ri 6a2i 8ai þ 3 () 1=2 6 ai 6 2=3; p ffiffi ffi Ri ¼ Ri ð1 ai Þ 3 () 2=3 P ai P 1 as shown in Fig. 4. That feature helps to decide about the collision possibility in a very efficient way. 2.2. Collision detection Collision detection among polyhedra is done in four steps. 2.2.1. Spheres interference The first and fast check is to compare distance between every possible diamond couple, dij, and to compare with the sum of their radii: Ri is the radius of the sphere which circumscribes the polyhedron, ith (truncated octahedron). If the circumscribed spheres do not touch each other, there is no possibility of collision and Eq. (6) is satisfied: Ri þ RJ < d ij :
ð6Þ
Fig. 4. Maximum diameter of the polyhedron after truncation of a regular octahedron. Those polyhedra with regular polygonal faces are indicated with arrows.
If Eq. (6) is not satisfied, there is a possible collision and the analysis should go much deeper. 2.2.2. Vertex within another polyhedron The second check is whether a vertex of polyhedron i is within the volume of the polyhedron j. If all the faces are well oriented, it is to say: their normals point outwards from the polyhedron, the vertex is within the polyhedron when it is to the left of every face, which means that when its coordinates are substituted into the equation of the plane containing a face, the result is always negative. If no vertex of polyhedron i is within the polyhedron j, next we should check the other way: Does a vertex of j lie within polyhedron i? If no vertex is found within the other polyhedron we might feel that there is no collision, and this is the case in a 2D world (if we were working with polygons), but two polyhedra may also collide with each other when one edge, belonging to one of them, cuts through any face of the second one. Therefore, a collision is only rejected when additionally no edge is found to cut any face. 2.2.3. One shortcut: split plane Of course the computation of edge/face collision checks consumes a lot of time (most of the time) and any trick to save this verification deserves the effort. The intermediate step between the vertices detection and edge to face computation is to find a split plane between both polyhedra. If it is possible to find a plane that keeps all the vertices of the first polyhedron to one side and all of the second polyhedron to the other, they certainly do not collide. If we work with convex polyhedra (as it is the case) the candidates for
J. Flaquer et al. / Computational Materials Science 41 (2007) 156–163
split planes may be reduced to only their faces. Therefore, this intermediate check is also fast and very efficient, saving a lot of computational time.
E¼
n1 X n X
159
maxðjxi xj j; jy i y j j;jzi zj jÞ þ maxðLx ;Ly ;Lz Þ;
j¼1 i¼jþ1
ð7Þ 2.2.4. Edge/face intersection Every edge (line) cuts a face (plane) at one point and it should be checked if this intersection point is within the polygon (actual face) or outside. A simple code in Matlab [14] is used to generate the computer programme. For every edge of the first polyhedron (i) the intersection with the plane defined by every face of polyhedron j is calculated. Now it is decided if the intersection point is located within the actual edge (between their vertices) or not; and afterwards if the intersection point is within the face or not. From the intersection point, two vectors are built towards two consecutive vertices of that face. Vertices are chosen in a clockwise sequence. The cross product of both vectors is positive if the intersection point is to the right of this face side, and negative otherwise. The point is only inside the face if all the cross products, from all the consecutives couples of vectors to all the edges of this face, are negative. 2.2.5. Note In theory, it is not necessary to check for the position of the vertices inside the other polyhedron, but only to check for the edge of one polyhedron intersecting a face of the second one. In fact it is not required to check the other way: the edges of the second intersecting the faces of the first. That edge/face check does all the required collision verifications. Problems arise when the intersection point is too close to some edge or vertex: the cross product tends to zero. So, the second step is to compute the distance of the intersection point to the vertices of the second polyhedron and its edges. Thus, the above described checking for vertices is most useful to avoid singular cases when doing the edge/face comparisons, for example, an edge lying on a face. (The closer the particle packing, the more frequently these border line situations are encountered.) So it is a good idea to clean up the edge/surface collision detection as much as possible of the difficult cases by means of those previous checks. 2.3. Packing As it is well known, a generalized movement is formed by a spatial movement (3 degrees of freedom) plus three additional twists in the space. Particle densification is performed by simulated annealing, a well-known technique [15] that will be explained in the next section. The main idea is to achieve a diamond distribution that minimizes the objective function, E, given in Eq. (7), which can be considered as the ‘‘energy’’ of the whole system.
where xi, yi and zi stand for the coordinates of the centre of the ith polyhedron; Lx, Ly and Lz represent the lengths of the minimum size box, containing all the particle’s centres. A generalized random movement is generated for a randomly selected particle; each moduli being normally distributed. The movement is always accepted when the value for the expression (7) decreases. The minimization of Eq. (7) clusters the initial distribution of diamonds into a cubic shape. Of course, other energy functions can be used, for example E¼
n X n 1X d3 2 j¼1 i¼1 ij
ð8Þ
or E¼
p 6
n X
ðx2i þ y 2i þ z2i Þ3=2 ;
i¼1
where dij represents the distance between the centres of particles i and j. In the second case, the particle distance to the absolute centre is penalized according to the volume of the sphere of this size. But both functions tend to form a sphere of particles and, if a cubic box is pursued, the corners will be empty spaces and the volume fraction not as large as with previous Eq. (7). Note that none of the energies proposed in Eqs. (7) or (8) changes when a particle is rotated, maintaining its centre; so, twists do not affect the system energy (but they may help to pack particles in a later displacement; so they are permitted and so are all movement without energy change). 2.3.1. Simulated thermal annealing Simulated annealing uses, in our case, generalized movements. The optimization procedure tries to replicate the minimization of the target function, as the energy of a system that cools down slowly and follows Boltzmann’s law. Initially the system is at a ‘‘temperature’’ T0 and the system’s total energy is E0. We, arbitrarily, take an initial value for T0, for example, T0 = 1001 K and the system ‘‘energy’’ is computed according to Eq. (7), as explained before. The algorithm proceeds as follows. A temperature drop DT is defined, which, also arbitrarily will be taken as 1 K. The actual ‘‘temperature’’ is set to its initial value: T = T0. A particle is chosen at random. A generalized displacement (displacements plus rotations: Dx, Dy, Dz, D/, Dh, Du) is chosen, also at random. This particle movement is checked for collision with the other particles. If the movement is not physically possible (there is a collision with another particle), the particle is returned to its original position and another particle is chosen for movement. If the proposed movement is physically possible (there is no collision) we
160
J. Flaquer et al. / Computational Materials Science 41 (2007) 156–163
Kelvin’s truncated octahedron (tetrakaidecahedron). Other size distributions, bimodal distributions and shape distribution are possible, but the results are being computed and the results will be shown in a future paper. 3.1. Particle volume fraction
Fig. 5. Cluster of 26 diamond-like polyhedra with different sizes and shapes.
shall compute the ‘‘energy’’ modification produced within the system, DE. If the energy of the system is reduced, DE < 0, the change is always accepted; the temperature is reduced by a further step: T = T DT; and another particle is chosen for movement. If the system energy is not reduced, DE P 0, a random probability (uniformly distributed between 0 and 1), p, is generated and compared with kT 0 DE q ¼ exp ; ð9Þ T
Volume fraction is computed within a parallelepiped or box, parallel to the coordinate axes, as shown in Fig. 6. This box contains all the particle centres and it is the smallest one. Its maximum dimension appears in the expression that calculates the ‘‘energy’’ of the system. The volume fraction is obtained by dividing the volume of the diamonds contained within this box by the volume of the box. Observe that for a small system, like the one shown in Fig. 6, most of diamonds intersect the counting box. Among 27 particles, there is usually only one – at the centre of the box – not intersecting the counting box, and the other 26 intersect the box. Even for a much larger system, for example 125 particles, only 27 will not cut the counting box; that is only a 22% of the particles are completely embedded in the box. 3.2. Thermal conductivity Thermal conductivity is computed, for a given diamond packing by line tracing. Diamonds and the aluminium metal matrix have different thermal conductivities: kd = 1800 W/m K [5,16], and km = 245 W/m K [17]. If a diamond volume fraction Vd = 0.64 is introduced, Maxwell’s equation [6]
where k is an arbitrary constant (we usually use k = 10). If p < q the new situation is accepted – even when there is no reduction in the energy; the temperature is updated to T = T DT, and another particle is chosen for movement. The algorithm ends when, after a significant consecutive number of trials, nothing changes (no movement is accepted). In deciding a figure for this number of trials it helps to observe the evolution of the system energy that usually, after this number of trials, flattens down into a horizontal line. The smallest temperature is kept to be 1 K and never drops to zero. Also, to avoid too many changes with energy increments, at the beginning of the process, the value of q is limited to 0.3. Thermal agitation allows the system to escape from local minima (traps) and, if appropriately performed, to get actual system minima. It is shown to allow for the densification of the particles in a very efficient way and to obtaining, in this way, very large particle volume fractions, see Fig. 5. 3. Results In the example a uniform size distribution is used for the diamonds (R = 100 lm), and their shape is a perfect
Fig. 6. Cluster of 27 diamond-like polyhedra with identical sizes and shapes. Volume fraction of diamonds is 0.64 within the indicated parallelepiped.
J. Flaquer et al. / Computational Materials Science 41 (2007) 156–163
2 kkmd 1 V d þ kkmd þ 2 kc ¼ km 1 kkmd V d þ kkmd þ 2
ð10Þ
predicts a thermal conductivity for the composite of kc = 810 W/m K, which is clearly an upper bound for the actual case. Maxwell assumes perfectly bonded spherical particles diluted in the matrix. According to Torquato [7], Maxwell’s approximation gives good estimates for non-dilute conditions (provided that the spheres are well separated from each other), and also for different types of spheres. Additionally two interfacial thermal conductances should be assumed for the contacts between the {0 0 1} diamond faces (very well bonded to the Al matrix) and the contacts between the {1 1 1} diamond faces (with a very poor bond) and the metal matrix. The following values are used: h001 = 1 · 108 W/m2 K and h111 = 1 · 107 W/ m2 K, so, about 10 times worse for the interfacially poorly bonded {1 1 1} diamond surfaces. These values have been estimated from typical values found in the literature for other material systems [18,19]. The model of Hasselman and Johnson [8] introduces the effect of a thermal barrier resistance at the interface between particles and matrix 2 kkmd 1 krhd V d þ kkmd þ 2 þ 2 krhd kc ¼ km ; ð11Þ 1 kkmd þ krhd V d þ kkmd þ 2 þ 2 krhd where r is the radius of the particles, which are assumed to be spherical and of identical size, and h the interface conductance. Of course, Eq. (11) was not derived for cubooctahedral-shaped particles, nor for different thermal barrier resistances at different crystal surface orientations. Actual measurements of the composite conductivity, by Ruch et al. [5] for the same volume fraction of diamonds (Vd = 0.62) provide values ranging from kc = 130 to 670 W/m K, depending on the production route: mechanically assisted infiltration (squeeze casting) or gas pressure infiltration, respectively. The reason for such large difference is that in the squeeze cast samples, no chemical bonding (no interfacial carbides) is observed at all, neither on {1 1 1} nor {1 0 0} faces, whereas there is some chemical bonding between aluminium and {1 0 0} diamond faces in the gas pressure infiltrated sample. Line tracing draws lines, parallel to the preferred axis, and computes the different lengths on the metal matrix and within diamonds, Lm and Ld, respectively. In addition, it counts how many square faces, N4, are cut and how many hexagonal, N6, are intersected by the line. The thermal resistance (in units of K/W) is estimated for this line as the serial addition: 1 Lm Ld N 4 N6 Ri ¼ þ þ þ ; ð12Þ si km kd h001 h111 where si represents the cross-sectional area assigned to ith line.
161
The volume fraction of diamonds within line i can be expressed as V d;i ¼
Ld Ld ¼ ; Lm þ Ld Li
ð13Þ
where Li = Lm + Ld is the total length of the ith line. In analogy the fraction of {0 0 1} faces over the total particle surface, S001, in line i can be written as S 001;i ¼
N4 N4 ¼ : N4 þ N6 Ni
ð14Þ
Accordingly, Eq. (12) can be modified as 1 1 V d;i V d;i S 001;i 1 S 001;i Ri ¼ þ þ Li þ Li ni ; si km kd h001 h111 ð15Þ where ni = Ni/li is the density of intersected interfaces per unit length. The thermal resistance of the whole composite is estimated as the parallel resistance for all the box, R, lines, in the present case a regular set of 20 · 20 = 400 lines. X 1 1 P ¼ : ð16Þ Ri s i R si i
i
A lower bound for the thermal conductivity of the composite, kc is then estimated as 400 1 1 X 1 ¼ Li R Li si R i si i¼1 1 400 X 1 V d;i V d;i S 001 1 S 001 ¼ þ þ ; þ ni km kd h001 h111 i¼1
kc ¼
ð17Þ where the Li are assumed to be equal for all lines. Eq. (17) can be evaluated with relatively low computational effort for given phase arrangements. For the one shown in Fig. 6, the obtained composite thermal conductivity is kc = 255 W/m K. The percentage of {0 0 1} faces over the total particle surface, S001 = N4/(N4 + N6), depends on the truncation ratio 3a2 pffiffiffi () 0 6 a 6 1=2; 3a2 þ 3ð1 3a2 Þ pffiffiffi 2 4 3ð2 3aÞ ¼ pffiffiffi 4 3ð2 3aÞ2 þ 6½4ð1 aÞ2 2ð2 3aÞ2
S 001 ¼ S 001
ð18Þ
() 1=2 6 a 6 2=3; S 001 ¼ 1 () 2=3 P a P 1; as it is shown in Fig. 7. Alternatively, a mean interface thermal conductance can be estimated as h ¼ h001 S 001 þ h111 ð1 S 001 Þ:
ð19Þ
This value, which depends on a, can be introduced into Hasselman and Johnson’s Eq. (11) to estimate the effect of the diamond shapes. The results are shown in Fig. 8.
162
J. Flaquer et al. / Computational Materials Science 41 (2007) 156–163
Fig. 7. Percentage of {0 0 1} faces against the truncation ratio of the octahedron.
Observe that line tracing ‘‘extrudes’’ the composite material in spaghetti-like lines and then joins their ends, with each ‘‘column’’ being perfectly isolated with respect to its neighbours. For sufficiently long lines through statistically homogeneous composites the first term on the righthand side of Eq. (17) can be recognized as the classical Wiener [20] lower bound for the thermal conductivity of composites with perfect interfaces. The second term accounts for the interfacial resistances in a consistent way. The Wiener bounds are one-point bounds [7] that are rather slack for composites with ideal interfaces, and Eq. (17) is very conservative for particle reinforced composites with interfacial resistances because heat flow around obstacles is not accounted for. This obviously can lead to very low predictions in the case of a conductive matrix containing particles that are weakly conducting or have interfaces of low conductance. In the limit of high defect densities, ni, or low conductances h001 and/or h111, the predictions of Eq. (17) are dominated by the second term on the right-hand side. Accordingly, there is a marked effect of the particle size. The predicted overall conductivity tends to vanish for very small particles with finite interfacial conductances despite the presence of a contiguous conducting matrix, i.e., the trivial lower bounds are attained. The Hasselman–Johnson model gives identical results to Mori-Tanaka estimates [21] and can, accordingly be viewed as being related to two-point (Hashin–Shtrikman) bounds. For radii exceeding a critical radius of rc ¼
Fig. 8. Effect of the diamond sizes and shapes (truncation ratio) on the thermal conductivity of the diamond/Al composite, for a constant diamond volume fraction Vd = 0.64 with a thermal conductivity, kd = 1800 W/m K. The short-dashed horizontal line represents the thermal conductivity of the aluminium matrix, km = 245 W/m K. The black star stands for the line tracing computation with interface conductances of h001 = 1 · 108 W/m2 K and h111 = 1 · 107 W/m2 K, a particle radius of r = 50 lm and truncation ratio a = 1/3. The solid and long-dashed lines represent predictions of the Hasselman–Johnson model for r = 50, 100, and 200 lm, respectively; evaluated with averaged conductances following Eqs. (18) and (19). The blue star pertains to the same conditions as the line tracing model. The vertical line indicates the range of experimental measurements [5]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
4. Discussion Different procedures are proposed to compute or bound the thermal conductivity of metal matrix composites. These procedures have different pros and cons, as follows.
1 kd km h kd km
ð20Þ
[22], its results are lower estimates and for r < rc upper estimates. When combined with Eq. (20) it can provide estimates at low computational cost, but cannot account for details of the actual particle shape and uses a coarse approximation for the effects of the distribution of the interfacial conductances on the particle surfaces. Details of the microgeometry may obviously be studied by combining the finite element method with an appropriate micromechanical approach [23]. Among these, periodic homogenization with unit cells requires periodic particle arrangements, making windowing methods and embedding models more appropriate choices for geometries as shown in Figs. 5 and 6. However, as the particle volume fraction increases, the average distance between particles decreases, which requires meshes with high number of small elements. This leads to models with high numbers of degrees of freedom that pose major demands in computer capacity. 5. Conclusions A fast and efficient way of generating microstructures composed of monocrystalline diamonds was described. A set of cubo-octahedral geometries has been used. The geometrical advantages of using uniform polyhedra have been exploited, and the way of computing possible collisions has been detailed. It may be used to model different properties
J. Flaquer et al. / Computational Materials Science 41 (2007) 156–163
of MMCs reinforced with single crystal synthetic diamonds, as used within ExtreMat project for improved thermal conductivity materials. A shape selection of diamonds, increasing the amount of {0 0 1} faces, improves the composite thermal conductivity from kc = 525 W/m K for Kelvin tetrakaidecahedra to 700 W/m K for cubical diamonds. A 33% improvement in the composite thermal conductivity may be expected in this way. The same improvement is obtained using larger diamonds (and reducing the interface area per unit of volume) to r = 180 lm, as follows from Fig. 8. Acknowledgements This work has been carried out within the framework of the European Integrated Project ‘‘ExtreMat’’ (New Materials for Extreme Environments). Thanks are given to the European Commission for funding this Integrated Project under FP6 (Contract No. NMP3-CT-2004-500253), as well as by the Spanish Ministry for Science and Technology (FTN2003-08228-C03-02). It only reflects the view of the authors and nor the European Community nor the Spanish Ministry for Science and Technology are liable for any use of the information contained therein. References [1] ExtreMat, Integrated European Project: New Materials for Extreme Environments, funded under FP6 (NMP3-CT-2004-500253), 2004. [2] Indus Global Superabrasives, India,
163
[4] S. Kleiner, F.A. Khalid, P.W. Ruch, S. Meier, O. Beffort, Scripta Mater. 55 (2006) 291–294. [5] P.W. Ruch, O. Beffort, S. Kleiner, L. Weber, P.J. Uggowitzer, Compos. Sci. Technol. 66 (2006) 2677–2685. [6] J.C. Maxwell, A Treatise on Electricity and Magnetism, 3rd ed., Oxford University Press, Oxford, 1904. [7] S. Torquato, Random Heterogeneous Materials Microstructure and Macroscopic Properties, Springer, New York, 2002. [8] D.P.H. Hasselman, L.F. Johnson, J. Compos. Mater. 21 (1987) 508– 515. [9] R. Rodrı´guez-Martı´n, I. Ocan˜a, A. Martı´n-Meizoso, Mater. Sci. Forum 482 (2005) 171–174. [10] D. Zwillinger, Standard Mathematical Tables and Formulae, 30th ed., CRC Press, Boca Rato´n, FL, 1996. [11] W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes The Art of Scientific Computing, Cambridge University Press, Cambridge, MA, 1986. [12] W. Thomson, Acta Math. 11 (1887) 121–134. [13] D. Weaire, R. Phelan, Phil. Mag. Lett. 69 (1994) 107–110. [14] Matlab 7.0, The MathWorks Inc., 1984. [15] P.S.M. van Laarhoven, E.H.L. Aarts, Simulated Annealing Theory and Applications, Kluwer Academic Publishers, Dordrecht, 1992. [16] J.E. Field, The Properties of Natural and Synthetic Diamond, 3rd ed., Academic Press, London, 1992. [17] CES Selector Version 4.5, Cambridge Engineering Selector, Granta Design Ltd., Cambridge, UK, 2004. [18] R.J. Stoner, H.J. Maris, Phys. Rev. B 48 (1993) 167373. [19] H.-K. Lyeo, D.G. Cahill, Phys. Rev. Lett. 96 (2006) 186101. [20] O. Wiener, Abh. Math.-Natur. Sa¨chs. Akad. Wiss. Leipzig 32 (1912) 509. [21] Y. Benveniste, J. Appl. Phys. 61 (1987) 2840. [22] D. Duschlbauer, Computational Simulation of the Thermal Conductivity of MMCs under Consideration of the Inclusion-Matrix Interface, VDI-Verlag, Du¨sseldorf, 2004. [23] H.J. Bo¨hm, A short introduction to continuum micromechanics, in: H.J. Bo¨hm (Ed.), Mechanics of Microstructured Materials, vol. 1, Springer, Vienna, 2004.