Experimental and numerical study of the thermal and hydraulic effect of EMC screens in radio base stations: detailed and compact models

(1)

(2)

(3)

Abstract

Today’s telecommunication cabinets use Electro Magnetic Compliance (EMC) screens in order to reduce electromagnetic noise that can cause some miss functions in electronic equipment.

Many radio base stations (RBSs) use a 90-degree building architecture: the flow inlet is perpendicular to the EMC screen, which creates a complex flow, with a 90-degree air turn, expansions, compressions, perforated plates and PCBs. It is of great interest to study how the EMC screen interacts with the rest of components and analyze the total pressure drop and how much the flow pattern changes due to the placement of the screen.

Velocity, pressure and temperature measurements as well as flow pattern visualizations have been carried out to gain good insight into the flow and heat transfer characteristics in a subrack model of an RBS. Furthermore, these measurements have been very useful for validating detailed CFD models and evaluating several turbulence models.

Nowadays, industrial competition has caused a substantial decrease in the time-to-market of products. This fact makes the use of compact models in the first stages of the design process of vital importance. Accurate and fast compact models can to a great extent decrease the time for design, and thus for production.

Hence, to determine the correlations between the pressure drop and flow pattern on the PCBs as a function of the geometry and the Reynolds number, based on a detailed CFD parametric study, was one objective. Furthermore, the development of a compact model using a porous media approach (using two directional-loss coefficients) has been accomplished. Two correlations of these directional loss coefficients were found as a function of the geometry and Reynolds number.

(4)

Acknowledgements

First I would like to thank my supervisor Professor Bahram Moshfegh and co-supervisor Dr Hans Jonsson for all their help with this thesis. Without you, this thesis would not have been accomplished. Thanks, Bahram, for employing me and giving me the opportunity to work in this beautiful area. Thanks also for all your comments and discussions about results and analyses. Although our offices were in different cities, you were always available to find some time to discuss the project wherever you were. Thanks, Hans, for all your supervising work and for the encouragement during these years.

I would like to thank Professor Carlos Bastero, Professor Miguel Angel Serna and Professor Tomás Gómez Acebo from TECNUN (University of Navarra) who trusted in me and decided to give me a scholarship to come to KTH. Also thanks to Professor Björn Palm and Associate professor Per Lundqvist, who accepted me in their division at KTH. Thanks, Björn, for helping me at the beginning of my stay in Sweden and because you always have your door opened and available to help anyone who comes by. And Per, for being always of good humor!!

Thanks to the KK Foundation, Sweden and Antonio Aranzabal Foundation, Spain for economically supporting this project. Thanks to Mr Mikael Ardvisson and Dr Jukka Rantala and their staff teams at Ericsson and Nokia, respectively, for participating in this project and for all the good comments during all the presentations at Kista and Helsinki. Thanks to everybody at the Applied Thermodynamic division at KTH: Prof. Eric Granryd, Tim, Jaime, Claudi, Peter, Joachim, Anders, Rathmat, Oxana, Wimol, Cecilia, Samer, Branko, Rashid, Jan-Erik, Erik, Åke, Nabil, Yang, Emilio, Getachew, Shota, Teclemariam, Primal, Richard, Wahib, and of course our technicians, Benny S. and Benny A., the computer support team, Tony and Birger, and our efficient secretaries, Inga and Susy and in a special way to Miguel Castiella for the help with the visualizations. It has been a pleasure to enjoy a lot of meals, coffee times, kart driving, excursions and all kinds of football with

(5)

Thanks also to everybody at the Department of Technology and Built Environment, University of Gävle: Daniel, Ulf, Samuel and Sofia. It always felt so familiar to be with you during my short visits there. Special thanks to Dr Mathias Cehlin, who always finds some time to help everybody with the experimental stuff and to the most efficient technician I have ever met: Larry Smith, who did such nice work with the test rig. Thanks to Hans Lundström for all the help with the hot-wire anemometers, Tommy Gaude for the good work with the perforated plates, Beth Chapple for the English corrections, Mrs Eva Wännström, our helpful secretary and again Mathias and Daniel for reading carefully the manuscript of the thesis.

I extend also my gratitude to Associate Professors Alejandro Rivas and Juan Carlos Ramos that lead the Fluid Flow and Thermal Engineering Division of TECNUN (University of Navarra) and taught me in my undergraduate studies, and to the rest of the nice people in that Division with whom I am working now: Gorka, Jon, Mireia and Juan.

I would like to thank everybody with whom I have lived during these years in the best residence hall in the world (it is not an exaggeration): Lärkstaden Studiecentrum. All of you have really made me enjoy this stay in Sweden. You are my family in Sweden.

At last but not least, thanks to my parents Luis Miguel and Rosa María, and my sisters Judith, María Elena and Beatriz for everything: for life, love and all the encouragement during these years and many other things! To you this thesis is dedicated!

Raúl Antón Remírez Stockholm, October 2006

(6)

Nomenclature

A, area m2

a, inlet height m

Ainlet, inlet area m2

Apcb, cross sectional area between two PCBs m2

b, gap between sub-rack inlet and EMC screen m

D, hole diameter m

d, sub-rack depth m

h, heat transfer coefficient W/m2_ºC

k, turbulent kinetic energy m2_/s2

n, number of holes between 2 PCBs -

p, pressure Pa

P, average pressure Pa

p', fluctuating pressure Pa

P, modified kinematic pressure m2_/s2

q, heat flux W/m2 Q, heat flow W S, hole spacing m t. temperature ºC T, average temperature ºC t', fluctuating temperature ºC tpcb, PCB thickness m ts, screen thickness m u, velocity m/s U, average velocity m/s u', fluctuating velocity m/s

Wslot, distance between 2 PCBs m

x, coordinate along the sub-rack depth m

z, coordinate along the sub-rack height m

'p, pressure difference Pa

'T, temperature difference ºC

Dthermal diffusivity m2_/s

Hscreen porosity -

Hdissipation rate of turbulence kinetic energy m2_/s3

Qdynamic viscosity m2_/s

(12)

1 Introduction

1 . 1 T h e r m a l M a n a g e m e n t o f E l e c t r o n i c

S y s t e m s

Thermal management of electronic systems is one of the bottlenecks in today’s electronic industry. The increasing clock speeds and miniaturization implies a huge increase in the heat dissipation per unit area in electronic systems. It is a real challenge to cool electronic components in an efficient way. Air is still the most common method for cooling electronics, due to safety reasons and availability. For some high-dissipation cases, liquid cooling cannot be avoided any more, since there is no way to cool very high heat fluxes with air. However, most electronics are still cooled by air.

Thermal management consists of reducing failures due to thermal-related problems in such a way that a good level of reliability is achieved. This means assuring the good performance and long life of the electronic product. The normal way to achieve this is to make sure that the junction temperature of all electronic components is not excessive.

Heat dissipation from a surface may be described by Newton’s law of cooling:

T hA

Q '

According to this law, there are four ways to try to improve thermal management. One way is to reduce the heat dissipation, for example reducing the voltage in the electronic components, as has been done during the last few decades (Kim and Kim, 2004). A second way is to increase the heat transfer coefficient (see Figure 1a). A lot of research has also been performed in this direction: air jet impingement (Rundström and Moshfegh, 2006), liquid cooling, etc. A third way is by increasing the dissipation area (Jonsson and Moshfegh, 2001), e.g. see

(13)

difference by reducing the surrounding temperature (see Figure 1c) or to use new materials inside the components that withstand higher temperatures.

To ensure that the junction temperature is not excessive, there must be an adequate heat transfer coefficient on the surfaces of the chips and a correct temperature surrounding the component. The correct component ambient temperature will be achieved only if there is good thermal management on the board level. In order to achieve a convenient inlet temperature on the board level, good thermal management should be achieved at the sub-rack level as well as at the cabinet level. This chain of temperature level requirements eventually leads to achieving the right temperature in the surroundings of the cabinet: room temperature. To set the room temperature to a particular value will normally be accomplished by using air conditioning units.

1 2 3 4 5 6 a ) S p r a y c o o l i n g b ) H e a t s i n k c ) A c t i v e c o o l i n g i n a d a t a c e n t e r F i g u r e 1 . S e v e r a l w a y s t o i n c r e a s e t h e c o o l i n g i n a n e l e c t r o n i c s y s t e m .

This thesis is focused on the system level: a sub-rack of a radio base station. However, the analysis that has been performed in this thesis is not only valid for a RBS but also for all the telecommunication sub-racks with a cooling architecture in which the airflow makes a 90 degrees turn.

1 . 2 R a d i o B a s e S t a t i o n : A C o o l i n g

B u i l d i n g A r c h i t e c t u r e

In a telecommunication network (e.g. a mobile telephone network), the RBSs are the nodes of the network, where wireless signals are received and transmitted. The RBSs have several shelves called sub-racks (see Figure 2b), filled with cards with electronic components (see Figure 2c). These boards are the power board assemblies (PBAs) where components are mounted on printed circuit boards (PCBs). In most of the cases, RBSs are cooled by forced air convection.

(14)

a ) C a b i n e t l e v e l b ) S u b - r a c k l e v e l c ) B o a r d l e v e l F i g u r e 2 . S e v e r a l l e v e l s i n t h e t h e r m a l m a n a g e m e n t o f a c a b i n e t .

In a typical RBS building architecture, the sub-racks are cooled individually, in such a way that the same amount of fresh air enters into each sub-rack. At the top of the sub-rack, several blowers (a radial fan tray) that suck the air (a pull configuration, see Figure 2b) are placed. Furthermore, two Electro Magnetic Compliance (EMC) screens are placed at the inlet and the outlet of the sub-rack slots. The function of this EMC screen is to decrease electromagnetic noise for safety and good performance reasons. In a sub-rack (see Figure 2b), the air first enters through the inlet and makes a 90-degree turn while impinging and going through the EMC screen. Then the air flows past the electronic components, thereby removing the heat. The air finally goes through the second EMC screen and is drawn out by the fans and ejected to the back chimney of the cabinet, ending in the cabinet surroundings.

This pull configuration is more typical than the push configuration (where axial fans are placed in the inlet plenum), since it creates a more uniform flow over the PCBs.

1 . 3 E l e c t r o m a g n e t i c C o m p l i a n c e

The desire to reduce the size of electronic systems results in closer placement of many electronic components. This proximity may create problems due to electromagnetic interference.

Noise is an undesired electrical signal in a circuit; interference is the effect of the noise. Noise should be reduced to a level at which it does not create interference.

Requirements for electromagnetic compatibility demand that electronic equipment works properly in its electromagnetic environment and that the equipment is not a source of electromagnetic interference.

(15)

A way to control the propagation of magnetic fields is by means of a shield. A shield may be a metallic perforated plate that efficiently hinders electromagnetic radiation going through it. In a sub-rack, this means to hinder electromagnetic radiation from leaving the sub-rack and electrical disturbances from entering it. The ideal shield would be to completely enclose the electronics in a metallic box. However, this prevents cooling, so there should be some apertures, although this means some electromagnetic leakage.

The level of leakage of electromagnetic radiation from the shield depends mainly on three parameters (Ott, 1988): the maximum linear dimension, the wave impedance (ratio of the electrical and magnetic fields) and the frequency of the source. A typical shield is a perforated plate (an EMC screen). It has small openings but enough free area that permits cooling of the electronic systems. For commercial products, a good rule of thumb is to avoid openings greater than 1/20 of a wavelength.

1 . 4 M o t i v a t i o n f o r T h i s S t u d y

There are three different motivations for this thesis: general reasons related to the selection of fans and to flow distribution problems, as well as a wish to study complex flow phenomena that appear in this building architecture and finally to develop compact models for a more efficient design process.

1 . 4 . 1 G e n e r a l r e a s o n s

It is important to know the pressure drops in the RBS cabinet in order to know which size of fans are required. In general, the pressure loss in a system increases with an increase in the flow rate, and it is important to match the pressure increase of the fan to the pressure losses caused by the flow through the system. It is therefore important to identify where those pressure losses occur.

For low screen porosities, the pressure drop is mainly due to the low porosity of the screen. In this case, the volume before the EMC screen has a tendency to act more as a plenum, and therefore the EMC screen behaves as a honeycomb flow straightener. This leads to a more uniform flow distribution (flow pattern) after the screen. On the other hand, high porosity EMC screens lack this feature, and thus a less uniform flow pattern is obtained between the PCBs.

(16)

The reason why information about the flow pattern is so important is that high dissipative components should be placed in locations with high velocity (i.e., in the bulk flow), and not in areas with very low velocities. This is a consequence of the strong dependence between the heat transfer coefficient and the velocity.

1 . 4 . 2 C o m p l e x f l o w p h e n o m e n a

According to a recent article (Joshi et al., 2002), there are still gaps in system level modeling. One of the main issues to be solved is to eliminate uncertainty in input parameters like the screen pressure loss coefficient. For complex flows as in an RBS sub-rack, it is difficult to calculate the pressure drop and to predict the flow pattern in between the PCBs. To simply add pressure loss coefficients for flow obstacles in the flow path will in general produce inaccurate results for flow obstacles placed close to each other, due to such factors as hydraulic interference between the obstacles (Smith, 1971; Fried and Idelchik, 1989).

Since the air entering the sub-rack makes a 90-degree turn, the flow does not impinge the screen at the same velocity. The impinge angle may also vary along the sub-rack depth, as shown in Figure 3. Furthermore, the flow will experience a sudden change in direction after colliding with the interior faces of the EMC screen holes.

F i g u r e 3 N o n u n i f o r m v e l o c i t y p r o f i l e ( c o l o u r e d b y v e l o c i t y m a g n i t u d e ) b e f o r e a n d a f t e r t h e E M C s c r e e n

1 . 4 . 3 C o m p a c t m o d e l l i n g

Nowadays, when the design time is cut due to industrial competition, thermal management engineers do not have as much time as they want to analyze the thermal problems. In order to have an efficient design process, computers performing simulations of the behavior of the

(17)

dynamics (CFD) simulations. Moreover, much more information can be extracted from the CFD simulations than would be possible with normal experimental investigations. However, it is a challenge to choose the proper turbulence model, near wall treatment and mesh strategy for a particular application, and it is also quite a time-consuming activity. Therefore, engineers tend to minimize their models by employing compact models of certain components, e.g. heat sinks, components or, as in this case, EMC screens. The importance and benefits of compact models are many. Primarily, the use of compact models reduces the time needed to get a good design for a given application and an accurate compact model of an EMC screen would be of great importance and benefit to the industry.

1 . 5 A i m o f t h e S t u d y

This thesis is focused on analyzing the heat and fluid flow inside an RBS sub-rack, and mainly analyzing the effect of the first EMC screen on the flow pattern, pressure drop and temperature field in a PCB sub-rack. The ultimate goal is to predict the total pressure drop and the flow details after the EMC screen in a PCB slot of a RBS sub-rack by means of a fast and accurate tool either based on correlations or by a CFD compact model based on a porous media approach.

1 . 6 M e t h o d s

There are two main methods to study an electronic system or a flow problem in general: experimental and theoretical. To use the experimental method is very time consuming; normally one does not get a whole field measurement, but local or average data are obtained. Many times, it is necessary to study the system at a reduced scale. Other times, in order to get more accurate measurements, a wind tunnel is used, or a simplified model is built. Anyway, experiments are still the most reliable method. Hence, the experimental method is not only used to gain insight into the heat and fluid flow in the system but it is very often used to check the validity of a numerical result.

A theoretical method is based on the laws of nature. Sometimes, it is possible to get an analytical solution, but many times a numerical solution is needed, due either to the mathematical complexity of the laws (e.g. nonlinear partial differential equations) or the complexity of the geometry where the laws are applied. This numerical approach implies an approximation of the analytical solution and thus it should be validated by experiments.

(18)

The positive feature of the theoretical approach is that we obtain detailed information about the whole domain (sometimes the discretized domain), it is relatively cheap and fast. Hence, it is interesting to work with both methods: to validate several CFD models by experiments and then a larger range in the analyzed parameters may be studied by CFD. This has been the approach in the present thesis.

1 . 7 L i m i t a t i o n s

In the present case, the reference values have been the experimental studies. These studies were not performed in a real RBS sub-rack but in a model inside a wind tunnel. The size of the sub-rack model is close to standard RBS sub-racks. The position of the model is horizontal, i.e., the PCBs are horizontal. In reality, the normal position is vertical. The reason why the sub-rack was placed horizontally was that an entrance duct of one meter was placed at the entrance of the sub-rack model, as can be seen in Figure 4. The entrance duct was used to get low turbulence intensity level (fully turbulent developed flow) in order to have reliable hot-wire measurements. Due to the 90-degree turn, the velocity profiles at the inlet of the sub-rack are neither uniform nor symmetric and it was important to measure these profiles in order to use accurate boundary conditions in the numerical study. This approach is beneficial, since the profiles at the inlet of the sub-rack are asymmetric. However, this is also a limitation, because in reality there is no entrance duct, and the flow would not always be perpendicular, which may affect the pressure drop. Anyway, it is thought to be a good approximation, and a repeatable method.

F i g u r e 4 . W i n d t u n n e l s e t - u p .

Another approximation is that the sub-rack model analyzed in this thesis lacks swirl boundary conditions at the outlet. It is supposed that the second screen (not studied in this work) will break the swirl produced by the blowers at the top of the sub-rack. In Nevelsteen et al. (2006) the

(19)

screen. It can be seen from their results that the swirl is actually reduced after the EMC screens. Although the building architecture studied in this thesis is a pull configuration, Nevelsteen et al.’s results can justify this approximation.

1 . 8 R e s e a r c h P r o c e s s

The process of the thesis is depicted in Figure 5. First an experimental parametric study was done (paper I, and the experimental part of paper II). In the experimental study, pressure and velocity measurements were performed at several locations before and after the screen for several velocities, sub-rack geometries, as well as screen porosities. The influence of these parameters on the flow pattern was also studied using the smoke-wire technique. All these experiments were used to validate detailed CFD models, and at the same time to evaluate different turbulence models and near wall treatments used in the numerical investigations in papers III and IV.

Experimental Parametric Study

To validate

Detailed CFD Parametric Study (velocity profile at the inlet of the subrack)

Detailed CFD Parametric Study (uniform velocity profile at the

inlet of the entrance duct)

Comprehensive Detailed CFD Parametric Study (uniform velocity profile at the

Compact model of the EMC screen by a porous media approach

Comprehensive Compact CFD Parametric Study (uniform velocity profile at the

Design Tools

To validate

To create

F i g u r e 5 . R e s e a r c h p r o c e s s .

In CFD simulations, accurate boundary conditions are needed. In order to use accurate boundary conditions, the velocity and turbulence

(20)

intensity profiles at the inlet of the sub-rack were used (paper III). In order to acquire reliable measurements, an entrance duct was placed before the sub-rack inlet in such a way that the profiles were fully developed turbulent profiles.

It was interesting to study the influence of different geometrical parameters and velocities for a greater range. However, the velocity profiles were needed at the sub-rack inlet and thus the study was not possible to perform. However, it was shown in paper IV that by increasing the CFD model using an entrance duct of 0.5 meters and using uniform velocity and turbulence intensity profiles at the inlet of the entrance duct, a good prediction of the velocity and turbulence intensity profiles at the inlet of the sub-rack was achieved.

A comprehensive parametric study could be performed with this approach, as was done in paper VI. Before performing this parametric study, an evaluation of turbulence models using an enhanced wall treatment (see details in section 3.5.1) was done in paper IV. The best of these models was also used for the analysis of the conjugate heat transfer as a function of several parameters (paper V). After that, a large number of isothermal detailed CFD simulations were done and correlations were developed to predict the total pressure drop. Furthermore, to characterize the flow pattern after the screen, two parameters were defined that account for the wetted area and the flow uniformity after the screen. Correlations for these two parameters were also developed. Since it can actually be said that these correlations are some kind of compact model, the objective of the thesis was partially fulfilled. However, a further step could be done: to try to create a compact model using a porous media approach (paper VII) and to validate the model by means of the comprehensive detailed CFD parametric study.

The aim of the compact model is to obtain the pressure loss coefficients of a volume hydraulical impedance that defines the EMC screen. The use of these pressure loss coefficients provides an accurate prediction of the pressure drop and the flow pattern after the screen. Correlations that predict the two directional loss coefficients as a function of the geometry and Reynolds number were developed. One positive aspect to this porous media approach is that the model can be completed by setting electronic components and other devices, but knowing that the flow pattern and pressure drop due to the 90-degree turn and screen porosity to be predicted faster and more accurately. The last step was to develop programs in a user-friendly environment that give all the information

(21)

way and as a function of the geometrical parameters: inlet height, PCB thickness, screen thickness, distance between inlet and screen, screen porosity, distance between two PCBs, and the Reynolds number.

1 . 9 L i t e r a t u r e S u r v e y

There are enormous numbers of papers that focus on thermal management of electronics equipment. There are also plenty of papers that analyze the pressure drop through perforated plates. However, very little has been written about the effect of an EMC screen in a building architecture in which the flow makes a 90-degree turn and impinges the screen at an angle from the normal plane. To the best of my knowledge, there are no compact models for EMC screens for these configurations.

1 . 9 . 1 P e r f o r a t e d p l a t e s

The pressure drop through perforated plates has been extensively studied by Idelchik (1994), who wrote a handbook for pressure loss coefficients that is very often used in the thermal management community. Other authors that have studied this issue are Carrothers and Baines (1965) and Smith (1971). Carrothers and Baines found the pressure drop through an inclined screen proportional to the square of the velocity component perpendicular to the screen. Smith published a book with the state of the art for correlations of the pressure drop through perforated plates, in which with the aid of experimental visualization he shows that porosity and screen thickness are the most important parameters for the prediction of pressure drop at high Reynolds numbers. Furthermore, Gan and Riffat (1997) analyzed the pressure drop through perforated plates experimentally and numerically. They also showed how the dimensionless pressure drop through the screens decreases as the ratio of screen thickness and hole diameter increases.

1 . 9 . 2 E M C s c r e e n s

There have been earlier attempts to model EMC screens for particular cases, as in Nevelsteen et al. (2006), where the screens close to axial fans were modelled by a compact model with three directional loss coefficients. Baelmans et al. (2003) analyzed experimentally the distance of influence by the screen, remarking on the difficulty of predicting the flow beside components placed close to the screen when the screen is modelled by means of hydraulic surface resistance. No correlations were developed in these works. Kordyban (2000) emphasized the need to use volume hydraulic impedance and not planar impedance to model EMC

(22)

screens or filters. He stated also that using a volume hydraulic impedance with only one coefficient is wrong. A one-coefficient model would reduce the velocity component perpendicular to the screen but not the other component, and thus the inclined flow would turn to the wrong direction. It was also stated that on some occasions the in-plane pressure loss coefficient may be larger than the perpendicular loss coefficient, however no correlations were developed in his study either.

Akella and Ortega (2005) analyzed the pressure loss coefficients in a slot (1 RU card passage). In this case the air was perpendicular to the screens. They measured the pressure loss coefficient for the first screen, for the board populated with electronic components, and finally for the second screen. They said that there may be interaction between the different components, but if these interactions are small a flow network modeling (FNM) approach would be viable for this case. They showed also that the pressure loss coefficients were influenced very little by the Reynolds number.

Bejan et al. (1995) used the intersection-of-asymptotes method applied to the cooling of stacks of plates shielded by porous screens, but in this case the air flow was also perpendicular to the EMC screens.

1 . 9 . 3 C F D a n d e x p e r i m e n t s i n e l e c t r o n i c s y s t e m s

w i t h s c r e e n s

Many studies have been performed numerically by means of CFD, for example Nakamura and Komura (1996), Biswas et al. (1999) and Tahmaspur and Berhe (2002). Nakamura and Komura studied the airflow in a sub-rack with and without EMC screens. The study used CFD and measurement techniques such as LDA (Laser Doppler Anemometry). Biswas et al. also performed experiments and CFD simulations in an electronic box in which there were grilles, however these grilles were modelled using pressure loss coefficients from a handbook. Tahmaspur and Berhe compared honeycomb and perforated vents in a TW-150 Passive Optical Network. They did a detailed model of the perforated plate in order to estimate the pressure loss coefficient of that screen in that application.

1 . 9 . 4 C o m p a c t m o d e l

Nakayama et al. (2001 and 2004) present a methodology in which, after spending time in developing detailed models, one can compress the information in a database and create compact models. The compact

(23)

temperature is not exceeded. Compact models have been developed not only for electronic packages like in Vinke and Lasance (1997), but also for subsystems like heat sinks. Narasimhan et al. (2003) developed a compact model of heat sinks using a porous media approach that was based on a boundary layer model between the plate fins of the heat sinks. Recently Kim and Kim (2004a) have developed a compact model based on a volume average technique, which has been applied for straight fin heat sinks as well as for pin fin heat sinks (Kim et al. 2004b). Another approach to shorten the design time is optimization methods, such as in Parry et. al. (2004), where a methodology is developed for getting the best design that fulfils an objective function and subject to some constraints, in which some of them may be integers like in the case of an optimized heat sink. Another way of compact modelling is flow network modelling, with early works like Butterbaugh and Kang (1995). Kowalski and Radmehr (2000) analyzed the flow in a cabinet in which they combine a flow network analysis together with some CFD simulations for some details (to get the impedance curve of the axial fans that they used). However, it should be emphasized that this method cannot be used when there is hydraulic interference between the different flow components. It is used between locations in which the pressure is uniform, without gradients, for some complex geometry may cause inaccuracies with this method.

(24)

2 Experiments

The experimental measurements in this thesis have two purposes. The first goal is to gain insight into the hydraulic and thermal behavior of the flow inside the sub-rack model. The second objective is to use the experimental measurements as a reference with which to validate the detailed CFD predictions.

2 . 1 E x p e r i m e n t a l S e t - u p

The experiments were performed in a sub-rack model inside a wind tunnel. In the following section the wind tunnel is presented.

2 . 1 . 1 W i n d t u n n e l

The general layout of the wind tunnel is shown in Figure 6.

Fan Laminar flow element Entrance duct Inlet Subrack Screen F i g u r e 6 . W i n d t u n n e l l a y o u t .

The wind tunnel walls are made of Plexiglass. The air velocity is controlled by frequency regulation of the electric power driving the fan.

(25)

A sub-rack was modelled by means of Plexiglass plates (4 mm thick) as is shown in Figure 7 and Figure 8 (all the dimensions in mm). There are a total of 14 PCB dummies. F i g u r e 7 . A d e t a i l e d v i e w o f t h e s u b - r a c k s l o t s . 265.5 280 15 50 Screen Side view Inlet d 1000 Top view x z F i g u r e 8 . S u b - r a c k m o d e l .

Two models were built with the same dimensions except for the depth,

(26)

The size of the commercial sub-racks follows the standards for electronic packaging that have been developed by the International Electrotechnical Commission (IEC). The employed unit for the height of the rack is U=44.45 mm (17.5 in). The typical values for the sub-rack height are 3U, 6U and 9U. The corresponding PCB heights for these sub-racks are (in mm): 100, 233 and 366. In the model for this thesis, the PCB height was always 280 mm. The typical PCB depths for these sub-racks are 160, 220 and 280 mm. In the model for this thesis, the PCB depth was either 200 mm or 260 mm. A typical distance between two PCBs is 3×5.08 mm (3×2 in); in the model for this work it was 14 mm. As can be seen, the tested sub-rack geometry is close to that of commercial sub-racks and based upon discussions with Nokia and Ericsson. A commercial sub-rack appears in Figure 9.

F i g u r e 9 . A t y p i c a l c o m m e r c i a l 6 U s u b -r a c k .

2 . 1 . 2 E M C s c r e e n s

Stainless steel screens (1 mm thickness) with 60° staggered circular holes were used (see Figure 10). The location of the screen is indicated in Figure 8. Three different porosities, H, have been used for the screen, shown in Table 1

The EMC screens were manufactured with a Computer Numerical Controlled (CNC) drilling machine. After the drilling a sand blaster was used to remove the remaining steel debris in a few holes.

(27)

T a b l e 1 . S c r e e n c h a r a c t e r i s t i c s

D (mm) S (mm) D+S (mm) _H(%)

3.25 1.25 4.5 46.9 3.50 1.00 4.5 54.3 3.75 0.75 4.5 62.2

The porosity in the above table is calculated according to:

d W D n /4 slot 2 S H (1) 4 S D 18 F i g u r e 1 0 . M o d e l o f t h e s c r e e n a n d P C B s l o t .

2 . 2 D a t a R e d u c t i o n

The following nomenclature designates the different geometries and velocities:

V1, V2, V3, V4, and V5 are the five velocities studied in the wind tunnel. V1 is the minimum (1.5 m/s) and V5 is the maximum (7.5 m/s). More details follow in Boundary conditions (section 3.6).

p1, p2, and p3 are the three screen porosities (see Table 1. Screen characteristics) studied. p1 is the maximum porosity (62.2%) and p3 is the minimum porosity (46.9%).

q1, q2, and q3 are the heat fluxes used in the PCBs. q1 is the minimum heat flux (266.5 W/m2) and q3 is the maximum one (431.5 W/m2). More details can be found in Boundary conditions (section 3.6).

S designates the short model, i.e., sub-rack depth equal to 200 mm; L refers to the large model, i.e., sub-rack depth 260 mm.

(28)

2 . 3 V e l o c i t y

2 . 3 . 1 E x p e r i m e n t a l p r o c e d u r e — i n l e t p r o f i l e

The system was always allowed to stabilize for approximately two minutes before taking the measurements with a cross hot-wire anemometer. The sample frequency used was 4 kHz and the sample time 20 sec. The average velocity and turbulence intensity were always obtained.

F i g u r e 1 1 . T r a v e r s i n g s y s t e m u s e d t o m e a s u r e t h e v e l o c i t y p r o f i l e s a t t h e i n l e t o f t h e s u b - r a c k .

Twenty-one measurements along the inlet width were measured in order to get an accurate profile. A Mitutoyo traversing system (see Figure 11) was used with 0.01 mm digital resolution and 0.254 mm accuracy. The accuracy for the point closest to the wall may be worse. These profiles were used through the project; the average of those velocity profiles agrees ±5% with the laminar flow element. The laminar flow element measures the average velocity in the whole sub-rack, but the hot-wire measurements were done at the middle of the sub-rack and thus some disagreement with the laminar element was expected.

2 . 3 . 2 E x p e r i m e n t a l p r o c e d u r e — i n t e r i o r p o i n t s

When measuring the local velocities, the positioning of the probe was done manually. First the probe was put in its position, and then the control system for the wind tunnel cycled through five preset wind tunnel velocities. The wind tunnel was then turned off, the probe moved

(29)

tunnel velocities were obtained as described earlier. The sample frequency and the sample time were 4 kHz and 20 sec, respectively. The measurements before the screen were done along half of the depth of the screen (locations b1, b2 and b3 in Figure 12).

z Entrance duct inlet x Flow direction EMC screen z _{velocity point} pressure tap x 10a 9a 8a 7a 6a 5a 4a 3a 2a 1a 10b 9b 8b 7b 6b 5b 4b 3b 2b 1b L R T b1 b2 b3 Subrack Inlet a1 a3 a5 a2 a4 F i g u r e 1 2 . P r e s s u r e a n d v e l o c i t y l o c a t i o n s .

The velocity measurements after the screen were done at 14.5 cm behind the screen (locations a1, a2, a3, a4 and a5 in Figure 12) making sure that the turbulence intensity was below 25%. The velocity measurements before the screen were performed at a distance about 21 mm perpendicular from the screen. Although it is not stated in paper I, some of the velocity measurements after the screen (a2, a3, a4 and a5) may not be trustworthy because of the obstacle that the hot-wire body probe produces on the flow. However, it is believed that the first reading (a1), in which the probe body is just 2 mm, can be trusted.

2 . 3 . 3 U n c e r t a i n t y A n a l y s i s

The Dantec calibration unit was used to calibrate the probe with a fourth-order polynomial curve and with linearization errors less than ±1%. During tests, the air temperature at the wind tunnel inlet was kept at around 21°C. Recalling that the wind tunnel is turned on and off for each position of a local velocity measurement, the repeatability of the wind tunnel should also be considered. The velocity was measured during several days at several occasions and the estimated uncertainty in the velocity is ±2.5%. In Table 2 the total estimated uncertainty is calculated as the square root of the sum of the squares of the maximum error of the hot-wire anemometer and wind tunnel repeatability.

(30)

T a b l e 2 . S u m m a r y o f t h e t o t a l e s t i m a t e d u n c e r t a i n t y f o r t h e v e l o c i t y m e a s u r e m e n t s

Maximum error from calibration ±1% Wind tunnel repeatability ±2.5%

Total ±2.7%

In papers I to V the same velocity boundary conditions were used. This boundary condition comes from the hot-wire anemometer measurements. As the thermal and pressure measurements were done on different occasions, the wind tunnel repeatability should be taken into account. A flow laminar element (see Figure 6) was used as a checking device, making sure that the velocity was within the wind tunnel repeatability, as it actually was.

2 . 4 P r e s s u r e

2 . 4 . 1 E x p e r i m e n t a l p r o c e d u r e

The pressure taps are placed flush at the middle of the sub-rack bottom wall at 65 mm before the screen, and flush in a PCB at 101.5 mm after the screen (see Figure 12). The pressure transducer is connected to a Campbell scientific 21X logger, which saves the data in the computer. The pressure measurement procedure described here is slightly different from the one used in paper I, where some valves were used. Some small leakage was found in some valves that were used for pressure measurements in paper I. Hence, the pressure measurements were repeated but without the valves.

A pressure reading consisted of an average of 500 samples that were taken over 40 seconds. This process was repeated 5 times, once every minute to make up one series. Three series were obtained every ten minutes. Then the average of all the 15 values was taken.

2 . 4 . 2 U n c e r t a i n t y A n a l y s i s

The pressure taps are connected to an Autotran 750 pressure transducer by approximately 0.8 meters of tubing. The accuracy (maximum error) according to the manufacturer is ±1% of the full-scale reading of 25 Pa. For the highest velocity (if pressure readings were larger than 25 Pa) a similar pressure transducer with a full-scale reading of 100 Pa was used; the accuracy in this case was also ±1% of the full-scale reading. According to the measurements the pressure uncertainty (twice the

(31)

velocity (V1§1.6m/s) to an average of ±0.5 Pa at the highest velocity (V5§7.5m/s).

Since the velocity and pressure drop measurements were not obtained simultaneously, the repeatability of the wind tunnel may also affect the pressure drop measurement uncertainty. The pressure drop is approximately proportional to the square of the velocity. Hence, there is an uncertainty in the pressure drop due to the wind tunnel repeatability being estimated to twice the wind tunnel repeatability, i.e., ±5% (as shown above). The total estimated uncertainty is calculated as the square root of the sum of the squares of the maximum error of the pressure transducer, twice the wind tunnel repeatability and an average uncertainty during a run of the wind tunnel, i.e., around 0.3 Pa for the lowest velocity and 1.5 Pa for the highest velocity.

2 . 5 T e m p e r a t u r e

In this section, temperature measurement on a PCB dummy used in paper V is discussed.

2 . 5 . 1 E x p e r i m e n t a l p r o c e d u r e

The temperature was measured using type T thermocouples connected to a Campbell scientific 21X logger via a Campbell multiplexer. Sixteen thermocouples were attached, eight at each side of the PCB (in order to duplicate the eight measurement locations). The thermocouple tips were soldered to the copper inside a 0.5u1mm channel. After the thermocouple attachment to the copper, the thermocouple channel was filled with conductive grease, leaving a smooth PCB surface. See paper V for more details and the location coordinates of the thermocouples. A PCB consisted of two copper plates of around 2 mm thickness in which had been milled a 0.4×1.5 mm channel (see Figure 13). In one of them, an electrically insulated heater cable (0.6 mm diameter) is placed and finally both plates are glued with epoxy and the help of a hydraulic press. In five of the fifteen PCBs the same heat flux was applied. The heat applied in each of them was the same, because the same length of heater cable was used and due to the serial connection, meaning the same current passed through each of the PCBs. Five heated PCBs were used to achieve heat symmetry conditions in the PCB in which the measurements took place (the one in the middle of the five).

(32)

The measurements were done automatically in such a way that a logger takes 10 samples (once every minute) if and only if the difference in temperature between two consecutive series is less than 0.05ºC in all the measurement locations at the same time. At this time steady state was assumed.

F i g u r e 1 3 . D e t a i l s o f a h a l f P C B w i t h t h e m i l l e d c h a n n e l a n d t h e c o p p e r P C B s i n s i d e t h e s u b - r a c k m o d e l .

2 . 5 . 2 U n c e r t a i n t y A n a l y s i s

The uncertainty of the thermocouples was estimated after calibration in an ice bath to ±0.1ºC. The uncertainty in the measurements is calculated by twice the standard deviation of the 10 temperature samples taken at each point; an average of this uncertainty for all the points and cases is equal to ±0.2ºC. As thermocouples were set at both sides of the PCB, at the same distance from the center plane of the PCB, the average difference between both sides’ measurements for all the locations and cases was ±0.2ºC. The total estimated temperature uncertainty is determined by the square root of the square of the previous uncertainties and is equal to ±0.3ºC.

The heat losses were calculated by a three-dimensional CFD conduction analysis. The heat transfer coefficient on the PCB surface and the bulk temperature inside the PCB slot were taken into account. These two parameters were obtained from the detailed PCB slot model. Furthermore, the temperature between the double walls at both sides of the model was measured by thermocouples and a heat transfer coefficient (natural convection and radiation) equal to 6 W/m2_{K was} assumed. In Figure 14 a detail of the model used to calculate the heat

(33)

losses is presented. The model size was 385,000 cells and the double precision Fluent solver was used.

F i g u r e 1 4 . 3 D c o n d u c t i o n p r e s s u r e l o s s e s m o d e l ( m o d e l i n f i g u r e i s m i r r o r e d a r o u n d t h e s y m m e t r y p l a n e a n d t h u s c o n t a i n s 3 8 5 , 0 0 0 x 2 c e l l s ) .

The estimated heat losses are defined as the heat leaving through the plexiglass out of the sub-rack model and expressed as a percentage of the heat load applied to the PCB. Initially, the radiation was neglected, and the error due to this assumption has been estimated and also expressed as a percentage of the heat load applied to the PCB. These losses (conduction and radiation) ranged between 1.4% for the maximum velocity cases and 3.8% for the minimum velocity cases. Some further simulations have been done that take into consideration these conduction and radiation heat losses in order to quantify the error in the PCB temperature prediction. As can be seen in Table 3, the predictions for the high velocity (V5) are even closer to the experimental measurements, but for the low velocity level a worse result is obtained. The average disagreement in Table 3 is defined as the average disagreement of the temperature measured and predicted at the eight locations. Anyway, these heat losses do not give rise to a large uncertainty in the prediction of temperature.

T a b l e 3 . T e m p e r a t u r e a v e r a g e d i s a g r e e m e n t ( % i n º C ) b e t w e e n e x p e r i m e n t s a n d C F D p r e d i c t i o n s ( R N G kH m o d e l ) , t a k i n g i n t o c o n s i d e r a t i o n t h e h e a t l o s s e s p2-V5-q1 (%) p2-V5-q3 (%) p2-V1-q1 (%) p2-V1-q3 (%) Nominal case 2.40 1.96 -4.60 -3.02 Nominal case

without heat losses

1.9 1.24 -6.69 -5.75

Plexiglass

Copper

Heater Air

(34)

Further, the temperature and velocity measurements were not taken simultaneously. This means that the repeatability of the wind tunnel should be considered. As stated before, the velocity was measured during several days at several occasions and the estimated repeatability is 2.5% of the velocity.

Further simulations have been performed in order to study the average disagreement in the eight temperature locations when the velocity is varied ±2.5% of its nominal value in 4 cases: maximum and minimum value of the analyzed velocities and heat fluxes with the intermediate screen porosity. In Table 4, the average temperature disagreement at eight locations between the experimental values and CFD predictions, in which the nominal velocity is increased or decreased by 2.5%, are shown. It can be inferred that this uncertainty in the velocity does not give rise to a large uncertainty in the prediction of the temperature.

T a b l e 4 . T e m p e r a t u r e a v e r a g e d i s a g r e e m e n t ( % i n º C ) b e t w e e n e x p e r i m e n t s a n d C F D p r e d i c t i o n s ( R N G kH m o d e l ) w i t h a c h a n g e i n t h e n o m i n a l v e l o c i t y b y 2 . 5 % p2-V5-q1 (%) p2-V5-q3 (%) p2-V1-q1 (%) p2-V1-q3 (%) Nominal velocity 2.40 1.96 -4.60 -3.02 Nominal velocity +2.5% 1.83 1.13 -5.57 -4.19 Nominal velocity -2.5% 3.11 2.83 -3.58 -1.79

The natural convection was ignored due to the position of the PCBs in the test rig. They were placed horizontally and that means that the Grashoff number is very small if we compare it with the square of the Reynolds number (Gr Re2was always less than 0.059); thus, natural convection was ignored.

Anyway, it was examined whether natural convection would be important if the PCBs were placed vertically as in reality. That was not the case studied in paper V, however it is interesting to know if large differences may be expected or not. In order to study that, the worst case (S-p1-V1-q3) was investigated, taking into consideration the gravity and

(35)

and without taking into consideration the buoyancy effect in the vertical PCB was equal to 1.3ºC (or 1.7%). Once the porosity increases, the input heat flux is reduced, or the inlet velocity is increased this difference is greatly reduced. At medium or high velocities the effect of buoyancy can totally be ignored.

2 . 6 F l o w P a t t e r n

2 . 6 . 1 S m o k e v i s u a l i z a t i o n

The smoke-wire technique consists of coating the stretched wires with oil and heating them up. After heating them, the oil evaporates and immediately after that the oil condenses, forming smoke-like micro-droplets. The smoke was illuminated by xenon light. The visualized flow pattern was recorded on video. The experimental set-up is presented in Figure 15.

(36)

Five wires were placed in the middle between two PCB dummies that are located in the middle of the sub-rack model (see Figure 16).

Bulk flow a b c wire 1 wire 2 wire 3 wire 4 wire 5 w 1 w 2 w 3 w 4 w 5 d Recirculation Bulk flow Frontier point z coordinate x coordinate F i g u r e 1 6 . L o c a t i o n o f t h e w i r e s a f t e r t h e s c r e e n .

The procedure for generating smoke is presented in Figure 17. First, the wires were inserted inside the sub-rack through 0.6 mm holes. Two cotton pieces with oil were attached at the plexiglass walls. The wire was long enough (more than twice the distance between the cotton pieces) to be easily moved and impregnated with oil. In this way, it was not necessary to open the wind tunnel in order to impregnate the wire with oil. Flow Weight Wire Syringe Cotton is stuck

on the wall Horizontal wire movement_{in order to wet the wire}

F i g u r e 1 7 . E x p e r i m e n t a l p r o c e d u r e i n t h e s m o k e - w i r e t e c h n i q u e .

The objective of the videos was to identify the frontier line between the bulk flow with a high velocity level and the reversed flow (recirculation) with a low velocity level (see Figure 16). This frontier line is a function of several parameters. Some of those parameters are analyzed: the velocity level, the ratio between the inlet and the depth of the sub-rack, and the screen porosity. The frontier line is identified by five points (see

(37)

component of the velocity equal to zero. In order to making it easier to see the frontier point, a black paper was placed on the PCB dummy. Another technique based on tufts (see Figure 18) was used, giving results very similar to CFD predictions. However, the uncertainty was difficult to estimate (Castiella, 2005), and the technique could only be used for high velocities.

F i g u r e 1 8 : F l o w v i s u a l i z a t i o n u s i n g t u f t s .

2 . 6 . 2 U n c e r t a i n t y A n a l y s i s

First it was studied whether the position of the camera influenced the measured frontier points. After 60 videos it was concluded that this factor was of little importance due to the fact that the wires are close to the PCB. Actually, the uncertainty of the camera position was estimated to be up to about 0.5 cm (twice the standard deviation of the frontier point values observed in the videos).

Furthermore, in order to decrease the uncertainty, each wire was recorded between two and five times. The uncertainty of wires 1 to 4 (and depth equal to 200 mm) in the observed values is around 0.5 cm (twice the standard deviation of the observed values of the five videos that were recorded per wire), for the wire 5; at the outlet and for the large case (d=260 mm), the estimated uncertainty increases up to 0.5–1 cm. The estimated total uncertainty for the first four wires and the short model (depth equal to 200 mm) is 0.71 cm and for the fifth wire and the large model is around 1 cm.

(38)

3 Computational Fluid

Dynamics

3 . 1 I n t r o d u c t i o n t o t h e C F D M e t h o d

Computational fluid dynamics (CFD) is a way to calculate the fluid and heat flow in any geometrical domain by numerical methods. It can be called a whole-field method because it provides all the information in the whole numerical domain. The numerical domain is the result of discretizing the continuum domain in a set of volumes or points.

Nowadays, CFD is a method that has vastly increased its use in industry due to the increase of computational capability in terms of memory and processors. Many papers that use CFD predictions for both chip level and system level of electronic equipment can be found. However, there is still a need to check the validity of the numerical models, turbulence models and mesh density distributions by experimental measurements. Most CFD software packages have three elements: the preprocessor, the solver and the postprocessor. The preprocessor is the software that is used to generate the geometry and the mesh (the discretized domain) and finally to set the boundary conditions. Then the solver applies numerical methods in order to convert the partial differential equations into a set of algebraic equations applied to the nodes of the mesh that can be solved by an iterative process. The third element is the postprocessor, which provides the ability to visualize the results. It also provides the tools to calculate e.g. mean values for a specified sub-domain.

Fluid flow is governed by conservation laws: mass and momentum, (for non-isothermal flow the energy equation is added) together with some relations between fluid properties (e.g. the ideal gas law) and other relations: e.g. a linear relation between stress and strain (Newtonian fluids). The momentum equations are also known as the Navier-Stokes equations. The problem is that these equations are nonlinear partial differential equations that need to be solved by numerical methods. The

(39)

equation for steady state incompressible flow can be written in the following form: 0 w w i i x u ₍₂₎

i i j i j u x P x u u ₂ w w w w Q (3)

T x T u j j 2 w w D (4)

3 . 2 T u r b u l e n c e

At high Reynolds numbers, instabilities appear in the laminar flow that ends up in the turbulent regime. Turbulence is characterized by three-dimensional chaotic flow with a high mixing capacity, thus increasing the momentum and heat exchange in the flow. It is also characterized by a huge spectrum of length, velocity and time scales (or eddies). The large scales are the energy-bearing eddies; they extract the turbulent kinetic energy from the main flow (due to main velocity gradients), and correlate well with the main flow and are thus anisotropic. Its characteristic length scale is related to the geometry of the domain.

The large eddies form smaller eddies by a vortex-stretching mechanism. These smaller eddies receive the energy from the large energy-bearing eddies and transfer it successively to smaller and smaller eddies. This phenomenon, known as energy cascade, ends up in the smallest eddies that are dissipated into heat (internal energy) by the action of molecular viscosity. The length scale of the smallest eddies, also known as the Kolmogorov scale, depends on the value of the molecular viscosity and is considered isotropic.

The production of turbulent kinetic energy depends on the generation of the large energy-bearing eddies (due to the main velocity gradients). The rate of dissipation of the turbulent kinetic energy into internal energy depends to a certain extent on the production and thus also of the large eddies. If the rate of production is equal to the rate of dissipation, the flow is in a state of equilibrium.

When performing design work, what is interesting to know is not the instantaneous flow variables but, often, the mean quantities. Normally

(40)

this is what is measured when doing experimental investigations. This idea ends up in the Reynolds decomposition. This consists of decomposing the instantaneous flow variables (velocity, pressure, temperature…) into a mean value and a fluctuating value, as is indicated in equation 5. The mean value is a time-averaged value or an ensemble average. The mean value of the fluctuation terms is zero, but the mean value of a product of fluctuation terms is not zero.

t T t p P p u U ui i i c c c (5)

If we apply equation 5 to equations 1–3, the Reynolds average Navier-Stokes (RANS) equation appears.

In the present study, the flow is assumed to be steady-state, three-dimensional, incompressible and turbulent. The buoyancy effect and radiation heat transfer are also assumed to be negligible. With these assumptions, and based on the Reynolds decomposition, the governing equations for the conservation of mass, momentum and energy are given by 0 w w i i x U (6)

j i j i i j i j u u x U x P x U U c c w w w w w w ₂ Q (7)

t u x T x T U j j j j c c w w w w ₂ D (8)

where P is the modified kinematic pressure and the unknowns, u cicuj and u_ictc, constitute the second-moment statistical correlation or so-called Reynolds stresses and turbulent heat fluxes.

3 . 3 T u r b u l e n c e m o d e l l i n g

The instantaneous real flow contains a very large length and time scale range. Hence, it requires a very fine mesh and time step in order to solve the Navier-Stokes equations. Nowadays, this is possible for simple

(41)

computers and a large amount of computational time are required, something that the industry still cannot afford. To resolve directly (without turbulence modeling) the instantaneous Navier-Stokes equations is called Direct Numerical Simulation (DNS). Some advanced turbulence models called Large Eddy Simulation (LES) also resolve the instantaneous flow variables, but only for the large eddies. For the small eddies, modeling is employed to resolve the time-averaged values of the flow variables. LES is also computationally very time demanding and thus is seldom used in industry.

The RANS equations are not a closed set of equations; due to the nonlinear term of the original set of equations, a new term called the Reynolds stresses appears and thus new unknowns are added to the equation system. Actually, the Reynolds stresses are a symmetry tensor and thus, six new unknowns appear. These Reynolds stresses represent the increase in the momentum transfer due to the turbulence (correlations between fluctuating velocities). We have a closure problem due to a larger number of unknowns (13) than number of equations (5), and thus turbulence modeling is required. Turbulence modeling is the way to solve the closure problem, i.e., how to specify the Reynolds stresses and turbulent heat fluxes.

One way to solve the closure problem is to use the Reynolds Stress Model (RSM). The RSM, (see Launder et al., 1975), or the so-called second-momentum model, is one of the most sophisticated tools currently used by engineers to predict turbulent flows. The aim is to solve the transport equations for the transport of the Reynolds stresses,

j iu

u cc . The stress components are solved by a transport equation. A linear pressure-strain model available in Fluent 6.1 is used in the simulations. This model is not as computationally time demanding as the LES.

The most popular model to approximate the Reynolds stresses is the eddy-viscosity turbulence model. The eddy-viscosity model is based on the Boussinesq assumption that makes an analogy between turbulent stresses and viscous stresses (given by Newton’s friction law). Therefore, the Reynolds stresses are proportionally related to the velocity gradients (strain) by an eddy-viscosity model for turbulence, and the turbulent heat fluxes are modelled using the eddy-diffusivity hypothesis in analogy with the eddy-viscosity hypothesis for the Reynolds stresses.

k S u ui j Qt ij Gij 3 2 2 c c (9)

(42)

i t t i x T t u w w c c V Q ₍₁₀₎

where Qt is the eddy viscosity, k is turbulent kinetic energy, Vt is the

turbulent Prandtl number (which is constant in the present study and equal to 0.85) and S_ij 0.5(wU_i/wx_jwU_j/wx_i).

The kinematic eddy viscosity by dimensional analyses can be expressed as the product of a turbulent length and velocity scale (using an analogy with the kinetic theory of gases).

Several models have been implemented depending on the way to supply the velocity and length scale. These models range from the so-called zero-equation model to the two-equation models. In zero-equation models, the length scale is based on geometry and the velocity scale is calculated by an algebraic relationship. In the one-equation models, one transport equation is solved in order to calculate the velocity scale (normally, an equation for k is solved) and the length scale is calculated by an algebraic relationship.

Finally, in the two-equation models, the platform of transport equations may be closed by two additional transport equations of two parameters related with the turbulent length and velocity scale. This model implies that the length and velocity scales of the mean flow and of the turbulence are proportional, and can be related by means of dimensional reasoning to turbulent kinetic energy, k, and one length-scale-related parameter like its dissipation rate, H, u=k0.5_{, l=k}1.5_H_{. Considering the}

above assumption, the turbulent eddy viscosity can be derived as,

vt=CPk2H, and it is valid only when local isotropy in the turbulence field

is assumed. H is given by j i j i x u x u w c w w c w Q H (11)

Thus, the two-equation eddy-viscosity model requires two additional transport equations for k and Hor another scale-related parameter like time scale Z Hk) to solve the spatial and temporal variation of the local velocity scale and the length scale.