Temperature Corrected Turbulence Model for Supersonic Oxygen Jet at High Ambient Temperature

Abstract

It is known that the two-equation turbulence models under-predict the turbulence mixing shear process for the compressible non-isothermal jet. In this study, a temperature corrected turbulence model on the basis of the SST k-ω model was proposed to predict the behavior of the cold supersonic oxygen jet injecting into the high-temperature environment. The results show that the corrected SST k-ω model is superior to previous corrected k-ε models when simulating the whole process of the supersonic oxygen jet flowing from a Laval nozzle into free space. Meanwhile, the calculation results of the outlet free jet at 285 K, 772 K and 1002 K are in good agreement with experimental data and empirical formulas in the literature. Furthermore, the behavior of the supersonic oxygen jet at 1873 K is predicted by the proposed model. In addition, the effects of the ambient temperature on the jet core length and the interaction between multiple jets were also studied.

1. Introduction

Supersonic oxygen jets are widely used for oxidizing the impurity elements and stirring the bath of liquid metal in steelmaking and other metal refining processes. The jet speed attenuates from supersonic to subsonic with the jet mixing with the surrounding gas. The study of the supersonic oxygen jet behavior and its attenuation law are essential to obtain the optimal kinetic parameters in blowing process, for example, in the converter steelmaking process.

The characteristic of the supersonic jet flow has been investigated by many researchers through cold model experiments.^1,2,3,4) However, there were few studies on the jet behavior under high ambient temperature condition. SUMI et al.⁵⁾ researched the supersonic oxygen jet behavior at 285 K, 772 K and 1002 K by hot model experiments. The results suggested that the velocity attenuation of the jet was restrained and the potential core length extended in high temperature environment. Moreover, some researchers^6,7,8) tried to predict the jets flow at high temperature using the existing turbulence models while the results need to be verified.

In numerical calculation process, the effect of compressible dissipation when computing high-speed jet has been taken into account. However, temperature fluctuations also affect the turbulence transport, which is not taken into consideration. This leads to poor prediction when simulating non-isothermal jet flow using the exiting two-equation turbulence models.^9,10,11) It is because these models lack the ability to predict the observed decrease in growth rate of mixing layer. Therefore, some researchers^12,13,14) have attempted to correct the k-ε turbulent model and the results are shown below:

Turbulent viscosity in standard k-ε turbulent model is defined as   

$${\mu }_t=C_{\mu }\frac{\rho k^2}{\epsilon },\mathrm{\hspace{0.17em} }{C}_{\mu }=0.09$$(1)

where ρ is the gas density, k is the turbulent kinetic energy, ε is the turbulent dissipation rate.

Sarkar¹²⁾ modified C_μ in compressive shear flow as a function of Mach number gradient defined as:   

$$M_g=\frac{Sl}{a}$$(2)

  

$$C_{\mu }=0.07exp\left(-0.4M_g\right)$$(3)

where S is the mean shear rate, l is the turbulence length scale, this modified model accurately predicts the behavior of the axial compressible jet at room temperature, but under-predict the jet behavior at high temperature.

For the high temperature jet flow into low-temperature environment, Hamid et al.¹³⁾ proposed a corrected k-ε turbulence model with increasing the turbulence viscosity. In this model, the temperature gradient was introduced to correct the turbulent viscosity coefficient C_μ as a function of local total temperature gradient normalized by the turbulence length scale T_g:   

$$T_g=\sqrt{{\left(\frac{\partial T_t}{\partial x_i}\right)}^2}\frac{\left(k^{3/2}/\epsilon \right)}{T_t}$$(4)

The turbulence March number defined as,   

$$M_{\tau }=\frac{\sqrt{2k}}{a}$$(5)

where a is the local speed of sound.

The corrected turbulent viscosity coefficient C_μ:   

$$C_{\mu }=0.09C_T$$(6)

the correction factor C_T:   

$$C_T=\left[1+\frac{C_1{T}_g^m}{1+C_{2}f\left(M_{\tau }\right)}\right]$$(7)

Alam et al.¹⁴⁾ predicted the behavior of the cold oxygen free jet ejected from the Laval nozzle flowing in high-temperature environment by a corrected k-ε model. The correction variable is same as Abdol-Hamid et al.¹³⁾ and the corrected C_μ is as follows:   

$$C_{\mu }=\frac{0.09}{C_T}$$(8)

C₁, C₂ and m are the constants in C_T, they were determined by matching experimental data.

The modification of C_μ by Abdol-Hamid et al.¹³⁾ and Alam et al.¹⁴⁾ can be expressed as:   

$$C_{\mu }=0.09C_T^n$$(9)

where   

$$n=\{\begin{array}{l}+1\\ 0\\ -1\end{array}\hspace{1em}\mathrm{\hspace{0.17em} }\begin{array}{c}{T}_{\text{jet}}>T_{\text{ambient}}\\ T_{\text{jet}}=T_{\text{ambient}}\\ T_{\text{jet}}<T_{\text{ambient}}\end{array}$$(10)

It is clear that the above corrected turbulence models^12,13,14) are all based on the k-ε model.¹⁵⁾ At the same time, it is known that the k-ε model is only suitable to simulate the fully development turbulence flow with high Reynolds number and lack the ability to predict the boundary layer flow near the wall. However, the SST k-ω model introduced by Menter¹⁶⁾ fills in this gap. It combines the k-ω model and k-ε model so that the k-ω is used in the inner region of the boundary layer and switches to the k-ε in the free shear flow. Recently, Zhao¹⁷⁾ has found that the calculation results of SST k-ω model have a better agreement with the cold experimental data¹⁸⁾ than other two-equation models when simulating the supersonic oxygen jet at room temperature. However, the effects of the large temperature gradient on non-isothermal jet have not yet been taken into account in the existing SST k-ω model.

In this study, the SST k-ω turbulent model was temperature corrected and used to predict the cold supersonic oxygen jet from the Laval nozzle injecting into the high-temperature environment. The simulation results were compared with the experimental results of SUMI et al.⁵⁾ and some empirical data under same conditions. In addition, the supersonic oxygen jet behavior at 1873 K was predicted.

2. Computational Model

2.1. Numerical Simulation Assumptions

In order to satisfy the transport equations, flowing assumptions were made:

(1) The gas, oxygen, is compressible and complies with the ideal gas law.

(2) The molecular viscosity of gas is a function of temperature and complies with the Suthurland’s formula.

(3) The flow, in Laval nozzle, is isentropic and non-isothermal in fluid domain.

2.2. Governing Equations

The mass, momentum and energy transport equations based on the above assumptions are satisfied the compressible Navier-Stokes equations expressed as:   

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

  

$$\frac{\partial \left(\rho u_j{u}_i\right)}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\frac{\partial {\sigma }_{ij}}{\partial x_j}$$(12)

  

$$\frac{\partial \left(\rho u_j{c}_{p}T\right)}{\partial x_j}=u_j\frac{\partial p}{\partial x_j}+{\sigma }_{ij}\frac{\partial u_i}{\partial x_j}-\frac{\partial q_j}{\partial x_j}$$(13)

where μ is the molecular dynamic viscosity; q_i is the heat conduction vector; σ_ij is the viscous stress tensor; c_p is the specific heat capacity.

2.3. Turbulent Model and Temperature Correction

In SST k-ω model, the transport equations of turbulence kinetic energy, k, and the specific dissipation rate, ω, are as follows:   

$$\frac{\partial \left(\rho k\right)}{\partial t}+\frac{\partial \left(\rho u_{j}k\right)}{\partial x_j}={\tau }_{ij}\frac{\partial u_i}{\partial x_j}-{\beta }^{*}\rho \omega k+\frac{\partial }{\partial x_j}\left[\left(\mu +{\sigma }_k{\mu }_t\right)\frac{\partial k}{\partial x_j}\right]$$(14)

  

$$\begin{array}{c}\frac{\partial \left(\rho \omega \right)}{\partial t}+\frac{\partial \left(\rho u_j\omega \right)}{\partial x_j}=\frac{\gamma }{{\nu }_t}{\tau }_{ij}\frac{\partial u_i}{\partial x_j}-\beta \rho {\omega }^2+\frac{\partial }{\partial x_j}\left[\left(\mu +{\sigma }_{\omega }{\mu }_t\right)\frac{\partial k}{\partial x_j}\right]\\ +2\left(1-F_1\right)\frac{\rho {\sigma }_{\omega 2}}{\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}\end{array}$$(15)

  

$$F_1=tanh\left({\mathrm{\Phi }}_1^4\right)$$(16)

  

$${\mathrm{\Phi }}_1=min\left[max\left(\frac{\sqrt{k}}{0.09\omega y},\frac{500\nu }{\rho y^2\omega }\right),\frac{4\rho k}{{\sigma }_{\omega ,2}{D}_{\omega }^{+}{y}^2}\right]$$(17)

  

$$D_{\omega }^{+}=max\left[2\rho \frac{1}{{\sigma }_{\omega ,2}}\frac{1}{\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j},{10}^{-20}\right]$$(18)

The turbulent viscosity is defined as:   

$${\mu }_t=\frac{\rho a_{1}k}{max\left(a_1\omega ,\mathrm{\hspace{0.17em} }\mathrm{\Omega }{F}_2\right)}$$(19)

  

$$F_2=tanh\left({arg}_2^2\right)$$(20)

  

$${arg}_2=max\left(2\frac{\sqrt{k}}{0.09\omega y};\mathrm{\hspace{0.17em} }\frac{500\nu }{y^2\omega }\right)$$(21)

where τ_ij is the Reynolds stress tensor, Ω is the vorticity tensor, ω is the specific dissipation rate, ν is the laminar flow viscosity, y is the distance to the wall, the model constants a₁=0.31, β^*=0.09, σ_k1=0.85, σ_ω1=0.5, σ_k2=1.0, σ_ω2=0.856.

The effects of the temperature gradient on turbulent mixing were taken into account in the SST k-ω model. A temperature corrected turbulence viscosity was proposed with local total temperature gradient as a variable by the local turbulence length scale T_g:   

$$T_g=\frac{\sqrt{{\left(\frac{\partial T_t}{\partial x_i}\right)}^2}\cdot \frac{k^{1/2}}{\omega }}{T_t}$$(22)

The function expression of T_g is similar to Eq. (4) used by Abdol-Hamid et al.,¹³⁾ but the formula of the turbulence length scale in SST k-ω model converts from k^3/2/ε to k^1/2/ω.

The variable C_T was determined by the similar function of Eq. (7), which were used by Abdol-Hamid et al.¹³⁾ and Alam et al.¹⁴⁾ The coefficients and constants of C_T were obtained by trial and error to match the hot experimental data of SUMI et al.⁵⁾ on the supersonic jet at high temperature. The coefficients and constants in Eq. (7) are shown in Table 1.

Table 1. Coefficients and constants of C_T.

Authers	Turbulence Model	m	C1	C2
Abdol-Hamid et al.13)	k-ε	3	24.39	24.39
Alam et al.14)	k-ε	0.6	1.2	1
Present study	SST k-ω	0.5	1.2	1

The C_T in the corrected SST k-ω model is expressed as:   

$$C_T=\left[1+\frac{T_g^{0.5}}{1+1.2f\left(M_{\tau }\right)}\right]$$(23)

where f(M_τ)=($M_{\tau }^2$−$M_{\tau 0}^2$)H(M_τ−M_τ0), H(x) is the Heaviside step function; M_τ0=0.1 and f(M_T)=0 for no compressibility correction.

The corrected turbulence viscosity in SST k-ω model is expressed as:   

$${\mu }_{t\mathrm{\hspace{0.17em} }\text{corrected}}=\frac{{\mu }_{t\mathrm{\hspace{0.17em} }\text{SST}k-\omega }}{C_T}$$(24)

The temperature modification occurs in the turbulent shear layer where the local total temperature changes. However, flow characteristics not related to shear layer mixing, such as flow expansion, compression, internal shocks and near-wall boundary layer flow will not be affected by this temperature corrected turbulence model.

Compared with the previous modified k-ε models,^13,14) the superiority of the corrected SST k-ω model proposed in present study is that it uses the k-ω model which can accurately predict the low Reynolds flow in the near-wall boundary layer and switches to the temperature corrected k-ε model for the free shear flow. Therefore, the corrected SST k-ω model is suitable to predict the jet flowing from the Laval nozzle into high temperature environments.

2.4. Computational Domain and Simulation Methodology

The computational domain of the single jet was divided into 41280 grids by the software GAMBIT ver.2.4 as shown in Fig. 1. Half of the computational domain and the axisymmetric boundary condition were selected to save the computation time. Geometry and boundary conditions are shown in Table 2, which is consistent with the experimental conditions of SUMI et al.⁵⁾

Fig. 1.

Grid configuration for the simulation of single jet. (Online version in color.)

Table 2. Geometric parameters and boundary conditions.

Parameter	Value for single jet	Value for multiple jets
Inlet stagnation pressure (Pa)	515491	768545
Outlet pressure (Pa)	101325	101325
Inlet oxygen temperature (K)	285	285
Nozzle diameter (mm)	Throat: 7.9; Exit: 9.2	Throat: 31.5; Exit: 40.5
Ambient temperature (K)	285, 772, 1002, 1873	1873
Inclination (°)	–	11.5
Mach number at the exit	1.72	1.98

The calculation process was performed by the software ANSYS FLUENT 16.0. For compressible supersonic oxygen jets, density-based solver was selected. The corrected turbulent viscosity, written in the C programming languages, was applied to the solving process in the form of User Define Function (UDF).

3. Results and Discussion

3.1. Velocity and Pressure Distribution of the Single Jet

Figure 2 shows the axial velocity distribution calculated by two-equation turbulence models at 285 K, 772 K and 1002 K. The reference experimental data and simulation results of all models are close at room temperature as shown in Fig. 2(a) and then the gaps appear as the ambient temperature increase as shown in Figs. 2(b) and 2(c). However, the calculation results of corrected SST k-ω model still have a good agreement with the experimental ones at 772 K and 1002 K. Besides, in Figs. 2(b) and 2(c), the standard k-ω model seems much better than SST k-ω model. The reason is as follows: when calculating the free shear flow, the SST k-ω model switches to the standard k-ε model which is known to under-predict the cold turbulence jet at high temperature. While, the standard k-ω model uses the specific dissipation rate(ω) rather than the dissipation rate(ε) is more sensitive to the change of turbulent kinetic energy in the turbulence mixing process than SST k-ω model.

Fig. 2.

The calculated axial velocity distribution with different turbulence models at different ambient temperatures; (a) 285 K; (b) 772 K; (c) 1002 K.

Figure 3 shows that the calculated axial velocity distribution with the corrected SST k-ω model is in a good agreement with the experimental values in the literature⁵⁾ at 1002 K, while the previous corrected k-ε models^13,14) under-predict it. It indicates that the corrected SST k-ω model is superior to the corrected k-ε models^13,14) when simulating the whole flow process of the oxygen jet which contains flow expansion, compression, internal shocks, near-wall boundary layer flow and free shear flow. Because the proposed model is based on the SST k-ω model which combines the k-ω model and k-ε model and has high accuracy and reliability in a wide range of flow fields.

Fig. 3.

The calculated axial velocity distribution with different corrected models at 1002 K.

The calculated axial velocity and dynamic pressure distribution of the oxygen jet with the proposed model at different ambient temperature are consistent with the experimental values as shown in Figs. 4 and 5. In addition, the effects of ambient temperature on the transmission of dynamic pressure weakens as the distance from the outlet increases. Given this good matching degree at 285 K, 772 K and 1002 K, the supersonic oxygen jet behavior, at 1873 K, was predicted. The calculation results suggest that the axial velocity and dynamic pressure decay at 1873 K will be more slow than those at 285 K, 772 K and 1002 K, meanwhile, the potential core length extends.

Fig. 4.

The calculated axial velocity distribution with the proposed model at different temperature.

Fig. 5.

The calculated axial dynamic pressure distribution with the proposed model at different temperature.

Ito and Muchi¹⁹⁾ proposed the empirical formula for supersonic oxygen jet, as shown in Eq. (25):   

$$-\frac{1}{2ln\left(1-U_{\text{m}}\right)}=\alpha \cdot {\left(\frac{{\rho }_a}{{\rho }_{\text{e}}}\right)}^{1/2}\cdot \frac{x}{d_{\text{e}}}-\beta ,\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }{U}_{\text{m}}=\frac{U}{U_{\text{e}}}$$(25)

where U is the axial velocity, U_e is the exit velocity, ρ_a is the ambient gas density, ρ_e is the jet density at the nozzle exit, d_e is the nozzle exit diameter, α=0.0841 and β=0.6305.

The ratio U/U_e is a function of (ρ_a/ρ_e)^0.5(x/d_e) as shown in Fig. 6. Both the computational and experimental results have a good agreement with the empirical formula at the four ambient temperatures.

Fig. 6.

The dimensionless velocity ratio as a function of the (ρ_a/ρ_e)^1/2(x/d_e) at different temperature.

Figure 7 shows the radial distribution of the velocity at different distances from nozzle exit. At the position of x/d_e=5, the jet has high speed with the potential core inside and the turbulence mixing layer at the external of the potential core is thin, so the radial distribution of the velocity is little affected by ambient temperature. At the position of x/d_e=20, the jet potential core only exists at 1873 K, and then the radial velocity distributions become different as the temperature increases. Compared with the exit velocity of the nozzle, the radial velocity at the position of x/d_e=50 significantly decreases at the four temperatures as the potential core disappears and the turbulent mixing fully develops.

Fig. 7.

The radial distribution of the velocity at different distances from the nozzle exit; (a) x/d_e=5; (b) x/d_e=20; (c) x/d_e=50.

3.2. Temperature Distribution of the Single Jet

In this corrected model, the Prandtl number was 0.715 obtained by Kleinstein²⁰⁾ and SUMI et al.⁵⁾ has further verified that this value is still suitable when the March number is 1.72.

Figure 8 shows the static temperature distribution on the center axis at different ambient temperatures. The calculation results are consistent with the experimental values at 285 K. However, the experimental measurement value becomes greater than the calculated one and the gap between the two gradually increases as the ambient temperature rises. The reasons may be as follows: (1) The Prandtl number varies with the temperature, while it is a fixed value in this calculation process, which may bring a certain influence on the accuracy of the calculation. (2) The thermocouples may receive an extra radiant heat from the ambient gas or the furnace wall at axial measuring points for the axial temperature of the jet is lower than the ambient temperature. Besides, the axial static temperature of the jet flow in the position of x/d_e=21.7 is the lowest among the four measurement points, and theoretically the radiant heat absorbed by the thermocouple in this measurement point should be the most, which may explain why the gap between the calculation results and the experimental measurements of this position at high-temperature is the largest in Fig. 8. For the above reasons, the static temperature distribution of the non-isothermal jet still needs further study.

Fig. 8.

The calculated and experimental axial static temperature at different temperature.

3.3. The Potential Core Length

The empirical formula of potential core length can be obtained when the U_m in Eq. (25) equals 1, and then the potential core length is expressed as:   

$$x_{\text{c}}/d_e=\left(\beta /\alpha \right){\left({\rho }_e/{\rho }_a\right)}^{1/2}$$(26)

Allemand et al.²¹⁾ proposed another empirical formula and its expression is as follows:   

$$x_{\text{c}}/d_{\text{e}}={\left({\rho }_{\text{e}}/{\rho }_a\right)}^{1/2}\left(4.2+1.1\left(M_{\text{e}}^2+1-T_{\text{e}}/T_a\right)\right)$$(27)

where x_c is the potential core length, T_e is the nozzle exit temperature, T_a is the ambient temperature, M_e is the exit March number.

Table 3 shows the potential core lengths at different ambient temperatures. The potential core length of the jet increases as the ambient temperature rises and its calculated value of at 1873 K is 2.63 times than that at 285 K. Because the density of ambient gas and the grow rate of turbulent mixing reduce in high-temperature environment which resulting in the jet attenuation being suppressed. The calculated potential core lengths by the proposed model are between the two results of Eqs. (26) and (27). In addition, the effects of ambient temperature T_a and the outlet temperature T_e on the potential core length were taken into account in Eq. (27). The potential core lengths calculated by the temperature corrected model at the four temperatures in Table 3 are close to the calculation results of the empirical formulas in the literature.

Table 3. The comparison of potential core lengths.

Potential core length	285 K	772 K	1002 K	1873 K
xc/de (Ito and Muchi19))	8.61	14.17	16.15	22.08
xc/de (Allemand21))	9.43	16.38	18.80	26.15
xc/de (Present study)	9.20	13.67	16.00	24.23

3.4. Multiple Oxygen Jets at High Ambient Temperature

The characteristic of multiple jets with 11.5 degrees inclination was studied using the corrected model. Figure 9 shows the calculation results of two SST k-ω models at different temperatures. Compared with the calculation results of the original SST k-ω model, the calculated velocity distribution by the corrected model, little changes at room temperature. However, the difference gradually emerges because the grow rate of turbulence mixing reduces at high temperature for the function of the temperature correction.

Fig. 9.

The axial velocity distribution of multiple jets at different temperature calculated with two models; (a) SST k-ω model; (b) corrected SST k-ω model. (Online version in color.)

Figure 10 shows the velocity distribution in the cross section perpendicular to the axis of the lance, at the distance of h=1.2 m from the outlet and Fig. 11 shows the paths of jets. The jets impingement area increases and jets core deflection toward the lance axis weakens as the ambient temperature increases. The greatest degree of deviation occurs at 285 K and the jet center disappears when the vertical distance from the outlet is greater than 1.6 m.

Fig. 10.

The velocity distribution in the cross-sectional of h=1.2 m calculated with the corrected model at different temperature. (Online version in color.)

Fig. 11.

The multiple jets core deflection toward the lance axis at different temperature.

4. Conclusions

In this study, a corrected SST k-ω model, taking the effect of temperature fluctuations on turbulence mixing, was proposed to predict the behavior of the cold supersonic oxygen jet entering into high temperature environments. Combined with the calculated results and the experimental data, the following conclusions are obtained:

(1) The effects of temperature fluctuations on turbulence mixing cannot be neglected when predicting the behavior of the cold oxygen jet flowing through Laval nozzle into high-temperature environment. There are differences between velocity distributions at high temperature gained by experiments in the literature and the calculated results with the previous corrected k-ε models. While the calculated axial velocity and dynamic pressure distributions with the proposed corrected SST k-ω model, at 285 K, 772 K, 1002 K, have a good agreement with the experiment data⁵⁾ and empirical formulas^19,21) in the literature. In addition, the potential core length extends as the ambient temperature increases, and the calculated length at 1873 K is 2.63 times than that at 285 K.

(2) The axial static temperature distributions of the jet at different temperatures remain to further research for the variable Prandtl number and the measurement under non-isothermal conditions.

(3) The interaction of multiple jets and the deflection of the jet core toward the lance axis weaken as the ambient temperature increases.

References

• 1)   K.  Naito,  Y.  Ogawa,  T.  Inomoto,  S.  Kitamura and  M.  Yano: ISIJ Int., 40 (2000), 23.
• 2)   Y.  Tago and  Y.  Higuchi: ISIJ Int., 43 (2003), 209.
• 3)   Y.  Higuchi and  Y.  Tago: ISIJ Int., 41 (2001), 1454.
• 4)   Y.  Higuchi and  Y.  Tago: ISIJ Int., 43 (2003), 1410.
• 5)   I.  Sumi,  Y.  Kishimoto,  Y.  Kikuchi and  H.  Igarashi: ISIJ Int., 46 (2006), 1312.
• 6)   W. J.  Wang,  Z. F.  Yuan,  H.  Matsuura,  H. X.  Zhao,  C.  Dai and  F.  Tsukihashi: ISIJ Int., 50 (2010), 491.
• 7)   H. J.  Odenthal,  U.  Falkenreck and  J.  Schlüter: Proc. European Conf. on Computational Fluid Dynamics, Eccomas Cfd, Netherlands, (2006).
• 8)   M. M.  Li,  Q.  Li,  S. B.  Kuang and  Z. S.  Zou: Steel Res. Int., 86 (2015), 1517.
• 9)   J. C.  Lau,  P. J.  Morris and  M. J.  Fisher: J. Fluid Mech., 93 (1979), 1.
• 10)   D.  Papamoschou and  A.  Roshko: J. Fluid Mech., 197 (1988), 453.
• 11)   S.  Sarkar,  G.  Erlebacher,  M. Y.  Hussaini and  H. O.  Kreiss: J. Fluid Mech., 227 (1991), 473.
• 12)   S.  Sarkar: J. Fluid Mech., 282 (1995), 163.
• 13)   K. S. A.  Hamid,  S. P.  Pao,  S. J.  Massey and  A.  Elimilgui: J. Fluids Eng., 126 (2004), 844.
• 14)   M.  Alam,  J.  Naser and  G.  Brooks: Metall. Mater. Trans. B, 41 (2010), 636.
• 15)   B. E.  Launder and  B. I.  Sharma: Heat Mass Transf., 1 (1974), 131.
• 16)   F. R.  Menter: AIAA J., 32 (1994), 1598.
• 17)   F.  Zhao,  Y. L.  Zhang,  R.  Zhu and  H.  Wang: J. Univ. Sci. Technol. Beijing, 36 (2014), 366.
• 18)   Z. P.  Cai,  Y. S.  Xie and  A. W.  Xia: Iron Steel, 15 (1980), 14.
• 19)   S.  Ito and  I.  Muchi: Tetsu-to-Hagané, 55 (1969), 1152.
• 20)   G.  Kleinstein: J. Spacecr. Rocket., 1 (1964), 403.
• 21)   B.  Allemand,  P.  Bruchet,  C.  Champinot,  S.  Melen and  F.  Porzucek: Microsc. Res. Tech., 98 (2001), 571.