Numerical analysis of non-uniform Cu(In, Ga)Se₂ growth in a selenization process on large-area substrates for mass production

ABSTRACT

Growth of a Cu(In, Ga)Se${}_2$ (CIGS) layer during a selenization process is numerically studied to understand mechanisms for formation of stains on large-area substrates batched. CIGS layers need to be uniformly deposited onto the substrates to obtain even conversion efficiency. However, it is difficult to control growth of large-area CIGS layers due to turbulent thermal-fluid flow leaving stains on the substrates. In the present research, the selenization process for an industrial-scale substrates of which sizes are order of square-meters is considered with integrated simulations of detailed key physical processes such turbulent convection, convective-radiative-conductive heat transfer, and chemical reactions. Ascending or descending gas generated by heaters is identified by the time-averaged velocity fields. Descending flow in the passages between substrates produces uneven flow rates across the substrates leading to inhomogeneous supply of heat energy and gas species to the surface chemical reactions. The uneven temperature distribution is the major cause for the stain formation on the substrates. Gross shapes of the stains are found to be well matched with the predicted velocity contour of gas flow above the substrate. The stains are expected to be alleviated by rectifying gas flow such that flow rates become uniform across substrates before entering the passages.

KEYWORDS:

• CIGS solar cell
• selenization process
• stain formation
• large-eddy simulation
• chemical reaction model
• efficiency prediction model

1. Introduction

CIGS solar cells have been receiving attention as a promising thin-film solar cell. The solar cells include a Cu(In, Ga)Se${}_2$ (CIGS) layer as a light absorber which has the highest light-absorption coefficient among compound solar cells (Chen et al., 2017). The high absorption coefficient of the layer guarantees efficient absorption of incident photons even with a thickness of a few micrometers due to the direct band-gap of a CIGS compound. The highest conversion efficiency of CIGS solar cells has been reported to reach 23% at a cell size of $1{\mathrm{cm}}^2$ (Nakamura et al., 2019) and 19% at a small-module size (Stölzel et al., 2019). The reported conversion efficiency is close to that of the most efficient but expensive silicon-based solar cells.

The CIGS layer should have an optimal thickness and chemical composition over a substrate surface to achieve such high conversion efficiency. A CIGS solar cell is composed of stacked multiple layers. Mostly, a soda-lime glass works as a substrate because sodium from the substrate diffuses to a light absorber and helps grain growth (W. Li et al., 2019). A molybdenum (Mo) layer sputtered on the substrate works as a back contact. Above the Mo layer, a CIGS layer is deposited as a P-type compound semiconductor, and on top of that, a CdS compound is deposited as a buffer layer. Among these layers, high-quality production of the CIGS layer is crucial for maximizing the generation of charge carriers and their current across the layers. Indeed, quantum efficiency of a CIGS layer, which represents a ratio of the number of collected charge carriers to the number of incident photons on a solar cell, is a function of a CIGS layer thickness and a light-absorption coefficient. Moreover, the light-absorption coefficient is determined by band-gap energy which is adjusted by chemical composition ratios of gallium (Ga), copper (Cu), and indium (In) compounds inside the CIGS layer. These facts inevitably raise the necessity of sophisticated manufacturing processes that can produce a uniform layer of a CIGS compound to keep the optimal thickness and chemical composition. Various processes have been researched and developed for producing a high-quality CIGS layer. Among these processes, a sputtering-selenization process has been considered as a method for mass production of large-area CIGS modules since it can minimize production time and material usage (Frantz et al., 2015; Gulkowski & Krawczak, 2020).

However, producing a large-area CIGS layer with a uniform thickness and homogeneous composition remains a challenge in industrial production. It is related to difficulties involved in optimization of a selenization process. In detail, the sputtering-selenization method consists of a sputtering process and a selenization process. At the sputtering process, metal precursors, Cu, In, and Ga, are uniformly sputtered on a Mo-coated substrate. Difficulties are mainly involved in the second process known as selenization. In the selenization process, the uniformly deposited precursor are annealed with hydrogen selenide/nitrogen gas (H${}_2$Se/N${}_2$) atmosphere up to about 770 Kelvin in a chemical vapor deposition (CVD) reactor. The H${}_2$Se gas highly reacts with the metal precursors and finally, forms Cu(In, Ga)Se${}_2$ grains. Therefore, maintaining a uniform temperature and homogeneous supply of the H${}_2$Se gas to the substrates is essential for achieving a uniform thickness and homogeneous composition in the CIGS layer. However, an industrial-size reactor is several meters in length and has many slots of narrow gaps to batch multiple substrates. Electric heaters are attached to the reactor walls in order to reach annealing temperature and develop turbulent flow of gas-species through natural convection. Process conditions and design parameters such as the reactor dimension, heater location, temperature profile, or any structure that changes the reactor shape may affect temperature distribution over the substrates and mass transport characteristics inside the reactor, such as flow velocity, flow pattern, and gas concentration. Therefore, it is difficult to optimize a selenization process to provide a uniform temperature and homogeneous supply of reactants across the whole substrates installed. Indeed, Yu et al. (2018) have reported that stain with uneven chemical composition appeared on meter-size solar modules produced in an atmospheric pressure CVD (APCVD) reactor for mass production. Consequently, as the area of the stained region is increased, the module efficiency is decreased.

In-situ evaluation of the influence of each design or process parameter on a selenization process presents a variety of challenges that hinder the use of experimental techniques. An APCVD reactor is a closed chamber to seal toxic gases inside and has an internal boat which is made of quartz or an opaque carbon-carbon material to withstand high temperature gradients and to make radiation effective. Unfortunately, flow visualization techniques mostly use optical systems, so there is a requirement that the system under observation must have transparent windows to transmit and receive light. For example, Liu et al. (2017) and Reinhold-López et al. (2012) used particle image velocimetry to visualize flow patterns in a thermal reactor. Two transparent windows were fitted on two ends of the cylindrical reactor. Similarly, Woods et al. (2004) developed a high-pressure CVD system with an integrated real-time optical characterization technique using laser light scattering. The characterization system could monitor gas-flow kinetics and precursor decomposition kinetics. They made slits that pass through the reactor wall to send and receive scattered light. Inevitably, the experimental techniques (Liu et al., 2017; Reinhold-López et al., 2012; Woods et al., 2004) require design changes in the reactor shapes or materials, increasing manufacturing costs and time, for transmitting light. Besides, the structural changes could increase the possibility of leakage of the extremely toxic and flammable H${}_2$Se gas inside. Due to such limitations, the experimental studies have been mainly focused on a laboratory-scale reactor where the flow regime is limited to laminar flow and thus can be easily characterized. Another challenge is measuring temperature distributions on module surfaces. Temperature sensors may suffer from thermal and chemical damage by the high temperature ambient, and their bundle of wires can disturb heat transfer and gas-flow in narrow gaps between installed substrates. For example, Yu et al. (2018) was only able to attach 9 temperature sensors on each substrate to avoid such problems, thus obtaining discontinuous temperature distributions during the selenization process.

In order to enhance a comprehensive understanding of physical phenomena during a selenization process, numerical simulations have been drawing attention as a method of investigating gas-flow, heat transfer, and chemical reactions, and as a method of modeling electrical properties of deposited layers. Previous numerical studies related to a deposition process can be classified into three groups by physical phenomena that are considered as follows: (i) gas-flow and heat transfer (Choi & Kim, 2014; Vahedein & Schrlau, 2016, 2018; Yu et al., 2018), (ii) gas-flow and heat transfer with detailed chemical reactions (Luo et al., 2021; Noh et al., 2019; Vorob'ev et al., 2002), and (iii) conversion efficiency of deposited layers (Decock et al., 2010; Delgado-Sanchez et al., 2017; Ghorbani et al., 2020; Houck et al., 2019; C. Li et al., 2021; Mohammadnejad & Parashkouh, 2017). Previous numerical studies of a selenization process mainly belong to the group (iii), where electrical characteristics of deposited layers are modeled and optimized according to band-gap energy, layer thickness, relative permittivity of CIGS-CdS-ZnO layers, doping rates, defect density, and material parameters through SCAPS-1D or AMPS-1D simulators (Decock et al., 2010; Gloeckler et al., 2003; Houck et al., 2019). In contrast, the number of studies in the group (i) and (ii) is relatively limited and they mainly studied the heat distribution during a small-area deposition process where the fluid motion is restricted to laminar flow. Recently, Yu et al. (2018) investigated turbulent thermal-fluid flow during a selenization process of batched large-area substrates in an industrial CVD reactor. They proposed the predicted non-homogeneity in the turbulent flow as the cause of staining that occurred on the fabricated CIGS modules. However, their study omitted reactions of gas- and surface-chemistry, only providing indirect evidence of the cause of the stain formation. Unfortunately, so far, it has been difficult to find numerical studies involving a selenization process with detailed chemical reactions that belong to the group (ii). The above mentioned previous studies have been mostly limited to a ‘cell-scale’ small substrate and considered only one or two aspects among the underlying detailed physical phenomena.

Differently from the previous studies, in the present research, the selenization process for an industrial-scale substrates of which sizes are order of square-meters for mass production is considered with integrated simulations of detailed key physical processes such turbulent convection of gas phase species, convective-radiative-conductive heat transfer through the substrates, and surface-chemistry including abstraction, recombination, and desorption. The goal of the present research is to gain a comprehensive understanding of underlying mechanism for stain formation on a large-area CIGS layer during the selenization process. In order to gain such an understanding, it is necessary to study detailed chemical reactions in the vicinity of a substrate region. Specifically, the chemical reactions considered are the thermal decomposition of H${}_2$Se and the growth of CuInSe${}_2$ and CuGaSe${}_2$ from Se elements, metal precursors, and their binaries. Conversion efficiency of the grown CIGS layer is evaluated through a diffusion equation of charge carriers. Various numerical methods have been developed to analyze fluid flow and heat transfer in the interaction of fluids and solid structures (Aghaei et al., 2021; Ghalandari et al., 2019; Karar et al., 2021; Salih et al., 2019). In this work, a large-eddy simulation (LES) method is applied for predicting turbulent thermal-fluid flow of gas-species. Heat transfer of conduction and radiation is simulated using a finite element method (FEM). The surface-chemistry is modeled by a site concept (Ern et al., 1996) with summarized reaction pathways of a selenization process as will be shown later. Various physical models are integrated into a single numerical framework to simulate the selenization process of the CIGS layer with full details. In addition, a numerical model based on the layer thickness and compound concentration, with which one can predict the conversion efficiency is integrated in this study. The numerical predictions are compared against previous numerical and experimental results. Through a quantitative analysis between non-homogeneity of turbulent thermal-fluid flow and non-uniform growth of the CIGS compound, major sources of the stain formation are elucidated.

The paper is organized as follows. Computational details such as a large-eddy simulation, a finite element method, gas- and surface-chemistry models, and an efficiency-prediction model are described in Section 2. In Section 3, the methods are validated for the prediction of a GaAs film growth rate and the prediction of the CIGS solar-cell efficiency. Numerical results of the selenization process and the details of the stain formation are explained and discussed in Section 4, followed by concluding remarks in Section 5.

2. Computational methodology

2.1. LES and FEM

An LES method is used for simulating fluid motions. The basic idea of the LES is that the Navier–Stokes equations directly solve an eddy motion larger than the grid cell size, whereas a subgrid-scale (SGS) model simulates the eddy motion smaller than the grid cell size. Governing equations for fluid motions are the incompressible Navier–Stokes equations as follows: (1) $\begin{array}{rl}\frac{\mathrm{\partial }{\overline{u}}_i}{\mathrm{\partial }{x}_i}& =0,\end{array}$(1) (2) $\begin{array}{rl}\frac{\mathrm{\partial }{\overline{u}}_i}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\overline{u}}_i{\overline{u}}_j}{\mathrm{\partial }{x}_j}& =-\frac{1}{\rho }\frac{\mathrm{\partial }\overline{p}}{\mathrm{\partial }{x}_i}+\nu \frac{{\mathrm{\partial }}^2{\overline{u}}_i}{\mathrm{\partial }{x}_j\mathrm{\partial }{x}_j}-\frac{\mathrm{\partial }{\tau }_{ij}}{\mathrm{\partial }{x}_j}\\ & +g_i\beta (\overline{T}-{\overline{T}}_{ref}),\end{array}$(2) where the spatial filter which filters out flow scales smaller than the grid cell size is denoted as $\overline{()}$, ρ is the fluid density, p is the pressure, ν is the kinematic viscosity of the fluid, ${\tau }_{ij}$ is the subgrid-scale tensor, $g_i$ is the gravitational acceleration, β is the thermal expansion coefficient. The fourth term on the right-hand side of Equation (2(2) $\begin{array}{rl}\frac{\mathrm{\partial }{\overline{u}}_i}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{\overline{u}}_i{\overline{u}}_j}{\mathrm{\partial }{x}_j}& =-\frac{1}{\rho }\frac{\mathrm{\partial }\overline{p}}{\mathrm{\partial }{x}_i}+\nu \frac{{\mathrm{\partial }}^2{\overline{u}}_i}{\mathrm{\partial }{x}_j\mathrm{\partial }{x}_j}-\frac{\mathrm{\partial }{\tau }_{ij}}{\mathrm{\partial }{x}_j}\\ & +g_i\beta (\overline{T}-{\overline{T}}_{ref}),\end{array}$(2) ) accounts for buoyancy resulting from temperature difference and is calculated using the Boussinesq approximation.

The governing equations are discretized by a finite-volume scheme with second-order accuracy in time and space. Details of the numerical algorithms and the solution methods are described in You et al. (2008). In the present study, a dynamic global-coefficient SGS model (You & Moin, 2007) which is suitable for flow in a complex configuration is employed to model subgrid-scale turbulence.

Scalar transport equations are solved to predict transport of energy and gas-phase chemical species as follows: (3) $\begin{array}{rl}\frac{\mathrm{\partial }\overline{T}}{\mathrm{\partial }t}+{\overline{u}}_j\frac{\mathrm{\partial }\overline{T}}{\mathrm{\partial }{x}_j}& =(\alpha +{\alpha }_t)\frac{{\mathrm{\partial }}^2\overline{T}}{\mathrm{\partial }{x}_j\mathrm{\partial }{x}_j},\end{array}$(3) (4) $\begin{array}{rl}\frac{\mathrm{\partial }\overline{Y_i}}{\mathrm{\partial }t}+{\overline{u}}_j\frac{\mathrm{\partial }\overline{Y_i}}{\mathrm{\partial }{x}_j}& =\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}((D_i+D_{t,i})\frac{\mathrm{\partial }\overline{Y_i}}{\mathrm{\partial }{x}_j})+\frac{W_i}{\rho }{\overline{\omega }}_i,\end{array}$(4) where the quantity of the gas-species i is denoted as $({)}_i$, T is the temperature, α is the thermal diffusivity, $Y_i$ is the mass fraction, $D_i$ is the mass diffusion coefficient, $W_i$ is the molar weight, and ${\omega }_i$ is the molar production rate. ${\alpha }_t$ and $D_{t,i}$ are the thermal and mass eddy diffusivity coefficients, respectively, computed from the dynamic subgrid-scale model.

Conductive and radiative heat transfer during a selenization process is calculated using an FEM. The solid domain of a reactor is divided into small elements where the heat transfer equation is expressed in an integral form (Dhondt, 2004). The governing equation is expressed as follows: (5) $\rho c\frac{\mathrm{\partial }T}{\mathrm{\partial }t}=-\frac{\mathrm{\partial }{q}_j}{\mathrm{\partial }{x}_j}+q_{rc},$(5) where c is the specific heat in a unit of J/(kg·K), $q_j$ is the heat flux of thermal conduction, and $q_{rc}$ is the heat flux from radiation-convection heat transfer. For standard FEM procedures, a weak form of the energy equation is considered, and then a backward Euler scheme is used for time integration.

Thermal conduction is the propagation of internal energy in a solid. Following the Fourier's law, the heat flux equation for thermal conduction is given as follows: (6) $q_j=-k\frac{\mathrm{\partial }T}{\mathrm{\partial }{x}_j},$(6) where k is the thermal conductivity. Thermal radiation is produced from emission or absorption of electromagnetic waves on a solid surface as it undergoes internal energy-state transition. Since a selenization process is a high-temperature process, the thermal radiation is occurred by the heaters attached to the reactor walls. The heat flux of the radiation is the sum of emitted energy and absorbed energy on each surface and can be expressed as follows: (7) $q_{rad}=\epsilon \sigma T^4-\alpha \sigma T_{\mathrm{\infty }}^4,$(7) where ε is the emissivity, α is the absorptivity of thermal radiation, and σ is the Stefan–Boltzmann constant. Thermal radiation is calculated with a view factor to consider a geometry-factor of the reactor. The view factor is a value related to the proportion of the radiation which leaves the surface i that strikes the surface j. The view factor $F_{i\to j}$ is expressed by a twofold integral method as follows: (8) $F_{i\to j}=\frac{1}{A_i}{\int }_{A_i}{\int }_{A_j}\frac{\cos {\psi }_i\cos {\psi }_j}{\pi {\mathcal{R}}^2}\mathrm{d}{A}_j\mathrm{d}{A}_i,$(8) where i and j are the surface element indices at the distance $\mathcal{R}$. $A_i$ and ${\psi }_i$ are the element area and the angle between the surface normal and a ray between the two surface elements, respectively.

In the present study, self-coded in-house programs (You et al., 2008; Yu et al., 2018) are employed. The LES solver and the FEM solver are coupled by imposing boundary conditions on shared interfaces between a solid domain and a fluid domain. For satisfying the boundary conditions, which account for energy balance of conduction, radiation, and convection, temperature values are calculated explicitly and imposed at the interface nodes using a coupling method (Yu et al., 2018). Specifically, applying conservation of energy at the interface wall results in the following equation (9) $-k\frac{T_{wall}-T_{wall-\mathrm{\Delta }x}}{\mathrm{\Delta }x}-(h_c+h_r)(T_{wall}-T_{\mathrm{\infty }})=0,$(9) where $T_{wall-\mathrm{\Delta }x}$ is the solid internal temperature at a distance $\mathrm{\Delta }x$ from the interface wall in the normal direction, $h_c$ is the convective heat transfer coefficient, and $h_r$ is the radiative heat transfer coefficient. Therefore, temperature at the interface wall is determined as follows: (10) $T_{wall}=\frac{T_{wall-\mathrm{\Delta }x}+\frac{h\mathrm{\Delta }x}{k}\cdot T_{\mathrm{\infty }}}{1+\frac{h\mathrm{\Delta }x}{k}},$(10) where $T_{\mathrm{\infty }}$ and $h=h_r+h_c$ are calculated explicitly from the solution at the previous time step. In the present multi-code coupling technique (Alonso et al., 2006), volume-weighted interpolation of temperature in fluid and structure domains is conducted to update the boundary conditions for computation at the next time step.

2.2. Chemical-reaction model

A chemical-reaction model predicts chemical phenomena during a selenization process. The governing equation of any gas-species is Equation (4(4) $\begin{array}{rl}\frac{\mathrm{\partial }\overline{Y_i}}{\mathrm{\partial }t}+{\overline{u}}_j\frac{\mathrm{\partial }\overline{Y_i}}{\mathrm{\partial }{x}_j}& =\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}((D_i+D_{t,i})\frac{\mathrm{\partial }\overline{Y_i}}{\mathrm{\partial }{x}_j})+\frac{W_i}{\rho }{\overline{\omega }}_i,\end{array}$(4) ) which is a spatially filtered transport equation of the mass fraction $Y_i$. In Equation (4(4) $\begin{array}{rl}\frac{\mathrm{\partial }\overline{Y_i}}{\mathrm{\partial }t}+{\overline{u}}_j\frac{\mathrm{\partial }\overline{Y_i}}{\mathrm{\partial }{x}_j}& =\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}((D_i+D_{t,i})\frac{\mathrm{\partial }\overline{Y_i}}{\mathrm{\partial }{x}_j})+\frac{W_i}{\rho }{\overline{\omega }}_i,\end{array}$(4) ), the gas-phase molar production rate ${\overline{\omega }}_i$ of the ith species is expressed in a unit of mol/(m${}^3\cdot$s). as follows: (11) $\begin{array}{rl}{\overline{\omega }}_i& =\sum _{j=1}^{m^g}({\nu }_{ij}^{\mathrm{\prime }\mathrm{\prime }}-{\nu }_{ij}^{\mathrm{\prime }})\\ & \times (k_j^f(T)\prod _{k=1}^{n^g}[C_k{]}^{{\nu }_{kj}^{\mathrm{\prime }}}-k_j^r(T)\prod _{k=1}^{n^g}[C_k{]}^{{\nu }_{kj}^{\mathrm{\prime }\mathrm{\prime }}}),\end{array}$(11) where $m^g$ is the number of gas-phase reactions, $n^g$ is the number of gaseous species, ${\nu }_{ij}^{\mathrm{\prime }}$ and ${\nu }_{ij}^{\mathrm{\prime }\mathrm{\prime }}$ are the forward and reverse stoichiometric coefficients of the ith species in the jth reaction, respectively. $[C_k]$ is molar concentration of the kth species in a unit of g/mol. $k_j^f$ is the forward rate constant modeled by the Arrhenius equation as follows: (12) $k^f(T)=k_0\exp (-E_a/RT),$(12) where $k_0$ is the pre-exponential (or frequency) factor, $E_a$ is the activation energy in a unit of kJ/mol, and R is the gas constant. The reverse reaction constant, $k^r$, can be calculated by dividing $k^f$ by the equilibrium constant. In the present work, the gas-phase reaction considered is the thermal decomposition of H${}_2$Se. The required $k_0$ and $E_a$ are obtained from the experimental study of Verma et al. (1991). They confirmed that the thermal decomposition of H${}_2$Se has an activation energy $E_a$ of 166.27 kJ/mol and a pre-exponent $k_0$ of $5.11\times {10}^{10}{\mathrm{s}}^{-1}$. They also found that most of H${}_2$Se are decomposed to Se${}_6$ at 723 K. Therefore, the present study assumes that H${}_2$Se is decomposed into Se${}_6$ during a selenization process.

Boundary conditions of Equation (4(4) $\begin{array}{rl}\frac{\mathrm{\partial }\overline{Y_i}}{\mathrm{\partial }t}+{\overline{u}}_j\frac{\mathrm{\partial }\overline{Y_i}}{\mathrm{\partial }{x}_j}& =\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}((D_i+D_{t,i})\frac{\mathrm{\partial }\overline{Y_i}}{\mathrm{\partial }{x}_j})+\frac{W_i}{\rho }{\overline{\omega }}_i,\end{array}$(4) ) at the substrate involve detailed surface-chemistry mechanisms. Specifically, mass balance of each species at the substrate surface yields (13) $\rho Y_i(V_i\cdot n)=W_i{\mathrm{\Omega }}_i,$(13) where n is the vector normal to the substrate surface, $V_i$ is the species diffusion velocity in a unit of m/s, $W_i$ is the molar weight in a unit of g/mol, and ${\mathrm{\Omega }}_i$ is the surface molar production rate with a unit of mol/(m${}^2\cdot$s). To obtain ${\mathrm{\Omega }}_i$, rates of the surface chemical reactions of the ith species must be calculated first.

To consider the surface chemical reactions, a site concept (Atkins & De Paula, 2009) is employed in the present study. A site is a virtual location of an atom or a molecule on the substrate surface. The present work assumes that the initial number of Cu, In, and Ga sites are proportional to the optimal chemical composition of a CIGS solar cell, which has $\mathrm{Cu}/(\mathrm{Ga}+\mathrm{In})=0.9$ and $\mathrm{Ga}/(\mathrm{Ga}+\mathrm{In})=0.3$ (Hanna et al., 2001; Rau & Schock, 1999). During the selenization process, the balance of site fraction satisfies the relationship as follows: (14) $\begin{array}{rl}{\sigma }_C& =1-\sum {\sigma }_{\mathrm{Cu}binaries},\\ {\sigma }_I& =1-\sum {\sigma }_{\mathrm{In}binaries},\\ {\sigma }_G& =1-\sum {\sigma }_{\mathrm{Ga}binaries},\end{array}$(14) where ${\sigma }_C$, ${\sigma }_I$, and ${\sigma }_G$ are the site fractions of Cu, In, and Ga, respectively. Cu and In react with Se at around 520 K and become CuInSe${}_2$ at about 730 K. Ga starts to react with Se at around 700 K. Finally those chemicals become Cu(In, Ga)Se${}_2$ at around 800 K. In the present paper, site fractions of produced binaries reach the maximum during the heat-up process. Therefore, site fractions of Cu, In, Ga can be negligible compared with site fractions of binaries during the heat-up process. Here, it is assumed that the ternary compounds of CuInSe${}_2$ and CuGaSe${}_2$ do not take any sites but incorporate in a bulk state.

Surface reactions are classified into three types: abstraction, recombination, and desorption. Abstraction is the reaction of a gaseous species with a surface species. For example, the reaction of H${}_2$Se gas with the precursors belongs to the abstraction. A rate of abstraction is expressed as follows: (15) ${\gamma }_{abs,i}=\frac{p}{\sqrt{2\pi W_{i}RT}}\exp (-E_a/RT)X_i{\sigma }_j,$(15) where the term $p/\sqrt{2\pi W_{i}RT}$ is referred to as a collision flux, $X_i$ is the mole fraction of the ith gas-phase species, and ${\sigma }_j$ is the free site fraction of the jth surface species. The recombination is the chemical reaction between surface species. A recombination rate between the ith and jth surface species can be expressed as follows: (16) ${\gamma }_{rec,i}=A\exp (-E_a/RT)[{\phi }_i][{\phi }_j],$(16) where [${\phi }_i$] is the concentration of the ith surface species with a dimension of mol/m${}^2$ and defined as (17) $[{\phi }_i]=\frac{S_0}{N}{\sigma }_i,$(17) where $S_0$ is the total surface concentration, and N is the Avogadro constant. The total surface concentration $S_0$ for CuIn${}_{0.7}$Ga${}_{0.3}$Se${}_2$, which has a Ga ratio of 0.3, is $4.83\times {10}^{18}\mathrm{atom}/{\mathrm{m}}^2$. Chemical reactions among Cu, In, Ga, and Se reach saturation so that the reaction rate will decay as the concentration of the ith species decreases. The desorption is a chemical reaction in which volatile substances evaporate, such as In${}_2$Se. A desorption rate is expressed as (18) ${\gamma }_{des,i}=A\exp (-E_a/RT)[{\phi }_i].$(18) Finally, using the above rates of progress for three types of the surface reactions, a surface production rate of any chemical species can be given as follows: (19) ${\mathrm{\Omega }}_i=\sum _{j=1}^{m^s}{\nu }_{ij}^s{\gamma }_j,$(19) where $m^s$ is the number of surface reactions and ${\nu }_{ij}^s$ is the stoichiometric coefficient of the ith species in the jth surface reaction.

Table 1 summarizes the gas-phase and surface reaction mechanisms considered in this study. G indicates a gas-phase reaction, and S indicates a surface chemical reaction. The reaction G1 indicates the decomposition of H${}_2$Se. The number of gas-phase species considered is five, including N${}_2$, which dilutes H${}_2$Se in the reactor, and In${}_2$Se, which evaporates during the surface chemical reactions. The number of surface reactions considered is 25. The reactions from S1 to S18 describe the abstraction of H${}_2$Se or Se${}_6$ to the precursors' sites with formed binary compounds. H${}_2$Se or Se${}_6$ that reaches the substrate surface by the diffusion is adsorbed as a function of the gas-phase species concentration and the site fraction. Every single element occupies a respective site and retains its position after forming the binary materials. In S19 to S24, the binary compounds formed by the abstraction become ternary materials, which are CuInSe${}_2$ and CuGaSe${}_2$. The present method assumes that these ternary compounds do not participate in the subsequent reactions. The reaction of S25 belongs to desorption, which is the evaporation of the surface species In${}_2$Se. The heat released on the substrate surface is canceled out by the energy balance.

Table 1. Gas-phase and surface reaction mechanisms in the selenization process.

Reaction ($k=k_0\exp (-E_a/RT)$)	$k_0$ (s${}^{-1}$)	$E_a$ (kJ/mol)	Ref.
G1. H${}_2$Se $\rightleftharpoons$ H${}_2$+1/6Se${}_6$	5.11E+10	166.27	Verma et al. (1991)
S1. Cu${}^{(C)}$+H${}_2$Se $\to$ CuSe${}^{(C)}$+H${}_2$	–	100.0	Verma (1993)
S2. Cu${}^{(C)}$+1/6Se${}_6$ $\to$ CuSe${}^{(C)}$	–	100.0	Verma (1993)
S3. 2Cu${}^{(C)}$+H${}_2$Se $\to$ Cu${}_2$Se${}^{(C)}$+H${}_2$	–	100.0	Verma (1993)
S4. 2Cu${}^{(C)}$+1/6Se${}_6$ $\to$ Cu${}_2$Se${}^{(C)}$	–	100.0	Verma (1993)
S5. Cu${}^{(C)}$+2H${}_2$Se $\to$ CuSe${}_2^{(C)}$+2H${}_2$	–	100.0	Verma (1993)
S6. Cu${}^{(C)}$+1/3Se${}_6$ $\to$ CuSe${}_2^{(C)}$	–	100.0	Verma (1993)
S7. In${}^{(I)}$+H${}_2$Se $\to$ InSe${}^{(I)}$+H${}_2$	–	100.0	Verma (1993)
S8. In${}^{(I)}$+1/6Se${}_6$ $\to$ InSe${}^{(I)}$	–	100.0	Verma (1993)
S9. 2In${}^{(I)}$+3H${}_2$Se $\to$ In${}_2$Se${}_3^{(I)}$+3H${}_2$	–	100.0	Verma (1993)
S10. 2In${}^{(I)}$+1/2Se${}_6$ $\to$ In${}_2$Se${}_3^{(I)}$	–	100.0	Verma (1993)
S11. Ga${}^{(G)}$+H${}_2$Se $\to$ GaSe${}^{(G)}$+H${}_2$	–	100.0	Verma (1993)
S12. Ga${}^{(G)}$+1/6Se${}_6$ $\to$ GaSe${}^{(G)}$+H${}_2$	–	100.0	Verma (1993)
S13. 2Ga${}^{(G)}$+3H${}_2$Se $\to$ Ga${}_2$Se${}_3^{(G)}$+3H${}_2$	–	100.0	Verma (1993)
S14. 2Ga${}^{(G)}$+1/2Se${}_6$ $\to$ Ga${}_2$Se${}_3^{(G)}$	–	100.0	Verma (1993)
S15. 2In${}^{(I)}$+H${}_2$Se $\to$ In${}_2$Se${}^{(I)}$+H${}_2$	–	100.0	Verma (1993)
S16. 2In${}^{(I)}$+1/6Se${}_6$ $\to$ In${}_2$Se${}^{(I)}$	–	100.0	Verma (1993)
S17. 4In${}^{(I)}$+3H${}_2$Se $\to$ In${}_4$Se${}_3^{(I)}$+3H${}_2$	–	100.0	Verma (1993)
S18. 4In${}^{(I)}$+1/2Se${}_6$ $\to$ In${}_4$Se${}_3^{(I)}$	–	100.0	Verma (1993)
S19. CuSe${}^{(C)}$+InSe${}^{(I)}$ $\to$ CuInSe${}_2^{(b)}$	4.46E+07	67.3	W. K. Kim et al. (2005)
S20. 2CuSe${}^{(C)}$+In${}_2$Se${}_3^{(I)}$ $\to$ 2CuInSe${}_2^{(b)}$+Se	2.60E+16	162.2	S. Kim et al. (2005)
S21. CuSe${}^{(C)}$+Ga${}^{(G)}$+1/6Se${}_6$ $\to$ CuGaSe${}_2^{(b)}$	8.99E+09	109.0	W. K. Kim et al. (2007)
S22. CuSe${}_2^{(C)}$+In${}^{(I)}$ $\to$ CuInSe${}_2^{(b)}$	4.09E+12	123.0	W. K. Kim et al. (2006)
S23. GaSe${}^{(G)}$+1/2Cu${}_2$Se${}^{(C)}$+1/12Se${}_6$ $\to$ CuGaSe${}_2^{(b)}$	8.18E+11	123.7	W. K. Kim et al. (2008)
S24. Cu${}_2$Se${}^{(C)}$+Ga${}_2$Se${}_3^{(G)}$ $\to$ 2CuGaSe${}_2^{(b)}$	5.43E+07	129.0	Purwins et al. (2006)
S25. In${}_2$Se${}^{(I)}$ $\to$ In${}_2$Se+${\sigma }_I$	1.00E+12	167.4	Ern et al. (1996)

From the surface production rates, the growth rate of each ternary compound is given as follows: (20) $\begin{array}{rl}{G}_{{\mathrm{CuInSe}}_2}& =V_{{\mathrm{CuInSe}}_2}{\mathrm{\Omega }}_{{\mathrm{CuInSe}}_2},\\ G_{{\mathrm{CuGaSe}}_2}& =V_{{\mathrm{CuGaSe}}_2}{\mathrm{\Omega }}_{{\mathrm{CuGaSe}}_2},\end{array}$(20) where $V_{{\mathrm{CuInSe}}_2}$ and $V_{{\mathrm{CuGaSe}}_2}$ are the molecular volume of each compound in a unit of cm${}^3$/mol. The value of the molecular volume is 58.49 and $52.31{\mathrm{cm}}^3/\mathrm{mol}$ (Beck et al., 2000), respectively. The layer thickness is then calculated by a time integration of the growth rate during the selenization process. The Ga ratio ζ is then given as (21) $\zeta =d_{{\mathrm{CuGaSe}}_2}/(d_{{\mathrm{CuInSe}}_2}+d_{{\mathrm{CuGaSe}}_2}),$(21) where the layer thickness is denoted as d. The band-gap energy of the CIGS solar cell varies according to the Ga ratio. The CIGS layer's chemical formula can be expressed as CuIn${}_{1-\zeta }$Ga${}_{\zeta }$Se${}_2$. The band-gap energy of the compound is then a function of ζ as follows (Dullweber et al., 2000): (22) $E_g(\zeta )=1.02+0.67\zeta +0.11\zeta (\zeta -1).$(22) The light absorption coefficient ${\alpha }_l$ is also a function of the band-gap energy and the incident wavelength. Under a monochromatic light, the light absorption coefficient of the CIGS layer is given as (23) ${\alpha }_l(\zeta )=C_1(1-\exp (C_2(h\nu -E_g(\zeta )))),$(23) where $h\nu$ is the photon energy and $C_1$ and $C_2$ are the model constants. $C_1$ and $C_2$ are set to a value of $7\times {10}^4$ and a value of $-1.45$, respectively, from the experimental data (Minoura et al., 2013).

2.3. Efficiency-prediction model

A cell efficiency η is defined as a ratio of the power output to the incident energy from the sun and defined as (24) $\eta =\frac{P_m}{A_{sol}\times E},$(24) where $P_m$ is the maximum power output, $A_{sol}$ is the area of a solar cell, and E is the amount of light energy into a cell. Cell performance is usually measured under the standard test conditions, which are a cell temperature of ${25}^{\circ }\mathrm{C}$, irradiation of 1000 W/m${}^2$, and an air mass (AM 1.5) spectrum. Presence of clouds in the sky reduces the incoming solar radiation so that irradiation can be reduced to 400 W/m${}^2$. A fill factor FF is one of important parameters to characterize the performance of a solar cell and is defined as (25) $FF=\frac{P_m}{P_{max}}=\frac{V_m\times I_m}{V_{oc}\times I_{sc}},$(25) where $P_{max}$ is the theoretical maximum power, $V_{oc}$ is the open-circuit voltage, and $I_{sc}$ is the short-circuit current. $V_{oc}$ is the voltage when the electric current becomes zero, and $I_{sc}$ is the current when the voltage becomes zero, specifically. Another performance parameter is the internal quantum efficiency (IQE) which can be calculated as a function of wavelength. The IQE represents a ratio of excited electrons to the number of photons absorbed by a solar cell.

To evaluate the performance parameters of a CIGS layer, a diffusion equation of charge-carriers in the CIGS layer is solved. The thickness of the CIGS layer is $1\sim 3\text{μ}\mathrm{m}$ and thus the diffusion of charge-carriers occurs very quickly within the layer reaching a steady-state. Because the present LES–FEM–Chemistry model simulates only the deposition process of a CIGS layer, the efficiency-prediction model assumes that the buffer and window layers are uniformly deposited onto the CIGS layer. The buffer layer is mostly made of CdS by the chemical bath deposition, and the quality of the layer is known to be good enough not to affect the solar cell performance according to the research conducted by Kaur et al. (1980). The window layer is composed of ZnO and protects the solar cell from outer conditions. Therefore, the quality of the light absorption layer is a key for determinig cell performance.

The steady-state diffusion equation of charge-carriers density is given as follows (Soedergren et al., 1994): (26) $D_e\frac{{\mathrm{d}}^{2}n(x)}{\mathrm{d}{x}^2}-(n(x)-n_0)/{\tau }_c+\mathrm{\Phi }{\alpha }_l{\mathrm{e}}^{-{\alpha }_{l}x}=0,$(26) where $n(x)$ is the excess charge density in a unit of cm${}^{-3}$, $n_0$ is the charge density in the dark, $D_e$ is the diffusion coefficient of the charge-carriers in the layer in a unit of cm/s, and ${\tau }_c$ is the mean charge-carrier lifetime, Φ is the incident light intensity, ${\alpha }_l$ is the light-absorption coefficient in a unit of cm${}^{-1}$. The solution of the diffusion equation is derived as follows: (27) $n(x)=C_1{\mathrm{e}}^{x/L}+C_2{\mathrm{e}}^{-x/L}+\frac{L^2{\alpha }_l\mathrm{\Phi }}{D_e(1-L^2{\alpha }_l^2)}{\mathrm{e}}^{-{\alpha }_{l}x}+n_0,$(27) where $C_1$ and $C_2$ are the constants, and L is the diffusion length.

Under illumination without an external voltage, the charge-carriers are efficiently drawn off at the back contact, i.e. $n(0)=n_0$. At the front contact, the charge-carriers are reflected, so the gradient of $n(x)$ becomes zero. These boundary conditions determine the constants $C_1$ and $C_2$. By calculating the gradient of $n(x)$ at x = 0 with a multiply of $q_e{D}_e$, the current density $J_L$ is derived in a unit of mA/cm${}^2$ as follows: (28) $J_L=q_e\mathrm{\Phi }\frac{\begin{array}{c}[-L{\alpha }_l\cosh (d/L)\\ +\sinh (d/L)+L{\alpha }_l{\mathrm{e}}^{-d{\alpha }_l}]L{\alpha }_l\end{array}}{[1-L^2{\alpha }_l^2]\cosh (d/L)}.$(28) The internal quantum efficiency can then be written as follows (Soedergren et al., 1994): (29) $IQE=J_L/q_e\mathrm{\Phi }.$(29) When an external voltage V is applied to the solar cell, the charge-carrier concentration, at the back contact, increases to $n(0)=n_0\exp (q_{e}V/k_{B}Tm)$, where $k_B$ is the Boltzmann constant, T is the temperature, and m is the ideality factor. At the front contact, the gradient of $n(x)$ is still zero. By solving Equation (26(26) $D_e\frac{{\mathrm{d}}^{2}n(x)}{\mathrm{d}{x}^2}-(n(x)-n_0)/{\tau }_c+\mathrm{\Phi }{\alpha }_l{\mathrm{e}}^{-{\alpha }_{l}x}=0,$(26) ) again with these boundary conditions, the current density is (30) $J=J_L-\frac{q_e{D}_e{n}_0\sinh (d/L)}{L\cosh (d/L)}({\mathrm{e}}^{q_{e}V/k_{B}Tm}-1).$(30) From Equations (24(24) $\eta =\frac{P_m}{A_{sol}\times E},$(24) ), (25(25) $FF=\frac{P_m}{P_{max}}=\frac{V_m\times I_m}{V_{oc}\times I_{sc}},$(25) ), (29(29) $IQE=J_L/q_e\mathrm{\Phi }.$(29) ), and (30(30) $J=J_L-\frac{q_e{D}_e{n}_0\sinh (d/L)}{L\cosh (d/L)}({\mathrm{e}}^{q_{e}V/k_{B}Tm}-1).$(30) ), the performance parameters η, FF, and IQE are obtained.

3. Validation of the computational methodology

In this section, validation cases for the turbulent thermal-fluid flow, the chemical reaction model and the efficiency-prediction model are provided. Firstly, validation for the prediction of convective heat transfer on a vertical plate is performed. In the second case, the chemical reaction model predicts the growth rate of a GaAs thin-film with comparison to the results of previous experimental and numerical studies. Lastly, the efficiency and the fill-factor of an 1-D CIGS cell are predicted by the efficiency-prediction model. The model accuracy is quantitatively identified compared to previous numerical and experimental studies.

3.1. Validation of the thermal-fluid flow model

Figure 1 shows a schematic diagram of the computational domain and boundary conditions. Air enters through the bottom surface and exits through the upper surface. The left side wall is a heater with temperature $T_{wall}$. The temperature difference between the heater and the ambient $T_{\mathrm{\infty }}$ causes natural convection on the vertical plate. The grid used in the present simulation is composed of 4 million hexahedral cells, and the size of the first grid cell from the heater wall is less than 0.5 in the wall unit. Periodic boundary conditions are applied in the z direction, and a slip boundary condition is applied to the opposite side of the heater wall.

Figure 1. Schematic diagram of the computational domain and boundary conditions for simulations of natural convection on a vertical flat plate.

Figure 2(a) shows profiles of the average temperature along the wall normal direction. Temperature is normalized by the temperature difference between $T_{wall}$ and $T_{\mathrm{\infty }}$. If the Rayleigh number is less than ${10}^9$, natural convection is considered to be laminar. Profiles of the predicted normalized temperature are found to agree well with the similarity solutions in the laminar flow regime. Since similarity solutions do not exist in the turbulent flow regime, temperature profiles which are numerically predicted are only presented. The profile of the time-averaged velocity is also compared with the experimental data as shown in Figure 2(b). The streamwise velocity is normalized by the free stream velocity along the dimensionless transverse coordinate defined as $\zeta =-y(\mathrm{\partial }\theta /\mathrm{\partial }y{)}_{y=0}$. The Grashof number in the present study is $G{r}_x=3.94\times {10}^{10}$, and the simulation result is in good agreement with the experimental data at $G{r}_x=4.83\times {10}^{10},6.81\times {10}^{10},6.68\times {10}^{13}$. In those Grashof number regimes, the buoyancy force dominates the viscous force causing natural convection.

Figure 2. (a) Profiles of the time-averaged temperature along the wall normal direction. Present simulation ▵, $Ra=2.89\times {10}^4$; °, $Ra=3.62\times {10}^6$; □, $Ra=2.89\times {10}^7$; – –, $Ra=8.23\times {10}^9$; $\cdot \cdot \cdot$, $Ra=2.11\times {10}^{10}$. Similarity solution –, $Ra=2.89\times {10}^4$; –, $Ra=3.62\times {10}^6$; –, $Ra=2.89\times {10}^7$. (b) Profile of the time-averaged velocity along the wall normal direction. –, present study ($G{r}_x=3.94\times {10}^{10}$); ▴, Cheesewright (1968) ($G{r}_x=4.83\times {10}^{10}$); •, Smith (1972) ($G{r}_x=6.81\times {10}^{10}$); ▪, Miyamoto et al. (1982) ($G{r}_x=6.68\times {10}^{13}$).

3.2. Validation of the chemical reaction model

The chemical reaction model is validated for predicting the growth rate of a GaAs film. Ern et al. (1996) conducted the numerical study on a CVD process to deposit a GaAs film by injecting trimethylgallium (Ga(CH${}_3$)${}_3$) and arsine (AsH${}_3$) gases in a horizontal channel reactor. Ern et al. considered detailed chemical reactions in gases and on substrate surfaces, which were also studied by Mountziaris and Jensen (1991).

A schematic illustration of the computational domain and boundary conditions for the horizontal channel reactor is depicted in Figure 3. A susceptor is attached to the right behind the substrate located on the reactor bottom. An operating gas flows into the domain at a constant velocity of 0.193 m/s with no-slip conditions on the side walls. The wall-normal gradients of temperature and species-concentration are zero on the wall boundaries. A constant temperature condition of 948 K is applied to the substrate and the susceptor. Ga(CH${}_3$)${}_3$ and AsH${}_3$ are diluted in the carrier gas H${}_2$ and fed into the reactor. The partial pressure of Ga(CH${}_3$)${}_3$ is $1.8\times {10}^{-4}\mathrm{atm}$ and the partial pressure of AsH${}_3$ is $3.3\times {10}^{-3}\mathrm{atm}$. The total pressure is 1 atm. Under these conditions, the Reynolds number Re is 190 based on the reactor length in the streamwise direction. The computational domain is filled with a full hexahedral mesh where a total of 430,920 grid cells ($n_x\times n_y\times n_z=38\times 90\times 126$) is used. At least 4 cells are placed within the turbulent wall-unit of unity near the walls. The gas-phase and surface chemical reactions considered are described in detail by Ern et al. (1996) and Mountziaris and Jensen (1991).

Figure 3. Schematic diagram of the computational domain and boundary conditions for the simulation of the horizontal channel reactor. The substrate and the susceptor temperature is maintained at 948 K.

The results of the present study are compared to the experimental (Reep & Ghandhi, 1983) and numerical ones (Ern et al., 1996). Figure 4 shows the temperature isolines at the middle plane ($y=3.6\mathrm{cm}$). The growth rates of the GaAs thin-film for varying Ga(CH${}_3$)${}_3$ partial pressures are also depicted in the figure. Specifically, Figure 4(a) shows the numerical result of Ern et al. (1996) for the temperature distribution in the reactor and (b) shows the result of the present study. High temperature gradients are shown along the substrate and the susceptor with the outlet temperature of 700–984 K. Figure 4(c) shows the growth rates of the GaAs thin-film on the substrate as a function of the partial pressure of Ga(CH${}_3$)${}_3$. The GaAs growth rates are linearly proportional to the Ga(CH${}_3$)${}_3$ pressure and therefore the film growth is significantly affected by the transport of the Ga(CH${}_3$)${}_3$ gas. The present results are consistent with those of Ern et al.. However, at a low Ga(CH${}_3$)${}_3$ pressure, the growth rate predicted by the present work is found to be slightly closer to the experimental result (Reep & Ghandhi, 1983). In general, the growth rates in the numerical studies are predicted to be higher than experimental results because the numerical methods assume that the surface reaction occurs in a single control volume, which is right on the surface. Therefore, the chemical reaction can be somewhat overly predicted.

Figure 4. Contours of isotherms in the reactor with the susceptor temperature of 948 K; partial pressure of Ga(CH${}_3$)${}_3$: $1.8\times {10}^{-4}\mathrm{atm}$; partial pressure of AsH${}_3$: $3.3\times {10}^{-3}\mathrm{atm}$; and total pressure: 1.0 atm. (a) The numerical result of Ern et al. (1996) and (b) the present result. (c) The growth rate of the GaAs layer as a function of Ga(CH${}_3$)${}_3$ partial pressure. – · –, experimental result of Reep and Ghandhi (1983); – –, numerical result of Ern et al. (1996); and —, the present result.

3.3. Validation of the efficiency-prediction model

The National Renewable Energy Laboratory (NREL) (Contreras et al., 1999) reported a CIGS solar cell with an efficiency of 18.8% fabricated by a three-stage co-evaporation process. The deposited CIGS layer had a thickness of $2\text{μ}\mathrm{m}$ and a band-gap energy of 1.2 eV. Song et al. (2004) made an attempt to numerically predict the electrical performance of the NREL's solar cell. They solved a Poisson equation for electrostatic potential and continuity equations for free holes and electrons using the AMPS–1D simulation tool. Similarly, the efficiency-prediction model, in the present work, is validated by estimating the performance of the NREL's solar cell in comparison to Song et al.. The simulation parameters used are the measured thickness ($d=2\text{μ}\mathrm{m}$) and Ga ratio ($\zeta =0.3$) from the NREL (Contreras et al., 1999). The charge-carrier concentration ($n_0=2.5\times {10}^{13}$) and the ideal factor (m = 1.5) are selected as optimal values to fitting the measured performance curves.

The predicted results of the present model are found to well agree with the experimental and the numerical results as shown in Table 2. The values of $V_{oc}$ predicted by the numerical simulations show a slight difference of $3\sim 4\mathrm{mV}$ compared to the experimental result. $J_{sc}$ is predicted as low as 0.8 mA/cm${}^2$ by Song et al. and as low as 0.58 mA/cm${}^2$ by the present model. The fill factor is predicted to be 0.75% higher in Song et al. while 0.23% higher in the present simulation. The conversion efficiency of the CIGS solar cell is predicted to be slightly lower than that in the experiment. The current density–voltage (J–V) curve and the IQE are presented in Figure 5. The present model well predicts the J–V curve as shown in Figure 5(a). Figure 5(b) shows the IQE as a function of wavelengths from 600 nm to 1200 nm with good agreements. Overall, the present model is found to well predict the performance of the CIGS solar cells.

Figure 5. (a) Current density and voltage (J–V) curve and (b) internal quantum efficiency curve as a function of the wavelength. •, Experimental results of Contreras et al. (1999); – –, numerical results of Song et al. (2004); and —, the present results.

Table 2. Comparison of a solar cell performance with the experiment (Contreras et al., 1999) and numerical study (Song et al., 2004) regarding the open-circuit voltage $V_{oc}$, short-circuit current density $J_{sc}$, fill factor FF, and efficiency η.

Parameters	Experiment (Contreras et al., 1999)	Numerical study (Song et al., 2004)	Present model
$V_{oc}$ (mV)	678	675	682
$J_{sc}$ (mA/cm${}^2$)	35.22	34.42	34.64
FF ($\mathrm{\%}$)	78.65	79.4	78.88
η ($\mathrm{\%}$)	18.8	18.39	18.64

4. Simulation of the selenization process

The selenization process is a chemical vapor deposition process. Cu, Ga, and In are sputtered on substrates and exposed to H${}_2$Se gas. The H${}_2$Se gas is directly absorbed by the sputtered substrates or decomposed into Se elements. As the absorbed or decomposed Se reacts with the sputtered atoms, Cu, Ga, and In produce CuSe${}_2$, CuSe, Cu${}_2$Se, GaSe, Ga${}_2$Se${}_3$, In${}_4$Se${}_3$, and InSe binaries. Those binary compounds become ternary compounds which are CuInSe${}_2$ and CuGaSe${}_2$ at around 730 K. Finally, the CIGS compound, Cu(In, Ga)Se${}_2$, is formed at around 800 K. Those surface reactions are very sensitive to the substrate temperature. In following sections, the velocity and temperature fields during the heat-up and stabilization steps are analyzed. Chemical reactions in the selenization process including dissociation of H${}_2$Se gas and growth of the CIGS layer are also considered. Based on the predicted layer thickness and the Ga ratio, the conversion efficiency is evaluated using an efficiency-prediction model.

4.1. Computational configurations

A comprehensive simulation of a selenization process for batched substrates of several square meters is provided in this section. The structure of an industrial CVD reactor considered is depicted in Figure 6. Due to some confidential issues, the exact shape of the reactor and the detailed selenization recipes can not be presented here. Briefly, as shown in Figure 6(a,b), heaters surround side-, top-, and bottom-walls, and two boats are located on the front side and the rear side in the reactor. 30 substrates are vertically batched with narrow intervals in the boats. Physical quantities are normalized by their reference values selected.

Figure 6. Schematic diagram of the thermal reactor. (a) The reactor is surrounded by heaters and (b) 30 substrates are placed in the reactor.

As the heater temperature increases, natural convection develops along vertical walls. The turbulent convection transports energy and chemical species to the substrates. Since the reactor is sealed during the process, all boundaries are treated with no-slip conditions for the momentum equations. On the heater surfaces, temperature profiles from the experiment (Yu et al., 2018) are directly imposed with Dirichlet conditions. Adiabatic conditions are imposed outside of the reactor. For chemical species, the mass fraction is assumed to be zero at the wall and the heater, whereas normal gradients of the mass fraction become zero on the substrate surface.

The H${}_2$Se gas diluted by N${}_2$ is injected into the thermal reactor and the total pressure becomes 0.92 atm. The partial pressure of H${}_2$Se gas is 0.02 atm, and the partial pressure of N${}_2$ is 0.90 atm. It is noted that a mixture of the working gases used in the selenization process does not participate in thermal radiation since the nitrogen gas, which makes up the majority of the total pressure, is transparent to infrared radiation.

The normalized process time τ is defined as a ratio of elapsed time to the entire process time. The normalized heater temperature θ is defined as follows: (31) $\theta =\frac{T-T_{min}}{T_{max}-T_{min}},$(31) where $T_{min}$ is the minimum heater temperature and $T_{max}$ is the maximum heater temperature. The selenization process is divided into two sub-steps. The first sub-step is referred to as a heat-up step, and the second step is a stabilization step. The heaters increase the temperature inside during the heat-up step and then are adjusted to maintain a constant temperature during the stabilization step. The material properties of the working gases are given as a function of temperature.

A solid domain is constructed with a full hexahedral mesh of 0.1 million grid cells. A grid for the fluid domain consists of a total of 7 million hexahedral cells. Since the turbulent convection occurs due to the high temperature-gradients near the heaters, at least 4 cells are allocated within the distance of 1 wall unit from the fluid boundaries. No meaningful differences in the simulations results were found by further refinement of grid resolution for fluid and solid domains to 12.7 million cells and 0.28 million cells, respectively.

4.2. Thermal-fluid flow during the selenization

Accurate prediction of fluid motions and heat transfer during the selenization is important because the concentration and reaction of chemical species are significantly influenced by the thermal-fluid flow inside. The time-averaged velocities and temperature fields at several locations are investigated at two normalized time points, $\tau =0.39$ and $\tau =0.78$. Specifically, $\tau =0.39$ is in the heat-up step with a ramping rate of 0.09 K/s, and $\tau =0.78$ is in the stabilization step where the heater temperature is kept constant. Based on the average temperature and the mean velocity magnitude, the heat-up step is in the turbulence regime of $Ra=1.7\times {10}^9$ and $Re=3.8\times {10}^3$.

Contours of the y-directional velocity and temperature are given in Figure 7 on an xy-plane slicing the middle of the reactor. The velocity and temperature are averaged for 5 minutes which corresponds to $\mathrm{\Delta }\tau =0.025$ at $\tau =0.39$. Figure 7(a) shows the magnitude of y-velocity. To develop circulating flow inside, vertical heaters on the side walls have higher temperature profiles than the bottom heater. In particular, due to the high temperature on the rear-side heater, strong upward flow forms near the rear wall. When the rising flow reaches the top, it starts to descend toward the bottom of the reactor. The descending flow passes through narrow passages between the substrates as shown in Figure 7(a). However, the descending flow is not uniform across the substrates. Furthermore, some reverse flow coming out from the passages is often observed at the top of the substrates. The time-averaged temperature distribution at $\tau =0.39$ is presented in Figure 7(c). As expected, the raised hot gas is extensively distributed at the top area of the reactor. The hot gas heats the top region of the substrates first and gradually raises the temperature of the entire region.

Figure 7. Contours of the y-directional velocity time-averaged for $\mathrm{\Delta }\tau =0.025$ in the mid-section x−y plane of the reactor at (a) $\tau =0.39$ and (b) $\tau =0.78$. Contours of the temperature time-averaged for $\mathrm{\Delta }\tau =0.025$ in the mid-section x−y plane of the reactor at (c) $\tau =0.39$ and (d) $\tau =0.78$.

At the stabilization step, $\tau =0.78$, similar flow patterns and similar temperature distributions appear. Figure 7(b) presents the magnitude of the y-velocity. Strong upward flow still exists near the rear heater whereas strong downdraft newly appears at the front wall. The downdraft blows the stagnant hot gas down as shown in Figure 7(d). In the stabilization step, relatively uniform distributions of the flow velocity and temperature are observed in the passages between the substrates compared to the heat-up step.

The surface temperature of the substrate is an important factor that determines the rates of surface chemical reactions. Figure 8 shows time-averaged temperature contours on 15th and 23rd substrates at $\tau =0.39$. The present results are averaged for 5 minutes ($\mathrm{\Delta }\tau =0.025$) and the experimental data from 9 temperature-sensors are also averaged for the same time interval. The temperature is normalized by the difference between the maximum and the minimum temperatures on the substrate. Considering that the temperature fields were averaged for only 5 minutes which are far from the averaging duration for full statistical convergence for both experimental and numerical results, the comparison is rather qualitative than quantitative. It is found that both measured and predicted temperature fields have radial patterns such that a low temperature zone appears in the center of the substrates. Symmetric temperature fields are observed in the simulation result while the temperature distributions are found to be anti-symmetric in the experiment as shown in Figure 8. Since the amount of heat energy generated by the left and right heaters is the same, it is considered that anti-symmetry of outgoing heat flux from the external surface of the reactor, which is difficult to precisely control in the experiment, caused anti-symmetric distributions of temperature on substrate surfaces in the experiment. Figure 9 shows the spatially averaged temperature for all substrates batched during the selenization process. The black line indicates the experimental result, whereas the red line indicates the predicted result.

Figure 8. Temperature distributions on the 15th and 23rd substrates at $\tau =0.39$. Temperature is averaged for 5 minutes ($\mathrm{\Delta }\tau =0.025$). The experimental results: (a) the 15th and (b) the 23rd substrates. Results of the present study: (c) the 15th and (d) the 23rd substrates.

Figure 9. Spatially averaged temperature profile during the selenization process. The black line indicates temperature measured by the experiment, while the red line indicates the temperature predicted by the present study.

4.3. Chemical reactions during the selenization

When temperature increases, H${}_2$Se is decomposed into H${}_2$ and Se${}_6$. The H${}_2$Se gas and the decomposed Se element participate in surface adsorption and react with the deposited precursor to form binary compounds. Therefore, their concentration is important for prediction of the chemical composition of the grown CIGS layer. The present study considers both forward and reverse reactions of the thermal decomposition. The gas concentration (H${}_2$Se and Se${}_6$), y-velocity, and temperature are analyzed on an yz-plane in the middle of a passage between two substrates. Figure 10 shows time-averaged contours on the yz-plane in the middle of the passage between the 1st and 2nd substrates for $\tau =0.42$. In Figure 10(a), strong ascending flow is developed along the left and right walls and then the flow goes down right after reaching the top wall. The temperature distribution is presented in Figure 10(c), and the mole fraction of H${}_2$Se is shown in Figure 10(b). As H${}_2$Se is rapidly decomposed into H${}_2$ and Se${}_6$ at a high temperature environment, the mole fraction of H${}_2$Se has the opposite distribution against the temperature field, i.e. a high-temperature region has a low mole fraction of H${}_2$Se. In the same context, the mole fraction of Se${}_6$ has the distribution opposite to the mole-fraction field of H${}_2$Se as shown in Figure 10(d). Since the decomposition rate is an exponential function for temperature, the concentrations of H${}_2$Se and Se${}_6$ are inhomogeneous across the substrate's surface resulting in the non-uniform surface-reaction rates.

Figure 10. Contours time-averaged for 5 minutes ($\mathrm{\Delta }\tau =0.025$) near the 1st substrate at $\tau =0.42$: (a) the y-directional velocity, (b) the mole fraction of H${}_2$Se, (c) normalized temperature, and (d) the mole fraction of Se${}_6$.

Figure 11 shows the predicted thickness and Ga ratio in the CIGS layer at the 1st and 16th substrates. The thickness and the Ga ratio are important parameters for evaluating band-gap energy and conversion efficiency of the CIGS solar cell. In the 1st substrate, the layer of CIGS compounds, approximately $2\text{μ}\mathrm{m}$ thick, is deposited as illustrated in Figure 11(a). The layer thickness is about $1.8\text{μ}\mathrm{m}$ on the bottom region but about $2.0\text{μ}\mathrm{m}$ on the top region, specifically. The thickness isolines present a wavy pattern on the top region. On the 16th substrate, the CIGS layer grows much thicker as shown in Figure 11(b). Due to the hot gas blowing from the top and bottom areas of the reactor, the CIGS layer grows much thicker than other substrates, however, the non-uniformity in thickness becomes more severer.

Figure 11. Predicted CIGS layer thicknesses on (a) the 1st substrate and (b) 16th substrate. Predicted Ga ratios on (c) the 1st substrate and (d) the 16th substrate.

Across the whole substrates, the Ga ratio ranges from 0.23 to 0.29, which is close to the optimum value of 0.3. Figure 11(c) shows the predicted Ga ratio on the 1st substrate. The ratio on the top region in the substrate is about 0.26, whereas the bottom region is as high as 0.28. It means that the CIGS layer grows more on the top region with a lower Ga content, whereas grows less at the bottom region with a higher Ga content. As shown in Figure 11(c), the contour presents a wavy pattern on the top and bottom edges where the Ga content is low. This pattern is more pronounced on the 16th substrate as shown in Figure 11(d). The Ga content on the left side of the substrate is about 0.23, showing a marked difference from that on the right side which is about 0.27. It means that CuInSe${}_2$ grows more on the left side, but less on the right side of the substrate and vice versa for CuGaSe${}_2$.

4.4. Conversion efficiency of the solar cell

Using the predicted CIGS thickness and Ga ratio, cell-efficiency can be evaluated from the efficiency-prediction model. The Ga ratio is used to calculate band-gap energy using Equation (22(22) $E_g(\zeta )=1.02+0.67\zeta +0.11\zeta (\zeta -1).$(22) ). The current density is then obtained by Equation (30(30) $J=J_L-\frac{q_e{D}_e{n}_0\sinh (d/L)}{L\cosh (d/L)}({\mathrm{e}}^{q_{e}V/k_{B}Tm}-1).$(30) ) using the predicted thickness and band-gap energy in the consideration of shunt resistance and series resistance measured in the experiment. Figure 12 shows several performance parameters for each substrate. Figure 12(a) shows the conversion efficiency of each substrate under the standard test conditions. The efficiency is predicted to be on average 0.67% higher than the experimental data. The fill factor is measured to evaluate cell performance as a function of the open-circuit voltage and the short-circuit current density. Figure 12(b) presents the open-circuit voltage ($V_{oc}$). The 15th substrate has a slightly higher value in $V_{oc}$ than the experimental results, whereas the 18th, 19th, and 20th substrates present lower values than the experimental data. The overall $V_{oc}$ is predicted to be 1.53 mV lower. The fill factor is given in Figure 12(c). The predicted fill factor has a more uniform distribution than the experimental results. The averaged fill factor is calculated to be 0.23% lower than the experimental values. Meanwhile, the short-circuit current density ($J_{sc}$) is predicted to be lower than the experimental ones on the 15th and 20th substrates as shown in Figure 12(d). The $J_{sc}$ is predicted to be on average 0.22 mA/cm${}^2$ lower. Overall, the performance parameters of each substrate are accurately predicted.

Figure 12. CIGS solar cell performance along the substrates: (a) the conversion efficiency η, (b) open-circuit voltage $V_{oc}$, (c) fill factor FF, and (d) short-circuit current density $J_{sc}$. •, experimental results and red• with a red line, the present results.

In the present work, there are two limitations that explain small deviation in the predicted parameters. First, the prediction model does not account for physics in a molecular scale such as grain boundaries, defects, and doping. The second limitation is that the present work ignores the diffusion of molecules within the CIGS layer by applying the concept of site fraction. By considering the molecular-scale physics and the diffusion of the molecules, it is expected that the prediction accuracy of the model can be further improved (Idris et al., 2015; W. Li et al., 2019).

4.5. Stain formation affecting the cell performance

Factors that affect the stain formation during the selenization are analyzed in this section. Figure 13(a) shows the conversion efficiency on the 16th substrate measured from the experiment. Figure 13(b) presents the electroluminescent image of this substrate. The dark stains in the image, where the efficiency is low, extensively appear on the left side of the substrate. The bright stains, where the efficiency is high, appear on the right side of the substrate. Figure 13(c) and (d) show the predicted temperature and Ga ratio of the substrate, respectively. The temperature on the left side of the substrate is predicted to be significantly higher than that on the right side. Similarly, Figure 13(d) shows the non-uniform distribution of the predicted Ga ratio. The Ga ratio is low on the left side but high on the right side. These distribution patterns are well matched with the efficiency map and the electroluminescent image as shown in Figure 13(a,b) from the experiment. The present work has found that, in the dark stains, the Ga ratio is low due to the high temperature indicating that CuInSe${}_2$ grows excessively compared to CuGaSe${}_2$ resulting in low conversion efficiency. It means that the dark stains will become more pronounced when an uneven distribution of hot gas-flow occurs through the passages between substrates reducing the cell-performance.

Figure 13. Experimental results on the 16nd substrate: (a) contour of the measured conversion efficiency and (b) the electroluminescence image. Numerical results: (c) contour of the temperature time-averaged for 5 minutes ($\mathrm{\Delta }\tau =0.025$) at $\tau =0.42$ and (d) contour of the predicted Ga ratio.

Figure 14. Experimental results on 19th substrate: (a) contour of the measured conversion efficiency and (b) the electroluminescence image. Numerical results: (c) contour of the instantaneous y-directional velocity at $\tau =0.42$ and (d) contour of the predicted Ga ratio.

In the present reactor, wavy stains appear on the top and bottom regions of the substrates. Figure 14(a,b) show the measured efficiency map and electroluminescence image of the 19th substrate, respectively. Notably, a distinct wavy pattern appears in the dark stains at the top region of the substrate. Ascending and descending gas flow near the top region is believed to be the cause of the wavy pattern. Figure 14(c) shows the predicted instantaneous-velocity distribution of the gas flow rising and falling near the top region. Similarly, the Ga ratio is predicted to be low on the top region, exhibiting a wavy pattern, and to be high at the bottom region as shown in Figure 14(d).

In summary, from the present analysis, it is found that the behaviors of the gas flow determine the shape of the stains occurring at the top or bottom of the substrates, while the brightness of the stains, i.e. efficiency, is strongly influenced by the temperature distribution on the substrate. From the binary phase diagram of CuInSe${}_2$ reported in experimental studies (Muzzillo et al., 2016; Regmi et al., 2020), the concentration of Se${}_6$ gas generated by thermal decomposition determines the ratio of binaries that determines the conversion efficiency of solar cells, which supports the present analysis.

5. Concluding remarks

A comprehensive analysis of the selenization on large-area substrates has been performed with the consideration of thermal-fluid flow and detailed chemical reactions in an industrial CVD reactor. The performance parameters, i.e. efficiency, fill-factor, open-circuit voltage, and short-circuit current, are evaluated from the solution of a diffusion equation for charge-carriers in the CIGS thin-layer.

Non-uniform thermal-fluid flow is identified through the vertical velocity fields and the temperature distributions on the substrates. The problem is that the descending flow is not uniformly rectified before entering the passages between substrates. Therefore, the supply of heat energy and gas-species becomes uneven across the substrates. From the instantaneous velocity field, in particular, a strong disturbance is identified near the top regions of the substrates where the descending flow first arrives. In terms of the chemical reactions, the uneven temperature distribution causes spatially non-uniform decomposition of H${}_2$Se gas inside the reactor resulting the uneven growth of the CIGS layer.

In addition, the stain formation on large-area substrates is investigated. The uneven temperature distribution on the substrate is identified as a major cause of the stain formation. Furthermore, the shape of the stain appears to be affected by the velocity pattern of the gas flow. This result suggests that the stain formation due to the uneven distribution of temperature can be alleviated by rectifying the thermal-fluid flow such that the flow rates become uniform before entering the passages between the substrates. Attaching appropriate baffles on the top of the substrates may be a feasible way of achieving this (Fauzi et al., 2018).

As a future direction of the present research, it is considered to incorporate models for molecular-scale dynamics (Idris et al., 2015; W. Li et al., 2019) which can explain the characteristics of micro-structures on the substrate.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the National Research Foundation of Korea (NRF) under the Grant Number NRF-2021R1A2C2092146.