Turbulence modeling of solidification process during continuous casting by Anurag Mahajan A thesis submitted in partial fulfillment of the requirements for the degree of Master of Science in Mechanical Engineering Montana State University © Copyright by Anurag Mahajan (2002) Abstract: Over the last two decades, there has been a tremendous increase in the manufacturing of semi finished products such as steel/aluminum billets, seamless pipes, sintered bronze bearings, copper/aluminum bar stock, hexbars, etc., by continuous casting processes. Continuous casting of metals involves a phase change coupled with fluid flow and heat transfer and is typically non-linear in nature, due to the presence of a moving boundary between the solid and liquid phase. The present study investigated the effects of turbulence on heat transfer and fluid flow during the solidification process of continuous castings. Pure aluminum was modeled as a Newtonian incompressible fluid with Boussinesq approximation. Two-dimensional turbulent fluid flow with submerged and non-submerged entry nozzles was modeled using Prandtl mixing length theory and the average specific heat method. To simulate turbulent flow, effective viscosity and thermal conductivity was evaluated using the Prandtl mixing length theory. In the phase change zone, the latent heat was absorbed into the specific heat of the material. Specific heat evaluation within the two-phase region was done using the average specific heat method. The governing equations were coupled and discretized using the Galerkin technique. A series of simulations were performed to study the effects of inlet speed, superheat, mold and post mold cooling rates on the process characteristics and were presented in terms of viscosity contours, location and slope of the solidification front, axial velocity profiles, velocity vectors, local heat flux, shell thickness and solid fraction. The effects of mold length were investigated to establish remedies for breakout due to higher inlet speed and superheat. The effects of the submergence depth were established by comparing the submerged and non-submerged entry nozzle continuous casting processes. The current study reveals that inlet speed, mold cooling rate and mold length has the strongest effect on the quality of the cast product. The effect of the superheat on the slope and therefore the quality of the end product is negligible. The results obtained indicate that the submergence depth is desirable to reduce the turbulence level at the free surface and to reduce severe casting defects such as slag/flux inclusions in the final product.  TURBULANCE MODELING OF SOLIDIFICATION PROCESS DURING CONTINUOUS CASTING by AnuiagMahajan A thesis submitted in partial MfiUment of the requirements for the degree of Master of Science in Mechanical Engineering Montana State University Bozeman, Montana May 2002 riri# WTlqcI APPROVAL of a thesis submitted by Anurag Mahajan This thesis has been read by each member of the thesis committee and has been found to be satisfactory regarding content, English usage, format, citations bibhograp ic style, and consistency, and is ready for submission to the College of Graduate Studies. Dr. Ruhul Amin (Signature) Date Approved for the Department of Mechanical and Industrial Engineering Dr. Vic Cundy (Signature) Date Approved for the College of Graduate Studies DateDr. Bruce R. McLeod (Signature) IH STATEMENT OF PERMISSION TO USE In presenting tins thesis in partial fulfillment of the requirements for a master s degree at Montana State University, I agree that the Library shall make it available to borrowers under rules of the Library. If I have indicated my intention to copyright this thesis by including a copyright notice page, copying is allowable only for scholarly purposes, consistent with “fair use” as prescribed in the U.S. Copyright Law. Request for the permission for extended quotation from or reproduction of this thesis in whole or in parts may be granted only by the copyright holder. Signature Date La ACKNOWLEDGEMENTS I would like to thank Dr. Ruhul Amin for his guidance in my research and thesis work. I would also like to thank Dr. Alan George and Dr. Thomas Reihman for then work as committee members. I would also like to thank the Sandia National Laboratories for the original version of the computer code used in the present work and to David Grief for the modified version that I used as a starting point for the present thesis work. I also acknowledge the Department of Mechanical and Industrial Engineering and my parents for providing me with the financial assistance. Finally, I would like to thank the staff of the Department of Mechanical and Industrial Engineering and my fellow graduate students for their support in make my work successful. TABLE OF CONTENTS LIST OF TABLES...................... ................................................................... LIST OF FIGURES......................................................................................... NOMENCLATURE........................................................................................ ABSTRACT.................................................................................................... L INTRODUCTION................................................................................. Motivation for the Present Research.................................................... Background......................................................................................... 2. BASIC CONCEPTS AND CLASSIFICATION OF TURBULENCE MODELS............................................................................................. Governing Equations.......................................................................... General Turbulence Models.............................................................. Correlations............................................................................. IntegralMethods..................................................................... One Point Closure.................................................................. Two Point Closure................................................................. Large Eddy Simulation.......................................................... One Point Closure Turbulence Models............................................. Zero Equation Models............................................................ One Equation Models............................................................ Two Equation Models............................................................ Theory of Je ts .................................................................................... Equations of Motion for Plane Turbulent Jet......................... Energy Equation..................................................................... Prandtl Mixing Length in Turbulent Jets.............................. 3. PROBLEM FORMULATION......................................................... Description......................................................................................... Governing Equations........................................................................ Continuity Equation.............................................................. Momentum Equation............................................................ Page ... vii . .viii .xiii ..xvi I ,5 .7 ..18 ...18 ...20 ....21 ....21 ...21 ...22 ...22 ....23 ....25 ...27 ...27 ....31 ....33 ....34 ....34 ...36 ....36 ......38 ......38 ......39 Vl Energy Equation................................................................................. Initial and Boundary Condition .................................................................... Non-Dimensionalization of Governing Equations and Associated Boundary Conditions...................................................................................................... Non-Dimensional Boundary Conditions........................................................ 4. NUMERICAL METHOD Introduction............................................................ Finite Element Formulation.................................. Average Specific Heat and Prandtl Mixing Length Grid Independence................................................ Code Validation.................................................... 47 48 52 57 59 5. 6. RESULTS AND DISCUSSION. Effects of Withdrawal Speed...................................... Effects of Superheat.................................................... Effects of Mold Cooling Rate.................................... Effects of Post Mold Cooling Rate............................ Effects of Mold Length.................................... ......... Effects of Turbulence from Submerged Entry Nozzle ..65 ..81 ..94 105 113 118 CONCLUSIONS AND RECOMMENDATIONS Conclusions......................................................... Recommendation................................ -.............. REFERENCES CITED.............. '•...................... APPENDIX A : Sample Calculations.................. 145 .145 .158 .149 .154 V ll LIST OF TABLES Table Page ...421. Thermophysical properties of aLurnmum......... 2. Values of CC parameters used by other authors 157 1. 2 . 3. 4. 5. 6 . 7 . 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. Figure viii LIST OF FIGURES Page A typical continuous casting process layout..................... ................................. Basic classification of one point closure models................................................ Definition sketch of plane turbulent jet.............................................................. Prandtl mixing length description...................................................................... Computational domain in terms of dimensional parameters.............................. Typical 9-noded iso-parametric element............................................................ Typical part of the finite element mesh............................................................. Grid independence results- Temperature distribution for Pe=6.0, O0 =2.0, Bip = 0.1, Bim=OA.......................................... .................................................. Comparison with the analytical results of Siegel (1984).................................. Physical model of Wolf and Viskanta (1988).................................................... Comparison with experimental data by Wolf and Viskanta (1988).................. Effects of Fe on turbulent viscosity contours for ®0=1.5., Bim=OA3 Bip=OT.... Effects of Pe on solidification front, ©0=1.5, Bim=OA, Bip=O.!....................... Effects of Pe on velocity vector field for ®0=1.5, Bim=OA, Bip=O. I ................ Axial velocity profiles for G0=LS, Bim=OA, Bip=O. I ...................................... Effects of Pe on shell thickness for ©o=1.5, Bim=OA, Bip=O. I ......................... Effects of Pe on solid fraction Q0=LS, Bim=OA, Bip=OT................................ Effects of Pe on shell and core temperature distribution for ©0=1.5, Bim-OA, Bip=OT ........................................................................... ............................... 2 .23 32 .34 .37 .48 .56 .58 ..60 ..61 .62 .67 .69 ..71 ..72 ...73 ...75 ...78 IX 19. 20 . 21 . 22. 23 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. Figure Effects of Pe on local heat flux for ©0=1.5, Bim=OA Bip=OT............................79 Effects of Pe on total heat flux for O0=TA Bim=OA Bip=OT............................ 79 Effects of Pe on percentage heat removed in the mold for O0= IA Bim=OA ............................................................................................................................. *0 Effect of superheat on turbulent viscosity for Pe=3.5, Bim=OT, Bip=O. 1.............84 Effect of superheat on solidification front (a) © o= 1.2, 1.5,2.0, Pe=3.5 Bim=OT, Bip =OT (b) Oo= 1.2, 1.5, 2.0, Pe=6.0 Bim=OT, Bip =OT.................................................... 85 Effects of superheat on local heat flux (a) Pe=3.5 Bim=OT, Bip =OT (b) Pe=6.0 Bim=OT, Bip =OT.............................................................................. Effects of superheat on temperature distribution for Pe=3.5, Bim=OT, Bip=OT....88 Effects of superheat on temperature distribution for Pe=6.0, Bim=OT, Bip=OT....89 Effects of superheat on average heat flux for Bim=OT, Bip =OT......................... 90 Effects of superheat on shell thickness for Pe=3.5, Bim=OT, Bip=O. 1................. 91 Effects of superheat on solidification fraction for Bim=OT, Bip =OT.................. 91 Effects of superheat on percentage heat extracted in the mold for Bim=OT, Bip=OT................................................................................................................ 93 Effect of Bim on turbulent viscosity contours for Pe=3.5, © o= 1.5, Bip=OTO.......96 Effect OfBim on solidification front (a) Pe=3.5, © o= lA Bip =OTO (b) Pe=6.0, @o=2.0, Bip =0.20................................................................................................... Effects of mold cooling rate on local heat flux for Pe=6.0, © o=2.0, Bip =OT......99 Effects of mold cooling on shell thickness for Pe=3.5, Q0=I .5, Bip=O. 1.......... 100 LIST OF FIGURES-continued Page 34. LIST OF FIGUKES-continued Figure Page 35. 36. 37. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48. 49. 50. 51. Effects OfBim on shell and core temperature distribution for Pe - 6.0, ©o=2.05Bip=O. I ................................................................................................ Effects OfBim on shell and core temperature distribution for Pe = 3.5, 0 o=T.5,Bip=O.l................................................................................................ Effects OfBim on local heat flux for Pe = 3.5, ©0=1.5, Bip=O. I ...................... Effects OfBim on local heat flux for Pe = 6.0, ©o=2.0, Bip=O. I ...................... Effects OfBim on percentage heat removed for Pe = 3.5, ©0=1.5, Bip=O l..... Effects OfBim on percentage heat removed for Pe = 6.0, ©o=2.0, Bip=O. I ..... Effects of Bip on turbulent viscosity for Pe=6.0, ©0=2.0, Bim =0.4................ Effects of Bip on solidification front.............................................................. Effects of Bip on local heat flux..................................................................... Effects of Bip on temperature distribution for Pe=6.0, ©0=2.0, Bim =0.4...... Effects of Bip on temperature distribution for Pe=6.1, ©0=1.2, Bim =0.25 Effects of mold length on the solidification front, Pe=6.0, ©o=2.0, Bim=0.3, Bip=O l.......................................................................................................... Effects ofLm* on the local heat flux, Pe = 6.0, @0=2.0, Bim= 0.3, Bip = 01. Effects OfLm* on the percentage heat removed for Pe = 6.0, @0=2.0, Bim= 0.3,Bip = 0 1 .......................................................................................... Comparison of Lm* Vs Bim on the percentage heat extracted in the mold.... SEN computational domain in terms of dimensional parameters................. Effects of ys* on turbulent viscosity contours for Pe=3.5, ©0 =1.5, Bim=OA, Bip=O l......................... - ............................................................... 101 .102 103 103 104 ■ 104 107 .108 .110 111 .112 114 .116 .117 .117 .120 ...121 Xl 52. 53. 54. 55. 56. 57. 58. 59. 60. 61. 62. 63. 64. 65. 66. Figure LIST OF FIGURES-continued Page Effects of ys* on velocity vector field for O0 = 1.5, Bim = 0.4, Bip= OT..... SEN Effects of Pe on solidification front, O0=LS, Bim=OA, Bip=OT.......... Effect of Pe on local heat Auxfor ys* =1, O0==LS, Bim=OA, Bip=OT.......... Variation of heat flux with Pe for SEN andNSEN, O0=LS, Bim=OA, Bip=OT.......................................................................................................... Effect of Pe on percentage heat removed in the mold for SEN and NSEN, ©0=1.5, Bim=OAJBip=OT............................................................................. Effect of Pe on shell thickness for ys* =1, O0= I .5, Bim=OA, Bip=OT.......... Effect of Pe on solid fraction for ys* =1, O0=LS, Bim=OA, Bip=OT............ Effects of©0 on solidification front for SEN andNSEN, Pe=3.5, Bim=OA, Bip=OT.......................................................................................... Effect of O0 on local heat flux for ys* =1, Pe=3.5, Bim=OA, Bip=O.I .......... Variation of heat flux with Pe for SEN and NSEN, Pe-3.5, Bim-OA, Bip=OT......................................................................................................... Effect OfO0 on percentage heat removed in the mold for SEN andNSEN, Pe=3.5, Bim=OA, Bip=OT................................................................;........... Effect of 0 O on shell thickness, ys* =1, Pe=3.5, Bim=OA, Bip=O. I ............. Effect OfBim on solidification front, SEN andNSEN, Pe=3.5,0O =1.5, Bip=OT.......................... .............................................................................. Variation of heat flux with Bim for SEN andNSEN,Pe=3.5,0 0=L5, Bip=OT....................................................................................................... Heat removed in the mold with Bim for SEN and NSEN,Pe—3.5,0 0- 1.5, Bip=OT............................. - ........................................................................ 122 124 .126 126 128 .129 .130 .132 .133 .134 .135 ..136 .138 .139 ...140 XU LIST OF FIGUKES-continued Figure Page 67. Effect of Bip on solidification front for SEN and NSEN, Pe=6.1, G0 =1.2, ....141 68. Effects of Bip on local heat flux, ys* =1, Pe=6.1, ©0=1.2, Bim -0.25.................143 69. Variation of heat flux with Bim for SEN andNSEN, Pe=3.5, ®0-1.5, Bip=O. I .......................................................................................................... ....143 Effects of Bip on temperature distribution for Pe=6.0, ©0-2.0, Bim -0.470. 144 XUl Symbol Ar Bi Cp g h H k L Im P Pe Pix Q q*" Ste t T Y* u Uo Description. Aspect ratio of cast material, LAV Biot number, h w / k Specific heat Acceleration due to gravity Convective heat transfer coefficient Latent heat of fusion Thermal conductivity Length of cast material Prandtl mixing length Pressure Peclet number, U0 Wn / ocs Turbulent Prandtl number Total dimensionless heat flux Dimensionless local heat flux, equation 3.17 Stefan number, Cs (Tf - T00) / H Time Temperature Dimensionless temperature, TZT0 X component of velocity Inlet speed NOMENCLATURE XlV NOMENCLATURE - continued Us Withdrawal speed V Volume of the solidified shell Vo Total volume of the computational domain V Y component Of velocity W Half width of the cast material, see Fig. 5 Wn Half width of the entry nozzle, see Fig. 5 %,y Spatial coordinates x*, y* Non-dimensional spatial coordinates, x/W, y/W ys* Non-dimensional submergence depth Greek Symbols a Thermal diffusivity, k /p C P Coefficient of volumetric expansion 8 Dirac-delta function M- Molecular viscosity M-T Turbulent viscosity M-T* Dimensionless turbulent viscosity, M-T / M- P Density ©o Superheat, (T0 - T00) / (Ts - T00) a StefanBoltzman constant, 5.67E-8 WZm2K4 NOMENCLATURE - continued Subscripts and Superscripts Effective value Phase change Liquid Mold property Nozzle property Inlet value Post mold or spray region property Sohdproperty Ambient value Non-dimensional value XVl ABSTRACT Over the last two decades, there has been a tremendous increase in the manufacturing of semi finished products such as steel/aluminum billets, seamless pipes, sintered bronze bearings, copper/aluminum bar stock, hexbars, etc., by continuous casting processes. Continuous casting of metals involves a phase change coupled with fluid flow and heat transfer and is typically non-linear in nature, due to the presence of a moving boundary between the solid and liquid phase. The present study investigated the effects of turbulence on heat transfer and fluid flow during the solidification process of continuous castings. Pure aluminum was modeled as a Newtonian incompressible fluid with Boussinesq approximation. Two-dimensional turbulent fluid flow with submerged and non-submerged entry nozzles was modeled using Prandtl mixing length theory and the average specific heat method. To simulate turbulent flow, effective viscosity and thermal conductivity was evaluated using the Prandtl mixing length theory. In the phase change zone, the latent heat was absorbed into the specific heat of the material. Specific heat evaluation within the two-phase region was done using the average specific heat method. The governing equations were coupled and discretized using the Galerkin technique. A series of simulations were performed to study the effects of inlet speed, superheat, mold and post mold cooling rates on the process characteristics and were presented in terms of viscosity contours, location and slope of the solidification front, axial velocity profiles, velocity vectors, local heat flux, shell thickness and solid fraction. The effects of mold length were investigated to establish remedies for breakout due to higher inlet speed and superheat. The effects of the submergence depth were established by comparing the submerged and non-submerged entry nozzle continuous casting processes. The current study reveals that inlet speed, mold cooling rate and mold length has the strongest effect on the quality of the cast product. The effect of the superheat on the slope and therefore the quality of the end product is negligible. The results obtained indicate that the submergence depth is desirable to reduce the turbulence level at the free surface and to reduce severe casting defects such as slag/flux inclusions in the final product. ICHAPTER I INTRODUCTION Over the last two decades, there has been a tremendous increase in the manufacturing of semi-finished products, such as steel/alummum billets, seamless pipes, sintered bronze bearings, aluminum bronze bearing, copper/alummum bar stock, hexbars, etc., by continuous casting (CC) processes. In a typical vertical continuous casting process, superheated metal flows from a ladle into an open-ended, water-cooled mold through a tundish and a ceramic nozzle. The superheated metal is protected from air exposure and oxidation by slug covers in each stage. Once in the mold, molten metal solidifies against the mold walls to form a solidified shell. Powdered mold flux is used to provide thermal/chemical insulation along with a uniform heat transfer and lubrication to prevent metal sticking to the mold. Figure I represents a typical CC process setup. The process of CC has many benefits over the conventional casting process. High volume productivity, energy savings, ability to cast pure metals and alloys into different shapes and sizes with ease of post mold working, such as cold/hot rolling or forging on the CC production lines, are some of its well known benefits. The scopes of high automation, due to elimination of intermediate processes such as ingot stripping and manual de-molding, are the additional benefits of the process. On the other hand, inability to cast products with close finished geometrical tolerances, complex configurations, and internal defects, such as longitudinal or transverse cracks, are some of its limitations. Lack of control on the process variables can result in a very serious problem known as 2“breakout.” Breakout occurs if the solidifying metal (due to insufficient shell thickness) cannot withstand the ferrostatic pressure, resulting in rupture during the casting process. Optimization of process parameters and establishing good control is therefore quite essential to reduce the scrap and provide a high performance finished product. Ladle , Molten \ Metal Tundish Submerged entry nozzle Liquid pool Solidification front Post mold spray cooling Post Processing Rollers Figure I . A typical continuous casting process layout. 3Continuous casting of metals involves a phase change coupled with fluid flow and heat transfer and is typically non-linear in nature, due to the presence of a moving boundary between the solid and liquid phases. The momentum of the pouring stream in conjunction with thermal buoyancy generates a turbulent motion of superheated metal in the solidifying pool. This turbulent motion, along with the superheat of the pouring stream, mold and post mold cooling rates, influences the critical process performance characteristics, such as rate of solidification, location and slope of the solidification front, heat and mass transfer and, therefore, the quality of the end product. Consequently, knowledge of the fluid flow and heat transfer in CC operation is a prerequisite for effective process analysis. Simulation of the CC process has been approached both analytically and numerically. Analytical techniques include heat balance integration, Goodman (1958), integral profile methods. Hills and Moore (1967), embedding method, (Boley, 1974), variation technique, Yeh and Chung (1975), isotherm migration. Crank and Gupta (1975), source and sink method, Keung (1980) and Cauchy boundary value method, Siegel (1984). Although analytical methods provide approximate solutions, they have the advantage of simpler mathematical formulations. Numerical modeling of the CC process primarily uses finite element and finite difference methods. Two different computational mesh variations are commonly applied to the problems; namely, time variant mesh (commonly known as moving grid technique) and fixed mesh technique. The moving grid technique is comparatively complicated and traces the solid-liquid interface with time and is limited to simpler geometries. Numerical 4simulations for temperature history and/or location of the solidification front lead to over/under predictions of the results, as well as numerical oscillations. The fixed grid technique, on the other hand, is considerably simpler; however, it requires accounting of latent heat into material properties, such as specific heat or enthalpy. The fixed grid technique typically uses equivalent heat capacity and enthalpy models. Hsiao (1985) developed an efficient finite difference scheme using the equivalent heat capacity method. Lee and Chiou (1995) further developed Hsiao’s method into an average specific heat method for finite element method to simulate the properties of the elements undergoing the phase change. Amin and Greif (1999) used this method in conjugate heat transfer in a continuously moving metal during solidification. Although numerical analysis of coupled fluid and heat transfer in the CC process dates back almost two decades [Asai and Szekely (1975)], only a few studies have been aimed at coupled turbulent flow and heat transfer with solidification in the pure metal CC process. The stream function vorticity method, Asai and Szekely (1975), the standard high Reynold number (Re) k -s method, Thomas et al. (1990), the low Re k-s turbulence model considering thermocapillary and buoyancy effects, Abutalbi et al. (1995), and the low Re k-s turbulence model using popular enthalpy-porosity technique, Seyedein and Hasan (1997) are some of the techniques adopted to simulate turbulent CC of steel systems. In the present study, the Prandtl mixing length method for turbulent jet flows in conjunction with the average specific heat method was used to simulate the coupled turbulent flow and heat transfer with solidification in the pure aluminum CC process. 5Motivation for Present Research In recent time, commercializing the newly developed processes or improving the existing processes has been the main focus of the research community. However, commercial realization or improving an existing process demands extensive knowledge and understanding of the processes and is rarely achieved without establishing a strong experimental database. Extending the experience derived from existing technology to establish various aspects of a new process does not serve the purpose very often. Therefore, a great deal of interest is shown in developing advanced numerical modeling techniques to establish qualitative and quantitative aspects of various processes. The process of continuous casting (CC) is one such process that has been the focus of both industry and the research communities for the last two decades. Production of various ferrous and non-ferrous products by CC in the United States of America alone has increased by about 20% in last two years. For the same reasons, there has been a tremendous investment of funds for research and development of the CC process at both industrial and academic levels. Continuous casting of metals involves a phase change coupled with fluid flow and heat transfer and is typically non-linear in nature, due to the presence of a moving solid and liquid interface. Typical CC parameters include the withdrawal speed, superheat of the pouring stream, mold/post mold cooling rates and submergence depth of the nozzle. These parameters influence the critical process performance characteristics, such as rate of solidification, location and slope of the solidification front, heat and mass transfer and therefore, the quality of the end product. Consequently, prediction of these parameters in 6the CC operation is highly desirable for effective process analyses and was the motivation for the present numerical modeling. Although the CC process is turbulent in nature, only a few investigations are done on the turbulent aspects of the process. Moreover a significant amount of information is available on the numerical analyses of alloys, especially steel; yet few reports can be found for pure metals. The studies of alloys demand extensive computational time and lacks the simplicity with which simple CC process of pure metals can be evaluated. Thus the process of CC of pure metals still requires simple and efficient models to simulate the phenomenon. Therefore, simplified formulations must be explored and compared for adequacy with the available complex techniques. The simplicity and ease with which the zero equation model, based on the Prandtl mixing length theory, can be used Was a major motivation for the current work. The finite element code NACHOS II developed by Gartling (1987) and modified by Amin and Grief (1999) to incorporate the average specific heat method [Lee and Chiou (1995)] was used in the present analysis. The ease of maintaining larger time steps with accuracy using this method was another motivation for the present work. Details of the Prandtl mbrimr length and average specific heat methods are discussed in the subsequent chapters. 7Backgromid Many industrial applications of heat conduction, such as continuous castings, welding, extrusion, anti-icing etc., involve solidification or melting with phase change. Heat conduction problems involving solidification and melting are also regarded as the Stefan problems. Generally, the Stefan problems are non-linear, due to the presence of the moving boundary whose position is not known a priori. Pioneering work on such problems dates back to Evans et al. (1950), in which they investigated exact solutions of such problems. Since then, there has been a tremendous focus on numerical solution of the Stefan problems. Rubinstein (1971) proposed a complex mathematical model to solve the Stefan problem. This model was based on an interface that moves either into the solid region (melting) or into the liquid region (solidification), in accordance with the relative magnitude of temperature gradients on either sides of the interface. Due to the complex nature of the model, the applicability of this model was limited to one-dimensional problems only. Bonacina (1973) proposed a three-time level implicit fixed mesh approach, which was unconditionally stable and convergent. This scheme was based on an analytical approach consisting of an approximation of the latent heat by a large heat capacity over a small range of temperature. In this approach, temperature dependent coefficients were evaluated at intermediate time level; therefore, the complication of solving a set of non­ linear equations at every time step was avoided. However, this method was limited to 8problems with a small interval of temperature, over which solidification takes place and also deals with a drawback of severe nonlinearity in material property. Bonnerot and Janet (1977) used a time-variant mesh approach to solve a phase change problem involving heat conduction only. In this explicit scheme, the mesh was deformed dynamically so that a continuous track of the solid/liquid front location could be achieved. This method offered high accuracy; however, it was limited to simpler problems and geometries, with the drawback of large computational time. Gartling (1977) approached both the free and forced convection phase change problems by using temperature dependent properties based on a finite element method. In this method, field equations were discretized using Galerkin method. In the current study, NACHOS H developed by Gartling (1987) and modified by Amin and Greif (1999) was used as a starting point. Morgan et al. (1978) simplified and improved the previously mentioned work of Bonacina et al. (1973) by reconsidering the latent heat effect accompanied by a phase change. In this approach, the latent heat effect accompanied by a phase change was approximated by integration of the terms involving heat capacity. Therefore, it required an accurate evaluation of heat capacity at integration points. Their results led to an important finding that the time step size should not exceed the temperature interval over which solidification takes place. Siegel (1978) used a conformal mapping technique to approach the continuous casting (CC) process. In this approach, the ingot interface shape was obtained either for the liquid metal at the solidification temperature or for the solid and the liquid interface 9by transferring a uniform heat flux to the interface. The complexity of the conformal transformations required in the mapping method made it difficult to treat other types of interface heating conditions. Siegel (1984) approached the CC process by assuming the spatial variation of the heating of the interface by the liquid phase to be a known function. In his two-region model, conduction was assumed to be the only mode of heat transfer. The results were obtained by a Cauchy boundary method. In the present work, laminar code validation was done against results obtained by this method. Hsiao (1985) approached phase change problems by using the equivalent heat capacity method. The latent heat of fusion was taken into account by using a linear interpolation of the nodal temperatures. This scheme was independent of solidification interval and, hence, could be used for both pure metals as well as alloys. Wolf and Viskanta (1988) performed a combined experimental and numerical study of solidification of a pure metal in the presence of liquid superheat. Temperature and position of the interface were used to deduce the importance of natural convection in the liquid pool on the solidification process. Their results indicated the need of a computational scheme that was capable of simultaneously tracking several moving boundaries (such as solid-liquid, solid-gas and liquid-gas) that occurs when a material shrinks upon solidification. The numerical results in the current work were compared with the experimental results of Wolf and Viskanta (1988). Zabaras et al. (1990) used a front tracking finite element method to calculate the temperature and stress field development in solidifying pure metal. In this approach, interface velocity and location were treated as primary variables and rate dependent 10 viscoplastic-hypoelastic models were employed to determine the stress state on the interface. Their results indicated that the melt pressure significantly alters the deformation at the early stages of solidification. The focus of study by Jaluria (1992) associated with a continuously moving material undergoing thermal processes was the thermal buoyancy, transient, and forced flow in the ambient medium. Three main approaches were outlined in this work. The first was directed at the heat transfer in the material, assuming convective coefficients at the surface. The second solved for the flow generated by an isothermal-moving surface. The third approach was the combination of the above two that considered the conjugate problem that couples the flow and transport in the fluid with those of the material. Kang and Jaluria (1993) used the enthalpy method, assuming a heat transfer coefficient at the surface of the material to solve I-D two zone and 2-D problems. Their results revealed that the interface shape and the resulting temperature fields were strongly dependent on the withdrawal speed and the cooling rate in the mold. For a small Peclet number, their results were in good agreement with the previously published analytical results. The investigations of Huang et al. (1992) on superheated dissipation in the CC process revealed that most of the heat is dissipated in the mold or just below the mold. Also, the amount of superheat and casting speed were found to be the dominating factors. The investigations on two-dimensional mathematical models, which included the effects of nozzle jet angle and submergence depth on the shell thickness and heat flux, were able to predict many three-dimensional results. 11 Lee and Tzong (1995) presented a modified latent heat method for a binary alloy system by solving species and momentum equations in the liquid and the mushy zones. In this method, an interpolation technique was proposed to trace the interface of the mushy zone and the pure liquid region. The results obtained by this technique revealed that a great amount of latent heat could be released from an area when the eutectic front sweeps through it. Furthermore, the total volume of the mushy zone predicted by this method showed a satisfactory agreement with the experimental results. Saitoh and Sato (1994) adopted a boundary-fixing method (BFM) to handle the moving boundary and verified the forced and natural convection effects in the two- dimensional CC process. In the BFM technique, they considered the arbitrary geometries of both the moving interface and the domain boundary via the change of an independent variable, thereby reducing the original problem into the one-dimensional fixed boundary problem. They also analyzed the effects of control parameters on uniform and nozzle flow models of CC process. Their results suggested that the location of the solidification front was strongly dependent on the Peclet number. Lee and Chiou (1995) developed an average specific heat method to simulate the properties of the elements undergoing the phase change, which was later used by Greif (1998) in the conjugate heat transfer in a continuously moving metal during solidification. This method was found to have little effect of the time step on the average and Tnayimum error. Therefore, larger time steps could be used in this method to save computational time. Insensitivity towards the temperature interval over which the 12 solidification takes place was another advantage of this method. The present study of CC process is based on this method. Although numerical analyses of coupled fluid and heat transfer in the CC process date back almost two decades, only a few studies have been aimed at coupled turbulent flow and heat transfer with solidification in pure metal CC process. Brimacombe et al. (1974) approached the phenomenon of CC castings by the heat transfer aspect of the process alone. In this approach, the authors used an equation similar to a conduction equation with a relatively large thermal conductivity in the liquid pool region to account for possible effects of turbulence on the heat transfer. Asai and Szekely (1975) used a stream function vorticity-based method in conjunction with a simplified turbulent model based on a mixing length approach in the coupled fluid and heat transfer in CC steel systems. While this approach enabled tracer dispersion to be calculated, in this regard, the agreement between measurements and predictions was less satisfactory. Ability to estimate the trajectory of inclusion particles within the mold pool was an attractive feature of the model. This method also established that a region of reverse flow is inherently associated with the use of straight nozzles. The general dimensional model of the CC process proposed by Indyk and Hadden (1979) emphasized the necessity of coupling heat transfer and turbulent fluid flow. Experimental results obtained by the authors indicated the need of slow startup of the process. The experiments established that the surface appearance of the continuously cast billets serves as a good indicator of the internal grain structure and provides an insight into the process of solidification. One of the interesting findings of their results was that 13 at too low casting speeds and temperatures, transverse cracks take the form of open cracks. Examination of the grain growth indicated that, for the greater part of the cycle, the temperature gradient is perpendicular to the axis of the billet. These results signified a rapid heat transfer rate in the radial and virtually none in the longitudinal directions in the bizonal growth section. The experimental study of Nakato et al. (1984) was aimed at the factors affecting the formation of shell and longitudinal cracks. Their experiments indicated that the withdrawal speed, stream velocity from the submerged nozzle, mold taper and mold oscillations were influential factors in the shell formation. They also found that the casting speed and the mold powder are the most influential factors to the mean heat flux in the mold. This work was also aimed at enhancing the uniform shell thickness and crack prevention on the slabs. Thomas et al. (1990) presented a turbulent model for steel by investigating the superheat dissipation on the fluid flow and temperature field in a continuous slab casting system using a standard high Reynold number (Re) k-z method. However, these studies did not couple the process of solidification within the caster. The effects of casting speed, nozzle angle and turbulence on the CC of steel were the main focus of this investigation. The overall flow field was found to be relatively insensitive to the process parameters. Farouk et al. (1992) performed numerical and experimental investigations on a twin - belt (Hazelett) caster. Their coupled heat transfer and fluid flow numerical model used a generalized energy equation that was valid for the solid, liquid and mushy zone in the cast. A zr-e, turbulence model was used to calculate the turbulent viscosity in the melt 14 pool. Their work revealed that the solidification process within the cast acts as a strong turbulence damper. The size of the mushy zone was found to vary with the process parameters. Their results also reflected that heat transfer coefficient variations along the belt might result in reheating within the cast. Thomas (1993) presented a comprehensive system of mathematical models for CC of steel slabs. This work summarized the different models, which calculate the turbulent flow through the bifurcated nozzle, fluid flow through the mold, mass transfer and heat transfer across the mold powder interface. The model calculations were compared with the observations using physical water models (copper mold with water cooling chamber) and measurement on actual operating casters to reproduce many known phenomena in the CC of steel. Aboutalebi et al. (1995) approached the steel CC phenomenon by a low Re k-z turbulence model to analyze turbulent flow, solidification, and evolution of macrosegregation in a continuous billet caster. A controlled volume based finite- difference procedure was employed to solve the conservation equations associated with appropriate boundary conditions. Solidification was modeled based on a continuum model using the enthalpy method. Due to high nonlinearity in the system of equations, a number of techniques were used to accelerate the convergence time. Their results concluded that an increase in the casting speed results in an increase in the length of the recirculation zone, as well as the depth of the liquid pool. The eddy-viscosity contours obtained Aom the results revealed that the turbulence effects on the process are appreciable even below the mold. 15 Choudhary and Mazumdar (1995) reported a coupled fluid flow and heat transfer model based on the cell porosity method and artificial enhanced melt viscosity for CC of steel. In this method, the mushy zone was modeled by enhancing the dynamic viscosity of liquid steel by twenty times in the mushy zone. The results obtained by this method were practically insensitive to the turbulence phenomena in the liquid pool region. Also, the thermal buoyancy phenomena and boundary conditions applied to the heat flow equation at the meniscus were shown to have negligible effects on the predicted results. Seyedein and Hasan (1997) explored the steady state transport phenomena of the coupled turbulent flow, heat transfer and macroscopic solidification in a CC steel caster. This model was based on a generalized transport equation applicable to liquid, mushy and solid zones. The turbulent effect was taken into account by using a Iow-Reynolds number (Re) - K-z turbulence model. An enthalpy-porosity technique was used to model solidification of the mushy zone. The analysis of turbulent viscosity contours indicated wide variations of the turbulent viscosity within the caster. Their results also indicated that models based on artificially enhanced viscosity or thermal conductivity or have uncoupled the solidification process from heat transfer and fluid flow could not provide realistic results. Most of the modem CC machines use submerged-entry nozzles (SEN) to maximize the productivity and quality of the cast product. The type of nozzle and submergence depth are the most critical features of the submerged-entry nozzles that affect the process of CC. Therefore, researchers have focused their studies on submerged-entry nozzles and the effect of submergence depth on CC process. Asai and 16 Szekely (1975) used a straight nozzle in their study and established that the region of reverse flow is inherently associated with the use of straight nozzles. They also suggested that such flow patterns tend to encourage the entrapment of the inclusions by advancing dendrite arms. Starting from the early 1990’s SEN CC process has once again been a main focus of numerical investigations. Hershey et al. (1993) used a standard at-s turbulence method to study the turbulent flow through bifurcated submerged-entry nozzles (SEN) used in CC of steel. Both two-dimensional and three-dimensional simulations were performed. Their results suggested that, in a SEN CC process, the flow pattern is relatively insensitive to the submergence depth The flow pattern was found to be more sensitive to the high velocity jet leaving the nozzle. Their results also indicated that a separate numerical model of the nozzle should provide reasonable boundary conditions for the flow pattern in the mold. Najjar et al. (1995) performed a numerical study of steady turbulent flow through bifurcated nozzles in CC of steel. Two and three dimensional simulations of steady turbulent at-s flow in bifurcated nozzles indicated that high surface turbulence and quality problems related to higher casting speeds can be avoided by an angled SEN port. The analysis also indicated that turbulence levels in the jet are higher with higher casting speed and smaller port. Aboutalebi et al. (1995) proposed a low Re at- s turbulence model considering thermocapillary and buoyancy effects to analyze a SEN type steel CC process. A bifurcated nozzle with a zero port angle was used in this analysis. The results obtained 17 from this analysis reported the presence of a small vortex near the upper meniscus, while a weaker but larger vortex existed in the lower region of the mold, a condition typically associated with a bifurcated nozzle. The shell thickness predicted by this model was in agreement with the experimental results. The previously mentioned work of Seyedein and Hasan (1997) also included the effects of submergence depth on the flow, temperature and solidification profiles. Their study revealed that the submergence depth has a significant effect on the quality of the cast product. Although higher casting rates were observed with lower submergence depth, unwanted entrapment of the mold flux was found to be a potential danger to the quality of the cast product. Higher submergence depth, resulting in reduced turbulence level in the upper part of the caster, was desirable for a better quality product. Modem commercial realization of the continuous casting process has been attracting both industrial and academic sectors to further explore the process. The numerical simulation techniques have been focused more on the practical opportunities and their solutions for improving the quality and the productivity of the industrial machines. The effects of turbulence, effects of submerged entry nozzle, thermal and mechanical stress analysis, mechanism of crack formation, intelligent molds, and microstmcture analysis have been the main focus of ongoing research worldwide. 18 CHAPTER 2 BASIC CONCEPTS AND CLASSIFICATION OF TURBULENCE MODELS Turbulence is a highly complex physical phenomenon that exists in flow problems of scientific and engineering interest. Due to the complexities of this phenomenon, universal mathematical models of the phenomenon are not available. The present review is limited to turbulence models from an engineering point of view. This means that the detailed resolution of a turbulent flow will be eliminated with some type of average flow analysis. Governing Equations In principle, the Navier-Stokes equations are capable of fully analyzing the turbulent flow. However, due to the inherently unsteady, three-dimensional, and wide length scale spectrum nature of the turbulent flow, direct numerical solution of Navier- Stokes equations is not feasible in most cases. The standard alternative approaches to direct numerical simulation involve the solution to some form of averaged Navier-Stokes equations. In most flow problems, it is the mean flow that is of most interest. The turbulent fluctuations are significant only in the evolution of the mean flow. By performing a suitable average on the instantaneous Navier-Stokes equations, a standard mean flow problem can be derived. The turbulent effects in such problems are expressed in a few terms that can be modeled. 19 To understand this approach, let us express the instantaneous fluid velocity and pressure, fields as the sum of a mean and fluctuating component as follows: Mz=IZz-+Mz- (2.1) P - P + P (2.2) where the upper case represents a mean quantity and superscript prime denotes a fluctuation. Substitution of these definitions in the standard Navier-Stokes equations results in the following equations: SUi Qxl I 8^ U eu^ v J dx a? a ---------------------- 1- - 'j j dxt dx. ^azz <%z^ ^axz a^ - P U iU :+/% (2.3) (2.4) The above equations are in Cartesian components. The average quantity is indicated by the horizontal overbar. Equation (2.3) and (2.4) describe the behavior of the mean fluid velocity and pressure fields. These equations are very similar in form to the laminar equations. The extra term in the equation (2.4) is often termed as Reynolds stress and represents the effects of turbulent velocity fluctuations on the mean flow. In the absence of any additional equations, the appearance of the unknown additional Reynolds stress leads to the well-known closure problem for turbulence. Through appropriate manipulation of the Navier-Stokes equations, six equations for the Reynolds stress variables may be derived. Unfortunately, these derived equations contain additional unknowns in the form of their higher order moments. For example, equations derived for the third-order moments contain fourth-order moments. Cascading the 2 0 unknowns is referred to as the closure problem. In order to achieve a solution to the boundary value problem, the equation derivation must be stopped at some level and a model introduced for the remaining unknowns. A number of turbulence models are discussed in the following section. General Turbulence Models In order to complete the mean flow boundary value problem described by equation (2.3) and (2.4), six additional equations for the Reynolds stresses must be specified. This specification constitutes a turbulent model. A number of turbulence models, ranging in complexity from simple algebraic expressions to descriptions involving multiple, non-linear partial differential equations, are available. Unfortunately, there is no universal method to classify such models. Here, the classification scheme groups the turbulence models according to the following categories: O Correlations. ❖ Integral methods. ❖ One point closure. ❖ Two point closure. ❖ Large eddy simulation. ❖ Direct numerical simulation. In this section, salient features of the above-mentioned model are discussed briefly. 21 Correlations Correlations are limited to the standard problems that have been extensively verified by experimental methods. Correlations are typically developed for prediction of friction factors and heat transfer coefficients for various types of channel flows. Correlations generally do not play any significant role in practical engineering simulations. t' Integral Methods Integral methods are useful when dealing with boundary layer flows or the problems that have parabolic nature. Typically, in this approach, a standard mean field equation integrated over one coordinate is utilized to reduce the mathematical complexity of the problem. The numerical procedure can easily be simplified and expedited by the use of empirical relationships. However, these models also lack generality and are limited to the problems that have been extensively studied using experimental techniques. These types of models are still in use in industry for repetitive design computations. One Point Closure One point closure models are very popular for computational work. In this approach, the correlation between the fluctuating velocity u / and Uj ’ components always use values taken at the same physical location. Despite this limitation, one point models are the basis of a wide variety of successful models. In the present work, Prandtl mixing length method based on the zero equation model (which is a sub-classification of the one point closure models) was used. This approach is further discussed in the subsequent sections. 22 Two Point Closure Two point closure models remove the limitation of one point closure and introduce the effect of length scale explicitly. This approach is still in the development stage and is restricted to homogeneous and isotropic turbulence. Large Eddy Simulation Large Eddy Simulation (LES), also known as Subgrid Scale Modeling (SGM), is used to simulate the large-scale motions of the flow and only models the very small scale eddies. From theory and experimental approaches, it is known that the smallest scale motions tend to be universal in their behavior. Therefore, LES seems to be the idealistic approach for a broadly applicable turbulence model. However, due to the prime consideration of grid refinement, this approach still remains the most expensive approach. 23 One Point Closure Turbulence Models As mentioned previously, one point closure models are the most popular for computational work. Even within this category, there is further classification of the models as shown in the Figure 2. ONE POINT CLOSURE MODELS REYNOLDS STRESS MODELSEDDY VISCOSITY MODELS (Based on Boussinesq hypothesis) Zero Equation Models Algebraic Stress Models One Equation Models Stress Equation Models Two Equation Models Figure 2. Basic classification of one point closure models 24 The one point closure models are broadly classified as eddy viscosity and Reynolds stress models. The eddy viscosity models are further classified as zero equation, one equation, and two equation models. The Reynolds stress models are classified as algebraic stress and stress equation models, as shown in Figure 2. The Reynolds stress models (RSM), although more sophisticated than the eddy viscosity approach, result in a large system of partial differential equations and empirical parameters. The extensive computational nature of these models, along with a lack of certainty of the parameters, makes the RSM unattractive for general engineering applications. In the present work, the zero equation model was used to simulate the turbulent process of continuous casting (CC). Therefore, the “so-called” eddy viscosity models are discussed in detail here. The eddy viscosity approach is based on the Boussinesq hypothesis. By analogy with the molecular diffusion of the momentum, the Boussinesq hypothesis relates the turbulent momentum transport to the mean velocity field gradients. The Reynolds stresses in equation (2.4) are then expressed as: pu,u, — Pf *- + —- Qxj Qxi (2.5) where juT is the eddy or turbulent viscosity. Unlike the molecular viscosity, jll, which is a fluid property, the eddy viscosity, jut, is a local property of the flow. When the definition of the Reynolds stress in equation (2.5) is applied into the momentum equation, the mean flow equations become: 25 QUi dxt O I # 2%/J ap a - + - dxt dx. {m + Vt) 'j Sx1 (2 6) (2.7) Once the form of the eddy viscosity is specified, the mean flow can be solved in the same manner as a laminar flow, because the equations are the same except for an additional viscosity term. Although the turbulent flow problem now has been reduced into a familiar set of partial differential equations, the eddy viscosity variation with the flow field still needs to be specified. Scaling arguments presented by many researchers show that eddy viscosity, Pt, is proportional to a characteristic eddy velocity, Ue, and eddy length, Ze, which means: Pt a: p Ue Ie (2.8) Based on the determination of pT, the available turbulent models are classified as zero equation, one equation and two equation models. Zero Equation Models These models evaluate the eddy viscosity in (2.8) by an algebraic expression, such as equation (2.9) of ue and 4, and, therefore, are the simplest of all one point closure models. Most of these models are based on PrandtTs mixing length, which specifies 4 as the length scale across which turbulent mixing takes place. The characteristic eddy velocity, Me is given by: - + J 1/2 Bxi ac/, (2.9) 26 Which, upon substitution in (2.8), gives the following: 1/2 - + - y Qxj Dxi 2(7. (2 10) Mixing length, 4, is an empirically determined parameter that is well known for different types of flows, such as jets, wakes, and boundary layers, and can be evaluated via simple formulae. Once the form of mixing length is known, then equation (2.10) allows the eddy viscosity to be derived, and the turbulent model can be completed. A simple algebraic description of the eddy viscosity, such as equation (2.10) makes the finite element implementation of a zero equation turbulence model very simple and straightforward. From equation (2.10), it is evident that the turbulent viscosity depends on velocity gradients, and for the best accuracy, the velocity gradients should be evaluated at the integration points. The dependence of eddy viscosity on the mean flow velocity results in another source of nonlinearity. This, however, is not a critical consideration for most of the iterative solution techniques. In the present work, eddy length is expressed as Prandtl mixing length, Im The characteristic velocity, Ue, is given by: f-T<2Xy A a.A 2 Za-A 2 + 5u + .2xJ + _ayj (Z l l ) Therefore, the turbulent viscosity is expressed as follows: Mt ~ Ph Za-A2 Za-A2 Za-A2 ^a-A2 -rl/2 Su + Sv v2%v + Sv (2.12) where Prandtl mixing length, Im = 0.0165(y) 27 One Equation Models The limitation of zero equation models to simple geometries and the need for complex flow simulations has led to the development of more advance turbulence models. First step towards this direction was to replace the algebraic specification of the mixing length with a more general transport equation. Because the characteristic velocity, %% is proportional to the square root of the turbulent kinetic energy, k, this implies: M, a t " : (2.13) Based on this relationship, the eddy viscosity in (2.8) can be expressed as: « p (2.14) A partial differential equation for k can be derived from Navier-Stokes equation and is expressed as: dk dk> — + U1 — dt 1 dx.\ J J ^/Jt dk^ + pG - p s (2.15) where G is the generation term, sis the turbulent dissipation, and Qk is a constant. Despite the apparent advancement involving the evaluation of k, the model is still dependent on an algebraic specification of the mixing. Like zero equation models, variations in Ie are problem specific and are characterized for only a few types of flows. Two Equation (k-e) Model A natural evolution of the one equation model was the replacement of the algebraic relation for the mixing length with a second transport equation. Dimensional arguments lead to the following proportionality: Z .a— ' (2.16) s 28 From (2.8), (2.13) and (2.16), the following expression for eddy viscosity can be derived: Ic2 or the Kolmogorov - Prandtl relation /uT — C„/> (2.17) This equation relates the eddy viscosity directly to the turbulence variables k and s. The turbulent kinetic energy in the two equation model is given by equation (2.15). P r dk dk X — + Ui — a JUt dk ^ dXj + pG - p s (2-18) and the turbulent dissipation is expressed by an equation of similar form P r \ds TT os — + Uj J j Pt ds_ (2.19) where G is again a shear generation term and Cl, c2, cik and a 8 are empirically derived constants. In order to obtain closure to the turbulence model, certain terms such as G are modeled by relating their exact form to some function of the mean flow. The two equation, k-z model described by equations (2.18) and (2.19) can be used in conjunction with mean flow equations and the description of pr given by the equation (2.17). The set of equations is highly nonlinear and has a strong coupling between the various transport equations. However, the k-z model is far from universal and has a number of weaknesses. It has been observed by a number of investigators [Reddy and Garthng (1994)] that the dissipation equation in k-z model may cause instabilities in the flow simulation 29 that result in poor or nonconvergence of the numerical solution. Another aspect of these instability problems is the prediction of negative values of k and s. This nonphysical occurrence is often attributed to the inaccurate modeling of the source terms for k and s. The k-z model described by equations (2.18) and (2.19) is known as the high Reynolds number model, because its derivation was predicted on a single length scale associated with a fully developed flow. Close to the solid boundaries, the nature of the turbulent flow changes to the low Reynolds number model, where small eddies associated with dissipation become important. Authors like Thomas et al. (1990), Hershey et al. (1993) and Najjar et al. (1995.) Jhave used the high Reynold number k-z model to predict the process of continuous casting, whereas, authors like Aboutalebi (1995) and Seyedein and Hasan (1997) have used low Reynolds number k-z model. Turbulent Heat Transfer The previous discussions have been focused on isothermal turbulent flow only. Most of the models for isothermal flows have a direct counterpart when the energy equation is included in the problem description. Subsequent discussion is focused on heat transfer aspects of the turbulent flow. Extending the assumption in equations (2.1) and (2.2), the instantaneous temperature can be written as: d = % + 9 (2.20) where 0 is the average or the mean temperature, and the primed quantity is the fluctuation. When equation (2.20) and equation (2.1) are substituted into the instantaneous energy equation, the result is an equation of the following form: 30 PCt —+u' j z J J 7 0© k— v. d x , J /OCp%)# (2 21) The averaging process introduces the three new variables, U j Q representing the components of the Reynolds heat flux vector. These quantities can be modeled in several ways, all of which are analogous to the techniques used for Reynolds stresses. Assuming one point closure models, the Reynolds heat flux can be related to the mean temperature field through Boussinesq type hypothesis [Launder (1988)] as follows: „ " T -T k T 5 0 Ut 5 0 pC„u,d = —------■=—------ C, 5x, Pr^ 5%, (2.22) Substituting equation (2.22) into equation (2.21) gives the following form of the energy equation: PC1 ^50 r7 50^ — + U'j j J _d_ dXj k + ^ - Prx 5 0 dXj (2.23) where kx is the eddy or turbulent conductivity, and Prj is the turbulent Prandtl number. Once the form of the eddy conductivity or the turbulent Prandtl number is specified, then the energy equation in (2.23) can be solved in conjunction with the momentum equations. Empirical specification of the turbulent Prandtl number is the most popular approach in the turbulence modeling [Reddy and Gartling (1994)]. ha the present work, the turbulent Prandtl number was considered to be unity as per Aboutalebi et. al. (1995) and White (1991). Clearly, the modeling of the turbulent heat transfer phenomenon is an active area of research, although its progress is dependent on the progress made in turbulent flow modeling. 31 Theory of Jets As mentioned previously, the liquid metal in a CC process enters the mold through a submerged or non-submerged entry nozzle. Therefore, it is very important to understand the behavior of jet flow. Inthe present work, a plane nozzle was used for the analysis purpose and, therefore, this discussion is limited to plane turbulent jets only. Figure 3 shows a schematic representation of a plane turbulent jet with thickness bo and uniform velocity Uq. Experimental observations on the mean turbulent velocity field indicate that, in the axial direction of the jet, the jet flow could be divided into two distinct regions. The first region, close to the nozzle, is commonly known as the flow development region. In this region, as the turbulence penetrates inwards towards the axis of the jet, there is a wedge like region of undiminished mean velocity, U0, as shown in Figure 3b. This wedge is known as the potential core and is surrounded by a mixing layer on top and bottom. The second region is known as the fully developed flow region. In this region, the turbulence penetrates to the axis and, as a result, the potential core disappears. For most of the practical applications such as continuous castings, it is the fully developed flow region, which is of greatest interest. Therefore, the current discussion is mainly focused on the fully developed flow only. In the fully developed region, the transverse distribution of the mean velocity in the x-direction, that is, the variation of u with y at different sections, has the same geometric shape as shown in Figure 3c. At every section, u decreases continuously from a maximum value of um on the axis to a zero value at some distance from the axis. A comparison of dimensionless velocity (u/um) at different sections of length, b (value of y 32 where u is the half of the maximum velocity) is represented in Figure 3b and 3c. The plane turbulent jet theory [Rajaratnam (1976)] indicates that the velocity profiles at different sections, which could be superposed in this manner, are similar. Potential Core Plane of symmetry Fully Developed FlowFlow Development Region (a) Figure 3. Definition sketch of plane turbulent jet 33 Equations of Motion for Plane Turbulent Jet For a typical turbulent plane jet, the Reynolds equations in Cartesian coordinate system is give by [Schlichting, 1968)]: du duu---- i-v— dx dy I dp d2u du'v' du2 ------—+ v- p dx dy I dp dv2 ay dy du dv n dx dy ay (2.24) (2.25) (2.26) Equation (2.25) relates the pressure in the jet to the ambient pressure, Pro, and fluctuating velocity, v'. Accounting for turbulence and behavior of a two-dimensional jet, equation (2.24) reduces to the following form [Rajaratnam, (1976)]: du du I dr dx dy p dy (2.27) where x represents the shear stress. The momentum equation (2.27) in integral form can be represented as: 2 j pu2dy = M0 (2.28) o where, momentum flux Mo is an important physical quantity controlling the behavior of the turbulent jet. It replaces the individual values of the thickness of the jet, 2 bo, and mean velocity, U0. That is, for a given value of Mo, the same jet behavior is obtained for different combinations of bo and Uo- 34 The energy equation of a plane turbulent jet relates the rate of change of kinetic energy to the rate of turbulence production. The integral form of the energy equation is given by [Rajaratnam, (1976)]: Energy Equation d dx (2.29) Prandtl Mixing Length in Turbulent Jets There are three unknowns in the equations of motion (equations [2.26] and [2.27]) namely u, v and x; but the number of equations is only two. Hence, a third equation is required to obtain a complete solution to the equations of motion. For this missing equation, Prandtl mixing length formula is used. The mixing length, lm, is defined in the following fashion. When a fluid lump traveling at its original mean velocity is displaced due to the turbulent motion in the transverse direction from yi to yz, its velocity differs from the surrounding mean velocity at y? by AU. The distance (yz-yi) at which AU is equal to the mean transverse velocity is the Prandtl mixing length lm. Figure 4. Prandtl mixing length description 35 Prandtl mixing length is related to the turbulent shear stress in the following manner [White (1991)]: (2.30) From dimensional consideration, Im can be related to the width of the jet as Im oc b and b oc x [White (1991)]. Therefore, Im oc x or Im = c x, where c is the constant that is established empirically by various researchers. Various authors have reported the following values of c: Forthmann (1934), c=0.0165; Reichardt (1942), c=0.0164; and Zijnen (1958), c=0.0223 and C=O.0234. In the present work, the value of c = 0.0165 was used. The theoretical velocity profile as predicted by Tollmien [Rajaratnam (1976), White (1991)] are in very good agreement with ForthmamTs value of c (c=0.0165). Several authors have analyzed various turbulent phenomena, using this value of c (c=0.0165) and those results are well documented by Rajaratnam (1976) and Pai (1954). 36 CHAPTERS PROBLEM FORMULATION Description Figure 5 represents the computational domain along with the boundary conditions used in the current research. Pure aluminum was modeled as Newtonian incompressible fluid with Boussinesq approximation. Two-dimensional turbulent fluid flow, with a non- submerged entry nozzle, was modeled, using Prandtl mixing length theory and the average heat capacity method. The symmetry of the problem was used to model only one half of the domain for computation purposes. An aspect ratio, Ar of 25, was used in the current work. The various dimensions shown in Figure 5 are Lm=O-IdL, Lp=O.84L, and Wn=O. 2 W. The finite element code NACHOS II, developed by Gartling (1987) and modified by Amin and Grief (1999) to incorporate the average heat capacity method, was used in the analysis. NACHOS II is a general-purpose finite element code used to solve two- dimensional continuity, momentum, and energy equations. The code can be used for both steady state and transient problems. Details of the code can be found in Gartling (1987). In the present work, the code was modified to incorporate turbulent modeling of the CC process. The internal mesh generator of the code was used to generate the grid. Appropriate boundary conditions were employed in terms of convenient parameters. The inlet speed of the liquid metal in the nozzle was expressed in terms of a non-dimensional Peclet number. The superheat (difference between the inlet and solidification temperature) was expressed in non-dimensional form as 0 O. Boundary conditions in the 37 mold and post mold regions were applied in terms of non-dimensional Biot numbers Bira and Bip respectively. Euler’s method with gradually increasing time steps was used to solve the governing equations. T0, U0 -I I - W n Withdrawal speed direction Figure 5. Computational domain in terms of dimensional parameters. 38 The following assumptions were made in the formulation of the turbulent CC process to simplify the governing equations: 1. The geometry is assumed to be two-dimensional. 2. Liquid aluminum behaves as a Newtonian incompressible fluid with Boussinesq approximation. 3. Aluminum is considered to be homogenous and isotropic. 4. Density change with change of the phase is neglected. 5. The latent heat effect is considered through appropriate modifications of specific heat. 6. The flow through the nozzle is considered to be the same as jet flow. 7. The caster is vertical with respect to the gravitational field and the strand curvature effect is neglected. 8. The effects of mold taper and oscillations are ignored. 9. Viscous work and viscous dissipation are neglected. 10. No slip boundary conditions prevail at the walls. Based on the above assumptions, the following governing equations are used for the present research. Continuity Equation du dv dx dy Governing Equations (3.1) 39 The momentum equations (Navier Stokes equations) in x and y directions are given by: Momentum Equation du du du dp ot ox dy dx d2u d2u dv dv dv dpP ---- (* pU---- 1- pV---—------- h U gt ch dy where, pe = p + Pt d2v d2v dx2 dy2 dy2 (3^) (3.3) (3.4) (J-t in equation (3.4) is given by [Reddy and Gartling (1994)]: Pt — pi, f fa ., Y A a. A2 fa . A2 + VuaV du + dv \d x / + dv ^dyy (3.5) Based on the experimental results by Forthmann (1934) for Prandtl mixing length, lm, the following expression for Im was used in equation (3.5): Im = 0.0165(y) (3.6) Energy Equation d r a r dTpLn -----Vu----vv— p dt dx dy - Ir % •% > I__ _ where, k = k + Prr (3.7) (3.8) In the energy equation, density p and specific heat Cp are constant throughout the liquid and solid phases. In the solidification region, specific heat was evaluated using the average heat capacity method. Effective thermal conductivity in equation (3.8) was 40 evaluated using turbulent viscosity and turbulent Prandtl number, PrT. In the present work, the turbulent Prandtl number was considered to be unity as per Aboutalebi et al. (1995) and White (1991). The thermal conductivity, k, in equation (3.8) was defined as follows [Bonacina et al. (1973)]: k = ks f o r r < r / -A r k = k, + ^ A [ r - ( 7 - AT)] for Tf - AT < 7 < 7} + AT k = kt for 7 > 7y + A7 (3.9) where 7f indicates the solidification temperature, and 2AT is the small temperature interval over which solidification takes place. 41 Initial and Boundary Conditions The following initial and boundary conditions in dimensional form were used: I . At the free surface (y = L, 0 < x< W): Inside the pouring stream (0 < x< Wn): V = U0 (3.10) u = 0 (3 11) H Il S3 . (3.12) Outside the pouring stream (Wn < x< W): v = 0 (3 13) u = 0 (3.14) dT/dy = 0 (3.15) At the axis of symmetry (x=0, 0 < y < L): dT/dx = 0 (3 16) u = 0 (3.17) At the exit or outflow boundary (y = 0 , 0 < x < W): V = Us (3.18) u = 0 (3/19) STJdy = 0 (3.20) At the wall or the mold side ( x = W, 0 < y < L): Mold region (x = W, Lm < y < L): u = 0 (3.21) h = hm (Convection) (3.22) 42 Post mold region (x = W, 0 < y < Lp): u = 0 (3.23) h = hp (Convection) (3.24) Note that in Figure 5, x = 0 defines the axis of symmetry, x = W defines the outside edge, Wn defines the nozzle width, y = 0 defines the outlet, y = L defines the inlet plane, Lm defines the mold length and Lp defines the post mold length. In the present work, Lm = 0.16L, Lp = 0.84L, Wn = 0.2W. Initially, the metal was considered to be superheated with initial velocity at the inlet and the withdrawal speed at the exit of the caster. As mentioned previously, the material under investigation was aluminum. The thermophysical properties of aluminum are listed in the following table: Table I. Thermophysical properties of aluminum. Property Liquid aluminum Solid aluminum Density (KgZmj) 25425 2542.5 Specific Heat (J/Kg K) 1080 1076 Viscosity (Kg/m s) 1.3xl0'j - Thermal conductivity (W/m K) 94.03 238 Thermal expansion ( 1/K) 1.2xl0"4 - Latent Heat of Fusion ( J/Kg) 3.95xl05 - Phase change temperature (K) 934.52 933.52 43 Non-Dimensionalization of Governing Equations and Associated Boundary Conditions ' , The fluid flow and thermal energy balance equations, the turbulence model, and the associated boundary conditions summarized above provide a complete description of the mathematical model developed in the present work. The governing equations were solved in terms of primitive variables. However, for the sake of generality, the results are presented in terms of non-dimensional variables shown below. I. Dimensionless space variables x* = x/W (3.25)IIi (3.26) 2. Dimensionless velocity components u* = u/Uo (3.27) v* = v/Uo (3.28) Us* = UsZU0 (3.29) 3. Dimensionless pressure P* = P/(pU02) (3.30) 4. Dimensionless temperature T* = TZT0 (3 31) 5. Time t* =AJ0ZW (3.32) 6. Dimensionless turbulent viscosity Ii (3.33) 44 7. Dimensionless local heat flux The dimensionless local heat flux along the outside edge of the cast aluminum is defined [Choudhary and Majumdar (1994)] as: q ^ a ^(r-c) (3 34) ps'-VJ The above dimensionless parameters result in the following governing equations: d2u d2u* Ck'2 """ du* * du* * du* dp* I — — + M — + V — — ---------— H--------- dt dx dy dx Re dv* * dv* * dv* dt* +u ~d7+v d / dp* I d2v* d2v* dy* Re dx*2 dy*2 + - ! ) / [ / : ar .ar .ar ----------- h l l ----------- h V --------- at ck' ay" I d2T* d2T* Re.Pr LaSc^ + (3-35) (3.36) (3.37) (3.38) The non-dimensional parameters appearing in the governing equations are as follows: Re (Reynolds number) = UoWpZpe (3.39) Pr (Prandtl number) ~ CpZ ke (3.40) Ste (Stefannumber) = Cps(Tf -T 00)ZH (3.41) ©o (Superheat) = (T0 - T«,)Z(TS - T00) (3.42) Pe (Peclet number) = UsWZas (3.43) Bi (Biot number) = h WZk (3.44) 45 Non-Dimensional Boundary Conditions The non-dimensionalized boundary conditions used in the present work are as follows: 1. At the free surface (y* = LAV, 0 < x*< 1): Inside the pouring stream (0 < x* < Wn/W): v* = I u* = 0 T* = I Outside the pouring stream (Wn*< x* < 1): v* = 0 u* = 0 dT*/dy* = 0 2. At the axis of symmetry (x*=0,0 < y* < L/W): dT*!dx* = 0 u* = 0 3. At the exit or outflow boundary ( y* = 0, 0 < x* < 1): v* = 0.2 u* = 0 (3.45) (3.46) (3.47) (3.48) (3.49) (3.50) (3-51) (3 52) (3.53) (3.54) dT*/dy* = 0 (3.55) 46 4. At the wall or the mold side ( x* = 1 ,0 < y* < L/W): Mold region (x* = I, Lm* < y* < L*) u* = 0 Bi = Bim (Biot number) 5. Post mold region (x* = I, 0 < y* < Lp*): u* = 0 Bi = Bip (Biot number) In the present investigation, Lm /L=O. 16, LpZL=O.84, L/W=25, Wn /W=0.2. (3 56) (3.57) (3.58) (3 59) 47 CHAPTER 4 NUMERICAL METHOD Introduction Numerical modeling of the continuous casting process primarily uses finite element and finite difference methods. Two different computational mesh variations are commonly applied to the problems; namely; time variant mesh (commonly known as moving grid technique) and fixed mesh technique. The moving grid technique is comparatively complicated and traces the solid-liquid interface with time and is limited to simpler geometries. Numerical simulations for temperature history and/or location of the solidification front result in over/under predictions, as well as numerical oscillations of the true response. Fixed grid technique, on the other hand, is considerably simpler; however, it requires accounting of the latent heat into material properties, such as specific heat or enthalpy and is restricted to simpler geometries. In the present work, the finite element method was used to solve the continuous casting problem. The fixed grid technique was used and appropriate boundary conditions were applied to account for phase change. The latent heat effect was considered in the specific heat of the material. The turbulent viscosity was determined based on the Prandtl mixing length method. 48 Finite Element Formulation In the present modeling method, the finite element procedure starts with the division of the continuum region of interest Q into a number of simple shaped regions called finite elements. Eulerian description of the elements was used; that is, the elements were assumed to be fixed in space. Sub-parametric quadrilateral nine node elements (Figure 6) were used to generate the finite element mesh. Within an element, dependent variables (u;, P and T) are interpolated in terms of values to be determined at a set of nodal points. In order to determine the equations for the unknowns at these nodal points, an individual element is separated from the global system. Node I Node4 Node Node 5 Node 2 Node3 Node 6 Node I Figure 6. Typical 9-noded iso-parametric element 49 Within each element, the approximations of the velocity, pressure and temperature fields are done as follows: where P and T are unknowns of the element nodal points and , 1P and 0 are vectors of interpolation (basis) functions and superscript T represents a vector transpose. Substitution of the above-mentioned approximations into the field equations results in the set of following functional equations: Momentum: Ul (X i J ) = ^>T (X i )U-Xt ) r(x„f) = 0r(x,)T(f) (4.1) (4.2) (4.3) /,i(0 ,'P ,0 ,% „f,T ) = ^ / , 2(@,Y,0 ,« „ f , r ) = ^ (4.4) (4.5) Incompressibility Zp(^Mi) — Rp (4.6) Energy ^ .(0 ,0 ,«„T) = ^ (4.7) where R denotes the errors (residual) resulting from the use of approximations in equations (4.1, 4.2 and 4.3). The Galerkin form of the method of weighted residual is used to reduce these errors to zero. This is achieved by making the residuals orthogonal 50 to the interpolation functions (that is, 0 , vP and 0 in equations 4.1,4.2 and 4.3) over each element expressed as follows: = 0 (4.8) = 0 (4.9) ( # , / , ) - = 0 (4.1(1) ( 0 , / r ) - ( 0 , ^ ) = O (4.11) where < a, b> denotes the inner product as follows: ( a , b } - 1 a .b d Q . a (412) Cl in equation (4.12) represents the area of the element. The detailed manipulations involving the integrals defined in equations (4.8) - (4.11) are available in Gaxtling (1987). In the present work, the penalty method was used to eliminate the terms involv ing pressure in the matrix. This reduces the overall size of the problem by considering an incompressibility condition of the following form: dU; Sx1 Zp-P (113) where sp is the penalty parameter and P is the pressure. The penalty parameter is typically a small constant of the order of IO"6, The solution of usual momentum equations with the above constraint converges to the solution of the true incompressible problem as ep-> 0 . 51 GRTt/ = - ^ .M pP (4.14) Using this approximation for the continuity equation and the Galerkin method of weighted residuals produces the following: MU+ C(U)U + KpU + Km(U3T)U + B(T)T = Fv(T) (4.15) where Kp = -G R M p1GRt (4.16) s p The equation (4.14) can be solved for pressure and substituted into the discretized momentum equations produce the following: M 0" U > _ |_ C(U)+KP+KM(U,T) B(T) T f j Fy (T) ' 0 N f 0 D(U) + Le(T)_ T tGs(T,U) (4.17) Pressure can be recovered from equation (4.14) by inverting the equation and using the known velocity as follows: P = - M p1GRtU. (4.18) In the above equations, C and D matrices represent the advection (convection) of momentum and energy respectively. Km and Le matrices represent the diffusion of momentum and energy. GR matrix is the gradient operator and GRt is the divergence operator. M and N matrices are the mass and capacitance terms in the field equations. The B matrix represents the buoyancy term while the Fy and Gs vectors provide the forcing functions for the system in terms of volume forces and surface forces. 52 Average Specific Heat Method and Prandtl Mixing Length Lee and Chiou (1995) developed the average specific heat method based on the work on the adjusted specific heat method by Hsiao (1985). Amin and Grief (1999) incorporated this method in NACHOS II to analyze the laminar behavior of the continuous casting process. To describe the computation of the adjusted specific heat at a node, let us consider two nodes at temperatures Ti and T2. Tm is the melting temperature with 2AT as the small temperature interval over which solidification takes place. Also assume that one of the two nodes is inside the solidification region and Ti > T2. In the adjusted specific heat method, Hsiao (1985) discretized the material within Ti and T2 into solid, liquid and molten phases. Adjusted specific heat, Cp (Ti, T2) at a node was calculated according to the following criteria: c,(%,%) = c „ (4.19) C W D = C,, (4.20) (4.21) 53 ? ; -A r <% < 7 ;< 7 ;+ A r c , + 0 2AT • + - Ps Pl (4.22) ^ < 7 ; - A r ^ 7 ; - A T < 7 ; < 2 ; + A T I k - % ) 2AT ( 7 ; - 7 ; -A r )+ c ^ ( 7 ; -A r -^ ) (4.23) y ? ; -A r<^<7;+AraW7;+Ar<7; ~ ^ — +C" W + Tm-T 2^ C ir (Tl -T m-A T )Cr( V 2) = I (j ; -T 2) ^2AT 2 / (4.24) I/2AT = 0 c , f f l . r 2) = + c „ , . ( T ,- r , ) + - r j ] (4 .25 ) In the solid and liquid regions, specific heat was evaluated as Cps and CpI respectively. Nodes in the two-phase zone were considered for the elements undergoing phase transition. The heat capacity at any node within the freezing interval depends on all of the four adjacent nodes, as well as the four comer nodes. The specific heat of a node undergoing change of phase can be expressed [Lee and Chiou (1995)] as follows: = + + + (4.26) where T av,i = - ( 7 Ib + T 2n + 7L + 7 V i ) (4.27) 54 where Tav, , represents the average temperature of the /th adjacent sub-element. Cp (Tlj, Tav,i), the modified specific heat, is a function of the latent heat (//), liquid specific heat (Cpi), and solid specific heat (Cps). Tjn, Tln, C3n, and T14n in equation (4.27) represent the nodal temperatures of four nodded sub-elements. Figure 7 represents a typical part of a finite element mesh. Liquid Region Node under consideration Xi+lj+1) Solidified region (4) (ij Solidification region (i+U-1) Figure 7. Typical part of the finite element mesh. 55 Once the average temperature is determined for a sub-element under consideration, T1 and T2 values are assigned so that T1 > T2 to use in equations (4.19- 4.25). Thereafter, the adjusted specific heat, Cp (Tlj, Tav, j), is calculated for that element and the average specific heat is finally evaluated, per equation (4.26). Details of this method can be found in Hsiao (1985) and Lee and Chiou (1995). The Prandtl mixing length for each node was evaluated by using the following expression [Forthmann (1934)]: Im = 0.0165(y) (4.28) where y is the distance of the node from the tip of the nozzle. Once the values of mixing length were known, turbulent viscosity was calculated by using the following expression [Gartling (1987)]: /4- ~ (4.29) h in equation (4.29) represents the shear rate (ref. equation [2.12]) and was evaluated by using the following expression: + + Sv + Sv v 4k/ (4.30) The effective viscosity / 4 was evaluated by adding the molecular viscosity to the turbulent viscosity. The effective thermal conductivity in the liquid region was determined per the following expression: K - k + - Prr (4 31) 56 PrT in equation (4.31) represents the turbulent Prandtl number. In the present work, turbulent Prandtl number was considered to be unity, per Aboutalebi et al. (1995) and White (1991). The momentum equation was coupled with the energy equation through the buoyancy term. Therefore, these two equations were solved simultaneously to calculate the velocity and temperature values. The continuity and momentum equations in the solidified region were eliminated by introducing a very high value of the viscosity essentially to solve a very stiff momentum equation. Thus, the energy equation was the only governing equation controlling the heat transfer in the solidified region. A zero value to the x-component of the velocity vector was specified upon solidification with a constant withdrawal speed in the y- direction. It should be noted here that the withdrawal speed of the cast metal is determined based on the inlet speed of the liquid metal through the nozzle and conserving the total mass. That means the mass flow rate of the liquid metal entering through the inlet nozzle should be the same as the mass flow rate of the solidified cast metal that is being withdrawn. 57 Grid Independence In order to minimize the computation time with a desired level of accuracy, the grid independence test was considered to be essential. The computational time increases with the increase in the number of elements. A coarser grid can result in steep gradients and give misleading results. The mesh was denser in the mold and nozzle region due to very high temperature and velocity gradients produced by the superheated fluid entering the nozzle. Grid independence tests with three different grid sizes (200, 400 and 800 elements) were conducted to establish the accuracy of the results. The variation of the temperature distribution for Pe=6.0, G0 =2.0, Bip =0.1, Bim =0.4 was investigated to establish the independence of the temperature distribution for a given number of elements in the computational domain. A case with higher process parameters ensured that the grid would be effective for lower values of the process parameters as well. The results of the grid independence tests are shown in Figure 8 . The maximum dimensionless temperature difference (between cases with 400 and 800 elements) along the wall and the centerline was less than 1% and 2% respectively. Therefore, the grid with 400 elements was chosen for the current study to ensure accuracy of the results with minimum computational time. 58 400 Elements 200 Elements Mold Location: 21 6 .1, the model becomes unstable, indicating a limiting value of Fe for a given superheat and cooling rates. 69 Pe=3.5 Pe=4.0 Pe=4.5 Pe=5.0 Pe=5.4 Pe=5.6 Pe=5.8 Pe=6.1 Mold Location: 21 < y* < 25 Figure 13. Effects o f Pe on solidification front, Q0=LS, Bim=0.4, Bip=O. I 70 Figure 14 represents the velocity vectors in the flow field for Pe = 3.5, 4.5, 5.6, 6.1, and ©o= 1.5, Bim= 0.4, Bip = 0.1. The vectors at the inlet (y* = 25) are uniform due to the boundary condition. Because the velocity was non-dimensionalized with inlet speed, the length of the inlet vectors is the same for different Peclet numbers. The magnitude of the vectors decreases as the superheated fluid cools down and approaches the solidification front, and eventually, it reaches a zero value at the solidification front. The results also show a back flow region near the solidification front. It can also be observed that as the Peclet number increases, the magnitude of the reverse flow decreases. This is due to the reason that as the Peclet number increase, the solidification front moves down the caster. This results in a decrease in the magnitude of the vectors in the downward direction near the solidification front. This in turn results in a weaker back circulation of the fluid from the solidification front in the caster. This is more evident from Figure 15(a), where the axial velocity (v*) profiles for the case Pe = 3.5, ®0 = 1.5, Bim= 0.4, Bip = 0.1 at various depths (y* = 24, 23, 21) are plotted across the width (x*) of the caster. Here the casting speed is represented by the horizontal uniform velocity profile at y* =21. It can be clearly observed that the axial velocity decreases with increasing lateral distance from the liquid core. For regions with reverse flow, the axial velocity becomes negative near the solidifying shell. With further increase in the lateral distance, the negative axial velocity starts to decrease, and eventually, it become positive. Upon solidification, the axial velocity becomes the same as the withdrawal speed. Also the magnitude of negative axial velocity decreases with increasing pool depth, which is more evident from Figures 15(b) and 15(c) for higher Peclet numbers. 71 Figure 14. Effects o f Pe on velocity vector field for ©0=1.5, Bim=0.4, Bip=O. I 72 a) Pe = 3.5 b) Pe =5.4 c) Pe = 6.1 Figure 15. Axial velocity profiles for 0 O= 1.5, Bim= 0.4, Bip= 0.1 73 Figure 16 represents the effect of Pe on the shell thickness for the cases with Pe=3-5, 4.0, 4.5, 5.0, 5.4, and ©o = 1.5, Bim = 0.4, Bip = 0.1. The lines in this figure represent the solidification front. Results indicate that, with increasing Fe, there is a decrease in the shell thickness. This can be attributed to the increased level of the turbulence with increasing Fe. The increased turbulence with increasing Fe results in more heat being convected down the stream, and hence, increased temperature throughout the computational domain yielding lower shell thickness. Fe = 3.5 Fe = 4.0 Fe = 4.5 Fe = 5.0 Fe = 5.4 Shell Thickness Figure 16. Effects o f Fe on shell thickness for 0 O= 1.5, Bim= 0.4, Bip= 0.1 74 Figure 17 represents the variation of solid fraction (VTV0) with Fe for the cases with Pe — 3.5, to 6.1, and ©o= 1.5, Bim= 0.4, Bip = 0.1. It is clear from the results that the solid fraction for the problem under consideration bears a good dependence on Fe, and it decreases with increasing value of Fe for the reasons mentioned earlier. A curve fit was obtained to describe a power variation of the solid fraction (VTV0) with Fe. The corresponding equation is given by: The statistical significance of the above mentioned correlation was checked by calculating the relative variance (r1) of the fitted data and the numerically obtained data as: In equation (5.2), the superscript “0” represents the fitted data (Eq. 5.1), (VNo)i numerically obtained data, and (VTVo)av represents the average of the numerically obtained results. A variance of 0.9859 for the curve fit indicates that the data is in good agreement with the curve fit. It may be noted that, the above mentioned correlation is valid only for the process parameters, and the physical domain size considered in the present work. VTVo = 1.6152(Pe)',-0.4102 (5.1) ZKvw.),- (W V .U (5.2) 75 V/Vo= 1.6152Pe"° 4,02 3 3.5 4 4.5 5 5.5 6 6.5 7 Pe Figure 17. Effects of Pe on solid fraction ©0= 1.5, Bim= 0.4, Bip= 0.1 The effects of the Peclet number on the solidification front (Figure 13) and velocity vectors (Figure14), have already illustrated that, with increasing the Peclet number, the solidification front moves down the caster, resulting in deeper penetration of the superheated metal into the liquid core. Therefore, the temperatures throughout the domain were expected to increase. Figures 18(a) and 18(b) represent the temperature distribution along the core (center of the caster) and the shell (the mold wall) in the caster respectively. Results obtained indicate a change in the average slope from mold to post mold regions due to different cooling rates. 76 Figure 19 shows the variation of local heat flux (q*") as a function of length (y*) along the outer edge of the cast metal for different Fe. From the results obtained, it can be concluded that increasing Fe results in increased heat flux throughout the domain. This is due to the reason that as the inlet speed increases, there is an increase in the down stream heat convected through the turbulent fluid. This results in increased temperature throughout the caster, and therefore, higher heat flux. These results indicate an abruptness at the mold exit (y*=21), due to the discontinuity in the Biot number (Bim^ Bip). Also, the average slopes in the mold and post mold regions are different due to different cooling rates. The variation of the average heat flux with the Peclet number (Figure 20) indicates that, with increasing the Peclet number, there is an increase in the average heat flux, regardless of the mold and post mold cooling rates. This is because with the increasing Peclet number, the turbulent fluid penetrates further down the caster, resulting in increased temperature throughout the domain. The increased outer boundary temperatures throughout the domain results in a higher temperature difference between the caster and the surrounding air. Under the assumption of the constant convective heat transfer coefficients, a large temperature difference is directly reflected in higher heat flux values. A curve fit was obtained to describe a power variation of the average heat flux (Q) with Fe. The corresponding equation is given by: Q= 13.402(Pe)1 1254 (5.3) A variance of 0.9905 for the curve fit indicates that the data is in good agreement with the curve fit 77 The effect of the Peclet number on the percentage heat extracted in the mold was another meaningful analysis performed in the present study. Figure 21 represents the variation of the percentage heat extracted in the mold with the Peclet number. The percentage heat extracted in the mold is defined as the percentage of the total heat flux (from y* = 0 to y* = 25) that is extracted in the mold (from y* = 21 to y* = 25) and is expressed as follow: (ej%= E ^ jllooy" (5'4) The results obtained indicate that the heat extracted in the mold decreases with the increasing Peclet number. This is because with increasing Peclet number, larger velocity gradients are induced in the caster. This effect, in turn, leads to higher turbulence and reduced residence time of the superheated metal in the mold; therefore, more heat is carried by the superheated metal further down the caster outside the mold. This results in decreased percentage heat removed in the mold with increasing Peclet Number. The present study indicates that the mold length has a significant effect on the percentage heat extracted in the mold. The subject is discussed in details in the mold length effect section of this chapter. 78 Pe=3.5 Pe=4.0 Pe=4.5 Pe=5.0 Pe=5.4 Pe=6 . 1 Mold Location: a) Shell Pe=3.5 Pe=4.0 Pe=^ 4.5 Pe=5.0 Pe=5.4 Pe=6 . 1 Mold Location. b) Core Figure 18. Effects of Pe on shell and core temperature distribution for 0 O= 1.5, Bim= 0.4, Bip=O-I 79 Pe=6.1 Pe=5.4 Mold Location. Pe=4.5 Pe=4.0 Pe=3.5 Figure 19. Effects of Pe on local heat flux for O0= 1.5, Bim= 0.4, Bip= 0.1 1.1254Q = 13.402(Pe) Figure 20. Effects of Pe on total heat flux for G0= 1.5, Bim= 0.4, Bip= 0.1 80 100 I .c 2 I I S x 80 60 40 20 0 3 6 Figure 21. Effects of Pe on percentage heat removed in the mold for O0= 15, Bim = 0.4, Bip = 0.1 81 Effects of Superheat To study the effects of superheat on the process of CC, the superheat (0O) was expressed as the difference between the solidification temperature and the inlet temperature. A number of simulations were performed and the results were analyzed in terms of superheat effects on various CC process characteristics. The numerical results were obtained for the following sets of input data, where ©0 was increased in steps, and all other parameters were kept constant: ©0= 1.2,1.5, 2.0, and Pe = 3.5 Bim= 0.4, Bip = 0.1 ®o= 1.2,1.5,2.0, and Pe = 6.0, Bim= 0.4, Bip = O.! The results obtained indicated that the productivity of the CC process (in terms of the location of the solidification) was greatly affected by the superheat, whereas the quality (in terms of the average slope of the solidification) was not significantly affected by the amount of the superheat. A combination of higher superheat and Peclet number resulted in breakout, indicating a limiting value of the inlet speed for a given superheat for a safe operation. Although increase in the superheat resulted in increased temperature throughout the domain, intensity of the turbulence was not significantly affected by the increased superheat. Also compared to inlet speed, superheat was found to have lesser effects on the solidification characteristics. 82 The effects of superheat (®o) on turbulent viscosity contours for cases (i) ®0 = 1.2, 1.5, 2.0, and Pe = 3.5 Bim= 0.4, Bip = 0.1 and (ii) ®0 = 1.2, 1.5, 2.0, and Pe = 6.0, Bini= 0.4, Bip = 0.1 were studied in the current research. Figure 22 represents the effects of superheat (O0) on the turbulent viscosity (pT*) for ®0=1.2, 1.5, 2.0, and Pe=3.5 Bim=OA, Bip =0.1. The turbulent viscosity was found to follow the same trend as the effect of Fe. That is, the turbulent viscosity increases with increasing ©0 and spreads over a larger area. However, it was observed that the superheat has less effect on the spread of turbulence compared to the Peclet number (Figure 12). This is because although superheat increases the heat throughout the domain, it does not have any significant effect on the velocity gradients in the superheated fluid. The turbulent viscosity is more influenced by the fluid flow.characteristics than the heat transfer and, hence, is affected to a lesser degree with increasing superheat compared to increasing the Peclet number. The results obtained here are in agreement with the results published by Hershey et al. (1993), Najjar et al. (1995) and Seyedien and Hasan (1997), where they studied the dominance of the fluid flow characteristics on the CC process. Figures 23a and 23b represent effects of superheat on the location and the slope of the solidification front for the cases with (i) ®0 = 1.2, 1.5, 2.0, and Pe = 3.5 . Bim = 0.4, Bip = 0.1 and (ii)© o= 1.2,1.5,2.0, and Pe = 6.0 Bim= 0.4 respectively. These results indicate that the position of the solidification front is greatly affected by the superheat, whereas the average slope of the front is not significantly affected by the amount of the superheat. This is because with increasing superheat, more heat needs to be extracted from the superheated fluid. Hence, for a given mold cooling rate, the front moves down 83 the caster. On the other hand, the slope of the front is a strong function of fluid velocity (as evident from Figure 13) and, hence, does not change much with superheat. From the above observations, we can conclude that, compared to the inlet speed, the superheat has a lesser effect on the quality of the end product. The results also demonstrate that, for lower Peclet number cases, the slope is less steep, which is highly desirable for a better quality of the end product. Therefore, for higher casting speeds, the amount of the superheat is not a suitable parameter for controlling the quality of the end product. It can be observed from Figure 23b that higher superheat (©0 = 2.0) in combination with a higher Peclet number (inlet speed) can result in breakout. This indicates that, for a larger amount of superheat, lower casting speeds should be used to avoid spills and other breakout related disasters. Mold length related studies (discussed in later section of this chapter) indicate that, if the process is designed to run with higher superheat and inlet speed, a pre-designed higher mold length can significantly reduce the scope of the breakout conditions. 25 24 23 22 21 20 19 18 17 16 15 0 84 1.2 b) ©o= 1.5 c) ©o= 2.0 on turbulent viscosity for Pe = 3.5 Bim 85 2 5 2 5 / / I 0 O = 1 . 2 / 0 O = 1 . 5 I © o = 2 . 0 11 2 0 2 0 - / / 0 O = 1 . 2 y / / mIl 0 / i 1 5 - 1 5 / O(NIl© / / * * > > > » 1 0 - 1 0 5 - 5 0 0 O x* I O x* I (a) (b) Mold Location: 21 . 10 - 10 - 5 - 5 0 0 0 X* I C X* I a) Pe = 3.5, G>0= 1.5, Bip = 0.10 b) Pe = 6.0, @o = 2. 0, Bip = 0.20 Figure 32 Effect OfBim on solidification front. 98 Figure 32b represents a case of breakout with a higher Peclet number (inlet speed) and superheat (inlet temperature). It is evident from the results that the higher mold cooling rate does not prevent breakout under such conditions. Therefore, a good combination of inlet speed, superheat and mold cooling rate is mandatory for a safe and productive process. Figures 33 and 34 represent mold cooling effect on the local heat flux for Bim = 0.4, 0.35, 0.30, Pe = 6.0, 0o = 2.0, Bip = O.land shell thickness for Bim= 0.4, 0.30, 0.20, Pe = 3.5, ©o = 1.5, Bip = 0.1 respectively. Evidently, with increasing the mold cooling rate, there is an increase in the local heat flux and an increase in the shell thickness. This effect is due to the higher heat extraction from the turbulent fluid with the increasing coefficient of convective heat transfer in the mold. This also results in a lower surface temperature of the cast material and is shown in Figure 35 and 36. This yields an overall higher local heat flux and a higher shell thickness. Variation of average heat flux with the mold cooling rate for Bim = 0.4, 0.30, 0.20, Pe = 3.5, ®o = 1.5, Bip = 0.1 and Bim = 0.4, 0.35, 0.30, Pe = 6.0, ©0 = 2.0, Bip = 0.1 is shown in Figures 37 and 38 respectively. The results are in agreement with the abover mentioned discussion. Figures 39 and 40 show the effect of mold cooling rate on the percentage heat extracted in the mold for Bim= 0.4, 0.35, 0.30, 0.25, 0.20, Pe = 3.5, ©o = 1.5, Bip = 0.1 and Bim = 0.4, 0.35, 0.30, Pe = 6.0, ©0 = 2.0, Bip = 0.1 respectively. It is clear from the results that irrespective of the Peclet number and the superheat, the percentage heat extracted in the mold increases with increasing mold cooling rate. It is also clear from Figure 40 that, for higher a Peclet number and higher superheat, the mold 99 cooling rate has comparatively less affect on the percentage heat removed in the mold. This is because a higher Peclet number and superheat result in higher velocity gradients, accompanied by deeper penetration of the superheated fluid into the caster. This, in turn, leads to the reduced residence time of the superheated metal in the mold, and more energy is carried by the superheated metal down the stream. Relatively cold back circulating fluid with the increased superheat of the incoming fluid in the mold reduces the effective heat transfer in the mold. Therefore, for a higher Peclet number and higher superheat, the mold cooling rate has comparatively small affect on the percentage heat removed in the mold. _ Bim = 0.35 Mold Location: * 6 — Figure 33. Effects of mold cooling rate on local heat flux for Pe=6.0, ©o=2.0, Bip =0.1 100 Shell Thickness Figure 34. Effects of mold cooling on shell thickness for Pe=3.5, 0 O=1.5, Bip=O-I 101 Mold Location: a) Shell Mold Location: b) Core Figure 35. Effects of Bim on shell and core temperature distribution for Pe = 6.0, ©o= 2.0, Bip= 0.1 102 Mold Location: a) Shell Mold Location: b) Core Figure 36. Effects of Bim on shell and core temperature distribution for Pe = 3 .5 ,0O= 1.5, Bip= 0.1 103 Figure 37. Effects of Bim on average heat flux for Pe = 3.5, 0 O= 1.5, Bip= 0.1 140 - 130 - 120 - HO - Figure 38. Effects OfBim on average heat flux for Pe = 6.0, 0 O= 2.0, Bip= 0 .1 104 X 2 0 - Figure 39. Effects O f B i m on percentage heat removed for Pe = 3.5, 0 O= 1.5, Bip= 0.1 100 0.45 Bim Figure 40. Effects o f Bim on percentage heat removed for Pe = 6.0, 0 O = 2.0, Bip= 0 .1 'I Effects of Post Mold Cooling The mode of heat transfer in the post mold region was assumed to be forced convection. The coefficient of convective heat transfer in the post mold region was expressed as Biot number (Bip). The following simulations were analyzed to investigate the effects of Bim on the CC process: Bip =0.1,0.15, 0.20, Pe = 6.0, ©0 = 2.0, Bim = 0.40 Bip= 0.1, 0.15, 0.20, Pe = 6.0, ©0 = 2.0, Bim = 0.35 Bip= 0,1, 0.15, 0.20, Pe = 6.0, ©0 = 2.0, Bim = 0.30 Bip= 0.4, 0.35, 0.30, Pe = 6.1, 0 O= 1.2, Bim = 0.25 The above mentioned conditions include the cases where the solidification front started in mold, as well as in the post mold region (breakout). This was helpful to investigate the effects of the post mold cooling rate in mold and post mold regions to validate the model for the entire length of the caster. One of the interesting findings of the analysis was that, for higher inlet speed (Pe) and lower mold cooling rate (Bim), breakout conditions can be avoided by increasing the post mold cooling rate (Bip). Although, the post mold cooling rate has a significant effect on the temperature distribution within the caster, intensity of turbulence was not significantly affected by the post mold cooling rate. Also compared to other CC parameters, the post mold cooling rate was found to have lesser effects on the quality of the cast product. 105 106 Figure 41 represents the effects of the post mold cooling rate on the non- dimensional turbulent viscosity for Bip = 0.2, 0.1 and Pe = 6.0, 0 O = 2.0, Bim = 0.4. It is clear from the results obtained that for a given inlet speed (Fe), superheat (Q0) and mold cooling rate (Bim), the post mold cooling rate does not have a significant effect on the turbulent viscosity. This is because the turbulent viscosity is mainly governed by the fluid flow characteristics not significantly affected by the post mold cooling rate. Therefore, the post mold convective heat transfer coefficient does not have a significant effect on the turbulent viscosity. Analysis of the effects of the post mold cooling rate on the location and the slope of the solidification front is represented in Figures 42a (Bip = 0.4, 0.35, 0.30, Pe = 6.1, 0 O = 1.2, Bim = 0.25), 42b (Bip = 0.1, 0.15, 0.20, Pe = 6.0, ©0 = 2.0, Bim = 0.4) and 2c (Bip= 0.1, 0.15, 0.20, Pe = 6.0, 0 O = 2.0, Bim = 0.35). It is clear from the results that, for higher inlet speed and lower superheat and mold cooling rate, increasing the post mold cooling rate can avoid a breakout condition. This is well demonstrated in the results of Figure 42a, where increasing the post mold cooling rate from Bip = 0.30 to Bip = 0.40 moved the solidification front up in the caster avoiding breakout. The case with Bip = 0.30, Pe = 6.1, ©o = 1.2, Bim = 0.25 is a case of breakout, where the solidification front is outside the mold. Figures 42b and 42c represent the effects of the post mold cooling rate on the location and the slope of the solidification front for a case with higher inlet speed, superheat and mold cooling rate. It is clear that, for such a case where the process is designed for higher inlet speed, superheat and mold cooling rate, the post mold cooling rate does not affect the location and the slope of the solidification front significantly. 25 107 Mold Location: 21 UsWZas = 6.0 (Reference equation [3.43]) where a s = IcsZpCps => UsWpCpsZks = 6.0 => Us= 6.0Z( WpCpsZks) where W = IxlO"2 m, p = 2542.5 KgZm3, Cps = 1076.0 JZKg K and ks = 238.0 WZm K. Substitution of these values in the above equation for Us gives the value of withdrawal speed as Us= 5.22x10"2 mZs. Calculations for inlet temperature: For the given set of parameters, ©0=1.5 => (T0 - Tro)Z(Ts - Too) = 1.5 (Reference Equation [3.42]) => T0= 1.5*(TS - Too) + T00 where Ts = 933.52 K, Tro = 15.77 K. Substitution of these values in the above equation for T0 gives the value of inlet temperature as T0 = 1392.4 K 156 Calculations for mold convective heat transfer coefficient For the given set of parameters, Bim= 0.4 => hm WZk = 0.4 (Reference Equation [3.44]) => hm = 0.4k/W where W = IxlO'2 m, and k = 238.0 W/m K. Substitution of these values in the above equation for hm gives the value of coefficient of convective heat transfer in the mold as hm = 9520 WZm2K. Calculations for post mold convective heat transfer coefficient For the given set of parameters, Bip= 0.1 hp WZk = 0.1 (Reference Equation [3.44]) => hp = OJkAV where W = IxlO"2 m, and k = 238.0 W/m K. Substitution of these values in the above equation for hp gives the value of coefficient of convective heat transfer in the post mold region as hp = 2380 WZm2 K. Table 2 represents values of CC parameters used by other authors. It can be observed that the range of the parameters used in the current study is within the range of the other researchers. 157 Table 2. Values of CC parameters used by other authors. Parameters Aboutalebi et al. (1995) Kang and Jaluria (1993) Inlet Temperature (0C) 1550 1119.25 Casting Speed (m/s) 0.03 0.0087 Mold Heat Transfer 1270 3570 Coefficient (W/m2 K) Post Mold Heat Transfer 1080 11,900 Coefficient(WZm2K) MONTANA STATE UNIVERSITY - BOZEMAN