Beware of symmetry breaking and periodic flow regimes in axisymmetric CVD reactor setups

https://doi.org/10.1016/j.compchemeng.2019.02.005Get rights and content

Highlights

  • 3D, nonlinear analysis of the solution space of an axisymmetric CVD reactor.

  • Buoyancy gives way to multiple solutions and symmetry-breaking.

  • Forced, free, non-axisymmetric and periodic solutions coexist.

  • Non-axisymmetric solutions appear only in the mixed convection regime.

  • Development of generic computational framework to expand CFD software capabilities.

Abstract

We compute the solution space in the laminar flow regime of a vertical, stagnation point, realistic chemical vapor deposition reactor. Compared to the solution space acquired in our previous work [N. Cheimarios, E.D. Koronaki, A.G. Boudouvis, Enabling a commercial computational fluid dynamics code to perform certain nonlinear analysis tasks, Comput. Chem. Eng. 35 (2011) 2632-2645] for a 2D axisymmetric case, 3D computations reveal a much richer solution space which if not investigated thoroughly crucial information for the process maybe lost. In particular, axisymmetric solutions can co-exist with non-axisymmetric and periodic solutions, due to buoyancy. It is found that non-axisymmetric flows can appear only in the region of mixed convection flows, which argues the use of only 2D axisymmetric models for this type of flows.

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,F(h,p)=0,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,Nu=Lkref(TwaferTinlet)·1AwaferAwaferkTzdA,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)

Cited by (0)

1

Current address: Scienomics SARL, Paris 75008, France.

View full text