Turbulent Flow Phenomena and Ce₂O₃ Behavior during a Steel Teeming Process

Abstract

Steel flow phenomena and Ce₂O₃ inclusion behavior are presented in this paper. A three-dimensional model was developed to describe the steel flow phenomena and the inclusion behavior during a teeming process. The Kim-Chen modified k-ε turbulent model was used to simulate the turbulence properties and the Height-of-Liquid model was used to capture the interface between gas and steel. A Lagrangian method was then used to track the inclusions and to compare the behaviors of different-size inclusions in the steel flow. In addition, a statistical analysis was carried out by the use of a stochastic turbulence model to investigate the behaviors of different-size inclusions at different nozzle regions. The results show that the steel flow was the most turbulent at the connection region of the straight pipe part and the expanding part of the nozzle. All inclusions with a diameter smaller than 20 μm were found to have a similar trajectory and velocity distribution in the nozzle. However, inertia force and buoyancy force were found to play an important role for the behaviors of large-size inclusions/clusters. The statistical analysis results indicate that the regions close to the connection region between different angled nozzle parts seem to be very sensitive with respect to deposition of inclusions.

1. Introduction

Inclusions in molten steel have received worldwide concern due to their serious influence on both the steel product quality and the steel production process. These inclusions may come from the deoxidation process, reoxidation by air and/or slag during steel transfer, slag at the top surface of steel, refractory of the steel containing vessel, dust of scrap, and so on. The influence of inclusions on a continuous casting process is quite serious. They can break up a casting process by clogging a tun dish nozzle when they deposit onto the wall. Overall, this nozzle clogging has been a problem for the steel continuous casting process, especially for Al-killed steel, Ca-treated steel, Rare Earth Metal treated stainless steel and Ti-killed steel. The reasons for the clogging and the possible countermeasures to prevent clogging have been extensively studied.^{1,2,3,4,5,6,7,8)} Some research has been carried out to remove inclusions in molten steel during a ladle treatment as well as during a tundish operation.^9,10,11,12) However, it is impossible to obtain a completely clean steel with the current steel production technology. Clogging is closely related to the inclusion behavior in steel. Therefore, the knowledge on the steel flow and the inclusion behavior is rather important for an understanding of the nozzle clogging process and for making a prediction on clogging situations.

Computer simulation has become an effective and inexpensive method to study the steel flow and inclusion behavior. Both of them are difficult to be measured or seen directly in production. A great amount of mathematical simulations, as is shown in Table 1,^{13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35)} have propelled our understanding of both the steel flow and inclusion behavior. As early as 1973, Szekely et al.¹³⁾ modeled fluid flow in a mold with a straight nozzle and a bifurcated nozzle, respectively. In addition, Thomas and Bai et al.^{15,16,17,18,19,20,21,22,23)} carried out systematic research on steel flows in nozzles. The effects of nozzle parameters and operating practice on the steel flow characteristics in a SEN were studied. Furthermore, steel flows in nozzles were also studied and compared by using of different turbulence models.^{30,31,32,33,34,35)} Some researchers also studied the inclusion behavior in a nozzle during casting. Wilson et al.¹⁴⁾ investigated the steel flow characteristics in a nozzle and tracked the trajectories of inclusions. The deposition of inclusions due to a centripetal force and turbulence was also studied. Yuan et al.²⁴⁾ predicted the fraction of inclusions with different densities and sizes entrapped by a lining in a stopper-rod nozzle. Zhang et al.^25,26,27) tracked the trajectories and entrapment locations of inclusions in a slide-gate nozzle. Long et al.³⁵⁾ studied the Al₂O₃ inclusion behavior in a turbulent pipe flow. Effects of factors, such as release location of inclusion, inclusion size, pipe diameter, casting speed, on the entrapment probability were investigated.

Table 1.Summary of mathematical models of fluid flow and inclusion behavior in nozzles.

Author	Year	Model	Heat	Argon	Dimension	Nozzle location	Inclusion				Verification
							Diameter, μm	Sphere	Density, kg/m3	Track method
Szekely et al.13)	1973	KPM	Y	N	2D	SEN	–	–	–	–	WM
Wilson et al.14)	1987	–	N	N	2D	GEN	5,15,25,35,50,70,90,100	Y	3500	LM	Experiment
Thomas et al.15)	1990	k-ε	N	N	2D	SEN	–	–	–	–	WM
Hershey et al.16)	1993	k-ε	N	N	2D,3D	SEN	–	–	–	–	WM
Najjar et al.17)	1995	k-ε	N	N	2D,3D	SEN	–	–	–	–	WM
Thomas et al.18)	1997	–	N	Y	2D	TN1	–	–	–	–	WM
Bai et al.19)	2000	k-ε	N	Y	3D	TN1	–	–	–	–	WM
Bai et al.20)	2001	k-ε	N	Y	3D	TN1	–	–	–	–	WM
Bai et al.21)	2001	k-ε	N	Y	3D	TN1	–	–	–	–	WM
Bai et al.22)	2001	k-ε	N	Y	3D	TN1	–	–	–	–	PM
Yuan et al.23)	2004	LES	N	N	3D	TN2	–	–	–	–	WM
Yuan et al.24)	2004	LES	N	N	3D	TN2	40,10	Y	2700,5000	LM	WM
Zhang et al.25)	2003	k-ε	N	N	3D	TN1	49,225	Y	5000	LM	PM
Zhang et al.26)	2006	k-ε	N	Y	3D	TN1	50,200,300	Y	3500	LM	PM
Zhang et al.27)	2008	k-ε	N	N	3D	TN1	5	Y	5000	LM	PM
Pfeiler et al.28)	2005	k-ε	N	Y	3D	SEN	10,500,1000	Y	3700	LM	–
Banderas et al.29)	2005	k-ε	N	Y	3D	SEN	–	–	–	–	WM
Sambasivam30)	2006	k-ε, RSM	N	N	3D	SEN	5	Y	Al2O3	LM	PT
Real et al.31)	2006	k-ε, LES	N	N	3D	SEN	–	–	–	–	WM
Real et al.32)	2007	LES	N	N	3D	SEN	–	–	–	–	WM
Zhang et al.33)	2011	k-ε, RSM	N	N	3D	TN1	–	–	–	–	WM
Chaudhary et al.34)	2011	k-ε, RANS, LES	N	N	3D	SEN	–	–	–	–	MAE
Long et al.35)	2011	–	N	N	2D	PIPE	10,20,40,60,80,100	Y	3500	LM	PM
Present study	2012	k-ε	N	N	3D	GEN	0.5,1,3,10,20,100,400	Y	6800	LM	experiment

GEN: Gradually Expanding Nozzle, KPM: Kolmogorov-Prandtl Model, LM: Lagrangian method, MAE: Metal Alloy Experiment, PM: Plant Measurement, PT: Plant Trials, SEN: Submerged Entry Nozzle, TN1: Slide-gate Controlled Tundish Nozzle, TN2: Stopper-rod Controlled Tundish Nozzle, WM: Water Model, –: didn’t study or not given, N: No, Y: Yes.

Previous modeling studies, in Table 1, mainly focused on the steel flow characteristics in a nozzle. Behaviors of different kinds of inclusions with different sizes released from different locations of a nozzle are not fully understood. In this paper, the Ce₂O₃ behavior and steel flow phenomena in a gradually expanding nozzle during a teeming process were studied. Ce₂O₃ is formed during Rare Earth Metal (REM)-alloyed stainless steel production. Ce₂O₃ is difficult to remove due to that the density is close to steel. Furthermore, it tends to stick on the nozzle wall when it comes close to the wall.^36,37,38,39) Based on experimental results, effect of Ce₂O₃ size on its behavior was studied at a teeming stage. A statistical analysis was also carried out to investigate the inclusion behavior in different nozzle regions. The aim of this paper is to investigate the Ce₂O₃ inclusion behavior in nozzle, which is not found in previous research as is indicated in Table 1. A systematic analysis combined with steel flow phenomena was expected to reveal a better understanding.

2. Experiment and Experimental Data

In the following sections, the experimental results to be used as validation for the model have been obtained from literatures and will be shown below.

2.1. Experimental Method

Steel was melted in a 600 Hz induction furnace (Fig. 1) of a 0.48 m inner diameter and with an Al₂O₃-lining, a 600 kg nominal melt size and an 800 kVA electrical output power. At the bottom of the furnace, a zirconia nozzle with a diameter of 5 mm was installed. The steel temperature was measured continuously by two Rh-Pt thermocouples. The nozzle was encased by graphite. A small induction coil outside the graphite heated the graphite. In addition, the nozzle temperature was also measured continuously by thermocouples. In this way, the temperature of ZrO₂ nozzle can be controlled to the desired level and the clogging due to the solidification of steel in nozzle can be avoided. About 350 kg steel was melted. The experimental temperature was 1480°C, while the nozzle temperature was 1620°C. After an adjustment of the steel temperature to the desired level, a de-oxidation metal with approximately 50 wt% Ce, 35 wt% La, 9 wt% Nd and 6 wt% Pr could be put into the molten steel. After inclusion separation for about 10 min, an alumina-graphite stopper-rod was elevated and the steel flowed through the nozzle and into a mold placed on a scale. The cumulative mass on the scale was sampled by a 1 Hz frequency. The steel composition is shown in Table 2.

Fig. 1.

Scheme of the experimental equipment.⁴⁰⁾

Table 2.Steel compositions.⁴¹⁾

Elements	C	Cr	Ni	Rare Earth	Al	S	Si	Mn	P
Content, wt%	0.04	20.00	11.80	0.08	0.003	0.001	1.64	0.92	0.023

2.2. Experimental Data

Firstly, the result of a teeming experiment without clogging is shown in Fig. 2. This result was used to verify the simulation of the teeming process.

Fig. 2.

Total teemed steel weight as a function of time in a teeming experiment without clogging.⁴²⁾

Bi et al.⁴¹⁾ investigated the inclusions from samples of an experiment with nozzle clogging. The typical SEM pictures of the inclusions and clusters from a film filter after electrolytic extraction of the samples are shown in Fig. 3, where the white points are REM-oxide inclusions and the grey color points are contaminations. From Fig. 3(a), it can be seen that single inclusions are very small, with a diameter less than 1 μm to a few microns. Clusters can be clearly seen in Figs. 3(b) and 3(c) with a diameter around 10 μm and 20 μm, respectively. The compositions of inclusions were also analyzed by using the EDS method, with the average compositions of about Rare Earth metal 70 wt% and O 25 wt%. REM-oxide clusters with a diameter more than 100 μm are difficult to find due to their very small number and only a few grams of steel was used in electrolytic extraction.

Fig. 3.

SEM pictures of inclusions and clusters from the electrolytic extraction of sample.⁴¹⁾

Kojola et al.⁴⁰⁾ showed the pictures of inclusion network near the nozzle wall and the nozzle center in their research. Inclusion compositions in steel and in clogged nozzle were also compared. It was found that inclusions in steel have similar compositions to the inclusions in the center of clogged nozzle, while the near clogged nozzle wall part contains somewhat higher silicon and aluminum. The picture of a part of the clogged nozzle is shown in Fig. 4(a). A similar clogged nozzle picture, as is shown in Fig. 4(b), was obtained in an experiment (Erik Roos, personal communication, January 10, 2013). It can be seen that a serious clogging region is located at the connection region of the straight pipe part and the gradually expanding part of the nozzle. The clogging of this region stops the teeming process. Although a lot of inclusions exist at the straight pipe nozzle wall, they are not the vital reason for the interruption of teeming. In Kojola et al. study,⁴⁰⁾ clogging at the upper part of the nozzle near the connection region between the nozzle and ladle bottom was also found.

Fig. 4.

Pictures of a part of the clogged nozzle that corresponds to region 5 and 6 in Figure 5(b).((a)⁴⁰⁾)

3. Mathematical Model Description

The ladle teeming process was mathematically modeled based on the same conditions as are valid for the pilot-plant experiments described earlier. The commercial software PHOENICS⁴³⁾ was used to solve the equations. A three-dimensional model of the ladle was developed. The Kim-Chen modified k-ε model⁴⁴⁾ was used to describe the turbulent flow properties. The Height-of-Liquid model, which is an option available in PHOENICS, was used to capture the moving interface between steel and air. The Log-Law wall function was used to bridge the near-wall and the outside fully-developed turbulent region. Hybrid was used as the differencing scheme. In the hybrid scheme, the diffusive flux is cut out for Peclet numbers in excess of 2.0. In addition, particle trajectory in steel was calculated by using a Lagrangian method.^45,46,47) The Lagrange approach originates from Taylor’s classic (1921) paper ‘Diffusion by Continuous Movements’.⁴⁸⁾ Taylor argued that the transport of species due to molecular diffusion is negligible compared to advective and turbulent transport and consequently the mean flow field can be evaluated from data on the motion of fluid particles.

The computational domain and 6 divided nozzle regions can be seen in Fig. 5. The radius and height of the cylindrical-polar domain is 0.24 m and 0.45 m, respectively. The time step for the calculation of the steel flow was set as 0.001 s considering the calculation stability and accuracy. The global convergence criterion was set to 0.01% for all variables.

Fig. 5.

Schematic diagram of: (a) calculation domain and particle release location and (b) nozzle regions.

3.1. Model Assumptions

The model was developed based on the following assumptions:

• Steel and air are incompressible Newtonian fluids;

• The steel/gas interface is frictionless and flat. No penetration or shear stresses or interfacial tension exist at the interface;

• Solidification and heat transfer are not considered;

• The roughness of the ladle and nozzle wall is set to zero;

• The slag at the top surface of steel is not considered;

• A fixed steel flow field is assumed during the tracking of inclusions. It is reasonable to use the fixed flow field due to that only a small change of the flow field occurs during the short time that the inclusions pass through the nozzle. In addition, the focus is to compare the inclusion behaviors in the same flow field;

• The inclusion shape is assumed to be a sphere;

• Interactions between inclusions are not considered;

• A one-way coupling between steel and inclusion was used, i.e. influence of inclusions on steel flow is not considered.

3.2. Transport Equations

Based on the model assumptions, the transport equations for the three-dimensional model of a ladle teeming process need to be solved to obtain the flow field. The particle tracking equation was thereafter solved based on the fixed flow field to obtain inclusion trajectories.

3.2.1. Governing Equation

The governing equation for the conservation of mass for a transient flow and an incompressible fluid is as follows:   

$$\frac{\partial }{\partial x_i}\left(\rho u_i\right)=0$$(1)

The conservation equation of momentum in i direction can be expressed by the following equation:   

$$\begin{array}{l}\frac{\partial }{\partial t}\rho u_i+\frac{\partial }{\partial x_j}(\rho u_i{u}_j)\\ =-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\left(\mu +{\mu }_t\right)\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\right]+\rho g_i \end{array}$$(2)

where μ is the molecular viscosity, μ_t is the turbulence viscosity and g_i is the gravity in i direction.

3.2.2. Turbulence Model

The standard high-Reynolds-number k-ε turbulence model employs a single timescale to characterize the dynamic processes occurring in turbulence. Turbulence, however, comprises fluctuating motions with a spectrum of time-scales. Therefore, a single-scale approach is unlikely to be adequate under all circumstances. In order to remedy this deficiency, Chen and Kim⁴⁴⁾ proposed a model which improves the dynamic response of the dissipation equation by introducing an additional timescale (k/P_k), where k is the turbulent kinetic energy and P_k is the volumetric production rate of k. An extra source term represents the energy transfer rate from large-scale to small-scale turbulence controlled by the production-range timescale and the dissipation-range timescale. This feature offers advantages in separated flows and also in other flows where the turbulence is removed from a local equilibrium situation. In this study, the Chen-Kim modified k-ε model has been used to calculate the turbulence viscosity:   

$${\mu }_t=\rho C_d{C}_{\mu }\frac{k^2}{\mathit{\epsilon }}$$(3)

where C_μ and C_d are constants, ε is the turbulence dissipation rate.

The equations for calculating k and ε are given below:   

$$\rho \frac{\partial }{\partial t}k+\frac{\partial }{\partial x_i}\left[\rho k{u}_i-\frac{{\mu }_t}{\text{PRT}(k)}\frac{\partial k}{\partial x_i}\right]=\rho (P_k+{\Gamma }_b-\mathit{\epsilon })$$(4)

  

$$\begin{array}{l}\rho \frac{\partial }{\partial t}\mathit{\epsilon }+\frac{\partial }{\partial x_i}\left[\rho \mathit{\epsilon }{u}_i-\frac{{\mu }_t}{\text{PRT}(\mathit{\epsilon })}\frac{\partial \mathit{\epsilon }}{\partial x_i}\right]\\ =\rho \frac{\mathit{\epsilon }}{k}(C_{1\mathit{\epsilon }}{P}_k+C_{3\mathit{\epsilon }}{\Gamma }_b+C_{2\mathit{\epsilon }}\mathit{\epsilon })-\rho C_{4\mathit{\epsilon }}\frac{P_k^2}{k}\end{array}$$(5)

where C_μ = 0.5478, C_d = 0.1643, C_1ε = 1.15, C_2ε = 1.9, C_3ε = 1.0, C_4ε = 0.25, PRT(k) = 0.75, PRT(ε) = 1.15, are constants in the turbulence model.^43,44) The parameters PRT(k) and PRT(ε) are the turbulent Prandtl number of variable k and ε, respectively.

3.2.3. Interface Tracking Model

For the Height-of-Liquid model in Phoenics,⁴³⁾ the position of the interface is determined by calculating the volume fraction of liquid, α_l, in each cell from:   

$${\alpha }_l=({{\displaystyle \sum }}^{}m-{{\displaystyle \sum }}^{}{m}_1)/({\rho }_l*V)$$(6)

The interface is defined to lie in the cells with an α_l value in the range: 0 < α_l < 1.

The density field can be obtained as follows:   

$$\rho ={\rho }_l*{\alpha }_l+{\rho }_g*(1-{\alpha }_l)$$(7)

The viscosity field is calculated as follows:   

$$\mu ={\mu }_l*{\alpha }_l+{\mu }_g*(1-{\alpha }_l)$$(8)

where ρ_l is the liquid density, ρ_g is the gas density, V is the cell volume, ∑m is the total mass of liquid in a vertical column of cells, ∑m₁ is the total mass of liquid below the interface cells.

3.2.4. Particle Tracking Model

Particle tracking is done by using a Lagrangian method.⁴⁷⁾ A particle position can be obtained by solving the following equation:   

$$\frac{d{x}_{pi}}{dt}=u_{pi}$$(9)

where x_pi is the particle position, u_pi is the particle velocity in i direction.

The particle velocity is obtained by solving the following particle momentum equation:   

$$m_{\text{P}}\frac{d{u}_{pi}}{dt}=D_{\text{P}}\left(U_i-u_{pi}\right)+m_{\text{P}}{g}_i\left(1-\frac{\rho }{{\rho }_{\text{P}}}\right)$$(10)

where m_P is the mass of the particle, D_P is a drag function, U_i is the continuous-phase velocity.

The influence of the lift force is excluded since it is questioned by many researchers regarding its significance, its region of validity and its formulation and it is not clear that it improves the results.^49,50,51) Furthermore, other hydrodynamic forces are excluded such as the Basset history, the virtual mass, the Faxen correction and the pressure gradient since they also are much smaller than the drag force for particle sizes of relevance in this study.⁵²⁾ Likewise Brownian motion only affects very small particles and is therefore not considered to be significant.⁵³⁾ The drag function, D_P, is as follows:   

$$D_{\text{P}}=\frac{1}{2}\rho \frac{\pi d_{\text{P}}^2}{4}{C}_D\left|U_i-u_{pi}\right|$$(11)

where C_D is the drag coefficient, which is given by Clift, Grace and Weber⁵⁴⁾ as follows:   

$$C_D=\frac{24}{R_{e\_p}}\left(1+0.15R_{e\_p}^{0.687}\right)+\frac{0.42}{1+4.25\times {10}^4{R}_{e\_p}^{-1.16}}$$(12)

where R_{e_p} is the particle Reynolds number. This drag coefficient function is valid for spherical particles with R_{e_p} < 3 × 10⁵.

3.2.5. Particle Stochastic Turbulence Model

In order to incorporate the effect of turbulent fluctuations on inclusion motion, a stochastic turbulent model can be used. The first approach that was based on the eddy lifetime⁵⁵⁾ spawned the eddy-interaction models in which fluid velocities (eddies) are taken to be stochastic quantities, which remains constant for the lifetime of the eddy Δt_e or, if shorter, the transit time of the particle through the eddy Δt_r.^56,57) Therefore, the continuous-phase velocity can be expressed by the following equation:   

$$U_i=u_i+{\mathit{ú}}_l$$(13)

where u_i and ${\mathit{ú}}_l$ are the continuous-phase average velocity and the fluctuating component, respectively. u_i is obtained from the Eulerian equations for the continuous-phase. Moreover, ${\mathit{ú}}_l$ is calculated by assuming that each component follows a normal distribution with a mean value of 0.0 and a standard deviation of   

$$\sigma =\sqrt{2k/\epsilon }$$(14)

The fluctuating component is assumed to act over a time interval, which is the minimum of Δt_e and Δt_r. Δt_e and Δt_r can be obtained from the following equations:   

$$\mathrm{\Delta }{t}_e=\frac{l_e}{\left|{\mathit{ú}}_l\right|}$$(15)

  

$$\mathrm{\Delta }{t}_r=\frac{l_e}{\left|U_i-u_{pi}\right|}$$(16)

where l_e is the eddy size which can be calculated by the following equation:   

$$l_e=\frac{C_{\mu }^{3/4}{k}^{3/2}}{\epsilon }$$(17)

3.3. Boundary Conditions

The detailed information of the boundary conditions of steel flow and inclusion motion is given in Table 3.

Table 3.Summary of boundary conditions.

Boundary	u	v	w	k	ε	Pgauge	D, m	Inclusion
Inlet	–	–	–	$\frac{\partial k}{\partial z}=0$	$\frac{\partial \epsilon }{\partial z}=0$	0	0.48	–
Outlet	–	–	–	$\frac{\partial k}{\partial z}=0$	$\frac{\partial \epsilon }{\partial z}=0$	0	0.005	escape
Wall	0	0	0	0	0	–	–	stick

Inlet

A constant pressure equal to the atmospheric pressure was imposed at the inlet of the ladle. Five release locations of inclusions are on the curve with a radius of 5 cm in Fig. 5(a). The exact information of these locations is shown in Table 4. Considering both the analysis effort and the aim to obtain representative result, twenty inclusions are released from each location for every case.

Table 4.Release locations of inclusions.

Locations	α,defined in Fig. 5(a)
1	5°
2	25°
3	45°
4	65°
5	85°

Ladle wall

A no-slip boundary condition was imposed for the steel and gas at the ladle wall. The inclusions were assumed to stick to the wall once they touched the ladle wall.

Outlet

A pressure-outlet condition was set at the outlet. Owing to that the nozzle outlet was exposed to air, its pressure was set to the atmospheric pressure. The inclusions escape from the calculation when they exit from the nozzle outlet.

3.4. Initial Conditions

The initial steel height from the outlet to the top steel surface is around 0.4 m for about 350 kg steel in the present ladle. Gravity is the driving force for the steel flow. The initial values of velocity and turbulence properties for the two phases are set to zero. The release velocities of particles are assumed to be zero considering a small steel flow velocity at the particle-release locations. In order to understand the inclusion behavior in molten steel as much as possible, Ce₂O₃ inclusions with different diameters, 0.5 μm, 3 μm, 10 μm, 20 μm and 100 μm, are simulated based on the sample analysis. Due to that a very small number of inclusions exist that are larger than 100 μm and to save some simulation effort, only 400 μm inclusions of those inclusions larger than 100 μm are considered in the simulation part.

3.5. Properties

According to reference data,⁵⁸⁾ the density and viscosity of the steel used in the calculation are assumed to 6998.5 kg/m³ and 6.4 × 10^–3 kg/(m·s), respectively. The air temperature is set to 20°C. The density and viscosity of air are 1.2 kg/m³ and 1.8 × 10^–5 kg/(m·s), respectively.⁴³⁾ The density of Ce₂O₃ inclusion at the steelmaking temperature is assumed to 6800 kg/m³, based on the data from reference.⁵⁹⁾

4. Results and Discussions

In order to obtain a grid independent solution, two grids, a fine grid with 16 cell layers and a coarse grid with 10 cell layers at the nozzle straight pipe part, were tested. It can be seen from Fig. 6 that the coarse grid is good enough to get a grid independent solution, with a difference of about 1.5% between the two cases. However, in order to get more accurate trajectories of inclusions, the fine grid was used in the remaining simulations. The number of elements in x, y and z direction is 20, 112 and 62, respectively. Simulation of the teeming process by using the fine grid was further verified by an experiment without clogging, which is also shown in Fig. 6. In order to save computational time, only the initial 100 s of the teeming process was calculated. A good agreement between experiment and simulation was obtained. The verified model was used to describe the steel flow phenomena. The obtained steel flow field was thereafter used to investigate the inclusion behaviors in molten steel.

Fig. 6.

Experimental results and simulation results of different grid numbers.

4.1. Steel Flow Phenomena

The steel flow field was firstly solved. The properties of the steel flow field in the nozzle at a point when around 300 kg steel is left in the ladle are shown in Figs. 7–10.

Figs. 7–10.

Predicted properties of the steel flow field in nozzle at a time when around 300 kg steel is left in ladle.

Figures 7 and 8 show the turbulent kinetic energy and turbulent dissipation rate of the steel flow field in the nozzle, respectively. It can be seen that the turbulent properties reach their maximum values approximately at the connection region of the straight pipe part and the expanding part of the nozzle (region 5). This illustrates that the steel flow is very chaotic in this region. The biggest shear stress value also exists around this region, as is shown in Fig. 9. More specifically, it illustrates that the steel flow has a large stress force on the nozzle wall. As reported by Uemura et al.,⁶⁰⁾ an inclusion that comes close to a nozzle wall sinters and forms a neck with the refractory of a nozzle wall or with an inclusion that already deposited on a nozzle wall. According to Sambasivam,³⁰⁾ the neck will be broken and the inclusion may be washed away when a shear stress at the neck is higher than the sinter bond strength. Therefore, in the regions with a high shear stress level, a deposition of an inclusion may be prevented even though the inclusion is present or moving close to the nozzle wall. Figure 10 shows the velocity distribution of the steel flow in the nozzle. The quick change of velocity contour in around the region 5 shows a high velocity gradient, which means that a turbulent flow developed very quickly. In the connection regions, where a nozzle geometry transition exists, a sudden change of the steel flow direction occurs. The velocity of radical steel flow greatly decreases due to its collision with a downwards directed steel flow at the nozzle core part. This kind of collision also increases the flow chaos, contributing to the turbulence level of the steel flow.

4.2. Inclusions Tracking Neglecting Stochastic Turbulent Motion of Inclusions

Six different sizes of inclusions, 0.5 μm, 3 μm, 10 μm, 20 μm, 100 μm and 400 μm, released from location 3 in Table 4 were tracked using a Lagrangian method under the previously obtained fixed flow field. In order to obtain a clear view on the flow abilities of different sizes of inclusions, a stochastic turbulent model for inclusion movement is not used at first. In this way, the uncertainty that a stochastic turbulent random motion leads to is avoided. Figure 11 shows the locations of inclusions at different times in the nozzle. It can be seen that inclusions with a diameter 0.5 μm, 3 μm, 10 μm, 20 μm and 100 μm have similar trajectories. However, for an inclusion with a 400 μm diameter, the trajectory is obviously different from other inclusions. It moves closer to the nozzle center and takes a much longer time before it reaches the nozzle region, as can be seen in Fig. 12(b), than the other inclusions. The behaviors of inclusions are mainly determined by three forces: i) an inertia force, ii) a drag force and iii) a buoyancy force due to a density difference between an inclusion and steel. In the current situation, the angle between the upwards buoyancy force and the downwards drag force is larger than 90°. Therefore, the drag force in the z direction needs to combat the buoyancy force to make inclusions move to the nozzle region. A larger-size inclusion has a bigger buoyancy force than a smaller-size inclusion. Meanwhile, the steel flow velocity is small at the release location 3, which leads to a smaller drag force due to a smaller magnitude of velocity difference between steel and inclusion. The relative magnitude of these two forces will determine the inclusion behavior. In order to explain the obviously different behaviors of 400 μm inclusions compared to other sizes of inclusions, the buoyancy force and the drag force for 100 μm and 400 μm inclusions at the release location 3 as well as at the straight pipe location were calculated, as is shown in Table 5. At the release location 3, it can be seen that the downwards drag force of 400 μm inclusions has a similar magnitude to the upwards directed buoyancy force. However, a much larger drag force than a buoyancy force was obtained for the 100 μm inclusions. This means that 100 μm inclusions can move much faster in the downwards z direction compared to 400 μm inclusions. The competition of these two forces in the z direction makes 400 μm inclusions to take a longer time than smaller inclusions, like 100 μm inclusions, to move to the nozzle region. This also gives them more time to move towards the nozzle center, as is shown in Fig. 11. The larger inertia of big inclusions than that of small inclusions also causes them to take a longer time to respond under the same conditions. At the straight pipe location, around 0.026 m from the nozzle outlet, the drag forces of both 100 μm and 400 μm inclusions are much larger than the buoyancy forces, which causes them to move fast in the nozzle pipe region.

Fig. 11.

Locations of inclusions at different times.

Fig. 12.

Velocities of inclusions in nozzle, (a) and (b) are y-direction and z-direction velocities of inclusions as a function of time, respectively; (c) and (d) are y-direction and z-direction velocity distributions of inclusions at different distance from the release location 3, respectively.

Table 5.Buoyancy force and drag force of inclusions in the z direction.

Location	Size, μm	Buoyancy force, N	Drag force, N	Acceleration a, m/s2
Release location 3	100	–1.02 × 10–9	1.50 × 10–8	3.93
	400	–6.53 × 10–8	6.70 × 10–8	7.66 × 10–3
Pipe location, 0.026 m from nozzle outlet	100	–1.02 × 10–9	1.50 × 10–7	36.79
	400	–6.53 × 10–8	2.09 × 10–5	91.70

The change of inclusion velocities in the y direction, parallel to the cross section of the nozzle, and the z direction, vertical to the cross section of the nozzle, as a function of time are shown in Figs. 12(a) and 12(b), respectively. The characteristics of the inclusion velocities in the x direction, which is not shown in this paper, are similar to those in the y direction, except with respect to the velocity magnitude. It can be seen that inclusions with a diameter of 0.5 μm, 3 μm, 10 μm and 20 μm have a similar velocity pattern with a minor difference of magnitude at the same elapsed time. This means that they have similar trajectories, which was previously shown in Fig. 11. With the increase of inclusion size, especially for inclusions larger than 100 μm, the maximum velocities of inclusions in both the y and z directions decrease. This can clearly be seen in Figs. 12(c) and 12(d). In the y direction, the main reason for that is the inertia of inclusions. The smaller fluctuation of the y-direction velocity of big inclusions in Fig. 12(c) illustrates the larger influence of inclusion inertia on their motion. For the z-direction inclusion velocity, both the inertia force and the buoyancy force should be responsible for a little bit smaller velocity magnitude for the larger inclusions than for the smaller inclusions. From Fig. 12(c), it can be seen that a sharp decrease of the y-direction inclusion velocities occurs at the locations of 3 cm and 7 cm away from the nozzle outlet, where the nozzle geometry transition exists. As previously mentioned, a rapid decrease of steel velocity in the y direction should be the reason for that. In Fig. 12(d), the inclusion velocity increases rapidly within a 2 cm distance, from 5 cm down to a 3 cm distance from the nozzle outlet. This location is situated just above the straight pipe part of the nozzle. This also illustrates that the turbulence intensity increases quickly in this region.

4.3. Inclusion Behavior Including Stochastic Turbulent Motions

In order to understand the inclusion behaviors at different nozzle regions, a statistical analysis was carried out to investigate the sensitivity of different nozzle regions on the possible inclusions deposition. Considering the result of 4.2 and a small number of large-size inclusions, e.g. larger than 100 μm, existing in steel, three sizes of inclusions, 1 μm, 10 μm and 100 μm, were tracked. Twenty inclusions were released from each location in Table 4.

The number of inclusions that touch the nozzle wall is shown in Fig. 13. Figure 13(a) shows the influence of inclusion release locations on the number of inclusions touching the wall. It can be seen that some inclusions released from location 1 and 2 touch the nozzle wall. However, there is no inclusion, released from location 3, 4 and 5, touching the nozzle. This illustrates that the deposition possibility of inclusions from location 3, 4 and 5 is very low. The inclusions from location 1 have the highest possibility of touching the nozzle wall among all the release locations. One possible reason is that location 1 is close to the nozzle wall, which provides an inclusion with a short distance to pass through to the nozzle wall. Figure 13(b) shows the number of different-size inclusions from all the release locations touching the nozzle wall. For each size of inclusions, 10 μm inclusions seem to have a little bit higher possibility to touch the nozzle wall than the other two sizes of inclusions. Figures 13(c) and 13(d) show the influence of inclusion sizes and release locations on the number of inclusions touching the wall of each nozzle region. It can be seen that for the 1 μm inclusions, the distribution of inclusions that touch the nozzle wall along the nozzle height is more uniform than the other two sizes of inclusions. For 10 μm and 100 μm inclusions, two nozzle regions, region 1 and region 5, have a higher number of inclusions touching the nozzle wall compared to the other regions. This illustrates that region 1 and region 5 have a higher possibility for nozzle clogging than the other regions. In region 1, the steel flow velocity is very small. The turbulence properties in Figs. 7 and 8 also show that turbulence intensity is not high in this region. Furthermore, Fig. 13(d) shows that all the inclusions that touch the nozzle wall within regions 1 to 4 are from release location 1. Therefore, the following reasons should be responsible for a large number of inclusions moving to the wall in nozzle region 1: 1) the release location 1 is close to the wall of nozzle region 1; 2) turbulence plays a positive role for the transport of inclusions to the wall; 3) centripetal force, which has also been investigated by Wilson,¹⁴⁾ should be helpful for the inclusion transport during flow direction changing in region 1. In region 5, as is previously shown in Fig. 10, steel velocity changes quickly. Therefore, the steel flow turbulence developed very fast. Furthermore, the collision between the radial steel flow and the downwards core flow contributes to the flow turbulence. Inclusions are transported to the nozzle wall by turbulent eddies. A smaller diameter of the nozzle in this region gives inclusions a shorter distance to pass through to reach the nozzle wall. Both the high turbulence intensity and the short moving distance cause inclusions easy to touch the nozzle wall. For the inclusions released from location 2, it can be seen from Fig. 13(d) that inclusions only touch the nozzle wall in region 5 and region 6. As previously mentioned, a high turbulence as well as a short moving distance should be the reason for that. However, a higher steel flow velocity in region 6 than that in region 5 causes inclusions to have a shorter residence time in the straight pipe nozzle part. This short residence time causes some inclusions not to have enough time to move to the nozzle wall in region 6, even though the flow turbulence is also high in region 6.

Fig. 13.

Number of inclusions touch the nozzle wall, (a) influence of all 5 release locations on the number of inclusions touching the nozzle wall, (b) influence of different-size inclusions from all release locations on the number of inclusions touching nozzle wall, (c) and (d) influence of inclusion sizes from all release locations and release locations for all sizes of inclusions on the number of inclusions touching the nozzle wall at each nozzle regions.

The experimental results, as is previously shown in Fig. 4, illustrate the reliability of this simulation work and the above explanation. It can be seen that a serious clogging is found in region 5, which is in a good agreement with the present simulation results. Region 1 is located at the connection part of the nozzle and ladle bottom. After experiment, the steel left in ladle must be tapped out making it difficult to get a sample of solidified steel at that region. However, supporting the model results is the comparable experiments reported by Kojola,⁴⁰⁾ demonstrating that clogging also frequently occurs in the upper part of the nozzle. Both the simulation and experiment show that the transition region of geometry, or flow field, is the sensitive region for an inclusion deposition, as well as for clogging.

The statistical analysis in this study gives some information on inclusion behaviors in nozzles. Even though it is not clear from the results of this study on whether inclusions deposit or not after touching the nozzle wall, it provides some information on where inclusions touch the nozzle wall and what places should be sensitive for nozzle clogging. Thus, a future deposition model should be developed to describe the inclusion behaviors when they come close to the nozzle wall. Also, when more complex production systems are to be studied using particle (inclusion) tracking, a stochastic model that is exactly compatible with properties of the crossing-trajectory effect as determined by Kolmogorov similarity scaling for high Reynolds-number turbulence could be used. As was pointed out by Reynolds et al.,⁶¹⁾ model predictions for particle dispersion in grid generated turbulence are shown to be in close agreement with experimental data of Snyder and Lumley⁶²⁾ and integral timescales are shown to be compatible with the much used parameterizations advocated by Csandy⁶³⁾ and by Frenkel.⁶⁴⁾ In the future, perhaps the most suitable model for the evolution of fluid velocities along particle trajectories in complex metallurgical production systems would be the Langevin equation⁶⁵⁾ used in conjunction with velocity fields calculated using the Eulerian approach.

5. Conclusions

Ce₂O₃ behavior and steel flow in a nozzle during teeming were investigated in this paper. Firstly, steel flow at the point of 300 kg steel left in ladle was studied. It was found that the turbulent properties reach their maximum values in around the connection region of the straight pipe part and the expanding part of the nozzle. At the connection part, the radial directed flow velocity rapidly decreases due to flow collision with the downwards directed steel flow, which increases the flow turbulence. Given the steel flow, single inclusions with different sizes were tracked. It was found that 0.5 μm, 3 μm, 10 μm and 20 μm inclusions have similar trajectories and velocity distributions in the nozzle. However, the trajectories of larger inclusions (400 μm) were quite different from the smaller ones. Both the inertia force and the buoyancy force play a very important role for the behaviors of large-size inclusions. The results of the statistical analysis performed for 1 μm, 10 μm and 100 μm inclusions show that inclusions that enter the nozzle inlet from a close-wall location have a high probability of touching the nozzle wall. The nozzle inlet region and the connection region of the straight pipe and the expanding part of the nozzle were found to be the sensitive regions for an inclusion deposition as well as for nozzle clogging. A future deposition model is required to describe the inclusion behaviors when they come close to the nozzle wall.

Acknowledgements

Peiyuan Ni wants to thank Chinese Scholarship Council for the financial support of his study at KTH-Royal Institute of Technology, Stockholm, Sweden. He would like to thank Anders Tilliander for the help regarding the software and computer. Erik Roos and Yanyan Bi are also appreciated for their help regarding experiment part of this paper.

Nomenclature

ρ density of fluid, kg/m³

ρ_l density of liquid phase, kg/m³

ρ_g gas density, kg/m³

ρ_P density of particle, kg/m³

u_i, u_j mean fluid velocity over a time step in particular directions, m/s

U_i instantaneous velocity of fluid phase, m/s

${\mathit{ú}}_l$ fluctuating component of fluid phase velocity, m/s

u_pi particle velocity in i direction, m/s

x_i, x_j particular directions

x_pi particle position

t time, s

Δt_e lifetime of the local eddy, s

Δt_r transient time taken for the particle to cross the eddy, s

μ molecular viscosity of fluid phase, Pa·s

μ_l molecular viscosity of liquid phase, Pa·s

μ_g molecular viscosity of gas phase, Pa·s

μ_t turbulence viscosity of fluid phase, Pa·s

∑m total mass of liquid in a vertical column of cells, kg

∑m₁ total mass of liquid below the interface cells, kg

m_P mass of particle, kg

d_P diameter of particle, m

g_i gravity in i direction, m/s²

p static pressure, N/m²

k turbulent kinetic energy, m²/s²

ε turbulent kinetic energy dissipation rate, m²/s³

α_l volume fraction of liquid phase

V cell volume, m³

l_e eddy size, m

C_μ , C_d, C_1ε, C_2ε, C_3ε, C_4ε, PRT(k), PRT(ε) constants in turbulence models

References

• 1)   L.  Jr. Trueba, K. D. Peaslee, J. D. Smith and M. Karakus: Steel Res. Int., 77 (2006), 37.
• 2)   J.  Szekely and  R. T.  Yadoya: Metall. Mater. Trans., 4 (1973),1379.
• 3)   S.  Ramachandran,  K. D.  Peaslee and  J. D.  Smith: ISS Trans., 30 (2003), 55.
• 4)   C.  Toulouse,  A.  Pack,  A.  Ender,  S.  Petry: Steel Res. Int., 79 (2008), 149.
• 5)   M.  Tsuda,  H.  Shinagawa,  R.  Kamata,  Y.  Hiraga and  T.  Hara: J. Tech. Assoc. Refract., 20 (2000), 138.
• 6)   W.  Fix,  H.  Jacobi and  K.  Wunnenberg: Steel Res., 64 (1993), 71.
• 7)   Y.  Fukuda,  Y.  Ueshima and  S.  Mizoguchi: ISIJ Int., 32 (1992), 164.
• 8)   R.  Dekkers,  B.  Blanpain and  P.  Wollants: Proc. of 60th Electric Furnace Conf., ISS, Warrendale, PA, (2002), 287.
• 9)   Y. J.  Kwon,  J.  Zhang and H. G. LEE: ISIJ Int., 48 (2008), 891.
• 10)   L.  Wang,  Q.  Zhang,  S.  Peng and  Z.  Li: ISIJ Int., 45 (2005), 331.
• 11)   H.  Arai,  K.  Matsumoto,  S. I.  Shimasaki and  S.  Taniguchi: ISIJ Int., 49 (2009), 965.
• 12)   G. S.  Diaz,  A. R.  Banderas,  J. j.  Barreto and  R. D.  Morales: Steel Res. Int., 80 (2009), 223.
• 13)   J.  Szekely and  R. T.  Yadoya: Metall. Trans., 4 (1973), 1379.
• 14)   F. G.  Wilson,  M. J.  Heesom,  A.  Nicholson and  A. W. D.  Hills Ironmaking Steelmaking, 14 (1987), 296.
• 15)   B. G.  Thomas,  L. J.  Mika and  F. M.  Najjar: Metall. Trans. B, 21B (1990), 387.
• 16)   D. E.  Hershey,  B. G.  Thomas and  F. M.  Najjar: Int. J. Numer. Methods Fluid., 17 (1993), 23.
• 17)   F. M.  Najjar,  B. G.  Thomas and  D. E.  Hershey: Metall. Mater. Trans. B, 26B (1995), 749.
• 18)   B. G.  Thomas,  A.  Dennisov and  H.  Bai: Proc. of ISS 80th Steelmaking Conf., Warrendale, PA, (1997), 375.
• 19)   H.  Bai and  B. G.  Thomas: Materials Processing in the Computer Age III, TMS Annual Meeting, eds. by V. Voller and H. Henein, TMS, Warrendale, PA, (2000), 85.
• 20)   H.  Bai and  B. G.  Thomas: Metall. Mater. Trans. B, 32B (2001), 253.
• 21)   H.  Bai and  B. G.  Thomas: Metall. Mater. Trans. B, 32B (2001), 269.
• 22)   H.  Bai and  B. G.  Thomas: Metall. Mater. Trans. B, 32B (2001), 707.
• 23)   Q.  Yuan,  B. G.  Thomas and  S. P.  Vanka: Metall. Mater. Trans. B, 35 (2004), 685.
• 24)   Q.  Yuan,  B. G.  Thomas and  S. P.  Vanka: Metall. Mater. Trans. B, 35 (2004), 703.
• 25)   L.  Zhang and  B. G.  Thomas: Proc. of XXIV National Steelmaking Symp., Morelia, Mich, Mexico, (2003), 184.
• 26)   L.  Zhang: 5th Int. Conf. on CFD in the Process Industries, CSIRO, Melbourne, (2006), 1.
• 27)   L.  Zhang,  Y.  Wang and  X.  Zuo: Metall. Mater. Trans. B, 39B (2008), 534.
• 28)   C.  Pfeiler,  M.  Wu and  A.  Ludwig: Mater. Sci. Eng. A, 413-414 (2005), 115.
• 29)   A. R.  Banderas,  R. D.  Morales,  R. S.  Perez,  L. G.  Demedices and  G. S.  Diaz: Int. J. Multiphas. Flow, 31 (2005), 643.
• 30)   R.  Sambasivam: Ironmaking Steelmaking, 33 (2006), 439.
• 31)   C.  Real,  R.  Miranda,  C.  Vilchis,  M.  Barron,  L.  Hoyos and  J.  Gonzalez: ISIJ Int., 46 (2006), 1183.
• 32)   C.  Real,  L.  Hoyos,  F.  Cervantes,  R.  Miranda,  M. P.  Pardave and  J.  Gonzalez: Mecanica Computational, XXVI (2007), 1292.
• 33)   X.  Zhang,  X.  Jin,  Y.  Wang,  K.  Deng and  Z.  Ren: ISIJ Int., 21 (2011), 581.
• 34)   R.  Chaudhary,  C.  Ji,  B. G.  Thomas and  S. P.  Vanka: Metall. Mater. Trans. B, 42 (2011), 987.
• 35)   M.  Long,  X.  Zuo,  L.  Zhang and  D.  Chen: ISIJ Int., 50 (2010), 712.
• 36)   A.  Memarpour,  V.  Brabie and  P. G.  Jönsson: Ironmaking Steelmaking, 38 (2011), 229.
• 37)   A.  Katsumata and  H.  Todoroki: ISS Trans., 29 (2002), 51.
• 38)   M.  Andersson,  J.  Appelberg,  A.  Tilliander,  K.  Nakajima,  H.  Shibata,  S.  Kitamura,  L.  Jonsson and  P.  Jönsson: ISIJ Int., 46 (2006), 814.
• 39)   J. W.  Farrell and  D. C.  Hilty: AIME Proc. Electric Furnace Conf., AIME, New York, (1971), 31.
• 40)   N.  Kojola,  S.  Ekerot,  M.  Andersson and  P. G.  Jönsson: Ironmaking Steelmaking, 38 (2011), 1.
• 41)   Y.  Bi,  E.  Roos,  A.  Karasev and  P.  Jönsson: Proc. of 8th Int. Conf. on Clean Steel, OMBKE, Budapest, Hungary, (2012).
• 42)   E.  Roos,  A.  Tilliander,  A.  Karasev,  P.  Jönsson and  S.  Ekerot: Pilot Scale Trials about Clogging Phenomena in Ce-alloyed Stainless Steel, ISRN, Stockholm, KTH/MSE-13/03-SE+APRMETU/TR, (2013).
• 43)  PHOENICS, PHOENICS-VR Reference Guide, in CHAM TR 326, CHAM, London, (2008).
• 44)   Y. S.  Chen and  S. W.  Kim: NASA Report, NASA CR-179204, George C. Marshall Space Flight Center, Alabama, (1987).
• 45)   C. T.  Crowe,  M. P.  Sharma and  D. E.  Stock: ASME J. Fluids Eng., (1977), 325.
• 46)   J.  Moffat and  K.  Pericleous: PHOENICS Journal of Computational Fluid Dynamics and Its Applications, 2 (1989), 21.
• 47)   J. C.  Ludwig,  N.  Fueyo and  M. R.  Malin: The GENTRA User Guide, in CHAM TR 211, CHAM, London, (2006).
• 48)   G. I.  Taylor: Proc. London Math. Soc., 20 (1921), 196.
• 49)   T.  Hussein,  J.  Smolik,  V. -M.  Kerminen and  M.  Kulmala: Aerosol Sci. Technol., 46 (2012), 44.
• 50)   M. R.  Sippola and  W. W.  Nazaroff: Lawrence Berkeley National Laboratory Report LBNL – 51432, Lawrence Berkeley National Laboratory, Berkeley, California, (2002).
• 51)   M. J.  Zufall,  W. P.  Dai,  C. I.  Davidson and  V.  Etyemezian: Atmos. Environ., 33 (1999), 4273.
• 52)   H.  Ounis and  G.  Ahmadi: J. Fluids Eng., 112 (1990), 114.
• 53)   L.  Amy and  A.  Goodarz: Aerosol Sci. Technol., 16 (1992), 209.
• 54)   R.  Clift,  J. R.  Grace and  M. E.  Weber: Drops and Particles, Academic Press, New York, (1978).
• 55)   A. D.  Gosman,  E.  Ioannides: AIAA 19th Aerospace Sciences Meeting, Paper 81-0323, Reston, VA, USA, (1981).
• 56)   A.  Ormancey: PhD. thesis, Ecole Nationale Supérieure des Mines, Paris, (1984).
• 57)   D. I.  Graham: Int. J. Multiphas. Flow, 24 (1998), 335.
• 58)   Y.  Kawai and  Y.  Shiraishi: Handbook of Physico-chemical Properties at High Temperatures, ISIJ, Tokyo, (1988).
• 59)   G. V.  Samsonov,  C. N.  Turton and  T. I.  Turton: The Oxide Handbook, IFI/Plenum, New York, (1973).
• 60)   K. I.  Uemura,  M.  Takahashi,  S.  Koyama and  M.  Nitta: ISIJ Int., 32 (1992), 150.
• 61)   A. M.  Reynolds and  G. L.  Iacono: Phys. Fluids, 16 (2004), 43
• 62)   W. H.  Snyder and  J. L.  Lumley: J. Fluid Mech., 48 (1971), 41.
• 63)   G. T.  Csandy: J. Atmos. Sci., 20 (1963), 201.
• 64)   F. N.  Frenkel: Adv. Appl. Mech., 3 (1953),61.
• 65)   S. B.  Pope: Turbulent Flows, Cambridge University Press, Cambridge, UK, (2000), 483.