Numerical Simulation of InGaSb Crystals Growth under Microgravity Onboard the International Space Station

InxGa1-xSb bulk crystal was grown using a GaSb(seed)/InSb/GaSb(feed) sandwich-structured system onboard the International Space Station (ISS). In order to investigate the transport phenomena especially in terms of interface shapes and dissolution heights, the dissolution process was simulated under a micro-gravity level of the ISS. Simulation results showed that the seed/melt interface was concave towards the seed due to the temperature distribution of the system. This prediction is in good agreement with the results of our previous experimental study.


Introduction
In x Ga 1-x Sb is a promising compound semiconductor which can be utilized for various applications such as optoelectronic and thermo-photovoltaic devices.These applications require high quality bulk homogeneous crystals of In x Ga 1-x Sb (uniform crystal compositions) because of its tunable lattice parameter and wavelength ranging from 6.09~6.48Å and 1.7~6.8m by adjusting its composition [1].It is however difficult to grow such crystals on Earth (under the effect of its gravity) because of the following two thermosolutal effects: (1) First is the thermosolutal convection in the melt, driven by density differences, that gives rise to fluctuations in crystal composition in the grown crystal, and consequently lowers crystal quality [2,3]; and (2) Second is the large separation between the liquidus and solidus curves of the system's phase diagram that causes large compositional differences (segregation) between the grown crystal (solid) and the solution (liquid).However, under microgravity environment such as the ISS, such adverse thermosolutal effects can be minimized.This is mainly due to the suppression of convective flow in the solution (melt).Therefore, several experimental and numerical studies were performed for a better understanding of the InGaSb crystal dissolution and growth process by comparing the terrestrial and microgravity conditions [3][4][5][6][7][8][9].Also recently Inatomi et al. have compared the InGaSb crystals grown on the ISS and on Earth, and demonstrated that the shape of the growth interface was a slightly concave towards the seed under microgravity [10].Numerical simulation studies reported thus far on InGaSb have considered smaller model domains than the actual experimental setups.This was mainly due to the high computational demand required by the boundary fitted coordinate (BFC) system used with a Lagrangian-Eulerian (ALE) method.This technique was also not suitable for handling the deformed interfaces [11].Since the volume-averaging model [12] was appropriate for the InGaSb growth ampoule geometry, the computational cost has been reduced by developing a new volume-averaging model [13].These show that there is a need for the development of a new and complete numerical simulation study to model the actual experimental set up to address the issues mentioned earlier.To this end, we have simulated the processes of dissolution of GaSb into InSb in the actual experimental size domain under microgravity condition.A schematic view of the growth/dissolution system (model geometry) is shown in Fig. 1.The rectangular GaSb/InSb/GaSb material sandwich system [5] was stacked in a quartz ampoule and then sealed with Boron Nitride (BN) and carbon sheet as seen in the figure.The stack of materials in the furnace was subjected to three constant vertical temperature gradients as shown in Fig. 1 (b) with top being hotter and bottom being cooler (the horizontal isothermal lines).The process is as follows.We first heat the sample at 1 °C/sec for 162.38 sec and then keep it at this temperature.

Numerical Methods
The melting points of InSb and GaSb are 525 °C and 712 °C, respectively.Thus, during the process InSb melts first.Then GaSb feed material (solid) dissolves into the molten InSb (the dissolution process) making a mixture of In-Ga-Sb (growth solution).Then, this mixture becomes supersaturated under the applied temperature profile and starts accumulating on the GaSb seed crystal at the bottom, and consequently the single crystal of InGaSb alloy begins to grow on the seed (the growth process).
The physical properties of In-Ga-Sb solution and GaSb solid given in Refs.14 and 15 are used.
Although the solution composition changes with time during the formation of the In-Ga-Sb solution and also during the dissolution of GaSb into InSb melt [6], we have used constant material properties of the In-Ga-Sb solution in the simulations for simplicity.The physical properties of BN and quartz are the same as given in Ref. 16.

Governing Equations
Governing equations of the liquid phase (In-Ga-Sb solution) are the those differential equations obtained from the well-known balance laws of thermofluids, namely conservation of mass, balance of momentum, balance of energy, and conservation of species of mass (mass transport).Only energy balance is written for the solid phases (boron nitride and quartz).In the derivation of these equations we have made the following assumptions:  Liquid phase is an incompressible and Newtonian fluid;  Densities of the liquid and solid phases are constant (except the terms in the momentum equations due to the Boussinesq approximation), and thus volume changes (shrinkage and expansion) due to phase changes are not considered;  Changes in physical properties (as a result of compositional changes) during the formation of the In-Ga-Sb solution and also during the dissolution of GaSb into InSb are negligible.
Under these assumptions, the governing equations of the liquid phase, namely continuity, momentum conservation, energy conservation, and mass transport equations become: where, v is the velocity, t time,  density, p pressure,  kinematic viscosity, g the gravitational acceleration,  T the thermal expansion coefficient, T the temperature,  C the solutal expansion coefficient, C concentration,  thermal diffusivity and D diffusion coefficient.The governing equation in the solid phases (BN and Quartz) is the energy balance given as: We solved the governing equations using the new volume-average continuum model [13] by applying OpenFOAM.

Numerical model and equations
Figure 1 (c) shows the computational domain of the growth ampoule used for the numerical simulation.In the crystal (solid) section, the volume averaging method was employed.This method utilizes volume fractions of the solid and liquid phases.The volume fractions were determined by using the applied temperature profile, the concentration profile, and the phase diagram.Details on this method can be found in our previous article [13].
The mathematical description of this method is summarized below: where h is enthalpy,  thermal conductivity, c heat capacity.Momentum source terms B and M in Eq. ( 7) represent the body forces due to gravity (in the Boussinesq approximation "B") and interfacial drag (like damping in porous mushy structure in terms of Kozeny-Carman equation "M") [12].In Eqs. ( 8) and ( 9), the energy and species source terms D and F are the diffusion-like and convection-like mixture sources [17] defined as: In this technique, it is important to identify the liquid and solid phases in the system and also how to compute temperature distribution from the enthalpy profile.This has been done by defining volume fractions that enable us to identify the liquid and solid phases in the domain.The volume (2) 011107-3 JJAP Conf.Proc.(2016) 011107 fractions of the solid and liquid phases ( s ,  l ) are related to the solid and liquid mass fractions (f s , f l ) via The solid and liquid mass fractions (f s , f l ) are determined from the binary phase diagram of the InSb-GaSb system.Details of this method of calculating fractions can be found in [13].Temperature distribution is obtained from the enthalpy profile by using latent heat L:

Boundary conditions and discretization
The boundary conditions for the top and bottom walls are adiabatic for the temperature field, no-slip conditions for the velocity field, and no flux condition (in the normal direction to the ampoule walls) for the concentration field.In addition, on the boundaries between the feed and seed crystals and the walls, and the seals (BN and Quartz) are assumed no-slip condition for the velocity field.Also, the same along the normal direction to the wall and seal boundaries; no mass and heat fluxes are allowed.On the outside of the ampoule along the outer walls, the temperature profile has been determined from the measured values from our experiment [10].The governing equations together with the boundary conditions were discretized by the finite volume method, and the pressure-velocity coupling was handled by the Pressure Implicit with Splitting of Operators (PISO) algorithm [18], and were then solved using an open-source CFD code: the OpenFOAM.For the numerical grid we used 20 and 40 uniform segments for the Quartz and BN sections, and 180 uniform segments for the internal domain in the horizontal direction (x-direction) while 120 uniform segments for the GaSb (feed and seed) section, and 180 uniform segments for the InSb section in the vertical direction (y-direction).The simulation was performed under the microgravity level of 10 -4 G that corresponds to a typical average gravity level observed on the ISS.The g-jitter (gravity fluctuation) on the ISS was not considered.The gravity direction was aligned with the axis of the sandwich system and was directed towards the seed crystal.As for the initial conditions, the complete molten state of InSb was regarded as the initial state in the simulation and therefore the initial GaSb concentration in the solution was taken zero.Figure 2 shows the computed temperature field, the liquid fraction and concentration conditions at 80th sec and 150th sec during the heating process.than the periphery (at the same height), and there is small dissolution at this time.However, at the 150th sec, there is a much cooler section in the center of the feed/melt and seed/melt interfaces as seen from the temperature distribution and isothermal lines in Fig. 2 (b).The feed/melt interface is almost flat while the seed/melt interface is slightly concave towards the melt according to the temperature distribution.For the concentration distribution, we can clearly see that dissolution started at the feed/melt and seed/melt interfaces.Figure 3 presents the computed temperature distribution, and the liquid fraction and concentration conditions at the 300th and 8000th sec of the process.After 162.38 sec, the temperature was kept constant as mentioned earlier.Since the highest and lowest temperature levels were kept stable, heat would only come from the vertical walls, and there would be no heat transfer through the top and bottom boundaries (adiabatic).As seen from Fig. 3 (a), as time progresses, the feed and seed sections are warming up.The computed temperature distribution is quite different from that at the 150th sec.The feed/melt interface becomes concave towards the melt and the seed/melt interface becomes concave towards the seed as indicated by the computed isothermal lines.As time progresses, the temperature field changes, and the center of the seed/melt interface becomes hotter than the periphery.This enhances dissolution in the center, making the feed/melt interface concave towards the melt.Concentration distribution also indicates that the dissolution still continues at this period of the process.The same goes for the seed/melt interface; getting hotter in the center leading to more dissolution in the center, and making the seed/melt interface concave towards the seed.The predictions of concave feed/melt and seed/melt interfaces as described above are in agreement with the experiments of Inatomi et al. [10].As seen from the computed values, the feed/melt interface was concave towards the melt all the time, while for the seed/melt interface was initially concave towards the melt and then became concave towards the seed.This observation is due to the change occurring in the temperature field during the process.The process is easier to understand based on the schematics in Fig. 4. Heat flows from the periphery towards the center during the heating period, thus, at a certain height, temperature is higher at the edges than the center.This was the reason that the periphery region dissolved faster, making the feed/melt interface concave towards the melt.During the constant temperature period, the process is as follows.Heat mainly flows from the top corners and transfers to the bottom corners, as seen from the isothermal lines in Fig. 4 (b).At this time, the center of the feed/melt interface was still cooler, while the seed/melt interface was hotter.So the central region of the seed/melt interface dissolved faster than the periphery, making the seed/melt interface concave towards the seed gradually.

Conclusion
The dissolution process of an InGaSb system was simulated under a micro-gravity level of the ISS.Simulation results showed that the seed/melt interface is concave towards the seed while the feed/melt interface is concave towards the melt.These results were due to the temperature distribution of the system.The numerical predictions are in good agreement with the results of our previous experimental and numerical studies.

Fig. 1
Fig. 1 Schematics of (a) the growth ampoule, and (b) the effective ambient temperature profile, and (c) grid system used for numerical simulation.

Fig. 2
Fig. 2 (i) Temperature distribution, (ii) liquid fraction and (iii) concentration for (a) 80 sec and (b) 150 sec during the heating period in the central section of the domain.
Figure2shows the computed temperature field, the liquid fraction and concentration conditions at 80th sec and 150th sec during the heating process.Figure2(a) indicates that the center is cooler

Fig. 3
Fig. 3 Computed temperature distribution (i), liquid fraction, (ii) and concentration (iii) at (a) 300 sec, (b) 8000 sec during the constant temperature period in the center.

Fig. 4
Fig. 4 Mechanism of the interface evolution: (a) heating period and (b) constant temperature period.