- Email: [email protected]
Large eddy simulation (LES) calculations of natural convection in cylindrical enclosures with rack and seeds for crystal growth applications
International Journal of Thermal Sciences 123 (2018) 42e57
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
Large eddy simulation (LES) calculations of natural convection in cylindrical enclosures with rack and seeds for crystal growth applications Hooman Enayati a, Abhilash J. Chandy b, *, Minel J. Braun a a b
Department of Mechanical Engineering, University of Akron, Akron, OH, 44325, USA Department of Mechanical Engineering, Indian Institute of Technology Bombay, Mumbai, 400076, India
a r t i c l e i n f o
a b s t r a c t
Article history: Received 31 January 2017 Received in revised form 27 July 2017 Accepted 24 August 2017
This paper presents three-dimensional (3D) large eddy simulation (LES) calculations of natural convection in laterally-heated cylindrical enclosures using ANSYS FLUENT. The problem considered here is at a Rayleigh number, ðRaÞ of 8:8 106 , based on a characteristic length scale defined by the ratio of the reactor's volume to the lateral area. The numerical simulations are carried out for three different geometrical configurations: an empty reactor, a reactor with a rack and 32 seeds, and thirdly, a reactor with a rack and 80 seeds mounted inside. The objective of the current study is to understand the effects of the rack and seeds on the flow patterns and thermal distribution inside the reactor. In order to investigate these effects, a comprehensive analysis of velocity and temperature on the selected planes and points inside the domain were performed. Results showed that with an increase in seeds within the reactor, the flow slowed down more, the turbulence intensity reduced, and the extent of temperature inversion increased. These have important implications for the choice of ammonothermal method in terms of acidic or basic, and also for species transport and reactions, which eventually affect the deposition rates of GaN crystals in autoclaves. © 2017 Elsevier Masson SAS. All rights reserved.
1. Introduction Natural or free convection, one of the modes of heat transfer, occurs due to density differences, resulting from concentration or temperature differentials. This phenomena can be observed both in the engineering and natural applications. For instance as an industrial application, natural convection heat transfer plays a critical role in a crystal growth reactor. High-quality Gallium Nitride (GaN) substrates grown in such reactors can be used in light-emitting diodes (LEDs), laser diodes (LDs), UV detectors and highfrequency/high-power electronic devices and automotive applications [1]. There are several methods including hydride vapor phase epitaxy (HVPE), solution growth, sodium flux and ammonothermal techniques, that are used to produce GaN crystals. These techniques have their own advantages and disadvantages [2,3]. An ammonothermal crystal growth technique is an effective
* Corresponding author. E-mail addresses: [email protected] (H. (A.J. Chandy), [email protected] (M.J. Braun).
Enayati),
http://dx.doi.org/10.1016/j.ijthermalsci.2017.08.025 1290-0729/© 2017 Elsevier Masson SAS. All rights reserved.
[email protected]
way to grow GaN bulk single crystals with a high quality, but slow growth rate (100 mm/day). This technique uses ammonia as the solvent, and is similar to the hydrothermal technique, which is widely utilized to grow quartz crystals using water as the solvent. There have been several studies in the past utilizing both these techniques [4e7]. In the ammonothermal growth technique, the reactor has two separated thermal zones. This thermal configuration is essential in order to create a buoyant flow inside the autoclave. The preloaded GaN crystals, located in the nutrient basket, are dissolved in the dissolution zone and transported by virtue of natural convection to the crystallization zone, where GaN crystals are crystallized on seeds. As GaN crystals are almost insoluble in pure ammonia, mineralizers are added in this process. Two types of mineralizers can be used in an ammonthermal technique inside the reactor; acid and basic ([8]). In the case of a basic mineralizer such as KNH2 , the solubility of GaN crystals decreases at elevated temperatures, also called retrograde solubility. On the other hand, in the case of an acidic mineralizer such as NH4 CL, there is a direct relation between the solubility and the temperature. There is a recognized relationship between the type of mineralizer, the position of the growth crystals in the reactor, and the distribution and
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
position of the hot and cold zones. Generally, the choice of mineralizer is dependent upon the design regarding the type of crystal to be grown and the thermal distribution inside the reactor, respectively. In the present paper, the position of the rack and seed crystals replicates the architecture of a prototype of a GaN reactor. Furthermore, a baffle can be used in order to control the flow direction, and possibly create a stronger temperature gradient inside the reactor. Experimental flow visualization and measurements of temperature inside a crystal growth reactor are near impossible due to the opaque walls and harsh environments, resulting from high pressures and temperatures. Also, optimization of the growth process by performing several new sets of experiments can be expensive and time-consuming. Through computational fluid dynamics (CFD), it is possible to understand the ammonothermal crystal growth process better, at a relatively cheaper cost and with a reasonable accuracy. Finally, it is also important to realize that a proper understanding of the flow and temperature fields is necessary in order to gain knowledge about the mass transfer mechanisms, which in turn is critical in optimizing the crystal growth process. During the past few decades, several studies have been conducted to numerically investigate natural convection in an enclosure. They included two-dimensional (2D) and three-dimensional (3D) studies in various kids of geometries including rectangular enclosures [9e14]. The thermal boundary conditions in those studies were mainly side-wall or top-bottom-wall heating configurations. In addition to constant heat flux or temperature applied on the boundaries, heat generation sources inside the domain have also been employed. The current study benefits from a temperature gradient between the top and bottom halves of the lateral walls of the reactor. Two different but constant temperatures are applied laterally on the walls onto two zones, separated by an adiabatic section. Such a thermal configuration is relatively new as it is different from the more-widely used side-wall heating and topbottom-wall heating configurations. Convective heat transfer in a cylindrical enclosure with different thermal boundary conditions have been studied numerically and experimentally over the past several years as well [15e21]. Yuan et al. [22] investigated natural convective flow in a horizontal concentric annuli with different inner shapes including cylindrical, square, triangular and elliptical cross sections. The outer and inner walls were considered to be at constant temperature. They studied the effects of different inner shapes and surface radiation on the flow field and heat transfer, and validated their results with available literature. They concluded that at high temperature levels, radiation had a very important role in the overall heat transfer behavior. They also established correlations for the average Nusselt number. Malik et al. [23] performed numerical and experimental studies of axial temperature gradient and the heat transfer within an enclosure. The experimental set up included different materials and diameters for the inner and outer cylinders, respectively. The thermal boundary condition included different bottom disc temperatures ranging from 353 to 433 K. They validated their numerical results with experimental data. While several investigations have been carried out focusing on the role of the natural convection heat transfer in a crystal growth process [24e30], more work is still needed. Li et al. [31] presented a 3D numerical simulation of a cylindrical autoclave with an aspect ratio (height/diameter) of 10, subjected to circumferentially nonuniform constant temperature on the upper and lower walls of the autoclave, with and without a baffle in the middle of the enclosure. They concluded that for Ra ¼ 3:84 1011, calculated
43
based on the autoclave diameter, this nonuniform heating condition had significant effects on the fluid flow and temperature distribution in the bulk flow. For instance, temperature deviation on the upper chamber wall, resulted in a distortion of temperature isotherms in the upper growing chamber. Li et al. [32] studied a rectangular enclosure both experimentally and numerically, the latter both in 2D and 3D. The sidewalls on the lower half of the enclosure were heated, while those on the upper half were cooled uniformly. They found a good agreement between the 3D numerical results and the experimental data with respect to the thermal and fluid flow distributions. They concluded that the 2D model predictions compared qualitatively well with experiments. However, there was a 11% difference of total heat transfer rate between the 2D and 3D results. Li et al. [33], numerically studied 3D simulations of flow and heat transfer in a cylindrical hydrothermal reactor and presented them along with an experimental set up. The thermal configuration included a cooling and heating procedure on the upper and lower halves, respectively. An unsteady jet-like flow pattern in each chamber of the reactor was observed due to the fluid exchange at the baffle opening. They concluded that the growth environment was mainly affected by the total heating rate. Erlekampf et al. [34], performed a numerical analysis of an ammonothermal process for the bulk growth of nitride crystals. They explored the impacts of various baffle shapes on the flow velocity and temperature fluctuations in the domain. The current authors [35,36], used a 2D axisymmetric RANS model to perform numerical simulations of a vertical cylinder with laterally-heated walls. The transition regime from laminar flow to turbulent flow was defined based on the Ra criteria. The approach followed in the current study in terms of the definition of the Ra is similar to [36], where the characteristics length is half of the reactor's radius or R=2. Furthermore, effects of a baffle inside an ammonothermal crystal growth reactor have also been investigated using 3D large eddy simulation (LES) [37]. Results demonstrated the presence of temperature inversion for specific baffle openings in the reactor. Several numerical and experimental investigations specifically with regard to GaN crystal growth using the ammonothermal technique have also been studied [38e44]. Yoshida et al. [45] developed autoclaves in order to grow high quality GaN crystals under a certain temperature ( < 800 C) and pressure ( < 150 MPa) conditions. A retrograde solubility of GaN crystals using an acidic mineralizer ðNH4 CLÞ, was achieved for an operating pressure of 110 MPa and temperature values of 650 C and higher. They concluded that a high-temperature acidic ammonothermal technique could yield low-etch pit density (around 104 cm2) and improve the quality of the GaN crystals. Ehrentraut et al. [46] investigated the solubility of GaN in an ammonothermal method containing an acidic mineralizer ðNH4 CLÞ. It was found that NH4 CL was a good choice among the acidic mineralizers for GaN crystal growth in terms of lattice distortion and its concentration played a key role in controlling the solubility of GaN. 2. Scope of work In spite of the several studies in the past on GaN crystal growth investigations, there has not yet been a comprehensive and fundamental analysis of the flow and heat transfer patterns focusing on the effects of different geometric configurations related to the number of seeds in an ammonothermal crystal growth reactor. ANSYS FLUENT numerical simulations of turbulent natural convection flow in an empty reactor, and a reactor with 32 and 80 seeds, are presented, in order to understand the influence of
44
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
internal furniture on convective flow, thermal maps, and potentially mass transfer mechanisms, all of which affect the crystal growth rates. Temperature distributions and flow patterns on several points and planes are analyzed and compared among the cases, all for a laterally heated cylindrical enclosure at Ra ¼ 8:8 106 . The characteristic length scale used here for the Ra definition is the ratio of the volume of the enclosure to its lateral area, resulting in half of the radius or R=2, where R is the radius of the enclosure. For the first time, LES, along with the dynamic Smagorinsky model for subgrid-scale (SGS) closure is employed for simulations of an ammonothermal crystal growth reactor. In all cases, the aspect ratio (total length/diameter) is 5.17 and Pr ¼ 5:85.
In the above equations, ð…Þ represents the LES filtering operation, ui is the velocity component, p is the pressure, r is the density, g is the acceleration due to gravity, b is the coefficient of thermal expansion, T is the temperature, Tref is the reference temperature, sij is the stress tensor due to molecular viscosity and tij is the subgrid-scale (SGS) stress tensor, defined in Eqs. (3) and (4):
3. Governing equations
The SGS turbulent viscosity, nt , in Eq. (4) is defined according to the classical Smagorinsky model as [47]:
The governing equations used in this study are the incompressible filtered Navier-Stokes equations subjected to the Boussinesq approximation. The continuity, momentum and energy equations are presented in Equations (1), (2) and (7), respectively, as follows:
vðui Þ ¼0 vxi vsij 1 vp vtij vðui Þ v ui uj þ ¼ þ g b T Tref : vt vxj r vxj vxj vxj
"
vui vuj sij ¼ m þ vxj vxi
!#
1 3
tij tkk dij ¼ 2nt Sij
(3)
(4)
nt ¼ ðCs DÞ2 S
(5)
In Eq. (5), D is the filter width, Cs is the Smagorinsky constant qffiffiffiffiffiffiffiffiffiffiffiffiffi and S ¼ 2Sij Sij . The strain-rate tensor ðSij Þ for the resolved scale
(1)
is defined by,
(2)
1 vui vuj Sij ¼ þ 2 vxj vxi
! (6)
The dynamic Smagorinsky model is employed in the study here
Fig. 1. 3D geometry for the (a) empty reactor (b) reactor with a rack and 32 seeds (c) reactor with a rack and 80 seeds (d) Y Z cross-section of the reactor with rack and 80 seed (Fig. 1c).
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
45
Fig. 3. Two-dimensional cross section of Fig. 1a on X Y plane and thermal boundary conditions.
where Prt is the subgrid Prandtl number. The non-dimensional parameters namely the Rayleigh ðRaÞ, Grashof ðGrÞ and Prandtl ðPrÞ numbers are given by:
Ra ¼
gbDTL3
na
; Gr ¼
Ra n ; Pr ¼ Pr a
(8)
Here L is the characteristic length scale, DT is the temperature difference, n is the kinematic viscosity, and a is the thermal diffusivity. 4. Numerical details 4.1. Numerical domain
Fig. 2. Dimensions of the rack and seeds in Fig. 1c (in mm) (a) Front view (b) Top view.
and according to this model, Cs is calculated based on Germono's identity [48]. The energy equation is given by:
vT v ui T v þ ¼ vt vxj vxj
nt vT Prt vxj
! (7)
The three-dimensional geometry of the empty reactor, and the reactor with a rack and 32 and 80 seeds are shown in Fig. 1aed. Fig. 2a and b presents the front and top views of Fig. 1c. A schematic of the cylindrical reactor along with the thermal boundary conditions on the different sections are presented in Fig. 3. Please note that all three cases in the current study have the same thermal boundary conditions. The dimensions of the different parts of the reactor (hot wall, cold wall, insulator, the radius of the reactor, rack and seeds), used in the numerical simulations and presented in Figs. 1e3, are based on an in-house experimental apparatus. All the seeds have a square cross section with a length of 50.84 mm. By comparing the top view of Fig. 1b and c, it is clear that their geometrical dimensions are similar, but the number of the seeds and their positions are different (e.g. 32 and 80 seeds). Also, note that the seeds and the rack (Fig. 1b and c) are located below the coordinate system indicated in Fig. 1a. This can also be understood
46
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
Fig. 4. Grids on two different planes for the empty reactor with 6:4 106 elements (a) X Y plane cut in Fig. 1(a), (b) Y Z plane cut in Fig. 1(a), (c) Exaggerated view near the wall in Fig. 4b to see the elements near the wall.
Table 1 Physical properties of the working fluid. Property
Value 3
Density (kg/m ) Heat capacity (J/kg-K) Thermal conductivity (W/m-K) Viscosity (kg/m-s) Thermal expansion coefficient (1/K)
997.1 4180.9 0.61028 0.00085384 0.000257
by comparing the dimensions shown in Figs. 2a and 3. The mesh distribution for the empty reactor on two different planes is shown in Fig. 4aec. The domain includes mainly hexadehdral elements with an average skewness and orthogonal quality of 6.06 102 and 0.985, respectively. The grid density near the wall in the empty
Table 2 Positions of the radial planes/points ðX ¼ Z ¼ 0Þ from the reference coordinate system in the Y direction in m. Point
1
2
3
4
5
Y position Point Y position Point Y position
0.467 6 0.172 11 þ0.1
0.408 7 0.113 12 þ0.2
0.349 8 0.054 13 þ0.3
0.290 9 þ0.025
0.2317 10 þ0.05
reactor case is such that there are 8 elements within 4 mm (See Fig. 4c), in order to be able to capture the velocity gradient in the boundary layers near the wall.
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
47
Fig. 5. Velocity vectors colored by the instantaneous velocity magnitude in the vicinity of walls; (a) X Y plane cut (b) Exaggerated view near the mixing zone.
4.2. Boundary and initial conditions 3D transient LES calculations of natural convection in an ammonothermal crystal growth reactor at Ra ¼ 8:8 106 are presented in the current study. For all the simulations, the entire computational domain (shown in Fig. 1) was initialized by considering the working fluid, i.e. water at rest along with a uniform temperature of 300 K. Constant temperatures 303.15 K and 296.15 K are imposed on the lower and upper sections of the reactor, respectively. The lower and upper caps of the reactor, and the midregion between the hot and cold sections are considered to be adiabatic (See Fig. 3). 4.3. Numerical set up The filtered governing equations (Eqs. (1)e(7)) along with the initial and boundary conditions, provided in Section 4.2, are numerically solved by utilizing the finite volume approach in the commercial software, ANSYS FLUENT [49]. The pressure implicit with splitting of operator (PISO) algorithm is considered for the pressure-velocity coupling and the bounded central differencing scheme is employed to spatially discretize the governing equations. Temporal discretization of the governing equation was accomplished using a second order implicit formulation. The properties of the working fluid, by considering the Boussinesq approximation, are provided in Table 1, and they were used in all three simulations presented in this study. 4.4. Mesh convergence study The current authors showed that for the numerical case
involving an empty reactor (with the same boundary conditions of the current study), total mesh elements of 3:4 106 showed a converged solution [50]. This conclusion was further verified by using experimental data. Although grid points of 3:4 106 was sufficient for the current calculations, the number of elements in this study for the empty reactor was increased such that all the three cases including those with the rack and seeds have a comparable type of mesh and resolution. In this study, the total number of elements for all three cases including empty reactor, reactor with rack and 32 seeds and the reactor with a rack and 80 seeds are 6,400,000, 9,640,000 and 8,564,000, respectively. 4.5. Time convergence study A time step of 0.005 s was used in all simulations here. A time convergence study was not performed for the 3D simulations, due to the computational effort involved with conducting multiple time-steps for different geometries. For instance, the simulation for the empty reactor with 6:4 106 elements took approximately 1300 CPU hours on 28 cores of the Intel Xenon processor, with a time step of 0.005 s. However, a time convergence study for a 2D simulation of the empty reactor with a same geometry and boundary conditions was indeed performed as part of a previous study by the current authors [36]. It was concluded that for the same Rayleigh ðRa ¼ 8:8 106 Þ and Prandtl ðPr ¼ 5:85Þ numbers, a time step of 0.01 s was appropriate for a converged solution in 2D. By extrapolation, half of this time step is considered in all simulations in the current study, thereby leading to a time step of Dt ¼ 0:005 s (See Ref. [50] as well). A convergence criteria of 5 105 for the solution of the continuity equation and 1 105 for the momentum and energy equations are considered in the
48
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
Fig. 6. Normalized mean velocity contours on the selected planes for the empty reactor over a time window of 139.5 s (a) Z ¼ 28.575 mm (b) Z ¼ - 28.575 mm (c) Z ¼ 9.525 mm (d) Z ¼ - 9.525 mm; Legend is valid for all the four plots.
numerical calculations here. 5. Results and discussion Numerical analysis of flow and heat transfer for each geometrical configuration, illustrated in Fig. 1, is presented in this section. Three geometrical configurations, including an empty reactor, a reactor with a rack and 32 seeds and a reactor with a rack and 80 seeds, are computationally modeled. Normalized velocity and temperatures are discussed in detail in Sections 5.1e5.2. For the natural convection problem in this paper, there isn't a pre-specified velocity scale, and hence, velocity is normalized using unity or 1 m/ s for convenience, whereas temperature is normalized such that L T ¼ TTT , where TH and TL are the high and low temperatures, H TL
respectively. The four axial planes of 1, 2, 3 and 4 shown in Fig. 2b, are used to study the flow patterns and temperature distributions for the analysis of each geometrical case. Planes 1, 2, 3 and 4 are located at þ28.575 mm, þ9.525 mm, 28.575 mm and 9.525 mm, respectively, all away the origin in the Z-direction and on X Y planes. Each plane is located in the middle of two consecutive seeds. Thirteen planes perpendicular to the Y-direction are also considered, in order to study the time-averaged velocity magnitude
for a time window of 139.5 s. The positions of these planes are presented in Table 2. Furthermore, thirteen points on the centerlines of these planes are also chosen to study the time-averaged Yvelocity for a similar time window. Table 2 also includes the position of all these points. It is expected that the presence of obstacles like rack and seeds affects the flow patterns and the temperature gradients inside a crystal growth autoclave. Hence, it is important to analyze the flow and thermal distributions inside the three geometrical configurations of the current study.
5.1. Flow patterns Figs. 6e8 present the contours of normalized mean velocity on the selected planes, shown in Fig. 2b, for the empty reactor, reactor with a rack and 32 seeds and a rack and 80 seeds over a time window of 139.5 s. Although the normalized mean velocities are in the same range of values on those planes for all the three cases (0e0.015), the velocity contours do vary qualitatively case by case. The velocity values in the central regions of the planes are relatively higher than other regions of the reactor (like regions close to the lower and upper caps), except in the vicinity of walls, where the high velocities are observed in the boundary layer. This
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
49
Fig. 7. Normalized mean velocity contours on the selected planes for the reactor with a rack and 32 seeds over a time window of 139.5 s (a) Z ¼ 28.575 mm (b) Z ¼ - 28.575 mm (c) Z ¼ 9.525 mm (d) Z ¼ - 9.525 mm; Legend is valid for all the four plots.
central region is referred to as the mixing zone, since the hot and cold fluids meet each other, and exchange momentum and enthalpy. In fact, the lower hot fluid layers move up and the upper cold fluid layers move down, and interact with each other in this region, near the lateral walls. This also results in random fluid motions at higher velocities and vortical formations at different scales, which in turn aids in the overall mixing and circulation within the reactor. Fig. 5 shows the fluid layers near the wall and the momentum boundary layer formation in the lower hot and upper cold sections of the empty reactor. Note that the velocity magnitude in the boundary layer increases away from the walls to the free stream portion of the flow and then decreases again. As the flow is turbulent in this cylindrical geometry, the thickness of the boundary layer varies as the free stream flow is non-uniform inside the reactor.
5.1.2. Reactor with rack and 32 seeds Fig. 7 shows the normalized mean velocity contours on the same planes, but for the reactor that has 32 seeds mounted on the rack. The presence of the seeds makes the flow with higher velocity values circulate mainly near the mixing zone or in the upper part of the reactor. This is far different from the case of the empty reactor. In fact, the main difference between Figs. 6 and 7 is the flow distribution in the lower section of the reactor. The flow is slower in this region for the reactor with seeds, due to the presence of the obstacle. This occurs as the flow cannot move freely enough in order to exchange momentum within the domain. However, the flow circulation and vortices in the mixing zone can move towards the upper section of the reactor easily, and can maintain higher velocities, since there are no obstacles in that section. In addition, the higher mean velocity values exist in the vicinity of the lateral walls, rather than the interior regions of the reactor in both cases.
5.1.1. Empty reactor Fig. 6aed shows the mean velocity distribution on several planes in the empty reactor. There are similar features, like the relatively higher velocities in the mixing zone and the decrease in values as the flow moves towards the caps on the top and bottom.
5.1.3. Reactor with rack and 80 seeds Similarly, Fig. 8 shows the normalized mean velocity contours on the corresponding planes for the case with 80 seeds mounted on the rack inside the reactor. The flow distribution on each plane is similar to the case with 32 seeds (see Fig. 7). However the
50
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
Fig. 8. Normalized mean velocity contour on the selected planes for the reactor and a rack and 80 seeds over a time window of 139.5 s (a) Z ¼ 28.575 mm (b) Z ¼ - 28.575 mm (c) Z ¼ 9.525 mm (d) Z ¼ - 9.525 mm; Legend is valid for all the four plots.
maximum mean velocity here is lower (0.014 m/s), although only slightly, compared to the 32 seeds case, and overall, the velocity distribution in the lower section of the reactor is different, due to the added seeds in the reactor, which further slows down the fluid in the lower part of the reactor. 5.1.4. Additional discussion More information about the effects of different numbers of seeds on the flow can be obtained by analyzing the area-averaged velocity over a time window of 139.5 s for each case of the current study on the selected planes of Table 2. As shown in Fig. 9, the empty reactor generally has the highest velocity magnitude on the planes considered here, in comparison to the other two cases, while the case with 80 seeds has the lowest. The presence of the seeds slows down the flow in the lower section of the reactor as shown in Fig. 9. However, they do not affect the flow upstairs to such an extent in all three cases, as there are no obstacles and in general, the flow is only affected by the buoyancy force in boundary layers and in the central core region by shear-induced entrainment. Note that for each case, the highest averaged velocity magnitudes are on the planes that are located near the mixing zone, for instance, Y ¼ 0:025 m. As the location of a plane moves towards the top and
bottom of the reactor, the velocity magnitude decreases on the corresponding plane. This is expected due to momentum exchange in the mixing zone, where the hot and cold fluid parcels interact with each other, and the chaotic movements and vortices ensue. As a result, mean velocity values on the planes near the mixing zones are higher. Furthermore, Fig. 10 shows the time-averaged Y-velocity component over a time window of 139.5 s for the thirteen points listed in Table 2, for each geometrical configuration. As shown in the figure, the presence of the seeds (32 or 80) reduces the Y-velocity values in the lower section of the reactor, as indicated by the values on Points # 1e8. This behavior is expected as the presence of an obstacle (seeds in this case) in the domain slows down the motions. The direction of the Y-velocity component in the lower section of the reactor also changes with the presence of the seeds. This is important from the perspective of species transport and reactions involving GaN crystals in the ammonothermal technique. Knowledge of fluid flow directions and magnitudes would help in designing an appropriate geometrical configuration for better deposition rates on the seeds inside the crystal growth autoclave. Fig. 10 also shows that the direction of Y-velocity component in the upper part of the reactor does not change with the geometrical
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
Fig. 9. Time-averaged velocity magnitudes on different planes over a time window of 139.5 s for three geometrical configurations.
Fig. 10. Time-averaged Y-velocities for different points over a time window of 139.5 s for three geometrical configurations.
51
52
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
Fig. 11. Resolved turbulent kinetic energy (TKE) on Z ¼ 28.575 mm for the (a) empty reactor (b) eeactor with the rack and 32 seeds (c) eeactor with the rack and 32 seeds.
configuration, and overall, higher Y-velocity values can be expected in comparison to the lower section of the reactor. Another interesting feature to analyze in such problems is the amount of turbulence within each configuration, since the latter will have an effect on the mixing and circulation and therefore, the etching and deposition rates. Turbulent kinetic energy (TKE) is used for that purpose, and the values of TKE and its distribution in each geometry are studied here. TKE is the mean kinetic energy per unit mass for the associated eddies in a turbulent flow. It is characterized by measured root-mean-square ðRMSÞ of the velocity fluctuations in a computational domain. For the current LES study, the resolved turbulence kinetic energy is quantified as follows:
TKE ¼ 0:5 RMS2x þ RMS2y þ RMS2z
(9)
where RMS is the root-mean-square of each velocity component. RMSy for instance can be defined as
RMSy ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðv 〈v〉Þ2
where v and 〈v〉 represent the instantaneous and mean Y-velocity components, respectively. Fig. 11 presents the TKE on a selected plane of Z ¼ 28:575 m for three different geometries. As is clear from the figure, the presence of the seeds decreases the maximum value of the TKE as a result of the decreased velocities due to the obstacles. The empty reactor has the highest TKE value, and in the regions around the mixing zone, a wider distribution of relatively higher values of TKE can be observed. Fig. 11, also shows that at the
current Ra, the flow is in the low-turbulent regime, since in all the three cases, the TKE is very low in the domain except in the mixing zone. 5.2. Temperature distribution As discussed earlier in Section 1, since GaN is insoluble in pure ammonia, mineralizers are added to strengthen the solubility of GaN crystals in an ammonothermal crystal growth reactor. Two main types of mineralizers, basic and acidic, can be chosen based on the temperature distribution inside the reactor. The thermal distribution inside an ammonothermal crystal growth reactor is a determining factor for the type of the mineralizer used in the nutrient basket that is filled with GaN crystals. For instance, for a retrograde solubility, the nutrient basket should be kept in a cooler zone inside the crystal growth autoclave. In fact, it is also important to know how the temperature distribution would develop inside the reactor as it could affect the etching and deposition rates of GaN crystals in the reactor. 5.2.1. Empty reactor Fig. 12 presents the velocity vectors superimposed on the normalized mean temperature distribution, averaged over 139.5 s, inside the empty reactor. The fluid layers adjacent to the lower lateral hot walls of the reactor are warm. These layers travel upward due to momentum and buoyancy forces, cross the mixing zone and fill the upper core of the reactor with a warm fluid. The cooler fluid layers in the vicinity of upper cold lateral walls travel down and fill up the lower section of the autoclave with the cooler fluid. This
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
53
Fig. 12. Normalized mean temperature contours on the selected planes for the empty reactor over a time window of 139.5 s (a) Z ¼ 28.575 mm (b) Z ¼ - 28.575 mm (c) Z ¼ 9.525 mm (d) Z ¼ - 9.525 mm; Legend is valid for all the four plots.
eventually results in a cooler fluid in the lower section of the reactor and hotter fluid in the upper section. This phenomena is like a temperature inversion, as the warmer flow exists in the upper section of the reactor, while the walls are colder. Similarly, the cooler fluid exists in the lower section of the reactor, where the walls are hotter. Note that in the problem studied here, as the temperature distribution is almost uniform away the walls, the buoyancy forces resulting from the temperature gradients inside the reactor, are weak. However, the boundary layers move up and down along the wall mainly due to the thermal boundary conditions, while the bulk fluid within the core of the reactor move mainly due to entrainment and shear. 5.2.2. Reactor with rack and 32 seeds Fig. 13 shows the velocity vectors superimposed over the normalized mean temperature over the same time window and for the same planes as presented earlier, except this time, for the case with 32 seeds mounted on a rack inside the rector. A discussion of the temperature distribution similar to the empty reactor is valid here too. However, the presence of 32 seeds causes a larger temperature gradient, especially away from the lateral walls, between the lower and upper sections of the reactor. Here, the upper section
of the reactor is slightly warmer in comparison to the empty reactor case. This happens because the rack and seeds in this case obstructs the flow significantly in the lower section of the reactor. So while the warmer fluid is indeed transported up towards the upper section, the colder fluid does not move as freely towards the lower section (compared to the empty reactor) due to the presence of the rack and seeds. This basically results in a lower overall circulation in the reactor, which leads to a less uniform temperature distribution in the core of the reactor, compared to the empty reactor case. 5.2.3. Reactor with rack and 80 seeds Fig. 14 illustrates the velocity vectors superimposed on the normalized mean temperature contours over a time window of 139.5 s, again for the same planes, but for the case with 80 seeds mounted on the rack inside the rector. A similar temperature inversion behavior is observed here too. The current geometrical configuration though provides the warmest upper zone of the reactor among all the three cases. Furthermore, unlike the case with 32 seeds mounted on the rack inside the reactor, the current geometrical configuration of the (80) seeds provides a slightly nonuniform temperature in the lower section of the reactor. Note that the region close to the walls on the lower section of the reactor has
54
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
Fig. 13. Normalized mean temperature contour on the selected planes for the reactor with a rack and 32 seeds over a time window of 139.5 s (a) Z ¼ 28.575 mm (b) Z ¼ 28.575 mm (c) Z ¼ 9.525 mm (d) Z ¼ - 9.525 mm; Legend is valid for all the four plots.
slightly higher temperatures. 5.2.4. Additional discussion As was earlier illustrated through Figs. 7 and 8, the presence of the seeds decreases the velocity magnitude in the lower section of the reactor. The velocity vectors in Figs. 12e14 shows a similar trend as well. The spatial density of these velocity vectors in the domain qualitatively represents the velocity magnitude. The density of the velocity vectors is higher near the walls due to the boundary layer formation for all the three cases. While the density of the velocity vectors in the empty reactor is comparable in the upper and lower sections of the reactor, it is not the case for seeds within the reactor. For instance, here, the velocity vectors are denser in the mixing zone and in the upper section of the reactor, compared to the lower section. This further confirms the existence of a slower flow inside the reactor for the cases with seeds. In addition, it can be seen from Figs. 12e14 that the velocity vectors in the bulk flow of the lower hot section are downward while in the upper cold section of the reactor are upward. The thickness of the thermal boundary layer and normalized mean temperature profiles across six radial lines on X Y plane (See Fig. 3) are shown in Fig. 15. The position of each line is tabulated in Table 3. It is clear from Fig. 15a and b that the temperature distribution is uniform far away the walls on the selected lines for
all the three cases of the study, as was mentioned earlier in the discussion of the temperature contours. The thermal boundary layer formation, where a sharp temperature gradient between the bulk and near wall temperatures exists, can also be noted in these figures. The thickness of the thermal boundary layer increases closer to the mixing zone, while away the mixing zone, the thermal boundary layer is smaller. As discussed earlier for Fig. 5, in the lower section of the reactor, the flow in the vicinity of walls is upwards, while in the upper section of the reactor, the flow is downwards. This results in a boundary layer development from bottom to top, and vise-versa in the lower and upper sections of the reactor, respectively. This would also further confirm that the thermal boundary thickness keeps developing up until the mixing zone. 6. Conclusion A study of natural convection in a laterally-heated cylindrical enclosure was performed by conducting a series of numerical simulations using the commercial CFD software, ANSYS FLUENT. LES calculations of three cases with different geometrical configurations, including an empty reactor, a reactor with a rack and 32 seeds and a reactor with a rack and 80 seeds, were analyzed. Distributions of velocity magnitude, temperature and resolved
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
55
Fig. 14. Normalized mean temperature contours on the selected planes for the reactor with a rack and 80 seeds over a time window of 139.5 s (a) Z ¼ 28.575 mm (b) Z ¼ 28.575 mm (c) Z ¼ 9.525 mm (d) Z ¼ - 9.525 mm; Legend is valid for all the four plots.
Table 3 Positions of the lines (in m) from the reference coordinate system in the Y direction (Z ¼ 0, 0.079375 m < x < 0.079375 m). Point
1
2
3
4
5
6
Y position
0.349
0.231
0.113
þ0.1
þ0.2
þ0.3
turbulent kinetic energy, along with velocity vectors, were presented for each case on selected axial and radial planes located in the domain. All the three cases of study have the same boundary conditions with a Ra ¼ 8:8 106 , based on the characteristic length scale of R2, where R is the radius of the cylinder. In the mixing zone, where the upper cold and lower hot boundary layers meet each other, the fluid exchanges momentum and heat, resulting in mixing and circulation. Vortices of different scales that are formed in and around this region, evolve, disappear and eventually reappear. Furthermore, this mixing zone exhibits the highest velocities in the domain, away from the walls. The presence of the rack and seeds in the lower part of the reactor, decreased the velocity in these regions, due to the added obstructions. The maximum velocity within the domain stayed the same, since it was a function of the velocity in the boundary layer, which in turn solely depends on the thermal boundary conditions that are identical for all the configurations. It was also shown that
the direction and values of the Y-velocity changed dramatically inside the reactor with the added furniture. For instance, the Yvelocity values in the lower section was highest for the empty reactor among all the geometries. It was also noted that the direction of the Y-velocity has potential implications on the deposition rates of GaN crystals inside the reactor, since a favorable flow towards the seeds can bring more GaN crystals from the nutrient basket. This further affirms the fact that a sufficient knowledge of the fluid flow is critical in order to understand species transport and reactions in a crystal growth reactor. In addition, TKE distributions indicated a reduction in turbulence as well with the addition of rack and increase in the number of seeds within the reactor. Due to the thin wall boundary layers that exist for the current problem at a relatively high- Ra, the thermal resistance between the bulk flow and the boundary layer region is weak. This also implies relatively uniform temperature values far away the walls in the reactor. The current thermal distribution inside the domain, results in a weaker flow circulation as the buoyancy force is weak due to the absence of strong temperature gradients. However, a strong buoyancy force does exist near the walls due to the thermal boundary conditions, which is the primary driver of the flow circulation in the reactor. The boundary layers entrain the flow in the core of the reactor via shear, thus creating some movements there as well. In all the three configurations studied here, a temperature
Fig. 15. Normalized mean temperatures on the selected lines for three geometrical configurations over a time window of 139.5 s (a) rhree lines in the lower section of the reactor (b) rhree lines in the upper section of the reactor.
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
inversion was also observed. This phenomenon is evident from the existence of warmer fluid in the upper cold section of the reactor and a colder fluid in the lower hot section of the reactor. This kind of temperature inversion can determine the type of a mineralizer (acidic or basic) in the nutrient basket of a crystal growth reactor. Finally, it was also shown that the thermal boundary layer was the thickest in the mixing zone. Acknowledgments This paper is based upon work supported by the National Science Foundation under Grant No. 1336700. References [1] Kachi Tetsu. Recent progress of gan power devices for automotive applications. Jpn J Appl Phys 2014;53(10):100210. [2] Ehrentraut Dirk, Meissner Elke, Bockowski Michal. Technology of gallium nitride crystal growth, vol. 133. Springer Science & Business Media; 2010. [3] Mamun Molla Md, Paul Sreebash C, Hossain Md Anwar. Natural convection flow from a horizontal circular cylinder with uniform heat flux in presence of heat generation. Appl Math Model 2009;33(7):3226e36. [4] Roux B, Louchart O, Terhmina O. Hydrodynamic aspect of hydrothermal synthesis of quartz bulk flow regimes. Le J de Physique IV 1994;4(C2):C2e3. ski R, Wysmołek A, Baranowski J, Kamin ska M, Doradzin ski R, [5] 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, [6] Dwilin Sierzputowski L. On gan crystallization by ammonothermal method. Acta Phys Pol A 1996;90(4):763e6. [7] Purdy Andrew P. Ammonothermal synthesis of cubic gallium nitride. Chem Mater 1999;11(7):1648e51. [8] Ehrentraut Dirk, Hoshino Naruhiro, Kagamitani Yuji, Yoshikawa Akira, Fukuda Tsuguo, Itoh Hirohisa, et al. Temperature effect of ammonium halogenides as mineralizers on the phase stability of gallium nitride synthesized under acidic ammonothermal conditions. J Mater Chem 2007;17(9):886e93. [9] Chen Q-S, Prasad V, Chatterjee A. Modeling of fluid flow and heat transfer in a hydrothermal crystal growth system: use of fluid-superposed porous layer theory. J heat Transf 1999;121(4):1049e58. [10] 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. [11] 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. [12] 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. [13] Kahwaji Ghalib, Ali Omar M. Numerical investigation of natural convection heat transfer from square cylinder in an enclosed enclosure filled with nanofluids. Int J Recent Adv Mech Eng (IJMECH) 2015;4(4). [14] Iyi Draco, Hasan Reaz. Natural convection flow and heat transfer in an enclosure containing staggered arrangement of blockages. Procedia Eng 2015;105:176e83. [15] Hsieh Shou-Shing. Thermal correlation of natural convection in bottomcooled cylindrical enclosures. J Thermophys Heat Transf 1990;4(1):123e6. [16] Li YH, Lin KC, Lin TF. Computation of unstable liquid metal convection in a vertical closed cylinder heated from the side and cooled from above. Numer Heat Transf Part A Appl 1997;32(3):289e309. [17] Choukairy Kh, Bennacer R, Beji H, Jaballah S, El Ganaoui M. Transient behavior inside a vertical cylindrical enclosure heated from the side walls. Numer Heat Transf Part A Appl 2006;50(8):773e85. [18] Barhaghi Darioush G. A study of turbulent natural convection boundary layers using large-eddy simulation. Chalmers University of Technology, Department of Applied Mechanics; 2007. [19] Malik Asif Hussain, Shah Ajmal, Khushnood Shahab. CFD analysis of heat transfer within a bottom heated vertical concentric cylindrical enclosure. J Phys: Conf Ser 2013;439:012004. [20] 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. [21] Moukalled F, Darwish M, Kasamani J, Hammoud A, Khamis Mansour M. Buoyancy-induced flow and heat transfer in a porous annulus between concentric horizontal circular and square cylinders. Numer Heat Transf Part A Appl 2016;69(9):1029e50. [22] Yuan Xing, Tavakkoli Fatemeh, Vafai Kambiz. Analysis of natural convection in horizontal concentric annuli of varying inner shape. Numer Heat Transf Part A
57
Appl 2015;68(11):1155e74. [23] Asif Hussain Malik, Shahab Khushnood, and Ajmal Shah. Experimental and numerical study of buoyancy driven flow within a bottom heated vertical concentric cylindrical enclosure. 2013. [24] Chen Q-S, Pendurti S, Prasad V. Effects of baffle design on fluid flow and heat transfer in ammonothermal growth of nitrides. J Cryst growth 2004;266(1): 271e7. [25] Li Hongmin, Braun Minel J. Numerical investigation on multi-hole baffle designs for industry hydrothermal synthesis of single crystals. Model Simul Mater Sci Eng 2005;13(8):1249. [26] Chen Q-S, Pendurti S, Prasad V. Modeling of ammonothermal growth of gallium nitride single crystals. J Mater Sci 2006;41(5):1409e14. [27] 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. [28] Jiang Yan-Ni, Chen Qi-Sheng, Prasad V. Numerical simulation of ammonothermal growth processes of gan crystals. J Cryst Growth 2011;318(1):411e4. [29] 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. [30] 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. [31] Li Hongmin, Wang Guo-Xiang, Evans Edward A. Three-dimensional flow of solution in an industry-size hydrothermal autoclave subjected to nonuniform heating effects of a baffle on flow and temperature separation. J Cryst growth 2004;271(1):257e67. [32] 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. [33] 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. [34] 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. [35] Enayati Hooman, Chandy Abhilash J, Braun Minel J. Numerical simulations of natural convection in a laterally-heated cylindrical reactor, vol. 10. NYC, DOI: American Society of Thermal and Fluids Engineers; 2015. [36] 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. [37] H. Enayati, A. J. Chandy, and M. J. Braun. Numerical simulations of turbulent natural convection in laterally-heated cylindrical enclosures with baffles for crystal growth. In IMECE 2016, Phoenix, Nov-2016. [38] Popov VN, S Tsivinskaya Yu, Bekker TB, Kokh KA, Kokh AE. Numerical investigation of heat-mass transfer processes in hydrothermal growth system. J Cryst growth 2006;289(2):652e8. [39] 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. [40] Ehrentraut Dirk, Fukuda Tsuguo. The ammonothermal crystal growth of gallium nitride, a technique on the up rise. Proc IEEE 2010;98(7):1316e23. [41] Ehrentraut Dirk, Pakalapati Rajeev T, Kamber Derrick S, Jiang Wenkan, Pocius Douglas W, Downey Bradley C, et al. High quality, low cost ammonothermal bulk gan substrates. Jpn J Appl Phys 2013;52(8S):08JA01. [42] Richter Theresia MM, Niewa Rainer. Chemistry of ammonothermal synthesis. Inorganics 2014;2(1):29e78. [43] 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. [44] 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. [45] Yoshida Kazuo, Aoki Kensuke, Fukuda Tsuguo. High-temperature acidic ammonothermal method for gan crystal growth. J Cryst Growth 2014;393: 93e7. [46] Ehrentraut Dirk, Kagamitani Yuji, Yokoyama Chiaki, Fukuda Tsuguo. Physicochemical features of the acid ammonothermal growth of gan. J Cryst Growth 2008;310(5):891e5. [47] Smagorinsky J. General circulation experiments with the primitive equations. I. The basic experiment. Mon Weather Rev 1963;91(3):99e164. [48] Sagaut P. Large eddy simulations of incompressible flows. Berlin: SpringerVerlag; 2005. [49] Ansys Co. Ansys fluent users guide. 2014. Version 15.0. Ansys. Ansys. [50] Enayati Hooman, Chandy Abhilash J, Braun Minel J, Horning Nicholas. 3d large eddy simulation (les) calculations and experiments of natural convection in a laterally-heated cylindrical enclosure for crystal growth. Int J Therm Sci 2017;116:1e21.
Contents lists available at ScienceDirect
International Journal of Thermal Sciences journal homepage: www.elsevier.com/locate/ijts
Large eddy simulation (LES) calculations of natural convection in cylindrical enclosures with rack and seeds for crystal growth applications Hooman Enayati a, Abhilash J. Chandy b, *, Minel J. Braun a a b
Department of Mechanical Engineering, University of Akron, Akron, OH, 44325, USA Department of Mechanical Engineering, Indian Institute of Technology Bombay, Mumbai, 400076, India
a r t i c l e i n f o
a b s t r a c t
Article history: Received 31 January 2017 Received in revised form 27 July 2017 Accepted 24 August 2017
This paper presents three-dimensional (3D) large eddy simulation (LES) calculations of natural convection in laterally-heated cylindrical enclosures using ANSYS FLUENT. The problem considered here is at a Rayleigh number, ðRaÞ of 8:8 106 , based on a characteristic length scale defined by the ratio of the reactor's volume to the lateral area. The numerical simulations are carried out for three different geometrical configurations: an empty reactor, a reactor with a rack and 32 seeds, and thirdly, a reactor with a rack and 80 seeds mounted inside. The objective of the current study is to understand the effects of the rack and seeds on the flow patterns and thermal distribution inside the reactor. In order to investigate these effects, a comprehensive analysis of velocity and temperature on the selected planes and points inside the domain were performed. Results showed that with an increase in seeds within the reactor, the flow slowed down more, the turbulence intensity reduced, and the extent of temperature inversion increased. These have important implications for the choice of ammonothermal method in terms of acidic or basic, and also for species transport and reactions, which eventually affect the deposition rates of GaN crystals in autoclaves. © 2017 Elsevier Masson SAS. All rights reserved.
1. Introduction Natural or free convection, one of the modes of heat transfer, occurs due to density differences, resulting from concentration or temperature differentials. This phenomena can be observed both in the engineering and natural applications. For instance as an industrial application, natural convection heat transfer plays a critical role in a crystal growth reactor. High-quality Gallium Nitride (GaN) substrates grown in such reactors can be used in light-emitting diodes (LEDs), laser diodes (LDs), UV detectors and highfrequency/high-power electronic devices and automotive applications [1]. There are several methods including hydride vapor phase epitaxy (HVPE), solution growth, sodium flux and ammonothermal techniques, that are used to produce GaN crystals. These techniques have their own advantages and disadvantages [2,3]. An ammonothermal crystal growth technique is an effective
* Corresponding author. E-mail addresses: [email protected] (H. (A.J. Chandy), [email protected] (M.J. Braun).
Enayati),
http://dx.doi.org/10.1016/j.ijthermalsci.2017.08.025 1290-0729/© 2017 Elsevier Masson SAS. All rights reserved.
[email protected]
way to grow GaN bulk single crystals with a high quality, but slow growth rate (100 mm/day). This technique uses ammonia as the solvent, and is similar to the hydrothermal technique, which is widely utilized to grow quartz crystals using water as the solvent. There have been several studies in the past utilizing both these techniques [4e7]. In the ammonothermal growth technique, the reactor has two separated thermal zones. This thermal configuration is essential in order to create a buoyant flow inside the autoclave. The preloaded GaN crystals, located in the nutrient basket, are dissolved in the dissolution zone and transported by virtue of natural convection to the crystallization zone, where GaN crystals are crystallized on seeds. As GaN crystals are almost insoluble in pure ammonia, mineralizers are added in this process. Two types of mineralizers can be used in an ammonthermal technique inside the reactor; acid and basic ([8]). In the case of a basic mineralizer such as KNH2 , the solubility of GaN crystals decreases at elevated temperatures, also called retrograde solubility. On the other hand, in the case of an acidic mineralizer such as NH4 CL, there is a direct relation between the solubility and the temperature. There is a recognized relationship between the type of mineralizer, the position of the growth crystals in the reactor, and the distribution and
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
position of the hot and cold zones. Generally, the choice of mineralizer is dependent upon the design regarding the type of crystal to be grown and the thermal distribution inside the reactor, respectively. In the present paper, the position of the rack and seed crystals replicates the architecture of a prototype of a GaN reactor. Furthermore, a baffle can be used in order to control the flow direction, and possibly create a stronger temperature gradient inside the reactor. Experimental flow visualization and measurements of temperature inside a crystal growth reactor are near impossible due to the opaque walls and harsh environments, resulting from high pressures and temperatures. Also, optimization of the growth process by performing several new sets of experiments can be expensive and time-consuming. Through computational fluid dynamics (CFD), it is possible to understand the ammonothermal crystal growth process better, at a relatively cheaper cost and with a reasonable accuracy. Finally, it is also important to realize that a proper understanding of the flow and temperature fields is necessary in order to gain knowledge about the mass transfer mechanisms, which in turn is critical in optimizing the crystal growth process. During the past few decades, several studies have been conducted to numerically investigate natural convection in an enclosure. They included two-dimensional (2D) and three-dimensional (3D) studies in various kids of geometries including rectangular enclosures [9e14]. The thermal boundary conditions in those studies were mainly side-wall or top-bottom-wall heating configurations. In addition to constant heat flux or temperature applied on the boundaries, heat generation sources inside the domain have also been employed. The current study benefits from a temperature gradient between the top and bottom halves of the lateral walls of the reactor. Two different but constant temperatures are applied laterally on the walls onto two zones, separated by an adiabatic section. Such a thermal configuration is relatively new as it is different from the more-widely used side-wall heating and topbottom-wall heating configurations. Convective heat transfer in a cylindrical enclosure with different thermal boundary conditions have been studied numerically and experimentally over the past several years as well [15e21]. Yuan et al. [22] investigated natural convective flow in a horizontal concentric annuli with different inner shapes including cylindrical, square, triangular and elliptical cross sections. The outer and inner walls were considered to be at constant temperature. They studied the effects of different inner shapes and surface radiation on the flow field and heat transfer, and validated their results with available literature. They concluded that at high temperature levels, radiation had a very important role in the overall heat transfer behavior. They also established correlations for the average Nusselt number. Malik et al. [23] performed numerical and experimental studies of axial temperature gradient and the heat transfer within an enclosure. The experimental set up included different materials and diameters for the inner and outer cylinders, respectively. The thermal boundary condition included different bottom disc temperatures ranging from 353 to 433 K. They validated their numerical results with experimental data. While several investigations have been carried out focusing on the role of the natural convection heat transfer in a crystal growth process [24e30], more work is still needed. Li et al. [31] presented a 3D numerical simulation of a cylindrical autoclave with an aspect ratio (height/diameter) of 10, subjected to circumferentially nonuniform constant temperature on the upper and lower walls of the autoclave, with and without a baffle in the middle of the enclosure. They concluded that for Ra ¼ 3:84 1011, calculated
43
based on the autoclave diameter, this nonuniform heating condition had significant effects on the fluid flow and temperature distribution in the bulk flow. For instance, temperature deviation on the upper chamber wall, resulted in a distortion of temperature isotherms in the upper growing chamber. Li et al. [32] studied a rectangular enclosure both experimentally and numerically, the latter both in 2D and 3D. The sidewalls on the lower half of the enclosure were heated, while those on the upper half were cooled uniformly. They found a good agreement between the 3D numerical results and the experimental data with respect to the thermal and fluid flow distributions. They concluded that the 2D model predictions compared qualitatively well with experiments. However, there was a 11% difference of total heat transfer rate between the 2D and 3D results. Li et al. [33], numerically studied 3D simulations of flow and heat transfer in a cylindrical hydrothermal reactor and presented them along with an experimental set up. The thermal configuration included a cooling and heating procedure on the upper and lower halves, respectively. An unsteady jet-like flow pattern in each chamber of the reactor was observed due to the fluid exchange at the baffle opening. They concluded that the growth environment was mainly affected by the total heating rate. Erlekampf et al. [34], performed a numerical analysis of an ammonothermal process for the bulk growth of nitride crystals. They explored the impacts of various baffle shapes on the flow velocity and temperature fluctuations in the domain. The current authors [35,36], used a 2D axisymmetric RANS model to perform numerical simulations of a vertical cylinder with laterally-heated walls. The transition regime from laminar flow to turbulent flow was defined based on the Ra criteria. The approach followed in the current study in terms of the definition of the Ra is similar to [36], where the characteristics length is half of the reactor's radius or R=2. Furthermore, effects of a baffle inside an ammonothermal crystal growth reactor have also been investigated using 3D large eddy simulation (LES) [37]. Results demonstrated the presence of temperature inversion for specific baffle openings in the reactor. Several numerical and experimental investigations specifically with regard to GaN crystal growth using the ammonothermal technique have also been studied [38e44]. Yoshida et al. [45] developed autoclaves in order to grow high quality GaN crystals under a certain temperature ( < 800 C) and pressure ( < 150 MPa) conditions. A retrograde solubility of GaN crystals using an acidic mineralizer ðNH4 CLÞ, was achieved for an operating pressure of 110 MPa and temperature values of 650 C and higher. They concluded that a high-temperature acidic ammonothermal technique could yield low-etch pit density (around 104 cm2) and improve the quality of the GaN crystals. Ehrentraut et al. [46] investigated the solubility of GaN in an ammonothermal method containing an acidic mineralizer ðNH4 CLÞ. It was found that NH4 CL was a good choice among the acidic mineralizers for GaN crystal growth in terms of lattice distortion and its concentration played a key role in controlling the solubility of GaN. 2. Scope of work In spite of the several studies in the past on GaN crystal growth investigations, there has not yet been a comprehensive and fundamental analysis of the flow and heat transfer patterns focusing on the effects of different geometric configurations related to the number of seeds in an ammonothermal crystal growth reactor. ANSYS FLUENT numerical simulations of turbulent natural convection flow in an empty reactor, and a reactor with 32 and 80 seeds, are presented, in order to understand the influence of
44
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
internal furniture on convective flow, thermal maps, and potentially mass transfer mechanisms, all of which affect the crystal growth rates. Temperature distributions and flow patterns on several points and planes are analyzed and compared among the cases, all for a laterally heated cylindrical enclosure at Ra ¼ 8:8 106 . The characteristic length scale used here for the Ra definition is the ratio of the volume of the enclosure to its lateral area, resulting in half of the radius or R=2, where R is the radius of the enclosure. For the first time, LES, along with the dynamic Smagorinsky model for subgrid-scale (SGS) closure is employed for simulations of an ammonothermal crystal growth reactor. In all cases, the aspect ratio (total length/diameter) is 5.17 and Pr ¼ 5:85.
In the above equations, ð…Þ represents the LES filtering operation, ui is the velocity component, p is the pressure, r is the density, g is the acceleration due to gravity, b is the coefficient of thermal expansion, T is the temperature, Tref is the reference temperature, sij is the stress tensor due to molecular viscosity and tij is the subgrid-scale (SGS) stress tensor, defined in Eqs. (3) and (4):
3. Governing equations
The SGS turbulent viscosity, nt , in Eq. (4) is defined according to the classical Smagorinsky model as [47]:
The governing equations used in this study are the incompressible filtered Navier-Stokes equations subjected to the Boussinesq approximation. The continuity, momentum and energy equations are presented in Equations (1), (2) and (7), respectively, as follows:
vðui Þ ¼0 vxi vsij 1 vp vtij vðui Þ v ui uj þ ¼ þ g b T Tref : vt vxj r vxj vxj vxj
"
vui vuj sij ¼ m þ vxj vxi
!#
1 3
tij tkk dij ¼ 2nt Sij
(3)
(4)
nt ¼ ðCs DÞ2 S
(5)
In Eq. (5), D is the filter width, Cs is the Smagorinsky constant qffiffiffiffiffiffiffiffiffiffiffiffiffi and S ¼ 2Sij Sij . The strain-rate tensor ðSij Þ for the resolved scale
(1)
is defined by,
(2)
1 vui vuj Sij ¼ þ 2 vxj vxi
! (6)
The dynamic Smagorinsky model is employed in the study here
Fig. 1. 3D geometry for the (a) empty reactor (b) reactor with a rack and 32 seeds (c) reactor with a rack and 80 seeds (d) Y Z cross-section of the reactor with rack and 80 seed (Fig. 1c).
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
45
Fig. 3. Two-dimensional cross section of Fig. 1a on X Y plane and thermal boundary conditions.
where Prt is the subgrid Prandtl number. The non-dimensional parameters namely the Rayleigh ðRaÞ, Grashof ðGrÞ and Prandtl ðPrÞ numbers are given by:
Ra ¼
gbDTL3
na
; Gr ¼
Ra n ; Pr ¼ Pr a
(8)
Here L is the characteristic length scale, DT is the temperature difference, n is the kinematic viscosity, and a is the thermal diffusivity. 4. Numerical details 4.1. Numerical domain
Fig. 2. Dimensions of the rack and seeds in Fig. 1c (in mm) (a) Front view (b) Top view.
and according to this model, Cs is calculated based on Germono's identity [48]. The energy equation is given by:
vT v ui T v þ ¼ vt vxj vxj
nt vT Prt vxj
! (7)
The three-dimensional geometry of the empty reactor, and the reactor with a rack and 32 and 80 seeds are shown in Fig. 1aed. Fig. 2a and b presents the front and top views of Fig. 1c. A schematic of the cylindrical reactor along with the thermal boundary conditions on the different sections are presented in Fig. 3. Please note that all three cases in the current study have the same thermal boundary conditions. The dimensions of the different parts of the reactor (hot wall, cold wall, insulator, the radius of the reactor, rack and seeds), used in the numerical simulations and presented in Figs. 1e3, are based on an in-house experimental apparatus. All the seeds have a square cross section with a length of 50.84 mm. By comparing the top view of Fig. 1b and c, it is clear that their geometrical dimensions are similar, but the number of the seeds and their positions are different (e.g. 32 and 80 seeds). Also, note that the seeds and the rack (Fig. 1b and c) are located below the coordinate system indicated in Fig. 1a. This can also be understood
46
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
Fig. 4. Grids on two different planes for the empty reactor with 6:4 106 elements (a) X Y plane cut in Fig. 1(a), (b) Y Z plane cut in Fig. 1(a), (c) Exaggerated view near the wall in Fig. 4b to see the elements near the wall.
Table 1 Physical properties of the working fluid. Property
Value 3
Density (kg/m ) Heat capacity (J/kg-K) Thermal conductivity (W/m-K) Viscosity (kg/m-s) Thermal expansion coefficient (1/K)
997.1 4180.9 0.61028 0.00085384 0.000257
by comparing the dimensions shown in Figs. 2a and 3. The mesh distribution for the empty reactor on two different planes is shown in Fig. 4aec. The domain includes mainly hexadehdral elements with an average skewness and orthogonal quality of 6.06 102 and 0.985, respectively. The grid density near the wall in the empty
Table 2 Positions of the radial planes/points ðX ¼ Z ¼ 0Þ from the reference coordinate system in the Y direction in m. Point
1
2
3
4
5
Y position Point Y position Point Y position
0.467 6 0.172 11 þ0.1
0.408 7 0.113 12 þ0.2
0.349 8 0.054 13 þ0.3
0.290 9 þ0.025
0.2317 10 þ0.05
reactor case is such that there are 8 elements within 4 mm (See Fig. 4c), in order to be able to capture the velocity gradient in the boundary layers near the wall.
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
47
Fig. 5. Velocity vectors colored by the instantaneous velocity magnitude in the vicinity of walls; (a) X Y plane cut (b) Exaggerated view near the mixing zone.
4.2. Boundary and initial conditions 3D transient LES calculations of natural convection in an ammonothermal crystal growth reactor at Ra ¼ 8:8 106 are presented in the current study. For all the simulations, the entire computational domain (shown in Fig. 1) was initialized by considering the working fluid, i.e. water at rest along with a uniform temperature of 300 K. Constant temperatures 303.15 K and 296.15 K are imposed on the lower and upper sections of the reactor, respectively. The lower and upper caps of the reactor, and the midregion between the hot and cold sections are considered to be adiabatic (See Fig. 3). 4.3. Numerical set up The filtered governing equations (Eqs. (1)e(7)) along with the initial and boundary conditions, provided in Section 4.2, are numerically solved by utilizing the finite volume approach in the commercial software, ANSYS FLUENT [49]. The pressure implicit with splitting of operator (PISO) algorithm is considered for the pressure-velocity coupling and the bounded central differencing scheme is employed to spatially discretize the governing equations. Temporal discretization of the governing equation was accomplished using a second order implicit formulation. The properties of the working fluid, by considering the Boussinesq approximation, are provided in Table 1, and they were used in all three simulations presented in this study. 4.4. Mesh convergence study The current authors showed that for the numerical case
involving an empty reactor (with the same boundary conditions of the current study), total mesh elements of 3:4 106 showed a converged solution [50]. This conclusion was further verified by using experimental data. Although grid points of 3:4 106 was sufficient for the current calculations, the number of elements in this study for the empty reactor was increased such that all the three cases including those with the rack and seeds have a comparable type of mesh and resolution. In this study, the total number of elements for all three cases including empty reactor, reactor with rack and 32 seeds and the reactor with a rack and 80 seeds are 6,400,000, 9,640,000 and 8,564,000, respectively. 4.5. Time convergence study A time step of 0.005 s was used in all simulations here. A time convergence study was not performed for the 3D simulations, due to the computational effort involved with conducting multiple time-steps for different geometries. For instance, the simulation for the empty reactor with 6:4 106 elements took approximately 1300 CPU hours on 28 cores of the Intel Xenon processor, with a time step of 0.005 s. However, a time convergence study for a 2D simulation of the empty reactor with a same geometry and boundary conditions was indeed performed as part of a previous study by the current authors [36]. It was concluded that for the same Rayleigh ðRa ¼ 8:8 106 Þ and Prandtl ðPr ¼ 5:85Þ numbers, a time step of 0.01 s was appropriate for a converged solution in 2D. By extrapolation, half of this time step is considered in all simulations in the current study, thereby leading to a time step of Dt ¼ 0:005 s (See Ref. [50] as well). A convergence criteria of 5 105 for the solution of the continuity equation and 1 105 for the momentum and energy equations are considered in the
48
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
Fig. 6. Normalized mean velocity contours on the selected planes for the empty reactor over a time window of 139.5 s (a) Z ¼ 28.575 mm (b) Z ¼ - 28.575 mm (c) Z ¼ 9.525 mm (d) Z ¼ - 9.525 mm; Legend is valid for all the four plots.
numerical calculations here. 5. Results and discussion Numerical analysis of flow and heat transfer for each geometrical configuration, illustrated in Fig. 1, is presented in this section. Three geometrical configurations, including an empty reactor, a reactor with a rack and 32 seeds and a reactor with a rack and 80 seeds, are computationally modeled. Normalized velocity and temperatures are discussed in detail in Sections 5.1e5.2. For the natural convection problem in this paper, there isn't a pre-specified velocity scale, and hence, velocity is normalized using unity or 1 m/ s for convenience, whereas temperature is normalized such that L T ¼ TTT , where TH and TL are the high and low temperatures, H TL
respectively. The four axial planes of 1, 2, 3 and 4 shown in Fig. 2b, are used to study the flow patterns and temperature distributions for the analysis of each geometrical case. Planes 1, 2, 3 and 4 are located at þ28.575 mm, þ9.525 mm, 28.575 mm and 9.525 mm, respectively, all away the origin in the Z-direction and on X Y planes. Each plane is located in the middle of two consecutive seeds. Thirteen planes perpendicular to the Y-direction are also considered, in order to study the time-averaged velocity magnitude
for a time window of 139.5 s. The positions of these planes are presented in Table 2. Furthermore, thirteen points on the centerlines of these planes are also chosen to study the time-averaged Yvelocity for a similar time window. Table 2 also includes the position of all these points. It is expected that the presence of obstacles like rack and seeds affects the flow patterns and the temperature gradients inside a crystal growth autoclave. Hence, it is important to analyze the flow and thermal distributions inside the three geometrical configurations of the current study.
5.1. Flow patterns Figs. 6e8 present the contours of normalized mean velocity on the selected planes, shown in Fig. 2b, for the empty reactor, reactor with a rack and 32 seeds and a rack and 80 seeds over a time window of 139.5 s. Although the normalized mean velocities are in the same range of values on those planes for all the three cases (0e0.015), the velocity contours do vary qualitatively case by case. The velocity values in the central regions of the planes are relatively higher than other regions of the reactor (like regions close to the lower and upper caps), except in the vicinity of walls, where the high velocities are observed in the boundary layer. This
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
49
Fig. 7. Normalized mean velocity contours on the selected planes for the reactor with a rack and 32 seeds over a time window of 139.5 s (a) Z ¼ 28.575 mm (b) Z ¼ - 28.575 mm (c) Z ¼ 9.525 mm (d) Z ¼ - 9.525 mm; Legend is valid for all the four plots.
central region is referred to as the mixing zone, since the hot and cold fluids meet each other, and exchange momentum and enthalpy. In fact, the lower hot fluid layers move up and the upper cold fluid layers move down, and interact with each other in this region, near the lateral walls. This also results in random fluid motions at higher velocities and vortical formations at different scales, which in turn aids in the overall mixing and circulation within the reactor. Fig. 5 shows the fluid layers near the wall and the momentum boundary layer formation in the lower hot and upper cold sections of the empty reactor. Note that the velocity magnitude in the boundary layer increases away from the walls to the free stream portion of the flow and then decreases again. As the flow is turbulent in this cylindrical geometry, the thickness of the boundary layer varies as the free stream flow is non-uniform inside the reactor.
5.1.2. Reactor with rack and 32 seeds Fig. 7 shows the normalized mean velocity contours on the same planes, but for the reactor that has 32 seeds mounted on the rack. The presence of the seeds makes the flow with higher velocity values circulate mainly near the mixing zone or in the upper part of the reactor. This is far different from the case of the empty reactor. In fact, the main difference between Figs. 6 and 7 is the flow distribution in the lower section of the reactor. The flow is slower in this region for the reactor with seeds, due to the presence of the obstacle. This occurs as the flow cannot move freely enough in order to exchange momentum within the domain. However, the flow circulation and vortices in the mixing zone can move towards the upper section of the reactor easily, and can maintain higher velocities, since there are no obstacles in that section. In addition, the higher mean velocity values exist in the vicinity of the lateral walls, rather than the interior regions of the reactor in both cases.
5.1.1. Empty reactor Fig. 6aed shows the mean velocity distribution on several planes in the empty reactor. There are similar features, like the relatively higher velocities in the mixing zone and the decrease in values as the flow moves towards the caps on the top and bottom.
5.1.3. Reactor with rack and 80 seeds Similarly, Fig. 8 shows the normalized mean velocity contours on the corresponding planes for the case with 80 seeds mounted on the rack inside the reactor. The flow distribution on each plane is similar to the case with 32 seeds (see Fig. 7). However the
50
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
Fig. 8. Normalized mean velocity contour on the selected planes for the reactor and a rack and 80 seeds over a time window of 139.5 s (a) Z ¼ 28.575 mm (b) Z ¼ - 28.575 mm (c) Z ¼ 9.525 mm (d) Z ¼ - 9.525 mm; Legend is valid for all the four plots.
maximum mean velocity here is lower (0.014 m/s), although only slightly, compared to the 32 seeds case, and overall, the velocity distribution in the lower section of the reactor is different, due to the added seeds in the reactor, which further slows down the fluid in the lower part of the reactor. 5.1.4. Additional discussion More information about the effects of different numbers of seeds on the flow can be obtained by analyzing the area-averaged velocity over a time window of 139.5 s for each case of the current study on the selected planes of Table 2. As shown in Fig. 9, the empty reactor generally has the highest velocity magnitude on the planes considered here, in comparison to the other two cases, while the case with 80 seeds has the lowest. The presence of the seeds slows down the flow in the lower section of the reactor as shown in Fig. 9. However, they do not affect the flow upstairs to such an extent in all three cases, as there are no obstacles and in general, the flow is only affected by the buoyancy force in boundary layers and in the central core region by shear-induced entrainment. Note that for each case, the highest averaged velocity magnitudes are on the planes that are located near the mixing zone, for instance, Y ¼ 0:025 m. As the location of a plane moves towards the top and
bottom of the reactor, the velocity magnitude decreases on the corresponding plane. This is expected due to momentum exchange in the mixing zone, where the hot and cold fluid parcels interact with each other, and the chaotic movements and vortices ensue. As a result, mean velocity values on the planes near the mixing zones are higher. Furthermore, Fig. 10 shows the time-averaged Y-velocity component over a time window of 139.5 s for the thirteen points listed in Table 2, for each geometrical configuration. As shown in the figure, the presence of the seeds (32 or 80) reduces the Y-velocity values in the lower section of the reactor, as indicated by the values on Points # 1e8. This behavior is expected as the presence of an obstacle (seeds in this case) in the domain slows down the motions. The direction of the Y-velocity component in the lower section of the reactor also changes with the presence of the seeds. This is important from the perspective of species transport and reactions involving GaN crystals in the ammonothermal technique. Knowledge of fluid flow directions and magnitudes would help in designing an appropriate geometrical configuration for better deposition rates on the seeds inside the crystal growth autoclave. Fig. 10 also shows that the direction of Y-velocity component in the upper part of the reactor does not change with the geometrical
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
Fig. 9. Time-averaged velocity magnitudes on different planes over a time window of 139.5 s for three geometrical configurations.
Fig. 10. Time-averaged Y-velocities for different points over a time window of 139.5 s for three geometrical configurations.
51
52
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
Fig. 11. Resolved turbulent kinetic energy (TKE) on Z ¼ 28.575 mm for the (a) empty reactor (b) eeactor with the rack and 32 seeds (c) eeactor with the rack and 32 seeds.
configuration, and overall, higher Y-velocity values can be expected in comparison to the lower section of the reactor. Another interesting feature to analyze in such problems is the amount of turbulence within each configuration, since the latter will have an effect on the mixing and circulation and therefore, the etching and deposition rates. Turbulent kinetic energy (TKE) is used for that purpose, and the values of TKE and its distribution in each geometry are studied here. TKE is the mean kinetic energy per unit mass for the associated eddies in a turbulent flow. It is characterized by measured root-mean-square ðRMSÞ of the velocity fluctuations in a computational domain. For the current LES study, the resolved turbulence kinetic energy is quantified as follows:
TKE ¼ 0:5 RMS2x þ RMS2y þ RMS2z
(9)
where RMS is the root-mean-square of each velocity component. RMSy for instance can be defined as
RMSy ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðv 〈v〉Þ2
where v and 〈v〉 represent the instantaneous and mean Y-velocity components, respectively. Fig. 11 presents the TKE on a selected plane of Z ¼ 28:575 m for three different geometries. As is clear from the figure, the presence of the seeds decreases the maximum value of the TKE as a result of the decreased velocities due to the obstacles. The empty reactor has the highest TKE value, and in the regions around the mixing zone, a wider distribution of relatively higher values of TKE can be observed. Fig. 11, also shows that at the
current Ra, the flow is in the low-turbulent regime, since in all the three cases, the TKE is very low in the domain except in the mixing zone. 5.2. Temperature distribution As discussed earlier in Section 1, since GaN is insoluble in pure ammonia, mineralizers are added to strengthen the solubility of GaN crystals in an ammonothermal crystal growth reactor. Two main types of mineralizers, basic and acidic, can be chosen based on the temperature distribution inside the reactor. The thermal distribution inside an ammonothermal crystal growth reactor is a determining factor for the type of the mineralizer used in the nutrient basket that is filled with GaN crystals. For instance, for a retrograde solubility, the nutrient basket should be kept in a cooler zone inside the crystal growth autoclave. In fact, it is also important to know how the temperature distribution would develop inside the reactor as it could affect the etching and deposition rates of GaN crystals in the reactor. 5.2.1. Empty reactor Fig. 12 presents the velocity vectors superimposed on the normalized mean temperature distribution, averaged over 139.5 s, inside the empty reactor. The fluid layers adjacent to the lower lateral hot walls of the reactor are warm. These layers travel upward due to momentum and buoyancy forces, cross the mixing zone and fill the upper core of the reactor with a warm fluid. The cooler fluid layers in the vicinity of upper cold lateral walls travel down and fill up the lower section of the autoclave with the cooler fluid. This
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
53
Fig. 12. Normalized mean temperature contours on the selected planes for the empty reactor over a time window of 139.5 s (a) Z ¼ 28.575 mm (b) Z ¼ - 28.575 mm (c) Z ¼ 9.525 mm (d) Z ¼ - 9.525 mm; Legend is valid for all the four plots.
eventually results in a cooler fluid in the lower section of the reactor and hotter fluid in the upper section. This phenomena is like a temperature inversion, as the warmer flow exists in the upper section of the reactor, while the walls are colder. Similarly, the cooler fluid exists in the lower section of the reactor, where the walls are hotter. Note that in the problem studied here, as the temperature distribution is almost uniform away the walls, the buoyancy forces resulting from the temperature gradients inside the reactor, are weak. However, the boundary layers move up and down along the wall mainly due to the thermal boundary conditions, while the bulk fluid within the core of the reactor move mainly due to entrainment and shear. 5.2.2. Reactor with rack and 32 seeds Fig. 13 shows the velocity vectors superimposed over the normalized mean temperature over the same time window and for the same planes as presented earlier, except this time, for the case with 32 seeds mounted on a rack inside the rector. A discussion of the temperature distribution similar to the empty reactor is valid here too. However, the presence of 32 seeds causes a larger temperature gradient, especially away from the lateral walls, between the lower and upper sections of the reactor. Here, the upper section
of the reactor is slightly warmer in comparison to the empty reactor case. This happens because the rack and seeds in this case obstructs the flow significantly in the lower section of the reactor. So while the warmer fluid is indeed transported up towards the upper section, the colder fluid does not move as freely towards the lower section (compared to the empty reactor) due to the presence of the rack and seeds. This basically results in a lower overall circulation in the reactor, which leads to a less uniform temperature distribution in the core of the reactor, compared to the empty reactor case. 5.2.3. Reactor with rack and 80 seeds Fig. 14 illustrates the velocity vectors superimposed on the normalized mean temperature contours over a time window of 139.5 s, again for the same planes, but for the case with 80 seeds mounted on the rack inside the rector. A similar temperature inversion behavior is observed here too. The current geometrical configuration though provides the warmest upper zone of the reactor among all the three cases. Furthermore, unlike the case with 32 seeds mounted on the rack inside the reactor, the current geometrical configuration of the (80) seeds provides a slightly nonuniform temperature in the lower section of the reactor. Note that the region close to the walls on the lower section of the reactor has
54
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
Fig. 13. Normalized mean temperature contour on the selected planes for the reactor with a rack and 32 seeds over a time window of 139.5 s (a) Z ¼ 28.575 mm (b) Z ¼ 28.575 mm (c) Z ¼ 9.525 mm (d) Z ¼ - 9.525 mm; Legend is valid for all the four plots.
slightly higher temperatures. 5.2.4. Additional discussion As was earlier illustrated through Figs. 7 and 8, the presence of the seeds decreases the velocity magnitude in the lower section of the reactor. The velocity vectors in Figs. 12e14 shows a similar trend as well. The spatial density of these velocity vectors in the domain qualitatively represents the velocity magnitude. The density of the velocity vectors is higher near the walls due to the boundary layer formation for all the three cases. While the density of the velocity vectors in the empty reactor is comparable in the upper and lower sections of the reactor, it is not the case for seeds within the reactor. For instance, here, the velocity vectors are denser in the mixing zone and in the upper section of the reactor, compared to the lower section. This further confirms the existence of a slower flow inside the reactor for the cases with seeds. In addition, it can be seen from Figs. 12e14 that the velocity vectors in the bulk flow of the lower hot section are downward while in the upper cold section of the reactor are upward. The thickness of the thermal boundary layer and normalized mean temperature profiles across six radial lines on X Y plane (See Fig. 3) are shown in Fig. 15. The position of each line is tabulated in Table 3. It is clear from Fig. 15a and b that the temperature distribution is uniform far away the walls on the selected lines for
all the three cases of the study, as was mentioned earlier in the discussion of the temperature contours. The thermal boundary layer formation, where a sharp temperature gradient between the bulk and near wall temperatures exists, can also be noted in these figures. The thickness of the thermal boundary layer increases closer to the mixing zone, while away the mixing zone, the thermal boundary layer is smaller. As discussed earlier for Fig. 5, in the lower section of the reactor, the flow in the vicinity of walls is upwards, while in the upper section of the reactor, the flow is downwards. This results in a boundary layer development from bottom to top, and vise-versa in the lower and upper sections of the reactor, respectively. This would also further confirm that the thermal boundary thickness keeps developing up until the mixing zone. 6. Conclusion A study of natural convection in a laterally-heated cylindrical enclosure was performed by conducting a series of numerical simulations using the commercial CFD software, ANSYS FLUENT. LES calculations of three cases with different geometrical configurations, including an empty reactor, a reactor with a rack and 32 seeds and a reactor with a rack and 80 seeds, were analyzed. Distributions of velocity magnitude, temperature and resolved
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
55
Fig. 14. Normalized mean temperature contours on the selected planes for the reactor with a rack and 80 seeds over a time window of 139.5 s (a) Z ¼ 28.575 mm (b) Z ¼ 28.575 mm (c) Z ¼ 9.525 mm (d) Z ¼ - 9.525 mm; Legend is valid for all the four plots.
Table 3 Positions of the lines (in m) from the reference coordinate system in the Y direction (Z ¼ 0, 0.079375 m < x < 0.079375 m). Point
1
2
3
4
5
6
Y position
0.349
0.231
0.113
þ0.1
þ0.2
þ0.3
turbulent kinetic energy, along with velocity vectors, were presented for each case on selected axial and radial planes located in the domain. All the three cases of study have the same boundary conditions with a Ra ¼ 8:8 106 , based on the characteristic length scale of R2, where R is the radius of the cylinder. In the mixing zone, where the upper cold and lower hot boundary layers meet each other, the fluid exchanges momentum and heat, resulting in mixing and circulation. Vortices of different scales that are formed in and around this region, evolve, disappear and eventually reappear. Furthermore, this mixing zone exhibits the highest velocities in the domain, away from the walls. The presence of the rack and seeds in the lower part of the reactor, decreased the velocity in these regions, due to the added obstructions. The maximum velocity within the domain stayed the same, since it was a function of the velocity in the boundary layer, which in turn solely depends on the thermal boundary conditions that are identical for all the configurations. It was also shown that
the direction and values of the Y-velocity changed dramatically inside the reactor with the added furniture. For instance, the Yvelocity values in the lower section was highest for the empty reactor among all the geometries. It was also noted that the direction of the Y-velocity has potential implications on the deposition rates of GaN crystals inside the reactor, since a favorable flow towards the seeds can bring more GaN crystals from the nutrient basket. This further affirms the fact that a sufficient knowledge of the fluid flow is critical in order to understand species transport and reactions in a crystal growth reactor. In addition, TKE distributions indicated a reduction in turbulence as well with the addition of rack and increase in the number of seeds within the reactor. Due to the thin wall boundary layers that exist for the current problem at a relatively high- Ra, the thermal resistance between the bulk flow and the boundary layer region is weak. This also implies relatively uniform temperature values far away the walls in the reactor. The current thermal distribution inside the domain, results in a weaker flow circulation as the buoyancy force is weak due to the absence of strong temperature gradients. However, a strong buoyancy force does exist near the walls due to the thermal boundary conditions, which is the primary driver of the flow circulation in the reactor. The boundary layers entrain the flow in the core of the reactor via shear, thus creating some movements there as well. In all the three configurations studied here, a temperature
Fig. 15. Normalized mean temperatures on the selected lines for three geometrical configurations over a time window of 139.5 s (a) rhree lines in the lower section of the reactor (b) rhree lines in the upper section of the reactor.
H. Enayati et al. / International Journal of Thermal Sciences 123 (2018) 42e57
inversion was also observed. This phenomenon is evident from the existence of warmer fluid in the upper cold section of the reactor and a colder fluid in the lower hot section of the reactor. This kind of temperature inversion can determine the type of a mineralizer (acidic or basic) in the nutrient basket of a crystal growth reactor. Finally, it was also shown that the thermal boundary layer was the thickest in the mixing zone. Acknowledgments This paper is based upon work supported by the National Science Foundation under Grant No. 1336700. References [1] Kachi Tetsu. Recent progress of gan power devices for automotive applications. Jpn J Appl Phys 2014;53(10):100210. [2] Ehrentraut Dirk, Meissner Elke, Bockowski Michal. Technology of gallium nitride crystal growth, vol. 133. Springer Science & Business Media; 2010. [3] Mamun Molla Md, Paul Sreebash C, Hossain Md Anwar. Natural convection flow from a horizontal circular cylinder with uniform heat flux in presence of heat generation. Appl Math Model 2009;33(7):3226e36. [4] Roux B, Louchart O, Terhmina O. Hydrodynamic aspect of hydrothermal synthesis of quartz bulk flow regimes. Le J de Physique IV 1994;4(C2):C2e3. ski R, Wysmołek A, Baranowski J, Kamin ska M, Doradzin ski R, [5] 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, [6] Dwilin Sierzputowski L. On gan crystallization by ammonothermal method. Acta Phys Pol A 1996;90(4):763e6. [7] Purdy Andrew P. Ammonothermal synthesis of cubic gallium nitride. Chem Mater 1999;11(7):1648e51. [8] Ehrentraut Dirk, Hoshino Naruhiro, Kagamitani Yuji, Yoshikawa Akira, Fukuda Tsuguo, Itoh Hirohisa, et al. Temperature effect of ammonium halogenides as mineralizers on the phase stability of gallium nitride synthesized under acidic ammonothermal conditions. J Mater Chem 2007;17(9):886e93. [9] Chen Q-S, Prasad V, Chatterjee A. Modeling of fluid flow and heat transfer in a hydrothermal crystal growth system: use of fluid-superposed porous layer theory. J heat Transf 1999;121(4):1049e58. [10] 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. [11] 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. [12] 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. [13] Kahwaji Ghalib, Ali Omar M. Numerical investigation of natural convection heat transfer from square cylinder in an enclosed enclosure filled with nanofluids. Int J Recent Adv Mech Eng (IJMECH) 2015;4(4). [14] Iyi Draco, Hasan Reaz. Natural convection flow and heat transfer in an enclosure containing staggered arrangement of blockages. Procedia Eng 2015;105:176e83. [15] Hsieh Shou-Shing. Thermal correlation of natural convection in bottomcooled cylindrical enclosures. J Thermophys Heat Transf 1990;4(1):123e6. [16] Li YH, Lin KC, Lin TF. Computation of unstable liquid metal convection in a vertical closed cylinder heated from the side and cooled from above. Numer Heat Transf Part A Appl 1997;32(3):289e309. [17] Choukairy Kh, Bennacer R, Beji H, Jaballah S, El Ganaoui M. Transient behavior inside a vertical cylindrical enclosure heated from the side walls. Numer Heat Transf Part A Appl 2006;50(8):773e85. [18] Barhaghi Darioush G. A study of turbulent natural convection boundary layers using large-eddy simulation. Chalmers University of Technology, Department of Applied Mechanics; 2007. [19] Malik Asif Hussain, Shah Ajmal, Khushnood Shahab. CFD analysis of heat transfer within a bottom heated vertical concentric cylindrical enclosure. J Phys: Conf Ser 2013;439:012004. [20] 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. [21] Moukalled F, Darwish M, Kasamani J, Hammoud A, Khamis Mansour M. Buoyancy-induced flow and heat transfer in a porous annulus between concentric horizontal circular and square cylinders. Numer Heat Transf Part A Appl 2016;69(9):1029e50. [22] Yuan Xing, Tavakkoli Fatemeh, Vafai Kambiz. Analysis of natural convection in horizontal concentric annuli of varying inner shape. Numer Heat Transf Part A
57
Appl 2015;68(11):1155e74. [23] Asif Hussain Malik, Shahab Khushnood, and Ajmal Shah. Experimental and numerical study of buoyancy driven flow within a bottom heated vertical concentric cylindrical enclosure. 2013. [24] Chen Q-S, Pendurti S, Prasad V. Effects of baffle design on fluid flow and heat transfer in ammonothermal growth of nitrides. J Cryst growth 2004;266(1): 271e7. [25] Li Hongmin, Braun Minel J. Numerical investigation on multi-hole baffle designs for industry hydrothermal synthesis of single crystals. Model Simul Mater Sci Eng 2005;13(8):1249. [26] Chen Q-S, Pendurti S, Prasad V. Modeling of ammonothermal growth of gallium nitride single crystals. J Mater Sci 2006;41(5):1409e14. [27] 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. [28] Jiang Yan-Ni, Chen Qi-Sheng, Prasad V. Numerical simulation of ammonothermal growth processes of gan crystals. J Cryst Growth 2011;318(1):411e4. [29] 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. [30] 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. [31] Li Hongmin, Wang Guo-Xiang, Evans Edward A. Three-dimensional flow of solution in an industry-size hydrothermal autoclave subjected to nonuniform heating effects of a baffle on flow and temperature separation. J Cryst growth 2004;271(1):257e67. [32] 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. [33] 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. [34] 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. [35] Enayati Hooman, Chandy Abhilash J, Braun Minel J. Numerical simulations of natural convection in a laterally-heated cylindrical reactor, vol. 10. NYC, DOI: American Society of Thermal and Fluids Engineers; 2015. [36] 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. [37] H. Enayati, A. J. Chandy, and M. J. Braun. Numerical simulations of turbulent natural convection in laterally-heated cylindrical enclosures with baffles for crystal growth. In IMECE 2016, Phoenix, Nov-2016. [38] Popov VN, S Tsivinskaya Yu, Bekker TB, Kokh KA, Kokh AE. Numerical investigation of heat-mass transfer processes in hydrothermal growth system. J Cryst growth 2006;289(2):652e8. [39] 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. [40] Ehrentraut Dirk, Fukuda Tsuguo. The ammonothermal crystal growth of gallium nitride, a technique on the up rise. Proc IEEE 2010;98(7):1316e23. [41] Ehrentraut Dirk, Pakalapati Rajeev T, Kamber Derrick S, Jiang Wenkan, Pocius Douglas W, Downey Bradley C, et al. High quality, low cost ammonothermal bulk gan substrates. Jpn J Appl Phys 2013;52(8S):08JA01. [42] Richter Theresia MM, Niewa Rainer. Chemistry of ammonothermal synthesis. Inorganics 2014;2(1):29e78. [43] 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. [44] 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. [45] Yoshida Kazuo, Aoki Kensuke, Fukuda Tsuguo. High-temperature acidic ammonothermal method for gan crystal growth. J Cryst Growth 2014;393: 93e7. [46] Ehrentraut Dirk, Kagamitani Yuji, Yokoyama Chiaki, Fukuda Tsuguo. Physicochemical features of the acid ammonothermal growth of gan. J Cryst Growth 2008;310(5):891e5. [47] Smagorinsky J. General circulation experiments with the primitive equations. I. The basic experiment. Mon Weather Rev 1963;91(3):99e164. [48] Sagaut P. Large eddy simulations of incompressible flows. Berlin: SpringerVerlag; 2005. [49] Ansys Co. Ansys fluent users guide. 2014. Version 15.0. Ansys. Ansys. [50] Enayati Hooman, Chandy Abhilash J, Braun Minel J, Horning Nicholas. 3d large eddy simulation (les) calculations and experiments of natural convection in a laterally-heated cylindrical enclosure for crystal growth. Int J Therm Sci 2017;116:1e21.