- Email: [email protected]
Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder
ARTICLE IN PRESS
JID: SPJPM
[m3+dc;June 22, 2017;11:9]
Available online at www.sciencedirect.com
St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12 www.elsevier.com/locate/spjpm
Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder✩ Sergey I. Smirnov∗, Evgueni M. Smirnov, Alexander A. Smirnovsky Peter the Great St. Petersburg Polytechnic University, 29 Politekhnicheskaya St., St. Petersburg 195251, Russian Federation Available online xxx
Abstract The obtained results of direct numerical simulation of the free mercury convection in a rotating cylindrical container heated from below are presented. Setting the Prandtl number equal to 0.025 and the height-to-diameter ratio equal to 1.0, effects of container rotation and heat transfer in horizontal solid walls have been studied. The effective Rayleigh number was close to 106 . The Navier–Stokes equations, written with the Boussinesq approximation, were solved using the fractional-step method. The instant and time-averaged flow fields, the pulsation spectra and the integral heat transfer data were analyzed. The in-house code SINF/Flag-S results were compared with the available experimental data, and with the data obtained using the commercial software ANSYS Fluent 15.0. Copyright © 2017, St. Petersburg Polytechnic University. Production and hosting by Elsevier B.V. This is an open access article under the CC BY-NC-ND license. (http://creativecommons.org/licenses/by-nc-nd/4.0/) Keywords: Rayleigh–Bénard convection; Rotating container; Liquid metal; Conjugate heat transfer.
Introduction Free convection with low Prandtl numbers is widespread in nature and technology (convection in the outer core of the Earth, steelmaking, semiconductor crystal growth from the melt, heat exchange processes in nuclear reactors, and so on). One of the model problems in this field is Rayleigh–Bénard convection in a vertically oriented cylindrical vessel (container) filled with low-Prandtl-number fluid (liquid
✩ Peer review under responsibility of St. Petersburg Polytechnic University. ∗ Corresponding author. E-mail addresses: [email protected] (S.I. Smirnov), [email protected] (E.M. Smirnov), [email protected] (A.A. Smirnovsky).
metal or semiconductor melt). In numerical simulation, this problem is typically considered in a nonconjugate setting where boundary conditions are imposed directly on the inner surfaces of the solid walls confining the container. It is most often assumed that the horizontal walls are isothermal and that the lateral cylindrical surface is adiabatic. Specifying the temperature directly at the inner surfaces of the confining walls is a good approximation to reality only if the thermal resistance of the walls is much less than the resistance of the fluid medium. However, in cases where Prandtl numbers are low, the thermal resistances of the fluid and solid regions may turn out to be comparable, which requires taking into account the conjugate heat transfer effects. This point is well illustrated by the results obtained in the experimental study of free convection in a
http://dx.doi.org/10.1016/j.spjpm.2017.05.009 2405-7223/Copyright © 2017, St. Petersburg Polytechnic University. Production and hosting by Elsevier B.V. This is an open access article under the CC BY-NC-ND license. (http://creativecommons.org/licenses/by-nc-nd/4.0/) (Peer review under responsibility of St. Petersburg Polytechnic University).
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
JID: SPJPM
2
ARTICLE IN PRESS
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
non-rotating cylinder [1], according to which varying the ratio of the physical parameters of the fluid and solid media has a significant impact on the flow structure and on the heat transfer. To date, numerous experimental computational studies have been carried on turbulent Rayleigh– Bénard convection developing in regions with different geometric shapes and in fluids with different Prandtl numbers (see, for example, the review paper [2] containing an extensive reference list). Among these studies, various experimental works were dedicated to turbulent fluid convection in cylindrical containers with small height-to-diameter ratios and for fluids with Prandtl numbers exceeding unity (these numbers are typical for water, for example). Funfschilling et al. [3,4] presented the results of the recent extensive studies on water convection in non-rotating containers. Water convection in rotating containers was studied in Refs. [5–7]. Considerably fewer studies were dedicated to the consideration of Rayleigh–Bénard convection at low Prandtl numbers. Turbulent convection of mercury in a non-rotating cylindrical container was experimentally studied in Refs. [8,9]. The results of the studies on liquid metal convection in rotating containers filled by mercury have been reported in Refs. [10–12]; the first two of these papers dealt with turbulent Rayleigh–Bénard convection in cavities whose diameter exceeded their height by an order of magnitude, and the third one presented the findings for the case of a rotating cylinder with the height equal to its diameter. The authors of all of the above-mentioned experimental studies have strived to simulate the conditions allowing to neglect the heat transport effects in the walls confining the container. The effect of thermal resistance of the walls on the convection structure and the integral heat transfer has been studied in relatively few papers. Some of the notable experimental works considering the effect of conjugate heat exchange on convection in a cylindrical (non-rotating) container are the detailed studies [1,13], carried out for the media with different Prandtl numbers: Pr = 0.7 [1] and Pr = 4.4 [13]. During the last three decades, turbulent Rayleigh– Bénard convection has been extensively studied by Direct Numerical Simulation (DNS), covering a wide range of Prandtl numbers. A horizontal layer with the periodicity conditions imposed on the vertical boundaries is often considered as the computational domain. The results of detailed computational studies on convection in a stationary horizontal layer at Prandtl num-
bers about unity are given in [14–17], for the case of a rotating layer. Convection in a stationary layer at low Prandtl numbers was studied by the DNS method in Ref. [14] for Pr = 0.07 with the Rayleigh number Ra ranging from 104 to 107 , as well as in Ref. [18] for Pr = 0.025 and Ra = 105 . The effect of superimposed global rotation was studied in Ref. [19] for Pr = 0.1 and 104 ≤ Ra ≤ 108 . The effects of conjugate heat exchange on the structure of turbulent convection and the heat transfer in a Non-rotating horizontal layer, confined by walls of finite thickness, were studied by the DNS method in Ref. [20] for the case of Pr = 0.025, Ra = 105 . No such investigations have been carried out thus far for the case of a rotating layer. A substantial amount of computational data obtained by the DNS method for turbulent Rayleigh– Bénard convection developing in cylindrical containers (in the non-conjugate problem setting) filled with fluids with significantly varying Prandtl numbers has been accumulated by now (see, for example, Ref. [21], citing numerous experimental and theoretical studies in this field). The convection in stationary cylindrical containers filled with fluids with the Prandtl number varying from 0.1 to 104 and the Rayleigh number varying from 105 to 109 was considered in Ref. [22]. The results of the DNS computations of turbulent convection in a rotating cylinder with the Prandtl number about unity and 108 ≤ Ra ≤ 109 were presented in Refs. [21,23]. The case of Rayleigh–Bénard convection in a cylindrical (non-rotating) container filled with fluid with a low Prandtl number was studied by the DNS method in Refs. [22,24–26]. Verzicco et al. [24,25] numerically simulated the convection of mercury (Pr = 0.025) for Rayleigh numbers up to 106 . The results of numerical simulations at higher Rayleigh numbers (up to 109 ) were presented in Ref. [22] for Pr = 0.1 and in Ref. [26] for Pr = 0.021. Abramov et al. [27] compared the data obtained in the conjugate and non-conjugate settings for the case of mercury convection at Ra = 106 . We could find no simulation studies on turbulent convection of fluids with low Prandtl numbers under rotation of a cylindrical container in the literature available to us. This paper presents the results of direct numerical simulation of turbulent free convection of mercury in a rotating cylinder heated from below, in the conjugate and non-conjugate problem settings. The main simulations were performed using software code of our own development. The results of these simulations have been compared with the experimental data [12] as well as with the results of the
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
ARTICLE IN PRESS
JID: SPJPM
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
3
Fig. 1. Computational domain with the problem parameters (a) and the cross-sections of a typical computational mesh by the vertical (b) and the horizontal (c) planes. The solid disks (cross-hatched) and the fluid region are shown; Th and Tc are the constant temperature values.
computations performed with the ANSYS Fluent 15.0 software package [28]. Problem setting Direct numerical simulation of natural mercury (Pr = 0.025) convection in a rotating container heated from below was performed in two settings: in the conjugate one taking into account the heat transfer effects in the horizontal walls of finite thickness, and the nonconjugate one corresponding to the idealized case with the horizontal walls of zero thickness. The computational domain for the conjugate problem (Fig. 1a) includes a cylindrical cell filled with fluid, with the diameter D and the height H (H/D = 1) and two solid-state disks with the thickness h. This paper presents the computations for the case h/D = 0.25. Impermeability and no-slip conditions are imposed on all boundaries of the inner cylindrical cell filled by the fluid. Constant temperature values (Th , Tc ) are specified at the outer horizontal surfaces of the solidstate disks as the thermal boundary conditions. The inner surfaces of the disks are regarded as the surfaces of thermal coupling between the fluid and the solid. The lateral cylindrical surfaces of the inner cell and the disks are assumed to be adiabatic. The thermal convective flow of the fluid with constant physical properties and heat transfer in a rotating container with walls of finite thickness are described by system of Eqs. (1)–(3) including the continuity
equation, unsteady Navier–Stokes equations, written in the Boussinesq approximation to account for buoyancy effects in the gravitational and the centrifugal force fields, and the unsteady equation of convective– diffusive heat transport (diffusive-only in the solid regions). The system has the following form: ∇ ·V=0
(1)
∂V 1 + (V · ∇ )V = − ∇ p + β(T0 − T ) g + ω2 R ∂t ρf −2ω × V + ν∇ 2 V
(2)
∂T λi 2 + ( V · ∇ )T = ∇ T, (3) ∂t ρiCi where V = (Vx , Vy , Vz ) is the velocity field; t is the time, p is the pressure, T is the temperature, ρ i is the medium density (i = f or s), ν is the kinematic viscosity, β is the thermal expansion coefficient, Ci is the medium heat capacity, λi is the thermal conductivity coefficient, g is the gravitational acceleration, T0 is the temperature under hydrostatic equilibrium, ω is the angular velocity of rotation, R is the distance from the current point of the computational domain to the rotation axis. The energy equation is solved simultaneously for the fluid (i = f) and the solid (i = s) subdomains. From this point on, we will exclude the buoyancy effects in the centrifugal force field (the term β(T0 – T)ω2 R in the equation of motion) from consideration, assuming them to be much less than the buoyancy effects in the gravitational force field.
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
ARTICLE IN PRESS
JID: SPJPM
4
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
In addition to the Prandtl number Pr = μCf / λf , two other dimensionless parameters governing unsteady convection can be constructed from the set of the physical properties of the media considered in solving the conjugate problem, namely, λf /λs and ρ f Cf /ρ s Cs . In this investigation, it was assumed that λf /λs = 1.641, ρ f Cf / ρ s Cs = 0.472; this corresponds to the ratios for the initial dimensional parameters, taken for mercury as the fluid and for steel as the wall material. Scale similarity criteria (the Rayleigh number and the rotation parameter, constructed from the scale quantities) specified in the course of the simulations, were determined for the conjugate and the nonconjugate settings as follows:
method (FSM) to advance in physical time. The FSM is widely used (in various variants) to solve numerically unsteady hydrodynamic problems (see, for example, Refs. [29–31]). The principle of the FSM is in separating the spatial operators forming the momentum equation, and in interpreting the role of the pressure gradient as a projection operator that converts an arbitrary velocity field to a solenoidal one. As in Ref. [29], the below-described method is based on the Crank– Nicolson scheme (which has a second-order accuracy in time) for the momentum equation, approximated as follows: Vn+1 − Vn 1 = − ∇ pn+1/2 + β T0 − T n+1/2 g t ρf
Ra = Pr · (gβT0 H 3 /ν 2 );
(4)
K = 2ω · (H/(gβT0 ))0.5 ,
(5)
(8)
where T0 is the temperature difference (T0 = Th – Tc ). In the conjugate setting, the convection in the cavity is actually determined by the values of the effective Rayleigh number Raeff [27] and the effective rotation parameter Keff (in the case of the rotating cavity), which can be constructed from the difference T between the space- and time-averaged temperatures at the fluid-solid interface. Thus, (6) Rae f f = Pr · gβT H 3 /ν 2 ;
where t is the time step; n, n + 1 are the time layers. We should note that the convective term [∇ · (VV )]n+1/2 in Eq. (8) can be computed by different methods. The most common one is extrapolation of the convective term as an ensemble from two preceding time levels by the Adams–Bashforth scheme [29]:
Ke f f = 2ω · (H/(gβT ) ) . 0.5
(7)
The effective values of the governing parameters are computed after the end of the computations carried out in the conjugate setting. To extract the actual conjugate heat exchange effects from the obtained simulation data, the scale similarity criteria are selected (iteratively) such that the effective Rayleigh number is almost equal to some reference value chosen when solving the non-conjugate problem. In this investigation, the reference value of the Rayleigh number was approximately equal to 106 . Numerical method The finite-volume ‘unstructured’ code SINF/Flag-S, developed by the team of the Fluid Dynamics, Combustion and Heat Transfer department of Peter the Great St. Petersburg Polytechnic University, was used for the main series of computations. As one of the options for solving Navier–Stokes equations, this code uses an original version of the implicit fractional-step
− 2ω × Vn+1/2 +ν∇ 2 Vn+1/2 −[∇ · (VV )]n+1/2 ,
[∇ · (VV )]n+1/2 = 1.5[∇ · (VV )]n − 0.5[∇ · (VV )]n−1 . (9) Since this scheme is explicit (with respect to the convective terms), the computational algorithm as a whole is stable only for Courant numbers less than unity. As noted in Ref. [31], using Eq. (9) for flows with high Reynolds numbers may require the introduction of a stabilizing term. Additionally, the experience of using this scheme in the SINF/Flag-S code revealed that for highly skewed cells, the Courant number should be significantly reduced for computational stability. In view of this, the implicit setting is preferable for computing the convective term. In particular, Ref. [30] also suggested, with respect to the Crank– Nicolson scheme, to compute this term by the semiimplicit scheme: (10) [∇ · (VV )]n+1/2 = ∇ · Vn Vn+1 . According to Ref. [30], scheme (10) provides computational stability at Courant numbers greater than unity while involving only two time levels. It is known, however, that the Crank–Nicolson scheme is neutrally stable and can lead to non-physical oscillations with time in a number of problems with strong nonlinear (convective) effects. A common way of
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
ARTICLE IN PRESS
JID: SPJPM
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
5
suppressing these oscillations is adopting a scheme weighing the contributions from the first-order implicit Euler scheme and the Crank–Nicolson scheme. This, however, leads to a reduction in the order of time accuracy. To settle these issues, the SINF/Flag-S code uses an original modification for computing the convective terms, which combines extrapolation from two preceding time levels by the Adams–Bashforth scheme with implicitness introduced into the scheme as follows:
manner. Eq. (8) is split into two equations:
[∇ · (VV )]n+1/2 = ∇ · (Vn+1/2 Vn+1/2 ),
Eq. (16) of the system is written relative to the ‘predictor’ velocity V∗ , which is found using the pressure field on the preceding time level pn –1/2 . The second equation relates the sought-for velocity and pressure in the new time level. The closing relation is the solenoidal condition imposed on the field of the sought-for velocity:
(11)
and the velocity in the intermediate (n + 1/2) level computed as Vn+1/2 = 0.5 Vn+1 + Vn . (12) The quantity V n+1/2 is regarded as the result of linear extrapolation from two lower time levels: Vn+1/2 = 1.5V − 0.5V n
n−1
.
(13)
The advantage of this approach for computing the convective term over the one proposed in Ref. [30] is that if the above-mentioned drawback of the Crank– Nicolson scheme emerges, it is possible to switch directly to approximating the equation of conservation with three-level discretization of the time derivative by the Euler ‘backward difference’ scheme (in this case the extrapolation is performed to the (n + 1) level instead of the intermediate (n + 1/2) one). This scheme does not generate non-physical oscillations in time and has a second order of accuracy. Since non-physical oscillations have not emerged in carrying out the turbulent convection computations (shown below) by the Crank–Nicolson scheme, the computations were performed by this scheme with the convective term in the equation of motion approximated by formula (11). The term in the right-hand side of this equation, representing the effect of the buoyancy force, includes the quantity T n+1/2 = 0.5 T n+1 + T n , (14) which is determined when solving the energy equation, also approximated by the Crank–Nicolson scheme: T
−T λi 2 n+1/2 = ∇ T − ∇ · (Vn+1/2 T n+1/2 ). t ρiCi (15)
n+1
n
The algorithm for computing the velocity field satisfying the equation of motion and the continuity equation in a new time level is constructed in the following
V∗ − Vn 1 = − ∇ pn−1/2 + β T0 − T n+1/2 g t ρf − ω × V∗ + Vn + 0.5ν∇ 2 V∗ + Vn − 0.5∇ · Vn+1/2 V∗ + Vn (16) Vn+1 − V∗ 1 = − ∇ pn+1/2 − pn−1/2 . t ρf
∇ · Vn+1 = 0.
(17)
(18)
Combining Eqs. (17) and (18) produces the Poisson equation for the pressure increment p = pn +1/ 2 – pn –1/2 in the form t ∇ 2 (p) = −ρ f ∇ · V∗ .
(19)
The sought-for velocity field is calculated from solving Eq. (19): Vn+1 = V∗ − t
∇ (p) . ρf
(20)
The solution of Eqs. (15), (16) and (19) in the SINF/Flag-S code is sought iteratively, introducing residuals and corrections; for example, the ‘predictor’ velocity is sought as ∗ V∗ → Vm+1 = Vm∗ + δV∗ ,
(21)
where m is the number of the previous iteration. The velocity Vn in the preceding time level is taken as the first approximation (m = 0) for the velocity V∗ . Eq. (16) is converted to the form ∗ n δV∗ ˜ s Vm∗ , ... + Vm − V ≡ −R, + L δV∗ = − R t t (22) ˜ s (Vm∗ , ... ) is the where L is the ‘stabilizing’ operator, R residual of the right-hand part of Eq. (16), R is the full residual of the momentum equation. Eqs. (15) and (19) are transformed similarly (however, the time derivative is absent in the latter equation). Writing the solved equations in the form (22) allows to use different stencils (sets of computational points) for approximating the spatial operators in the
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
ARTICLE IN PRESS
JID: SPJPM
6
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
left- and right-hand sides. The stabilizing operator in the SINF/Flag-S code is generally approximated with the most compact stencil, which, in particular, makes computing parallelization by the Domain Decomposition technology considerably easier, especially in case of unstructured grids; the convective terms in the operator L are approximated by the first order upwind scheme. The stencil is extended when computing the residual in the right-hand side of Eq. (22), which allows, on the one hand, to use high-order schemes for the convective terms (the QUICK scheme that nominally has the third order of accuracy was used in our computations), and, on the other hand, to take into account the skewness of the computational cells. Similar remarks apply to the method used for approximating Eqs. (15) and (19). Since the SINF/Flag-S code is based on using collocated grids, the Rhie–Chow correction [32], computed by the pressure field in the preceding time level, is introduced in order to suppress the development of ‘even-odd’ spatial oscillations in the pressure and velocity fields, similarly to the computations performed in Ref. [33]. Computational aspects All computations were carried out on the Polytechnic RSC Tornado cluster of the Peter the Great St. Petersburg Polytechnic University’s supercomputer center, with peak performance of about a thousand PFLOPS. As noted above, for solving the conjugate problem, the effective Rayleigh number had to be maintained equal to the reference Rayleigh number corresponding to the non-conjugate case. To do this, several simulations were performed in the conjugate setting for each selected value of the scale rotational parameter. The Raeff values were computed after the computations had been completed. The obtained Raeff values were then interpolated to the point corresponding to the reference value Ra = 106 (selected for the non-conjugate problem), and the scale Rayleigh number that had to be specified in the conjugate setting was evaluated to provide the effective Rayleigh number equal to 106 . The final stage of the simulation was performed with the estimated value of the scale Rayleigh number for the scale rotational parameter that had not been varied in the course of the above-described procedure of selecting Ra. Some of the results obtained by using the SINF/Flag-S code were compared with the data of the computations carried out with the commercial soft-
ware ANSYS Fluent 15.0, where the convective terms were also discretized using the QUICK scheme (the diffusive terms were approximated by a second-order central difference scheme). The time advance was also performed by the fractional step method, with three subiterations specified for each equation (the “Max. Corrections” parameter was equal to 3). The rest of the parameters of the ANSYS Fluent 15.0 solver were default. We used a computational grid consisting of hexagonal elements: about 5·105 cells belonged to the fluid region and 2·105 cells to each of the horizontal walls. The grid was clustered to the walls and to the interface surfaces (the vertical cell dimension near the interface amounted to about 1.5·10–4 H). The distribution of the cells in various cross-sections is shown in Fig. 1b,c. The time step was of the order of one thousandth of the characteristic convection time of the problem, defined as tconv = (H/(gβT ))0.5 .
(23)
The local Courant number did not exceed unity. The samples used for averaging amounted to 1000– 3000 dimensionless times. Computing a process lasting 10 dimensionless times (about 20,000 time steps) took about three hours on twenty computing cores for both codes. The velocity components and the temperature in the points located within the fluid region, near the adiabatic wall in the central horizontal plane and near the interface between the media, were monitored during the computations. The plane-averaged heat flux on the outer walls and at the interfaces was also monitored.
Computational results As noted above, the scale quantities were chosen so as to ensure the effective Rayleigh number close to 106 . In order to maintain continuity with the previous study [27], the reference value of the Rayleigh number, selected for solving the non-conjugate problem, was 9.64·105 . In solving the conjugate problem, the scale Rayleigh numbers providing Raeff = 9.64·105 were selected for three variants with different values of the specified scale rotational parameter, namely, for K = 0; 0.499 and 5.823. The selected values of the scale Rayleigh number were, respectively, 2.72·106 , 2.66·106 and 1.28·106 . The resulting values of the effective rotational parameter Keff amounted to 0, 0.834 and 6.706.
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
JID: SPJPM
ARTICLE IN PRESS
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
7
Fig. 2. Isosurfaces of the vertical velocity component at three values of the effective rotational parameter: Keff =0 (a); 0.834 (b); 6.706 (c). The dark structures correspond to the downward flow; the light structures correspond to the upward flow.
Fig. 2 shows the convective flow structure, simulated in the conjugate setting with three different rotational intensities. At low values of the rotational parameter (Keff = 0 and 0.834) the flow can be generally characterized as a global convective cell (Fig. 2a, b) whose visual images were earlier presented in [21,23,27]. The flow structure changes with increasing rotation intensity: the convective cell disintegrates into a group of smaller structures (Fig. 2c). Similar changes in the flow structure with increasing rotational parameter were noted by the authors of Ref. [6], who performed the DNS computations with Pr = 0.7, and in Ref. [23] with Pr = 6.4. We should also note that the vertical velocity isosurfaces obtained for the same values of the governing parameters in the non-conjugate computations are similar to those shown in Fig. 2. Fig. 3 shows the pattern of the temporal changes in the dimensionless vertical velocity component, computed for different rotation rates of the container in the conjugate setting (the buoyancy rate Vb =(HgβT)0.5 serves as the velocity scale). The time sample (which was also used for averaging) is much longer in the absence of rotation than for the variants with the rotating container. This is due to the fact that the abovementioned global circulation cell oscillates in space from time to time. However, it can be seen from Fig. 3a that two relatively stable states are observed for this structure. Similar behavior of the numerical solution is also observed in simulations using the AN-
SYS Fluent 15.0 software package. We should also note that similar phenomena have been observed in the experiments described in Ref. [9], where they were explained by the high sensitivity of the global convection cell to small defects of the experimental setup and by non-ideal isothermality of the horizontal walls. With applied rotation of the container, the oscillation of the global convection cell ceases to be random in character. The action of the Coriolis force under moderate rotation leads to a precession of the convection cell with a constant angular velocity; this can be clearly observed in Fig. 3b, which shows the periodic alternation of the sections corresponding to the positive and negative values of the vertical velocity. Strong rotation results not only in the disintegration of the global convection cell into several structures, but also in a sharp decrease in the fluctuation amplitude of the vertical velocity component (see Fig. 3c), which indicates the general suppression of convection. The computational results in a non-conjugate setting allow to come to similar conclusions for the cases of moderate and strong rotation. Fig. 4 shows the dependences of the dimensionless RMS (Root Mean Square) fluctuations of the velocity and temperature components on the dimensionless vertical coordinate (measured from the lower interface, along the axis of the container) for different variants of computing the convection in a rotating cavity. Here u, v, w are the fluctuations of the velocity components Vx ,
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
JID: SPJPM
8
ARTICLE IN PRESS
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
Fig. 3. Dynamics of the vertical velocity component near the adiabatic wall in the central horizontal cross-section at three values of the effective rotational parameter: Keff = 0 (a), 0.834 (b), 6.706 (c). The computations were performed using the ANSYS Fluent 15.0 software package (gray curves) and the SINF/Flag-S code (black curves).
Fig. 4. Vertical distributions of various functions according to the results of the computations in the conjugate (curves 1с, 2с) and non-conjugate (1n, 2n) problem settings for two values of the effective rotational parameter: Keff = 0.834 (1с, 1n) and 6.706 (2с, 2n). The functions shown in the figure: RMS fluctuations u, w (a) and v (b) of the Vx , Vz (a) and Vy (b) velocity components, respectively; RMS fluctuations of the temperature Т(c); mean temperature profiles (d).
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
ARTICLE IN PRESS
JID: SPJPM
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
9
Fig. 5. Energy spectra of the vertical velocity component according to the computational results in the conjugate setting at three values of the effective rotational parameter: Keff = 0 (a), 0.834 (b), 6.706 (c). The dotted lines are the curves of the spectral density decreasing by the Kolmogorov ‘–5/3’ law.
Vy , Vz , respectively. The profiles of the mean dimensionless temperature are also shown. We should stress that in solving the conjugate problem the temperature was counted from the averaged value of the temperature of the lower interface and transformed into a dimensionless quantity by dividing this value by the difference T between the averaged temperatures of the lower and upper interfaces. It can be seen from the data shown in Fig. 4a–c that in case of moderate rotation (Keff = 0.834), the computations in the conjugate setting predict slightly higher fluctuation intensities for the velocity and the temperature than those obtained for the non-conjugate problem: the differences in the vicinity of the central horizontal cross-section are about 10%. Strong rotation (Keff = 6.706) leads to a sharp decrease in the fluctuation level, and the difference in the profiles of the RMS values obtained for the conjugate and the non-conjugate settings practically disappears. Fig. 4d shows that the transition to the conjugate setting has no pronounced effect on the behavior of the mean temperature profile in the region occupied by the fluid both under moderate and relatively strong rotation. In the latter case, the mean temperature profile becomes linear, pointing to the prevalence of diffusive heat transfer over the convective one for the entire fluid region. The effect of rotation rate on the energy spectrum of the vertical velocity component is illustrated by Fig. 5 (the monitoring point is located in the central horizontal plane, at a distance of 0.05Н from the adiabatic wall). The dimensionless oscillation frequency f = f · tconv ,
(24)
where f’ is the dimensional frequency, is plotted on the horizontal axis. The energy spectra obtained for the cases of zero and moderate rotation rates (Fig. 5a, b), are sufficiently filled and indicate that the simulated convective flow is turbulent. The dashed lines show the dependences corresponding to a decrease in the spectral density by the Kolmogorov ‘–5/3’ law. At the same time, a spectral energy density peak can be observed (see Fig. 5b for Keff = 0.834); it corresponds to the precession frequency of the global convection cell, the Strouhal number for this precession is Sh = fpeak = 0.016. With increasing rotation rate, the turbulence is suppressed and the energy spectrum changes drastically (Fig. 5c). Moving on to discussing the results characterizing the integral heat transfer, let us first consider the data of the verification computations carried out in the nonconjugate setting for the case of the non-rotating cavity. Table 1 contains a comparison of the space- and time-averaged Nusselt numbers obtained in the course of these computations with the numerical and experimental data available in literature for the Rayleigh number equal to 106 . To bring the data obtained in this study with Ra = 9.64·105 to the Ra = 106 , we used the relation Nu1 /Nu2 = (Ra1 /Ra2 )n ,
(25)
where Nu1 , Nu2 were the Nusselt numbers at Ra1 =106 and Ra2 = 9.64·105 , and the exponent n was taken to equal 0.25, in accordance with the correlation dependence established, in particular, in the experimental study in Ref. [12]. Table 1 shows that all the computed data are mutually consistent and in good agreement with the ex-
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
ARTICLE IN PRESS
JID: SPJPM
10
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
Table 1 Comparison of the computed (DNS) and experimental Nusselt numbers for the case of non-conjugate problem setting and a nonrotating cavity. DNS computation
Experiment
Parameter Nu , %
[m3+dc;June 22, 2017;11:9]
This study
[27]
[24]
[25]
[8]
[9]
[12]
5.58 0
5.66 1.4
5.59 0.2
5.55 –0.5
6.46 14
5.08 –9.8
5.85 4.6
Notes. The computed data obtained by direct numerical simulation (DNS) in the present study at Ra = 9.64·105 , brought to the case Ra = 106 (K = 0) using the relation Nu1 /Nu2 =(Ra1 /Ra2 )n , where Nu1 , Nu2 are the Nusselt numbers for Rayleigh numbers Ra1 =106 and Ra2 = 9.64·105 ; n = 0.25. The Nu values were space- and time-averaged; the value characterizes the percentage deviation from the value obtained by us.
Table 2 Effect of the rotational parameter and heat transfer in the horizontal walls of the cavity on the Nusselt number values, computed for different cases of the problem setting and with two different codes. Computational code
Problem setting
Ra
K
Keff
Nu
SINF/Flag-S SINF/Flag-S SINF/Flag-S SINF/Flag-S SINF/Flag-S SINF/Flag-S ANSYS Fluent 15.0 ANSYS Fluent 15.0
Conjugate Conjugate Conjugate Non-conjugate Non-conjugate Non-conjugate Conjugate
2.72·106 1.28·106 2.66·106 9.64·105 9.64·105 9.64·105 2.72·106
0 5.823 0.499 0 6.706 0.834 0
0 6.706 0.834 0 6.706 0.834 0
5.975 1.064 5.845 5.527 1.062 5.320 5.980
Conjugate
1.28·106
5.823
6.706
1.070
Note. The computational results for the Nusselt number, averaged over surface and time, were obtained with Raeff = 9.64·105 .
perimental results of Ref. [12]. Earlier experimental studies [8,9] provide the Nusselt number values deviating from the DNS results by 10–15% in the opposite directions. The effects of the rotational parameter and the heat transport in the horizontal walls on the integral heat transfer in thermal convective flow can be determined by analyzing values of the space- and time-averaged Nusselt number Nu , presented in Table 2. It can be seen that the values of this number obtained for the conjugate problem under zero and moderate rotation of the container are 8–10% higher than in the case of the non-conjugate problem (such an estimate was earlier obtained for a non-rotating cylinder in Ref. [27]). This agrees with the above-noted increase in the fluctuation intensity of velocity and temperature observed in solving the conjugate problem. As a result of a significant increase in the rotation rate (Keff = 6.706), when convection is suppressed, the mean Nusselt number Nu falls to near unity for both the conjugate and the non-conjugate problems.
The computations using the ANSYS Fluent 15.0 software package were also performed for the variants corresponding to the effective rotational parameter values Keff = 0 and 6.706. In both cases, the obtained results agree very well with the data of the computations by the SINF/Flag-S code. Fig. 6 illustrates the temporal variation of the surface-averaged heat flux (Nusselt number Nu) at the interface. Notably, in case of the conjugate problem, the fluctuation amplitude of the number Nu, computed for the variant with moderate rotation of the container (Fig. 6a), is markedly less than the amplitude of its fluctuations on the isothermal wall in the case of the non-conjugate problem. A possible explanation for this is that temperature fluctuations are capable of permeating the solid region through the interfacial surfaces. Due to the general suppression of convection, the fluctuation intensity of the surface-averaged heat flux at
Fig. 6. Time fluctuations of the surface-averaged Nusselt number at the interface at two values of the effective rotational parameter: Keff = 0 (a) and 6.706 (b). The black curves are the results of the computations in the conjugate setting; the gray curves are the results of the computations in the non-conjugate setting.
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
ARTICLE IN PRESS
JID: SPJPM
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
the interface is also strongly reduced under strong rotation (Fig. 6b). Conclusions We have performed direct numerical simulation (DNS) of turbulent mercury convection in a rotating and a non-rotating cylindrical container heated from below, with the height equal to the diameter, with zero and finite conductivity of the solid horizontal walls. The computational data presented was obtained with the effective Rayleigh number close to 106 . We have established that the flow that develops in the container under small and moderate values of the rotational parameter can be generally classified as a global convection cell for the cases of both the conjugate and the non-conjugate settings of the heat transfer problem. In the absence of rotation, the cell oscillates from time to time. Precession of the convection cell with a constant angular velocity emerges in a rotating cavity. The flow pattern changes considerably with a further increase in the rotation rate: the convection cell disintegrates into a group of smaller structures, while the fluctuation amplitude of the vertical velocity component decreases sharply. Up to the values of the effective rotational parameter close to unity, the integral heat transfer through the fluid layer under comparable conditions (the equality of the effective Rayleigh number) can be higher by 10% if the heat transport in the solid horizontal walls is taken into account. With strong rotation practically suppressing convection, the heat transfer is mainly governed by the diffusion mechanism in the case of both the conjugate and the non-conjugate heat exchange problems. References [1] R. Verzicco, Effects of nonperfect thermal sources in turbulent thermal convection, Phys. Fluids 16 (6) (2004) 1965–1979. [2] G. Ahlers, S. Grossmann, D. Lohse, Heat transfer and large scale dynamics in turbulent Rayleigh–Bénard convection, Rev. Mod. Phys. 81 (2) (2009) 503–537. [3] D. Funfschilling, E. Brown, A. Nikolaenko, G. Ahlers, Heat transport by turbulent Rayleigh–Bénard convection in cylindrical samples with aspect ratio one and larger, J. Fluid Mech. 536 (2005) 145–154. [4] S. Weiss, G. Ahlers, Turbulent Rayleigh–Bénard convection in a cylindrical container with aspect ratio =0.50 and Prandtl number Pr = 4.38, J. Fluid Mech. 676 (2011) 5–40. [5] J.-Q. Zhong, G. Ahlers, Heat transport and the large-scale circulation in rotating turbulent Rayleigh–Bénard convection, J. Fluid Mech. 665 (2010) 300–333.
11
[6] J.-Q. Zhong, R.J.A.M. Stevens, H.J.H. Clercx, et al., Prandtl-, Rayleigh-, and Rossby-Number dependence of heat transport in turbulent rotating Rayleigh–Bénard convection, Phys. Rev. Lett. 102 (4) (2009) 044502. [7] Y. Liu, R.E. Ecke, Heat transport measurements in turbulent rotating Rayleigh–Bénard convection, Phys. Rev. E 80 (3) (2009) 036314. [8] T. Takeshita, T. Segawa, J.A. Glazier, M. Sano, Thermal turbulence in mercury, Phys. Rev. Lett. 76 (9) (1996) 1465–1468. [9] S. Cioni, S. Ciliberto, J. Sommeria, Strongly turbulent Rayleigh–Bénard convection in mercury: comparison with results at moderate Prandtl number, J. Fluid Mech. 335 (1997) 111–140. [10] H.T. Rossby, A study of Bénard convection with and without rotation, J. Fluid Mech. 36 (2) (1969) 309–335. [11] J.M. Aurnou, P.L. Olson, Experiments on Rayleigh–Bénard convection, magnetoconvection and rotating magnetoconvection in liquid gallium, J. Fluid Mech. 430 (2001) 283–307. [12] E.M. King, J.M. Aurnou, Turbulent convection in liquid metal with and without rotation, PNAS USA 110 (17) (2013) 6688–6693. [13] E. Brown, A. Nikolaenko, D. Funfschilling, G. Ahlers, Heat transport in turbulent Rayleigh–Bénard convection: effect of finite top- and bottom-plate conductivities, Phys. Fluids 17 (2005) 075108. [14] R.M. Kerr, J.R. Herring, Prandtl number dependence of Nusselt number in direct numerical simulations, J. Fluid Mech. 419 (1) (2000) 325–344. [15] T. Hartlep, A. Tilgner, F.H. Busse, Transition to turbulent convection in a fluid layer heated from below at moderate aspect ratio, J. Fluid Mech. 544 (2005) 309–322. [16] K. Julien, S. Legg, J. McWilliams, J. Werne, Rapidly rotating turbulent Rayleigh–Bénard convection, J. Fluid Mech. 322 (1996) 243–273. [17] R.P.J. Kunnen, B.J. Geurts, H.J.H. Clercx, Direct numerical simulation of turbulent rotating Rayleigh–Benard convection, in: Proceedings of the Sixth International ERCOFTAC Workshop on Direct and Large-Eddy Simulation VI, 2006, pp. 233–240. [18] I. Oti´c, G. Grötzbach, Direct numerical simulation and RANS modeling of turbulent natural convection for low Prandtl number fluids, in: Proceedings of the American Society of Mechanical Engineers, 2, 2004, pp. 159–165. [19] H.K. Pharasi, R. Kannan, K. Kumar, J.K. Bhattacharjee, Turbulence in rotating Rayleigh–Bénard convection in low-Prandtl-number fluids, Phys. Rev. E 84 (4) (2011) 047301. [20] N. Ivanov, S. Smirnov, Numerical simulation of turbulent conjugate free convection in horizontal low-Pr fluid layer, in: Proceedings of the Seventh Baltic Heat Transfer Conference, 2015, pp. 29–34. [21] R.J.A.M. Stevens, H.J.H. Clercx, D. Lohse, Heat transport and flow structure in rotating Rayleigh–Bénard convection, Eur. J. Mech. B Fluids 40 (2013) 41–49. [22] G. Silano, K.R. Screenivasan, R. Verzicco, Numerical simulations of Rayleigh–Bénard convection for Prandtl numbers between 10−1 and 104 and Rayleigh numbers between 105 and 109 , J. Fluid Mech. 662 (2010) 409–446. [23] G.L. Kooij, M.A. Botchev, B.J. Geurts, Direct numerical simulation of Nusselt number scaling in rotating Rayleigh–Bénard convection, Int. J. Heat Fluid Flow 55 (2015) 26–33.
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
JID: SPJPM
12
ARTICLE IN PRESS
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
[24] R. Verzicco, R. Camussi, Transitional regimes of low-Prandtl thermal convection in a cylindrical cell, Phys. Fluids 9 (5) (1997) 1287–1295. [25] A. Abramov, A. Korsakov, Direct numerical modeling of mercury turbulent convection in axisymmetric reservoirs including magnetic field effects, Heat Transf. Res. 35 (1–2) (2004) 76–84. [26] J.D. Scheel, J. Schumacher, Global and local statistics in turbulent convection at low Prandtl numbers, J. Fluid Mech. 802 (2016) 147–173. [27] A. Abramov, E. Smirnov, A. Smirnovsky, Numerical simulation of turbulent Rayleigh–Bénard conjugate convection of low Pr fluid in a cylindrical container, in: Proceedings of the Seventh Baltic Heat Transfer Conference, 2015, pp. 11–16. [28] ANSYS Fluent 15.0, User Guide (2013). [29] J. Kim, P. Moin, Application of a fractional-step method to incompressible Navier–Stokes equations, J. Comput. Phys. 59 (1985) 308–323.
[30] Y.-J. Jan, T.W.-H. Sheu, A quasi-implicit time advancing scheme for unsteady incompressible flow, Part I: validation, Comput. Methods Appl. Mech. Eng. 196 (45–48) (2007) 4755–4770. [31] R. Bevan, E. Boileau, R. Loon, R. Lewis, R. Nithiarasu, A comparative study of fractional step method in its quasi-implicit, semi-implicit and fully-explicit forms for incompressible flows, Int. J. Numer. Methods Heat Fluid Flow 26 (2016) 595–623. [32] C.M. Rhie, W.L. Chow, Numerical study of the turbulent flow past an airfoil with trailing edge separation, AIAA J. 21 (11) (1983) 1525–1532. [33] S.W. Armeld, N. Williamson, M.P. Kirkpatrick, R. Street, A divergence free fractional-step method for the Navier–Stokes equations on non-staggered grids, ANZIAM J. 51 (2010) C654–C667 (EMAC2009).
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
JID: SPJPM
[m3+dc;June 22, 2017;11:9]
Available online at www.sciencedirect.com
St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12 www.elsevier.com/locate/spjpm
Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder✩ Sergey I. Smirnov∗, Evgueni M. Smirnov, Alexander A. Smirnovsky Peter the Great St. Petersburg Polytechnic University, 29 Politekhnicheskaya St., St. Petersburg 195251, Russian Federation Available online xxx
Abstract The obtained results of direct numerical simulation of the free mercury convection in a rotating cylindrical container heated from below are presented. Setting the Prandtl number equal to 0.025 and the height-to-diameter ratio equal to 1.0, effects of container rotation and heat transfer in horizontal solid walls have been studied. The effective Rayleigh number was close to 106 . The Navier–Stokes equations, written with the Boussinesq approximation, were solved using the fractional-step method. The instant and time-averaged flow fields, the pulsation spectra and the integral heat transfer data were analyzed. The in-house code SINF/Flag-S results were compared with the available experimental data, and with the data obtained using the commercial software ANSYS Fluent 15.0. Copyright © 2017, St. Petersburg Polytechnic University. Production and hosting by Elsevier B.V. This is an open access article under the CC BY-NC-ND license. (http://creativecommons.org/licenses/by-nc-nd/4.0/) Keywords: Rayleigh–Bénard convection; Rotating container; Liquid metal; Conjugate heat transfer.
Introduction Free convection with low Prandtl numbers is widespread in nature and technology (convection in the outer core of the Earth, steelmaking, semiconductor crystal growth from the melt, heat exchange processes in nuclear reactors, and so on). One of the model problems in this field is Rayleigh–Bénard convection in a vertically oriented cylindrical vessel (container) filled with low-Prandtl-number fluid (liquid
✩ Peer review under responsibility of St. Petersburg Polytechnic University. ∗ Corresponding author. E-mail addresses: [email protected] (S.I. Smirnov), [email protected] (E.M. Smirnov), [email protected] (A.A. Smirnovsky).
metal or semiconductor melt). In numerical simulation, this problem is typically considered in a nonconjugate setting where boundary conditions are imposed directly on the inner surfaces of the solid walls confining the container. It is most often assumed that the horizontal walls are isothermal and that the lateral cylindrical surface is adiabatic. Specifying the temperature directly at the inner surfaces of the confining walls is a good approximation to reality only if the thermal resistance of the walls is much less than the resistance of the fluid medium. However, in cases where Prandtl numbers are low, the thermal resistances of the fluid and solid regions may turn out to be comparable, which requires taking into account the conjugate heat transfer effects. This point is well illustrated by the results obtained in the experimental study of free convection in a
http://dx.doi.org/10.1016/j.spjpm.2017.05.009 2405-7223/Copyright © 2017, St. Petersburg Polytechnic University. Production and hosting by Elsevier B.V. This is an open access article under the CC BY-NC-ND license. (http://creativecommons.org/licenses/by-nc-nd/4.0/) (Peer review under responsibility of St. Petersburg Polytechnic University).
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
JID: SPJPM
2
ARTICLE IN PRESS
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
non-rotating cylinder [1], according to which varying the ratio of the physical parameters of the fluid and solid media has a significant impact on the flow structure and on the heat transfer. To date, numerous experimental computational studies have been carried on turbulent Rayleigh– Bénard convection developing in regions with different geometric shapes and in fluids with different Prandtl numbers (see, for example, the review paper [2] containing an extensive reference list). Among these studies, various experimental works were dedicated to turbulent fluid convection in cylindrical containers with small height-to-diameter ratios and for fluids with Prandtl numbers exceeding unity (these numbers are typical for water, for example). Funfschilling et al. [3,4] presented the results of the recent extensive studies on water convection in non-rotating containers. Water convection in rotating containers was studied in Refs. [5–7]. Considerably fewer studies were dedicated to the consideration of Rayleigh–Bénard convection at low Prandtl numbers. Turbulent convection of mercury in a non-rotating cylindrical container was experimentally studied in Refs. [8,9]. The results of the studies on liquid metal convection in rotating containers filled by mercury have been reported in Refs. [10–12]; the first two of these papers dealt with turbulent Rayleigh–Bénard convection in cavities whose diameter exceeded their height by an order of magnitude, and the third one presented the findings for the case of a rotating cylinder with the height equal to its diameter. The authors of all of the above-mentioned experimental studies have strived to simulate the conditions allowing to neglect the heat transport effects in the walls confining the container. The effect of thermal resistance of the walls on the convection structure and the integral heat transfer has been studied in relatively few papers. Some of the notable experimental works considering the effect of conjugate heat exchange on convection in a cylindrical (non-rotating) container are the detailed studies [1,13], carried out for the media with different Prandtl numbers: Pr = 0.7 [1] and Pr = 4.4 [13]. During the last three decades, turbulent Rayleigh– Bénard convection has been extensively studied by Direct Numerical Simulation (DNS), covering a wide range of Prandtl numbers. A horizontal layer with the periodicity conditions imposed on the vertical boundaries is often considered as the computational domain. The results of detailed computational studies on convection in a stationary horizontal layer at Prandtl num-
bers about unity are given in [14–17], for the case of a rotating layer. Convection in a stationary layer at low Prandtl numbers was studied by the DNS method in Ref. [14] for Pr = 0.07 with the Rayleigh number Ra ranging from 104 to 107 , as well as in Ref. [18] for Pr = 0.025 and Ra = 105 . The effect of superimposed global rotation was studied in Ref. [19] for Pr = 0.1 and 104 ≤ Ra ≤ 108 . The effects of conjugate heat exchange on the structure of turbulent convection and the heat transfer in a Non-rotating horizontal layer, confined by walls of finite thickness, were studied by the DNS method in Ref. [20] for the case of Pr = 0.025, Ra = 105 . No such investigations have been carried out thus far for the case of a rotating layer. A substantial amount of computational data obtained by the DNS method for turbulent Rayleigh– Bénard convection developing in cylindrical containers (in the non-conjugate problem setting) filled with fluids with significantly varying Prandtl numbers has been accumulated by now (see, for example, Ref. [21], citing numerous experimental and theoretical studies in this field). The convection in stationary cylindrical containers filled with fluids with the Prandtl number varying from 0.1 to 104 and the Rayleigh number varying from 105 to 109 was considered in Ref. [22]. The results of the DNS computations of turbulent convection in a rotating cylinder with the Prandtl number about unity and 108 ≤ Ra ≤ 109 were presented in Refs. [21,23]. The case of Rayleigh–Bénard convection in a cylindrical (non-rotating) container filled with fluid with a low Prandtl number was studied by the DNS method in Refs. [22,24–26]. Verzicco et al. [24,25] numerically simulated the convection of mercury (Pr = 0.025) for Rayleigh numbers up to 106 . The results of numerical simulations at higher Rayleigh numbers (up to 109 ) were presented in Ref. [22] for Pr = 0.1 and in Ref. [26] for Pr = 0.021. Abramov et al. [27] compared the data obtained in the conjugate and non-conjugate settings for the case of mercury convection at Ra = 106 . We could find no simulation studies on turbulent convection of fluids with low Prandtl numbers under rotation of a cylindrical container in the literature available to us. This paper presents the results of direct numerical simulation of turbulent free convection of mercury in a rotating cylinder heated from below, in the conjugate and non-conjugate problem settings. The main simulations were performed using software code of our own development. The results of these simulations have been compared with the experimental data [12] as well as with the results of the
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
ARTICLE IN PRESS
JID: SPJPM
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
3
Fig. 1. Computational domain with the problem parameters (a) and the cross-sections of a typical computational mesh by the vertical (b) and the horizontal (c) planes. The solid disks (cross-hatched) and the fluid region are shown; Th and Tc are the constant temperature values.
computations performed with the ANSYS Fluent 15.0 software package [28]. Problem setting Direct numerical simulation of natural mercury (Pr = 0.025) convection in a rotating container heated from below was performed in two settings: in the conjugate one taking into account the heat transfer effects in the horizontal walls of finite thickness, and the nonconjugate one corresponding to the idealized case with the horizontal walls of zero thickness. The computational domain for the conjugate problem (Fig. 1a) includes a cylindrical cell filled with fluid, with the diameter D and the height H (H/D = 1) and two solid-state disks with the thickness h. This paper presents the computations for the case h/D = 0.25. Impermeability and no-slip conditions are imposed on all boundaries of the inner cylindrical cell filled by the fluid. Constant temperature values (Th , Tc ) are specified at the outer horizontal surfaces of the solidstate disks as the thermal boundary conditions. The inner surfaces of the disks are regarded as the surfaces of thermal coupling between the fluid and the solid. The lateral cylindrical surfaces of the inner cell and the disks are assumed to be adiabatic. The thermal convective flow of the fluid with constant physical properties and heat transfer in a rotating container with walls of finite thickness are described by system of Eqs. (1)–(3) including the continuity
equation, unsteady Navier–Stokes equations, written in the Boussinesq approximation to account for buoyancy effects in the gravitational and the centrifugal force fields, and the unsteady equation of convective– diffusive heat transport (diffusive-only in the solid regions). The system has the following form: ∇ ·V=0
(1)
∂V 1 + (V · ∇ )V = − ∇ p + β(T0 − T ) g + ω2 R ∂t ρf −2ω × V + ν∇ 2 V
(2)
∂T λi 2 + ( V · ∇ )T = ∇ T, (3) ∂t ρiCi where V = (Vx , Vy , Vz ) is the velocity field; t is the time, p is the pressure, T is the temperature, ρ i is the medium density (i = f or s), ν is the kinematic viscosity, β is the thermal expansion coefficient, Ci is the medium heat capacity, λi is the thermal conductivity coefficient, g is the gravitational acceleration, T0 is the temperature under hydrostatic equilibrium, ω is the angular velocity of rotation, R is the distance from the current point of the computational domain to the rotation axis. The energy equation is solved simultaneously for the fluid (i = f) and the solid (i = s) subdomains. From this point on, we will exclude the buoyancy effects in the centrifugal force field (the term β(T0 – T)ω2 R in the equation of motion) from consideration, assuming them to be much less than the buoyancy effects in the gravitational force field.
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
ARTICLE IN PRESS
JID: SPJPM
4
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
In addition to the Prandtl number Pr = μCf / λf , two other dimensionless parameters governing unsteady convection can be constructed from the set of the physical properties of the media considered in solving the conjugate problem, namely, λf /λs and ρ f Cf /ρ s Cs . In this investigation, it was assumed that λf /λs = 1.641, ρ f Cf / ρ s Cs = 0.472; this corresponds to the ratios for the initial dimensional parameters, taken for mercury as the fluid and for steel as the wall material. Scale similarity criteria (the Rayleigh number and the rotation parameter, constructed from the scale quantities) specified in the course of the simulations, were determined for the conjugate and the nonconjugate settings as follows:
method (FSM) to advance in physical time. The FSM is widely used (in various variants) to solve numerically unsteady hydrodynamic problems (see, for example, Refs. [29–31]). The principle of the FSM is in separating the spatial operators forming the momentum equation, and in interpreting the role of the pressure gradient as a projection operator that converts an arbitrary velocity field to a solenoidal one. As in Ref. [29], the below-described method is based on the Crank– Nicolson scheme (which has a second-order accuracy in time) for the momentum equation, approximated as follows: Vn+1 − Vn 1 = − ∇ pn+1/2 + β T0 − T n+1/2 g t ρf
Ra = Pr · (gβT0 H 3 /ν 2 );
(4)
K = 2ω · (H/(gβT0 ))0.5 ,
(5)
(8)
where T0 is the temperature difference (T0 = Th – Tc ). In the conjugate setting, the convection in the cavity is actually determined by the values of the effective Rayleigh number Raeff [27] and the effective rotation parameter Keff (in the case of the rotating cavity), which can be constructed from the difference T between the space- and time-averaged temperatures at the fluid-solid interface. Thus, (6) Rae f f = Pr · gβT H 3 /ν 2 ;
where t is the time step; n, n + 1 are the time layers. We should note that the convective term [∇ · (VV )]n+1/2 in Eq. (8) can be computed by different methods. The most common one is extrapolation of the convective term as an ensemble from two preceding time levels by the Adams–Bashforth scheme [29]:
Ke f f = 2ω · (H/(gβT ) ) . 0.5
(7)
The effective values of the governing parameters are computed after the end of the computations carried out in the conjugate setting. To extract the actual conjugate heat exchange effects from the obtained simulation data, the scale similarity criteria are selected (iteratively) such that the effective Rayleigh number is almost equal to some reference value chosen when solving the non-conjugate problem. In this investigation, the reference value of the Rayleigh number was approximately equal to 106 . Numerical method The finite-volume ‘unstructured’ code SINF/Flag-S, developed by the team of the Fluid Dynamics, Combustion and Heat Transfer department of Peter the Great St. Petersburg Polytechnic University, was used for the main series of computations. As one of the options for solving Navier–Stokes equations, this code uses an original version of the implicit fractional-step
− 2ω × Vn+1/2 +ν∇ 2 Vn+1/2 −[∇ · (VV )]n+1/2 ,
[∇ · (VV )]n+1/2 = 1.5[∇ · (VV )]n − 0.5[∇ · (VV )]n−1 . (9) Since this scheme is explicit (with respect to the convective terms), the computational algorithm as a whole is stable only for Courant numbers less than unity. As noted in Ref. [31], using Eq. (9) for flows with high Reynolds numbers may require the introduction of a stabilizing term. Additionally, the experience of using this scheme in the SINF/Flag-S code revealed that for highly skewed cells, the Courant number should be significantly reduced for computational stability. In view of this, the implicit setting is preferable for computing the convective term. In particular, Ref. [30] also suggested, with respect to the Crank– Nicolson scheme, to compute this term by the semiimplicit scheme: (10) [∇ · (VV )]n+1/2 = ∇ · Vn Vn+1 . According to Ref. [30], scheme (10) provides computational stability at Courant numbers greater than unity while involving only two time levels. It is known, however, that the Crank–Nicolson scheme is neutrally stable and can lead to non-physical oscillations with time in a number of problems with strong nonlinear (convective) effects. A common way of
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
ARTICLE IN PRESS
JID: SPJPM
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
5
suppressing these oscillations is adopting a scheme weighing the contributions from the first-order implicit Euler scheme and the Crank–Nicolson scheme. This, however, leads to a reduction in the order of time accuracy. To settle these issues, the SINF/Flag-S code uses an original modification for computing the convective terms, which combines extrapolation from two preceding time levels by the Adams–Bashforth scheme with implicitness introduced into the scheme as follows:
manner. Eq. (8) is split into two equations:
[∇ · (VV )]n+1/2 = ∇ · (Vn+1/2 Vn+1/2 ),
Eq. (16) of the system is written relative to the ‘predictor’ velocity V∗ , which is found using the pressure field on the preceding time level pn –1/2 . The second equation relates the sought-for velocity and pressure in the new time level. The closing relation is the solenoidal condition imposed on the field of the sought-for velocity:
(11)
and the velocity in the intermediate (n + 1/2) level computed as Vn+1/2 = 0.5 Vn+1 + Vn . (12) The quantity V n+1/2 is regarded as the result of linear extrapolation from two lower time levels: Vn+1/2 = 1.5V − 0.5V n
n−1
.
(13)
The advantage of this approach for computing the convective term over the one proposed in Ref. [30] is that if the above-mentioned drawback of the Crank– Nicolson scheme emerges, it is possible to switch directly to approximating the equation of conservation with three-level discretization of the time derivative by the Euler ‘backward difference’ scheme (in this case the extrapolation is performed to the (n + 1) level instead of the intermediate (n + 1/2) one). This scheme does not generate non-physical oscillations in time and has a second order of accuracy. Since non-physical oscillations have not emerged in carrying out the turbulent convection computations (shown below) by the Crank–Nicolson scheme, the computations were performed by this scheme with the convective term in the equation of motion approximated by formula (11). The term in the right-hand side of this equation, representing the effect of the buoyancy force, includes the quantity T n+1/2 = 0.5 T n+1 + T n , (14) which is determined when solving the energy equation, also approximated by the Crank–Nicolson scheme: T
−T λi 2 n+1/2 = ∇ T − ∇ · (Vn+1/2 T n+1/2 ). t ρiCi (15)
n+1
n
The algorithm for computing the velocity field satisfying the equation of motion and the continuity equation in a new time level is constructed in the following
V∗ − Vn 1 = − ∇ pn−1/2 + β T0 − T n+1/2 g t ρf − ω × V∗ + Vn + 0.5ν∇ 2 V∗ + Vn − 0.5∇ · Vn+1/2 V∗ + Vn (16) Vn+1 − V∗ 1 = − ∇ pn+1/2 − pn−1/2 . t ρf
∇ · Vn+1 = 0.
(17)
(18)
Combining Eqs. (17) and (18) produces the Poisson equation for the pressure increment p = pn +1/ 2 – pn –1/2 in the form t ∇ 2 (p) = −ρ f ∇ · V∗ .
(19)
The sought-for velocity field is calculated from solving Eq. (19): Vn+1 = V∗ − t
∇ (p) . ρf
(20)
The solution of Eqs. (15), (16) and (19) in the SINF/Flag-S code is sought iteratively, introducing residuals and corrections; for example, the ‘predictor’ velocity is sought as ∗ V∗ → Vm+1 = Vm∗ + δV∗ ,
(21)
where m is the number of the previous iteration. The velocity Vn in the preceding time level is taken as the first approximation (m = 0) for the velocity V∗ . Eq. (16) is converted to the form ∗ n δV∗ ˜ s Vm∗ , ... + Vm − V ≡ −R, + L δV∗ = − R t t (22) ˜ s (Vm∗ , ... ) is the where L is the ‘stabilizing’ operator, R residual of the right-hand part of Eq. (16), R is the full residual of the momentum equation. Eqs. (15) and (19) are transformed similarly (however, the time derivative is absent in the latter equation). Writing the solved equations in the form (22) allows to use different stencils (sets of computational points) for approximating the spatial operators in the
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
ARTICLE IN PRESS
JID: SPJPM
6
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
left- and right-hand sides. The stabilizing operator in the SINF/Flag-S code is generally approximated with the most compact stencil, which, in particular, makes computing parallelization by the Domain Decomposition technology considerably easier, especially in case of unstructured grids; the convective terms in the operator L are approximated by the first order upwind scheme. The stencil is extended when computing the residual in the right-hand side of Eq. (22), which allows, on the one hand, to use high-order schemes for the convective terms (the QUICK scheme that nominally has the third order of accuracy was used in our computations), and, on the other hand, to take into account the skewness of the computational cells. Similar remarks apply to the method used for approximating Eqs. (15) and (19). Since the SINF/Flag-S code is based on using collocated grids, the Rhie–Chow correction [32], computed by the pressure field in the preceding time level, is introduced in order to suppress the development of ‘even-odd’ spatial oscillations in the pressure and velocity fields, similarly to the computations performed in Ref. [33]. Computational aspects All computations were carried out on the Polytechnic RSC Tornado cluster of the Peter the Great St. Petersburg Polytechnic University’s supercomputer center, with peak performance of about a thousand PFLOPS. As noted above, for solving the conjugate problem, the effective Rayleigh number had to be maintained equal to the reference Rayleigh number corresponding to the non-conjugate case. To do this, several simulations were performed in the conjugate setting for each selected value of the scale rotational parameter. The Raeff values were computed after the computations had been completed. The obtained Raeff values were then interpolated to the point corresponding to the reference value Ra = 106 (selected for the non-conjugate problem), and the scale Rayleigh number that had to be specified in the conjugate setting was evaluated to provide the effective Rayleigh number equal to 106 . The final stage of the simulation was performed with the estimated value of the scale Rayleigh number for the scale rotational parameter that had not been varied in the course of the above-described procedure of selecting Ra. Some of the results obtained by using the SINF/Flag-S code were compared with the data of the computations carried out with the commercial soft-
ware ANSYS Fluent 15.0, where the convective terms were also discretized using the QUICK scheme (the diffusive terms were approximated by a second-order central difference scheme). The time advance was also performed by the fractional step method, with three subiterations specified for each equation (the “Max. Corrections” parameter was equal to 3). The rest of the parameters of the ANSYS Fluent 15.0 solver were default. We used a computational grid consisting of hexagonal elements: about 5·105 cells belonged to the fluid region and 2·105 cells to each of the horizontal walls. The grid was clustered to the walls and to the interface surfaces (the vertical cell dimension near the interface amounted to about 1.5·10–4 H). The distribution of the cells in various cross-sections is shown in Fig. 1b,c. The time step was of the order of one thousandth of the characteristic convection time of the problem, defined as tconv = (H/(gβT ))0.5 .
(23)
The local Courant number did not exceed unity. The samples used for averaging amounted to 1000– 3000 dimensionless times. Computing a process lasting 10 dimensionless times (about 20,000 time steps) took about three hours on twenty computing cores for both codes. The velocity components and the temperature in the points located within the fluid region, near the adiabatic wall in the central horizontal plane and near the interface between the media, were monitored during the computations. The plane-averaged heat flux on the outer walls and at the interfaces was also monitored.
Computational results As noted above, the scale quantities were chosen so as to ensure the effective Rayleigh number close to 106 . In order to maintain continuity with the previous study [27], the reference value of the Rayleigh number, selected for solving the non-conjugate problem, was 9.64·105 . In solving the conjugate problem, the scale Rayleigh numbers providing Raeff = 9.64·105 were selected for three variants with different values of the specified scale rotational parameter, namely, for K = 0; 0.499 and 5.823. The selected values of the scale Rayleigh number were, respectively, 2.72·106 , 2.66·106 and 1.28·106 . The resulting values of the effective rotational parameter Keff amounted to 0, 0.834 and 6.706.
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
JID: SPJPM
ARTICLE IN PRESS
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
7
Fig. 2. Isosurfaces of the vertical velocity component at three values of the effective rotational parameter: Keff =0 (a); 0.834 (b); 6.706 (c). The dark structures correspond to the downward flow; the light structures correspond to the upward flow.
Fig. 2 shows the convective flow structure, simulated in the conjugate setting with three different rotational intensities. At low values of the rotational parameter (Keff = 0 and 0.834) the flow can be generally characterized as a global convective cell (Fig. 2a, b) whose visual images were earlier presented in [21,23,27]. The flow structure changes with increasing rotation intensity: the convective cell disintegrates into a group of smaller structures (Fig. 2c). Similar changes in the flow structure with increasing rotational parameter were noted by the authors of Ref. [6], who performed the DNS computations with Pr = 0.7, and in Ref. [23] with Pr = 6.4. We should also note that the vertical velocity isosurfaces obtained for the same values of the governing parameters in the non-conjugate computations are similar to those shown in Fig. 2. Fig. 3 shows the pattern of the temporal changes in the dimensionless vertical velocity component, computed for different rotation rates of the container in the conjugate setting (the buoyancy rate Vb =(HgβT)0.5 serves as the velocity scale). The time sample (which was also used for averaging) is much longer in the absence of rotation than for the variants with the rotating container. This is due to the fact that the abovementioned global circulation cell oscillates in space from time to time. However, it can be seen from Fig. 3a that two relatively stable states are observed for this structure. Similar behavior of the numerical solution is also observed in simulations using the AN-
SYS Fluent 15.0 software package. We should also note that similar phenomena have been observed in the experiments described in Ref. [9], where they were explained by the high sensitivity of the global convection cell to small defects of the experimental setup and by non-ideal isothermality of the horizontal walls. With applied rotation of the container, the oscillation of the global convection cell ceases to be random in character. The action of the Coriolis force under moderate rotation leads to a precession of the convection cell with a constant angular velocity; this can be clearly observed in Fig. 3b, which shows the periodic alternation of the sections corresponding to the positive and negative values of the vertical velocity. Strong rotation results not only in the disintegration of the global convection cell into several structures, but also in a sharp decrease in the fluctuation amplitude of the vertical velocity component (see Fig. 3c), which indicates the general suppression of convection. The computational results in a non-conjugate setting allow to come to similar conclusions for the cases of moderate and strong rotation. Fig. 4 shows the dependences of the dimensionless RMS (Root Mean Square) fluctuations of the velocity and temperature components on the dimensionless vertical coordinate (measured from the lower interface, along the axis of the container) for different variants of computing the convection in a rotating cavity. Here u, v, w are the fluctuations of the velocity components Vx ,
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
JID: SPJPM
8
ARTICLE IN PRESS
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
Fig. 3. Dynamics of the vertical velocity component near the adiabatic wall in the central horizontal cross-section at three values of the effective rotational parameter: Keff = 0 (a), 0.834 (b), 6.706 (c). The computations were performed using the ANSYS Fluent 15.0 software package (gray curves) and the SINF/Flag-S code (black curves).
Fig. 4. Vertical distributions of various functions according to the results of the computations in the conjugate (curves 1с, 2с) and non-conjugate (1n, 2n) problem settings for two values of the effective rotational parameter: Keff = 0.834 (1с, 1n) and 6.706 (2с, 2n). The functions shown in the figure: RMS fluctuations u, w (a) and v (b) of the Vx , Vz (a) and Vy (b) velocity components, respectively; RMS fluctuations of the temperature Т(c); mean temperature profiles (d).
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
ARTICLE IN PRESS
JID: SPJPM
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
9
Fig. 5. Energy spectra of the vertical velocity component according to the computational results in the conjugate setting at three values of the effective rotational parameter: Keff = 0 (a), 0.834 (b), 6.706 (c). The dotted lines are the curves of the spectral density decreasing by the Kolmogorov ‘–5/3’ law.
Vy , Vz , respectively. The profiles of the mean dimensionless temperature are also shown. We should stress that in solving the conjugate problem the temperature was counted from the averaged value of the temperature of the lower interface and transformed into a dimensionless quantity by dividing this value by the difference T between the averaged temperatures of the lower and upper interfaces. It can be seen from the data shown in Fig. 4a–c that in case of moderate rotation (Keff = 0.834), the computations in the conjugate setting predict slightly higher fluctuation intensities for the velocity and the temperature than those obtained for the non-conjugate problem: the differences in the vicinity of the central horizontal cross-section are about 10%. Strong rotation (Keff = 6.706) leads to a sharp decrease in the fluctuation level, and the difference in the profiles of the RMS values obtained for the conjugate and the non-conjugate settings practically disappears. Fig. 4d shows that the transition to the conjugate setting has no pronounced effect on the behavior of the mean temperature profile in the region occupied by the fluid both under moderate and relatively strong rotation. In the latter case, the mean temperature profile becomes linear, pointing to the prevalence of diffusive heat transfer over the convective one for the entire fluid region. The effect of rotation rate on the energy spectrum of the vertical velocity component is illustrated by Fig. 5 (the monitoring point is located in the central horizontal plane, at a distance of 0.05Н from the adiabatic wall). The dimensionless oscillation frequency f = f · tconv ,
(24)
where f’ is the dimensional frequency, is plotted on the horizontal axis. The energy spectra obtained for the cases of zero and moderate rotation rates (Fig. 5a, b), are sufficiently filled and indicate that the simulated convective flow is turbulent. The dashed lines show the dependences corresponding to a decrease in the spectral density by the Kolmogorov ‘–5/3’ law. At the same time, a spectral energy density peak can be observed (see Fig. 5b for Keff = 0.834); it corresponds to the precession frequency of the global convection cell, the Strouhal number for this precession is Sh = fpeak = 0.016. With increasing rotation rate, the turbulence is suppressed and the energy spectrum changes drastically (Fig. 5c). Moving on to discussing the results characterizing the integral heat transfer, let us first consider the data of the verification computations carried out in the nonconjugate setting for the case of the non-rotating cavity. Table 1 contains a comparison of the space- and time-averaged Nusselt numbers obtained in the course of these computations with the numerical and experimental data available in literature for the Rayleigh number equal to 106 . To bring the data obtained in this study with Ra = 9.64·105 to the Ra = 106 , we used the relation Nu1 /Nu2 = (Ra1 /Ra2 )n ,
(25)
where Nu1 , Nu2 were the Nusselt numbers at Ra1 =106 and Ra2 = 9.64·105 , and the exponent n was taken to equal 0.25, in accordance with the correlation dependence established, in particular, in the experimental study in Ref. [12]. Table 1 shows that all the computed data are mutually consistent and in good agreement with the ex-
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
ARTICLE IN PRESS
JID: SPJPM
10
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
Table 1 Comparison of the computed (DNS) and experimental Nusselt numbers for the case of non-conjugate problem setting and a nonrotating cavity. DNS computation
Experiment
Parameter Nu , %
[m3+dc;June 22, 2017;11:9]
This study
[27]
[24]
[25]
[8]
[9]
[12]
5.58 0
5.66 1.4
5.59 0.2
5.55 –0.5
6.46 14
5.08 –9.8
5.85 4.6
Notes. The computed data obtained by direct numerical simulation (DNS) in the present study at Ra = 9.64·105 , brought to the case Ra = 106 (K = 0) using the relation Nu1 /Nu2 =(Ra1 /Ra2 )n , where Nu1 , Nu2 are the Nusselt numbers for Rayleigh numbers Ra1 =106 and Ra2 = 9.64·105 ; n = 0.25. The Nu values were space- and time-averaged; the value characterizes the percentage deviation from the value obtained by us.
Table 2 Effect of the rotational parameter and heat transfer in the horizontal walls of the cavity on the Nusselt number values, computed for different cases of the problem setting and with two different codes. Computational code
Problem setting
Ra
K
Keff
Nu
SINF/Flag-S SINF/Flag-S SINF/Flag-S SINF/Flag-S SINF/Flag-S SINF/Flag-S ANSYS Fluent 15.0 ANSYS Fluent 15.0
Conjugate Conjugate Conjugate Non-conjugate Non-conjugate Non-conjugate Conjugate
2.72·106 1.28·106 2.66·106 9.64·105 9.64·105 9.64·105 2.72·106
0 5.823 0.499 0 6.706 0.834 0
0 6.706 0.834 0 6.706 0.834 0
5.975 1.064 5.845 5.527 1.062 5.320 5.980
Conjugate
1.28·106
5.823
6.706
1.070
Note. The computational results for the Nusselt number, averaged over surface and time, were obtained with Raeff = 9.64·105 .
perimental results of Ref. [12]. Earlier experimental studies [8,9] provide the Nusselt number values deviating from the DNS results by 10–15% in the opposite directions. The effects of the rotational parameter and the heat transport in the horizontal walls on the integral heat transfer in thermal convective flow can be determined by analyzing values of the space- and time-averaged Nusselt number Nu , presented in Table 2. It can be seen that the values of this number obtained for the conjugate problem under zero and moderate rotation of the container are 8–10% higher than in the case of the non-conjugate problem (such an estimate was earlier obtained for a non-rotating cylinder in Ref. [27]). This agrees with the above-noted increase in the fluctuation intensity of velocity and temperature observed in solving the conjugate problem. As a result of a significant increase in the rotation rate (Keff = 6.706), when convection is suppressed, the mean Nusselt number Nu falls to near unity for both the conjugate and the non-conjugate problems.
The computations using the ANSYS Fluent 15.0 software package were also performed for the variants corresponding to the effective rotational parameter values Keff = 0 and 6.706. In both cases, the obtained results agree very well with the data of the computations by the SINF/Flag-S code. Fig. 6 illustrates the temporal variation of the surface-averaged heat flux (Nusselt number Nu) at the interface. Notably, in case of the conjugate problem, the fluctuation amplitude of the number Nu, computed for the variant with moderate rotation of the container (Fig. 6a), is markedly less than the amplitude of its fluctuations on the isothermal wall in the case of the non-conjugate problem. A possible explanation for this is that temperature fluctuations are capable of permeating the solid region through the interfacial surfaces. Due to the general suppression of convection, the fluctuation intensity of the surface-averaged heat flux at
Fig. 6. Time fluctuations of the surface-averaged Nusselt number at the interface at two values of the effective rotational parameter: Keff = 0 (a) and 6.706 (b). The black curves are the results of the computations in the conjugate setting; the gray curves are the results of the computations in the non-conjugate setting.
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
ARTICLE IN PRESS
JID: SPJPM
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
the interface is also strongly reduced under strong rotation (Fig. 6b). Conclusions We have performed direct numerical simulation (DNS) of turbulent mercury convection in a rotating and a non-rotating cylindrical container heated from below, with the height equal to the diameter, with zero and finite conductivity of the solid horizontal walls. The computational data presented was obtained with the effective Rayleigh number close to 106 . We have established that the flow that develops in the container under small and moderate values of the rotational parameter can be generally classified as a global convection cell for the cases of both the conjugate and the non-conjugate settings of the heat transfer problem. In the absence of rotation, the cell oscillates from time to time. Precession of the convection cell with a constant angular velocity emerges in a rotating cavity. The flow pattern changes considerably with a further increase in the rotation rate: the convection cell disintegrates into a group of smaller structures, while the fluctuation amplitude of the vertical velocity component decreases sharply. Up to the values of the effective rotational parameter close to unity, the integral heat transfer through the fluid layer under comparable conditions (the equality of the effective Rayleigh number) can be higher by 10% if the heat transport in the solid horizontal walls is taken into account. With strong rotation practically suppressing convection, the heat transfer is mainly governed by the diffusion mechanism in the case of both the conjugate and the non-conjugate heat exchange problems. References [1] R. Verzicco, Effects of nonperfect thermal sources in turbulent thermal convection, Phys. Fluids 16 (6) (2004) 1965–1979. [2] G. Ahlers, S. Grossmann, D. Lohse, Heat transfer and large scale dynamics in turbulent Rayleigh–Bénard convection, Rev. Mod. Phys. 81 (2) (2009) 503–537. [3] D. Funfschilling, E. Brown, A. Nikolaenko, G. Ahlers, Heat transport by turbulent Rayleigh–Bénard convection in cylindrical samples with aspect ratio one and larger, J. Fluid Mech. 536 (2005) 145–154. [4] S. Weiss, G. Ahlers, Turbulent Rayleigh–Bénard convection in a cylindrical container with aspect ratio =0.50 and Prandtl number Pr = 4.38, J. Fluid Mech. 676 (2011) 5–40. [5] J.-Q. Zhong, G. Ahlers, Heat transport and the large-scale circulation in rotating turbulent Rayleigh–Bénard convection, J. Fluid Mech. 665 (2010) 300–333.
11
[6] J.-Q. Zhong, R.J.A.M. Stevens, H.J.H. Clercx, et al., Prandtl-, Rayleigh-, and Rossby-Number dependence of heat transport in turbulent rotating Rayleigh–Bénard convection, Phys. Rev. Lett. 102 (4) (2009) 044502. [7] Y. Liu, R.E. Ecke, Heat transport measurements in turbulent rotating Rayleigh–Bénard convection, Phys. Rev. E 80 (3) (2009) 036314. [8] T. Takeshita, T. Segawa, J.A. Glazier, M. Sano, Thermal turbulence in mercury, Phys. Rev. Lett. 76 (9) (1996) 1465–1468. [9] S. Cioni, S. Ciliberto, J. Sommeria, Strongly turbulent Rayleigh–Bénard convection in mercury: comparison with results at moderate Prandtl number, J. Fluid Mech. 335 (1997) 111–140. [10] H.T. Rossby, A study of Bénard convection with and without rotation, J. Fluid Mech. 36 (2) (1969) 309–335. [11] J.M. Aurnou, P.L. Olson, Experiments on Rayleigh–Bénard convection, magnetoconvection and rotating magnetoconvection in liquid gallium, J. Fluid Mech. 430 (2001) 283–307. [12] E.M. King, J.M. Aurnou, Turbulent convection in liquid metal with and without rotation, PNAS USA 110 (17) (2013) 6688–6693. [13] E. Brown, A. Nikolaenko, D. Funfschilling, G. Ahlers, Heat transport in turbulent Rayleigh–Bénard convection: effect of finite top- and bottom-plate conductivities, Phys. Fluids 17 (2005) 075108. [14] R.M. Kerr, J.R. Herring, Prandtl number dependence of Nusselt number in direct numerical simulations, J. Fluid Mech. 419 (1) (2000) 325–344. [15] T. Hartlep, A. Tilgner, F.H. Busse, Transition to turbulent convection in a fluid layer heated from below at moderate aspect ratio, J. Fluid Mech. 544 (2005) 309–322. [16] K. Julien, S. Legg, J. McWilliams, J. Werne, Rapidly rotating turbulent Rayleigh–Bénard convection, J. Fluid Mech. 322 (1996) 243–273. [17] R.P.J. Kunnen, B.J. Geurts, H.J.H. Clercx, Direct numerical simulation of turbulent rotating Rayleigh–Benard convection, in: Proceedings of the Sixth International ERCOFTAC Workshop on Direct and Large-Eddy Simulation VI, 2006, pp. 233–240. [18] I. Oti´c, G. Grötzbach, Direct numerical simulation and RANS modeling of turbulent natural convection for low Prandtl number fluids, in: Proceedings of the American Society of Mechanical Engineers, 2, 2004, pp. 159–165. [19] H.K. Pharasi, R. Kannan, K. Kumar, J.K. Bhattacharjee, Turbulence in rotating Rayleigh–Bénard convection in low-Prandtl-number fluids, Phys. Rev. E 84 (4) (2011) 047301. [20] N. Ivanov, S. Smirnov, Numerical simulation of turbulent conjugate free convection in horizontal low-Pr fluid layer, in: Proceedings of the Seventh Baltic Heat Transfer Conference, 2015, pp. 29–34. [21] R.J.A.M. Stevens, H.J.H. Clercx, D. Lohse, Heat transport and flow structure in rotating Rayleigh–Bénard convection, Eur. J. Mech. B Fluids 40 (2013) 41–49. [22] G. Silano, K.R. Screenivasan, R. Verzicco, Numerical simulations of Rayleigh–Bénard convection for Prandtl numbers between 10−1 and 104 and Rayleigh numbers between 105 and 109 , J. Fluid Mech. 662 (2010) 409–446. [23] G.L. Kooij, M.A. Botchev, B.J. Geurts, Direct numerical simulation of Nusselt number scaling in rotating Rayleigh–Bénard convection, Int. J. Heat Fluid Flow 55 (2015) 26–33.
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009
JID: SPJPM
12
ARTICLE IN PRESS
[m3+dc;June 22, 2017;11:9]
S.I. Smirnov et al. / St. Petersburg Polytechnical University Journal: Physics and Mathematics 000 (2017) 1–12
[24] R. Verzicco, R. Camussi, Transitional regimes of low-Prandtl thermal convection in a cylindrical cell, Phys. Fluids 9 (5) (1997) 1287–1295. [25] A. Abramov, A. Korsakov, Direct numerical modeling of mercury turbulent convection in axisymmetric reservoirs including magnetic field effects, Heat Transf. Res. 35 (1–2) (2004) 76–84. [26] J.D. Scheel, J. Schumacher, Global and local statistics in turbulent convection at low Prandtl numbers, J. Fluid Mech. 802 (2016) 147–173. [27] A. Abramov, E. Smirnov, A. Smirnovsky, Numerical simulation of turbulent Rayleigh–Bénard conjugate convection of low Pr fluid in a cylindrical container, in: Proceedings of the Seventh Baltic Heat Transfer Conference, 2015, pp. 11–16. [28] ANSYS Fluent 15.0, User Guide (2013). [29] J. Kim, P. Moin, Application of a fractional-step method to incompressible Navier–Stokes equations, J. Comput. Phys. 59 (1985) 308–323.
[30] Y.-J. Jan, T.W.-H. Sheu, A quasi-implicit time advancing scheme for unsteady incompressible flow, Part I: validation, Comput. Methods Appl. Mech. Eng. 196 (45–48) (2007) 4755–4770. [31] R. Bevan, E. Boileau, R. Loon, R. Lewis, R. Nithiarasu, A comparative study of fractional step method in its quasi-implicit, semi-implicit and fully-explicit forms for incompressible flows, Int. J. Numer. Methods Heat Fluid Flow 26 (2016) 595–623. [32] C.M. Rhie, W.L. Chow, Numerical study of the turbulent flow past an airfoil with trailing edge separation, AIAA J. 21 (11) (1983) 1525–1532. [33] S.W. Armeld, N. Williamson, M.P. Kirkpatrick, R. Street, A divergence free fractional-step method for the Navier–Stokes equations on non-staggered grids, ANZIAM J. 51 (2010) C654–C667 (EMAC2009).
Please cite this article as: S.I. Smirnov et al., Endwall heat transfer effects on the turbulent mercury convection in a rotating cylinder, St. Petersburg Polytechnical University Journal: Physics and Mathematics (2017), http://dx.doi.org/10.1016/j.spjpm.2017.05.009