3D large eddy simulation (LES) calculations and experiments of natural convection in a laterally-heated cylindrical enclosure for crystal growth

International Journal of Thermal Sciences 116 (2017) 1e21

Contents lists available at ScienceDirect

International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts

3D large eddy simulation (LES) calculations and experiments of natural convection in a laterally-heated cylindrical enclosure for crystal growth Hooman Enayati, Abhilash J. Chandy*, Minel J. Braun, Nicholas Horning Department of Mechanical Engineering, The University of Akron, Akron, OH, 44325-3903, USA

a r t i c l e i n f o

a b s t r a c t

Article history: Received 10 October 2016 Received in revised form 26 December 2016 Accepted 26 January 2017

The ammonothermal crystal growth technique is one of the most effective methods used for growing Gallium Nitride (GaN) crystals in terms of the final product quality and manufacturing efficiency. Popular applications of GaN crystals include light emitting diodes (LEDs) and high-frequency electronic devices. A laterally-heated cylindrical enclosure with a top to bottom temperature gradient is considered in this study to investigate the fluid mechanics and heat transfer in an ammonothermal crystal growth reactor. Three-dimensional (3D) large eddy simulations (LES) results of natural convection in a laterally heated cylindrical reactor, using the commercial computational fluid dynamics (CFD) software, ANSYS FLUENT, are presented in this paper for a Rayleigh number (Ra) of 8:8
106 . The Ra is defined based on a length scale that is equal to the ratio of volume to the lateral area of the cylinder (R=2). In addition, an experiment of a geometrically- and dynamically-similar geometry is developed, and particle image velocimetry (PIV)-based flow visualizations are carried out for the purpose of validating the numerical model. Comparisons between experiments and numerical simulations showed that flow patterns were qualitatively similar, and Fourier transforms of velocity magnitudes at selected points in the domain matched reasonably well. An added interesting observation in the simulation was the existence of temperature inversion, which has potential implications on the choice of mineralizer (acidic/basic). © 2017 Elsevier Masson SAS. All rights reserved.

Keywords: Natural convection Lateral heating Vertical cylinder Low turbulent flow LES Three dimensional modeling Boussinesq approximation Crystal growth

1. Introduction Natural convection occurs due to the density differences resulting from concentration or temperature gradients. It is observed in several natural flows such as in the ocean and in the atmosphere, where a temperature gradient results in transport of nutrients and pollutants. The effectiveness of heat transfer due to natural convection also has widespread practical applications in energy storage, solar heating configurations, nuclear reactors and crystal growth industry. With regard to the crystal growth industry, Gallium nitride (GaN) crystals are manufactured via an ammonothermal process that utilizes natural convection in order to circulate the nutrients that create GaN in crystallized form in the reactor. The ammonothermal crystal growth technique is similar to the hydrothermal technique, that is utilized to grow the quartz crystals,

* Corresponding author. E-mail addresses: [email protected] (H. Enayati), [email protected] (A.J. Chandy). http://dx.doi.org/10.1016/j.ijthermalsci.2017.01.025 1290-0729/© 2017 Elsevier Masson SAS. All rights reserved.

with an important distinction. The difference is the type of the solvent in the reactor, which is ammonia for ammonothermal, and water for hydrothermal techniques [1e3]. The two distinct thermal zones in an ammonothermal crystal growth reactor would create a buoyant force in the reactor. The pre-loaded GaN crystals in the nutrient basket are dissolved in the dissolution zone and, by virtue of natural convection, they are transported to the crystallization zone, where they deposit on the seeds. In the ammonothermal method, as GaN crystals are almost insoluble in pure ammonia, mineralizers (acidic or basic) are added to the autoclave. Other alternative crystal growing methods are hydride vapor phase epitaxy (HVPE), solution growth, sodium flux, and the stoichiometric melting methods [4,5]. Natural convection has been studied both experimentally and numerically by several researchers for the past few decades [6e26]. Numerical studies included both 2D and 3D, and rectangular and cylindrical enclosures with different boundary conditions. The thermal configuration has a very significant role in determining the characteristics of temperature and flow distribution in the reactor. They primarily include side-wall heating, and top and bottom wall

2

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

heating with either constant temperature or heat flux. In addition, previous investigations have also looked at the effect of heat generation sources within the domain. The current study tries to introduce and analyze a relatively new thermal configuration, where the lower hot wall at a constant temperature in a cylindrical enclosure is separated from the upper colder wall, again at a constant temperature, by an insulated section, with no heat source. With regard to natural convection specifically in a cylindrical enclosure, there have been quite a few 2D and 3D numerical and also experimental studies with a variety of thermal configurations [27e31]. For instance, Chakir et al. [32] performed a 2D numerical calculation of turbulent natural convection of mixed gases in a horizontal annulus. Several simulations of thermal configurations involving a constant heat flux and constant temperature on the inner cylinder were carried out, while the outer cylinder was maintained at a constant temperature. The heat transfer and fluid motion were analyzed for Rayleigh number (Ra) values ranging from 105 to 1010 , and numerical results were compared to experiments. Barhaghi [33], studied the structure of turbulent natural convection boundary layers, and the effect of the buoyancy on a mixed convection boundary layer using direct numerical simulations (DNS) and large eddy simulations (LES). Flow and turbulence parameters in all cases were studied, and the results were compared to experiments as well. More recently, Malik et al. [34], performed a CFD analysis of conjugated heat transfer within a bottom-heated cylindrical enclosure. The study investigated the influence of the inner cylinder material and the outer cylinder geometric configurations on the heat transfer mechanism within the enclosure. Kang et al. [35], studied experimentally the effects of the Prandtl number and the curvature on the natural convection in a cylindrical geometry. In order to study the effect of curvature, they measured the Nusselt numbers for Ra of 1:4
109 - 3:7
1012 at Pr ¼ 2094. They showed a good agreement with laminar natural convection correlation data on the vertical plates. Several investigations have also been performed about the role of the natural convection in a crystal growth process [5,36e42]. Experimental and numerical studies in a rectangular cuboid were carried out by Li et al. [37]. The upper half of the sidewalls was cooled uniformly, while the lower half of the sidewalls was heated. The Ra, based on the width of the chamber as the characteristic length, was 1.77
108 . They showed a good agreement between the experimental data and the 3D numerical results, with respect to the fluid flow and thermal distributions. Li et al. [43] also carried out experiments and 3D numerical simulations of flow and heat transfer in a cylindrical hydrothermal reactor, again cooled on the upper-half and heated on the lower-half, separated by a baffle. Due to the fluid exchange at the baffle opening, an unsteady jet-like flow pattern in each of the chambers of the reactor was observed. They also concluded that the time-averaged flow profiles and temperature were axially symmetric. Erlekampf et al. [44], conducted numerical simulations of an ammonothermal autoclave. The effect(s) of various baffle shapes on the flow pattern, fluctuations of the temperature as well as flow velocities were studied. They came to the conclusion that 3D numerical simulations were necessary to accurately model the heat and mass transport in the ammonothermal autoclaves, rather than a 2D model. The current authors previously performed a 2D axisymmetric numerical study in a vertical cylinder with laterally-heated walls, using ANSYS FLUENT. A Ra criteria based on flow transition from fully laminar flow to transitional flow was introduced [45]. They considered the cylindrical volume of the reactor to its lateral area as the characteristic length in calculating the Ra. There have been several experimental and numerical investigations with regard to ammonothermal-based GaN crystal growth [46e61]. D'Evelyn et al. [62] used high-pressure

ammonothermal crystal growth to produce high-quality GaN single crystals. The pressures and temperatures were varied between 5 and 20 kbar and 600e1000 C, respectively. They managed to reduce the level of impurity in the GaN crystals and successfully demonstrated homoepitaxial laser diodes. They concluded that this technique could produce high-quality bulk GaN substrates. Bao et al. [63] studied different types of acidic mineralizers such as NH4 Cl, NH4 Br, NH4 I and NH4 F, for the ammonothermal GaN crystal growth method. It was concluded that NH4 F was a promising mineralizer, since it could afford a negative temperature gradient in supercritical ammonia at a temperature range of 550e650 C. Furthermore, it could increase the growth rate and quality of GaN crystals. Mirzaee et al. [64] presented an in-depth 2D axisymmetric numerical study for the ammonothermal GaN crystal growth technique in a Rene-41 autoclave, including fluid flow, heat transfer, and dissolution and crystallization rates. The operating temperature and pressure were 500 C and 1.5 bar, respectively. The top and bottom parts of the external sidewalls had a constant temperature of 525 C and 475 C, respectively. They proposed to create a gap between the nutrient basket and the sidewall in order to have a stronger flow in the nutrient basket. Although, the presence of the gap resulted in a better mixing and transportation of the nutrient particles, it adversely affected the crystal growth rate near the seeds. The proposed model (with the gap) showed a steady growth rate of 92mm per day. There is still a need for detailed experimental observations and numerical analysis of natural convection in a cylindrical enclosure with laterally-heated boundary conditions. The current study presents 3D LES calculations of turbulent natural convection in a laterally-heated cylindrical enclosure at a Ra number of 8:8
106 , Pr ¼ 5:85, and an aspect ratio (H=D) of 5.17, using the commercial CFD software, ANSYS FLUENT. LES uses filtering operations in order to solve the filtered governing equations such that, the turbulent scales are separated. As a result, the large scales are resolved, while the sub-filter scales are modeled using subgrid-scale (SGS) models [65]. Such a methodology is superior to a Reynolds-averaged Navier-Stokes (RANS) approach, which models all scales. Also, LES is computationally cheaper than DNS, which resolves all the length scales. Over the last two decades, LES has demonstrated excellent comparisons with experiments for nonreacting flows. See the review paper by Pope [66] for comprehensive list of LES investigations. Flow patterns and thermal maps are analyzed through contours of instantaneous velocity and temperature on selected planes. An experiment, with a similar setup as the numerical simulations, is also carried out, and flow patterns and velocities are compared to the computational results. In order to assess the effectiveness of the heating in terms of the consistency of the wall temperatures, which has direct implications on the correlation between experiments and computations, several thermocouples are placed on the inner wall of the reactor. For the numerical results validation purposes, the thermal boundary conditions (temperatures of the walls) are recorded and then directly extended to the numerical computations. For that purpose, Type-K surface thermocouples are installed on the inner wall of the enclosure. It was found during the experiments that these thermocouples recorded constant temperatures along the walls (variation of less than 0:3+ C) both in the azimuthal and longitudinal directions, and in the cold and hot sections. Therefore, the simulations are conducted using constant wall temperatures germane to the experiment. This work attempts to provide more information about the flow patterns and the thermal distribution in a cylindrical crystal growth reactor. The unique aspect of the current work is the usage of 3D LES in conjunction with experiments for investigating relatively large laterally-heated cylindrical reactors. It should be noted that the size

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

3

of such enclosures plays a very important role in the overall flow patterns. For instance for a large reactor, as will be demonstrated in the results below, the core flow is entrained by shear rather than buoyancy. 2. Formulation The governing equations employed in this study are the filtered incompressible filtered Navier-Stokes equations subjected to the Boussinesq approximation. Specifically, Equations (1), (2) and (7) refer to the continuity, momentum and energy equations, respectively.

vðui Þ ¼0 vxi

(1)


 d   vðui Þ v ui uj 1 vsij 1 vp 1 vtij þ ¼   þ g b T  Tref vt r vxj r vxj r vxj vxj

(2)

In Equation (2), sij is the stress tensor due to molecular viscosity and tdij is the deviatoric SGS stress tensor, defined in Equations (3) and (4):

"

vui vuj sij ¼ m þ vxj vxi

!# (3)

1 3

tdij ¼ tij  tkk dij ¼ 2mt Sij

(4)

The SGS turbulent viscosity (mt ) in Eq. (4) is defined according to the classical Smagorinsky model as [67]:

 

mt ¼ rðCs DÞ2 S

(5)

In Equation (5), D is the filter width, Cs is the Smagorinsky   qffiffiffiffiffiffiffiffiffiffiffiffiffi   constant and S ¼ 2Sij Sij . The strain-rate tensor (Sij ) for the resolved scale is defined by,

1 vui vuj Sij ¼ þ 2 vxj vxi

! (6)

The dynamic Smagorinsky model is employed in the study here and according to this model, Cs is calculated based on Germono's identity. More details about this model can be found in Ref. [68]. The energy equation is defined by,


 vT v ui T v þ ¼ vt vxj vxj



nt vT

!

Prt vxj

(7)

where Prt is the subgrid Prandtl number. The non-dimensional parameters namely the Rayleigh, Grashof and Prandtl numbers are given by

Ra ¼

g bDTL3

na

; Gr ¼

Ra n ; Pr ¼ a Pr

(8)

3. Computational details 3.1. Geometry and corresponding numerical domain The 3D geometrical configuration is shown in Fig. 1 Fig. 2 shows a 2D cross-section along the axis of the 3D cylindrical reactor on the

Fig. 1. 3D geometry.

X  Y plane. The internal part (core) of the cylinder and the radius of the caps are 53.975 mm and 25.4 mm, respectively. The heights of the cold section, insulator and hot section are 296.32 mm, 12.7 mm and 462.28 mm, respectively. These values used in the numerical simulation, represent an in-house experimental apparatus designed on the basis of an actual industrial crystal growth autoclave. The working fluid (water) is considered to be Newtonian, incompressible and with constant transport properties except for density in the body force term (Eq. (2)), by virtue of the Boussinesq approximation in the Navier-Stokes equation. Fig. 3 shows two axial planes in the domain chosen to analyze the flow patterns and the thermal distribution (in Section 5). Fig. 4 shows the corresponding top view of these planes. Two additional planes are shown in this figure mainly for the purpose of grid-dependence analysis presented in Section 3.4. In this figure, the angle between the planes 2 and 3 is equal to 60 . Planes 1 and 7 are perpendicular to each other as well. Fig. 5(a) and (b) show the structured and unstructured mesh in the X  Y plane and in the B  B plane, respectively. The number of elements in this mesh is around 3.4 million. It should be noted that the grid is refined closer to the walls, in order to accurately resolve the boundary layers. The mesh density close to the wall is such that there are 8 elements within 6 mm (See Fig. 6(a)). Fig. 6(b) shows the boundary layer flow near the wall. Note that the velocity vectors in this figure are sized and colored by magnitude. The velocity increases across and away from the wall up to the free stream flow indicated by the black line in Fig. 6(b), and then starts to decrease. Obviously, the free stream flow here will be non-uniform due to the cylindrical geometry and turbulent flow characteristics. Also, note that the decrease in

4

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

Fig. 4. Four different planes from the top view.

a laterally-heated cylindrical enclosure are presented here. The initial conditions are such that the numerical domain presented in Fig. 1 was initialized with a uniform temperature of 300 K while the fluid was considered at rest. The two sections of the reactor, the lower hot and upper cold sections, have constant temperature of 303.15 K and 296.15 K, respectively. The upper and lower caps, and the insulator are considered to be adiabatic (see Fig. 2). 3.3. Numerical set up

Fig. 2. 2D cross section of geometry in X  Y plane with boundary conditions; Also shown are positions of two points in mm for mesh convergence study and line and plane locations for flow and thermal analysis.

As illustrated in Fig. 5(a) and (b) and described in Section 3.1, a non-uniform structured and non-structured mesh was employed. Equations (1)-(7) along with the boundary and initial conditions, were numerically solved by utilizing ANSYS FLUENT. This commercial software employs a finite volume method to discretize the governing equations [69]. The pressure-velocity coupling was applied by using the PISO algorithm (pressure implicit with splitting of operator). A second order implicit formulation was utilized for the unsteady flow computation. The bounded central differencing scheme was employed to discretize the momentum and energy equations. 3.4. Mesh convergence study

Fig. 3. Positions of different planes in the reactor.

velocity in the core is due to the fact the finite amount of momentum contributed by the boundary layer is imparted to the entire central core of the fluid in the reactor. 3.2. Boundary and initial conditions Transient calculations of natural convection at Ra ¼ 8:8
106 in

A mesh convergence study was performed for three different grids (547000, 1139000 and 3434000 elements). Two separate points (A and B in Fig. 2) were chosen in the domain and the fast Fourier transform (FFT) of the velocity component in the Y direction was calculated, and the results are shown in Fig. 7(a) and (b) for the two finer grid sizes. Both points A and B are studied in the X  Y plane (see Fig. 4). The difference between the amplitudes of the two grid sizes for a random frequency value, as seen in Fig. 7(a) and (b), is small. The mean velocity values in the Y direction and for the points A and B within the time frame of 200 s were also compared with the corresponding experimental values. The difference between the numerical and experimental values for the mesh of 1:14
106 was 26% and 22% for points A and B, respectively. The corresponding value for the mesh size of 3:43
106 was 14.7% and 6.1% for the points A and B, respectively. Based on these results and the FFT plots, the mesh of 3:43
106 was selected for further study. It can also be seen from Fig. 7(a) and (b) that in the high frequency zone a 5/3 slope is observed in the log-log scale as shown in Fig. 7(a) and (b), which implies isotropic turbulence at the small scales [65]. The turbulent kinetic energy magnitude theoretically is a function of velocity components. In the specific physical problem studied here, the values of X and Z velocity components were (not shown here) much lower compared to the Y velocity components. From Fig. 7(a) and (b), it appears that the experimental data is closer to the numerical results obtained with the coarser mesh rather than the finer one. This behavior can be explained as follows: The experimental data indicates that there is a buildup of energy at the high frequencies, due to insufficient resolution of the sampling.

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

5

Fig. 5. Mesh on two different planes (a) Plane 1 (X  Y) (b) Plane B  B in Fig. 2.

Fig. 6. (a) Exaggerated view of the mesh from Fig. 5(b), near the wall, (b) boundary layer flow near the wall from Fig. 5(a).

This will also be observed and explained in figures later on in the paper. Based on this aspect of the experimental data, it is expected that the coarser mesh numerical simulations would match the data better than the finer mesh results.

the Intel Xeon processor, depending on the grid density of the calculations.

3.5. Time convergence study

4.1. Experimental overview

A time step of 0.005 s was used for the transient computations presented here. The same time step was used for the gridindependent studies involving the different grid densities as well. It should be noted that a time convergence study was not carried out for 3D simulations due to the computational effort involved with conducting multiple 3D calculations. However, a comprehensive time convergence analysis was indeed performed for 2D simulations using the same geometry and boundary conditions in a previously published paper [70]. It was found in the previous study that, for the same Ra and Pr studied here, (Ra ¼ 8:8
106 ; Pr ¼ 5:85), a time step of 0.01 s was appropriate for successful convergence. By inference and extrapolation, a more conservative time step of 0.005 s was chosen for the present study. A convergence criterion of 106 was employed for the solutions of the continuity, momentum and the energy equations. Typical simulation times ranged from 800 to 1100 CPU hours on 28 cores of

The experimental reactor geometry is the same as the one used in the numerical experiments. It should be noted that the purpose of the experimental facility is to provide direct comparison of flow and temperature results with the numerical model, thereby contributing to the validation of the latter. The temperature boundary conditions implemented during the experiment were also used for the activation of the numerical model for the natural convection environment. In the experiment, certainly the environment is fully 3D. However, procedure used for the non-intrusive flow visualization of the full flow field, isolated 2D longitudinal sections in the central symmetry region of the reactor. These are sections of the reactor that are used to establish the twodimensionality of the flow, and provide qualitative and quantitative validation of the numerical predictions. The test apparatus (Fig. 8) consists of an acrylic cylindrical tube representing the reactor vessel, which in turn is encased in two

4. Experiments

6

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

Fig. 7. Fast Fourier transform of velocity component in Y direction for two grids of 1.14 mil and 3.43 mil and in X  Y plane (a) Point A (b) Point B.

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

7

Fig. 8. Details of the test section with the cylindrical reactor and surrounding constant temperature bath jackets (a) Picture of the reactor tube inside the heating jackets (b) Exploded solid drawing view showing the relationship between the jackets, reactor tube and the jackets separator.

parallelepipedic enclosures separated by a Teflon insulator. These external enclosures are designed to function as constant temperature baths, the lower one containing the hot temperature fluid source, while the upper containing the cooler fluid, thereby engendering a natural circulation environment in the reactor. These two separated baths are meant to provide isothermal boundary conditions for the cylindrical test reactor along the portions of its

Fig. 10. The flow in each subsection (marked A through E) was measured in a separate instance in time. The flow outside of these sections was concealed from the camera, and therefore could not be measured.

inner wall exposed to the respectively, hot and cold surfaces. The numerical simulations were conducted to reproduce this particular

Fig. 9. Test loop showing the reactor, constant temperature hot and cold bath jackets and the sink reservoirs connected to each one of the jackets. Sink reservoirs are connected through constant volume pumps to their corresponding hot and cold jackets.

8

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

experimental environment. The fluid (water) in each of the constant temperature baths is circulated from two larger thermal sink reservoirs (hot and cold) by means of two separate circulating pumps (Fig. 9). These volumetric flow pumps are meant to provide high flow rates to their respective parallelipipedic enclosures, such that nearly constant respectively, hot and cold temperatures are maintained along the reactor wall surfaces. The temperatures in the two sink reservoirs are maintained constant by means of programmable logic controllers (PLC) controls. A series of Type-K thin foil, fast response thermocouples were used to measure the wall surface temperatures of the reactor cylinder during each run. These temperature measurements consisted of a total of 60 points on the outside surface of the tube (15 axial positions and 4 circumferential positions) and 5 points on the inside surface. In addition Type-K, 0.25 mm bid thermocouples were used to measure the inlet and exit temperatures from the constant temperature enclosures encasing the test reactor. The designated measurement error of the individual thermocouples factory was ± 1 K. However, this error was limited by calibrating all thermocouples at the same time against a consistent standard. Bias error was eliminated through the adjustment of the calibration curves in the data acquisition system and the subsequent overall

error in the measurement of the temperature difference between the hot and cold inner surfaces was less than 0.55 K. The flow visualization method requires seeding of the fluid with light reflecting neutrally buoyant tracer particles that were in average 10 mm in size. A Lexel laser in combination with an optical train of spherical and cylindrical lenses is used to transform the cylindrical beam (2 mm in diameter and wavelength of 514 nm) of the laser into a thin sheet of light. When passing through the test section, this light sheet isolates a longitudinal cross section that is approximately 0.5 mm thick. Such a thin section has the ability of displaying particles that execute their motion predominantly in a 2D plane, thus giving a 2D rendition of the flow. A tomographic procedure, not followed in this work, can identify longitudinal slices throughout the entire test section thus creating a succession of 2D cuts in a 3D flow. The laser used in these experiments is a continuous wave laser yielding 5 W. The high speed camera used to acquire images was a Photron APX RS powered by Photron FASTCAM Viewer software. The exposure time for an individual 1024
1024 pixel image was 0.008 s (125 Hz) and 50 sequential frames (0.4 s) were used and added digitally to create one image that reconstructed the trajectories of the particle paths, which rendered a qualitative image of

Fig. 11. Visualization of flow structures in subsections A and B of the reactor from Fig. 10; (a),(c) PIV images of the flow at two different time instances b) Flow trajectories for a period of 0.4 s; Temperature difference between hot and cold sections is DT ¼ 7  C and Ra ¼ 8:8
106 .

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

the flow. Since the laser yielded a continuous wave, the camera frequency of acquiring images was used as the de facto strobe of the flow motion. Using the camera action, a series of sequential images of the change in the particles ’positions with time was captured at a rate of 125 fps. These images provide a picture of the flow, thus able to visualize movements of both the bulk flow and that of the boundary layers. When processed using particle image velocimetry (PIV), these sequences of qualitative images yield a quantitative 2D vector field that can resolve the individual particle velocities in an x  y system of coordinates.

9

For the present set of experiments, the high aspect ratio of the reactor combined with the limited length of the longitudinal sheet of light, makes it necessary to capture the images in a succession for five subsections, in order to retain a high image resolution of the flow velocity vectors in the queried sections. Furthermore, the current experimental setup, due to the construction of the reactor and heating jackets, prevents the uninterrupted visual interrogation of the entire cross section. The small sections at the top, bottom and center of the reactor, as well as the edge of the velocity field along the surface of the reactor were not accessible to the camera due to curvature. Additionally, some of the thermocouples used to

Fig. 12. Visualization of flow structures in subsections C, D and E of the reactor from Fig. 10; Temperature difference between hot and cold sections is DT ¼ 7  C and the Ra ¼ 8:8
106 ; a) PIV processing of the flow b) flow trajectories for a period of 0.4 s.

10

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

gather temperature information for the boundary conditions obstructed the flow imaging process. Details of the subsections and concealed areas are displayed in Fig. 10. The consequence of this limitation, is that the presented velocity field is discontinuous in some regions. However, the reader will be able to obtain a very good overall understanding of both qualitative and quantitative features in this overall reactor. The qualitative information was extracted using a PIV procedure embedded in the Insight 4G software. A deformation grid was chosen with 16
16 pixel interrogation regions. This grid engine algorithm allows the greatest accuracy by allowing the image to deform to accommodate the velocity gradients within the windows. A standard FFT algorithm was used to correlate the paired

images fit with a Gaussian curve. Here, one pixel represented 149 mm. A time step of 56 ms was chosen between the paired images. This time step allowed for the faster moving particles to utilize the full scale resolution without too many pairs being lost due to movement out of the window of consideration. The velocity fields calculated using the PIV procedure were validated using the corresponding pathlines that were created by summing sequences of raw images in a certain time span. The length of these particle trails represented a velocity of the flow in that region. When compared at several points within the field, the deviation between the PIV calculated- and path length measuredvelocities were within 5 percent. The process was then repeated with a smaller time step to verify that error due to out of plane

Fig. 13. Qualitative flow in subsection C at four different relative, but successive instances in time (t ¼ 0 is an arbitrary starting point): a) t ¼ 0 b) t ¼ 4 c) t ¼ 8 and d) t ¼ 12. Temperature difference between the hot and cold sections is DT ¼ 7  C and Ra ¼ 8:8
106 .

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

Fig. 14. Fast Fourier transform of the velocity component in Y direction for the mesh size of 3,430,000 on the six planes for (a) Point A, and (b) Point B.

11

12

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

particle movement was controlled. 4.2. Experimental results The experimental flow patterns presented in Figs. 11e13 are segregated according to the A, B, C, D, and E frames defined in Fig. 10. The results presented in Figs. 11 and 12 were collected at 1200 s from the start of the experiment. At this time, it was established that thermal equilibrium was achieved throughout the enclosure, based on the steadiness of the temperatures at the wall, and the observed flow patterns, which presented a quasi-steady structure. Note that even though the results of the experiments and simulations are not (and cannot) be compared at the same time instance, the comparison is meaningful due to the fact that the flow reached a quasi-steady structure which keeps repeating itself,

presenting the same flow scales and formations, even though the latter are in a state of flux. Also, further note that this time discrepancy situation is due to the fact that the initialization of the heating process in the experiments take a longer time to die out and reach constant wall temperature, while in the numerical simulations we are starting directly with constant wall temperatures. Fig. 11(a) presents the quantitative PIV reconstruction of the flow and in the A and B sections of the reactor, while Fig. 11(b) presents the digital summation of approximately 50 images acquired at intervals of 0.008 s. The pathlines presented in Fig. 11(b) help reconstruct the flow and provide locations of vortices as well as streams of fluid circulating between the vortices. The interesting locations, 1, 2, 3 and 4 of such flow formations are marked with red circles. The snaking flow straddling both A and B subsections is also marked. The corresponding images reconstructed using PIV are

Fig. 15. (a) Velocity contours and (b) Velocity vectors on the X  Y plane (Z ¼ 0) at time, t ¼ 849 s; Vectors are sized uniformly and colored by velocity magnitude.

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

indicated in Fig. 11(a). The color of the vectors (and not their length) is indicative of their magnitude. The low velocity vectors are consistent with the indicated darker regions, Fig. 11(b), while the longer and more luminous streaks are indicative of larger flow velocities. The punctuated luminous dots mark regions of highly visible low velocity particles. Fig. 12 presents the flow reconstruction in sections C, D and E of the hot section, using the same procedure described above for sections A and B. Vortex regions of high and low velocity are identified along with streams of flow that circulate between the vortices. Like in Fig. 11, a quantitative PIV vector is obtained using Insight 4G. The presentation in parallel of the qualitative flow reconstruction (Fig. 12(b)) and the PIV vector map (Fig. 12(a)) is meant for offering concomitantly an image of the raw flow and its corresponding PIV digital embodiment. In subsequent sections, it will be the PIV vector map that will be compared

13

to the numerically simulated flow. In both Figs. 11 and 12, the flow structures presented in the separate sub-sections fit together relatively well, considering that each has been obtained at slightly different time instances (but within the same experimental session). It is important to note that while the locations of various flow structures fluctuate, the overall structure remains relatively constant. That is due to the fact that the flow is relatively slow and more importantly quasi-steady. Some of the structures evolve and even repeat at a large time scale, even though they may change slightly in size and location. Such an example is exemplified in Fig. 13 which presents subsection C in isolation, during a sequence of four time steps at a 4 s interval, spanning a total of 12 s. Several eddies appear with seemingly arbitrary movement across time steps. However, the general pattern of the flow remains

Fig. 16. (a) Velocity contours and (b) Velocity vectors on the X  Y (Z ¼ 0) plane at time, t ¼ 917 s; Vectors are sized uniformly and colored by velocity magnitude.

14

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

consistent. These results point to the conclusion that although much of the flow is seemingly random in nature, the general spatial characteristics of the flow stay the same, and can be described. 5. Numerical simulations This section includes numerical analysis and investigations of flow patterns and temperature distribution in a cylindrical domain as seen in Fig. 1 subjected to the prescribed boundary conditions described in Section 3.2. Results are also compared with the experimental data. Six azimuthal planes 60 apart (see planes 1, 2, 3, 4, 5 and 6 in Fig. 4), are chosen for the purpose of comparing predictions with data. A FFT of the Y- velocity values from the experiments and the

numerical simulations was calculated for two points, A and B (see Fig. 2), on each of the planes. The FFT outcomes are shown in Fig. 14. Fig. 14(a) displays the FFT of the Y-velocity values for point A in all the planes (planes 1, 2, 3, 4, 5 and 6 in Fig. 4) for a time window of 200 s in a log-log scale for both experiments and computations. This time window was chosen such that it was large enough to capture the flow features taking place across multiple flow cycles both in the numerical and experimental results. Note that as it was observed both in experiments and corresponding computations that, as the flow reached a quasi-steady formation for the respective Ra, various vortices or recirculation zones with their respective length scales disappeared and reappeared periodically. Due to the repetitive nature of the flow, qualitative and meaningful comparisons can be made between the results of the experiments and

Fig. 17. (a) Velocity contours and (b) Velocity vectors on the X  Y (Z ¼ 0) plane at time, t ¼ 980 s; Vectors are sized uniformly and colored by velocity magnitude.

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

numerical calculations. For the FFT analysis conducted here, the time windows, and not the start times are meaningful, with the former chosen to be the same for the experiments and simulations. In the numerical simulations, the time window of 780 se980 s was selected to study the FFT results. In Fig. 14(a) for instance, P2A represents point A in Fig. 2, on the plane 2 as seen in Fig. 4. Similarly, P5A represents point A in Fig. 2, on the plane 5 as seen in Fig. 4, and so on. Firstly, as it can be seen in Fig. 14, the predictions from the corresponding planes compare well with each other indicating that overall, the flow fluctuations are similar on all the planes. Comparison with experiments shows that there is a buildup of energy at the high frequencies for the experiments, indicating an insufficient resolution of the sampling of experimental data. Similar information is provided in Fig. 14(b) for point B. Here though, there is a better comparison between experiments and simulations. On quantitatively comparing velocity values between one of the observed planes in the experimental setup and 6 different planes of the numerical calculations, the maximum difference between the mean Y velocity values in a time window of 200 s, was 14.7% and 6.1% for points A and B, respectively. These results in addition to the FFTs presented earlier imply that the velocity fluctuations could be higher for point A compared to point B. This is the reason for the insufficient resolution observed in the experimental data for the high frequencies in the point A results. For further numerical analysis presented in the next section, two planes, the X  Y and Y  Z planes (see Fig. 4), were chosen and results are presented in Section 5.1.

15

5.1. Flow distribution Figs. 15e17 show the velocity magnitude contours and velocity vectors on the X  Y plane at three different times corresponding to t ¼ 849, 917, and 980 s, respectively. The velocity magnitudes in Figs. 15e17 are normalized by the value of 1 m/s. It is clear that at all times, near the walls, the flow moves downward in the upper section and upward in the lower section of the reactor. In the vicinity of the insulator in the domain, eddies form and the velocity values are higher in comparison to the lower or higher parts of the reactor far from the walls. Here, the velocity magnitudes are higher, because this is the region where the hot and cold fluids meet each other and mixing occurs, resulting in more recirculation and randomness in the flow. Figs. 15e17 clearly exhibit the wall boundary layer regions (for example see encircled regions in Fig. 15). These boundary layers with the help of additional shear are instrumental in moving the fluid core of the cold and hot sections. These results are qualitatively similar to the experimental results. As shown in Figs. 15(b), 17(b) and 11(c), the velocity vectors are mainly from bottom to top in the cold section of the reactor. However, as the flow has randomness, the directions of the velocity vectors can change in time as shown in Figs. 16(b) and 11(a). The velocity vectors on the upper right corner are from top to bottom, but in the lower center, they are from bottom to top. The maximum velocities in both the experiments and numerical simulations are similar to each other (approximately 0.018 m/s). Earlier, it was stated that the velocity magnitude away from the mixing zone (Note: mixing zone is the region close to the insulator, where the mixing between hot and cold fluids occurs) and towards

Fig. 18. Instantaneous velocity contours at time ¼ 980 s on two different X  Y planes corresponding to locations (a) Z ¼ 0:03 m and (b) Z ¼ 0:03 m, and (c) five different X  Z planes.

16

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

the two caps of the reactor were smaller than those in the mixing zone itself. This can be also be seen in the instantaneous velocity contours shown on two planes in Fig. 18 for a time of 980 s in the simulation. Fig. 18 (a) and (b) display the velocity contours on the X  Y planes at Z ¼ 0:03 m and Z ¼ 0:03 m, respectively. Contours of instantaneous velocity are shown in Fig. 18 (c) for X  Z planes corresponding to Y ¼ 0:15 m, Y ¼ 0:3 m, Y ¼ 0:45 m, Y ¼ 0:6 m and Y ¼ 0:75 m from the origin, and again they show higher velocity regions being present closer to the mixing zone. When one compares the flow in the X  Y symmetry plane (Fig. 18) to the flow in the X  Y planes at Z ¼ 0:03 m (Fig. 18), it is clear that the flow in

the symmetry plane is more organized and regular, in comparison to the Z ¼ 0:03 m, X  Y planes, where many additional small scales seem to be present. This also implies that the turbulent flow presented here is not axisymmetric, and a 2D solution of the problems might not be appropriate in order to describe the physical flow accurately. The area averaged mean Y-velocity and X-velocity profiles across five different planes are shown in Fig. 19. The five planes are corresponding to Y ¼ 0:4 m, Y ¼ 0:45 m, Y ¼ 0:5 m, Y ¼ 0:55 m and Y ¼ 0:6 m. Firstly, the Y-velocities are higher in magnitude than the X-velocities. This is expected due to the prescribed wall

Fig. 19. Mean radial (a) Y-velocity and (b) X-velocity profiles across the planes corresponding to Y ¼ 0:4 m, Y ¼ 0:45 m, Y ¼ 0:5 m, Y ¼ 0:55 m and Y ¼ 0:6 m.

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

temperature gradients in the Y-direction. Secondly, at all locations, the velocities are higher, close to the boundary layer than in the center of the domain, especially in the case of the Y-velocities, further reaffirming the observations from the velocity contours

17

presented earlier. In the case of the X-velocities, this is true only for the Y ¼ 0:5 m plane, whereas for the rest of the planes (Y ¼ 0:45 m, Y ¼ 0:55 m and Y ¼ 0:6 m), the flow stays quite uniform. This nonuniformity on Y ¼ 0:5 m can be expected, since this location is near

Fig. 20. Isosurfaces of vorticity magnitude for values corresponding to (top to bottom) 1.25 (1/s), 1.5 (1/s), 1.75 (1/s) and 2 (1/s), at time t ¼ 980 s.

18

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

the mixing zone where the vortices form and move toward the upper and lower sections of the reactor, resulting in more chaotic movements and thereby larger X-velocity fluctuations. Fig. 20 presents the different iso-surfaces of the vorticity magnitude in the domain at the time instance of 980 s, corresponding to values of 1.25 s1, 1.5 s1, 1.75 s1 and 2 s1. The isosurfaces illustrate the low-turbulent nature of the flow, and interestingly this is restricted to the mixing region near the interface between the hot and cold zones. This again is further evidence that high velocity and recirculation zones are mostly present only in a very small region near the insulator of the reactor.

5.2. Thermal distribution Fig. 21 shows the instantaneous contours of the normalized temperature in the domain on two different planes, the X  Y and Y  Z planes. The temperature was normalized as follows:

T ¼

T  TL TH  TL

(9)

where TL is the lowest temperature or the temperature applied on the upper cold wall. As it is clear on both the planes, the temperature is almost uniform in the entire hot section in the bottom, whereas there are some variations in the upper cold section. In addition, there is an inversion of the temperature field in the upper and lower parts of the enclosure. In other words, overall, the upper section is warmer than the lower section. This happens, because

the hot fluid on the boundary of the lower half has enough momentum to penetrate the upper cold section and heat it, and the reverse occurs with the cold fluid from the top. The implications of an almost uniform temperature in the entire domain are that the buoyancy forces are considerably weakened, potentially leading to slower flow in the interior regions far away the walls in the domain. So it is evident that flow occurring in the reactor is primarily driven by the temperature gradients on the boundary and not in the interior of the reactor. This further indicates that conduction plays a negligible role in heat transfer within the domain. Fig. 22 shows the normalized mean temperature profiles across four radial lines (see lines A-D in Fig. 2) on the two planes X  Y and Y  Z. The lines correspond to Y ¼ 0:2 m, Y ¼ 0:4 m, Y ¼ 0:6 m and Y ¼ 0:8 m. As it can be seen from the figures, there is a sharp temperature gradient near the walls and is positive and negative in the upper and lower sections, respectively. However, in the interior regions of the reactor the temperature is quite uniform as was observed in the contours. 5.3. Comparison of the results in LES and RANS models In a previous study by the current authors [70], a RANS-based k-

u SST model was employed to perform 2D axisymmetric calculations in order to understand the transition from laminar to turbulent flow, specifically identifying the transition Ra for a fixed thermal and geometrical configuration. The velocity and temperature contours for Ra ¼ 8:8
106, using a 2D axisymmetric model, are presented in Fig. 23(a) and (c), respectively. Also shown here are

Fig. 21. Instantaneous normalized temperature contours on the (a)X  Y ðZ ¼ 0Þ and (b)Y  Z ðX ¼ 0Þ planes at time, t ¼ 980 s.

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

19

Fig. 22. Normalized mean temperature contours on the (a) X  Y and (b)Y  Z planes at time, t ¼ 980 s.

velocity and temperature contours from 3D RANS calculations of the same geometry, except for the curvatures on the corners, in Fig. 23(b) and (d). It should be noted that the 2D configurations of Fig. 23(a) and (c) represent one-half of the reactor, based on the axis of symmetry, while for the 3D configurations of Fig. 23(b) and (d), one uses the X  Y plane of symmetry applied to the full reactor. A comparison between the maximum turbulence intensity magnitudes in 2D and 3D RANS simulations, showed relatively similar values (0.42% and 0.43%). A comparison between Figs. 17(a), 23(a) and 23(b) shows that velocity distribution in the 2D simulation is different, although the maximum velocity magnitudes are similar. The same kind of behavior is observed by comparing Figs. 21(a), 23(c) and 23(d). These comparisons lead the authors to conclude that a 3D simulation is essential, in order to accurately study the heat and mass transport in this type of reactor. Similar conclusions were also noted in the study by Erlekampf et al. [44] as well. Furthermore, a comparison between Figs. 21(a) and 23(d), shows that temperature distribution by the RANS model is not as accurate, since it illustrates more uniform temperatures away from the walls in the domain. However, in an inarguably more accurate calculation like LES, the lower section is now cooler than the upper section of the reactor. Li et al. [37] reported a similar thermal distribution in their experimental results. All of this is further justification for using LES to accurately predict the thermal distributions in the domain as well. In a separate study by the current authors [71], it was shown that by introducing a baffle with a specific opening area, between the hot and cold zones, such a thermal distribution could be flipped, i.e. the lower hot section stays

warmer than the upper cold section. These two different thermal distributions have potential implications for the design of a crystal growth reactor, for instance, retrograde vs. forward solubility in the ammonothermal crystal growth. 6. Conclusions 3D numerical simulations of natural convection in a laterallyheated vertical cylindrical reactor, with lower hot and upper cold zones, were conducted using LES (for the first time, based on the authors' knowledge) in the commercial CFD software, ANSYS FLUENT. The problem investigated here consists of a temperature gradient driven buoyant flow, characterized by Pr ¼ 5:85 and Ra ¼ 8:8
106 , based on the characteristic length defined as the ratio of the volume of the reactor to its lateral area. The reactor had an aspect ratio of 5.17. Velocity vectors, vorticity iso-surfaces and contours of mean temperature and velocity on several axial and transverse planes were presented in this paper. An in-house cylindrical reactor was used in an experimental effort, in order to validate the numerical results with the identical geometry and boundary conditions. Flow patterns in the cylindrical reactor were captured using a PIV method, and the velocity magnitude was calculated from the experimental data. Both the experiments and numerical simulations showed that mixing of the fluids primarily happened in the region close to the interface between the hot and cold zones (mixing zone). In this region, velocity magnitudes were higher, as evident from results on the orthogonal planes, in comparison to the upper and lower

20

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21

Fig. 23. (a) and (c) 2D RANS velocity and temperature contours, respectively; (b) and (d) 3D RANS velocity and temperature contours, respectively; in the central symmetry region, all at t ¼ 980 s.

sections of the reactor away from the walls and the interface. Overall, the results of the numerical simulations were qualitatively similar to those of the experiments in terms of flow patterns and maximum velocity magnitudes. In addition, two physical points were selected in the cylindrical enclosure, and the mean values of Y-velocity were calculated for both the numerical simulation and the experiment. For the two points, the differences between the numerical and the corresponding experimental values were 6% and 14%. FFT of the Y-velocity component for the two specific points were also calculated both from the experiment and the numerical calculation, and comparisons showed reasonable qualitative agreement between the two. A uniform temperature distribution was observed away the walls inside the reactor. This happens as the thin wall boundary layers exist for the current reactor size and Ra, near the walls. In fact, the thermal resistance between the bulk flow and the boundary layer region is weak. Although the buoyancy force is weak due to the absence of large temperature gradients in the core regions of the reactor, the shear force is primarily responsible for entraining the flow in these areas. Note that the buoyant flow is dominant in the boundary layer regions and it is the main driver of the flow inside the reactor. An other interesting thermal phenomena observed here which was the temperature inversion inside the

reactor. The presence of the warmer flow in the upper cold section of the reactor and cooler fluid in the lower hot section of the reactor can affect the choice of the mineralizer. Finally, 2D and 3D RANS calculations were also presented as part of an effort to demonstrate the necessity of employing accurate 3D modeling of turbulence natural convection in laterally-heated cylinders. References  ski R, Wysmołek A, Baranowski J, Kamin  ska M, Doradzin  ski R, [1] Dwilin  ski J, et al. GaN synthesis by ammonothermal method. Acta Phys Pol A Garczyn 1995;88(5):833e6.  ski R, Baranowski JM, Kamin  ska M, Doradzin  ski R, Garczyn  ski J, [2] Dwilin Sierzputowski L. On GaN crystallization by ammonothermal method. Acta Phys Pol A 1996;90(4):763e6. [3] Purdy Andrew P. Ammonothermal synthesis of cubic gallium nitride. Chem Mater 1999;11(7):1648e51. [4] Ehrentraut Dirk, Meissner Elke, Bockowski Michal. Technology of gallium nitride crystal growthvol. 133. Springer Science & Business Media; 2010. [5] Moldovan Stefan Ilie. Numerical simulation and experimental validation of fluid flow and mass transfer in an ammonothermal crystal growth reactor. PhD thesis. University of Akron; 2013. [6] Bodoia JR, Osterle JF. The development of free convection between heated vertical plates. J Heat Transf 1962;84(1):40e3. [7] Flack RD. The experimental measurement of natural convection heat transfer in triangular enclosures heated or cooled from below. J Heat Transf 1980;102(4):770e2.

H. Enayati et al. / International Journal of Thermal Sciences 116 (2017) 1e21 [8] Poulikakos D, Bejan A. Natural convection experiments in a triangular enclosure. J heat Transf 1983;105(3):652e5. [9] Poulikakos Dimos, Bejan Adrian. The fluid dynamics of an attic space. J Fluid Mech 1983;131:251e69. [10] Ostrach S. Natural convection in enclosures. J Heat Transf 1988;110(4b): 1175e90. [11] Fusegi Toru, Min Hyun Jae, Kuwahara Kunio. A numerical study of 3D natural convection in a cube: effects of the horizontal thermal boundary conditions. Fluid Dyn Res 1991;8(5e6):221. [12] Andrew E, McLeod, Bishop Eugene H. Turbulent natural convection of gases in horizontal cylindrical annuli at cryogenic temperatures. Int J Heat Mass Transf 1989;32(10):1967e78. [13] Barakos G, Mitsoulis E, Assimacopoulos Do. Natural convection flow in a square cavity revisited: laminar and turbulent models with wall functions. Int J Numer Methods Fluids 1994;18(7):695e719. [14] Lee TL, Lin TF. Three-dimensional natural convection of air in an inclined cubic cavity. Numer Heat Transf Part A Appl 1995;27(6):681e703. [15] Salmun Haydee. Convection patterns in a triangular domain. Int J Heat Mass Transf 1995;38(2):351e62. [16] Salmun Haydee. The stability of a single-cell steady-state solution in a triangular enclosure. Int J heat mass Transf 1995;38(2):363e9. [17] Asan H, Namli L. Numerical simulation of buoyant flow in a roof of triangular cross-section under winter day boundary conditions. Energy Build 2001;33(7):753e7. [18] Guo Yanhu, Bathe Klaus-Jürgen. A numerical study of a natural convection flow in a cavity. Int J Numer methods fluids 2002;40(8):1045e57. [19] Muhammad AR, Sharif, Mohammad Taquiur Rahman. Natural convection in cavities with constant flux heating at the bottom wall and isothermal cooling from the sidewalls. Int J Therm Sci 2005;44(9):865e78.  [20] Xam an J, Alvarez G, Lira L, Estrada C. Numerical study of heat transfer by laminar and turbulent natural convection in tall cavities of facade elements. Energy Build 2005;37(7):787e94. [21] Badr HM, Habib MA, Anwar S, Ben-Mansour R, Said SAM. Turbulent natural convection in vertical parallel-plate channels. Heat mass Transf 2006;43(1): 73e84. [22] Li Hongmin, Xing Changhu, Braun Minel J. Natural convection in a bottomheated top-cooled cubic cavity with a baffle at the median height: experiment and model validation. Heat mass Transf 2007;43(9):895e905. [23] Lo DC, Young DL, Tsai CC. High resolution of 2D natural convection in a cavity by the DQ method. J Comput Appl Math 2007;203(1):219e36. [24] Hussain Salam Hadi, Hussein Ahmed Kadhim. Natural convection heat transfer in a differentially heated square enclosure with a heat generatingconducting circular cylinder at different diagonal locations. In: 6th international advanced technologies symposium; 2011. p. 13e8. , Menezo Christophe, Bouia Hassan, [25] Daverat Christophe, Pabiou Herve Xin Shihe. Experimental investigation of turbulent natural convection in a vertical water channel with symmetric heating: flow and heat transfer. Exp Therm Fluid Sci 2013;44:182e93. [26] Iyi Draco, Hasan Reaz. Natural convection flow and heat transfer in an enclosure containing staggered arrangement of blockages. Procedia Eng 2015;105:176e83. [27] Hsieh Shou-Shing. Thermal correlation of natural convection in bottomcooled cylindrical enclosures. J Thermophys Heat Transf 1990;4(1):123e6. [28] Edwards DK, Catton Ivan. Prediction of heat transfer by natural convection in closed cylinders heated from below. Int J Heat Mass Transf 1969;12(1):23e30. [29] Huang Durn-Yuan, Hsieh Shou-Shing. Analysis of natural convection in a cylindrical enclosure. Numer Heat Transf Part A Appl 1987;12(1):121e35. [30] Lin YS, Akins RG. Pseudo-steady-state natural convection heat transfer inside a vertical cylinder. J heat Transf 1986;108(2):310e6. [31] Papanicolaou E, Belessiotis V. Transient natural convection in a cylindrical enclosure at high Rayleigh numbers. Int J Heat Mass Transf 2002;45(7): 1425e44. [32] Chakir A, Souli Mhamed, Aquelet Nicolas. Study of a turbulent natural convection in cylindrical annuli of gas-insulated transmission lines 400 kV. Appl Therm Eng 2003;23(10):1197e208. [33] Barhaghi Darioush G. A study of turbulent natural convection boundary layers using large-eddy simulation. Chalmers University of Technology, Department of Applied Mechanics; 2007. [34] Malik Asif Hussain, Shah Ajmal, Khushnood Shahab. CFD analysis of heat transfer within a bottom heated vertical concentric cylindrical enclosure. In: Journal of physics: conference seriesvol. 439. IOP Publishing; 2013. p. 012004. [35] Kang Gyeong-Uk, Chung Bum-Jin, Kim Hyoung-Jin. Natural convection heat transfer on a vertical cylinder submerged in fluids having high prandtl number. Int J Heat Mass Transf 2014;79:4e11. [36] Jiang Yan-Ni, Chen Qi-Sheng, Prasad V. Numerical simulation of ammonothermal growth processes of GaN crystals. J Cryst Growth 2011;318(1):411e4. [37] Li Hongmin, Braun Minel J, Evans Edward A, Wang G-X, Paudel Govind, Miller Jason. Natural convection flow structures and heat transfer in a model hydrothermal growth reactor. Int J heat fluid flow 2005;26(1):45e55. [38] Wang Buguo, Callahan MJ, Rakes KD, Bouthillette LO, Wang S-Q, Bliss DF, et al. Ammonothermal growth of GaN crystals in alkaline solutions. J Cryst Growth 2006;287(2):376e80. [39] Masuda Yoshio, Suzuki Akira, Ishiguro Tohru, Yokoyama Chiaki. Numerical simulation of heat and fluid flow in ammonothermal GaN bulk crystal growth process. Jpn J Appl Phys 2013;52(8S):08JA05.

21

[40] Masuda Y, Suzuki A, Mikawa Y, Kagamitani Y, Ishiguro T, Yokoyama C, et al. Numerical simulation of GaN single-crystal growth process in ammonothermal autoclaveeeffects of baffle shape. Int J Heat Mass Transf 2010;53(5): 940e3. [41] Li Hongmin, Xing Changhu, Braun Minel J. Flows in a lower half heated upper half cooled cylindrical model reactor loaded with porous media. Heat mass Transf 2007;43(11):1201e11. [42] Masuda Yoshio, Sato Osamu, Tomida Daisuke, Yokoyama Chiaki. Convection patterns and temperature fields of ammonothermal GaN bulk crystal growth process. Jpn J Appl Phys 2016;55(5S):05FC03. [43] Li Hongmin, Braun Minel J, Xing Changhu. Fluid flow and heat transfer in a cylindrical model hydrothermal reactor. J Cryst growth 2006;289(1):207e16. [44] Erlekampf J, Seebeck J, Savva P, Meissner E, Friedrich J, Alt NSA, et al. Numerical time-dependent 3D simulation of flow pattern and heat distribution in an ammonothermal system with various baffle shapes. J Cryst Growth 2014;403:96e104. [45] Enayati Hooman, Chandy Abhilash J, Braun Minel J. Numerical simulations of natural convection in a laterally-heated cylindrical reactor. In: American society of thermal and fluids engineers. NYC; 2015. [46] Nernst Walther. Theoretische chemie vom standpunkte der Avogadro’schen regel und der thermodynamik. Stuttgart: Verlag Von Ferdinand Enke; 1907. [47] Douglas R, Ketchum, Kolis Joseph W. Crystal growth of gallium nitride in supercritical ammonia. J Cryst Growth 2001;222(3):431e4. [48] Yoshikawa A, Ohshima E, Fukuda T, Tsuji H, Oshima K. Crystal growth of GaN by ammonothermal method. J Cryst growth 2004;260(1):67e72. [49] Popov VN, Tsivinskaya Yu S, Bekker TB, Kokh KA, Kokh AE. Numerical investigation of heat-mass transfer processes in hydrothermal growth system. J Cryst growth 2006;289(2):652e8. [50] Pendurti S, Chen Q-S, Prasad V. Modeling ammonothermal growth of GaN single crystals: the role of transport. J Cryst growth 2006;296(2):150e8. [51] Ehrentraut Dirk, Kagamitani Yuji, Fukuda Tsuguo, Orito Fumio, Kawabata Shinichiro, Katano Kizuku, et al. Reviewing recent developments in the acid ammonothermal crystal growth of gallium nitride. J Cryst Growth 2008;310(17):3902e6. [52] Hashimoto Tadao, Wu Feng, James S, Speck, Nakamura Shuji. Growth of bulk GaN crystals by the basic ammonothermal method. Jpn J Appl Phys 2007;46(10L):L889. [53] Fukuda Tsuguo, Ehrentraut Dirk. Prospects for the ammonothermal growth of large GaN crystal. J Cryst Growth 2007;305(2):304e10. [54] Ehrentraut Dirk, Kagamitani Yuji, Yokoyama Chiaki, Fukuda Tsuguo. Physicochemical features of the acid ammonothermal growth of GaN. J Cryst Growth 2008;310(5):891e5. [55] Ehrentraut Dirk, Fukuda Tsuguo. The ammonothermal crystal growth of gallium nitride, a technique on the up rise. Proc IEEE 2010;98(7):1316e23. [56] Chen Qi-Sheng, Jiang Yan-Ni, Yan Jun-Yi, Li Wei, Prasad V. Modeling of ammonothermal growth processes of GaN crystal in large-size pressure systems. Res Chem Intermed 2011;37(2e5):467e77. [57] Chen Qi-Sheng, Yan Jun-Yi, Jiang Yan-Ni, Li Wei. Modeling on ammonothermal growth of GaN semiconductor crystals. Prog Cryst Growth Charact Mater 2012;58(2):61e73. [58] Theresia MM, Richter, Niewa Rainer. Chemistry of ammonothermal synthesis. Inorganics 2014;2(1):29e78. [59] Bao Quanxi, Saito Makoto, Hazu Kouji, Kagamitani Yuji, Kurimoto Kouhei, Tomida Daisuke, et al. Ammonothermal growth of GaN on a self-nucleated GaN seed crystal. J Cryst Growth 2014;404:168e71. [60] Kucharski R, Zajac M, Puchalski A, Sochacki T, Bockowski M, Weyher JL, et al. Ammonothermal growth of GaN crystals on HVPE-GaN seeds prepared with the use of ammonothermal substrates. J Cryst Growth 2015;427:1e6. [61] Bockowski M, Iwinska M, Amilusik M, Fijalkowski M, Lucznik B, Sochacki T. Challenges and future perspectives in HVPE-GaN growth on ammonothermal GaN seeds. Semicond Sci Technol 2016;31(9)(093002). [62] Evelyn MPD, Hong HC, D-S Park, Lu Hai, Kaminsky E, Melkote RR, et al. Bulk GaN crystal growth by the high-pressure ammonothermal method. J Cryst growth 2007;300(1):11e6. [63] Bao Quanxi, Saito Makoto, Hazu Kouji, Furusawa Kentaro, Kagamitani Yuji, Kayano Rinzo, et al. Ammonothermal crystal growth of GaN using an nh4f mineralizer. Cryst Growth & Des 2013;13(10):4158e61. [64] Mirzaee Iman, Charmchi Majid, Sun Hongwei. Heat, mass, and crystal growth of GaN in the ammonothermal process: a numerical study. Numer Heat Transf Part A Appl 2016;70(5):460e91. [65] Pope Stephen B. Turbulent flows. Cambridge university press; 2000. [66] Pope Stephen B. Ten questions concerning the large-eddy simulation of turbulent flows. New J Phys 2004;6(1):35. [67] Smagorinsky J. General circulation experiments with the primitive equations. I. The basic experiment. Mon Weather Rev 1963;91(3):99e164. [68] Sagaut P. Large eddy simulations of incompressible flows. Berlin: SpringerVerlag; 2005. [69] Co Ansys. Ansys fluent users guide. Version 16.2. Ansys. Ansys. 2015. [70] Enayati Hooman, Chandy Abhilash J, Braun Minel J. Numerical simulations of transitional and turbulent natural convection in laterally heated cylindrical enclosures for crystal growth. Numer Heat Transf Part A Appl 2016;70(11): 1195e212. [71] Enayati Hooman, Chandy Abhilash J, Braun Minel J. Numerical simulations of turbulent natural convection in laterally-heated cylindrical enclosures with baffles for crystal growth. In: IMECE 2016, phoenix; Nov-2016.