Abstract
Thermal ground testing is an accepted and frequently used method for simulating the aerodynamic heating of high speed flight vehicles. A numerical method based on a finite volume method for a quartz lamp heating system, used in thermal testing, is proposed. In this study, the unstructured finite-volume method (UFVM) for radiation has been formulated and implemented in a fluid flow solver GTEA on unstructured grids. For comparison and validation of the proposed method, a 2D furnace with cooling pipes was chosen. The results obtained from the proposed FVM agreed well with the exact solutions. Numerical results show that the quartz lamp can be simplified as a slat with the same temperature radiation source, and a simplified 2D thermal testing case was then simulated with the coupling effects of radiation, convection, and conduction heat transfer. Different temperature loading curves and ratios of intervals between the lamps and lamp length () were studied using the proposed method. The radiation heat flux on the metal surface was a wave-shaped curve. Comparing the different interval ratios, we found that the smaller the interval ratio, the larger the maximum value and the smaller the difference between the maximum and minimum heat flux.
1. Introduction
The quartz lamp heating system is widely used in industries and laboratories as a key part of rapid thermal processing (RTP) equipment or thermal testing. RTP is essential in the manufacturing of semiconductor devices such as integrated circuits, memories, or solar cells. They correspond to key stages in wafer production operations such as annealing (RTA), oxidation (RTO), or chemical vapor deposition (RTCVD) [1–3]. A quartz lamp heating system is also the most used heating method for thermal testing because of its high efficiency (approximately 90%). Thermal ground testing is an accepted and frequently used method to simulate the aerodynamic heating of high speed flight vehicles. Coupling radiation and conduction heat transfer are important in these problems. To model the combination of mechanical and thermal loads, the time history effects of aerodynamic heating cannot be overemphasized. Consequently, it is appropriate to study the structure fluid coupling transient heat transfer. Therefore, some attention has been focused on this field in the past decades [1]. In this study, a numerical algorithm based on the finite volume method is proposed for coupling the radiation and transient conjugate heat transfer of quartz lamp heating systems and structures.
The prediction of fluid flow has been realized by applying computational fluid dynamics (CFD). There are many CFD methods including stream function vorticity [4], original variable method [5], and finite difference (FD) and finite volume (FV) methods. The finite volume method is particularly popular because it leads to exact compliance with the conservation laws and direct solutions of the original variables. Owing to the direct extent of this method to convective and conduction heat transfer, the finite volume method is extensively used in solving fluid flow and convective heat transfer problems. In some cases, the effect of radiation heat transfer is important for natural convection problems. The mechanism of radiation heat transfer is quite different from convection and conduction approaches; therefore, there are many methods for radiation heat transfer, such as Hottel’s zone method, the Monte Carlo method, the discrete exchange factors method, the discrete ordinates method, and the P-N method, which are difficult to couple with the fluid flow solver. To bridge the gaps, Raithby and Chui proposed a finite volume method to solve the radiation transport equation in 1990 [6]. Kim et al. studied the radiation in an enclosure with obstacles using the unstructured finite volume method and the influence of types of manipulation of the control angle overlap [7]. Asllanaj et al. proposed a cell vertex based finite volume method for radiation heat transfer [8]. Kim et al. extended the unstructured finite volume method to polygonal grids [9]. Talukdar et al. have simulated radiation heat transfer on a multiblock grid using the finite volume method [10].
In this study, a modified FVM method based on unstructured grids is used to solve the radiation intensity and calculate the radiation heat transfer source term [11]. A complex geometry furnace with embedded cooling and quartz lamp heating system are used to evaluate the method. The quartz lamp heating system is then studied using the proposed numerical method. From the numerical results, the temperature history of the temperature control points is sometimes a very complex function of the loading curve. The goal of the study is to determine the optimized parameters and for a uniform flux and temperature distribution along the aluminum block. The control of the lamps can be manipulated to match a temperature required for the block.
2. Mathematical Model
2.1. Fluid Flow and Energy Control Equation
The governing equations for natural convection are set for nonlinear mass, momentum, and energy conservation equations:
The Boussinesq approximation is used to correlate the buoyant force and temperature variation:
2.2. Radiation Transfer Equation
For a gray medium, the radiation transfer equation (RTE) may be written, for a ray traveling along the direction, in the following form:
In this study, we ignore the scattering of the participating medium, and (3) is then rewritten in following form:
The boundary condition for a diffusely emitting and reflecting wall is written as
Both the emissivity and absorption coefficients are assumed to be constant and independent with wavelength.
2.3. Energy Equation in a Solid
In a solid zone, only conduction is considered and the control equation is written in the following form:
3. Numerical Methods
3.1. Finite Volume Method
The finite volume method for the N-S equation is well covered in the literature, so we only give the discretization process of the RTE. As mentioned above, we integrate the RTE over the control volume and solid angle, so (4) is written in the following form:
The Gauss theorem is used to transfer the volume integral to the surface integral, and the left-hand side of (7) can then be derived. Here, we suppose that the radiation intensity is uniform on the cell faces, so the following formulation holds. By assuming that the intensity is constant in the control volume and solid angle, the right-hand side of (7) is written in the following form:
Up to now an assumed distribution has been necessary to relate a control-volume facial intensity to a cell center intensity. Among many others, a step scheme is used in which a downstream facial intensity is set to be equal to the upstream cell value. It not only ensures positive intensity, but it is also simple and convenient. The details of the derivation can be found in the literature [11].
3.2. Hybrid Grid Method and Solid Angle Discretization
The computational domain is divided into a series of triangular, quadrilateral, or polygonal nonoverlap small cells as in Figure 1. All the solved variables and properties are stored at the cell center, such as in Figure 1. The control equation for each variable is integrated over a control volume and then written in discretized form (the details of this procedure are given below).
The solid angle, steradians, is discretized into () directions as given by Kim et al. [7], where is a polar angle ranging from 0 to and is an azimuthal angle ranging from 0 to in Figure 2. Here the control angle is equally divided such that and .
Once the angular limit for each control angle is fixed, the directional weight, which determines an inflow or outflow of radiant energy across the control-volume face depending on its sign, is evaluated as follows: where and .
3.3. Boundary Conditions
The boundary condition is implemented as follows for the gray medium:
3.4. Natural Convection and Radiation Interaction
Once the radiation intensity is obtained, we calculate the radiation source term with the following formulation:
For the adiabatic boundary condition, the wall temperature is the unknown of the problem. It is part of the solution obtained from the condition that the heat is balanced between radiation and conduction at the wall: where and are the radiation heat flux and unit outward normal vector of the wall.
3.5. Conjugate Heat Transfer
For the fluid cell and solid cell next to the interface as shown in Figure 3, the formulations hold according to Fourier’s law:
Because the heat flux through the interface is identified, we obtain the following equation:
From (13) to (15), we derive the temperature on the interface: where and .
Substituting (16) into the energy equations (13) and (14), the following formulations are reached:
3.6. Solution Algorithm
The RTE is solved using the FVM, and the heat source or sink and heat flux on the fluid solid interface is obtained by integrating the radiation intensity. Meanwhile, a SIMPLE based fluid solver is activated and energy conservation equation is then simultaneously solved both in the solid and fluid domains. As we assume that the time scale of radiation is much shorter than that of convection and heat conduction, steady radiation heat transfer is used in this study. The flow chart for every time step is given in Figure 4.
4. Numerical Results and Discussion
4.1. Radiation in a Furnace with Embedded Cooling Pipes
A geometric model of a furnace with cooling pipes which was used by Kim et al. [7] to validate their method was considered. The furnace was simplified as a two-dimensional enclosure 2.0 m wide and 6.29 m long. There were 22 cooling pipes along the top wall between 2.22 m and 6.29 m. The pipe, 0.03 m in diameter, was 0.025 m away from the top wall. The furnace wall and pipe emissivities were 0.5 and 0.8, respectively. The furnace wall and gas medium temperatures were obtained from the literature [12], while the cooling pipe temperature was assumed to be 300 K. The gas absorption coefficient was set at . The constant temperature distribution was adopted from the literature [12] as in Figure 5, which was given on four different lines, . The temperature was symmetric along the central line, while and corresponded to the central line and furnace wall, respectively. A ninth order polynomial fitting was used for the temperature along the given lines as in Figure 5 and a linear interpolation was then used in the direction from these data. The results are shown in Figure 6.
Figure 7 shows the incident radiation heat flux distribution along the top wall with and without cooling pipes. It was established that the cooling pipes marginally affected the heat flux. The incident heat flux on the top wall with no cooling pipes was a smooth curve, but with cooling pipes it fluctuated dramatically. Because the cooling pipes intercept the incident radiation coming from the hot medium toward the wall, shadows appear on the heat flux distribution. The incident heat flux was not zero directly behind the pipes because of the nonperpendicular radiation which was consistent with the literature [7]. With the present method, a detail about the radiation heat transfer is given and an optimization on the distance between the pipes and wall and pipe intervals can be done with this tool.
4.2. Quartz Lamp Heating System
Fields provided the basic quartz heater system characteristics and design and two specific applications, the space shuttle and flight, are also given [13]. The sample case sketch of thermal ground testing is shown in Figure 8 and generally includes the test piece, quartz lamp and loading equipment, and some measurement apparatus, which is not given in the figure. The quartz lamp structure shown in Figure 8 is a helical tungsten filament covered by quartz glass. The diameters of tungsten filament and quartz glass are 2 and 12 mm, respectively. The model of the quartz lamp was studied before we simulated the thermal ground testing case. We assumed that the quartz glass is transparent for radiation and the tungsten filament is blackbody, therefore the quartz lamp can be replaced by a circular with 2 mm in diameter. The temperature of the quartz lamp was 1200 K and the wall was assumed to be a blackbody, which absorbs all the radiation from the system and the radiation was neglected. The reflection plate was assumed to be an adiabatic boundary condition. There are many quartz lamps used in a quartz lamp heating system and the radiation heat transfer characteristics are studied before being used in a practical problem. A case of three lamps is used combined with a periodical assumption. The cyclic condition was used on the top and bottom boundaries. All boundaries ware assumed to be blackbodies and the absorption coefficient of air was .
The quartz lamp can be simplified as a slat [14]. Here numerical experiments with a slat and cylinder model are conducted to appreciate the accuracy of this assumption. Two models of the quartz lamp, with grids as shown in Figure 9, were simulated in the proposed method; the first was a circular model, 2 mm in diameter, and the other was a rectangular model. The computational domain was and there were 8 quartz lamps at intervals. The center of a lamp was away from the left reflection plate and from the test piece surface. The wall heat flux on the test piece surface is shown in Figure 10. We found that the results for cases 1 and 3 were quite approximate, as were the results for cases 2 and 4. The quartz lamp was modeled with a narrow slat radiation source whose length was equal to the lamp diameter. We also noted that the grid had a great influence on the results. Consequently, a rectangular radiation source was subsequently used instead of the quartz lamp.
(a)
(b)
(c)
(d)
The computational domain was a enclosure as in Figure 11(a) and the coordinates were constructed as in Figure 11(b). There were an aluminum block with on the left side and four quartz lamp heaters placed at intervals on the right side.
(a)
(b)
An adiabatic boundary condition was used on all boundaries except and , to which a Dirichlet boundary condition and coupling condition were applied, respectively. The initial temperature of the gas medium was and the heat temperature varied with time as in Figure 12. Constant material properties were used in this research as shown in Table 1.
The numerical results are shown in Figures 13–15. Figure 13 presents the temperature distribution in the solid domain at different times. There were two hot zones in the solid domain which was directly opposite the quartz lamp heater. The radiant heat flux was nonuniform because of the gap between the heaters, which caused a thin nonuniform temperature layer.
(a)
(b)
(c)
(a)
(b)
The temperature, velocity vector, and vorticity at time 70 s are shown in Figure 14, and five large eddies formed because of the buoyant force. The air near the quartz lamp surface was heated and the density reduced, so the air was driven upward. It was similar for the metal surface where the flow was also upward. This led to a high pressure near the top wall, which drove the fluid downward. Therefore a complex fluid flow occurred in the closure.
There was not only a complex fluid flow but also a complex heat transfer mechanism. The three heat transfer modes interacted with each other and were quite different at different times. Initially, the lamp temperature was low, the radiation was ignored, and the air velocity was very small, so the conduction was dominant in the first period. With increasing lamp temperature, the radiation increased dramatically and the large temperature differences created a buoyant flow. The convection and radiation combined heat transfer was dominant in this period. The temperature histories at different monitor points are shown in Figure 15(a). The temperature was low on MP E and K because there was no radiation incident on these two corner points, while the temperatures of the other five points were different. Because the MP G and I corresponded to the center of the quartz lamp and MP F H and I were at the center of the interval, the temperature was highest on MP G and I. The mechanism of the heat transfer during this process was complex, which was verified by the temperature histories of the monitor points shown in Figure 15(b). The monitor points G, C and I, B are symmetric about the x coordinate. The difference in the temperature history was caused by convection heat transfer.
Three different temperature loading curves as in Figure 12 were activated with the quartz lamp heater and the corresponding temperature histories are shown in Figure 16. Referring to Figure 16, the temperature history at the monitor points (MP), which are also called temperature control points, was sometimes a very complex function of the loading curve. When the quartz lamp temperature was less than 500 K, all temperature histories were the same. Case 3 was different to cases 1 and 2 after 130 s because the lamp temperature over the last 30 s was greater than 1200 K. Referring to Figure 17, the temperature history was approximately a cubic function of time when the quartz temperature was constant at 1200 K.
The effect of the interval width between the lamps was studied. Three different ratios of the width and lamp length were simulated using the proposed method. Case 1 as in Figure 12 used the lamp temperature loading curve. The radiation heat flux at different times is shown in Figure 18. The radiation heat flux on the metal surface was a wave-shaped curve. The maximum value was exactly opposite the center of the quartz lamp and the minimum value opposite the interval center. Comparing the different interval ratios, we found that the smaller the interval ratio, the larger the maximum value and the smaller the difference between the maximum and minimum heat flux.
| (a) |
| (b) |
| (c) |
The normalized radiation heat flux on the interface at time 120 s is shown in Figure 19. Increasing with the , the high heat flux points moved toward the center , and the heat flux at the central point increased. From the figure, we also found the ratio of maximum and minimum heat flux increased with increasing . An optimization ratio is to get uniform incident heat flux with uniform lamp temperature.
5. Conclusions
In this study a finite volume method (FVM) on hybrid grids for radiation and transient conjugate heat transfer is proposed and implemented using the fluid flow solver GTEA. To test the accuracy of the proposed method, a complex geometry furnace with cooling pipes is simulated. The numerical method was adopted to study the quartz lamp heating system. From the numerical results, the temperature histories at temperature control points were generally a very complex function of the loading curve. When the quartz lamp temperature was less than 500 K, all temperature histories were the same. With a quartz lamp temperature constant at 1200 K, case 3 was different to cases 1 and 2 after 130 s because the lamp temperature lasted for more than 30 s, but the trend was similar for the other two cases. The radiation heat flux on the metal surface was a wave shaped curve. The maximum value was exactly opposite the center of the quartz lamp and the minimum value opposite the interval center. Comparing the results for different interval ratios, we found that the smaller the interval ratio, the larger the maximum value and the smaller the difference between the maximum and minimum heat flux.
Nomenclature
| : | Body force in momentum equation |
| : | Directional weights (—) |
| : | Radiation flux on cell face |
| : | Grashof number (—) |
| : | Gravity acceleration |
| : | Enthalpy |
| : | Radiation intensity |
| : | Blackbody radiation intensity, |
| : | Unit vector in directions |
| : | Total number of polar (azimuthal) angles (—) |
| : | Pressure |
| : | Prandtl number (—) |
| : | Heat flux |
| : | Temperature |
| : | Reference temperature |
| Velocity vector | |
| : | Fluid. |
| : | Thermal diffusivity |
| : | Coefficient of thermal expansion |
| : | Azimuthal angle |
| : | Temperature ratio (—) |
| : | Absorption coefficient, |
| : | Dynamic viscosity |
| : | Kinetic viscosity |
| : | Polar angle |
| : | Emissivity (—) |
| : | Stefan-Boltzmann constant |
| : | Scattering coefficient |
| : | Discrete control angle |
| : | Density . |
| in out: | Radiation flux direction on the cell face |
| : | Index of radiation direction. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The financial support from the National Natural Science Foundation of China under Grant 51206031 is gratefully acknowledged. The authors thank the reviewers for their helpful advice to improve their paper.