ISSN 0021-8944, Journal of Applied Mechanics and Technical Physics, 2016, Vol. 57, No. 7, pp. 1264–1275. © Pleiades Publishing, Ltd., 2016. Original Russian Text © A.S. Teimurazov, P.G. Frick, 2015, published in Vychislitel’naya Mekhanika Sploshnykh Sred, 2015, Vol. 8, No. 4, pp. 433–444.

Numerical Study of Molten Magnesium Convection in a Titanium Reduction Apparatus A. S. Teimurazov* and P. G. Frick** Institute of Continuous Media Mechanics, Russian Academy of Sciences, Ural Branch, Perm, 614013 Russia e-mail: *[email protected], **[email protected] Received October 5, 2015; in final form, December 30, 2015

Abstract—The structure of a convective flow of molten magnesium in a metallothermic titanium reduction apparatus has been studied numerically for various retort heating and cooling configurations. The mathematical model is based on the thermogravitational convection equations for a singlephase fluid in the Boussinesq approximation. A nonuniform computational grid with a total of 5 million grid points was used. The LES (Large Eddy Simulation) technique was applied to take into account the turbulence. The problem was considered in a three-dimensional nonstationary formulation, which allowed one to construct the instantaneous and average characteristics of the process and to analyze the velocity and temperature pulsation fields. It has been found that steady axisymmetric flows occur at moderate Grashof numbers (Gr ~ 107–108), while unsteady turbulent flows take place at Grashof numbers corresponding to the actual titanium reduction process (Gr ~ 1012). The influence of the degree of nonhomogeneity of the heat release due to the titanium reduction reaction that runs mainly on the magnesium surface has been studied. The computations have been performed for two configurations of the system for maintaining the apparatus heating regime: with furnace heaters operating at full power and with switched-off furnace heaters. Fundamental differences in convective flow structure for these two heating methods have been revealed. It has been established that the most intense velocity and temperature pulsations arise in the region adjacent to the interface between the cooled and heated parts of the retort side surface. Key words: convection, turbulence, low Prandtl numbers, liquid metal, magnesium, metallothermic titanium reduction DOI: 10.1134/S0021894416070129

1. INTRODUCTION The molecular heat transfer in liquid metals with a good thermal conductivity is generally so efficient that the convective heat transfer may be neglected when technological devices are simulated. Turbulent convection in fluids with a low Prandtl number Pr = ν χ (here, ν is the kinematic viscosity and χ is the thermal diffusivity) has a number of peculiarities and is being actively studied by fundamental science, although its experimental studies involve great difficulties and are very rare (see, e.g., the review [1]). However, there are situations where convective flows in a metal can affect significantly the operation of equipment, and either experimental or numerical data are needed to take into account their role in technological processes. For example, in recent years, the tasks of designing the liquid-metal cooling systems of reactors have stimulated the development of computational codes to simulate turbulent convective flows of liquid metals [2, 3] and the execution of experiments on turbulent convection in liquid sodium in long cylinders with various orientations relative to gravity [4–7]. The metallothermic titanium reduction process is another example of a technological process in which thermogravitational convection plays an important role [8, 9]. The retort for titanium reduction is a cylindrical vessel with a diameter up to 2 m and a height up to 4 m containing liquid magnesium at a temperature of 850°C. Titanium tetrachloride (TTC) is supplied to the magnesium surface. The chemical reaction runs predominantly on the surface to produce spongy titanium and magnesium dichloride that sink to the retort bottom:

2Mg + TiCl 4 = Ti + 2MgCl 2 + Q . The titanium reduction process takes more than two days and largely remains a “black box,” because the possibilities for any measurements during the process are very limited (only the temperature of the outer 1264

NUMERICAL STUDY OF MOLTEN MAGNESIUM CONVECTION

1265

retort surface is measured at several points). A stable reaction should be provided by controlling the supply of titanium tetrachloride, draining the magnesium chloride, and cooling and heating the furnace. At present, up to 5% of the production cycles are rejected due to the contingencies associated with the disruption of the magnesium chloride settling process, the emergence of titanium sponge on the magnesium surface, and local retort overheating. For decades, titanium sponge production engineers have associated the main hopes for controlling the reaction progress with the magnesium level measurements in the reactor, for which purpose various indirect methods were proposed (see, e.g., [10]). However, attempts to use them in actual production have shown that they work only at the initial stages of the reaction without providing a stable reaction before cycle completion. It follows from a generalization of the results of these works that the reaction control is not possible without a thorough understanding of not only the chemistry but also the hydrodynamics of the process in the reactor. Since the reaction is accompanied by the release of a large amount of heat (Q = 1707 kJ per 1 kg of titanium tetrachloride [11, 12]), the radial temperature gradient can reach hundreds of degrees per meter. The retort is placed in a furnace that together with a blower creates the thermal regime of the reaction. Temperature gradients cause convective flows in the body of magnesium whose influence on the titanium reduction process remains the subject of discussion by specialists. It is a change in the pattern of convective flow that may be responsible for the so-called “non-separation,” whereby the reaction by-product (the magnesium salt whose density is close to the magnesium density) ceases to sink to the reactor bottom and the process turns out to be spoiled. Direct measurements of process characteristics involve considerable difficulties due to the high operating temperature, large masses and sizes of the facility. This determines the interest in the numerical simulation of flows in a reactor. Knowing their structure at various stages of the process, a metal flow can be formed so as to provide a more stable reaction. The first attempts to determine the intensity of convective flows in a reactor based on numerical simulations of magnesium convection in a cylindrical vessel were made almost 40 years ago [13, 14]. The simulations were performed for greatly lowered Grashof numbers, 10 5 instead of 1012 , which forces us to treat the velocity estimates made in a real reactor with caution. Recently, the results of numerical simulations of magnesium convection in a cylinder were obtained for realistic control parameters [15], but the computations were performed for the stationary axisymmetric case. The goal of this paper is to construct a hydrodynamic model of the process in a three-dimensional formulation of the problem with a proper allowance for the turbulent flow and to investigate the convection in the reactor for various heating and cooling regimes using this model. 2. FORMULATION OF THE PROBLEM AND THE MATHEMATICAL MODEL We consider a convective flow in the reactor of only one phase, liquid magnesium, while the contribution of the second liquid phase, magnesium chloride, and the porous medium, titanium sponge, which sink to the reactor bottom, is disregarded. Our three-dimensional mathematical model is based on the thermogravitational convection equations for a single-phase fluid in the Boussinesq approximation. We use the LES (Large Eddy Simulation) approach to take into account the turbulence, namely Smagorinsky’s model [16]. Within this approach, the large turbulent flow scales are calculated explicitly, while the contribution of smaller-scale eddies is taken into account implicitly through subgrid-scale closure. According to the latter, small-scale turbulence is parameterized via large-scale flow characteristics using empirical coefficients. The equations for a large-scale flow are obtained by applying the procedure of filtering the thermal convection equations. The length scale Δ for the filter is determined by the computational grid size, i.e., 13 Δ = ( hx hy hz ) , where hx , hy , hz is the grid step size along each of the axes of the Cartesian coordinate system. The eddies with a size greater than Δ are deemed “large,” while the eddies with a size less than Δ are deemed “small,” the contribution from which should be simulated. The system of equations describing a large-scale flow is

∂v + v ⋅ ∇ v = − 1 ∇ P + ν Δ v + gβ Te , eff z ∂t ρ ∇ ⋅ v = 0, ∂T + v ⋅ ∇ T = χ Δ T . eff ∂t JOURNAL OF APPLIED MECHANICS AND TECHNICAL PHYSICS

Vol. 57

No. 7

2016

1266

TEIMURAZOV, FRICK

(a)

z

(b)

z

(c)

ΓU ΓSU

h H

ΓSL 0

R

x

x

0

y

y

ΓL

Fig. 1. Scheme of the computational domain: a general view (a), the vertical xz section with the designated boundaries (b); the computational grid structure (c).

Here, t is the time, v is the fluid velocity vector, P is the pressure (the deviation from the hydrostatic pressure P0), T is the temperature (the deviation from the mean T0 ), ρ is the mean density, g is the acceleration of gravity, β is the thermal volume expansion coefficient, and e z is a unit vector directed along the z axis. Instead of the kinematic viscosity ν and thermal diffusivity χ , the equations contain the effective viscosity ν eff and effective thermal diffusivity χeff , respectively. The effective viscosity is defined by the expression 2 2 ν eff = ν + ν t , where ν t = C S Δ S is the turbulent viscosity (CS is Smagorinsky’s constant, S is the norm

(

of the strain rate tensor S = ∇ v + (∇ v )

T

)

2 ). The effective thermal diffusivity is found from the formula

χ eff = ν t Pr t + ν Pr , where Prt is the turbulent Prandtl number. The following values are taken in our computations: C S = 0.14 and Prt = 0.9 . The retort for titanium reduction has a radius of 0.75 m and is filled with liquid magnesium to a level of 2.5 m. The upper part of the side surface is blown with cold air; the height of the cooling zone is 0.7 m. The scheme of the computational domain and the coordinate system are shown in Fig. 1a. The computational domain is a cylinder with radius R , height H , and a cooled upper part of the side surface with height h. Figure 1b designates the following: Γ U is the upper surface with area S U , Γ SU is the part of the cooled side surface with area S SU , Γ SL is the lower part of the side surface with area S SL , and Γ L is the base with area S L . All boundaries of the domain are assumed to be solid; the no-slip conditions are specified for the velocity at all boundaries ( v = 0 ). The form of the condition for the velocity at the free upper boundary Γ U deserves a separate discussion. Since films emerge on a free surface in liquid metals [17], as a rule, it is assumed in calculations that there is sticking here. However, the fulfilment of this condition in the case of a chemical reaction is by no means obvious. To clarify the question about the influence of the form of boundary conditions for the velocity on the flow structure, we performed a number of computations with slip at the upper boundary. It turned out that, in this case, the flow near the free surface is slightly more intense than the flow with sticking, but the flow structure does not change fundamentally. For this reason, our computations here were performed using the sticking conditions for the velocity at the upper boundary. In our computations we consider three different forms of thermal boundary conditions (TBCs), three heating configurations. A situation whereby a reaction with a thermal power Q R runs on the metal surface, the metal body in the upper part is cooled with air through the side surface, and the furnace heaters are switched off corresponds to the first configuration, TBC-1 (Fig. 2a). In this case, a constant heat flux qU is specified at the upper boundary of the cylinder Γ U , a constant heat flux − qSU is removed from the side JOURNAL OF APPLIED MECHANICS AND TECHNICAL PHYSICS

Vol. 57

No. 7

2016

NUMERICAL STUDY OF MOLTEN MAGNESIUM CONVECTION

(a)

qU

(b)

1267

(c)

qU

qU −qSU

−qSU

∇T = 0

∇T = 0

−qSU

∇T = 0

qSL

∇T = 0

qL

Fig. 2. Schemes of the heat flux distribution at the boundaries of the computational domain for various heating configurations: (a) TBC-1, uniform heating from the reaction; (b) TBC-2, nonuniform heating from the reaction according to a parabolic law; (c) TBC-3, heating from the reaction and from the furnace heaters.

surface Γ SU , the lower part of the side surface Γ SL and the base ΓL are assumed to be thermally insulated. Thus, TBC-1 are

QR QR –k ∂ T = qU = , – k ∂T = −qSU = − , ∂T = 0. SU S SU ∂n ΓU ∂ n ΓSU ∂ n ΓSL , Γ L The second configuration, TBC-2 (Fig. 2b), is analogous to the first one, except that the heat flux at the upper boundary Γ U is distributed not uniformly with parabolic law with a maximum at the center and a minimum on the periphery. TBC-2 are taken to investigate the influence of a nonuniform heat flux distribution at the surface of molten metal Γ U on the structure of the convective motion in the retort. Such a temperature nonuniformity can emerge in the reactor during the actual titanium reduction process due to nonuniform titanium tetrachloride sputtering on the magnesium surface. Consequently, TBC-2 are –k ∂ T ∂n

ΓU

= qU =

2QR ⎛ ( x 2 + y 2 ) ⎞ ⎜1 − ⎟, SU ⎝ R2 ⎠

– k ∂T ∂n

Γ SU

= −qSU = −

QR , S SU

∂T ∂n

Γ SL , Γ L

= 0.

The third configuration, TBC-3 (Fig. 2c), combines the heating from above at the boundary Γ U with a constant heat flux qU from a chemical reaction with a thermal power Q R , the heating from below at the boundary Γ L with a heat flux qL from the lower heater with a power QL , and the heating of the lower part of the side surface Γ SL with a heat flux qSL from the side heaters with a total power QSL . All of this heat is removed from the retort by cooling the boundary Γ SL . Thus, TBC-3 are

(Q + QSL + QL ) = −qSU = − R – k ∂T , ∂ n ΓSU S SU ΓU Q Q = qSL = SL , –k ∂ T = qL = L . –k ∂ T ∂ n ΓSL ∂n ΓL S SL SL The first and third configurations correspond to two extreme situations: the first is the heat release only as a result of the reaction with switched-off furnace heaters and the third is the heat supply from all of the simultaneously operating furnace heaters. The configurations considered are most revealing for qualitative estimates of the flow structure, although it should be noted that other heating regimes are also possible in the actual titanium reduction process; in particular, only separate sections of the furnace heater can be switched on simultaneously. In all our computations we checked whether the integral thermal balance in the reactor was kept. The governing parameters in our problem are the Grashof number Gr and the Prandtl number Pr . Just as in [13], the Grashof number was found via the heat flux at the upper boundary: Gr = gβ qU R 4 ( ν 2k ) , where the characteristic temperature difference is qU R k , k is the thermal conductivity of the liquid. This method of determining the Grashof number is justified for TBC-1 and TBC-2 presented in Figs. 2a and 2b, but it disregards the contribution from the furnace heaters in the configuration corresponding to –k ∂ T ∂n

= qU =

QR , SU

JOURNAL OF APPLIED MECHANICS AND TECHNICAL PHYSICS

Vol. 57

No. 7

2016

1268

z/H

TEIMURAZOV, FRICK

z/H

(а)

0.8 0.6

V 1000 800

0.4

500 300

0.6

V 1190 900

0.4

600 300

0.8

0 0.2

0.8

z/H T 0.474 0.341 0.210 0.078

0.6 0.4

0.6

1400 900 450 0

0.4

−0.8 −0.4 0 0.4 0.8 x/R z/H

(e)

0.8

T 0.676

0.6

0.329 0.197

0.4

(f)

0.8

T 0.312

0.6

0.143

0.4

0.121

0.091 0.2

−0.8 −0.4 0 0.4 0.8 x/R

V 1930

0.2

0.083 0.2

0.8

−0.8 −0.4 0 0.4 0.8 x/R

(d)

(c)

0 0.2

−0.8 −0.4 0 0.4 0.8 x/R z/H

z/H

(b)

0.410

0.2

−0.8 −0.4 0 0.4 0.8 x/R

−0.8 −0.4 0 0.4 0.8 x/R

Fig. 3. Velocity (a–c) and temperature (d–f) fields in the xz section at Gr = 2.2 × 10 7 corresponding to various heating configurations: TBC-1 (a, d), TBC-2 (b, e), and TBC-3 (c, f).

TBC-3 (Fig. 2c). Below, we will use both dimensional and dimensionless quantities when discussing the results. The equations are discretized by the finite volume method. The PISO (Pressure Implicit with Splitting of Operators) algorithm [18] based on the pressure correction procedure is used to solve the derived system of equations. The computational grid contains collocated grid points, i.e., the values of all variables are calculated at the same grid points [19]. The terms with time derivatives are discretized using an implicit Euler scheme. The convective terms are calculated based on the TVD (Total Variation Diminishing) scheme and the flux limiter proposed by Sweby [20]. This scheme possesses the transportiveness and boundedness properties and has a second-order approximation accuracy [21]. The discrete analogs of the diffusion terms are constructed on the basis of a computational scheme with central differences. In total, the implicit conservative computational scheme represents the original system of differential equations with a first-order accuracy in time and a second-order accuracy in coordinate. The Biconjugate Gradient Method (BiCG) [22] is used to solve the system of linear algebraic equations. In all of the main computations the grid has a regular block structure with a total of 5 million grid points and is nonuniform (Fig. 1c) with the minimum step near the boundary (3 mm) and a large spatial step in the inner part of the retort (up to 10 mm). Grid refinement near the boundary is needed to properly resolve the temperature and velocity boundary layers. The Courant number in our computations does not exceed 0.5. The free and open source software package OpenFOAM Extend 3.1 is used to implement the described numerical scheme. The computations were performed on the “Triton” computer cluster at the Institute of Continuous Media Mechanics of the Ural Branch of the Russian Academy of Sciences (Perm) and on the “Uran” supercomputer of the IMM of the Ural Branch of the Russian Academy of Sciences (Yekaterinburg). The acceleration of the operation of the parallel algorithm relative to the computing time on a single processor was: a factor of 5.5, 14.2, 43.8, and 62.6 on 8, 32, 128, and 160 processors, respectively. It took 240 s of the computer time on 160 processors to compute 1 s of the time of the physical process at Gr = 2.2 × 1012 . An JOURNAL OF APPLIED MECHANICS AND TECHNICAL PHYSICS

Vol. 57

No. 7

2016

NUMERICAL STUDY OF MOLTEN MAGNESIUM CONVECTION

z/H

(а)

(b) V 2337 1850

0.8 0.6

1200 0.4 0.2

y/R 0.8 0.4 0

600

−0.4

0

−0.8

z/H

V 1950 1700 1080

0.8

550 0

0.4

−0.8−0.4 0 0.4 0.8 x/R

1269

(c) T 0.394

0.6 0.183 0.064 0.134 0.2

−0.8 −0.4 0 0.4 0.8 x/R

−0.8 −0.4 0 0.4 0.8 x/R

Fig. 4. Results for the TBC-1 heating configuration at Gr = 2.2 × 10 8 : the instantaneous velocity fields in the xz section (a) and the xy section at height z = 0.8H (b); the temperature field in the xz section (c).

averaging interval of at least 2000 s is needed to properly estimate the mean fields and statistical characteristics of the pulsations. 3. RESULTS OF COMPUTATIONS In our computations the main specified parameters were the heating configurations (Fig. 2) and the input and output heat fluxes to which the Grashof number is related. The computations were performed for a wide range of Grashof numbers, from Gr = 2.2 × 10 7 (under low heating) to Gr = 2.2 × 1012 (for the actual process in the reactor). Such heating regimes were chosen to trace how the flow structure changed qualitatively with heating rate and to determine the boundaries of applicability of the axisymmetric stationary mathematical models for the problem under consideration. The convective parameters of the fluid corresponded to liquid magnesium at a temperature of 850°С and were the following: the Prandtl number Pr = 0.008, the kinematic viscosity ν = 4.87 × 10 −7 m2 s–1, and the thermal expansion coefficient −4 β = 1.7 × 10 K–1 [23]. Figure 3 presents the results of our computations for the smallest Grashof number considered, Gr = 2.2 × 10 7 , for the specified heating configurations. Under uniform heating from above (TBC-1) a steady axisymmetric flow with a ring vortex in the upper part is established in the retort (Fig. 3a, here and below, ν R is chosen as a characteristic unit for the velocity). It should be noted that the vortex is not immediately adjacent to the upper boundary of the domain but is slightly offset from it. It can be seen from Fig. 3d, which displays the temperature field for the case in question, that a significant horizontal temperature gradient emerges in the boundary layer, while the isotherms in the central part of the cylinder are nearly horizontal. Such a structure of the temperature field qualitatively agrees with the results from [13] for the heating from above. Under nonuniform heating (TBC-2) the vortex in the upper third of the metal body is adjacent to the upper boundary and is more intense (Fig. 3b). This difference is explained by the fact that nonuniform heating produces an additional horizontal temperature gradient near the upper boundary (Fig. 3e). The situation changes significantly if in addition to the heating from above, there is heating in the lower part and on the side surface of the retort (TBC-3). The motion ceases to be steady and, therefore, Figs. 3c and 3f present the time-averaged velocity and temperature fields. The flow acquires a two-vortex structure (Fig. 3c), the upper and lower vortices have similar sizes, the maximum velocities increase twofold compared to the velocities in the case without any additional heaters. Let us trace the change of the flow structure as the Grashof number gradually increases by considering only the heating from above (TBC-1). The flow becomes unsteady and loses its axial symmetry with increasing Grashof number. Figures 4a and 4b show examples of the instantaneous velocity fields for Gr = 2.2 × 10 8 . The vortex intensity in the upper part of the cylinder rose, the vortex itself slightly shifted downward, and the temperature boundary layer became thinner (Fig. 4c). It should be noted that finding JOURNAL OF APPLIED MECHANICS AND TECHNICAL PHYSICS

Vol. 57

No. 7

2016

1270

TEIMURAZOV, FRICK

z/H

(b)

(а)

0.8 V 10 014

0.6

7700 4600 3080 0

0.4 0.2

z/H

V 6322 4600 3080 0

−0.8 −0.4 0 0.4 0.8 x/R (с)

(d) z/H 0.8

0.8 0.6

T 0.364

0.4

0.193 0.087

0.2

−0.8 −0.4 0 0.4 0.8 x/R

0.6

T 0.008

0.4 −0.8 −0.4 0 0.4 0.8 x/R

–0.124

0.071 0.045 0.124

Fig. 5. Results for the TBC-1 heating configuration at Gr = 1.1 × 1010 in the xz section: the instantaneous velocity field (a) and its magnified fragment (b); the temperature field (c) and its magnified fragment (d).

z, m 2.0 1.5

z, m

(а)

V, m/s 0.059 2.0 0.040 1.5 0.030 0.010 1.0 0

1.0

(b)

(c) V, m/s 0.05 0.03 0.01 0

0.5

0.5

y, m

T, °C 1172 1157 1142 1127 1116

0.6 0.3 0 −0.4 −0.6 −0.6−0.3 0 0.3 0.6 x, m

−0.6 −0.3 0 0.3 0.6 x, m

−0.6 −0.3 0 0.3 0.6 x, m

(d) z, m 2.4

T, °C 1172

2.0

1067 977 887 798

1.6 1.2 −0.6 −0.3 0 0.3 0.6 x, m

Fig. 6. Results for the TBC-2 heating configuration at Gr = 2.2 × 1012 : the instantaneous (a) and time-averaged (b) velocity fields in the xz section; the temperature field near the upper boundary of the computational domain in the xy section (c); the temperature field in the xz section in the upper part of the magnesium body (d). JOURNAL OF APPLIED MECHANICS AND TECHNICAL PHYSICS

Vol. 57

No. 7

2016

NUMERICAL STUDY OF MOLTEN MAGNESIUM CONVECTION z, m

z, m

(а) T, °C 945 917 882 847

2.0 1.5 1.0

2.0 1.5

1.0

m /s 0.011 0.009

0.003

2

0.5

2.0

−0.6 −0.3 0 0.3 0.6 x, m

(e)

(f) PI Ex ,

2

m /s

2

y, m PI PI 2 2 0.00157 0.6 Ex + Ey , m /s 0.00120 0.00254 0.3 0.00180 0.00080 0 0.00090 0.00040 −0.3 0 0 −0.6 −0.6−0.3 0 0.3 0.6 x, m

1.5 1.0

0 0.5

−0.6 −0.3 0 0.3 0.6 x, m

0

0

0.5

z, m

0.5

0.05 1.0

0.05

0.006 1.0

1.5

−0.6 −0.3 0 0.3 0.6 x, m 2

V, m/s 0.12 0.10

0.10

(d) P Ez ,

(c)

2.0

0.18 0.15

1.5

−0.6 −0.3 0 0.3 0.6 x, m z, m

V, m/s

2.0

789

0.5

z, m

(b)

1271

−0.6 −0.3 0 0.3 0.6 x, m

Fig. 7. Results for the TBC-3 heating configuration at Gr = 2.2 × 1012 : (a) the instantaneous temperature field, (b) the instantaneous velocity field, (c) the time-averaged velocity field, (d) the vertical velocity pulsation component E zP , (e) the horizontal velocity pulsation component E xP in the xz section; (f) the velocity pulsation field ( E xP + E yP ) in the horizontal xy plane at z = 0.72H.

the exact boundary at which the flow becomes unsteady and loses its axial symmetry was not among the main goals of our paper, because this boundary is in a region with control parameters that are definitely lower than those in the actual titanium reduction apparatus. An increase in the Grashof number is accompanied by a further growth in the flow velocity and the vortex size. Figure 5 presents the results of our computations for Gr = 1.1 × 1010 . The motion acquired a distinct boundary-layer pattern. The convective flow in the upper third of the cylinder, except the thin boundary layer, virtually disappeared (Figs. 5a, 5b). This is because a nearly uniform vertical temperature gradient without any nonuniformity in the horizontal direction was established here (Fig. 5c). A significant horizontal temperature gradient emerged only near the interface between the cooled and thermally insulated parts of the side surface (Fig. 5d). At the same time, the flow in the lower part of the retort remained fairly intense and unsteady. A Grashof number Gr ≈ 2.2 × 1012 corresponds to the actual titanium reduction process in the reactor. Figures 6a and 6b present the velocity fields for such a flow regime under nonuniform heating from above (TBC-2). The flow structure does not change qualitatively compared to the case of Gr = 1.1 × 1010 ; the flow in the upper third of the computational domain is still suppressed, but the intensity of the motion in the lower part increases, the turbulent flow occupies the entire volume below the cooling zone. In dimensional units for a retort with radius R = 0.75 m filled with a metal to height H = 2.5 m and with height h = 0.7 m of the cooled part of the side surface, the flow velocities in the boundary layer reach 6 cm s–1, while the velocity in the central part of the retort is ≈ 4 cm s–1. The ring vortex in the lower part of the retort is clearly discernible in Fig. 6b, which presents the time-averaged velocity field for the regime under conJOURNAL OF APPLIED MECHANICS AND TECHNICAL PHYSICS

Vol. 57

No. 7

2016

1272

TEIMURAZOV, FRICK

(a)

T, °C 840

z 3

1

820 800 780 760

1

2

2

3

2600

2800

3000

3200

0

x z (b)

845 835

1

1

2

825

810

2 2600

3

3 2800

3000

0

3200

x

t, s Fig. 8. Results for the TBC-3 heating configuration at Gr = 2.2 × 1012 : the time dependences of the temperature at specific points of the computational domain located at a distance of 1 cm from the side wall (a) and on the axis (b) (see the schemes on the right); lines 1 , 2 , and 3 correspond to the points at z = 2.3 , 1.8, and 0.3 m, respectively.

sideration. The velocities reach their maximum values in the boundary layer near the side wall at a height of 0.7–1.7 m. It is important that the nonuniformity of the heat flux distribution at the upper boundary of the computational domain in the regime under consideration has no fundamental influence on the flow structure: the isotherms in its upper part are virtually horizontal already at a distance of 5 cm from the surface (Figs. 6c and 6d). It is worth noting that it is natural to expect the emergence of the most intense vortices in the region of the greatest temperature gradients (in the upper third of the metal body). However, this is not observed at the actual parameters of the process. The temperature gradient caused by the surface reaction was shown to actually maintain the ring vortex in the upper part of the cylinder at low Grashof numbers and to suppress it at Grashof numbers corresponding to the actual titanium reduction process. Turning on additional heating of the lower part of the retort changes significantly the flow pattern. Figure 7 presents the results of our computations with the TBC-3 heating configuration at the actual parameters of the process in the reactor, when the heaters operate at full power equal to 423 kW, out of which the heating from below is QL = 94 kW and the heating on the side surface in its lower part is QSL = 329 kW; the thermal power of the surface reaction is QR = 205 kW. No uniform vertical temperature gradient is established in the upper part of the cylinder for such a thermal regime. The temperature field (Fig. 7a) reflects the unsteady pattern of the motion caused by the near-wall upflow from the heated lower part of the reactor and by the downflow arising in the cooled part of the retort. The boundary between these flows is not stationary and shifts upward-downward. It can be seen from the instantaneous velocity fields (Fig. 7b) that the turbulent flow occupies the entire volume of the molten metal. A two-vortex structure of the convective flow, one ring vortex in the upper part of the computational domain and one in its lower part, is clearly discernible in the time-averaged velocity field (Fig. 7c). The interface between the vortices is at a height of 0.72 H, i.e., at the boundary between the cooled and heated parts of the cylinder side surface. It can be seen that the vortex centers are shifted to the side boundary; the vortices reach their maximum intensity (12 cm s–1) near the boundaries. The flow intensity is low in the central part of the upper third of the cylinder, while an upflow with an average velocity of ≈ 7 cm s–1 is observed in the remaining volume. The maximum instantaneous velocity V max near the boundary reaches 20 cm s–1. JOURNAL OF APPLIED MECHANICS AND TECHNICAL PHYSICS

Vol. 57

No. 7

2016

NUMERICAL STUDY OF MOLTEN MAGNESIUM CONVECTION

EVz, m2 s−1 10−6 1 10−7 10−8

(a) −5/3

2

ET, K2 s 10−1 10−2

10−4

10−10

10−5

10−11

10−6 0.02 0.05 0.10 0.20 0.50 1.00 2.00 ω, Hz

(b) −5/3

z

1 4

10−3

10−9

1273

2

2

4

1

3

3

0

x

0.02 0.05 0.10 0.20 0.50 1.00 2.00 ω, Hz

Fig. 9. Results for the TBC-3 heating configuration at Gr = 2.2 × 1012 : the power spectral densities of the vertical velocity V z (a) and temperature (b) pulsations at specific points of the computational domain (see the scheme on the right): the points located at 1 cm from the side wall at heights z = 1.8 and 2.3 m correspond to lines 1 and 2; the points on the axis at heights z = 1.8 and 2.3 m correspond to lines 3 and 4 .

The distribution of velocity pulsations against a background of the mean flow is also of interest. Figure 7d demonstrates the vertical velocity V z pulsation field E zP in the vertical yz section of the computational domain. It follows from the figure that the largest pulsations take place in a narrow ring-like region ≈ 7 cm in thickness near the interface between the cooled and heated parts of the side surface. Figure 7e presents the horizontal velocity V x pulsation field E xP in the vertical yz section. The horizontal velocity pulsations also have maxima near the interface, but the regions with a significant pulsation amplitude are distributed over a larger volume. It should be noted that the vertical velocity pulsations are more intense than the horizontal ones (the maximum values differ by a factor of 7). Figure 7f presents the velocity pulsation field ( E xP + E yP ) in the xy plane at height z = 0.72H . In the horizontal section the pulsations are distributed along the radius relatively uniformly; their maximum values lie within the ring 25 cm < r < 65 cm. Figure 8 gives examples of the time dependence of the temperature at specific points of the computational domain. The most intense temperature pulsations arise in the vicinity of the side wall near the interface between the cooled and heated parts of the side surface. At the point located at a distance of 1 cm from the side surface and at a height of 2.3 m, the maximum temperature pulsation amplitudes reach 50°С (Fig. 8a). The pulsation intensity depends on time: the periods with large and small amplitudes alternate (Fig. 8a, line 2). The amplitudes in the central part of the cylinder were found to be considerably smaller than those near the side surface. At points on the retort axis the largest temperature pulsation amplitude is observed near the upper boundary (Fig. 8b, line 1), and the farther the point from the upper boundary, the smaller the pulsation amplitude at it. Figure 9a presents the power spectral densities of the vertical velocity V z pulsations at various points of the computational domain. Its pulsations near the side surface are seen to be more intense than the pulsations at the center in the entire frequency range. Figure 9b shows the power spectral densities of the temperature pulsations. At a height of 1.8 m the temperature pulsations near the side boundary and at the center differ by two orders of magnitude (lines 1 and 2), while at a height of 2.3 m the pulsation spectra near the boundary and at the center virtually coincide (lines 3 and 4). In both parts of Fig. 9 the straight line indicates a slope of − 5/3, which corresponds to the Kolmogorov law. The slope is shown for reference, because the interpretations in terms of homogeneous and isotropic turbulence, especially near the boundaries of the computational domain, require great caution, despite the fact that such a slope is observed in a frequency range exceeding one decade. The presented characteristics of the turbulent pulsations show that the mean flow is formed against a background of well-developed turbulence, and the instantaneous velocity and temperature fields differ significantly from the mean ones. In these conditions using a full three-dimensional model becomes crucial. Although the mean fields turn out to be axisymmetric, the three-dimensionality of the problem is JOURNAL OF APPLIED MECHANICS AND TECHNICAL PHYSICS

Vol. 57

No. 7

2016

1274

TEIMURAZOV, FRICK

needed to properly take into account the turbulence, because the LES models were developed precisely for the three-dimensional case. Computations in a stationary axisymmetric formulation [15] also showed that with switched-on heaters the most intense convective motion is caused by a temperature difference near the interface between the cooled and heated parts of the cylinder side surface, while the overall structure of the averaged flow is a two-vortex one. However, the shape and sizes of the vortices in computations in a full formulation differ significantly, showing that the flow has a distinct boundary-layer pattern at large Gr , and the vortex centers shift toward the retort side boundary. 4. CONCLUSIONS The structure of a convective flow of molten magnesium in a titanium reduction reactor was numerically studied for a wide range of Grashof numbers (from Gr = 2.2 × 10 7 to 2.2 × 1012 ). The problem was solved in a three-dimensional nonstationary formulation, which allowed us to calculate the instantaneous and average characteristics of the process and to analyze the velocity and temperature pulsation fields. Steady and axisymmetric flows in the reactor were found to be possible only at considerably lowered values of the control parameters. At parameters corresponding to the actual process, unsteady turbulent flows arise in the retort. The temperature gradient caused by the chemical reaction on the surface maintains a vortex ring in the upper part of the metal body at low Grashof numbers but suppresses it at Grashof numbers corresponding to the reduction process in the reactor. The computations were performed for two heating configurations in the reactor: with furnace heaters operating at full power and with switched-off ones. The convective flow structures in these two cases were shown to differ fundamentally. Consequently, by changing the retort heating regimes, we can influence the flow to create the most favorable conditions for the reaction to run. The maximum flow velocity in the reactor was estimated. The flow velocities reach the greatest values near the side boundaries, with the maximum instantaneous and mean velocities being V max ≈ 20 cm s–1 and V max ≈ 12 cm s–1, respectively. The most intense velocity and temperature pulsations arise in the region near the interface between the cooled and heated parts of the retort side surface with switched-off heaters. ACKNOWLEDGMENTS This work was financially supported by the Basic Research Program of the Ural Branch of the Russian Academy of Sciences (project no. 15-10-1-9). REFERENCES 1. Ahlers, G., Grossmann, S., and Lohse, D., Heat transfer and large scale dynamics in turbulent Rayleigh– Bénard convection, Rev. Mod. Phys., 2009, vol. 81, pp. 503–537. 2. Zaitsev, A.M., Semenov, V.N., and Shvetsov, Y.E., Simulation of mixing different temperature jets using CABARET method, Vychisl. Mekh. Sploshnykh Sred, 2013, vol. 6, no. 4, pp. 430–437. 3. Rogozhkin, S.A., Aksenov, A.A., Zhluktov, S.V., Osipov, S.L., Sazonova, M.L., Fadeev, I.D., Shepelev, S.F., and Shmelev, V.V., Development and verification of a turbulent heat transport model for sodium-based liquid metal coolants, Vychisl. Mekh. Sploshnykh Sred, 2014, vol. 7, no. 3, pp. 306–316. 4. Frick, P., Khalilov, R., Kolesnichenko, I., Mamykin, A., Pakholkov, V., Pavlinov, A., and Rogozhkin, S., Turbulent convective heat transfer in a long cylinder with liquid sodium, Europhys. Lett., 2015, vol. 109, no. 1, p. 14002. 5. Kolesnichenko, I.V., Mamykin, A.D., Pavlinov, A.M., Pakholkov, V.V., Rogozhkin, S.A., Frick, P.G., Khalilov, R.I., and Shepelev, S.F., Experimental study on free convection of sodium in a long cylinder, Therm. Eng., 2015, vol. 62, no. 6, pp. 414–422. 6. Vasil'ev, A.Yu., Kolesnichenko, I.V., Mamykin, A.D., Frick, P.G., Khalilov, R.I., Rogozhkin, S.A., and Pakholkov, V.V., Turbulent convective heat transfer in an inclined tube filled with sodium, Tech. Phys., 2015, vol. 60, no. 9, pp. 1305–1309. 7. Mamykin, A., Frick, P., Khalilov, R., Kolesnichenko, I., Pakholkov, V., Rogozhkin, S., and Vasiliev, A., Turbulent convective heat transfer in an inclined tube with liquid sodium, Magnetohydrodynamics, 2015, vol. 51, no. 2, pp. 329–336. JOURNAL OF APPLIED MECHANICS AND TECHNICAL PHYSICS

Vol. 57

No. 7

2016

NUMERICAL STUDY OF MOLTEN MAGNESIUM CONVECTION

1275

8. Garmata, V.A., Gulyanitskiy, B.S., Kramnik, V.Yu., Lipkes, Ya.M., Seryakov, G.V., Suchkov, A.B., and Homyakov, P.P., The Metallurgy of Titanium, USA: National Technical Information Service, 1970. http://www.dtic.mil/dtic/tr/fulltext/u2/717831.pdf. 9. Garmata, V.A., Petrun’ko, A.N., Galitskiy, N.V., Olesov, Yu.G., and Sandler, R.A., Titan (Titanium), Moscow: Metallurgiya, 1983. 10. Khalilov, R.I., Khripchenko, S.Yu., Frik, P.G., and Stepanov, R.A., Electromagnetic measurements of the level of a liquid metal in closed volumes, Meas. Tech., 2007, vol. 50, no. 8, pp. 861–866. 11. Sergeev, V.V., Galitskiy, N.V., Kiselev, V.P., and Kozlov, V.M., Metallurgiya titana (The Metallurgy of Titanium), Moscow: Metallurgiya, 1971. 12. Mal’shin, V.M., Zavadovskaya, V.N., and Pampushko, N.A., Metallurgiya titana (The Metallurgy of Titanium), Moscow: Metallurgiya, 1991. 13. Tarunin, E.L., Shihov, V.M., and Yurkov, Yu.S., Free convection in a cylindrical vessel with a given heat flow at the upper boundary, Uch. Zap. Perm. Univ., Gidrodinam., 1975, no. 327 (6), pp. 85–98. 14. Zimin, V.D., Lyakhov, Yu.N., and Sorokin, M.P., Convection in a vertical cylinder under heating from above, Uch. Zap. Perm. Univ., Gidrodinam., 1975, no. 327 (6), pp. 73–84. 15. Tsaplin, A.I. and Nechaev, V.N., Numerical modeling of non-equilibrium heat and mass transfer processes in a reactor for the production of porous titanium, Vychisl. Mekh. Sploshnykh Sred, 2013, vol. 6, no. 4, pp. 483–490. 16. Smagorinsky, J., General circulation experiments with the primitive equations. I. The basic experiment, Mon. Weather Rev., 1963, vol. 91, no. 3, pp. 99–164. 17. Kolesnichenko, I. and Khripchenko, S., Surface instability of the plane layer of conducting liquid, Magnetohydrodynamics, 2003, vol. 39, no. 4, pp. 427–434. 18. Issa, R.I., Solution of the implicitly discretised fluid flow equations by operator-splitting, J. Comput. Phys., 1985, vol. 62, no. 1, pp. 40–65. 19. Ferziger, J.H. and Peric, M., Computational Methods for Fluid Dynamics, Berlin: Springer, 2002. 20. Sweby, P.K., High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Numer. Anal., 1984, vol. 21, no. 5, pp. 995–1011. 21. Versteeg, H.K. and Malalasekera, W., An Introduction to Computational Fluid Dynamics: The Finite Volume Method, Harlow: Pearson Education, 2007. 22. Fletcher, R., Conjugate gradient methods for indefinite systems, Lect. Notes Math., 1976, vol. 506, pp. 73–89. 23. Eidenzon, M.A., Magnij (Magnesium), Moscow: Metallurgiya, 1969.

Translated by V. Astakhov

JOURNAL OF APPLIED MECHANICS AND TECHNICAL PHYSICS

Vol. 57

No. 7

2016