A 2D Conduction-Convection Model of a Jet Impinging on a Flat Plate by MATTHEW ROBERT BOLTON MMath., The University of East Anglia, 1998 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES (Department of Mathematics; Institute of Applied Mathematics; Applied Mathematics) WTe accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA May 2002 ©Matthew Robert Bolton, 2002 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for refer ence and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mathematics The University of British Columbia Vancouver, Canada Abstract A two dimensional model of a jet impinging on a heated wafer is constructed. After neglecting buoyancy, a Navier-Stokes solver provides the flow velocity such that the energy equation can be used to solve for the temperature in the gas. Using a coupling condition, a heat conduction solver in the wafer is coupled to the energy equation solver in the gas, and as such the heat transfer across the surface of the wafer can be computed. A predictive formula is constructed for the surface heat transfer coefficient in terms of the Reynolds number and the ratio of the jet height to the nozzle width. ii Contents Abstract ii Table of Contents iiList of Tables v List of Figures v 1 Heat Transfer in a Thin Semi-Conductor Wafer 1 1.1 Introduction 1 2 Formulation of the Problem 2 3 Modeling 4 3.1 Geometry3.2 Governing Equations3.3 Boundary and Initial Conditions 5 3.4 Wafer Model 8 3.4.1 Governing Equations3.4.2 Discretisation of Domain 8 3.4.3 Discretising the Boundary Conditions 11 3.5 Gas Flow Model 14 3.5.1 Governing Equations 5 3.5.2 Navier Stokes Discretisation: Centered flux evaluations, Implicit Euler time advance . 16 3.5.3 Energy Equation Discretisation: Centered flux evaluations, Trapezoidal time advance 18 3.6 Coupling the Two Models 9 4 Results 20 5 Conclusion and Proposed Future Study 31 A Validation 36 A.l Conduction in Wafer 3A.l.l Test Case 1A.1.2 Test Case 2 8 A.1.3 Test Case 3 9 A.2 Energy Equation 42 A.2.1 Test Case 1B Data Plots 3 References 7 iii List of Tables 1 Approximate values of parameters used 6 2 Ranges of parameters used for obtaining results 20 3 Values of parameters used for example calculation 33 iv List of Figures 1 Chamber side view 1 2 Showerhead top down view (Showerhead is the same size as the wafer) 2 3 Surface Wafer temperature image for a gas flow of 0.2Zmm-1. The lighter the color, the cooler the surface (Actual temperature data is unavailable due to commercial sensitivity; moreover the actual temperature values are not important for the purposes of this discussion.) 3 4 Surface Wafer temperature image for a gas flow of Q.hlmin-1. The rings of jet nozzles create colder rings on the wafer that are just visible in the image 3 5 The full geometry and boundary conditions for a jet impinging on a flat plate 6 6 Boundary conditions on the wafer. We have replaced the coupling condition at the interface with a simple convection term; this makes it is easier to validate the conduction model .... 8 7 Boundary conditions for the gas flow domain 14 8 Boundary conditions on channel inlet9 Gas Vector Field Characteristics 21 10 The heat transfer coefficient for H = 17mm, d — 5mm, Red — 9.33 22 11 The temperature in the gas above the wafer with the jet inlet at the top left of the figure (H = 17mm, d = 5mm, Red = 9.33) 3 12 The temperature in the wafer with the jet impinging at the top left of the figure (H = 17mm, d = 5mm, Red = 9.33) 213 (a) Scaled HTC for Re = 9.33 5 14 (b) Scaled HTC for Re = 21.3315 (a) Error for Re = 9.33 7 16 (b) Error for Re = 11.9917 (c) Error for Re = 15.11 28 18 (d) Error for Re = 18.2219 (e) Error for Re = 21.33 9 20 Maximum absolute error in predictive formula 30 21 Overhead view of one jet from a row of impinging axisymmetric jets 31 22 Sketch of the flow due to two neighboring jets as a possible future investigation 32 23 (a) Composite heat transfer function for example showerhead calculation 34 24 (b) Steady state temperature on wafer surface for example showerhead calculation 34 25 The boundary conditions for the temperature problem 1 36 26 The boundary conditions for the temperature problem 2 8 27 The boundary conditions for the temperature problem 3 9 28 The boundary conditions for the energy problem 1 42 29 (a) aPk = -0.6940 + 0.02337i?e - 0.0003527i?e2 4 30 (b) Apk = exp (1.7406 + 0.08516fle - 0.0016294fle2)31 (a) c*oo = -1.1799 + 0.02651Jte - 0.0008423i?e2 45 32 (b) Aoo = exp (1.8876 - 0.04313iRe + 0.001608fle2)33 (a) a. = .89427 - .015997i?e + .47219 x 10~5i?e2 6 34 (b) A* = exp (-5.66657 4- 0.029919i?e - 0.0002652i?e2) 4v 1 Heat Transfer in a Thin Semi-Conductor Wafer 1.1 Introduction In the production of micro chips it is necessary to produce semi-conductor wafers in a carefully controlled manner, such that structural degradation of the material is kept to a minimum. This structural degradation has a direct impact on the performance of the micro chips. The wafer is subjected to a series of processes (For a more detailed description of the manufacturing of Silicon Wafers, see [7].) Anealing is the process by which a wafer is heated to a high temperature (RJ 1050°C7) so that dopants can be introduced to the wafer, and thermally activated diffusion can take place. Rapid Thermal Processing (RTP) is the term commonly used to describe these types of process. Most of the techniques involved in RTP originate from ideas that arose from experiments conducted in the field of materials science. The first devices to resemble RTP were developed in the 50s and 60s. These devices used a convection furnace to heat the wafers in batches. A furnace creates a heat transfer environment that is difficult to control, due to the convection currents occurring. Although advances in furnace design have enabled these kinds of RTP device to heat up quickly, the speed is still appreciably slower than an equivalent device using a radiative influx of heat. The first commercial RTP chamber to use radiation was developed in 1999; unfortunately it was not a commercial success. In 1994, scientists at a Vancouver based company, Vortek Industries (www.vortek.com) began to research RTP as a possible application for their high powered arc lamps. Vortek is world renowned as a leading manufacturer of extremely high powered lamps that operate between 50 and 300kW. They are now at the stage in which a prototype RTP chamber has been built and adjustments are being made in order to refine their patented RTP tool. The key advantages to this process are: • The wafer can be heated to the peak temperature (w 1050°C7) very rapidly. This lowers the thermal budget (integral of temperature with time) thus makes for a more efficient production line. • The temperature across the wafer can be accurately controlled and thus superior wafers can be pro duced. It is desirable to maintain a uniform temperature throughout the wafer such that structural abnormalities and surface degradation is kept to a minimum, and thus the tracks which are to be placed on the wafer can be placed in closer proximity. There are many reasons to produce smaller micro chips, the most important being increased speed, reduced power consumption, and the advent of nanotechnology. A target temperature uniformity has been set at ±2°K due to market demand. Figure 1: Chamber side view The configuration of the equipment is depicted in fig. 1. An array of jets are used to cool the wafer with a flow of inert gas on the top side. The wafer must be held in a horizontal configuration for heating, as internal stresses within the thin wafer prevent a vertical configuration (which would perhaps be better suited 1 to convection cooling). The horizontal configuration also facilitates a smooth production production line, as machinery to hold and turn the wafer is not necessary and the wafer can simply sit in a recess. Figure 2: Showerhead top down view (Showerhead is the same size as the wafer) The wafer is housed in a black chamber that serves to absorb unwanted reflected radiation, and the light shines through a quartz window onto the bottom surface of the wafer. The space in the chamber is occupied with inert gas. The CCD camera uses the emitted radiation to measure temperature across the wafer surface and thus experimental data can be obtained. The top surface of the plate is subjected to a flow of inert gas that passes over the wafer via a showerhead that sits above the wafer. Fig 2 depicts a schematic diagram of this showerhead; note that the arrangement of holes shown describes an approximate layout for a typical showerhead configuration. Currently, much experimentation is underway concerning the placement of the holes. This showerhead can be adjusted from lcm to 2.5cm above the surface of the wafer. The wafer is about 15cm in diameter. When a wafer is in position and ready to be treated, the influx from the lamp is turned on. The influx from the lamp takes a few seconds to reach its maximum value, and as such the temperature of the wafer 'ramps' up to the peak temperature (about 1050°C) in about 5 seconds. The wafer is held at the peak temperature for a time period of around 15s. The lamp is then switched off and the temperature of the wafer ramps down, the most dominant method of cooling being the radiative heat loss. 2 Formulation of the Problem Vortek would like to know how best to arrange the apparatus such that the heat transfer across the surface of the wafer remains as uniform as possible. Put most simply, for a given flow rate, how high does the shower head need to be and how close together should the holes in it be placed to get a heat transfer coefficient that is as uniform as possible. Setting up a full 3D computational model would be the best way to obtain an accurate answer, but such a model would take a long time to configure, and computations are likely to take a long time (perhaps 4-5 hours per run, so to obtain the necessary range of data would take months.) As usual when solving applied mathematics problems the best approach is usually to simplify the problem as much as possible; this will provide insight into how best to formulate a more complicated model. Proceeding in this manner, we reduce the problem to that of a single (2D) jet impinging on a flat plate, which is heated by a constant influx of heat. Since the problem has now been made 2D, the information obtained will be qualitative rather than quantitative. However, we should hope that the qualitative relationships obtained will still be applicable to an equivalent axisymmetric 3D jet, since the domain we consider is effectively a cross section of such a jet. In simplifying the model we have also removed the effects of neighboring jets, and this probably represents the greatest deviation from the actual physical representation. Appealing to the pictorial data obtained from Vortek, it would seem that a jet introduces a localized peak to the heat transfer coefficient on the surface of the wafer. It seems likely that interference between the jets and the corresponding effect on the heat transfer coefficient is confined to the regions between the jets. We are hopeful then, that our model of a single jet gives a useful qualitative description of the heat transfer coefficient in the vicinity 2 of a single jet. We will be unable to provide quantitative data (with regards to the actual setup) using this model, but we aim to gain a qualitative understanding of how changing the height of the showerhead (H) and the jet nozzle width (d) affect the local heat transfer. Some attempt is made, however, to use parameters that will give variable values in the right kind of range (e.g. the peak temperature on the wafer is kept to around 1050°C which is around the same temperature used in the annealing process ). For the purposes of qualitatively describing the local HTC (Heat transfer coefficient) with regard to unifor mity, the influx parameter was modeled using its peak value and disregarding the ramping stages. However, in order to obtain quantitative data when implementing a 3D model, these effects would have to be included. Figure 3: Surface Wafer temperature image for a gas flow of 0.2Zmm-1. The lighter the color, the cooler the surface (Actual temperature data is unavailable due to commercial sensitivity; moreover the actual temperature values are not important for the purposes of this discussion.) Figure 4: Surface Wafer temperature image for a gas flow of 0.5/mm 1. The rings of jet nozzles create colder rings on the wafer that are just visible in the image. Note: The pictures presented above in no way reflect the quality of Silicon wafers produced by Vortek. The examples shown here are merely plots obtained during testing and are used here to illustrate how the impinging jets create changes in the heat transfer. The heat transfer patterns from the latest Vortek model are far superior to the ones presented here. 3 3 Modeling 3.1 Geometry As discussed before, the problem has been reduced that of a 2D jet impinging on a flat plate. It is hoped that the results obtained will be useful in describing (qualitatively at least) the heat transfer due to an equivalent 3D axially symmetric jet. Please refer to Table 1 for a description of parameters and their values. 3.2 Governing Equations The temperature T(x, y) in the wafer is modeled by the 2D heat conduction equation: The gas flow is modeled using the non-dimensional Navier-Stokes equations in artificial compressibility form. Artificial compressibility means that we have added an aphysical time derivative of pressure to the continuity equation. This derivative is scaled by a parameter /3 that sets the pseudo-compressibility of the fluid. With the addition of this time derivative of pressure, we are able to advance pressure and velocity in time together. Note that this approach can sometimes affect the time accuracy; care must be taken in computing transient processes. (Since we are computing steady state solutions this will not be of great concern.) The aim in the gas flow domain is to compute the flow field using the Navier-Stokes equations, and then to use the Energy equation to compute the temperature. In decoupling the velocity and pressure from the temperature we have neglected the effects due to buoyancy that may be of concern near the heated wafer, where large temperature gradients may be possible, especially during the initial (transient) stages of heating. Buoyancy is caused by the fact that hot air is less dense and thus a pressure difference exists between the hot and cold gas, the latter tending to sink under gravity, the former rising upwards. The jets' function is to blow away these rising and falling gas currents so that a uniform temperature can be established on the wafer surface. We hope then that in the actual prototype RTP device, the speed of the gas next to the wafer is large enough to blow away these unstable gas currents. For the purposes of our model we are convinced that in the vicinity of the jet any pressure gradients due to temperature differences will be far outweighed by the pressure gradient caused by the impinging gas. As we move away from the impinging zone, and the situation can be likened to that of a uniform flow past a heated plate, we may have to pay closer attention to the effects of buoyancy. Indeed if the flow is slow enough, it may be possible that convection cells are created in the form of a rolling vortex made up of a flow loop of rising and then falling gas close to the plate surface. This is the so called Poiseuille-Benard flow which is studied in the 2D case by [5]. In practice, it may be possible for the regions between the gas jets to produce flow fields in which convection cells can exist. However, since our model's physical significance is limited to the region near a single jet, it will not be instructive to investigate convection cells in the context of this model. We leave this area of investigation for future study. Here then, are the non-dimensional Navier-Stokes equations in artificial compressibility form: dP _1_ (cPu cPu dx + Re \dx 2 + dy 2 ) (b + H) d 4 The equations have been non-dimensionalised via the following: P-JL T-H X~d' U~U' pU*' T* - y ~ v y=d> v = u-where Ud v 2 TjjVrnazI is the mean velocity across the inflow, 1°K is a typical temperature fluctuation Hats were subsequently dropped, and the domain ranges become: L H 0<x<-, 0<y<—. In order to solve for the temperature in the gas flow domain, we use the flow field given by the solution of the Navier-Stokes equation to compute the energy equation for incompressible flow (i.e. the flow velocity v satisfies V.tf = 0, at least at steady state. This assumption is valid for the low Reynolds numbers encountered in this model). Re = U = T* = In which _ Ud _ inertia ^ _ fxcp _ dissipation ^ _ cpT* _ enthalpy v viscosity' k conduction' U2 kinetic energy 3.3 Boundary and Initial Conditions l.(a) Vertical Direction - Wafer A picture of the geometry, including the boundary conditions imposed is shown in fig. 5. At the top surface of the wafer we have the coupling condition corresponding to temperature continuity through the wafer-gas interface. We have also included the radiative heat loss term which is represented by the fourth power relationship on T. For a better explanation of how this term is derived see [3]. Basically, radiative heat loss occurs due to energy being released from the material in the form of an electromagnetic wave, usually with a frequency towards the red end of the electromagnetic spectrum. This is also why objects can appear to glow red when they are hot. This heat is radiated to the chamber walls which absorb the heat energy. The walls are specially painted black and the material carefully selected such that the reflected radiation is kept to a minimum. For more information on radiative heat transfer in the context of RTP see [7]. We hope that the heat radiated to the chamber walls is quickly dissipated by convection, and that our assumption of room temperature being maintained at the wall throughout is a realistic one. If an accurate quantitative model were being constructed, this effect may have to be reconsidered as it may well be expected that the chamber walls will heat up to some degree during the course of operation. We observe that since the radiation heat loss term is a fourth power relationship on T, it will be the dominant mechanism of cooling once the temperature is sufficiently large. At the bottom surface 5 Tx = 0 ux = v = 0 P = 0 T* = O u = v = 0 P„ = 0 ^ = fear, - oE (T4 - 1%, TT = 0 kSiTy = -a(t) + aE(Ti-T4Room) Figure 5: The full geometry and boundary conditions for a jet impinging on a flat plate Parameter Value Unit Wafer Thermal Conductivity (ksi) 149.0 Wm^K'1 Gas Thermal Conductivity (koas) 0.025 Wm~lK-1 Wafer Thermal Diffusivity (ust = ^) 9.0 x 10-5 Wafer Stefan Boltzman const, (a) 5.67 x 10"8 Wm-2K~4 Wafer Emissivity (E) 0.7 -Influx (a) 2.33 x 105 Wm-2 Room Temp. (TRoom) 297.0 K Wafer Length(L) 0.10 m Wafer Width (6) 0.0007 m Typical Shower Head Height (H) 0.02 m Typical Peak Incoming Velocity (|Vmo:i:|) 0.042 ms"1 Typical Peak Outgoing Velocity (Vout) 0.0105 ms-1 Typical velocity scale (U = %\Vmax\) 0.028 ms-1 Typical Reynolds No. (Re = ^) 9.33 -Table 1: Approximate values of parameters used we have a (constant) influx of radiant heat provided by the lamp, and an outflux of radiant heat. Since the gas is prevented from circulating around the wafer, and the flow field near the bottom of the plate is dynamically stable (i.e. the flow will be carried along the plate and then upwards by natural convection) , convection effects at the underside were ignored. Again, if we were looking for a quantitative representation, this may have to be reconsidered. (b) Vertical Direction - Gas Pressure, velocity and Temperature are prescribed on all sides of the boundary. At the Wafer-Gas interface we have a coupling condition based on temperature continuity. At the top of the domain (showerhead surface) we have prescribed the temperature to be at room temperature. As described before this neglects the radiant heat exchanged between the showerhead and the wafer surface, see [7] where this process has been studied. We do, however, include the radiation heat loss to the surroundings. 2.(a) Horizontal Direction - Wafer Since the wafer is so thin, the heat loss through the edges will be minimal, and thus we can use a no flux boundary condition. 6 (b) Horizontal Direction - Gas At the left hand side of the domain we have a symmetry plane, hence there is no flow across the boundary (u = 0) and the other quantities have x derivatives zero. At the outlet we assume the flow is fully developed and take the pressure to be zero (we are free to measure pressure from whatever reference value we choose.) 3. (a) Initial Condition - Wafer The initial temperature distribution in the wafer is taken to be constant (room temperature). 3.(b) Initial Condition - Gas The velocity, and pressure are taken to be zero everywhere. The temperature is taken to be room temperature, as in the wafer. Note that convergence times could probably be reduced by using an initial velocity field that was more similar to an impinging jet e.g. using data from computation with similar parameters. This method was not utilized however, to ensure the numerical solution was correctly converged. With the conditions as stated we have a well posed problem, and thus the solution to the problem is unique. 7 3.4 Wafer Model To simplify validation of the model, we break the model into two parts, that of the conduction in the wafer and the gas flow model. To separate the model we need to specify a new boundary condition at the wafer-gas interface. The most obvious choice is to use a linear outflux which corresponds to amalgamating the heat conduction (Newton's law of cooling) and convection effects into one term by a judicious choice of constant of proportionality, the surface heat transfer coefficient (h). For the purposes of validating this part of the model, we can take h as some constant. Later, we will use this boundary condition in comparison with the full coupled boundary condition to determine the heat transfer coefficient for the full problem. 3.4.1 Governing Equations The temperature in the wafer is modeled by the 2D heat conduction equation: dT fd2T d2T\ n _ n 0<X<L' 0<y<b-kSiTy = -h(T- TRoom) - aE(T* ks,Tv = -a(t) + ffE(T*-T4Room) Figure 6: Boundary conditions on the wafer. We have replaced the coupling condition at the interface with a simple convection term; this makes it is easier to validate the conduction model 3.4.2 Discretisation of Domain In order to solve the problem numerically we must divide the space into control volumes by forming a grid over the domain with spacings Ax and Ay in the x and y directions respectively. In this analysis, each cell is labeled by a subscript which refers to the ith column and jth. row. A bar over the variable denotes the average value. It can be shown (see [2]) that the central value of the cell is second order accurate to the average value over the cell. In this paper we shall simply refer to the variable as, say Tij and implement this as the central value of the (i,j)th cell, but remember that it actually means the average value of the (i, j)th cell. We must also divide the time into steps of size At. Some analysis of the combined space-time discretisation scheme is necessary to ensure that the time step used is compatible with the space step sizes. Space A second order centered flux was used to discretise the Laplacian: 1 dT _ (Tj+ij — 2Titj + Ti-ij Tjj+i — 2Tjj + Tj-i \ ~dti~K V AP Ay2 )' Using Von Neumann analysis (see [2]) it can be shown that the eigenvalues are 2K 2K Xkl = -7r-2(l~C0S<Pl') ~ -r-2 (1-COS0/). For more information on discretisation, including a proof that this centered discretisation is second order accurate see [2] 8 Time A second order centered implicit (Trapezoidal) time advance scheme was employed. It was important that the scheme be at least 2nd order since we have a transient process with time dependent boundary conditions, although in this paper we present results only for steady state solutions we anticipate that it may be use ful to have a time accurate model. Also, probably more importantly we must not damage the space accuracy.) The amplification factor G of a solution is defined as the growth rate of the solution from time level t — nAt to time level t + At = (n + l)At. As an example, consider dv — = Aw. dt The exact solution is v{t) = Aext, and thus the amplification factor is given by v(t) vn 2 6 By comparing a numerical schemes amplification factor with this Gexact we can determine a scheme's time accuracy. \G\ < 1 gives the condition for stability. Thus when A is expressed as a function of the space stepping Arc, Ay we can examine the space/time discretisation together. The trapezoidal scheme is advantageous because of its stability properties: The amplification factor is ^ , , A (AA*)2 (AA*)3 G = 1 + XAt + v ; + v 1 , 11 -I- ^1 |G|= 2 1 XAll • \x 2~l Thus we can see by comparing with the exact amplification factor that the scheme is 2nd order accurate. Now consider A : A = A + IB, , . A2 B2 A2 B2 l + A+ — + —<l-A+ — + — 4 4 4 4 => A < 0 for stability. This means that there is no (theoretical) restriction on the time step in order for the scheme to be stable because the eigenvalues for the space discretisation lie in the stable space of the time advance scheme. Applying this time advance scheme for each point in the space discretisation (3.4.2) yields, in matrix form, the following system of equations: At + [S*] 2 + [5»] 2 = °' where 9 T21 T = Ti-ij Ti,j Ti+l,j and [S,] = Ax2 (2 -1 -1 2 -1 -12 1 and similarly [Sy] Ay2 V (2 -1 -1 2 -1 -1 2 1 V \ / \ J Hence we obtain [I] + ^ + ^) sf = ([Sx] + [Sy]) f». = -At(FluxIntegral)ij Approximate Factorisation In order to solve the resulting matrix problem, broadly speaking we need to separate the problem into its x and y solution components. Then we may apply boundary conditions in one direction and solve that sub-problem, and then apply boundary conditions in the other direction to fully solve the problem at that time step. One way of delivering the matrices into the desired form is the method of approximate factorisation (for more information see [2]). This is demonstrated below; the matrix equation is written as a product. It may be checked that when this product is expanded out, the extra term is of order (At)2 and thus will not affect the time accuracy. Applying the method of approximate factorisation gives This introduces an error of 0(At3) and does not affect the time accuracy which is 0(At2). Specifically, if we denote the operator Ea$ by We obtain the final result nAt 1 -2Aa;2 (Eifi —2 + E-ito) Ea,bSTitj — 6Tj+aj+b. 1 — ——^ (-^O,! — 2 + EQ^-I) 2Ay2 STij = At{Flux Integral) Ij 10 In practise this equation can now be solved in two steps (we introduce 6Q{j as an intermediate matrix): K,At 1 1 - ^XZt - 2 + £-i>°) 6®i,j = ^{Fluxlntegral)^, 2Ax2 which gives the information necessary to update the temperature for each time step. 3.4.3 Discretising the Boundary Conditions A simple way to prescribe boundary conditions on a mesh is to use cells on the exterior side of the boundary, so called 'ghost' cells. Thus updating the value of the ghost cell based on a calculation involving the interior cells provides a quick method of imposing the boundary condition. This means that the same flux calculation can be used everywhere, which is efficient and can be easily implemented. Bottom of Wafer At the bottom of the wafer we have the radiative influx of heat provided by the lamp, and the heat loss radiated out to the chamber walls. We need to find an equation that will determine how to set the ghost cells for both the space and the time discretisations. Space Restating the boundary condition in continuous form:2 BT On y = 0, k— = -a(t)+R{T*-TRoom), where T(x,y,t) is the temperature, Tw{x, y, t) is the temperature on the bottom surface of the wafer, TRoom is the (constant) room temperature. R is the product oE the radiation parameter, used for brevity. In discrete form this becomes =• T,-_! = Tj + ^ (a(i) - R(Tt - TRoomj) . Implicit Time Advance In discrete form we have (using the notation Tn to represent T at time level n): k Ay = -an+R 1 Room If we denote ST — Tn+1 - Tn, writing the above equation at time level (n +1) and then subtracting gives 2the influx term a(t) is shown as a function of t to reflect the fact that the model would accommodate such a parameter. In computing the results for this study however, its peak value was used. (Also note that a(x,t) can also be accommodated.) 11 k Ay^ 0 = ~Q"+1 +an+R t((T-+1)2 + CO'JCC"1 + TZ){TZ+l - O] . The problem here is that the temperature at the bottom surface of the wafer (Tw) is unknown at time level (n + 1), but we can approximate it by Tw at time level n (a more detailed analysis of the error introduced here is given further down the page) « an - an+1 + 2R(T^)2(2T^)STW, 5Tj (k - 2R(TZ)3Ay) = ST^ (k + 2R(T^)3Ay) - (an+1 - an)Ay, where W ~ 2 Error Analysis of Approximation used in Radiation Boundary Condition (Note that this analysis is only pertinent to transient computations, the results obtained during this study were all obtained for steady state solutions. The analysis included here is to show that the model is capable of time accurate results should the model be required to simulate transient cases.) Consider the Taylor series expansion for T™+1 dTn T£+1 = + -^-At + 0(At2), (Tn+1)2 = (Tn)2 + 2T^At + 0(At2), (r„+l)2 + (rn)2 = 2(Tn)2 + 2T^At. so to second order dTn dTn (2(T£)2 + 2(r;)2^At)(2T^ + -jjS-At) ~ «T£+1)2 + (T^)2)(T^+l + T") dTn =• 4(T£)3 + 6(T%)2^At ~ (K+1)2 + (T:)2)(T^+1 + T"). (1) So in approximating the radiation boundary condition we have introduced an error of size |6(T")2^"-At|. Looking at the left hand side of (1), we must impose the condition dTn The error is largest when the temperature at the wall is varying quickly. As the temperature approaches steady state, the error will decrease to zero. This suggests that it may be necessary to take smaller time steps in certain regions where the temperature is varying rapidly. Refining the time step should provide a simple method to determine whether the error introduced is important. The analysis performed here gives a good indication of when this error is likely to be largest, and therefore over which time period the time convergence should be verified. Top of Wafer At the top of the wafer we have radiation heat loss as well as the linear outflux corresponding to the convection cooling provided by the impinging jet. The parameter h is in reality a function of x the distance along the plate. For the purposes of validating this model we take ft to be a constant. We need to find an equation that will determine how to set the ghost cells for both the space and the time discretisations. 12 Space In continuous form the boundary condition is dT k— = -h(x)(Tw - TRoom) - R(T* - TRoom), where T{x, y, i) is the temperature, Tw(x, y, t) is the temperature on the top surface of the wafer, TRoom. is the (constant) room temperature. R is the product aE the radiation parameter, used for brevity. Writing this in discrete form we obtain k (^ir1) = ~h (p^1) ~ TR°°m) ~R{T™ - T«°°™}' 1j — (T,--! (k - + L\yhTRoom - AyR(T* - TRoom)) (k + ^) Implicit Time Advance Using the notation T™ meaning Tw at time level n, and 5T = Tn+1 — Tn , we have: STj-STj-A^ u fSTj+STj^ Anfmn,3 fSTj-i+STj h^ , 1 T>ITn\i A.\ _ XT. (l, h^ _0OfTn\3, 5Tj[k + —±+ 2R(TZYAy = 8T^ [k - —^ - 2R(T^)iAy 13 3.5 Gas Flow Model We now consider the gas flow above the wafer. Fig. 7 below depicts a 2D impinging jet on a plane surface. Using a symmetry boundary condition along the jet axis, only half the jet needs to be considered. The gas flowing in is assumed to be at room temperature, and has a peak velocity |Vmai|. Since the inlet channel is narrow we may assume that the flow is fully developed in the channel and emerges as a quadratic velocity profile. This profile can be used to compute the pressure gradient at the entrance to the gas flow region. u = 0 Vln(x) = (x2-d2)j2\Vmax\ u = v = 0 Pv = 0 UX = V-P = 0 Tx = 0 -Ty = kSlTy - CJE T4 - Tj Figure 7: Boundary conditions for the gas flow domain Channel Inlet: Fully developed profile X Symmetry P Jet Inlet x ane V (x) = (x2 - d) j2\Vm Figure 8: Boundary conditions on channel inlet For a fully developed profile one may assume the following velocity in the y direction (see [8] or [1]): 14 V (x) = (x2 - d2) ^\Vmax\, and from the y momentum equation dP (32V\ 2 The mean velocity across the inflow can be related to the maximum of the velocity profile as follows: »=y; Vin(x)dx d3 rd Jo d2)dx 3.5.1 Governing Equations Under the assumption that the Reynolds number of the flow is large enough such that buoyancy effects can be ignored, temperature can be decoupled from velocity. Henceforth the problem is modeled by first solving the Navier-Stokes equations yielding the flow field at steady state and then solving the energy equation. As discussed in section 3, the non-dimensional Navier-Stokes equations in artificial compressibility form can be written as dP ldu 1 dv dt + /3dx + fidy ~ ' du du2 duv dt dx dy du duv dt dx df_ dy dP_ J_ dx Re dP_ J_ dy Re d2u d2u dx2 dy2 d2v d2v dx2 dy2 The incompressible energy equation is dT dvT dvT dt dx dy dvT _ 1 fd2T + Ec Re - Pr \ dx2 dy Re du dx + 2 'dv du\' ^dx dy J The equations have been non-dimensionalised via the following: P x X=d' « V V=d' u p = pU2 T= — T T* Thus _ Ud inertia Re = — = —, v viscosity Pr = pcp _ dissipation k conduction' Ec = cpT* enthalpy U2 kinetic energy' 15 where 2 U = -\Vmax\ is the mean velocity across the inflow, o T* = \°K is a typical temperature fluctuation Hats were subsequently dropped, and the domain ranges become: L „ H 0 < z < -, 0<y <—. a a 3.5.2 Navier Stokes Discretisation: Centered flux evaluations, Implicit Euler time advance In matrix form we have dU 9F^ + dG_Q dt dx dy ' where U =\u \ v i ' 0 \ G = uv — V J^du Re dy Re dy j Integrating over a computational cell, applying Gauss's theorem and dividing by the cell volume gives dt Ax Using centered flux evaluations and a functional notation Ay '+5J F (Uij, Ui+ij) I Ui,,-+Ui + l,j \ I 2/3 1 ("i.j+Hi + l.j ^ , Pi.j+fi + l.j 1 Ui + l.j-Uj,, \ 2 ) "•" 2 fle Aft 2 2 fle Ax and : G {Uij,Ui,j+i) ( 2/3 ui,j+ui,J + l vi,i+vi,i + l 1 Ui,3 + 1-Ui,j 2 2 fle Ati 2 1 / VjJ + Vjj+A Pi,J+Pi,i + l 1 Vi,i + i-Vi,i J \V 2 y 2 Re Ay 'J Using an implicit Euler time advance, the fully-discrete formulation is ui,j ui,j _ OUi,j _ »+;.J «-;,J At At Ax Aw In order to evaluate the fluxes at time level (n + 1), a Taylor series expansion is used to express these fluxes in terms of data at time level n. F^.^FOJlfW^) = F(UPtj+SUiJ,U^1<j + SUi+ld) 6F"k . 5Fn.l . , . SUi,j 16 The derivative terms are flux Jacobians: dUi,j dF"1 . •+5,3 at/, 8Gn., ! x 2/3 II.. j+Ui + l,j I in <+vi+i,i ~ 4 1 —o ^ 'J Aa: Ke 0 iJij+Biiil. + 0 0 "<,J+"i,j + l 4 Ay Re Ay Re 0 0 Aziie-1 0 0 Ui.j+tXj + l.j _ 1 4 As ile-I X 2/3 4 + 2/3 "•,i+"»,J+i 4 fi,J+^i,J + l _ 1 2 Aj/fle If we substitute the expanded fluxes into the fully-discrete equation, we obtain: I 1 9FP^ j dF?_k 1 dG" L 1 AG" A • t+2 X 1 2 ^ I x '..7+2 ___1£ZJL I A/7 At Aa; 3(7;,Aa; 917^ Ay at/*,,- Ay dUitj I lJ Ax at/j+ij Ax oUi-ij Ay oUij+i 'J Ay dUij-i TTn Z?n /^*n /^n Ax Ay Multiplying by At and relabeling terms we obtain {I + AtBx + AtBy) 5Uij + AtCx6Ui+ij + AtAx6Ui-itj + AtCySUitj+i + AAy6Uij-i Fn ! . — FJ1 t . G" , - Gn. l Ax Ay Approximate Factorisation and Block Tri-diagonal Solution Applying an approximate factorisation (see section 3.4.2) we get [I + AtAxE_ifi + AtBx + AtCxElfi]x F.™ i - -F-1i • G" A - Gn._L [I + AtAyE0^ + AtBy + AtCyEo^SUij = — *'\x ' ^ - \y * where EifiSUij = 5Ui+ij. As described at the end of section 3.4.2, this equation can be solved in two steps, using the Thomas algorithm to solve in each direction. 17 3.5.3 Energy Equation Discretisation: Centered flux evaluations, Trapezoidal time advance Applying Gauss's theorem over a control volume we arrive at: dt v H)Ain Ji'J~i Ay Re-Pr\dxi-$,jAx dyij_LAy Ec ((8u dv\ fdv du\ \iJ+? ' PrAxAy VV^x flW ' ' 8y, , + PrAxAy \\dx + flu/ "+ Vflw dx)V)iA_,_ Using the notation Xto mean the central value of the (i,j)th cell and defining Sij as Ecffdu dv\ ,(dv.du^ //«•• «-N 5^ = -sr((^:-^:)w+(^ + Ec //flu flu\ /flu du\ v 3 Pr \_ V Ax dy J " ' \dx ' flw, we can apply a Trapezoidal time advance scheme (centered implicit) to obtain fTi>U 1 ( RT AT i /fr<+1,J--2^iJ-+Ari_1J-\\ , W ™ ™ 1 fSTu+1-26Ti,j+6TiJ-1 2Ay V 'J /?e-PrV Aj/ _ J_ | „. .fn . *n 1 f~ 2^3 + A., I Vl'31i,j v*,3-i1i,j-l Aj iJg-Prl Ay So we can rewrite this in matrix form as (3) The right hand side of eqn. 3 is the flux integral at time level n. The left hand side can be approximately factorised as per section 3.4.2 to give [I] + Y \DX) + Y lDy]) & * (W + Y [DX]) + T [DY]) 5F And thus we can solve separately in each direction using the Thomas algorithm as described at the end of section 3.4.2. 18 3.6 Coupling the Two Models We are now in a position to put the two models together by coupling them at the interface using a temperature continuity condition. This gives us the boundary condition at the interface as seen in Fig. 5 : This equation is in dimensional form, so we must be a little careful in the implementation to ensure that the correct scaling is used as we move from the wafer domain (unsealed) to the gas domain (scaled). In practice this condition gives us a relationship to update the ghost cells at each time step, as follows: 1. Compute the heat conduction in the wafer using an initial condition in the wafer and gas domain of room temperature everywhere. This means the boundary condition for the top of the wafer is just a Dirichlet boundary condition, with the value as given by room temperature. 2. Once the conduction in the wafer is computed, we can compute the new surface wafer temperature using the old gas temperature data above the wafer (assuming temperature continuity across the interface we can just use T''J+2ri|'~' = TijWan). We can now use this wafer surface temperature as a boundary condition to advance the energy equation. 3. Once we know the new temperature distribution in the gas flow domain, we can use the coupling equation to obtain the new ghost cells such that the heat equation can be advanced. Steps 2-3 can now be repeated until the solution is within a tolerance value of steady state. For the purposes of this study, we only wish to know the heat transfer at steady state as this will enable us to extract the heat transfer coefficient by comparing with the boundary condition we assumed in section 3.4 for the heat loss through the wafer, 19 4 Results With the volume flux through the domain held constant, and the jet radius d fixed, results were obtained for a range of H. This was repeated for varying Reynolds numbers (Red = ^) corresponding to differing rates of flow. The values used for the parameters are tabulated below. A mesh size of 80 x 80 was found to give a reasonable accuracy / computation time trade off, the time taken for each test case was about Admins on a P4 1.8GHz system. Red U H/m d/m H/d 9.33 0.02800 0.017 0.005 3.4 9.77 0.02933 0.018 3.6 10.22 0.03066 0.019 3.8 10.66 0.03200 0.020 4.0 11.11 0.03333 0.021 4.2 11.55 0.03466 0.022 4.4 11.99 0.03600 0.023 4.6 12.44 0.03733 0.024 4.8 12.88 0.03866 0.025 5.0 13.33 0.04000 0.026 5.2 13.77 0.04133 0.027 5.4 14.22 0.04266 0.028 5.6 14.66 0.04400 0.029 5.8 15.11 0.04533 0.030 6.0 15.55 0.04666 0.031 6.2 15.99 0.04800 0.032 6.4 16.44 0.04933 0.033 6.6 16.88 0.05066 0.034 6.8 17.33 0.05200 0.035 7.0 17.77 0.05333 0.036 7.2 18.22 0.05466 0.037 7.4 18.66 0.05600 0.038 7.6 19.11 0.05733 0.039 7.8 19.55 0.05866 0.040 8.0 19.99 0.06000 20.44 0.06133 20.88 0.06266 21.33 0.06400 Table 2: Ranges of parameters used for obtaining results 20 The figure below shows the streamline pattern for a domain with H = 17mm, d = bmm, Red - 9.33. The Jet inlet is at the top left of the figure. The darker color shading corresponds to larger velocity magnitude. (Note that the vertical is exaggerated by ~ 5L.) 21 h \ • \ ' \ \ \ \ \ \ \ \ \ \ \ 0 0.2 0.4 0.6 0.8 1 X Figure 10: The heat transfer coefficient for H = 17mm, d = 5mm, Red = 9.33 The following figures (figs. 11 and 12) show how the temperature varies in the gas flow and wafer, for the same domain (H = 17mm, d = 5mm, Red = 9.33). However these figures probably do not represent the actual temperature observed experimentally, due to the limitations outlined in section 2. It may be noted that a change of lWm^1K~1 in h produces a temperature deflection of approximately 0.6° K. This estimate is not refined further however due to the qualitative nature of the study, i.e. we are really investigating the qualitative behaviour of h in order to asses its uniformity. 22 T/°K L x Figure 12: The temperature in the wafer with the jet impinging at the top left of the figure (H = 17mm, d = 5mm, Ra = 9.33) 23 We aim to relate the heat transfer coefficient (HTC) to the Reynolds number (Re) and the ratio (£). Firstly suitable scalings were chosen to represent the data: for the vertical scale (h - /loo) (hpk — hoc)' and for the horizontal scale x x7' where h is the heat transfer coefficient, /ioo is the value far from the jet (for each Re), hpk is the peak value (for each Re), x is the distance along the plate, (non-dimensionalised with respect to L, the plate length) x* is the x value where ) = \ (for each Re). Providing that we can find suitable relationships for the parameters that depend on Re, we are hopeful that a generalized formula for the HTC can be postulated. A sech2 (pr) curve seemed to be the most obvious choice of profile to fit the data, matching the curves at x = x* i.e. (h - /loo) _ ,,^2 (Ax' (hpk — /ioo) = sech2 (^^j with A = cosh"1 (v^) ^ 0.88137... The two plots below show the scaled HTC data along with the sech2 (^f) for Re = 9.33 and Re = 21.33. The fit is still pretty good for the larger Reynolds number. Larger values of H (corresponding to larger ^) give profiles that spread out a little further, the effect being more pronounced for larger Re, as one would expect. 24 (h - /ioo) 1 X* Figure 13: (a) Scaled HTC for Re = 9.33 (h - hoc) 1 H = 0.04 7 x 8 Figure 14: (b) Scaled HTC for Re = 21.33 25 In order to correlate these plots we need to find relationships for hpk,hoo, x* in terms of Re. By assuming relationships of the form hpk = APk (Re) aPk{Re) we can plot log hPk v apk {Re) log + log Apk, and thus obtain apk as the gradient, and log(Apfc) as the y intercept of the best fit line for each Re. This gives data for apk and Apk with Re so that plotting best fit lines to quadratic order through these curves gives us our relationships. This process can be repeated for hoo and x* and the resulting plots are given in appendix B. Putting all this together, we obtain the following formula for the (dimensional) HTC: where \x\ h = seen'* ( — I (hpu - hoc) + hoc, X* J X = cosh-1 (v^) ~ 0.88137..., Hjapk d HI hpk — Apk hoc — ^-oo x* — A* H with apk APK AQQ a* A, : -0.6940 + 0.02337.Re - 0.0003527Pe2, : exp (1.7406 + 0.08516Pe - 0.0016294i?e2), : -1.1799 + 0.02651fle - 0.0008423i?e2, : exp (1.8876 - 0.04313i?e + 0.001608.Re2), : 0.89427 - 0.0159963Ee + 0.4722 x 10_5i?e2, : exp (-5.66657 + 0.029919i?e - 0.0002652i?e2). In order to assess the accuracy of the formula in comparison to the computed HTC the difference between the computed value and the formulated value are plotted for increasing Re (figs. 15 to 19). In these plots x is the (dimensional) distance along the wafer (in meters). 26 Figure 15: (a) Error for Re = 9.33 Figure 16: (b) Error for Re = 11.99 27 8 0.1 Figure 17: (c) Error for Re = 15.11 Figure 19: (e) Error for Re = 21.33 29 If we now plot the maximum error (Emax) for each Re we obtain fig. 20 75 15 20 25 30 Re Figure 20: Maximum absolute error in predictive formula The least squares approximation to the curve is given by Emax = 0.2224 - 0.04273i?e + 0.002700.Re2. This relationship is useful for estimating the error introduced when extrapolating the formula to obtain results outside the range of the computed data. Looking at figs. 15 to 19 it appears that increasing ^ is not going to change the error behaviour substantially. For low Reynolds number, low values of ^ are perhaps prone to error when using the formulated result. 30 5 Conclusion and Proposed Future Study A 2D second order accurate computational model of the heat transfer processes occurring in a jet cooled Silicon wafer was constructed. The model was used to find a relationship between the surface heat transfer coefficient, the Reynolds number and the ratio of the jet height to the nozzle width. The formula was compared against the computed data and found to be in good agreement with the computational model. This provides a very useful relationship for calculating the heat transfer distribution. However, the validity of this result in the light of experimental data has not been investigated, which is perhaps the most significant failing of this paper. Most experimental data that was encountered catered only for high Reynolds numbers. Given more time for the study this is obviously something worth checking. Another point of contention that remains is whether the model created is actually of any use to engineers at Vortek. If we look at an enlarged overhead view of one jet from a row of jets (the neighbors having the same flow patterns), we may expect the following picture, fig. 21 Figure 21: Overhead view of one jet from a row of impinging axisymmetric jets In the vicinity of the jet it would seem that the flow radiating out from the jet in any particular direction is an analogue of the 2d flow described by the model. The velocity will be different, however, as the flow is spreading out radially. Further investigation of the flow field due to an impinging axisymmetric jet would be instructive here. The results obtained in this paper indicate that a 3D model using a similar coupling condition would be effective in describing the flow due to one or more impinging jets. Using a fully 3D model would allow experiments in jet configuration and should provide information on how best to arrange the jets on the shower head to obtain a uniform heat transfer coefficient. Results obtained in this study would be very useful when validating the model. Depending on how accurate the results turn out to be (after comparing with experimental and 3d numerical data) it may be useful to investigate interference between jets e.g. by modifying the boundary conditions to include, say, two neighboring jets. It is a little uncertain as to what the flow field would look like, and it is likely that the different ranges of Reynolds numbers will cause qualitatively different flow fields. As long as the gas flow can be computed, it should be a straight forward task to repeat the method used to postulate a predictive formula for the heat transfer due to two neighboring jets. This may provide much needed insight on how close the jets should be placed together in any particular configuration. In order to present some results that might be useful to the actual setup used at Vortek, an example showerhead was converted to an equivalent row of slot nozzles. The results are intended to be used to asses how the temperature is offset by a varying heat transfer coefficient. The example showerhead has holes placed in concentric rings across the 10cm disc. By taking a line drawn from the centre outwards, the rings were translated to slots (the radial distance is now just the distance along the 2D plate). The speed of the 31 Jet Inlet jet Met Outlet Figure 22: Sketch of the flow due to two neighboring jets as a possible future investigation gas through each hole was assumed to be the same, and this gives the same effective Reynolds number for each jet. Thus using the predictive model, the heat transfer coefficient may be formed via a juxtaposition of single (slot3) jets. This of course neglects effects due to cross flow, which would probably offset the data by a constant amount across the plate. The data pertinent to the example is tabulated in table 3 The jets on the example showerhead disc were placed as follows: 1 hole at the centre 8 holes at 0.5cm 12 holes at lew 16 holes at 2cm 32 holes at 4cm 72 holes at 8cm This configuration was converted to a slot nozzle configuration by having 1 slot jet at the centre 1 slot jet at 0.5cm 1 slot jet at 1cm 1 slot jet at 2cm 1 slot jet at 4cm 1 slot jet at 8cm The total volume flux through the example showerhead is QT = 0.5 l/min = 8.33 cm3^1. 3the word 'slot' is used to verify the fact that the jet is now 2D 32 Parameter Value Unit Wafer Thermal Conductivity (ksi) 149.0 Wm^K-1 Wafer Thermal Diffusivity (KSi = *£) Wafer Stefan Boltzman const, (a) 9.0 x 10"5 m^s-1 5.67 x 10-8 Wm^K-4 Wafer Emissivity (E) 0.7 -Influx (a) 2.33 x 105 Wm~2 Air kinematic viscosity (v) 1.5 x 10"5 i —i m s Room Temp. (TRoom) 297.0 K Wafer Length(L) 0.10 m Wafer Width (b) 0.0007 m Shower Head Height (H) 0.01 m Slot Jet radius (d) 0.00075 m Peak Jet Incoming Velocity (|Vmaa:|) 0.049 ms"1 Jet velocity scale (U = %\Vmax\) 0.0327 ms~1 Jet Reynolds No. (Re = ^) 3.27 -Table 3: Values of parameters used for example calculation The total area of holes is AT = 145 x 7T x 0.0752 cm2 » 2.56 cm2. Thus assuming the flow through each hole is the same, the average jet velocity is U = ~ =3.25 cms-1. AT This gives the jet Reynolds number ite = ^ = 1.63. v Based on this data (Re = 1.63, H = 0.01 m, d = 0.00075 m, thus f = 13.33), the model predicts that the heat transfer coefficient is h = .865395 sech2 (25.617484a;) + .323462. and the maximum error in this is predicted as Emax = ±0.165 Thus the composite heat transfer function is h = (.865395sech2 (25.617484a;) + .3120719659) + (.865395sech2(25.617484(a;- 0.005)) + .312072) + (.865395sech2(25.617484(a; - 0.01)) + .312072) + (.865395sech2(25.617484(a;- 0.02)) + .312072) + (.865395sech2(25.617484(a; - 0.04)) + .312072) + (.865395sech2(25.617484(x - 0.08)) + .312072). With an error4 Emax — ±0.99 Thus by computing the conduction in the wafer using the uncoupled wafer model (sec. 3.4), the temperature in the wafer can be computed at any time. For this example the temperature on the surface of the wafer at steady state along with the corresponding heat transfer function is plotted in figs. 23 and 24. 4This error is just a theoretical value formed from a sum of the errors in the composed heat transfer functions. The error due to the cross flow now introduced is likely to make this error estimate physically unrealistic. 33 Figure 24 shows the deflection in heat transfer coefficient predicted by the model for the arrangement of slot jets as described. It would be interesting here to compare the results found with the data obtained from experiment, unfortunately at the time of writing this information was not available to the author. Using the method described here, any arrangement of slot nozzles (and with differing Reynolds numbers) can be investigated. It is also worth noting that the spread of the jet (which is likely to be smaller for a round jet) could most easily be adjusted by raising the exponent of the sech function. There is also more tweaking that could be done in the light of experimental results e.g. changing the Reynolds number for each slot nozzle to better represent the flow. It is expected that neighboring jets will aid each other creating lower pressure areas where the holes in the showerhead are concentrated. 35 A Validation To ensure that the results obtained from the computations were correct, some exact solutions were formu lated and the L2Norm of the difference between the exact and computed meshes was calculated. We expect this L2Norm to decrease by a factor of 4 as the mesh is doubled in size for second order accuracy. Unless specified otherwise, the following data was used for the various physical parameters. The actual values of the parameters are not required to be exact (in fact the parameters will vary with temperature) but the values were chosen to be of the same order of magnitude so as best to validate the numerical scheme. The values used for the parameters are as tabulated in table 1. The temperature scale is in Kelvin. A.l Conduction in Wafer A.1.1 Test Case 1 [ This test case is designed to check the validity of the solution in the x direction.] In the interior of the rectangle we have dT _ (cPT 82T\ ~dt~K\dx2 + dy2 J Subject to the initial condition T {x, y, 0) = T0 (const) and the following boundary conditions y Figure 25: The boundary conditions for the temperature problem 1. Firstly we decompose T into the sum of a steady state component and a transient component. T = u (x, y) + v (x, y, t) With u (the steady component) satisfying Laplace's equation with the same boundary conditions as for T but with no initial condition (there is no t dependence, and v satisfying the full heat conduction equation with homogeneous boundary conditions and the initial condition v (x,y,0) = T0 — u (x,y) We obtain u {x>y) = i^W)cos (?)cosh ({Y -B) T) For v we try the following expansion 36 2^ Anm cos (—— J sin (o-„y) e ^ " > n=0 m=0 which satisfies the boundary conditions at the left, right and bottom of the domain. Then the boundary condition at y = 6 yields a cos (anb) = 0 => an = (2n + 1) 7T 2b n = 0,1,2. 2^ Anmcos(^-j—Jsml i ^—le V ^) n=0 m=0 ^ ' ied by 2 2 f" ,m , ^^ (m-KX\ . ((2n + 1) TlV The constants are determine bL r = LJ0 bJ0ToCOS{~L-)Sm{ 2b—)dydx-Lj0 bj0 «(*•*)«»(—)»»( SI )*dx dydx (2n+ 1)773/ 26 ® With Now And u {x'y) = co^ycos (?)cosh ((2/ -6) T ~ i iJl9n if m = 0, (J) _ J (2n+l)7r ' 10 if m ^ 0 fL /m-KX\ (2KX\ I C0S(—)C°S(—j 0 if m #2, \L )~\\ ifm=2 4L2 (2n + 1) ((2n + l)2 L2 + 1662) 7T And so 4T0 Sln , e-W)« 26 4L2(2n + l) /2TTX % [pn + l)2 L2 + 1662) 7T V A-(2n + l)7ry\ -*(.*+#) sm | — ) e y " > lb The steady state is reached very quickly, due to the high thermal conductivity and the very small thickness of the wafer. For the steady state case, the values for the L2 Norms were as follows: Mesh Size Z/2 Norm Factor 20 x 20 40 x 40 80 x 80 1.648e-006 4.129e-007 1.033e-007 3.991 3.997 For the transient case the L2 Norms obtained were Time Step Time Level Mesh Size L2 Norm Factor 0.0025 0.5s 20 x 20 0.1589 -0.00125 0.5s 40 x 40 0.03961 4.012 0.000625 0.5s 80 x 80 0.009896 4.003 37 A.1.2 Test Case 2 Again we are solving the heat conduction equation, this time with a prescribed flux on the bottom of the wafer, and Newton cooling (convection) at the top surface: dT %-{x,0,t)=a(x) dy dT — (x,b,t)=P(T-TRoom) Firstly, we can scale the temperature to simplify the convection term slightly by letting T — T — TRoom and dropping the hat, remembering to rescale the temperature at the end. dT — = a[x) dy Figure 26: The boundary conditions for the temperature problem 2. For steady state problem, the eigenfunctions are given by: /TWIN Xn = cos (-j—J YN = AN cosh (J^) + Bn sinh (^) The bottom boundary condition gives oo a w = Bn cos {—) 77 n—l So as long as a(x) is sufficiently smooth, we have 2 fL /nnx\ Bn = — / a fx) cos —— dx nir Jo \ L J On y — b: dT (Ansinh ("77") + Bncosh ) 77cos (~7~) = 13 Yl (Ancosh ("77") + Bnsinh ("77")) cos ("77") [sinh (2g^) /? - 2f C0Sh (n*J>)] [sinh (22*) 2i -/3cosh(2^)] 38 In particular if a (ar) = a cos {^jf-), for the steady state problem we obtain (after reseating the temperature by TRoom.) (/3sinh(^) -cosh(^) (sinh(^)^-^coSh(^)) cosh + sinh cos For the steady state case with /? = 2 x 10s, the values for the L2 Norms were as follows: Time Step Time Level Mesh Size Li Norm Factor 0.005 1.0s 20 x 20 1.795e-005 -0.0025 1.0s 40 x 40 4.499e-006 3.990 0.00125 1.0s 80 x 80 1.125e-006 3.999 A.1.3 Test Case 3 The transient solutions obtained so far have boundary conditions that are discontinuous across t = 0, and therefore very small time steps were required near t — 0 in order to get an accurate computation. To better validate the time accuracy, another solution was obtained that had zero influx at time t = 0 and then increased linearly. The boundary conditions are shown below. 7 a r dT Figure 27: The boundary conditions for the temperature problem 3 Since the boundary conditions are ID, (a(t) does not depend on x) the solution will not vary in the x direction. Letting T(y,t) = u(y,t) +v(y,t) we obtain two problems to solve for u(y,t) and v(y,t). du d2u m=Kdy^ 0<y<b u(y, 0) - TRoom un{0,t) =u„(6,i) =0 39 dv d2v -dt=KcV °<V<b v(y,0)=0 kvn(0,t) = -a{t) vn(b,t)=0 Clearly, with a constant initial temperature and insulated sides, u(y,t)=TRo For v(y,t) we try the expansion v(y,t) a0(t) OO + ^Tan(t) cos an(t) = lfo v(r1, t) cos (^p) dri dan{t) 2K fb d2v dt 2 f zv (mtr}\ 2K ~~b 2K ( ajfy ~k~ fb /mrr)\ n2n2 , J0 VC0S[^)-vr-dT> dan(t) fmr\2 If we suppose a(t) is linear dt • K(TTan{t) = ^a{t) an{t)=bkJ0 a(t) = ct 1 .A'(l-r) e 5* a(r)dr , 2KC _«n'.'(.-r) 2KC 2cfe kn2ir2 t + 27T2 (' e —^ r=0 JO KTL2N a0(<) ,9 OO v{y't)= 2bk +22a^)^{-r) n=l T = Tfioom + —+ V^-^ <+—^(e 1) cos(— n=l 40 For the transient solution, with a ramp up time of 5.0s, the following L2 Norms were obtained Time Step Time Level Mesh Size Li Norm Factor 0.01 1.0s 20 x 20 1.144e-005 -0.005 1.0s 40 x 40 2.874e-006 3.980 0.0025 1.0s 80x80 6.933e-007 4.145 41 A.2 Energy Equation A.2.1 Test Case 1 The boundary conditions are shown below. T = 0 u = 0 v = 0 T = 0 u = 0 v = 0 Figure 28: The boundary conditions for the energy problem 1 r2 exp (r\x + r2L) - r*i exp (r\L + r^x) r2 exp (r2L) — r\ exp (riL) where T(x, y) = sin (ny) wo , / «o , 2 ri'2=2^±V4^+7r and the velocity was simply specified by u0 = 1.0 v = 0 A Reynolds number of 666.666 was used, and the following L2Norms of the difference between the exact and computed solutions were obtained: Time Step Mesh Size L2 Norm Factor 0.005 0.0025 0.00125 20 x 20 40 x 40 80 x 80 1.738e-005 4.343e-006 1.0669e-006 4.0017 4.0710 42 B Data Plots The following plots show how the various parameter relationships were found. The trends for Apk, apjt and A*, a, show the best correlation (These parameters are related to the peak heat transfer and the spread in the heat transfer respectively) . The trends for Aoo, (related to the heat transfer far from the impingement point) are more sensitive to the convergence of the gas flow. In order to get a better trend here it may be necessary to converge the gas velocity next to the wall to within a smaller tolerance. 43 apk 10 12 14 16 18 20 Re Figure 29: (a) aPk = -0.6940 + 0.02337i?e - 0.0003527i?e2 InApk 10 12 18 20 Re Figure 30: (b) APk = exp (1.7406 + 0.08516i?e - 0.0016294i?e2) 44 10 12 14 16 18 20 Re Figure 31: (a) <*«, = -1.1799 + fJ.02651.Re - 0.0008423.Re2 10 12 14 16 18 20 Re Figure 32: (b) = exp (1.8876 - 0.04313i?e + 0.001608Pe2) 45 References [1] G.K. Batchelor. An introduction to Fluid Dynamics. Cambridge University Press, 1967. [2] J. Tannehill D. Anderson and R. Pletcher. Computational Fluid Mechanics and Heat Transfer. McGraw-Hill, 1984. [3] J. Hill and J. Dewynne. Heat Conduction. Blackwell Scientific Publications, 1987. [4] H. Martin. Heat and mass transfer between impinging gas jets and solid surfaces. Advances in Heat Transfer, 13:1-60, 1977. [5] A. Mojtabi and X. Nicolas. Two-dimensional numerical analysis of the poiseuille-benard flow in a rect angular channel heated from below. Phys. Fluids, 9(2):337-348, Feb 1997. [6] D. Phares R. Flagan and G. Smedley. The inviscid impingement of a jet with arbitrary velocity profile. Phys. Fluids, 12(8):2046-2055, 2000. [7] A. Wacher. A radiation model of a rapid thermal processing system. Master's thesis, University of British Columbia, 2001. [8] F. White. Viscous Fluid Flow. McGraw-Hill, 1974. 47
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A 2D conduction-convection model of a Jet impinging...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A 2D conduction-convection model of a Jet impinging on a flat plate Bolton, Matthew Robert 2002
pdf
Page Metadata
Item Metadata
Title | A 2D conduction-convection model of a Jet impinging on a flat plate |
Creator |
Bolton, Matthew Robert |
Date | 2002 |
Date Issued | 2009-08-20 |
Description | A two dimensional model of a jet impinging on a heated wafer is constructed. After neglecting buoyancy, a Navier-Stokes solver provides the flow velocity such that the energy equation can be used to solve for the temperature in the gas. Using a coupling condition, a heat conduction solver in the wafer is coupled to the energy equation solver in the gas, and as such the heat transfer across the surface of the wafer can be computed. A predictive formula is constructed for the surface heat transfer coefficient in terms of the Reynolds number and the ratio of the jet height to the nozzle width. |
Extent | 4336950 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | Eng |
Collection |
Retrospective Theses and Dissertations, 1919-2007 |
Series | UBC Retrospective Theses Digitization Project |
Date Available | 2009-08-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080028 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2002-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
URI | http://hdl.handle.net/2429/12402 |
Aggregated Source Repository | DSpace |
Download
- Media
- ubc_2002-0344.pdf [ 4.14MB ]
- Metadata
- JSON: 1.0080028.json
- JSON-LD: 1.0080028+ld.json
- RDF/XML (Pretty): 1.0080028.xml
- RDF/JSON: 1.0080028+rdf.json
- Turtle: 1.0080028+rdf-turtle.txt
- N-Triples: 1.0080028+rdf-ntriples.txt
- Original Record: 1.0080028 +original-record.json
- Full Text
- 1.0080028.txt
- Citation
- 1.0080028.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
United States | 18 | 0 |
France | 6 | 0 |
China | 5 | 8 |
Japan | 4 | 0 |
India | 3 | 1 |
Tanzania | 1 | 0 |
Germany | 1 | 0 |
City | Views | Downloads |
---|---|---|
Unknown | 13 | 0 |
Mountain View | 8 | 0 |
Nanjing | 4 | 2 |
Tokyo | 4 | 0 |
Ashburn | 3 | 0 |
Buffalo | 2 | 0 |
University Park | 1 | 0 |
Shenzhen | 1 | 5 |
Redmond | 1 | 0 |
Bangalore | 1 | 1 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080028/manifest