Numerical Simulation of Temperature Fields in a Three-Dimensional SiC Crystal Growth Furnace with Axisymmetric and Spiral Coils

: Three-dimensional numerical simulation platform for silicon carbide crystal growth furnace was established using C programing language, where a physical model of the furnace was built based on cylindrical coordinates; governing equations for electromagnetic and temperature ﬁelds were discretized by ﬁnite volume method; radiation characteristics were studied with the help of S2S model (surface to surface radiation model); and the least distance method was proposed to check radiation faces visibility efﬁciently. LU decomposition algorithm based on graphic processing unit (GPU) technology was developed to accelerate the solving process of surface to surface radiation. Then the radiation heat transfer in silicon carbide crystal (SiC) growth chamber and temperature ﬁeld of silicon carbide growth furnace were studied quantiﬁcationally at I = 1250 A and F = 16 kHz. The effects of coil structures (axisymmetric and spiral) on temperature ﬁeld and its gradient distributions were investigated by standard deviation method. The simulation results demonstrate that spiral electromagnetic coil generates non-axisymmetric temperature ﬁeld easily; the radiation heat ﬂux is 10 2 ~10 3 times more than conduction heat ﬂux, radiation heat transfer is helpful to increase temperature evenness; the spiral temperature ﬁeld on the SiC crystal cross-section reduces the poor homogeneity of temperature gradient, which will cause crystal to generate large defects.


Introduction
Silicon carbide crystal (SiC) crystals are used to product semiconductor devices under higher-temperature, higher-frequency, and higher-power conditions [1]. Compared with silicon, silicon carbide crystal has better performance on chemical and thermal properties, such as higher thermal conductivity, chemical stability, radiation stability, and mechanical stability [2]. The physical vapor transport (PVT) method is the most successful and common method to grow SiC crystals at pressure ≤ 50 mbar, increasing crystal growth rate [3][4][5]. In this method, graphite crucible is heated by induction coils and the temperature in the furnace reaches more than 2300 K. The temperature distribution is fatal to SiC growth process, which decides the SiC growth speed and the product quality. However, it is difficult to measure temperature in a confined silicon carbide crystal growth furnace at both high temperature and low pressure, and so far, the infrared thermometry method is adopted to monitor temperature through a gap cooler hole. The detailed study of temperature distribution in the furnace is beneficial to the modification of the furnace structure and operation system.
With the development of computer technology, numerical simulation is becoming a powerful tool to analyze temperature fields [3,[6][7][8][9], crystal growth processes [10][11][12], and thermal stresses [13] as the SiC grows. These parameters are meaningful to optimize the process of crystal growth. Both conduction and radiation are the main two heat transfer processes in 2~4 inches SiC crystal growth furnaces. The weaker convection heat transfer is always ignored. For radiation heat transfer, species in the growth chamber are not participated in radiation process, so the surface to surface (S2S) radiation model is used to calculate radiation heat flux [2,[14][15][16]. The axisymmetric assumption is adopted in the majority of numerical simulation research, so the model is two-dimensional, the effect of non-axisymmetric coils on electromagnetic field and temperature field are always missed artificially. Moreover, some graphite crucible surfaces cannot see other surfaces because of the nontransparent crystal. Therefore, some methods should be developed to check whether radiation surfaces are visible to other surfaces or not, ensuring the accurate calculation of radiation heat flux.
Herein, a three-dimensional numerical simulation algorithm was developed based on cylindrical coordinate system and C programing language to study the effects of axisymmetric model and spiral coils on electromagnetic and temperature fields at coil current I = 1250 A and current frequency F = 16 kHz. The least distance method is proposed to check visibility of radiation facing to others, and a parallel LU decomposition algorithm for radiation view factor matrix based on graphic processing unit technology is developed to accelerate the solution of surface to surface radiation heat transfer.

Modeling of Induction Heating and Heat Transfer
Radio frequency (RF) induction heating technology is used to heat the silicon carbide crystal growth furnace. The reasonable temperature field is premised on an accurate electromagnetic field. Generally, mathematical models for the electromagnetic field are developed based on the Maxwell where E and B are electric field intensity and magnetic induction intensity, respectively; µ 0 (4π × 10 −7 H/m) is permeability of free space; ε 0 (8.855 × 10 −12 F/m in vacuum) is permittivity of free space; J c is current flux in coil; and J i (A/m 2 ) is induction current flux in graphite crucible. To simplify Maxwell equations, magnetic potential, A (B = ∇ × A), is introduced and the simplified equation is proposed in 3D cylindrical coordinates 1 r ∂ ∂r r ∂A ∂r A can be expressed by (0, A 0 , 0), and A 0 = A 0r + iA 0i , so Equation (2) can be rewrote as [7,17] where ω is angular velocity (rad/s, ω = 2πF) and F is current frequency (Hz); J c is current flux in coil (A/m 2 , J c = I/A coil ) and I is current (A), A coil is area of coil cross section; σ (1.1 × 10 5 Ω −1 ·m −1 ) is coil electric conductivity.
For the different subsystems-such as crystals, silicon carbide powder, graphite crucible, insulation layer, and air of outside furnace-the heat transfer process can be described as the following energy Equation (4) where Q eddy is induction heating source term in the graphite crucible; Q IF is heating source term acting on cell interfaces, such as radiation and convection heat exchange (Q insu ) on felt walls, radiation heat exchange (Q rad ) on growth chamber walls.
Equations (1) and (2) are discretized by the finite volume method, and solved by the SIMPLE algorithm.
Q eddy as the induction heating source is calculated by [2,18] The third boundary condition is used to describe heat release on the felt walls, and the heating source can be added into Equation (2) with the additional source terms method where n is vector normal to outer wall. h is convection heat transfer coefficient between the insulation layer and air, σ is the Boltzmann constant (5.67 × 10 −8 W/(m 2 ·K 4 ), e is surface radiation coefficient, h eff is the equivalent heat transfer coefficient of convection and radiation heat transfer, d is the grid length in the radial direction, M is the circumference grid area of the volume mesh on the insulation layer, ∆V is the grid volume, and T f is the fluid temperature. S C,ad,insu and S P,ad,insu are source terms for solving algebraic equations of wall temperatures. Figure 1 shows a three-dimensional geometry and meshes for silicon carbide crystal growth furnace. To decrease calculation capacity, the electromagnetic and temperature fields are decoupled. Their physical models are presented in Figure 1a,b, respectively. Figure 1c indicates the mesh for physical model in Figure 1a. As shown in the figure, more fine meshes are used for silicon carbide crystal, silicon carbide powder and graphite crucible to get more accurate simulation results. It can be seen from Figure 1a that the coil is quasi-sprial, and the size of its cross section is 10 × 10 mm. The silicon carbide crystal is 50 mm in diameter and 5 mm in height. Graphite crucible is 125 mm in inner diameter, 6 mm thick in r-direction, and 20 mm thick in z-direction. Silicon carbide powder and growth chamber are 55 mm and 25 mm in height, respectively. The furnace is cooled by both top hole with diameter of 20 mm and bottom hole with diameter of 10 mm in diameter. The convection heat transfer coefficient of cooling water (T = 300 K) is 1000 W/(m 2 ·K). The felt is 37 mm thick in r-direction and 40 mm thick in z-direction. The convection heat transfer coefficient of air outside of the furnace (T = 300 K) is 50 W/(m 2 ·K). Table 1 lists properties of several materials, such as density, heat capacity, heat conductivity coefficient. Properties of porous silicon carbide powder are calculated by ϕ SiC_powder = εϕ SiC + (1 − ε)ϕ fluid , where ε is void fraction and equals to 0.45.  Table 1 lists properties of several materials, such as density, heat capacity, heat conductivity coefficient. Properties of porous silicon carbide powder are calculated by φSiC_powder = εφSiC + (1 − ε)φfluid, where ε is void fraction and equals to 0.45.

Surface Visibility Checking Algorithm for S2S Radiation Model
The classical surface to surface (S2S) radiation model is where A and B are coefficient matrixes of view factors (F); when j = k, δjk = 1, otherwise δjk = 0; T is wall temperature; ηjk = 1 in Fij means face j can see face k, otherwise ηjk = 0. So ηjk should be determined before calculating Fij. The radiation source term, Qrad, is added into energy equation by qr as the following forms [13] :

Surface Visibility Checking Algorithm for S2S Radiation Model
The classical surface to surface (S2S) radiation model is where A and B are coefficient matrixes of view factors (F); when j = k, δ jk = 1, otherwise δ jk = 0; T is wall temperature; η jk = 1 in F ij means face j can see face k, otherwise η jk = 0. So η jk should be determined before calculating F ij .
The radiation source term, Q rad , is added into energy equation by q r as the following forms [13] where δ sf and δ ff are distances from the center of solid mesh cell and fluid mesh cell to the fluid-solid interface center. Least distance method (LDM) is proposed to determine the value of η jk between two wall faces. interface center.
Least distance method (LDM) is proposed to determine the value of ηjk between two wall faces. Figure 2 exhibits a schematic diagram and a flowchart of LDM. An example (ith cell face) is given to illustrate the principle of LDM: (1) give the mathematical forms (F(R)) of block faces, such as growth surface and cylindrical surface of silicon carbide crystal; (2) search all other jth cell faces and save their face information as shown in Figure 2a, such as normal vector n, distance vector d, position coordinates of face center, line equation L(i, j) from ith face center to jth face center, angle θi and θj between n and d; (3) calculate and save all intersection points (P(L, F)) between line L(i, j) and face F(R); (4) calculate all distances from ith face center to P(L, F) and record the point owning a minimum distance; (5)  Two wall faces (located on SiC powder surface 1# mesh face and graphite crucible inner surface (2# mesh face)) are selected to verify LDM as shown in Figure 3. The red regions as shown in Figure  3a,c are ideal visible region for 1# mesh face and 2# mesh face, respectively. By comparing of Figure  3a,b, it can be seen that calculated visible region is as same as ideal visible region which locates on graphite cylindrical wall and crystal growth wall. Similar results are shown in Figure 3c,d. These simulation results can be demonstrated that LDM is an efficient and accurate method to check whether radiation faces can see other faces or not.

Acceleration to Solve S2S Radiation Heat
Matrix A and B in Equation (7) are only relevant with radiation view factors which are determined by physical geometry and meshes. With the increase of mesh number in growth chamber, the calculation capacity increases sharply, not only do matrix A and B need more memory, but also

Acceleration to Solve S2S Radiation Heat
Matrix A and B in Equation (7) are only relevant with radiation view factors which are determined by physical geometry and meshes. With the increase of mesh number in growth chamber, the calculation capacity increases sharply, not only do matrix A and B need more memory, but also it will take more time for a CPU (central processing unit) to calculate q r . For purpose of solving this problem, constant matrix C is introduced to simplify Equation (7), and matrix C is solved by matrix A and B q r = C·(σT 4 ) AC = B Here, the matrix C was solved by LU decomposition algorithm (which was an algorithm to solve linear algebraic equations) based on parallel GPU codes, and the models of Equations (3)-(8) were solved by a serials of CPU codes. All the simulation works were finished on a workstation with 8 CPUs and 32G memory. The GPU codes were run on a GeForce GTX 580: NVIDIA, Santa Clara, CA, USA graphics card, which owns 512 stream processors and 1536 MB video memory.
For the same mesh and the same physical model, the pre-solved and pre-saved matrix C can be read into the program which will save much time to solve Equation (7). To further improve calculation efficiency, a parallel LU decomposition algorithm based on GPU technology is used. Four cases are set up in order to compare calculation efficiency by CPU and GPU with different mesh numbers, and the time consuming results are listed as shown in Table 2. It can be found that GPU costs less time to get the matrix C than CPU. GPU technology also indicated the faster calculation efficiency with the increase of mesh size. In addition, the temperature difference (∆T) between crystal surface and SiC powder surface is compared with different mesh size. It can be seen that ∆T for the four cases are close to each other, so 56 × 43 × 71 is chosen to get accurate temperature field and decrease calculation capacity.

Test of Mesh Independence
In addition, four mesh sizes listed in Table 3 are compared with each other in order to investigate the relationship between mesh size and both calculation accuracy and calculation efficiency. 2304, 4320, 5184, and 6912 radiation faces are generated in the growth chamber for the four cases, respectively.  Figure 4 shows profiles of magnetic potential and temperature along with radius. For A 0r and A 0i presented in Figure 4a,b, the results of case 2, 3 and 4 are in good agreement with each other, which illustrates that meshes for electromagnetic field are fine enough. As seen in Figure 4c, the similar temperature profiles are presented for case 1, 2 and 3, temperature profile of case 4 is slightly lower than those of case 1, 2 and 3. This is because more radiation surfaces reduce to bigger calculation error.
Based on the above analysis, 56 × 64 × 131 for electromagnetic field and 56 × 43 × 71 for temperature field are selected. 4 74 × 64 × 131 74 × 43 × 71 Figure 4 shows profiles of magnetic potential and temperature along with radius. For A0r and A0i presented in Figure 4a,b, the results of case 2, 3 and 4 are in good agreement with each other, which illustrates that meshes for electromagnetic field are fine enough. As seen in Figure 4c, the similar temperature profiles are presented for case 1, 2 and 3, temperature profile of case 4 is slightly lower than those of case 1, 2 and 3. This is because more radiation surfaces reduce to bigger calculation error. Based on the above analysis, 56 × 64 × 131 for electromagnetic field and 56 × 43 × 71 for temperature field are selected.    Figure 5 demonstrates the contours of A 0r of rz cross section and the iso-surface for axisymmetric coil and spiral coil at I = 1250 A and F = 16 kHz. Axisymmetric A 0r are illustrated in Figure 5a for the axisymmetric coil, and it can be seen that it is meaningful 2D axisymmetric model substitute for 3D model for axisymmetric coils. However, it also can be found obviously that A 0r is non-axisymmetric as shown in Figure 5b. In addition, as presented in Figure 6, the profile of A 0r and A 0i in r-direction for axisymmetric and spiral coil also can prove the above viewpoint. Therefore, it demonstrates that the coil structure has a great effect on the electromagnetic field.  Figure 5 demonstrates the contours of A0r of rz cross section and the iso-surface for axisymmetric coil and spiral coil at I = 1250 A and F = 16 kHz. Axisymmetric A0r are illustrated in Figure 5a for the axisymmetric coil, and it can be seen that it is meaningful 2D axisymmetric model substitute for 3D model for axisymmetric coils. However, it also can be found obviously that A0r is non-axisymmetric as shown in Figure 5b. In addition, as presented in Figure 6, the profile of A0r and A0i in r-direction for axisymmetric and spiral coil also can prove the above viewpoint. Therefore, it demonstrates that the coil structure has a great effect on the electromagnetic field. In SiC growth chamber, SiC powder surface, graphite crucible inner surface, and SiC crystal surface radiate strongly. To study the relationship between radiation and conduction heat transfer quantitatively, Figure 7 presents the radiation heat flux (qr) and conduction heat flux (qcz) of SiC powder surface and SiC crystal surface along with radius. It can be seen that: (1) both qr and qcz on SiC powder surface are bigger than 0 as shown in Figure 7a, which illustrates that the surface owns higher temperature and release heat; (2) both qr and qcz on SiC crystal surface are less than 0 as shown in Figure 7b, which illustrates that the surface owns lower temperature and absorption heat;  In SiC growth chamber, SiC powder surface, graphite crucible inner surface, and SiC crystal surface radiate strongly. To study the relationship between radiation and conduction heat transfer quantitatively, Figure 7 presents the radiation heat flux (q r ) and conduction heat flux (q cz ) of SiC powder surface and SiC crystal surface along with radius. It can be seen that: (1) both q r and q cz on SiC powder surface are bigger than 0 as shown in Figure 7a, which illustrates that the surface owns higher temperature and release heat; (2) both q r and q cz on SiC crystal surface are less than 0 as shown in Figure 7b, which illustrates that the surface owns lower temperature and absorption heat; (3) the temperature in the furnace is non-axisymmetric which can be demonstrated by the non-axisymmetric profile of q r and q cz on both surfaces; (4) q r = 160~260 KJ/m 2 and q cz = 240~390 J/m 2 on powder surface, q r = −160~−80 KJ/m 2 and q cz = −270~−240 J/m 2 on crystal surface, so q r is about 2~3 orders of magnitude greater than q cz .  Figure 8 shows profiles of temperature of powder surface and crystal surface along with r axis with and without radiation heat transfer, respectively. It can be found that: (1) the average temperature at −0.9 < r/R < 0.9 is 3247.9 K without radiation and it decreases to 2651.7 K with radiation on powder surface; (2) the average temperature is 2212.6 K without radiation and it increases to 2555.3 K with radiation on powder surface; (3) the temperature grandient is 373.1 K/cm without radiation, while it decreases to 34.76 K/cm, which in the suitable temperature ranges (1-50 K/cm) for SiC crystal growth introduced in reference [14]. The above simulation results illustrate that the reasonable temperature difference can be gotten if radiation heat transfer is considered [7] and the radiation is helpful to uniform the temperature field in the chamber.   Figure 8 shows profiles of temperature of powder surface and crystal surface along with r axis with and without radiation heat transfer, respectively. It can be found that: (1) the average temperature at −0.9 < r/R < 0.9 is 3247.9 K without radiation and it decreases to 2651.7 K with radiation on powder surface; (2) the average temperature is 2212.6 K without radiation and it increases to 2555.3 K with radiation on powder surface; (3) the temperature grandient is 373.1 K/cm without radiation, while it decreases to 34.76 K/cm, which in the suitable temperature ranges (1-50 K/cm) for SiC crystal growth introduced in reference [14]. The above simulation results illustrate that the reasonable temperature difference can be gotten if radiation heat transfer is considered [7] and the radiation is helpful to uniform the temperature field in the chamber. Figure 8 shows profiles of temperature of powder surface and crystal surface along with r axis with and without radiation heat transfer, respectively. It can be found that: (1) the average temperature at −0.9 < r/R < 0.9 is 3247.9 K without radiation and it decreases to 2651.7 K with radiation on powder surface; (2) the average temperature is 2212.6 K without radiation and it increases to 2555.3 K with radiation on powder surface; (3) the temperature grandient is 373.1 K/cm without radiation, while it decreases to 34.76 K/cm, which in the suitable temperature ranges (1-50 K/cm) for SiC crystal growth introduced in reference [14]. The above simulation results illustrate that the reasonable temperature difference can be gotten if radiation heat transfer is considered [7] and the radiation is helpful to uniform the temperature field in the chamber.  Figure 9 shows temperature fields of both axisymmetric coil and spiral coil in the furnace at x = 0 cross section, and some similar results with two dimensional axisymmetric model [7,14,18] are obtained: (1) SiC powder and graphite crucible regions own the highest temperature (about 2700~2980 K) which is helpful for SiC to sublimate steadily, this is because the induction coil and induction heat locate in the areas surrounding graphite crucible. The result is coincided with that in the reference [14,18], and the highest growth temperature is about 3000~3100 K on the similar condition; (2) furnace top region owns lower temperature because of continuous heat release from top cooling hole, and then good temperature gradient is formed. Moreover, it can be found obviously that the temperature field of the spiral coil is non-axisymmetric, especially in the powder high temperature region. So it is necessary to simulate the temperature field in the SiC crystal growth furnace by three-dimensional models.  Figure 9 shows temperature fields of both axisymmetric coil and spiral coil in the furnace at x = 0 cross section, and some similar results with two dimensional axisymmetric model [7,14,18] are obtained: (1) SiC powder and graphite crucible regions own the highest temperature (about 2700~2980 K) which is helpful for SiC to sublimate steadily, this is because the induction coil and induction heat locate in the areas surrounding graphite crucible. The result is coincided with that in the reference [14,18], and the highest growth temperature is about 3000~3100 K on the similar condition;

Results and Discussion
(2) furnace top region owns lower temperature because of continuous heat release from top cooling hole, and then good temperature gradient is formed. Moreover, it can be found obviously that the temperature field of the spiral coil is non-axisymmetric, especially in the powder high temperature region. So it is necessary to simulate the temperature field in the SiC crystal growth furnace by three-dimensional models. Defects in SiC crystal-such as dislocation, microtubules, etc.-are generated mainly by nonuniform temperature gradient [2,13,19]. Figure 10 presents fields of temperature and temperature gradient (T_R, K/m) in r-direction on rθ cross section in the crystal. A good axisymmetric temperature field is shown in Figure 10a compared with non-axisymmetric temperature field of the spiral coil as shown in Figure 10b. It can be seen from Figure 10b that the highest temperature region deviates from the crystal center, which is as same as the results of Böttcher et al. [13]. The axisymmetric T_R is also gotten for the axisymmetric coil. It is beneficial to grow SiC crystal with fewer defects. As seen from Figure 10d, the T_R is non-axisymmetric for the spiral coil, and the worse T_R are presented in the region far from the center of crystal which may generate more defects. Defects in SiC crystal-such as dislocation, microtubules, etc.-are generated mainly by non-uniform temperature gradient [2,13,19]. Figure 10 presents fields of temperature and temperature gradient (T_R, K/m) in r-direction on rθ cross section in the crystal. A good axisymmetric temperature field is shown in Figure 10a compared with non-axisymmetric temperature field of the spiral coil as shown in Figure 10b. It can be seen from Figure 10b that the highest temperature region deviates from the crystal center, which is as same as the results of Böttcher et al. [13]. The axisymmetric T_R is also gotten for the axisymmetric coil. It is beneficial to grow SiC crystal with fewer defects. As seen from Figure 10d, the T_R is non-axisymmetric for the spiral coil, and the worse T_R are presented in the region far from the center of crystal which may generate more defects. Defects in SiC crystal-such as dislocation, microtubules, etc.-are generated mainly by nonuniform temperature gradient [2,13,19]. Figure 10 presents fields of temperature and temperature gradient (T_R, K/m) in r-direction on rθ cross section in the crystal. A good axisymmetric temperature field is shown in Figure 10a compared with non-axisymmetric temperature field of the spiral coil as shown in Figure 10b. It can be seen from Figure 10b that the highest temperature region deviates from the crystal center, which is as same as the results of Böttcher et al. [13]. The axisymmetric T_R is also gotten for the axisymmetric coil. It is beneficial to grow SiC crystal with fewer defects. As seen from Figure 10d, the T_R is non-axisymmetric for the spiral coil, and the worse T_R are presented in the region far from the center of crystal which may generate more defects. To study the effect of coil type on the temperature and its gradient evenness (ST, ST_R) quantificationally, the date at r/R = 0.38 are extracted and processed by the sample standard deviation method To study the effect of coil type on the temperature and its gradient evenness (S T , S T_R ) quantificationally, the date at r/R = 0.38 are extracted and processed by the sample standard deviation method where T i is the temperature at mesh cell center, T is the average temperature, and S means bad temperature or its gradient evenness. For the axisymmetric coil, both S T and S T_R equal 0; for the spiral coil, S T and S T_R are 1.59 and 170.93, respectively. So it is suggested that the coil is closed to axisymmetric structure.

Conclusions
A tool to simulate electromagnetic and temperature fields in the SiC crystal growth furnace is developed by using finite volume method with C programing language. LU decomposition algorithm based on GPU technology is developed to accelerate solving process of radiation heat flux, and the least distance method (LDM) is proposed for S2S radiation model to check whether mesh face can see others or not. The electromagnetic and temperature fields are simulated for axisymmetric and spiral coils at I = 1250 A and F = 16 kHz. The radiation and conduction heat transfer are compared with each other quantificationally, and the uniforms of temperature and its gradient are also discussed. The following conclusions are obtained: (1) The coil structure has a great effect on temperature field of SiC crystal growth furnace, spiral coil will cause non-axisymmetric temperature field; (2) Radiation heat flux is about 2~3 orders of magnitude bigger than conduction heat flux in the growth chamber, and radiation is helpful to make more uniform temperature field and decrease temperature gradient between powder surface and crystal surface; (3) The non-axisymmetric temperature and its gradient in the crystal caused by the spiral coil are adverse to grow crystal with fewer defects.
Author Contributions: Chunzhen Yang and Guangxia Liu wrote the paper; Chengmin Chen and Yanjin Hou organized figures; Min Xu and Yongxian Zhang contributed to the paper structure and sentences revise.