Strojniški vestnik - Journal of Mechanical Engineering 51(2005)7-8, 391-398 UDK-UDC 536.2 Izvirni znanstveni članek - Original scientific paper (1.01) Thermal Characterisation of Rectangular Cooling Shapes in Heat Generating Mediums - A Three-Dimensional Investigation Jaco Dirker1, Arnaud G. Malan2, and Josua .P. Mey 2 er 1Department of Mechanical Engineering, Rand Afrikaans University, Johannesburg, South Africa 2 Department of Mechanical and Aeronautical Engineering, University of Pretoria, Pretoria, South Africa, * jmeyer@up.ac.za Abstract The optimum aspect ratios of uniformly distributed embedded rectangular cross-sectioned solid cooling inserts in a heat-generating medium are investigated. Numerical investigations were performed to determine and characterise how various geometric and thermal parameters determine these optimum shapes. Introduction The current trend in power electronics is to increase the power conversion capability of circuits while reducing their size. One such a way is by means of three-dimensional integration of discrete components into multifunctional modules and the creation of standardized building blocks [1]. This can be done by using planar type structures or modules with various layers [2]. In order to satisfy future thermal demands associated with, for instance, integrated power electronic devices, the focus is starting to shift toward innovative design of the internal structure of power modules to assist in heat extraction, while maintaining high levels of electromagnetic performance and efficiency. Such an integrated power electronics system requires advances in different technologies, which depend upon finding solutions to deal with the multi-disciplinary issues in materials, electromagnetic compatibility and thermal management. Purpose Due to the manufacturing method proposed for constructing the new type of integrated power electronic modules, the cooling inserts would need to have rectangular cross-sections. It is the aim of this paper to describe the optimum aspect ratio of the cooling insert’s cross section such that the peak temperature within the heat-generating medium would be minimised. Little reference material was found in literature that describes the optimisation of such geometries. By characterising the influence of various parameters on the optimum cooling structure shape, the design of such embedded cooling systems are enabled. Procedure A numerical approach was followed in the study, as it would have been too time-consuming and costly to do the optimisation processes experimentally. Consider Fig. 1 that gives a representation of a proposed method of assisting heat flow from within a heat-generating medium to the surroundings. The sectioned view shows a generalised distribution of evenly spaced identical rectangular cross-sectioned heat extraction inserts running parallel to the z direction. Arrows indicate the flow of heat from the heat-generating medium via the heat extraction inserts to the surrounding. If such a set-up is large enough in the x and y directions, or insulated from the ambient, certain assumptions can be made. If each heat extraction structure is exposed to an isothermal uniform ambient in an identical way, rectangular regions drawn around each cooling structure such as those in Fig. 1 (dotted lines), would have identical temperature distributions. Heat-generating material Cooling insert Figure 1 Cross-sectioned representation of a heat generating material with embedded cooling structures 391 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 391-398 Nomenclature A x-y view cross-sectional area [m2] 7 thermal conductivity ratio [dimensionless] a x-y view aspect ratio [dimensionless] A half centre-to-centre x-directional offset of aC,rel relative x-y view aspect ratio of cooling insert neighbouring cooling inserts [m] [dimensionless] a half x directional dimension of cooling insert [m] C heat equation coefficient B half centre-to-centre x directional offset of CG vector of heat gain coefficients neighbouring cooling inserts [m] GM ration between volumetric heat generation and b half y directional dimension of cooling insert [m] thermal conductivity [K/m2] Z half z directional dimension of structure [m] k thermal conductivity [W/mK] M matrix of heat equation coefficients Subscripts M2 number of nodes in the x direction B relative rear directional (negative z direction) [dimensionless] C cooling insert N2 number of nodes in the y direction D whole domain [dimensionless] E relative easterly direction (positive x direction) q heat flux at exposed face of cooling insert [W/m ] F relative frontal direction (positive z direction) G heat gain or heat loss qM volumetric heat generation density [W/m3] max maximum R thermal interface resistance [m2K/W] min minimum T temperature [K] N relative northerly direction (positive y direction) x Cartesian coordinate S relative southerly direction (negative y direction) y Cartesian coordinate T relative current node indication z Cartesian coordinate W relative westerly direction (negative x direction) 0 reference cooling insert temperature Greek and special symbols a fraction of volume occupied by cooling inserts [dimensionless] Furthermore due to the even spacing of the cooling structures and the equal exposure to the surroundings in both sides in the z direction, the temperature field about one such cooling structure is symmetric in all three Cartesian directions. The rectangular region around such a heat extraction structure can thus be divided into quadrants, each of which would be sufficient to represent the temperature field about such a heat extraction structure. Such a three dimensional representative quadrant domain is shown in Fig. 2 where the heat extraction insert is exposed to the surroundings on its positive z side face. In order to obtain the steady state three-dimensional temperature distribution in such a quadrant domain, a cell-centred finite volume numerical approach was used. The schematic representation of the grid used for this is shown in Fig. 3 . The number of nodes shown here are not necessarily the number used during simulations. An important point to note is that no nodes were defined on the interface between the cooling structure and the heat-generating medium. This is due to discontinuities that appear in solutions for cases with small or no contact resistance at the interface. The grid used was localised uniform, which means for each block shown in Fig.3 , the grid spacing was uniform, but not necessarily uniform across the entire mesh. Uniform grid spacing was also used in the z direction. All boundaries were defined as being adiabatic, except the positive z side face of the heat extraction insert, which had a fixed reference temperature of T0 at x = 0 and y = 0, and a uniform heat flux to the surrounding. By applying the principle of conservation of energy, the steady state heat flux, qC [W/m2], on this surface can be obtained from the dimensions of, and volumetric heat generation density, qM [W/m3] within the representative domain: q C = - q ''M(AB-ab)Z ab (1) Here A, a, B, b, and Z are the dimensions of the representative domain as defined in Fig. 2 . Heat generating material y U—>x z Z Cooling ^^ structure Figure 2 Schematic of the representative domain A B a b 392 Dirker J. - Malan A.G. - Meyer J.P. Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 391-398 interface with thermal contact resistance (0 ;0 ;0) TC (A ;B ;Z) (A ;b ;Z) (A ;0 ;Z) cooling structure with uniform heat flux out Figure 3 Meshing scheme used z heat-generating medium Variable Definition The aspect ratios of the overall domain cooling structure in the xy-plane were respectively as: A aD = — >1 B aC = — b and the defined (2) (3) Only aD > 1 were investigated, as the mirror behaviour for aD<1 is equivalent. For instance ad = 2 would give the same results as aD = 0.5 by just rotating the x and y axes by 90°. The cross-sectional area of the entire representative domain and cooling structure in the xy view can be written respectively as: AD = AB AC = ab = W (4) (5) The fraction of the entire domain volume used for heat extraction purposes are denoted by a. Five variables namely AD [m2], aD [dimensionless], aC [dimensionless], a [dimensionless], and Z [m], define the geometry of the representative domain fully. Due to geometric constraints, the cooling insert aspect ratio, aC, has maximum and minimum allowed values namely: a C,min a C,max D = aD a (6) (7) The ratio between the thermal conductivity of the cooling structure, kC [W/mK] and thermal conductivity of heat-generating material kM [W/mK] was defined as: r= kC M (8) The ration between the volumetric heat-generation and thermal conductivity of the heat-generating material [K/m2] can be written as: GM = qm (9) Numerical Method Numerical simulation involves three separated stages namely pre-processing, processing or solution, and post-processing. Pre-processing include the definition of domain dimensions, thermal properties and conditions, mesh creation, and solution method specification. Due to the large number of anticipated simulation cases that would be needed to optimise cooling structure shapes and distribution, the conventional pre-processing stage would become very time consuming when using commercially available numerical simulation software packages. For this reason it was opted to create computer code that would automatically do the pre-processing functions and would then feed the problem information straight into a solution procedure or algorithm. A fully implicit solution algorithm approach was followed. This approach has the advantage of being computational efficient in the case of relatively small meshing schemes, and highly predictable concerning the running time needed to solve the temperature field. It is only dependent on the number of nodes in the domain and is not influenced by the physical conditions of the domain, as would be the case if a fully explicit procedure were used. In order to use an implicit method, it is necessary to express the temperature at a node in terms of its neighbouring temperature values by means of different coefficients [3]: CTT = CFTF + CSTS + CWTW + CETE + CNTN + CBTB + CG …(10) Subscripts F, S, W, E, N, and B refer to the six reference directions namely “front”, “south”, “west”, “east”, “north”, and “back” respectively. Subscript T refer to temperature node under consideration, and G indicates the heat gain or heat loss of the particular node. If there is N number of nodes, this will result in N number of equations. It is possible to express this system of equations as: MT = CG (11) Here M is a N-by-N coefficient matrix, CG is a vector containing the heat gain or loss coefficient for each node, and T is a vector containing the temperature values of all the nodes in the domain. k M y Thermal characterisation of rectangular cooling shapes in heat generating mediums 393 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 391-398 In order to construct such a matrix it is necessary to decide before hand on a numbering scheme for the nodes such that each node has a unique address. For a structured numbering scheme (shown in Fig. 4 for a 3x3x3 mesh) where there are M2 number of nodes in the x direction and N2 number of nodes in the y direction, M has a banded structure as shown in Fig. 5 with seven diagonals that contain non-zero entries. Each of these diagonals contains the relevant coefficients defined in Eq. (10). The coefficient matrix is thus a very sparsely populated matrix with the majority of entries being zero. The conventional classic method of solving the temperature vector would be to use Gauss-Jordan elimination to obtain the inverse of the matrix, but due to the high density of zero-entries there are more efficient ways of solving the temperature vector. One such a method is LU-decomposition of the coefficient matrix into an upper and lower triangular matrix. From this back- and forward-substitution can be used to obtain the temperature solution. Due to the banded structure of M, the traditional LU-decomposition algorithm can be adjusted to reduce the number of computations needed. This is done by only applying the decomposition algorithm within the band of M ranging from diagonals B to F. This is a valid adjustment as the entries of lower and upper triangular matrices corresponding to the regions outside this band only contain zero entries anyway. The number of computations needed can further be reduced if the bandwidth of M is decreased. This can be done by using a different numbering scheme for the nodes in the domain according to the reverse Cuthill-McKee algorithm. Refer to Fig. 6 for the percentage reduction in the bandwidth (compared with the structured numbering scheme) for three different size matrices at different number of nodes in the z direction after this algorithm was applied. It was found that for a structured two-dimensional mesh, where there is only one node in the z direction, there is no reduction in the bandwidth of the matrix. It is the most advantageous to apply the reverse Cuthill-McKee algorithm where the mesh has a width of two nodal point in the z direction. In the current study 10 nodes were used in the z direction. Using more than 10 nodes in the z direction, resulted in differences of less than 0.1% in solved temperature distributions. Numerical Results By running a sequence of simulations, the optimum cooling structure shape can be found for a particular representative domain geometry and thermal condition. It was found that for aD ?1, the lowest peak temperature is always associated with a cooling structure aspect ratio, aC, in the range of [aD; aC,max]. A convenient way of normalising aC is by defining a new relative cooling aspect ratio, aC,rel: aC,rel = ac -aD aC,max -aD (12) When aC,rel = 1, it means that the cooling insert is at its maximum aspect ratio as shown in Fig. 7, while when aC,rel = 0 the cooling insert has the same aspect ratio as the representative domain. It was found that the maximum temperature in the domain, Tmax, is directly proportional to GM and that an increase or decrease in the cooling structure temperature, TC, is translated directly into an identical increase or decrease in Tmax. Mathematically this can be expressed as: Tmax -TC ?GM (13) This means that neither GM nor TC has any influence on the shape of the optimum cooling geometry. 7 8 9 16 4] 5 6 13 17 14 18 15 11 26 23 27 24 12 * 1 2 3 10 x k=1 (front face) k =2 Figure 4 Structured numbering scheme used 19 20 21 k =3 (back face) TE W NB 1 y Figure 5 Banded structure of the coefficient matrix 70 60 50 40 30 20 10 0 3x3xkMesh 4x4xkMesh 5x5xkMesh 1 23 45 k (Number of nodes in z direction) Figure 6 Reduction in matrix bandwidth after applying the Reverse Cuthill-McKee algorithm S F 394 Dirker J. - Malan A.G. - Meyer J.P. Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 391-398 Optimisation Results for aD = 1 A domain with a square cross-sectional area is applicable to an arrangement where the cooling structure spacing in both Cartesian directions is identical. Fig. 8 shows that for the case where there is no thermal contact resistance, the physical size of the cooling insert has little influence on its optimum shape. Similarly, in Fig. 9 , the same is shown to be true of the influence of the depth of the domain. For cases where less than approximately 70% of the domain is occupied by cooling, the optimum cooling insert shape is found to be a flat continues “plate” as shown in Fig. 7. Also, it was found for cases with no thermal interface resistance, that the actual values of the thermal conductivities did not play any role, but rather that the ratio between the thermal conductivities, y, was of concern. In Fig. 10 the optimum cooling shapes for various ratios, where no thermal contact resistance is present, is given. When however thermal contact resistance, R [m K/W] is introduced, it is found that none of the above mentioned trends are valid. From Figs. 11-13 it is shown that the presence of contact resistance results in the relative optimum aspect ratio deviating from the aC,rel = 1 line at an earlier a value than in the case without the presence of contact resistance. This means that the optimum cooling shape for larger cooling fraction is no longer that of a flat continues plate. Also as is shown in Fig. 14 , that the physical sizes of the thermal conductivities influence the optimum cooling aspect ratio when interface thermal resistance is present. ¦Heat generating medium •Merged cooling structures x IHIIIIIMIIi a D =1 Z =1 m Figure 7 Physical meaning for the case where aC is at its maximum 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 : 1 W/mK = 1 W/mK R =0 m2K/W AD values [m2]: 0 0.2 0.4 0.6 0.8 1 Fraction of Domain used for Cooling (a) Figure 8 The influence of the physical size of the cooling structures on the optimum cooling insert shape for a mesh thickness of 1 m. 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 IIIIIIMIIIH , values [m]: -1 -0.5 -0.1 -0.05 -0.01 -0.05 -0.001 aD =1 R =1 m2 =1 W/mK = 1 W/mK 0m2K/W 0 0.2 0.4 0.6 0.8 1 Fraction of Domain used for Cooling (or) Figure 9 The influence of the z direction depth on the optimum cooling insert shape for a domain area of 1m2. 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 - 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Fraction of Volume Occupied by Heat Extractors (a) Figure 10 Influence of y on aC,rel, optimum for AD 1m2 Z = 0.001 m, and R = 0 m2K/W 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 - 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Fraction of Domain used for Cooling (a) Figure 11 Influence of thermal contact resistance on the optimum cooling shape for AD = 1 m2 and Z= 0.1 m A D k C k C M Thermal characterisation of rectangular cooling shapes in heat generating mediums 395 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 391-398 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 R =0.1 m2K/W R =0.01 m2K/W R =0.001 m2K/W R =0.0001 m2K/W R =0 m2K/W------- aD =1 AD = 0.001m2 Z =0.1 m ,kC = 1W/mK kM = 1 W/mK 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Fraction of Domain used for Cooling (a) Figure 12 Influence of thermal contact resistance on the optimum cooling shape for AD = 0.001 m2 and Z 0.1 m Optimisation Results for aD = 2 A domain with an aspect ratio of 2 represents a situation where the centre-to-centre distance between two adjacent cooling structures in one Cartesian coordinate direction is twice the centre-to-centre distance of two neighbouring cooling structures in the other Cartesian direction. Fig. 15 gives the behaviour of the optimum cooling shape for aD = 2 at various R-values. It was found that as before the presence of thermal contact resistance influences the optimum shape. It was also found that unlike aD = 1, the optimum cooling shape never has the same aspect ratio as the representative domain. It has a higher tendency to have a flat plate geometry as an optimum cooling shape. For smaller cooling structures it was found that the optimum geometry graphs fan out a bit as indicated in Fig. 16 for different R values. The influence of the conductivity ratio was found to be small. 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Fraction of Domain used for Cooling (a) Figure 13 Influence of thermal contact resistance on the optimum cooling shape for AD = 0.001 m2 and Z= 0.01 m 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 a D= 2 A D =0.001 m2 Z= 1 m kC = 140 W/mK k=2.1 W/mK M R values: [m2K/W] 0 0.0001 A 0.001 X 0.01 O 0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Fraction of Volume Occupied by Heat Extractors (a) Figure 15 Influence of R on the optimum cooling shape for aD = 2. 1 0.9 0.8 0.7 0.6 0.5 0.3 0.2 0.1 0 0 0.2 0.4 0.6 0.8 1 Fraction of Domain used for Cooling (a) Figure 14 Non-zero R for a constant ratio of thermal conductivity for AD= 1 m2 and Z = 0.001 m 0.00001 0.0001 Representative Domain Area (AD ) [m2] 0.001 Figure 16 Influence of AD for aD = 2 at a = 0.524 396 Dirker J. - Malan A.G. - Meyer J.P. Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 391-398 Optimisation Results for aD = 5 As was found for aD = 2, the case where aD = 5, has an even higher tendency toward optimum cooling shapes of a flat plate. From Fig. 17 it can also be seen that the influence of the thermal contact resistance is also not as severe as before. Similarly the physical size of the cooling structure, AD, has a lesser influence. 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 ¦ * m*- R values [m2K/W]: -0 -0.001 - 0.002 0.005 0.01 -0.02 0.05 0.1 aD =5 A D = 0.001 m2 z =1m k C kM 140 W/mK = 2.1 W/mK 0 0.2 0.4 0.6 0.8 1 Fraction of Domain used for Cooling (a) Figure 17 Optimum cooling shape aspect ratios for aD = 5 and various R values. 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 ¦ >¦¦¦¦¦=**» R values [m2K/W]: -0 -0.001 -0.002 0.005 0.01 -0.02 0.05 0.1 a D = 10 AD = 0.001 m2 z =1m kC = 140 W/mK kM = 2.1 W/mK 0 0.2 0.4 0.6 0.8 1 Fraction of Domain used for Cooling (a) Figure 18 Optimum cooling shape aspect ratios for aD = 10 and various R values. Optimisation Results for aD = 10 A similar trend is found to be true for higher representative domain aspect ratios as is shown in Fig. 18. In this case it was observed that almost through the entire spectrum of domain fractions used by cooling the optimum cross-sectional shape was that of a flat plate. Dimensional Scaling It is also important to determine when it is advantageous to optimise the cross sectional aspect ratio of the extractor inserts in order to accommodate higher heat generation densities while maintaining a certain maximum temperature within the representative domain. Fig. 19 shows, for a case with no internal interface thermal resistance, that the effectiveness of the heat extractor becomes less dependent on its cross sectional shape as the ratios of AD to Z2 decreases. For low ratios of AD to Z2, where the domain of a single insert becomes increasingly slender, the advantage of optimising the cross section in order to support higher heat-generation levels diminishes. The advantage of inserting heat extractors however does not diminish, but rather the effectiveness thereof for slender domains are more dependent on the fraction of the total volume occupied by the heat extraction system a, rather than the cross sectional shape thereof. It would thus be sensible in such cases to use a cross-sectional geometry that is easier to manufacture or a geometry that would conform to possible other restrictions. Similar trends were obtained for beryllium oxide and synthetic diamond inserts. However, for the case where thermal contact resistance is present it was found that unlike the trend shown in Fig. 19, the influence of the cross-sectional shape remains present for low AD to Z ratios. 50 45 40 35 30 25 20 15 10 5 0 0.000001 0.0001 0.01 AD /Z 2 1 100 Figure 19 Additional allowable increase in heat generation density after geometric optimisation Conclusions It was found that there are 7 variables that determine what the optimum embedded cooling insert cross-sectional aspect ratio will be, namely aD, AD, a, Z, R, kM, and kC . When no thermal contact resistance is presents it seems as if this list can be reduced to aD, «, andy. The actual temperature values within the heat-generating medium will however be dependent on all the above given parameters as well as the reference cooling temperature, T0 and the volumetric heat generation, q''' . The exact influence of the thermal contact resistance is difficult to describe from the presently available data, and more investigation will be needed in this regard. Thermal characterisation of rectangular cooling shapes in heat generating mediums 397 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 391-398 Indications are that as the domain aspect ratio is increased, the tendency for the optimum cooling structure to be a continues flat plate increases. References [1]Van Wyk JD, Strydom JT, Zhao L, Chen R, Review of the development of high density integrated technology for electromagnetic power passives, Proceedings of the 2nd International Conference on Integrated Power Systems (CIPS), 2002, pp. 25-34 [2] Barbosa P., Lee F.C., Van Wyk J.D., Boroyevich D., Scott E., Thole K., Odendaal H., Liang Z., Pang Y., Sewell E., Chen J., and Yang B., An overview of IPEM-based modular implementation for distributed power systems, Center for Power Electronics Systems (CPES) Power Electronics Seminar, Virginia Tech, 674 Whittemore Hall, Campus Mail Code 0179, Blacksburg, VA 24061, 2002, pp 70-76 [3] Patankar S.V., Numerical heat transfer and fluid flow, Hemisphere (Washington D.C.) , 1980 398 Dirker J. - Malan A.G. - Meyer J.P.