Beware of symmetry breaking and periodic flow regimes in axisymmetric CVD reactor setups
Introduction
Buoyancy driven flows – or free/natural convection flows - are flows associated with density variations, which originate from temperature gradients in non-isothermal flow domains. Because of its nonlinear character (see Eq. (2) – right hand side, third term), buoyancy is one of the classical hydrodynamic topics in stability theory, dynamical systems, and nonlinear science (Hoyas et al., 2005, Mancho and Herrero, 2000, Zhou et al., 2004). Apart from buoyancy's fundamental/theoretical importance, buoyancy driven flows are met in numerous industrial applications, which involve temperature gradients (Zhou et al., 2004); such an industrial process is Chemical Vapor Deposition (CVD).
CVD is widely used for the production of thin solid films from gaseous reactants. A stream of (reactive) gases is injected inside the CVD reactor which is equipped with heated surfaces, the substrates or wafers upon which the deposition of the solid film is performed. The stream of gases carries the precursor i.e. the compound which ideally, when reaching the substrate will decompose through a series of thermally activated surface reactions, which result in the solid film deposition. Due to the temperature differences between the inlet of the reactor where the stream of gases is injected and the temperature of the substrate, free convection flows are commonly met in CVD processes. However, free convection is an undesirable operating flow regime (Cheimarios et al., 2012) because key properties of the deposited film, such as thickness uniformity, do not meet the criteria for industrial use. For this reason, the reactor is preferable to operate in the forced convection flow regime.
The nonlinear nature of the physical/chemical phenomena inside the CVD reactor can generate a multiplicity of flow fields i.e. different flow regimes for exactly the same operating conditions. From a physical point of view, this so-called solution multiplicity is attributed to the different physical mechanisms which are mirrored on the set of partial differential equations describing the transport phenomena. For that, selecting and controlling the flow regime that the reactor will operate is far from trivial and becomes much more complex when designing a new CVD process and/or reactor (Huang and Weng, 2015, Skibinski et al., 2016; Varshney et al., 2009).
From an experimental point of view, smoke tracer experiments (Gadgil, 1993, Wang et al., 1986) or laser holography techniques (Jensen, 1989) although they provide valuable information about the flow field, they are limited for certain cases (Jensen, 1989). Furthermore experimental analysis can cause flow disturbances since it introduces measuring probes (Skibinski et al., 2016). On the other hand, the early methods of scaling analysis based on the dimensionless Reynolds, Grashof and Rayleigh numbers are not adequate since they can only provide a qualitative description (Vanka et al., 2004). Numerical simulations can provide in-depth understanding of the flow (Jensen, 1989, Vanka et al., 2004) inside the reactor but they fail in terms of the “big picture”, i.e. stand-alone numerical simulation limit the analysis only to stable solutions. Advantageous operating conditions can be hidden in regions where multiple solutions co-exist. Thus, in order to have a full understanding and control of the process, only numerical simulations along with a proper nonlinear analysis tool can provide an in-depth understanding of the entire process. Characteristic examples of solution multiplicity in CVD processes are mixed convection flows (Cheimarios et al., 2011, Van Santen et al., 2001) where free and forced convection mechanisms interact, symmetry breaking (Fotiadis et al., 1990, Jensen, 1989, Koronaki et al., 2014, Van Santen et al., 2000) and periodic solutions (Koronaki et al., 2016).
In this work, we compute the solution space – including both stable and unstable solutions - of a CVD reactor via 3D computations. In a region of realistic operating conditions of a CVD reactor we trace all the aforementioned attainable co-existing solutions and we prove - through comparison with computations neglecting the buoyancy term (Gr = 0) - that the solution multiplicity is due to existence of buoyancy forces. Furthermore, the 3D computations allow us to capture nonlinear phenomena such as “symmetry breaking”; the symmetric flow pattern imposed by the symmetric geometric characteristics of the reactor “breaks” and non-axisymmetric solutions - which greatly affect the properties of the deposited film (Chiu et al., 2000, Koronaki et al., 2014) – appear in a region where the flow is laminar. Such non-axisymmetric laminar flows are computed in regions where mixed convection flows are observed and cease to exist when free or forced convection dominates. The latter is in contrast to the findings presented in a previous work of Van Santen et al. (2000), in which they report that symmetry breaking is due to buoyancy (free convection) effects only.
The use of 2D vs 3D computations is a long-lasting debate from the early years of the computational analysis of CVD processes (de Keijser et al., 1988, Fotiadis et al., 1990). The computed solution space in this work is far more complex compared to the one constructed via 2D axisymmetric computations in our previous work (Cheimarios et al., 2011) and the work of others (Van Santen et al., 2001) for the same reactor and operating conditions. This raises major doubts for the use of 2D axisymmetric models especially in regions of mixed convection flows (Vanka et al., 2004) where symmetry breaking phenomena are observed.
From a computational point of view, computing the solution space is far from trivial, especially when the analysis concerns tracking of unstable solutions while working with commercial CFD software. The capabilities and tools of commercial CFD software (user friendly interfaces, advanced mesh generators, state of the art model integration, high performance solvers etc.) make them a valuable tool in process analysis, design and optimization. Still, standalone commercial CFD software such as Comsol Multiphysics (hereafter called Comsol) (“Comsol,” 2018), Fluent (“Ansys/Fluent,” 2018) and CFD-ACE++ (“CFD-ACE++,” 2018) are unable to converge to unstable solutions and/or perform nonlinear analysis tasks in an efficient and robust way. However, in order to compute a full and reliable solution space of the process i.e. a solution space that will not shadow parameter values which may be advantageous to work with, these unstable solutions (if any present) ought to be traced. The computational framework developed in this work expands the capabilities of the commercial CFD software Comsol, through an outer shell module based on the arc-length continuation technique. The latter creates a generalized –here applied in a CVD process- robust nonlinear analysis environment by exploiting all the advanced tools of Comsol, while enabling it to track the solution space accurately, efficiently and fast.
The rest of the paper is structured as follows: We first present the governing equations and boundary conditions followed by the computational framework employed in this work, which enables the computation of all attainable steady state solutions. In the results section we present the computed solution space and discuss the multiple solutions characteristics for a range of reactor operating conditions. The paper ends up with the conclusions.
Section snippets
Problem formulation
The geometry of the problem is that of a single-wafer, cold-wall CVD, stagnation point, axisymmetric reactor. A schematic representation of the CVD reactor along with its dimensions is shown in Fig. 1.
Details for the geometrical characteristics of the reactor can be found in Koronaki et al. (2014) and Van Santen et al. (2000). To obtain the velocity and pressure fields inside the reactor, the continuity (mass) along with the momentum equations are solved. Due to the different thermal conditions
Computer-aided nonlinear analysis
In order to compute the solution space, we implement the so-called pseudo-arc-length continuation algorithm (Kavousanakis et al., 2009, Keller, 1977) in conjunction with the Finite Element Method (FEM) as implemented in the commercial CFD (Multiphysics) software Comsol. The nonlinear system created by the discretization of the set of equations via FEM, in a general case of a N-dimensional space, can be written as,where F is the residual vector, h is the discretized form of the solution
Results and discussion
By implementing the aforementioned computational framework, we compute the solution space, which is presented in Fig. 2; in particular, we show the dependence of the average Nusselt number on the Reynolds number. The Nusselt number, Nu, is computed from,where Awafer, is the surface of the wafer. In all cases the studied flow is laminar with the maximum value of Re, being less than 5. The computations presented below are performed for Twafer = 700 K
Conclusions
We develop a generalized computational framework based on the arc-length continuation method which expands the capabilities of a commercial CFD code and enables it to perform non-linear analysis tasks in an efficient, reliable and robust way. This computational framework is used to construct the solution space of a CVD process in a laminar flow regime. The 3D computations reveal the existence of an interval of mass inflow rates values, in which buoyancy induced, forced convection, symmetry
References (31)
- et al.
Illuminating nonlinear dependence of film deposition rate in a CVD reactor on operating conditions
Chem. Eng. J.
(2012) - et al.
Enabling a commercial computational fluid dynamics code to perform certain nonlinear analysis tasks
Comput. Chem. Eng.
(2011) - et al.
Multiscale modeling in chemical vapor deposition processes: coupling reactor scale with feature scale computations
Chem. Eng. Sci.
(2010) - et al.
Transport phenomena in vertical reactors for metalorganic vapor phase epitaxy. I. Effects of heat transfer characteristics, reactor geometry, and operating conditions
J. Cryst. Growth
(1990) Optimization of a stagnation point flow reactor design for metalorganic chemical vapor deposition by flow visualization
J. Cryst. Growth
(1993)Transport phenomena and chemical reaction issues in OMVPE of compound semiconductors
J. Cryst. Growth
(1989)- et al.
Peculiar asymmetric flow pattern in a vertical axisymmetric VPE reactor
J. Cryst. Growth
(1988) - et al.
A novel free boundary algorithm for the solution of cell population balance models
Chem. Eng. Sci.
(2009) - et al.
Efficient tracing and stability analysis of multiple stationary and periodic states with exploitation of commercial CFD software
Chem. Eng. Sci.
(2016) - et al.
Simultaneous solution of large-scale linear systems and eigenvalue problems with a parallel GMRES method
J. Comput. Appl. Math.
(2009)
Bifurcation detection with the (un)preconditioned GMRES(m)
Comput. Methods Appl. Mech. Eng.
On multiple stability of mixed-convection flows in a chemical vapor deposition reactor
Int. J. Heat Mass Transfer
Symmetry breaking in a stagnation-flow CVD reactor
J. Cryst. Growth
Flow visualization studies for optimization of OMVPE reactor design
J. Cryst. Growth
Effect of the carrier gas on the metal-organic chemical vapor deposition of TiN from tetrakis-dimethyl-amido-titanium
Thin Solid Films
Cited by (0)
- 1
Current address: Scienomics SARL, Paris 75008, France.