- Email: [email protected]
Robust methodology for determination of heat transfer coefficient distribution in convection
Applied Thermal Engineering 30 (2010) 2815e2821
Contents lists available at ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
Robust methodology for determination of heat transfer coefficient distribution in convection Bowang Xiao a, *, Qigui Wang b, Gang Wang a, Richard D. Sisson Jr. c, Yiming Rong a a
Manufacturing Engineering, Worcester Polytechnic Institute, 100 Institute Road, Worcester, MA 01609, USA General Motor Company, Materials Technology, Global Powertrain Engineering, 823 Joslyn Avenue, Pontiac, MI 48340, USA c Materials Science and Engineering, Worcester Polytechnic Institute, 100 Institute Road, Worcester, MA 01609, USA b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 12 April 2010 Accepted 13 August 2010 Available online 20 August 2010
To predict the microstructures, residual stresses and distortions in the heat treated metal components, it is important to accurately know the heat transfer coefficients (HTCs) between the hot work piece and cooling media. In this paper, a new method is presented to accurately determine the node-based HTC distribution by coupling computational fluid dynamics (CFD) with optimal weight functions and scale factors. With this new method, the predicted temperature profile of the work piece during quenching (rapid cooling) is in excellent agreement with experimental measurements. This new method can be also applied to accurately predict convection heat transfer in thermal equipment such as heat exchangers and refrigerators, building thermal design and other heat transfer related situations. Ó 2010 Elsevier Ltd. All rights reserved.
Keywords: Heat transfer coefficient Convection Quenching Finite element analysis CFD
1. Introduction To improve mechanical properties, many heat treatable metal products are usually subject to heat treatment that involves a quenching process to produce supersaturated solid solution or martensitic phase transformation. A significant amount of residual stress can be developed in the quenched parts [1,2]. The existence of residual stresses can have a significant detrimental influence on the performance of a structural component. In many cases, the high tensile residual stresses can also result in a severe distortion of the component, and they can even cause cracking during quenching or subsequent manufacturing processes [3,4]. The magnitude of residual stress and distortion produced in the parts during quenching significantly depends on the cooling rate and the extent of non-uniformity of the temperature distribution in the parts. The heat transfer during quenching involves conduction as expressed in Equation (1), convection as expressed in Equation (2), and radiation as expressed in Equation (3). Experimental and numerical simulation results have shown that the convection HTC between the hot component and the quenching media plays an important role in resultant distortion, residual stress and hardness distribution of the quenched object [5e8].
* Corresponding author. Tel.: þ1 508 831 6160; fax: þ1 508 831 6412. E-mail address: [email protected] (B. Xiao). 1359-4311/$ e see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.applthermaleng.2010.08.017
vT Q_ cond ¼ KA vx
(1)
where Q_ cond ¼ heat transfer rate via conduction, W; K ¼ thermal conductivity of a material, W/m C; A ¼ cross section area, m2; vT/vx ¼ thermal gradient in the direction of the heat flow, C/m
Q_ conv ¼ hc AðT TN Þ
(2)
where Q_ conv ¼ heat transfer rate via convection, W; hc ¼ convection heat transfer coefficient, W/m2 C; A ¼ surface area, m2; T ¼ component surface temperature, C; TN ¼ quenchant tempera ture, C
4 Q_ rad ¼ 3sA T 4 TN
(3)
where Q_ rad ¼ heat transfer rate via radiation, W; s ¼ universal Stefan Boltzman constant, 5.6704 108 W/m2 K4; 3 ¼ emissivity of the body. To simulate the quenching process and predict microstructure, residual stresses and distortions in the as-quenched parts, accurate HTC distribution data are needed. Accurate HTC distribution data are also required for optimizing designing thermal equipment such as heat exchanges and refrigerators and thermal controls in buildings and electrical products [9e11]. There are many classical empirical equations reported in the literature for calculating convection HTC data by using dimensionless numbers (i.e. Reynolds number, Prandtl number and Nusselt
2816
B. Xiao et al. / Applied Thermal Engineering 30 (2010) 2815e2821
Fig. 1. Comparison of measured and calculated timeetemperature curves using a CFD code at two locations of an aluminum casting quenched in water [22].
number) [12]. However, their applications are limited since almost all of these dimensionless numbers are calibrated for specific experimental conditions which can be significantly different from the actual production process. An experimental approach with a small probe can be utilized to determine HTC. With the so-called lumped heat capacity method [13,14], the cooling curves of the probe are acquired and used to inversely calculate HTC in terms of probe surface temperature [12,14e18]. This method assumes that the temperature distribution in the small probe is uniform during cooling in order to calculate HTC inversely, which is also uniform over the probe surface [13]. Another method, known as the direct conduction method, assumes the temperature near a surface is a linear variation so that the HTC can be obtained based on the equilibrium of heat flux at the surface by conduction and convection [13,19]. It is, however, difficult to obtain HTC distribution data with this experimental method. In recent years, some methods have been reported to obtain HTC distributions for the entire heat transfer interface. One method called
iterative modification method, divides the part surfaces into several zones and optimizes the zone-based HTC distribution data iteratively by minimizing the errors between the predicted temperature distributions and the measured cooling curves at different locations of the part [14,20]. Although this method can provide accurate prediction of the temperature distributions of the locations where the thermocouples are instrumented, it cannot assure that the predicted temperature profiles at other locations in the part are matching the actual temperature distributions during quenching. Another limitation of the method is that temperature measurements are needed for the same part quenched at different conditions. With the recent advances in CFD, the code may be used to calculate the node-based HTC distribution [21]. However, no currently commercially available CFD code is designed to calculate accurate HTC distribution data. The current CFD prediction of heat transfer and temperature distribution of a quenched part is not accurate particularly in liquid quenching because the complicated interaction and heat transfer phenomena between liquid and the
START FILE: HTC distribution data exported from CFD simulation INITIAL Scale factors
MODIFY
Thermal simula tion in ABAQUS Read Temperature from ABAQUS results FILE: Experimental temperature vs. time curves
Calculate temperature differences ΔT
ΔT >= 5 oC and
YES
Iteration # <=10
NO
END
Fig. 2. Flow chart showing the procedure to optimize the HTC distribution (NOTE: the scale factors can vary with work piece temperature or quenching time) [23].
B. Xiao et al. / Applied Thermal Engineering 30 (2010) 2815e2821
2817
hot part are not fully understood and represented in the state-ofthe-art fluid flow and heat transfer code. Fig. 1 shows examples of the significant discrepancy observed in the thermal simulation using a state-of-the-art fluid flow and heat transfer code in comparison with experimental measurements in a water-quenched aluminum part [22]. In this paper, a new method to accurately determine node-based HTC distributions by coupling computational fluid dynamics (CFD) with optimal weight functions and scale factors is presented. 2. Integrated CFD and iterative modification method Fig. 2 summarizes the integrated CFD and iterative modification method for accurate determination of HTC distribution during quenching [23]. An initial set of node-based HTC distribution data are first obtained from the CFD simulation based on the part geometry, quenching setup, and quenching conditions including part initial temperature, material properties, quenchant flow velocity, quenching gas flow direction relative to the part, quenchant temperature and physical properties. The initial HTC values for the entire surface of the part calculated from the CFD simulation are then optimized by multiplying scale factors to minimize the errors between the predicted and the measured temperatures for a given quenching condition. Equations (4) and (5) illustrate how the scale factors are modified based on the temperature differences between the simulation and measurement. In the optimization process, the scale factors are modified iteratively until either acceptable temperature differences such as 5 C are reached or iteration number reaches the preset maximal iteration number (to prevent endless loop). It should be noted that the scale factors can vary with the surface temperatures of the part or quenching time. If the part surfaces are grouped into surface zones, different scale factors can be applied to various surface zones.
ScaleFactornew ¼ ScaleFactorold þ DScaleFactor
(4)
DScaleFactor ¼ k$DT
(5)
Fig. 4. A schematic illustration of the aluminum alloy casting with 14 thermocouples at selected locations.
modify the optimized baseline HTC distribution data for different quenching conditions (i.e., variations of quenching conditions such as flow velocity from the baseline) without performing complete heat transfer and optimization calculations, as shown in Fig. 3. For instance, a heat treating company has acquired some timeetemperature curves for a specific part under current production condition and has performed the optimization presented in Fig. 2 to obtain the HTC distribution data. Now the company increases the air velocity for a higher cooling rate. With the weight function for air velocity, the company can quickly determine the new HTC distribution by scaling up the old one. In this way, no quenching experiment is required for determining the new HTC distribution. Neither is the CFD simulation for the new condition.
where k is a constant and its value only affects the convergence rate. The calculation time in this optimization process varies with the complexity of the target part surface, the initial guess and the value of the constant k in Equation (5). For the optimization in the following validation section, the HTC data converge usually after about 10 iterations with a reasonable k. For a series of optimizations with similar conditions, the optimized HTC data for the previous condition can be used as an initial guess of the next optimization and the optimization process can even be completed after only 1 or 2 iterations. The calculation time for each iteration depends on the mesh used for the solid part, the time increment and the total cooling time. For the aluminum casting shown in Fig. 4, it takes less than 1 h to complete one iteration on a 4-CPU server by using Abaqus. When the node-based HTC data are optimized for a baseline quenching condition, a set of semi-empirical equations (or weight functions) as expressed in Equation (6) can be used to quickly
Scale factors CFD HTC simulation distribution
Experiment
New New HTC production Semi-empirical distribution condition equations
Fig. 3. A schematic illustration of the integrated CFD, iterative modification method and semi-empirical equations in determining node-based HTC distribution.
Fig. 5. The experimental system for the forced air quenching of the aluminum alloy casting.
2818
B. Xiao et al. / Applied Thermal Engineering 30 (2010) 2815e2821
modification factor and the air velocity for the velocity ranging from 4 m/s to 20 m/s, which was determined by experiments [24].
Kvelocity ¼ A
Vel Vel0
þB
(7)
where A ¼ 0.57; B ¼ 0.41; and Velo ¼ 10.5 m/s. 3. Validation and application of the method
Fig. 6. A sample timeetemperature profile at the vertical orientation and 18 m/s gas velocity.
HTC ¼ K1 K2 K3 /Kn $HTCo
(6)
where HTCo is the standard HTC at a baseline quenching condition, unit in W/m2 C; K1, K2, K3, ., Kn are modification factors. These modification factors are governing the effects of various influencing factors such as fluid velocity, fluid flow direction, fluid temperature, part surface quality and material. These modification factors can be either constant or temperature-dependent. They can be determined experimentally or numerically by CFD simulations. For instance, Equation (7) shows the relationship between the velocity
To validate the method, a mixed leg thickness picture-frame shape aluminum casting was utilized in the forced air quenching experiments. Fig. 4 shows the casting instrumented with 14 thermocouples. The dots and lengths of dash lines illustrate the locations of the thermocouples in three dimensions. In air quenching experiments, the casting was heated in a furnace to about 485 C and then placed on the fixture and quenched in the forced air, as shown in Fig. 5. The cooling curves at 14 locations of the casting during quenching were acquired by a data acquisition system. Temperature measurements were taken 5 times per second. Fig. 6 shows one exemplary timeetemperature profile for the aluminum casting during quenching at the vertical orientation as shown in Figs. 7 and 8. The air temperature, pressure and humidity were measured by a weather station (La Crosse Technology WS-7394U-IT-CH Wireless Forecast Station). The air velocity at the quenching area was calibrated as 18 m/s with an anemometer. A CFD simulation was conducted based on the experimental setup, as shown in Fig. 7. A big air box around the casting was created. Boundary conditions were applied to the air box. One face was defined as the air velocity inlet where normal air flow at specified air velocity was blown in. All other faces were defined as
Fig. 7. CFD modeling.
B. Xiao et al. / Applied Thermal Engineering 30 (2010) 2815e2821
2819
Fig. 10. A comparison of temperatureetime curves between the measurements and FEA simulation using the optimized HTC data (one simulation using original HTC for comparison purpose). Fig. 8. A color contour showing the CFD-predicted HTC distribution in the aluminum test casting during forced air quenching [23].
opening boundary where air can go in and out depending on the pressure difference. Keepsilon turbulence model and P1 radiation model were applied in this CFD simulation. Simulation results showed that radiation slightly affects the temperature cooling curve. In this simulation, air is treated as ideal gas. In the CFD simulation, stagnant air (initial status) was forced to flow and a steady air flow field was formed after a few seconds. In the simulation, fresh air was blown to the cooling zone and thus the inlet air temperature remained constant during simulation. Since the air flow (18 m/s, 1 atmospheric pressure) was set steady and air temperature remained constant, the calculated HTC distribution is considered steady. Fig. 8 presents the predicted node-based HTC distribution from the CFD simulation when the air flow has become steady. As expected, the predicted node-based HTC distribution
varies from location to location even on the same surface of the casting. However, the predicted node-based HTC distribution directly from the CFD simulation is not accurate. As shown in Fig. 9, the simulated cooling curves from the finite element analysis using the original CFD-determined HTC data significantly differ from the experimentally measurements at the locations as shown in Fig. 4. Accordingly, the original CFD-determined HTC data are modified and optimized with the scale factor as described in Equations (4) and (5). After the HTCs are optimized with scale factors, the nodebased HTC distribution represents better heat transfer phenomena during quenching. As shown in Fig. 10, the predicted cooling curves using the optimized HTC data agree very well with the experimentally measured temperature distributions. The maximum error (or difference) between the predicted temperatures and the measurements is þ/5 C. Fig. 10 also shows how the optimized HTC data drag the timeetemperature curve at location TC2 towards the experimental one and reduce the temperature differences from above 50 C to less than 5 C. As presented in Fig. 11, the corresponding scale factor varies with casting surface temperature and the original HTC magnitude is scaled up by over 30%. The optimization is further validated by comparing the cooling rates. Fig. 12 compares the experimentally determined cooling rates at TC2, TC3, TC5 and TC9 (on four different walls), predicted cooling rate using CFD-determined HTC data at TC2, and predicted cooling rate using optimized HTC data at TC2. The
Scale Factor Scale Factor
1.4 1.3 1.2 1.1 1 0
Fig. 9. A comparison of temperatureetime curves between the measurements and FEA simulation using the original CFD-determined HTC data (The differences between experimental TC3 and TC5 are very small).
100
200
300
400
Temperature (oC ) Fig. 11. The scale factor as a function of part surface temperature.
500
2820
B. Xiao et al. / Applied Thermal Engineering 30 (2010) 2815e2821
Cooling Rate Comparison 0 -1
Cooling rate ( oC/s)
-2 2 measured
-3
2 simulated (after optimization)
-4
2 simulated (before optimization)
-5
3 measured
-6 5 measured
-7
9 measured
-8 0
100
200 300 Surface Temperature ( oC )
400
500
Fig. 12. A comparison of cooling rates among experimentally obtained and numerically predicted ones.
experimental cooling rates are generated from the experimental cooling curves, which are moving averaged to remove abrupt oscillations. It is apparent that the predicted cooling rate using optimized HTC data agrees very well with the experimental one, but the predicted cooling rate using original CFD-determined HTC data differs from the experimental one significantly. In addition, the cooling rate plot validates two other things. First, at the beginning, because of the transportation of casting from the furnace to the cooling zone, there are big changes of cooling rates at the beginning of cooling. Second, the cooling rates in the main range are very linear, indicating that HTC data during gas cooling remain quite constant after the gas flow becomes steady. This conclusion can be drawn from the following derivation. In the quenching process, the heat loss rate of the hot part can be measured by Equation (8).
Q_ loss ¼ Cp m$dT=dt
(8)
where Q_ loss ¼ heat loss rate during cooling, W; Cp ¼ specific heat of the casting, J/kg C; m ¼ mass, kg; dT/dt ¼ cooling rate, C/s. If we simplify the problem and assume the heat loss of a small portion of the part is transferred by convection, since radiation is very small in this case, we can approximate Equation (8) as Equation (2). Thus, after substituting Equation (8) into Equation (2), we have the cooling rate as a function of surface temperature as expressed in Equation (9). The near linear relationship between cooling rate and surface temperature implies that HTC is close to a constant.
dT=dt ¼
hc A hc A T TN Cp m Cp m
(9)
complicated. As shown in Fig. 9, the temperature differences between the measurements and the predictions based on the HTC obtained directly from CFD simulation can be as high as 50 C for the relatively simple shape aluminum test casting. Therefore, the thermal calculation from CFD simulation cannot be directly used in structural analysis particularly for the residual stress and distortion prediction. To accurately simulate thermal history and predict residual stress and distortion in the quenched part, the node-based HTC distribution data determined from CFD simulation need to be optimized. The CFD-determined HTC distribution can be utilized as an initial set of HTC data for further optimization. In the HTC optimization, the node-based HTC distribution data are modified iteratively with scale factors or weight functions (semi-empirical equations). The scale factors and weight functions are not only temperature-dependent but also part surface geometry dependent. With the measured cooling curves from the quenched part, the scale factors and/or weight functions can be calibrated and optimized for a given quenching condition. As shown in Fig. 10, the temperature difference between the measurements and the predictions based on the optimal HTC distribution can be less than 5 C. Although this new method was initially put forward for quenching processes, it can be also applied to other heat treating processes like heating, annealing, tempering, etc., thermal design and management and others, whenever there is convection heat transfer and accurate prediction of temperature distribution in a part is required. 5. Conclusion In this paper, a new method integrating CFD simulation and an iterative modification and optimization method was developed to determine accurate node-based HTC distribution for a part with complicated geometry during rapid cooling. With this method, the accurate HTC distribution data can be obtained and thus transient heat transfer in the cooled part can be accurately predicted. When further integrated with semi-empirical equations (weight functions), this method can be used to determine new HTC distributions for the same part under different quenching conditions at no cost of experiments and simulations. This method can be applied in other numerical predictions of convection heat transfer, such as optimization of heat exchanger design, thermal management of electrical products, etc. Acknowledgements This work was supported financially by GM Powertrain in Pontiac, Michigan. The authors thank GM R&D Foundry for producing the aluminum alloy castings. The gratitude is also sent to CHTE (Center for Heat Treating Excellence) members for their sponsored quenching system. References
4. Discussion Although much commercial CFD software can solve thermal and fluid problems to a certain degree of accuracy, the software cannot accurately calculate the heat transfer coefficients and temperature distributions during highly transient heat transfer processes like quenching. This is not only because many simplifications and assumptions are assumed in the CFD calculations but also because the physics and the phase transformation phenomena involved at the interface between the hot part and cooling media are
[1] I. Elkatatny, Y. Morsi, A.S. Blicblau, S. Das, E.D. Doyle, Numerical analysis and experimental validation of high pressure gas quenching, Int. J. Therm. Sci. 42 (4) (2003) 417e423. [2] K. Li, B. Xiao, Q. Wang, Residual stresses in as-quenched aluminum castings, SAE Int. J. Mater. Manuf. 1 (1) (2009) 725e731. [3] P. Li, D.M. Maijer, T.C. Lindley, P.D. Lee, Simulating the residual stress in an A356 automotive wheel and its impact on fatigue life, Metall. Mater. Trans. B. 38 (4) (2007) 505e515. [4] Y.L. Lee, J. Pan, R. Hathaway, M. Barkey, Fatigue Testing and Analysis: Theory and Practice. Elsevier Butterworth-Heinemann, 2005. [5] A. Rose, O. Kessler, F. Hoffmann, H.W. Zoch, Quenching distortion of aluminum castings-improvement by gas cooling, Materialwiss. Werkst. 37 (1) (2006) 116e121.
B. Xiao et al. / Applied Thermal Engineering 30 (2010) 2815e2821 [6] Z. Li, R.V. Grandhi, R. Srinivasan, Distortion minimization during gas quenching process, J. Mater. Process. Technol. 172 (2) (2006) 249e257. [7] H. Li, G. Zhao, S. Niu, Y. Luan, Technologic parameter optimization of gas quenching process using response surface method, Comput. Mater. Sci. 38 (4) (2007) 561e570. [8] B. Xiao, G. Wang, Y. Rong, Hardenability and distortion control in high pressure hydrogen quenching, Int. J. Manuf. Res., in press. [9] B. Sahin, Y. Ust, I. Teke, H.H. Erdem, Performance analysis and optimization of heat exchangers: a new thermoeconomic approach, Appl. Therm. Eng. 30 (2e3) (2010) 104e109. [10] G.N. Xie, B. Sunden, Q.W. Wang, Optimization of compact heat exchangers by a genetic algorithm, Appl. Therm. Eng. 28 (8e9) (2008) 895e906. [11] Z. Luo, H. Cho, X. Luo, K. Cho, System thermal analysis for mobile phone, Appl. Therm. Eng. 28 (14e15) (2008) 1889e1895. [12] J.P. Holman, Heat Transfer, ninth ed. McGraw-Hill, New York, 2002. [13] M. Sedighi, C.A. McMahon, The influence of quenchant agitation on the heat transfer coefficient and residual stress development in the quenching of steels, Proc. Inst. Mech. Eng. Pt. B: J. Eng. Manuf. 214 (7) (2000) 555e567. [14] A. Sugianto, M. Narazaki, M. Kogawara, A. Shirayori, A comparative study on determination method of heat transfer coefficient using inverse heat transfer and iterative modification, J. Mater. Process. Technol. 209 (10) (2008) 4627e4632. [15] M. Maniruzzaman, J. Chaves, C. McGee, S. Ma, R. Sisson Jr., CHTE quench probe system: a new quenchant characterization system, in: 5th International Conference on Frontiers of Design and Manufacturing (ICFDM), 2002, pp. 619e625.
2821
[16] D.K. Funatani, M. Narazaki, M. Tanaka, Comparisons of probe design and cooling curve analysis methods, in: 19th ASM Heat Treating Society Conference Proceedings Including Steel Heat Treating in the New Millennium, 1999, pp. 255e263. [17] G.E. Totten, C.E. Bates, N.A. Clinton, Handbook of Quenchants and Quenching Technology. ASM International, Materials Park, OH, 1993. [18] N.I. Kobasko, A.A. Moskalenko, G.E. Totten, G.M. Webster, Experimental determination of the first and second critical heat flux densities and quench process characterization, J. Mater. Eng. Perform. 6 (1) (1997) 93e101. [19] C. Laumen, S.I. Midea, T. Lubben, F. Hoffman, P. Mayr, Measured heat transfer coefficients by using hydrogen as quenchant in comparison with helium and nitrogen, Accelerated Cooling/Direct Quenching of Steels, 1997, pp. 199e206. [20] M. Li, J.E. Allison, Determination of thermal boundary conditions for the casting and quenching process with the optimization tool OptCast, Metall. Mater. Trans. B 38B (4) (2007) 567e574. [21] G. Wang, B. Xiao, L. Zhang, K. Rong, Q. Wang, Heat transfer coefficient study for air quenching by experiment and CFD modeling, in: 137th Annual Meeting & Exhibition. Materials Characterization, Computation and Modeling, vol. 2, 2008, pp. 271e276. [22] O. Dahlberg, Z. Wu, Feasibility study of three dimensional quenching analysis, CD-adapco Technical Report to GM. 071-186-003, 2008. [23] Q. Wang, B. Xiao, G. Wang, Y. Rong, R.D.J. Sisson, Systems and methods for predicting heat transfer coefficients during quenching (Pending Patent), 2010. [24] B. Xiao, G. Wang, Q. Wang, M. Maniruzzaman, R.D. Sisson Jr., Y. Rong, An experimental study of heat transfer during forced air convection, J. Mater. Eng. Perform., in press, doi: 10.1007/s11665-010-9745-7.
Contents lists available at ScienceDirect
Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng
Robust methodology for determination of heat transfer coefficient distribution in convection Bowang Xiao a, *, Qigui Wang b, Gang Wang a, Richard D. Sisson Jr. c, Yiming Rong a a
Manufacturing Engineering, Worcester Polytechnic Institute, 100 Institute Road, Worcester, MA 01609, USA General Motor Company, Materials Technology, Global Powertrain Engineering, 823 Joslyn Avenue, Pontiac, MI 48340, USA c Materials Science and Engineering, Worcester Polytechnic Institute, 100 Institute Road, Worcester, MA 01609, USA b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 12 April 2010 Accepted 13 August 2010 Available online 20 August 2010
To predict the microstructures, residual stresses and distortions in the heat treated metal components, it is important to accurately know the heat transfer coefficients (HTCs) between the hot work piece and cooling media. In this paper, a new method is presented to accurately determine the node-based HTC distribution by coupling computational fluid dynamics (CFD) with optimal weight functions and scale factors. With this new method, the predicted temperature profile of the work piece during quenching (rapid cooling) is in excellent agreement with experimental measurements. This new method can be also applied to accurately predict convection heat transfer in thermal equipment such as heat exchangers and refrigerators, building thermal design and other heat transfer related situations. Ó 2010 Elsevier Ltd. All rights reserved.
Keywords: Heat transfer coefficient Convection Quenching Finite element analysis CFD
1. Introduction To improve mechanical properties, many heat treatable metal products are usually subject to heat treatment that involves a quenching process to produce supersaturated solid solution or martensitic phase transformation. A significant amount of residual stress can be developed in the quenched parts [1,2]. The existence of residual stresses can have a significant detrimental influence on the performance of a structural component. In many cases, the high tensile residual stresses can also result in a severe distortion of the component, and they can even cause cracking during quenching or subsequent manufacturing processes [3,4]. The magnitude of residual stress and distortion produced in the parts during quenching significantly depends on the cooling rate and the extent of non-uniformity of the temperature distribution in the parts. The heat transfer during quenching involves conduction as expressed in Equation (1), convection as expressed in Equation (2), and radiation as expressed in Equation (3). Experimental and numerical simulation results have shown that the convection HTC between the hot component and the quenching media plays an important role in resultant distortion, residual stress and hardness distribution of the quenched object [5e8].
* Corresponding author. Tel.: þ1 508 831 6160; fax: þ1 508 831 6412. E-mail address: [email protected] (B. Xiao). 1359-4311/$ e see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.applthermaleng.2010.08.017
vT Q_ cond ¼ KA vx
(1)
where Q_ cond ¼ heat transfer rate via conduction, W; K ¼ thermal conductivity of a material, W/m C; A ¼ cross section area, m2; vT/vx ¼ thermal gradient in the direction of the heat flow, C/m
Q_ conv ¼ hc AðT TN Þ
(2)
where Q_ conv ¼ heat transfer rate via convection, W; hc ¼ convection heat transfer coefficient, W/m2 C; A ¼ surface area, m2; T ¼ component surface temperature, C; TN ¼ quenchant tempera ture, C
4 Q_ rad ¼ 3sA T 4 TN
(3)
where Q_ rad ¼ heat transfer rate via radiation, W; s ¼ universal Stefan Boltzman constant, 5.6704 108 W/m2 K4; 3 ¼ emissivity of the body. To simulate the quenching process and predict microstructure, residual stresses and distortions in the as-quenched parts, accurate HTC distribution data are needed. Accurate HTC distribution data are also required for optimizing designing thermal equipment such as heat exchanges and refrigerators and thermal controls in buildings and electrical products [9e11]. There are many classical empirical equations reported in the literature for calculating convection HTC data by using dimensionless numbers (i.e. Reynolds number, Prandtl number and Nusselt
2816
B. Xiao et al. / Applied Thermal Engineering 30 (2010) 2815e2821
Fig. 1. Comparison of measured and calculated timeetemperature curves using a CFD code at two locations of an aluminum casting quenched in water [22].
number) [12]. However, their applications are limited since almost all of these dimensionless numbers are calibrated for specific experimental conditions which can be significantly different from the actual production process. An experimental approach with a small probe can be utilized to determine HTC. With the so-called lumped heat capacity method [13,14], the cooling curves of the probe are acquired and used to inversely calculate HTC in terms of probe surface temperature [12,14e18]. This method assumes that the temperature distribution in the small probe is uniform during cooling in order to calculate HTC inversely, which is also uniform over the probe surface [13]. Another method, known as the direct conduction method, assumes the temperature near a surface is a linear variation so that the HTC can be obtained based on the equilibrium of heat flux at the surface by conduction and convection [13,19]. It is, however, difficult to obtain HTC distribution data with this experimental method. In recent years, some methods have been reported to obtain HTC distributions for the entire heat transfer interface. One method called
iterative modification method, divides the part surfaces into several zones and optimizes the zone-based HTC distribution data iteratively by minimizing the errors between the predicted temperature distributions and the measured cooling curves at different locations of the part [14,20]. Although this method can provide accurate prediction of the temperature distributions of the locations where the thermocouples are instrumented, it cannot assure that the predicted temperature profiles at other locations in the part are matching the actual temperature distributions during quenching. Another limitation of the method is that temperature measurements are needed for the same part quenched at different conditions. With the recent advances in CFD, the code may be used to calculate the node-based HTC distribution [21]. However, no currently commercially available CFD code is designed to calculate accurate HTC distribution data. The current CFD prediction of heat transfer and temperature distribution of a quenched part is not accurate particularly in liquid quenching because the complicated interaction and heat transfer phenomena between liquid and the
START FILE: HTC distribution data exported from CFD simulation INITIAL Scale factors
MODIFY
Thermal simula tion in ABAQUS Read Temperature from ABAQUS results FILE: Experimental temperature vs. time curves
Calculate temperature differences ΔT
ΔT >= 5 oC and
YES
Iteration # <=10
NO
END
Fig. 2. Flow chart showing the procedure to optimize the HTC distribution (NOTE: the scale factors can vary with work piece temperature or quenching time) [23].
B. Xiao et al. / Applied Thermal Engineering 30 (2010) 2815e2821
2817
hot part are not fully understood and represented in the state-ofthe-art fluid flow and heat transfer code. Fig. 1 shows examples of the significant discrepancy observed in the thermal simulation using a state-of-the-art fluid flow and heat transfer code in comparison with experimental measurements in a water-quenched aluminum part [22]. In this paper, a new method to accurately determine node-based HTC distributions by coupling computational fluid dynamics (CFD) with optimal weight functions and scale factors is presented. 2. Integrated CFD and iterative modification method Fig. 2 summarizes the integrated CFD and iterative modification method for accurate determination of HTC distribution during quenching [23]. An initial set of node-based HTC distribution data are first obtained from the CFD simulation based on the part geometry, quenching setup, and quenching conditions including part initial temperature, material properties, quenchant flow velocity, quenching gas flow direction relative to the part, quenchant temperature and physical properties. The initial HTC values for the entire surface of the part calculated from the CFD simulation are then optimized by multiplying scale factors to minimize the errors between the predicted and the measured temperatures for a given quenching condition. Equations (4) and (5) illustrate how the scale factors are modified based on the temperature differences between the simulation and measurement. In the optimization process, the scale factors are modified iteratively until either acceptable temperature differences such as 5 C are reached or iteration number reaches the preset maximal iteration number (to prevent endless loop). It should be noted that the scale factors can vary with the surface temperatures of the part or quenching time. If the part surfaces are grouped into surface zones, different scale factors can be applied to various surface zones.
ScaleFactornew ¼ ScaleFactorold þ DScaleFactor
(4)
DScaleFactor ¼ k$DT
(5)
Fig. 4. A schematic illustration of the aluminum alloy casting with 14 thermocouples at selected locations.
modify the optimized baseline HTC distribution data for different quenching conditions (i.e., variations of quenching conditions such as flow velocity from the baseline) without performing complete heat transfer and optimization calculations, as shown in Fig. 3. For instance, a heat treating company has acquired some timeetemperature curves for a specific part under current production condition and has performed the optimization presented in Fig. 2 to obtain the HTC distribution data. Now the company increases the air velocity for a higher cooling rate. With the weight function for air velocity, the company can quickly determine the new HTC distribution by scaling up the old one. In this way, no quenching experiment is required for determining the new HTC distribution. Neither is the CFD simulation for the new condition.
where k is a constant and its value only affects the convergence rate. The calculation time in this optimization process varies with the complexity of the target part surface, the initial guess and the value of the constant k in Equation (5). For the optimization in the following validation section, the HTC data converge usually after about 10 iterations with a reasonable k. For a series of optimizations with similar conditions, the optimized HTC data for the previous condition can be used as an initial guess of the next optimization and the optimization process can even be completed after only 1 or 2 iterations. The calculation time for each iteration depends on the mesh used for the solid part, the time increment and the total cooling time. For the aluminum casting shown in Fig. 4, it takes less than 1 h to complete one iteration on a 4-CPU server by using Abaqus. When the node-based HTC data are optimized for a baseline quenching condition, a set of semi-empirical equations (or weight functions) as expressed in Equation (6) can be used to quickly
Scale factors CFD HTC simulation distribution
Experiment
New New HTC production Semi-empirical distribution condition equations
Fig. 3. A schematic illustration of the integrated CFD, iterative modification method and semi-empirical equations in determining node-based HTC distribution.
Fig. 5. The experimental system for the forced air quenching of the aluminum alloy casting.
2818
B. Xiao et al. / Applied Thermal Engineering 30 (2010) 2815e2821
modification factor and the air velocity for the velocity ranging from 4 m/s to 20 m/s, which was determined by experiments [24].
Kvelocity ¼ A
Vel Vel0
þB
(7)
where A ¼ 0.57; B ¼ 0.41; and Velo ¼ 10.5 m/s. 3. Validation and application of the method
Fig. 6. A sample timeetemperature profile at the vertical orientation and 18 m/s gas velocity.
HTC ¼ K1 K2 K3 /Kn $HTCo
(6)
where HTCo is the standard HTC at a baseline quenching condition, unit in W/m2 C; K1, K2, K3, ., Kn are modification factors. These modification factors are governing the effects of various influencing factors such as fluid velocity, fluid flow direction, fluid temperature, part surface quality and material. These modification factors can be either constant or temperature-dependent. They can be determined experimentally or numerically by CFD simulations. For instance, Equation (7) shows the relationship between the velocity
To validate the method, a mixed leg thickness picture-frame shape aluminum casting was utilized in the forced air quenching experiments. Fig. 4 shows the casting instrumented with 14 thermocouples. The dots and lengths of dash lines illustrate the locations of the thermocouples in three dimensions. In air quenching experiments, the casting was heated in a furnace to about 485 C and then placed on the fixture and quenched in the forced air, as shown in Fig. 5. The cooling curves at 14 locations of the casting during quenching were acquired by a data acquisition system. Temperature measurements were taken 5 times per second. Fig. 6 shows one exemplary timeetemperature profile for the aluminum casting during quenching at the vertical orientation as shown in Figs. 7 and 8. The air temperature, pressure and humidity were measured by a weather station (La Crosse Technology WS-7394U-IT-CH Wireless Forecast Station). The air velocity at the quenching area was calibrated as 18 m/s with an anemometer. A CFD simulation was conducted based on the experimental setup, as shown in Fig. 7. A big air box around the casting was created. Boundary conditions were applied to the air box. One face was defined as the air velocity inlet where normal air flow at specified air velocity was blown in. All other faces were defined as
Fig. 7. CFD modeling.
B. Xiao et al. / Applied Thermal Engineering 30 (2010) 2815e2821
2819
Fig. 10. A comparison of temperatureetime curves between the measurements and FEA simulation using the optimized HTC data (one simulation using original HTC for comparison purpose). Fig. 8. A color contour showing the CFD-predicted HTC distribution in the aluminum test casting during forced air quenching [23].
opening boundary where air can go in and out depending on the pressure difference. Keepsilon turbulence model and P1 radiation model were applied in this CFD simulation. Simulation results showed that radiation slightly affects the temperature cooling curve. In this simulation, air is treated as ideal gas. In the CFD simulation, stagnant air (initial status) was forced to flow and a steady air flow field was formed after a few seconds. In the simulation, fresh air was blown to the cooling zone and thus the inlet air temperature remained constant during simulation. Since the air flow (18 m/s, 1 atmospheric pressure) was set steady and air temperature remained constant, the calculated HTC distribution is considered steady. Fig. 8 presents the predicted node-based HTC distribution from the CFD simulation when the air flow has become steady. As expected, the predicted node-based HTC distribution
varies from location to location even on the same surface of the casting. However, the predicted node-based HTC distribution directly from the CFD simulation is not accurate. As shown in Fig. 9, the simulated cooling curves from the finite element analysis using the original CFD-determined HTC data significantly differ from the experimentally measurements at the locations as shown in Fig. 4. Accordingly, the original CFD-determined HTC data are modified and optimized with the scale factor as described in Equations (4) and (5). After the HTCs are optimized with scale factors, the nodebased HTC distribution represents better heat transfer phenomena during quenching. As shown in Fig. 10, the predicted cooling curves using the optimized HTC data agree very well with the experimentally measured temperature distributions. The maximum error (or difference) between the predicted temperatures and the measurements is þ/5 C. Fig. 10 also shows how the optimized HTC data drag the timeetemperature curve at location TC2 towards the experimental one and reduce the temperature differences from above 50 C to less than 5 C. As presented in Fig. 11, the corresponding scale factor varies with casting surface temperature and the original HTC magnitude is scaled up by over 30%. The optimization is further validated by comparing the cooling rates. Fig. 12 compares the experimentally determined cooling rates at TC2, TC3, TC5 and TC9 (on four different walls), predicted cooling rate using CFD-determined HTC data at TC2, and predicted cooling rate using optimized HTC data at TC2. The
Scale Factor Scale Factor
1.4 1.3 1.2 1.1 1 0
Fig. 9. A comparison of temperatureetime curves between the measurements and FEA simulation using the original CFD-determined HTC data (The differences between experimental TC3 and TC5 are very small).
100
200
300
400
Temperature (oC ) Fig. 11. The scale factor as a function of part surface temperature.
500
2820
B. Xiao et al. / Applied Thermal Engineering 30 (2010) 2815e2821
Cooling Rate Comparison 0 -1
Cooling rate ( oC/s)
-2 2 measured
-3
2 simulated (after optimization)
-4
2 simulated (before optimization)
-5
3 measured
-6 5 measured
-7
9 measured
-8 0
100
200 300 Surface Temperature ( oC )
400
500
Fig. 12. A comparison of cooling rates among experimentally obtained and numerically predicted ones.
experimental cooling rates are generated from the experimental cooling curves, which are moving averaged to remove abrupt oscillations. It is apparent that the predicted cooling rate using optimized HTC data agrees very well with the experimental one, but the predicted cooling rate using original CFD-determined HTC data differs from the experimental one significantly. In addition, the cooling rate plot validates two other things. First, at the beginning, because of the transportation of casting from the furnace to the cooling zone, there are big changes of cooling rates at the beginning of cooling. Second, the cooling rates in the main range are very linear, indicating that HTC data during gas cooling remain quite constant after the gas flow becomes steady. This conclusion can be drawn from the following derivation. In the quenching process, the heat loss rate of the hot part can be measured by Equation (8).
Q_ loss ¼ Cp m$dT=dt
(8)
where Q_ loss ¼ heat loss rate during cooling, W; Cp ¼ specific heat of the casting, J/kg C; m ¼ mass, kg; dT/dt ¼ cooling rate, C/s. If we simplify the problem and assume the heat loss of a small portion of the part is transferred by convection, since radiation is very small in this case, we can approximate Equation (8) as Equation (2). Thus, after substituting Equation (8) into Equation (2), we have the cooling rate as a function of surface temperature as expressed in Equation (9). The near linear relationship between cooling rate and surface temperature implies that HTC is close to a constant.
dT=dt ¼
hc A hc A T TN Cp m Cp m
(9)
complicated. As shown in Fig. 9, the temperature differences between the measurements and the predictions based on the HTC obtained directly from CFD simulation can be as high as 50 C for the relatively simple shape aluminum test casting. Therefore, the thermal calculation from CFD simulation cannot be directly used in structural analysis particularly for the residual stress and distortion prediction. To accurately simulate thermal history and predict residual stress and distortion in the quenched part, the node-based HTC distribution data determined from CFD simulation need to be optimized. The CFD-determined HTC distribution can be utilized as an initial set of HTC data for further optimization. In the HTC optimization, the node-based HTC distribution data are modified iteratively with scale factors or weight functions (semi-empirical equations). The scale factors and weight functions are not only temperature-dependent but also part surface geometry dependent. With the measured cooling curves from the quenched part, the scale factors and/or weight functions can be calibrated and optimized for a given quenching condition. As shown in Fig. 10, the temperature difference between the measurements and the predictions based on the optimal HTC distribution can be less than 5 C. Although this new method was initially put forward for quenching processes, it can be also applied to other heat treating processes like heating, annealing, tempering, etc., thermal design and management and others, whenever there is convection heat transfer and accurate prediction of temperature distribution in a part is required. 5. Conclusion In this paper, a new method integrating CFD simulation and an iterative modification and optimization method was developed to determine accurate node-based HTC distribution for a part with complicated geometry during rapid cooling. With this method, the accurate HTC distribution data can be obtained and thus transient heat transfer in the cooled part can be accurately predicted. When further integrated with semi-empirical equations (weight functions), this method can be used to determine new HTC distributions for the same part under different quenching conditions at no cost of experiments and simulations. This method can be applied in other numerical predictions of convection heat transfer, such as optimization of heat exchanger design, thermal management of electrical products, etc. Acknowledgements This work was supported financially by GM Powertrain in Pontiac, Michigan. The authors thank GM R&D Foundry for producing the aluminum alloy castings. The gratitude is also sent to CHTE (Center for Heat Treating Excellence) members for their sponsored quenching system. References
4. Discussion Although much commercial CFD software can solve thermal and fluid problems to a certain degree of accuracy, the software cannot accurately calculate the heat transfer coefficients and temperature distributions during highly transient heat transfer processes like quenching. This is not only because many simplifications and assumptions are assumed in the CFD calculations but also because the physics and the phase transformation phenomena involved at the interface between the hot part and cooling media are
[1] I. Elkatatny, Y. Morsi, A.S. Blicblau, S. Das, E.D. Doyle, Numerical analysis and experimental validation of high pressure gas quenching, Int. J. Therm. Sci. 42 (4) (2003) 417e423. [2] K. Li, B. Xiao, Q. Wang, Residual stresses in as-quenched aluminum castings, SAE Int. J. Mater. Manuf. 1 (1) (2009) 725e731. [3] P. Li, D.M. Maijer, T.C. Lindley, P.D. Lee, Simulating the residual stress in an A356 automotive wheel and its impact on fatigue life, Metall. Mater. Trans. B. 38 (4) (2007) 505e515. [4] Y.L. Lee, J. Pan, R. Hathaway, M. Barkey, Fatigue Testing and Analysis: Theory and Practice. Elsevier Butterworth-Heinemann, 2005. [5] A. Rose, O. Kessler, F. Hoffmann, H.W. Zoch, Quenching distortion of aluminum castings-improvement by gas cooling, Materialwiss. Werkst. 37 (1) (2006) 116e121.
B. Xiao et al. / Applied Thermal Engineering 30 (2010) 2815e2821 [6] Z. Li, R.V. Grandhi, R. Srinivasan, Distortion minimization during gas quenching process, J. Mater. Process. Technol. 172 (2) (2006) 249e257. [7] H. Li, G. Zhao, S. Niu, Y. Luan, Technologic parameter optimization of gas quenching process using response surface method, Comput. Mater. Sci. 38 (4) (2007) 561e570. [8] B. Xiao, G. Wang, Y. Rong, Hardenability and distortion control in high pressure hydrogen quenching, Int. J. Manuf. Res., in press. [9] B. Sahin, Y. Ust, I. Teke, H.H. Erdem, Performance analysis and optimization of heat exchangers: a new thermoeconomic approach, Appl. Therm. Eng. 30 (2e3) (2010) 104e109. [10] G.N. Xie, B. Sunden, Q.W. Wang, Optimization of compact heat exchangers by a genetic algorithm, Appl. Therm. Eng. 28 (8e9) (2008) 895e906. [11] Z. Luo, H. Cho, X. Luo, K. Cho, System thermal analysis for mobile phone, Appl. Therm. Eng. 28 (14e15) (2008) 1889e1895. [12] J.P. Holman, Heat Transfer, ninth ed. McGraw-Hill, New York, 2002. [13] M. Sedighi, C.A. McMahon, The influence of quenchant agitation on the heat transfer coefficient and residual stress development in the quenching of steels, Proc. Inst. Mech. Eng. Pt. B: J. Eng. Manuf. 214 (7) (2000) 555e567. [14] A. Sugianto, M. Narazaki, M. Kogawara, A. Shirayori, A comparative study on determination method of heat transfer coefficient using inverse heat transfer and iterative modification, J. Mater. Process. Technol. 209 (10) (2008) 4627e4632. [15] M. Maniruzzaman, J. Chaves, C. McGee, S. Ma, R. Sisson Jr., CHTE quench probe system: a new quenchant characterization system, in: 5th International Conference on Frontiers of Design and Manufacturing (ICFDM), 2002, pp. 619e625.
2821
[16] D.K. Funatani, M. Narazaki, M. Tanaka, Comparisons of probe design and cooling curve analysis methods, in: 19th ASM Heat Treating Society Conference Proceedings Including Steel Heat Treating in the New Millennium, 1999, pp. 255e263. [17] G.E. Totten, C.E. Bates, N.A. Clinton, Handbook of Quenchants and Quenching Technology. ASM International, Materials Park, OH, 1993. [18] N.I. Kobasko, A.A. Moskalenko, G.E. Totten, G.M. Webster, Experimental determination of the first and second critical heat flux densities and quench process characterization, J. Mater. Eng. Perform. 6 (1) (1997) 93e101. [19] C. Laumen, S.I. Midea, T. Lubben, F. Hoffman, P. Mayr, Measured heat transfer coefficients by using hydrogen as quenchant in comparison with helium and nitrogen, Accelerated Cooling/Direct Quenching of Steels, 1997, pp. 199e206. [20] M. Li, J.E. Allison, Determination of thermal boundary conditions for the casting and quenching process with the optimization tool OptCast, Metall. Mater. Trans. B 38B (4) (2007) 567e574. [21] G. Wang, B. Xiao, L. Zhang, K. Rong, Q. Wang, Heat transfer coefficient study for air quenching by experiment and CFD modeling, in: 137th Annual Meeting & Exhibition. Materials Characterization, Computation and Modeling, vol. 2, 2008, pp. 271e276. [22] O. Dahlberg, Z. Wu, Feasibility study of three dimensional quenching analysis, CD-adapco Technical Report to GM. 071-186-003, 2008. [23] Q. Wang, B. Xiao, G. Wang, Y. Rong, R.D.J. Sisson, Systems and methods for predicting heat transfer coefficients during quenching (Pending Patent), 2010. [24] B. Xiao, G. Wang, Q. Wang, M. Maniruzzaman, R.D. Sisson Jr., Y. Rong, An experimental study of heat transfer during forced air convection, J. Mater. Eng. Perform., in press, doi: 10.1007/s11665-010-9745-7.