Get documents about "Chapter 4 Two-Dimensional Analysis
Document Sample
Chapter 4
Two-Dimensional Analysis
Table of Contents
4.1. Introduction......................................................................................................36
4.2. Cooling Structure Distribution........................................................................36
4.3. Thermal Governing Equation and Conditions ...............................................37
4.4. Geometric Definitions and Restrictions .........................................................39
4.4.1. Definitions......................................................................................................39
4.4.2. Restrictions on the Cooling Structure Aspect Ratio ........................................40
4.5. Analysis Approaches ......................................................................................41
4.6. Analytical Analysis – Generalised Cooling Structure ...................................41
4.7. Numerical Analysis – Generalised Cooling Structure...................................44
4.7.1. Meshing .........................................................................................................45
4.7.2. Discretised Governing Equation.....................................................................46
4.7.3. Tri-Diagonal-Matrix-Algorithm (TDMA) Solution Method.................................48
4.7.4. Numerical Results..........................................................................................50
4.7.5. Optimisation Results ......................................................................................51
4.7.5.1. Domain with a Square Cross-Sectional area (aD = 1) ...................................................52
4.7.5.2. Domain with aD = 2 ........................................................................................................53
4.7.5.3. Domain with aD = 5 ........................................................................................................55
4.7.5.4. Domain with aD = 10 ......................................................................................................56
4.7.5.5. Behaviour of the Optimum Cooling Structure Aspect Ratio ..........................................57
4.7.6. Maximum Temperatures for Optimum Geometry ...........................................58
4.7.7. Summary of Results up to This Point .............................................................59
4.8. Numerical Analysis – Solid Cooling Structure ..............................................60
4.8.1. Meshing and Numerical Solution Method.......................................................60
4.8.2. Discretised Governing Equations ...................................................................62
4.8.2.1. Discretised Coefficients for Heat-Generating Medium Nodes ......................................62
4.8.2.2. Discretised Coefficients for Cooling Structure Nodes ...................................................63
4.8.3. Numerical Results..........................................................................................64
4.8.4. Verification of Results ....................................................................................65
4.8.5. Summary .......................................................................................................68
4.9. Conclusion .......................................................................................................68
4.10. Nomenclature...............................................................................................69
4.10.1. General Symbols ...........................................................................................69
4.10.2. Greek Symbols ..............................................................................................70
4.10.3. Subscripts......................................................................................................70
4.11. References ...................................................................................................70
Chapter 4 – Two-Dimensional Analysis 35
4.1. Introduction
In this chapter the foundation is laid down for the theoretical investigation of the cooling of a
homogeneous material with uniform heat-generation via rectangular cross-sectioned embedded
cooling structures. The aim of this section is to formulate general geometric restrictions and
definitions. Also, different methods of investigating and describing the temperature field within
the heat-generating medium and or the cooling structure inserts are considered.
Ultimately it is the purpose of these definitions to assist in the process of optimising the
distribution of the geometric shape of cooling insert to obtain the lowest possible peak
temperature within the heat-generating medium for a parallel-configured heat-extraction system.
4.2. Cooling Structure Distribution
Consider Figure 4.1 showing a representation of a sectioned heat-generating medium with a
general distribution of evenly spaced identical inserted rectangular cross sectioned cooling
structures running parallel to the z direction:
Heat
generating
material
Rectangular
embedded
y cooling
structures
x
Figure 4.1 Cross-section of a heat-generating material with embedded cooling structures
If such a composite structure consisting of uniformly spaced cooling inserts and heat-generating
material is large enough in the x and y directions, or insulated from the surrounding, certain
assumptions may be made concerning the temperature distribution. If each cooling structure has
identical material properties and dimensions, and is exposed to an isothermal ambient in an
identical way, control volumes drawn around each cooling structure, such as those shown in
Figure 4.1, would each have identical temperature distributions.
Due to the even spacing of the cooling structures, the temperature field about one such cooling
structure is symmetric in both the x and y Cartesian directions. The control volume around such a
cooling structure can thus be sub-divided into quadrants, each of which would be sufficient to
represent the temperature field about such a cooling structure. Such a three-dimensional quadrant
domain is shown in Figure 4.2 indicating the dimensions, which will be referred to in this and
subsequent chapters.
36 Chapter 4 – Two-Dimensional Analysis
2
2
Representative
domain
2
y
x
z Cooling
structure Heat-generating
material
Figure 4.2 A representative domain of the temperature distribution about a cooling structure
insert.
Dimensions and represent half the centre-to-centre distances between neighbouring cooling
structures in the x and y directions respectively, while dimensions and represent the half -
width and -height of a cooling structure in the x and y directions respectively. represents the
half-depth of the domain or control volume of interest.
4.3. Thermal Governing Equation and Conditions
For the application of interest certain approximations are made to ease the theoretical
investigation approach. Two of these have already been mentioned, namely that the heat-
generating medium is assumed to be homogenous in respect to material properties and uniform in
terms of thermal heat-generation.
The governing equation describing the three-dimensional temperature distribution in a heat-
generating solid can be written in differential form as:
∂ 2T ∂ 2T ∂ 2T q '''
+ + + =0 (4.1)
∂x 2 ∂y 2 ∂z 2 k
Here k [W/mK] is the thermal conductivity of the material and q ''' [W/m3] is the thermal
volumetric heat-generation rate within the material. For a solid cooling structure this equation is
also valid but with q ''' equal to zero (no heat-generation).
Due to the symmetric nature of the temperature profile in the control volume around a cooling
structure insert in the x and y directions, as explained earlier, adiabatic (zero heat transfer)
boundaries can be applied on certain domain faces when analysing the representative corner
quadrant domain theoretically. Refer to Figure 4.3.
Chapter 4 – Two-Dimensional Analysis 37
∂T
=0
∂y
∂T
=0
∂x
∂T
=0
∂x
y
x
z
∂T
=0
∂y
Figure 4.3 Adiabatic boundaries in the x and y directions
These boundary conditions can be written as:
∂T
=0 (4.2)
∂x x =0
∂T
=0 (4.3)
∂x x=
∂T
=0 (4.4)
∂y y=0
∂T
=0 (4.5)
∂y y=
From the principle of conservation of energy is follows that for steady-state conditions the rate of
heat-generation within the medium is equal to the rate at which heat is carried by the cooling
structure to the environment. If it is then assumed that the heat-generating medium is insulated
on the outside in the z-direction (being the worst-case-scenario), all heat would need to be
transported across the exposed xy faces of the cooling structure to the environment.
If it is assumed that all the heat transported to the environment from the exposed cooling structure
faces occurs at a uniform rate per unit area, the boundary conditions in the z direction can thus be
described as:
∂T
=0 (4.6)
∂z x> ∩ y>
∂T
= qC''
'
(4.7)
∂z x< ∩ y<
The uniform heat flux at which heat leaves the exposed cooling structure surface, qC' [W/m2],
'
follows from the principle of conservation of energy and can be written as:
38 Chapter 4 – Two-Dimensional Analysis
qM' (
''
− )
qC' = −
'
(4.8)
Here qM' [W/m3] refers to the volumetric heat-generation within the heat-generating medium. In
''
the rest of the test, M will refer to the heat-generating medium, and C to the cooling structure
(whether being a solid state heat-extractor, or a microchannel). Graphically, for a solid cooling
structure insert where heat is transferred to the environment in two directions, this can be
represented as shown in Figure 4.4.
∂T
=0
∂y
∂T heat-generating ∂T
=0 medium =0
∂z ∂z
∂T ∂T
= qC'
'
cooling structure = qC'
'
∂z ∂z
y
z ∂T
=0
∂y
Figure 4.4 Thermal boundary conditions in the y and z directions
Up to this point no mentioning has been made of the condition at the interface between the heat-
generating medium and cooling structure. If the cooling structure is solid, the heat transfer,
q [W], across the interface can be described with the aid of a thermal interface resistance, Rif
[m2K/W] between the two mediums:
Aif ∆T
q= (4.9)
Rif
Here Aif [m2] refers to the contact interfacial area between the mediums. If the cooling structure
is a channel filled with a fluid, the heat transfer can be described by a convective heat transfer
coefficient, h [W/m2K]:
q = hA∆T (4.10)
In the first case T represents the temperature difference across the solids-interface while in the
second case it refers to temperature difference between the channel wall and fluid temperature in
the centre of the channel.
4.4. Geometric Definitions and Restrictions
4.4.1. Definitions
Even though some dimensions have been defined in Figure 4.2, it is possible to define additional
variables which would be more geometrical meaningful to the investigation.
The aspect ratio of the domain which describes the relative distribution of cooling structures in
the x and y directions can be defined as:
Chapter 4 – Two-Dimensional Analysis 39
aD = (4.11)
For an aD value of 1, it follows that the centre-to-centre distance between neighbouring cooling
structures in the x an y directions are identical. For an aD value of 2 the centre-to-centre distances
in the x direction is twice that in the y direction, and so forth.
The xy aspect ratio of the cooling structure itself can be defined in a similar fashion as:
aC = (4.12)
For both these, the aspect ratios are defined as the relevant x dimension to the relevant y
dimension. The cross-sectional area, AD [m2], of the domain, which gives an indication of the
physical size of the region about an inserted cooling structure, can be defined as:
AD = (4.13)
The cross-sectional area of the cooling structure insert can be related to AD by the following
equation:
AC = AD = (4.14)
Here is the fraction of the heat-generating medium that is occupied by the cooling system. A
value of 0.1 represents a case where 10% of the total medium volume is occupied by the cooling
system.
It is shown by the below equations that the newly defined variables, aD, aC, AD, and are
sufficient to define the xy view of Figure 4.2.
= aD AD (4.15)
AD
= (4.16)
aD
= aC AD (4.17)
AD
= (4.18)
aC
4.4.2. Restrictions on the Cooling Structure Aspect Ratio
With the new geometric definitions and variables given in the previous sub-section, it is important
to consider the restrictions that exist on especially the cooling structure aspect ratio. When
referring to Figure 4.2, it may be noted that the defined cooling structure has to fit within the
representative domain for the current analysis to be meaningful. There thus exists a restricted
range of cooling structure aspect ratio values. The minimum and maximum cooling aspect ratios
can be expressed in terms of the representative domain aspect ratio and the fraction of the domain
that is used for cooling purposes:
aC ,min = a D (4.19)
aD
aC ,max = (4.20)
40 Chapter 4 – Two-Dimensional Analysis
Graphical representation of the minimum and maximum value aC cases are given in Figure 4.5.
Cooling structure
for = 0.25
y
Maximum aC Minimum aC
x case case
Figure 4.5 Maximum and minimum aspect ratio values
4.5. Analysis Approaches
Once a foundation has been laid down for the definition of geometric relations and restrictions, it
was possible to investigate different analysis methods with which the steady-state temperature
distribution of the representative domain can be found. One of the first options is that of the use
of either a two-dimensional or three-dimensional approximation of the temperature distribution.
Being the more simplified approach of the two, the two-dimensional approximation approach was
dealt with first.
A two-dimensional analysis is in actual fact a unit-depth three-dimensional analysis where the
influence of the third dimension on the variables being solved is ignored. In the following
sections different methods of performing two-dimensional analyses are explored. The two main
options for determining the two-dimensional temperature profile within the representative domain
discussed here are, analytical and numerical approaches.
Different methods of representing the cooling structure also exist including a generalised cooling
structure theory or a solid-state cooling structure theory. With a generalised cooling structure
theory it was attempted to describe the temperature distribution about a cooling structure
associated with either cooling channel or embedded solids configurations with the same set of
equations and relations.
4.6. Analytical Analysis – Generalised Cooling Structure
Obtaining an analytical solution for the temperature profile has the advantage that a less complex
procedure would be needed when attempting to optimise the distribution and geometric shape of
cooling structure inserts.
As will be demonstrated, obtaining an analytical solution for the representative domain is not
always possible, or is very difficult. For this reason it is necessary as a first approximation to
make certain assumptions to simplify the problem even further.
Consider the case where the type of the cooling structure is not specified. (Thus where it can
either be a solid or a channel filled with a fluid.) In general, the heat transfer from the heat-
generating medium to the cooling structure can then be described in terms of a heat transfer
coefficient, similarly defined as a convective heat transfer coefficient.
Chapter 4 – Two-Dimensional Analysis 41
In the hypothetical case where this average heat transfer coefficient on the interface or wall
between the cooling structure and the heat-generating region, h [W/m2K], and the core
temperature, TC, of the cooling structure is known, it would be possible to obtain the average
temperature at the interface or wall between the cooling structure and heat-generating medium,
Tif , by applying the principle of conservation of energy:
(
qtotal = h Aif Tif − TC ) (4.21)
Here, qtotal [W] refers to the total heat-generation rate within the representative domain with unit
depth while Aif [m2] refers to the area of contact between the cooling structure and the heat-
generation region.
If the interface temperature can be linked to the core temperature of the cooling structure, it
would be of use if the temperature profile within the heat-generating medium can be expressed in
terms of the interface temperature. An amended domain is shown in Figure 4.6 that incorporates
this idea. Note that the governing equation has been altered to be applicable to a two-dimensional
domain. The partial derivative in the z direction is equal to zero.
∂T
=0
∂y
∂T ∂ 2T ∂ 2T
=0 + +G =0
∂x ∂x 2 ∂y 2
T=T(x,y)
∂T
1 =0
T = Tif ∂x
T =T
Cooling
structure T = Tif
y
T =T
x
∂T
=0
∂y
Figure 4.6 Two dimensional domain of interest for an analytical approximation
For a two-dimensional space the governing equation can be rewritten as:
∂ 2T ∂ 2T
+ + GM = 0 (4.22)
∂x 2 ∂y 2
Here GM is defined as:
'''
qM
GM = (4.23)
kM
In this equation qM [W/m3] and kM [W/mK] refer to the uniform heat-generation density and
'''
thermal conductivity respectively of the heat-generating region.
42 Chapter 4 – Two-Dimensional Analysis
It may be noted that the differential equation is non-homogenous which makes it particularly
difficult to solve. It is possible however to transform it into a homogenous partial differential
equation by the applying the following transformation:
U ( x, y ) = T ( x, y ) + ψ ( y ) (4.24)
By considering the boundary conditions, it would be possible to determine the expression of
ψ ( y ) . Due to the discontinuity in the boundary conditions in the x and y directions at point 1 (see
Figure 4.6), it places the solution of the differential equation in this particular format outside the
scope of this study.
It is therefore opted to amend the analytical approach once more by simplifying the domain such
that there is no discontinuity on the boundaries (no obtuse internal angles). Refer to Figure 4.7
for the amended representative domain. It should be noted that an additional assumption is
introduced for this first order approximation, namely that the heat flux perpendicular to the
diagonal boundary is zero. This is only true for the case where the domain has exactly the same
dimensions in the x and y directions, which is not in general true for the domain under
investigation.
y ∂T
=0
x ∂y ( ;0)
(0;0)
∂ 2T ∂ 2T
∂T + +G=0 ∂T
=0 ∂x 2 ∂y 2 =0
∂x T=T(x,y) ∂n
ˆ
n
(0; ) T = Tint ( )
Figure 4.7 Amended domain for an analytical approximation
Once again the same transformation can be applied as before. After considering the new
boundary conditions the following expression is obtained:
2
(
ψ ( y ) = 1 GM y 2 − ( − )2 ) (4.25)
After applying the separation of variables method and substituting the four boundary conditions
of the simplified domain, the following equation can be obtained:
∞
−
Ci N ( − )sinh (Nx )cos(Ny ) + ( − )cosh N y+ sin (Ny ) = −( − )GM y (4.26)
i =1 −
1 − 2i
with N = π. (4.27)
2( − )
Here Ci is an integration constant.
Obtaining this integration constant has proved to be considerably difficult. The solution of this
constant unfortunately falls outside the scope of this study and thus so does obtaining an
analytical solution with this method for the temperature profile within the heat-generating
medium. An alternative analytical solution approach to find the steady-state temperature
distribution could also not be found and other means of predicting the steady-state temperature
distribution were thus needed.
Chapter 4 – Two-Dimensional Analysis 43
4.7. Numerical Analysis – Generalised Cooling Structure
Due to unfortunate complexity of obtaining an analytical solution for the temperature distribution,
the focus was shifted to numerical methods. Numerical methods require the discretisation of the
governing equation(s) over numerous control volumes or cells of the system or domain under
investigation. By solving the governing equation(s) for al the control volumes, an approximation
of the relevant variable(s) can be obtained.
Although many good and reliable numerical software packages are available, the approach here
involves building numerical solution algorithms from scratch. This may seem to be a case of
reinventing the wheel but due to various factors it was necessary to move away from the
commercially available packages for the current section of study as will be described shortly.
Normally, numerical simulation involves three separated stages namely pre-processing,
processing, and post-processing. During pre-processing, the governing equations, solution
methods, and aspect such as the domain and mesh layout is defined. During processing, the
problem that was defined in pre-processing is solved via various iterative, implicit or other
methods. Post-processing is the stage where the obtained solution is arranged in an
understandable format, whether it is the creation plots or other outputs.
Due to the large number of anticipated simulation cases that would be needed to optimise cooling
structure shapes and distribution involving various variables for a two-dimensional domain, the
traditional pre-processing stage would become very time consuming when using commercially
available software packages. For instance, for a single material, broad geometrical and thermal
condition case study, at least 10 different mesh cases, would have to be constructed, analysed,
post-processed, and results compared to each in order to approximate the optimum cooling
geometry for the particular set of input values. Due to the large number of input value
combinations of interest to the scope of power electronic modules, several million case studies
would be required, each of which would contain a number of simulation runs.
Not only would each test case require a large amount of computer disc space, the effort required
administrating the large number of input combinations required would be immense. A method
was thus required with which the procedure could be automated, speeding up the pre-processing
stages, which traditionally require the largest fraction of time. Different commercially available
numerical packages were considered for use, but due to the lack of sufficient automation offered
by these none where found suitable for the purposes of the study.
For this reason it was opted to create computer code that would automatically do the pre-
processing functions and would feed the problem information straight into a solution procedure or
algorithm. With such a tool it would thus be possible to pre-program multiple simulation cases in
advance.
In this study, the finite difference method proved to be very useful for finding approximations of
the two-dimensional temperature distribution of the representative domain. This method involves
the creation of a mesh of nodes across the domain and then solving the temperatures at these
nodes. Various solution methods such as iterative and implicit procedures exist.
44 Chapter 4 – Two-Dimensional Analysis
4.7.1. Meshing
Similar as before in the previous sub-section, the cooling structure is represented here in general
as either being a solid or a channel filled with fluid. The local heat transfer across the interface is
defined by the following equation:
(
qif ,local = h Alocal Tif ,local − TC ) (4.28)
Note that in this case there is no longer made use of the uniform ‘constant’ averaged temperature
values along the interface. Rather with this approach, the local interface temperature, Tif,local, is
given more freedom and is determined according to the local heat transfer rate at that particular
spot on the interface or wall. The cooling structure behaviour and characteristics are still
modelled via single (constant) averaged heat transfer coefficient.
The heat-generating medium, which in this case forms the numerical representative domain, is
divided into three uniformly spaced blocks with N nodes in each direction per block. One row or
column is overlapped at the border between two blocks. See Figure 4.8.
one
overlapping
node column
N nodes N nodes
N nodes 1 2
N nodes
one
overlapping
node row
Cooling
3
Structure N nodes
y
x
N nodes
Figure 4.8 Schematic representation of domain structure
The spacing between nodes in a particular block was uniform while the node spacing across the
entire domain might not have been strictly uniform, as each block could have had its own x and y
node spacing depending on the domain dimensions. A graphical representation of the nodal mesh
of the domain and thermal boundary conditions is given in Figure 4.9.
Chapter 4 – Two-Dimensional Analysis 45
Adiabatic 2N-1
boundary number
of nodes
N number
of nodes
2N-1
number
h of nodes
h
TC
N number
Adiabatic
of nodes
boundary
Figure 4.9 Two-dimensional grid used for numerical analysis of a generalised cooling structure
It should be noted that in this approach, the cooling structure temperature is represented by the
core temperature, TC [˚C], only, and that the temperature distribution of the cooling structure is
not being solved at this stage.
All boundaries are adiabatic except for at the interface between the heat-generating medium and
the cooling structure, which is defined by equation (4.28). The ratio between the equivalent heat
transfer coefficient on the interface and the thermal conductivity, γ [m-1] is defined as:
h
γ= (4.29)
kM
4.7.2. Discretised Governing Equation
The discretised governing equation for different parts of the domain are given in this section.
Refer to Figure 4.10 for the definition of the relative numbering scheme used for neighbouring
nodes. North, south, east, and west were used as reference directions as indicated by subscripts
N, S, E, and W respectively.
TN control
volume
yN
TW T TE
yS
TS
xW xE
Figure 4.10 Relative numbering scheme for neighbouring nodes
46 Chapter 4 – Two-Dimensional Analysis
By applying the principle of conservation of energy to a control volume about a generalised node
at a specific location in the domain, an equation can be obtained which expresses the nodal
temperature in terms of the temperatures at the neighbouring nodes.
For nodes not falling on the interfaces between the blocks shown in Figure 4.8 the mesh spacing
was uniform in both the x and y directions such that:
x E = xW = x (4.30)
y N = ys = y (4.31)
For such cases the general formats of the discretised governing equation are given below. Refer to
the small diagram next to each equation for the location type where the equation is valid.
GM x 2 y 2 + (TE + TW ) y 2 + 2TN x 2 + 2γx 2 yTC
T=
Cooling ( )
2 x 2 + y 2 + 2γx 2 y
Structure
Insu- G M x 2 y 2 + 2(TE y 2 + T N x 2 ) + 2γx 2 yTC
lated T=
Cooling
2(x 2 + y 2 ) + 2γx 2 y
Structure
G M x 2 y 2 + 2TE y 2 + (T N + TS )x 2
T=
2(x 2 + y 2 )
Insu-
lated
G M x 2 y 2 + 2(TE y 2 + TS x 2 )
Insulated T=
2(x 2 + y 2 )
G M x 2 y 2 + 2TS x 2 + (TE + TW ) y 2
Insulated
T=
2(x 2 + y 2 )
GM x 2 y 2 + 2(TW y 2 + TS x 2 )
Insulated
T=
2(x 2 + y 2 )
G M x 2 y 2 + 2TW y 2 + (T N + TS )x 2
T=
2( x 2 + y 2 )
Insu-
lated
G M x 2 y 2 + 2(TW y 2 + T N x 2 )
T=
Insulated
2(x 2 + y 2 )
G M x 2 y 2 + 2T N x 2 + (TE + TW ) y
T=
Insulated
2(x 2 + y 2 )
Cool- G M x 2 y 2 + 2(T N x 2 + TE y 2 ) + 2γxy 2 TC
T=
2(x 2 + y 2 ) + 2γxy 2
ing
Struc- Insu-
ture lated
Cool- G M x 2 y 2 + (T N + TS )x 2 + 2TE y 2 + 2γxy 2TC
T=
2(x 2 + y 2 ) + 2γxy 2
ing
Struc-
ture
Chapter 4 – Two-Dimensional Analysis 47
G M x 2 y 2 + (TE + TW ) y 2 + (T N + TS )x 2
T=
2(x 2 + y 2 )
If however a nodal point falls on the interface between two blocks indicated in Figure 4.8, the
grid spacing in the x and y directions are not necessarily the same. This required some adjustment
to the discretised equations listed above.
In such cases, the total length of the control volume in the x and y directions can be written
respectively as:
X = x E + xW (4.32)
Y = y N + yS (4.33)
The adjusted discretised equations are given below:
G M XxW x E y 2 + (T N + TS ) XxW x E + 2(TE xW + TW x E ) y 2
T=
E W
2 X (x E xW + y 2 )
Insulated G M XxW x E y 2 + 2TS XxW x E + 2(TE xW + TW x E ) y 2
T=
2 X (x E xW + y 2 )
E W
N G M x 2Yy N y S + (TE + TW )Yy N y S + 2(T N y S + TS y N )x 2
T=
2Y (x 2 + y N y S )
S
N G M x 2Yy N y S + 2TW Yy N y S + 2(T N y S + TS y N )x 2
T=
Insu-
lated 2Y (x 2 + y N y S )
S
A special case is found at the corner of the cooling structure where both the x and y directional
grid spacing might not be uniform:
xW xE
yN
Cooling
Struc- yS
ture
T=
(
G M ( Xy N + x E y S )x E xW y N y S + 2 TE xW Yy N y S + TW x E y N y S + T N Xx E xW y S + TS x E xW y N + γ ( x E + y S )x E xW y N y S TC
2 2
)
(
2 xW Yy N y S + x E y N y S + Xx E xW y S + x E xW y N + γ ( x E + y S )x E xW y N y S
2 2
)
4.7.3. Tri-Diagonal-Matrix-Algorithm (TDMA) Solution Method
In this first numerical approach the temperature values for the nodes within the heat-generating
medium were solved iteratively. Various iterative methods for obtaining the temperature
distribution for a mesh exist. Normally an initial estimate is applied to all nodes in the domain.
In this case all initial nodal temperature values were set equal to TC (the cooling structure central
temperature).
48 Chapter 4 – Two-Dimensional Analysis
After each iteration run, the temperature values at each node are closer to the ultimate converged
values. Once the difference in temperature values from one iteration run to the next, is
satisfactory small, convergence is assumed.
A line-by-line method combining the Gauss-Seidel and the Tri-Diagonal-Matrix-Algorithm
(TDMA) methods was used. This method has the advantage that thermal boundary information
is simultaneously transferred across an entire line in the domain, something that would have taken
several iterations if a pure Gauss-Seidel approach were followed [1]. The TDMA method is
described shortly below.
Consider a line of nodes in the domain either being a row or a column running from one border to
the other shown in Figure 4.11:
T1 Ti Ti+1 TN
Figure 4.11 Row of node on which the TDMA method can be applied on
From the discretisation equations it is possible to express the temperature at on one node on the
line in terms of its neighbours on the same line as:
aTDMA,i Ti = bTDMA,iTi +1 + cTDMA,i Ti −1 + d TDMA,i (4.34)
Here aTDMA through to dTDMA are coefficients aiding in this process and are dependent on
geometric, material properties and thermal conditions.
From this, equation (4.34) can be rewritten as:
Ti = PTDMA,i Ti +1 + QTDMA,i (4.35)
With
bTDMA,i
PTDMA,i = (4.36)
aTDMA,i − cTDMA,i PTDMA,i −1
d TDMA,i + cTDMA,i QTDMA,i −1
QTDMA,i = (4.37)
aTDMA,i − cTDMA,i PTDMA,i −1
It can be deduced from the fact that, cTDMA,1 = 1, that
bTDMA,1
PTDMA,1 = (4.38)
aTDMA,1
and that
d TDMA,1
QTDMA,1 = (4.39)
aTDMA,1
Also, because bTDMA,N = 0, it can be found that
TN = Q N (4.40)
Chapter 4 – Two-Dimensional Analysis 49
In order to solve the temperature of the nodes on a particular line in one iteration, the following
steps should be taken:
1. Determine PTDMA,1 and QTDMA,1 with equations (4.38) and (4.39)
2. Use equations (4.36) and (4.37) to calculate PTDMA,i and QTDMA,i for i = 2,3,…N
3. Use equation (4.40) to determine TN
4. Use equation (4.35) to obtain Ti for i = N-1 to 1
With use of the line-by-line method in conjunction with Gauss-Seidel method, once one line of
nodes is solved, a neighbouring line can be solved with a similar method until the entire domain
has been covered. The influence of neighbouring nodes which do not fall on the line under
consideration is incorporated in dTDMA,i. This whole process of sweeping across the domain is
repeated until the convergence criterion has been satisfied.
In a two-dimensional space the progression to convergence is more rapid if a column-by-column
sweep is repeated by a row-by-row sweep.
4.7.4. Numerical Results
As a first step, a set of Visual Basic macros was added into Microsoft Excel 2000 spreadsheets to
do all computational work. A typical two-dimensional temperature distribution in the heat-
generating medium is shown in Figure 4.12. The region of the cooling structure, since it not is
modelled, has a uniform temperature of TC.
aD = 1
aC = 1
= 0.5
2
G M = 1.0 K/m
2
AD = 1 m
1.05
-1
=1m
1.00 o
TC = 0 C
0.95
Temperature [o C]
0.90
Boundary of the
0.85 heat generating
medium with the
0.80
cooling structure
This area of the
0.75 domain is at a
constant
0.70 temperature of
T C = 0 oC
Figure 4.12 A typical temperature distribution in the representative domain1
For the thermal boundary of a uniform average heat transfer coefficient on the interface of the
cooling structure and the heat-generating medium, it was found that the maximum temperature,
Tmax, in the domain is always in the corner furthest away from the cooling structure.
1
Refer to Microsoft Excel 2000 spreadsheet “2D-2ndApproach-BasicMacros.xls” for data and details of
the visual basic macro used to perform all calculations.
50 Chapter 4 – Two-Dimensional Analysis
Since the location of the maximum or peak temperature was known, it is possible to determine
what influence the aspect ratio of the cooling structure, aC, has on the maximum temperature for
different domain aspect ratios, aD. See Figure 4.13.
While changing the value of aC, the percentage of the original heat-generating medium volume
occupied by the cooling structure, , is kept constant; in this case 50% ( = 0.5) of the medium is
used for cooling purposes. In Figure 4.13 the horizontal axis is defined in such a way that results
for different aD values could be placed on one graph.
2.50 = 0.5
2
G M = 1 K/m
Maxim um Tem perature (Tmax ) [ 0C]
-1 a D values:
2.00 =1m
0 1
TC = 0 C 1.25
1.50 1.5
2
5
10
1.00
0.50
0.00
a C = a C,min aC=aD a C = a C,max
Cooling Structure Aspect Ratio (a C )
Figure 4.13 Influence of aC on the maximum temperature in the domain2
For aD > 1, it was found that the aspect ratio of the cooling structure corresponding to the lowest
peak temperature in the domain always fell between aC = aD and aC = aC,max. This was tested for
various conditions. From this, it is clear that a thermal optimum cooling structure aspect ratio
exist for a particular set of geometric and thermal conditions such that the peak temperature are
its lowest possible level.
As it is known in which interval of aC an optimum geometric shape for the cooling structure can
be found it is possible to define a new variable which represents the aspect ratio shape of the
cooling structure:
ac − aD
aC ,rel = (4.41)
aC ,max − aD
This new definition normalises aC. If the relative cooling structure aspect ratio, aC,rel, has a value
of 1, it represents the case where the cooling structure has the highest possible aspect ratio, while
when the aC,rel = 0 it represents the case where aC is equal to aD.
4.7.5. Optimisation Results
With the new variable which represents the relative aspect ratio of the cooling structure in terms
of its maximum allowed value and the aspect ratio of the domain, it is possible to optimise aC for
2
See Microsoft Excel 2000 spreadsheet “2D-2_ac_InfluenceOnTmax.xls”
Chapter 4 – Two-Dimensional Analysis 51
various conditions and compare the results. Some optimisation results are given in the following
sub sections.
Only aD values greater than 1 was investigated as mirror behaviour for aD less than one are
evident. Thus, aD = 2 would have the same maximum temperature as aD = 0.5 as these two case
are mirror images of each other about the line x = y. The same principle is valid for all inverse
pairs of aD.
It was found that 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 ∝ G M (4.42)
This means that neither GM nor TC has an influence on the shape of the optimum cooling
geometry.
4.7.5.1. Domain with a Square Cross-Sectional area (aD = 1)
An investigation3 was conducted on a domain case with a square cross-sectional area. This is
applicable to conditions where the cooling structure spacing in both Cartesian directions is
identical.
The optimum relative aspect ratio for an entire domain area of 1 m2 is shown in Figure 4.14 for
various values. For low values of , the optimum aC is close to its maximum value, while as
is increased (or a greater fraction of the domain is used for cooling) the optimum aC drops to be
equal to aD. Thus, the cooling structure and the representative domain has the same aspect ratio.
It was found that by increasing (thus increasing the convective heat transfer coefficient between
the heat-generating medium and the cooling structure), aC,rel,optimum decreases.
1
= 0.1
0.9 = 0.5
0.8 =1
=2 -1
0.7 = 5 (unit of : m )
= 10
aC,rel,optimum
0.6
0.5
= 20
= 50
0.4
0.3 aD =1
0.2 A D = 1 m2
0.1
0
0 0.2 0.4 0.6 0.8 1
Fraction of Domian used for Cooling ( )
Figure 4.14 Optimum relative cooling structure aspect ratios for a square domain of area of
1m2 with various values of and .
3
See Microsoft Excel 2000 spreadsheet “2D-2_Optimum_aCrel(1).xls”
52 Chapter 4 – Two-Dimensional Analysis
In Figure 4.15 to Figure 4.17 the influence, which different domain areas have on the optimum
aspect ratio of the cooling structure, is shown for values of 1 m-1, 10 m-1, and 50 m-1
respectively. It can be seen that as the domain area is decreased, aC,rel,optimum increases. A smaller
domain area represents a case where more cooling structures of smaller sizes are used without
increasing the fraction used for cooling.
1
= 0.1
0.8
aC,rel,optimum
= 0.2
0.6
= 0.3
0.4
aD =1 = 0.4
-1 0.2
=1m
0
1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 1.E+02
2
Area of Entire Domain (A D ) [m ]
Figure 4.15 Influence of the entire domain area for = 1m-1
1 aD =1
-1
0.8 = 10 m
aC,rel,optimum
0.6
= 0.1 0.4
= 0.2
= 0.3 0.2
= 0.4
0
1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 1.E+02
2
Area of Entire Domain (A D ) [m ]
Figure 4.16 Influence of the entire domain area for = 10m-1
1 aD =1
-1
0.8 = 50 m
aC,rel,optimum
0.6
= 0.1
0.4 = 0.2
= 0.3
= 0.4 0.2
0
1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 1.E+02
2
Area of Entire Domain (A D ) [m ]
Figure 4.17 Influence of the entire domain area for = 50m-1
4.7.5.2. Domain with 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. This
case was investigated4 and some results are given below.
4
See Microsoft Excel 2000 spreadsheet “2D-2_Optimum_aCrel(2).xls”
Chapter 4 – Two-Dimensional Analysis 53
Figure 4.18 shows the behaviour of aC,rel,optimum for various values of for a domain area, AD of
1 m2. The influence of the domain area on the optimum aspect ratio is shown in Figure 4.19 to
Figure 4.21. Similar trends are observed as with a case where aD = 1 except that the optimum
aspect ratio is never the same as the domain aspect ratio and that in general higher optimum
relative cooling structure aspect ratios are obtained.
1 (unit of
-1
:m ) = 0.1
0.95 = 0.5
=1
0.9 =2
0.85 =5
= 10
aC,rel,optimum
0.8
= 20
0.75 = 50
0.7
0.65
0.6 aD = 2
0.55 A D = 1 m2
0.5
0 0.2 0.4 0.6 0.8 1
Fraction of Domain used for Cooling ( )
Figure 4.18 Optimum relative cooling structure aspect ratios for a domain with aD = 2 with an
area of 1 m2
1
= 0.2
0.9
aC,rel,optimum
0.8 = 0.4
aD = 2 0.7 = 0.6
=1m -1 = 0.8
0.6
1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 1.E+02
Area of Entire Domain (A D ) [m2]
Figure 4.19 Influence of the entire domain area for = 1 m-1
1 aD = 2
0.9 = 10 m-1
aC,rel,optimum
= 0.2
0.8
0.7 = 0.4
0.6 = 0.8
= 0.6
0.5
1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 1.E+02
2
Area of Entire Domain (A D ) [m ]
Figure 4.20 Influence of the entire domain area for = 10 m-1
54 Chapter 4 – Two-Dimensional Analysis
1 aD = 2
0.9 = 50 m-1
aC,rel,optimum
0.8 = 0.2
0.7 = 0.4
0.6 = 0.6
0.5
= 0.8
1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 1.E+02
2
Area of Entire Domain (A D ) [m ]
Figure 4.21 Influence of the entire domain area for = 50 m-1
4.7.5.3. Domain with aD = 5
The optimum aC values5 for aD = 5 is given in Figure 4.22. As before, an increase in , has as a
result a decrease in aC,rel,optimum. The influence of the entire domain area magnitude is given in
Figure 4.23 to Figure 4.25. Similar trends are observed as with aD = 1 and aD = 2.
1
0.98
0.96
0.94 Series1
= 0.1
aC,rel,optimum
0.92 Series2
= 0.5
0.9 Series3
=1
1
0.88 Series4 (unit of : m- )
=2
Series5
=5
0.86
Series6 aD = 5
= 10
0.84 A D = 1 m2
Series7
= 20
0.82 Series8
= 50
0.8
0 0.2 0.4 0.6 0.8 1
Fraction of Domain used for Cooling ( )
Figure 4.22 Optimum relative cooling structure aspect ratios for a domain with aD=5 with an
area of 1 m2
5
See Microsoft Excel 2000 spreadsheet “2D-2_Optimum_aCrel(5).xls”
Chapter 4 – Two-Dimensional Analysis 55
1
0.98 = 0.2
aC,rel,optimum
0.96 = 0.4
0.94 = 0.6
aD = 5
0.92 = 0.8
= 1 m-1
0.9
1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 1.E+02
2
Area of Entire Domain (A D ) [m ]
Figure 4.23 Influence of the entire domain area for = 1 m-1
1 aD = 5
0.98 = 10 m-1
aC,rel,optimum
= 0.2
0.96
= 0.4
0.94
= 0.6
0.92
= 0.8
0.9
1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 1.E+02
2
Area of Entire Domain (A D ) [m ]
Figure 4.24 Influence of the entire domain area for = 10 m-1
1
= 0.2
aC,rel,optimum
0.95
= 0.4
0.9 = 0.6
aD = 5
= 50 m-1 = 0.8
0.85
1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 1.E+02
2
Area of Entire Domain (A D ) [m ]
Figure 4.25 Influence of the entire domain area for = 50 m-1
4.7.5.4. Domain with aD = 10
Figure 4.26 inicates the value of the optimum aspect ratios6 of the cooling structure for a domain
where the overall representative domain aspect ratio is 10. It may be noted that there is
significantly less difference in aC,rel,optimum for different values than before. The influence of the
overall domain area, AD, as shown in Figure 4.27 is also not as great as for lower aD values, and
aC,rel,optimum fell into a smaller range.
6
See Microsoft Excel 2000 spreadsheet “2D-2_Optimum_aCrel(10)updated.xls”
56 Chapter 4 – Two-Dimensional Analysis
1
1
(unit of : m- )
0.99
0.98
= 0.1
Series1
aC,rel,optimum
0.97 = 0.5
Series2
0.96 =1
Series3
=2
Series4
0.95 =5
Series5
0.94 =10
Series6 a D = 10
= 20
Series7 2
0.93 AD = 1m
Series8
=50
0.92
0 0.2 0.4 0.6 0.8 1
Fraction of Domain used for Cooling ( )
Figure 4.26 Optimum relative cooling structure aspect ratios for a domain with aD = 10 with
an area of 1 m2
1
= 0.2
0.99
aC,rel,optimum
0.98 = 0.4
0.97 = 0.6
a D = 10
0.96 = 0.8
= 1 m-
0.95
1.E-03 1.E-02 1.E-01 1.E+00 1.E+01 1.E+02
2
Area of Entire Domain (A D ) [m ]
Figure 4.27 Influence of the entire domain area for = 1 m-1
4.7.5.5. Behaviour of the Optimum Cooling Structure Aspect Ratio
For all the cases shown above it may be seen that an increase in the fraction of the entire domain
used for cooling, there is a decrease in the relative optimum aspect ratio of the cooling structure.
Also, from the above graphs it is evident that with an increase in aD, the optimum aspect ration of
the cooling structure tends to the maximum allowed cooling structure aspect ratio even for a wide
range in values.
This trend is also true for a decrease in the domain area. A decrease in the domain area represents
the case where the heat-generating medium is cooled by a greater number of cooling structures -
each with smaller cross sectional areas.
When the aspect ratio of the cooling structure is at its maximum allowed value (aC,rel,optimum = 1),
this results in neighbouring cooling structures along a particular Cartesian direction to be merged
into one continues strip or layer as shown in Figure 4.28.
Chapter 4 – Two-Dimensional Analysis 57
Heat generating
medium
y
Merged cooling
structures
x
Figure 4.28 Physical meaning of a case with the maximum allowed aC
4.7.6. Maximum Temperatures for Optimum Geometry
With the known optimum aspect ratios of the cooling structure for a wide variety of cases, the
behaviour of the maximum temperature in terms of GM, TC, , , aD and AD for these optimum
geometries were investigated7.
The maximum temperature in the domain for the optimum geometric shapes obtained in the
previous sections can now be plotted taking the influence of TC and GM into account according to
the relation given in equation (4.42).
Figure 4.29 shows that there is a significant drop in the maximum temperature as aD and is
increased. It would thus be more advantageous if the centre-to-centre distance between adjacent
cooling structures in a particular Cartesian direction were much greater than the distance between
adjacent cooling structures in the other Cartesian direction. This was found to be the optimum
geometry as aD is increased.
10.00
aD = 1
aD = 2
(Tmax -TC )/GM [m 2 ]
1.00 aD = 5
aD = 10
0.10
A D = 1 m2
= 0.4
0.01
0.1 1 10 100
[m -1 ]
Figure 4.29 Influence of aD and on the maximum temperature
In Figure 4.30 it is shown that a decrease in AD and an increase in will lower the maximum
temperature for a particular value of GM. These trends were found to be true for all cases of , ,
aD and AD.
7
See Microsoft Excel 2000 spreadsheet “2D-2_Temp BehaviouratOptimum.xls”
58 Chapter 4 – Two-Dimensional Analysis
10
aD = 2
(Tmax -TC )/GM [m 2 ]
-1
1 =1m
0.1
= 0.2
= 0.4
0.01 = 0.6
= 0.8
0.001
0.0001 0.001 0.01 0.1 1 10
Area of Entire Domain (A D ) [m2]
Figure 4.30 Influence of AD and on the maximum temperature
Unexpectedly it is seen that a very exponentially based linear relation exists between the
maximum temperature in the domain and the domain area. This can in future be useful to draw
up generalised graphs with which the maximum temperature in the domain can be estimated.
4.7.7. Summary of Results up to This Point
It was found that for this numerical scheme the optimum cooling structure aspect ratio is
dependent on physical geometric conditions, material properties, as well as the heat transfer
coefficient between the heat-generating medium and the cooling structure.
Increases in the aspect ratio of the domain, the fraction of the domain used for cooling, and the
ratio between the heat transfer coefficient and the thermal conductivity results in a decreased
maximum temperature.
From a pure thermal point of view the results indicate that continuous cooling strips or layers are
the most efficient in cooling a heat-generating medium where a relatively small volume fraction
is occupied by cooling inserts. Whether these optimum geometries can be applied in practice will
depend on physical, material, and manufacturing restrictions.
Even though this method was the first tangible step towards characterising the temperature profile
in a corner quadrant domain, it has many shortfalls. The use of an average heat transfer
coefficient links easily with conventional convective heat transfer theory for the case where the
cooling structure is a cooling channel, but the same might not be said for the case of a solid
cooling structure. A great amount of investigation would be needed to link a solid structure’s
thermal behaviour with a single heat transfer coefficient.
Also, for the case of a fluid filled cooling channel, it is not always accurate to model the
convective heat transfer coefficient as a constant along the whole perimeter of the channel. In
practice it is known that for rectangular cross-sectioned channels the convective heat transfer
coefficient in the corners are much less than in middle of the sides.
A more accurate numerical approach than the one used up to this point was needed to
approximate the two-dimensional temperature distribution in the representative corner quadrant
domain.
Chapter 4 – Two-Dimensional Analysis 59
4.8. Numerical Analysis – Solid Cooling Structure
It was found in the previous section that the use of a generalised approach for the cooling
structure is not very representative of a solid cooling insert. In this sub-section focus is placed on
solid cooling inserts and the numerical approach modified accordingly.
As a first attempt in embedding a cooling system into integrated power electronic modules, the
use of solid state cooling structures has a great manufacturing and operating advantage above the
use of fluid cooling schemes as described in Chapter 2. With this in mind, a separate numerical
model for determining the temperature in and around a solid state cooling structure was needed.
When considering heat conduction in a homogeneous solid in a one-dimensional space, it is easy
to describe the linear temperature profile across a length, x, by the following elementary
equation:
∆T q
= (4.43)
∆x kA
Unlike the one-dimensional case, it becomes significantly more complicated to describe the
temperature profile in a two-dimensional space as it involves the solution of the following second
order partial differential equation:
∂ 2T ∂ 2T
+ =0 (4.44)
∂x 2 ∂y 2
The solution of this equation is heavily dependent on conductor shape and all boundary
conditions, both of which is not readily available in the current corner quadrant domain type
when considering the interface between the cooling structure insert and the heat–generating
medium.
It was thus evident that the use of the first mentioned generalised numerical approach would be
difficult to implement. Rather than attempting to link the separate temperature profile solutions
in the cooling structure and heat-generating medium, it might be advantageous to find a solution
of the entire representative corner quadrant domain that includes both these regions.
4.8.1. Meshing and Numerical Solution Method
A similar numerical meshing scheme was used as before but with the solid-state cooling structure
or heat-extractor, the mesh was extending to include the cooling structure region as well. See
Figure 4.31 for more information. It should be noted that there were no nodes on the interface of
the heat-extraction and the heat-generating mediums. This was done due to mathematical
discontinuities that arise when the thermal contact resistance on the interface is zero or very
small. In the computer code, used for the current analysis, an option was given for the use of
either a piecewise uniform grid spacing or a non-uniform grid spacing.
60 Chapter 4 – Two-Dimensional Analysis
(0 ; ) ( ; )
( ; )
interface with
thermal contact ( ; )
resistance
TC
y ( ;0)
cooling structure with heat-generating
x
uniform heat flux out medium
Figure 4.31 Full two-dimensional nodal mesh across the entire domain
One of the greatest differences between the current meshing scheme and the previous one is that
with current one a ‘heat flux boundary’ is used on one of the cooling structure region’s xy faces to
simulate the axial conduction of heat from the domain into the z direction. The heat flux is
defined for a unit depth of the domain as:
qM (
' ''
− )
''
qC = − (4.45)
This is obtained by applying the principle of conservation of energy to the entire domain. The
negative sign indicates the loss of heat to the surroundings. Defining a heat flux in the z direction
for a two-dimensional xy analysis may seem contradictory, but it is easy to incorporate into the
governing equation used:
∂ 2T ∂ 2T
+ + GC = 0 (4.46)
∂x 2 ∂y 2
with
''
qC
GC = (4.47)
kC
Here kC represents the thermal conductivity of the solid cooling structure.
This is similar to the governing equation used for the heat-generating medium except that the G
term is now negative instead of positive indicating that heat is extracted at that location.
Another difference is the inclusion of internal thermal interface resistance, Rint [m2K/W] on the
planes where the heat-generating medium and the solid-state cooling structure meet. Due to the
presence of two solids being present, the ration is also adapted and redefined as being:
kC
γ= (4.48)
kM
Unlike before, it is now dimensionless.
The core cooling temperature, TC, was not solved for but rather fixed and used as a reference
temperature for the rest of the domain nodes.
As before, the Tri-Diagonal Matrix Algorithm (TDMA) combined with the Gauss-Seidel method
was used to speed up convergence. Various improvements were made to the solution method to
increase the time efficiency of the process.
Chapter 4 – Two-Dimensional Analysis 61
4.8.2. Discretised Governing Equations
For the TDMA method, it is convenient to express the temperatures of neighbouring nodes in
terms of each other by means of a coefficient equation:
CT T = +C S TS + CW TW + C E T E + C N T N + C G (4.49)
Subscripts S, W, E, and N refer to the four reference directions namely “south”, “west”, “east”
and “north” respectively as shown in Figure 4.10. Subscript T refer to temperature node under
consideration and G indicates the heat gain or heat loss of the particular node. The coefficients
for this equation for different node locations are given below.
4.8.2.1. Discretised Coefficients for Heat-Generating Medium Nodes
Adiabatic boundary to the south CT: 2 X [xE xW + yS 2 ]
CS: 0
CE: 2 xW yS 2
CW: 2 xE y N 2
CN: 0
CT: 2 X [xE xW + y N 2 ] CG: GM XxE xW yS 2
CE: 2 xW y N 2 Adiabatic boundaries to the north and the
CN: 2 XxE xW west
CG: GM XxE xW y N 2 CS: 2 xE 2
Adiabatic boundaries to the south and east CW: 0
CS: 0 CT: 2[xE 2 + y S 2 ]
CW: 2 y N 2 CE: 2 yS 2
CT: 2[xW 2 + y N 2 ] CN: 0
CE: 0 CG: GM xE 2 yS 2
CN: 2 xW 2 Adiabatic boundary to the west
CG: GM xW 2 y N 2 CS: 2 xE 2 y N
Adiabatic boundary to the east CW: 0
CS: 2 xW 2 y N CT: 2Y [y N yS + xE 2 ]
CW: 2YyN yS CE: 2YyN yS
CT: 2Y [y N y S + xW 2 ] CN: 2 xE 2 y S
CE: 0 CG: GM xE 2YyN yS
CN: 2 xW 2 yS Interface to the south and adiabatic
CG: GM xW 2YyN yS boundary to the west
Adiabatic boundaries to the north and east CS: 4 xE 2ζ
CS: 2 xW 2 CW: 0
CT: 2[YyN + xE 2 (1 + 2ζ )]
CW: 2 yS 2
CE: 2Yy N
CT: 2[xW 2 + y S 2 ]
CN: 2 xE 2
CE: 0
CN: 0 CG: GM xE 2YyN
CG: G M xW 2 y S 2 with ζ =
yN
Adiabatic boundary to the north yS (1 + γ ) + 2k M Rint
−1
CS: 2 XxE xW
CW: 2 xE yS 2
62 Chapter 4 – Two-Dimensional Analysis
Interface to the south Interface to the west and adiabatic boundary
CS: 4 XxE xW ζ to the south
CW: 2 xEYyN CS: 0
CT: 2 X [YyN + xE xW (1 + 2ζ )] CW: 4 y N 2ζ
CE: 2 xW YyN CT: 2[XxE + y N 2 (1 + 2ζ )]
CN: 2 XxE xW CE: 2 y N 2
CG: GM YyN z 2 XxE xW YyN CN: 2 XxE
yN CG: GM Xx E y N 2
with ζ =
yS (1 + γ −1 ) + 2k M Rint xE
with ζ =
Interface to the west xW (1 + γ ) + 2k M Rint
−1
CS: 2 XxE y N Central node
CW: 4YyN yS ζ CS: 2 XxE xW y N
CT: 2Y [ XxE + y N y S (1 + 2ζ )] CW: 2 xW YyN y S
CE: 2YyN yS CT: 2 XY [xE xW + y N yS ]
CN: 2 XxE yS CE: 2 xEYyN yS
CG: GM XxEYy N yS CN: 2 XxE xW y S
xE CG: GM XxE xW YyN y S
with ζ =
xW (1 + γ ) + 2k M Rint
−1
4.8.2.2. Discretised Coefficients for Cooling Structure Nodes
Adiabatic boundary to the south Interface to the east
CS: 0 CS: XxW y N
CW: xE y N 2 CW: YyN y S
CT: X [x E xW + y N 2 ] CT: Y [ XxW + y N y S (1+ 2ζ )]
CE: xW y N 2 CE: 2YyN yS ζ
CN: XxE xW CN: XxW y S
CG: GC XxE xW y N 2
' CG: GC XxW YyN yS
'
Interface to the east and adiabatic boundary xW
With ζ =
to the south xE (1 + γ ) + 2kC Rint
CS: 0 Interface to the north and east
CW: y N 2 CS: XxW
CT: XxW + y N 2 (1 + 2ζ ) CW: YyS
CE: 2 y N 2ζ CT: z 2 [YyS (1 + 2ζ x ) + XxW (1 + 2ζ y )]
CN: XxW CE: 2YySζ x
CG: G XxW y N
'
C
2
CN: 2 XxW ζ y
xW CG: GC XxW YyS
'
With ζ =
xE (1 + γ ) + 2kC Rint xW
With ζ x =
xE (1 + γ ) + 2kC Rint
yS
And ζ y =
y N (1 + γ ) + 2kC Rint
Chapter 4 – Two-Dimensional Analysis 63
Interface to the north Adiabatic boundary to the west
CS: XxE xW CS: xE 2 y N
CW: xEYyS CW: 0
CT: X [YyS + xE xW (1+ 2ζ )] CT: Y [y N y S + xE 2 ]
CE: xW YyS CE: YyN y S
CN: 2 XxE xW ζ CN: xE 2 yS
CG: GC XxE xW YyS
'
CG: GC x E 2YyN y S
'
yS Central node
With ζ =
y N (1 + γ ) + 2kC Rint CS: XxE xW y N
Adiabatic boundary to the west with CW: xEYyN yS
interface to the north CT: XY [ y N yS + xE xW ]
CS: x E 2 CE: xW YyN yS
CW: 0 CN: XxE xW y S
CT: YyS + xE 2 (1 + 2ζ )
CG: GC XxE xW YyN yS
'
CE: YyS
CN: 2 xE 2ζ
CS: GC xE 2YyS
'
yS
With ζ =
y N (1 + γ ) + 2kC Rint
4.8.3. Numerical Results
A Visual Basic Program was written into Microsoft Excel 2000 spreadsheets and tested. Figure
4.32 shows the converged steady state temperature solution obtained for a case where = 1 m,
= 0.9 m, = 1 m, = 0.5 m, kM = 10 W/mK, kC = 100 W/mK, and Rint = 0.2 m2K/W for a 40 by
40 node mesh8.
8
See Microsoft Excel 2000 spreadsheet “BasicMacros-2D-3rdApproach.xls” for data and the Visual Basic
Program
64 Chapter 4 – Two-Dimensional Analysis
=1m
= 0.9 m
=1m
1.032 = 0.5 m
k M = 10 W/mK
0.956
k C = 100 W/mK
0.88
R in t = 0.2 m 2K/W
0.804 Mesh: 40 by 40
0.728
Temp. Scale:
0.652
0.956-1.032
0.576 0.88-0.956
0.5 0.804-0.88
0.424 0.728-0.804
0.652-0.728
0.348
0.576-0.652
0.272
0.5-0.576
0.196 0.424-0.5
0.12 0.348-0.424
0.044 0.272-0.348
y -0.032
0.196-0.272
0.12-0.196
x 0.044-0.12
-0.032-0.044
Figure 4.32 Steady state temperature solution for an arbitrary case
An interesting observation, which could be made from this solution example, is that the
maximum temperature is not always located at the point furthest away from the cooling structure.
In this case the peak temperature is on the left hand rear end of the figure instead of in the rear
right hand corner.
4.8.4. Verification of Results
The solutions obtained from the program were tested by amongst other methods, comparing
solutions of transposed geometries to each other. If a geometry is rotated by 90˚, (represented by
switching the and , and and values around), the same temperature solution were obtained
for both cases. Other checks included the cases where either or . For such cases the
temperature solution can be approximated by a one-dimensional equation. Obtained numerical
solutions approached the expected one-dimensional results.
The obtained solution shown in Figure 4.32 was also verified by a commercially available
numerical solution package, STAR CD. This package is mainly aimed at computational fluid
dynamics simulation, but also has thermal solution capabilities.
The obtained solutions obtained from the Visual Basic code written and STAR CD compared
well with each other as shown in Figure 4.33.
Chapter 4 – Two-Dimensional Analysis 65
a) Solution obtained from own code b) Solution obtained from STAR CD
Figure 4.33 Comparison of solutions obtained with different methods
The commands given to STAR CD is listed in Table 4.1. It should be noted that the interfacial
resistance defined in STAR CD is only half the required value. This is due to the fact that the
package makes use of baffle cells to include interfacial resistances and that the user defined
thermal resistance is applied to both sides of these cells. The effective thermal resistance is thus
twice the specified value.
66 Chapter 4 – Two-Dimensional Analysis
Table 4.1 STAR CD Command line entries
! Domain with heat flux in, heat flux out, corner ctype 4
cellfaces = 0Celcius, Baffles with contact resistance vc3d 0 0.025 1 0 0.025 1 0 1 1
itle
TestHF5: A=1,a=0.9,B=1,b=0.5,Q=10,kM=10,kC=100,R=0.1 ! Create Flux-In Wall Boundary on Material 1
!clearing the grid structure cset news type 1
blkdel all view 0 0 1
cdel all $ ccom all $y term,,,rast $plty,ehid $replot
bdel all $ bcom all $y cplo
spldel all zoom,off $replot
vdel all bzone 3 all
!Switch on Conjugate Heat Transfer in order to use solid cset news type 2
cells view 0 0 1
conj on term,,,rast $plty,ehid $replot
cplo
!Define cell types zoom,off $replot
ctab 1 Solid 4 0 1 bzone 3 all
ctna 1 Solid_H
ctab 2 Solid 5 0 1 rdef,3,wall,standard
ctna 2 Solid_H nosl,stand,9,
ctab 3 solid 6 0 2 0,0,0,1,0
ctna 3 Solid_C flux,10
ctab 4 solid 7 0 2 rname,3,Heat_Flux_in
ctna 4 Solid_C_FixedTemp
ctab 5 baffle 8 0 1 ! Create Fixed Temp Wall BoundaRIES on one cell of Material 2
ctna 5 Interface
cset news type 4
ctab 6 baffle 9 0 1
view 0 -1 0
ctna 6 Interface
cplo
zoom, off $replot
!Gen Hotcells - Celltypes 1 and 2
bzone 4 all
ctype 1
vc3d 0 0.9 36 0.5 1 20 0 1 1
cset news type 4
view -1 0 0
!Create Baffle cells on the one side (xz face)
cset all cplo
view 0 -1 0 zoom, off $replot
term,,,rast $plty,ehid $replot bzone 4 all
cplo
zoom,off $replot rdef,4,wall,standard
czone baff 5 all nosl,stand,9,
0,0,0,0,0,
! Create Baffle boundary with resistance of 0.1 m2K/W fixed 0
cset news type 5 rnam 4 Fixed_Temp
cplo
bzone 1 all !Create Flux-Out Wall Boundary on Material 2
cset news type 3
rdef,1,baffle,standard view 0 0 1
nosl,stand,9, term,,,rast $plty,ehid $replot
0,0,0,1,0 cplo
1.e+30,1.e+30,0 zoom,off $replot
cond,0.1 bzone 5 all
SAME
cset news type 4
ctype 2 view 0 0 1
vc3d 0.9 1 4 0 0.5 20 0 1 1 term,,,rast $plty,ehid $replot
! Create baffle cells on the on side (yz face) cplo
cset news type 2 zoom,off $replot
view -1 0 0 bzone 5 all
cplo
zoom,off $replot rdef,5,wall,standard
czone baff 6 all nosl,stand,9,
0,0,0,1,0
! Create Baffle boundary with resistance of 0.1 m2K/W flux,-12.222222
cset news type 6 rname,5,Heat_Flux_OUT
cplo
bzone 2 all vmerg all
n
rdef,2,baffle,standard
nosl,stand,9, pmat 2 solid Solid_C
0,0,0,1,0 dens cons 100
1.e+30,1.e+30,0 spec cons 1000
cond,0.1 cond cons 100
SAME moni 1500
ctype 2 !solve only for temp:
vc3d 0.9 1 4 0.5 1 20 0 1 1 solve,n,n,n,n,n,n,y,n,n,n,n,n,
pmat 1 solid Solid_H wdata,restart,100,0,,
dens cons 100 powall,n,y,n
spec cons 1000 prfield,none,,,nouser
cond cons 10 prwall,n,y,n
moni 50
iter,500,0.001
!Gen Coldcells - Celltypes 3 and 4
ctype 3 save,,
vc3d 0.025 0.9 35 0 0.5 20 0 1 1 geom,,
vc3d 0 0.025 1 0.025 0.5 19 0 1 1 prob,,
Chapter 4 – Two-Dimensional Analysis 67
4.8.5. Summary
Comparing the numerical scheme described here, to that of the first numerical scheme, various
improvements for the case of a solid cooling insert are made. Many of the uncertainties were
removed and fewer were made. It however still only solves for the temperature distribution in a
two-dimensional slice and does not take into account the influence of the third dimension.
Until this influence is taken into account and results from two-dimensional analyses compared
with results from three-dimensional analyses, the accuracy of the two-dimensional approximation
cannot be guaranteed.
The two-dimensional code making use of the TDMA method generated thermal solutions, were
believed to be correct. The current two-dimensional simulation model can be extended into a
three-dimensional model. Because the ultimate analysis aim was to work with three-dimensional
temperature distributions, the results obtained from the current two-dimensional method is not
discussed. It was felt to be more important to develop a three-dimensional solid-state numerical
model, which could be used to do thermal optimisation of a heat-extraction system.
The use of a plain iterative solution method however still does not have the required performance
needed to do three-dimensional simulation time efficiently, and some improvements were needed
to speed up the solution process.
4.9. Conclusion
A representative corner quadrant domain was defined which could be used to describe the
temperature profile about a cooling structure insert. Two main theoretical approaches were
explored to find the temperature profile in this representative domain, namely analytical and
numerical solutions.
Obtaining an analytical solution for even the most simplified two-dimensional cases was believed
to fall outside the scope of this study and it was opted to use numerical approximations to obtain
the required temperature distributions. Due to the large number of anticipated simulations needed
to describe the temperature profiles and optimise the cooling structure distribution and shape, it
was decided to develop numerical computer code from scratch. Doing such large numbers of
simulations on commercially available numerical packages would have been very time
consuming.
A numerical scheme was developed using a generalised cooling structure of either a fluid filled
channel or a solid-state insert. From an investigation that was conducted with the aid of this
scheme, generalised trends were obtained for the behaviour of the optimum cross sectional shape
of the rectangular cooling structure as well as the temperature response.
It was found that the optimum cooling structure aspect ratio is dependent on physical geometric
conditions, material properties, as well as the heat transfer coefficient between the heat-
generating medium and the generalised cooling structure.
Increases in the aspect ratio of the domain, aD, the fraction of the domain used for cooling, , and
the ratio between the heat transfer coefficient and the thermal conductivity, [m-1], results in a
decreased maximum temperature. From a pure thermal point of view the results indicate that
Chapter 4 – Two-Dimensional Analysis 68
continuous cooling layers are the most efficient in reducing the peak temperature within a heat-
generating medium if a relatively small volume fraction is occupied by a cooling system.
It was found that solid cooling structures and fluid filled channel cooling structure have to be
investigated separately to achieve higher accuracies and reliability. A two-dimensional numerical
scheme was developed for the case where a solid-state cooling structure is used. A similar type
of temperature distributions was obtained from this scheme as before, with similar trends.
With the two-dimensional foundation been laid for a solid state cooling structure, it was possible
to extend this scheme into the third dimension and investigate the thermal behaviour of such a
configuration.
4.10. Nomenclature
4.10.1. General Symbols
Symbol Unit Description
A m2 area
a dimensionless rectangular aspect ratio defined as the x dimension over the y
dimension
m half of the centre-to-centre distance between neighbouring cooling
inserts in the x direction
m half of the x-direction dimension of the cooling insert
m half of the centre-to-centre distance between neighbouring cooling
inserts in the y direction
m half of the y-direction dimension of the cooling insert
Ci dimensionless element of a series of integration constants
G K/m ratio of volumetric heat gain over thermal conductivity
h W/m2K heat transfer coefficient
i dimensionless index number
k W/mk thermal conductivity
N dimensionless number of nodes or index terminator
q W heat transfer rate
q '' W/m2 heat flux over a surface
q ''' W/m3 volumetric thermal heat-generation
2
R m K/W thermal interfacial resistance
T ˚C or K temperature
U ˚C or K analytical transformation variable
X m full dimension of control volume in the x direction
x m Cartesian coordinate or a control volume dimension in x direction
Y m full dimension of control volume in the y direction
y m Cartesian coordinate or a control volume dimension in y direction
z m Cartesian coordinate or a control volume dimension in z direction
m z-direction dimension of domain of interest
Chapter 4 – Two-Dimensional Analysis 69
4.10.2. Greek Symbols
Symbol Unit Description
dimensionless fraction of representative domain used for cooling purposes
m-1 ratio between thermal conductivities of the heat-extraction material
and the heat-generating medium, or the ratio between the heat
transfer coefficient on the interface to the thermal conductivity of the
heat-generating medium
ζ dimensionless variable group used in the discretised equations
ψ integration function
4.10.3. Subscripts
Symbol Description
C cooling structure
D representative domain
E easterly reference direction
if interface between the heat-generating medium and the cooling structure
int internal (referring to the interface between the heat-extraction and heat-generation
mediums)
M heat-generating medium
N northerly reference direction
max maximum
min minimum
optimum optimum
rel relative
S southerly reference direction
W westerly reference direction
4.11. References
[1] Patankar S.V., Numerical heat transfer and fluid flow, Hemisphere (Washington D.C.) ,
1980
Chapter 4 – Two-Dimensional Analysis 70
Related docs
