# Modeling and Simulation of  a Chemical  Vapor Deposition

Authors: J. Geiser; M. Arab
Journal: Journal of Applied Mathematics
 (2011)

---

## Abstract

We are motivated to model PE-CVD (plasma enhanced chemical vapor deposition) processes for metallic bipolar plates, and their optimization for depositing a heterogeneous layer on the metallic plate. Moreover a constraint to the deposition process is a very low pressure (nearly a vacuum) and a low temperature (about 400 K). The contribution of this paper is to derive a multiphysics system of multiple physics problems that includes some assumptions to simplify the complicate
process and allows of deriving a computable mathematical model without neglecting the real-life processes.
To model the gaseous transport in the apparatus we employ mobile gas phase streams, immobile and mobile phases in a chamber that is filled with porous medium (plasma layers). Numerical methods are discussed to solve such multi-scale and multi phase models and to obtain qualitative results for the delicate multiphysical processes in the chamber. We discuss a splitting analysis to couple such multiphysical problems. The verification of such a complicated model is done with real-life experiments
for single species. Such numerical simulations help to economize on expensive physical experiments and obtain control mechanisms for the delicate deposition process.


---

## Body

## 1. Introduction

We  motivate our research on producing high-temperature films by low-pressure deposition processes. Such deposition can only be achieved with plasma-enhanced processes; see [1]. In standard applications based on depositing binary systems such as TiN and TiC, new material classes arose. More delicate to deposit are nanolayered ternary metal-carbide or metal-nitride materials, including MAX-phase materials, while we have to control three  species in the deposition process see [2].To understand the growth of a thin film done by PE-CVD (plasma-enhanced chemical vapor deposition) processes,  we have to model a complicated multiphysical process; see [3, 4].Based on the process constraints, which are very low-pressure (nearly vacuum) and a low temperature (about 400 K), we can simplify the complicated model.The delicate modeling part is to model the plasma, in which the deposition species, given as precursor gases, are transported; see [5]. We concentrate on plasma as an ionized medium, known as “cold” plasma or glow discharges; see [6, 7]. To model the gaseous transport of our deposition process in the cold plasma, a porous medium model with different permeable layers is an attractive simulation model; see [4]. In the chamber filled with an ionized plasma that is given as a porous medium, the gaseous transport involves several phases: mobile gas phase streams, and both immobile and adsorbed phases. For such processes, we present a multiphase and multispecies model; see [3, 8]. The species of the gaseous transport are MAX-phase materials; see [9–11], and we concentrate on Ti3SiC2 molecules, so we have Ti, Si, and C in the gas phases.Further a porous medium model with permeable layers is an attractive simulation model also for the electric field in the chamber that is used to control the gaseous transport to the deposition area. Here the expertise of similar simulation models are used; see [12].So the transport, chemical and sorption processes in a heterogeneous medium can be used to simulate species transport in a plasma-enhanced environment, controlled by pressure, temperature and by additional electric fields. Such modeling methods have been developed in recent years and are focused on producing high-temperature films; see [4].Fundamental models of CVD processes employ detailed descriptions of fluid flow, mass transfer, heat transfer and chemical kinetics; see [13, 14].While having a powerful deposition process in low-temperature and low-pressure regimes, we have also taken into account some drawbacks of such a process, which is driven with ionized plasma. Here, we deal with multicomponent problems with heavier and lighter species. Such conditions have consequences in slower (strongly adsorbed) or faster traveling (weakly adsorbed) species.The contributions of this paper are to couple numerical methods that are developed to solve such multiscale and multiphase models and to obtain qualitative results on the delicate multiphysical processes in the chamber. Here splitting schemes are applied as an attractive tool to decouple delicate processes; see [15]. We apply and analyze splitting schemes, in a way that combines discretization methods for partial differential and ordinary differential equations; see [12, 15].The immobile, adsorbed and reaction terms can be treated with fast Runge-Kutta solvers, whereas the mobile terms are convection-diffusion equations and are solved with splitting semi-implicit finite volume methods and characteristic methods, [16].Such a sequential treatment of the partial differential equations and ordinary differential equations allow of saving computational time, while expensive implicit Runge-Kutta methods are reduced to the partial operators and fast explicit Runge-Kutta methods are for the ordinary operators of the multiphase model.With various source terms we control the required concentration at the final deposition area.Different kinetic parameters allow of simulating the different time scales of the heavier and lighter gaseous species.The numerical results are verified with physical experiments and we discuss the applications to the production of so-called metallic bipolar plates.This paper is outlined as follows. In Section2, we present our mathematical model based on the multiphases. In Section 3, we present the underlying physical experiments for the model equations. In Section 4, we discuss discretization and solver methods with respect to their efficiency and accuracy. The numerical experiments are given in Section 5. In Section 6, we briefly summarize our results.


## 2. Mathematical Modeling

In the model we have included the following multiple physical processes, related to the deposition process:(i)
flow field of the ionized plasma: Navier-Stokes equation; see [17],(ii)
transport system of the species: mobile and immobile gaseous phases with Kinetically controlled adsorption (heavier and lighter species); see [4],(iii)
electric field: Poisson equation; see [13, 18].In the following we discuss the three models separately and combine all the models into a multiple physical model. We assume a two-dimensional domain of the apparatus with isotropic flow fields; see [17].


### 2.1. Flow Field

The conservation of momentum is given by (flow field: Navier-Stokes equation)(2.1)∂∂tv+v⋅∇v=-∇p,inΩ×[0,t],v(x,t)=v0(x),onΩ,v(x,t)=v1(x,t),on∂Ω×[0,t],
where v is the velocity field, p the pressure, v0 the initial velocity field and the position vector x=(x1,x2)t∈Ω⊂ℝ2,+. Furthermore, we assume that the flow is divergence free and the pressure is predefined.


### 2.2. Transport Systems (Multiphase Equations)

We model the ionized plasma as an underlying medium in the chamber with mobile and immobile phases. Here transport in the plasma with gaseous species contain of mobile and immobile concentrations, [3]. For such a heterogenous plasma, we applied our expertise in modeling multiphase transport through a porous medium; see [2].To amplify the modeling of the gaseous flow to the chamber which is filled with ionized plasma, we deal with the so-called far-field model based on a porous medium. Here the plasma can be modeled as a continuous flow [17], that has mobile and immobile phases; see [8].We assume nearly a vacuum and a diffusion-dominated process, derived from the Knudsen diffusion, [19]. In such viscous flow regimes, we deal with small Knudsen numbers and a pressure of nearly zero. We discuss a modified model, including new system parameter spaces such as porosity and permeability, which describe the plasma flow through the reactor.In Figure1, the chamber of the CVD apparatus is shown. It is modeled as a porous medium, where the porosity is given as the ratio between the void space and the total volume of the ceramic or bulk material:
(2.2)ϕ=VvVb,
where Vv is the void space and Vb is the volume of the ceramic or bulk material.Figure 1
Chamber of the CVD apparatus.For the one-phase model, the velocityv of the underlying fluid is calculated by Darcy's law and given by:
(2.3)v=-κμ∇pϕ,
where κ is the permeability, ∇p is the pressure gradient, ϕ is the porosity and μ the dynamic viscosity.To model a retarded gas concentration in the porous medium, we have to consider multiphases of the underlying fluid. The new model includes immobile and adsorbed phases, where the velocity of the fluid is zero and impermeable layers are given.In the model, we consider both absorption and adsorption taking place simultaneously and with given exchange rates. Therefore we consider the effect of the gas concentrations' being incorporated into the porous medium.We extend the model to two more phases:(i)
immobile phase,(ii)
adsorbed phase.In Figure2, the mobile and immobile phases of the gas concentration are shown in the macroscopic scale of the porous medium. Here the exchange rate between the mobile gas concentration and the immobile gas concentration control the flux to the medium.Figure 2
Mobile and immobile phase.In Figure3, the mobile and adsorbed phases of the gas concentration are shown in the macroscopic scale of the porous medium. To be more detailed in the mobile and immobile phases, where the gas concentrations can be adsorbed or absorbed, we consider a further phase. Here the adsorption in the mobile and immobile phase is treated as a retardation and given by a permeability in such layers.Mobile-adsorbed phase and immobile-adsorbed phase.
(a)(b)The model equations for the multiple phase equations are(2.4)ϕ∂tci+∇⋅Fi=g(-ci+ci,im)+kα(-ci+ci,ad)-λi,iϕci+∑k=k(i)λi,kϕck+Q̃i,inΩ×[0,t],(2.5)Fi=vci-De(i)∇ci+ηEci,(2.6)ϕ∂tci,im=g(ci-ci,im)+kα(ci,im,ad-ci,im)-λi,iϕci,im+∑k=k(i)λi,kϕck,im+Q̃i,im,inΩ×[0,t],(2.7)ϕ∂tci,ad=kα(ci-ci,ad)-λi,iϕci,ad+∑k=k(i)λi,kϕck,ad+Q̃i,ad,inΩ×[0,t],(2.8)ϕ∂tci,im,ad=kα(ci,im-ci,im,ad)-λi,iϕci,im,ad+∑k=k(i)λi,kϕck,im,ad+Q̃i,im,ad,inΩ×[0,t],(2.9)ci(x,t)=ci,0(x),ci,ad(x,t)=0,ci,im(x,t)=0,ci,im,ad(x,t)=0,onΩ,(2.10)ci(x,t)=ci,1(x,t),ci,ad(x,t)=0,ci,im(x,t)=0,ci,im,ad(x,t)=0,on∂Ω×[0,t],
where the  initial value is given as ci,0 and we assume a Dirichlet boundary conditions with the function ci,1(x,t) sufficiently smooth, all other initial and boundary conditions of the other phases are zero.The parameters are given as:ϕ: effective porosity [-],ci: concentration of the ith gaseous species in the plasma chamber,ci,im: concentration of the ith gaseous species in the immobile zones of the plasma chamber phase [mol/cm3],v: velocity through the chamber and porous substrate [20] [cm/nsec],De(i): element-specific diffusions-dispersions tensor [cm2/nsec],λi,i: decay constant of the ith species [1/nsec],Q̃i: source term of the ith species [mol/(cm3nsec)],g: exchange rate between the mobile and immobile concentration [1/nsec],kα: exchange rate between the mobile and adsorbed concentration or immobile and immobile adsorbed concentration (kinetic controlled sorption) [1/nsec],E: stationary electricfield in the apparatus [V/cm],η: the mobility rate in the electricfield; see [13] [cm2/nsec],withi=1,…,M and M denotes the number of components.The parameters in (2.4) are further described; see also [12].The four phases are treated in the full domain, such that we have a full coupling in time and space.The effective porosity is denoted byϕ and declares the portion of the porosities of the aquifer that is filled with plasma, and we assume a nearly fluid phase. The transport term is indicated by the Darcy velocity v, that presents the flow direction and the absolute value of the plasma flux. The velocity field is divergence free. The decay constant of the ith species is denoted by λi. Thereby, k(i) denotes the indices of the other species which react with species i.


### 2.3. Electric Field (Distribution)

In addition, in order to find the distribution of the electric field, it is necessary to solve the Poisson equation which relates the divergence of the local electric fields to the charge density [13], and is given by
(2.11)∇⋅E=eϵ0∑i=1nZici,
where e is the charge of the electron, ϵ0 the permittivity of the gas mixture, ci the particle density of species i and Zi the relative charge of species i; see [21].A general model for PECVD reactors is very complicated and we therefore simplify the model with some assumptions to obtain tractable problems.(i)
Predefined electric fieldE respecting different predefined areas Ωj, we apply:
(2.12)∇⋅E≈eϵ0∑i=1nZicconst,i|Ωj,
where we assume that ci≈cconst,iforΩj is a constant particle density of species i in domain Ωj with j∈1,…,m and m is the number of domains.(ii)
The plasma chemistry in the reactor is neglected.


### 2.4. Multiphysics Model

We derive the multiphysics model based on the multiphase and electric field model.We start by developing the multiphysics model in the following steps:(i)
transport system (multiphase model),(ii)
electrical field is approximated by a permeable layer (field model),(iii)
multiphysics model with multiphase transport and embedded electrical field.We assume that the number of gas particles in the plasma atmosphere Ar (Argon) is low and the density of the particles can be treated as nearly constant.Therefore, we can derive a constant field in each sector of the apparatus(2.13)E≈EΩj,
where EΩj is constant and Ωj⊂Ω are sectors in the domain.Such an electric field, see Figure4, can be assumed in near-field experiments to be nearly homogeneous,  where the field intensity is decreasing and can be assumed to be constant in small areas; see [4].Figure 4
Electrical sectors in the domain of the apparatus.Next the idea is to embed the electric field into a porous medium model. A simplified model can be derived with given domain-dependent parameters and a retardation factor that includes the influence of the electric field.First, the mobile phase equation (2.4) can be reformulated with retardation factor given by: (2.14)Ri∂∂tci+∇F̃i-RiRg,i=0,inΩ×[0,t],F̃i=vci-De(i)∇ci,(2.15)ci(x,t)=ci,0(x),onΩ,ci(x,t)=ci,1(x,t),on∂Ω×[0,t],
where ci is the particle density and F̃i the flux of species i. Rg,i is the reaction term of species i, meaning the right-hand side of (2.4), all other parameters as in (2.4)–(2.8). Ri is defined as a specific and constant Henry isotherm, called the retardation factor; see [22]. Such factors are known as sorption coefficients and model the rate between the mobile and immobile concentration, they are derived from experiments; see [23].In a second step, we calculate the retardation factor, which is based on the influence of the electric field.We subtract (2.4) from (2.14) and obtain (2.16)∇Fi=∇F̃iRi.
We include the assumption of the electric field dependence of the domain parts and consider (2.17)(1Ri-1)∇(v-De(i)∇)ci≈eϵ0∑k=1nZkcconst,k|Ωj,forj=1,…,m,
so we obtain the derivation of Ri for a specific domain Ωj⊂Ω and we have: (2.18)Ri=(eϵ0∑k=1nZkcconst,k|Ωj1∇(v-De(i)∇)ci+1)-1,forΩj.
In the third step we add the equations (2.6)–(2.8) to the modified mobile phase equation (2.14) and we obtain the full coupled system.Here the novel contribution is to derive a full coupled system of a multiphase and multiflow problem to include all the physical processes in the chamber. Such ideas have been developed and considered in fluid dynamical models for many years [24].Remark 2.1.
For the flow through the chamber, for which we assume a homogeneous medium and nonreactive plasma, we have considered a constant flow [1]. A further simplification is given by the very small porous substrate, for which we can assume the underlying velocity to a first approximation as constant [4].Remark 2.2.
For an in-stationary medium and reactive or ionized plasma, we have to take into account the relations of the electrons in thermal equilibrium. Such spatial variation can be considered by modeling the electron drift. Such modeling of the ionized plasma is done with Boltzmanns relation, [3].


## 2.1. Flow Field

The conservation of momentum is given by (flow field: Navier-Stokes equation)(2.1)∂∂tv+v⋅∇v=-∇p,inΩ×[0,t],v(x,t)=v0(x),onΩ,v(x,t)=v1(x,t),on∂Ω×[0,t],
where v is the velocity field, p the pressure, v0 the initial velocity field and the position vector x=(x1,x2)t∈Ω⊂ℝ2,+. Furthermore, we assume that the flow is divergence free and the pressure is predefined.


## 2.2. Transport Systems (Multiphase Equations)

We model the ionized plasma as an underlying medium in the chamber with mobile and immobile phases. Here transport in the plasma with gaseous species contain of mobile and immobile concentrations, [3]. For such a heterogenous plasma, we applied our expertise in modeling multiphase transport through a porous medium; see [2].To amplify the modeling of the gaseous flow to the chamber which is filled with ionized plasma, we deal with the so-called far-field model based on a porous medium. Here the plasma can be modeled as a continuous flow [17], that has mobile and immobile phases; see [8].We assume nearly a vacuum and a diffusion-dominated process, derived from the Knudsen diffusion, [19]. In such viscous flow regimes, we deal with small Knudsen numbers and a pressure of nearly zero. We discuss a modified model, including new system parameter spaces such as porosity and permeability, which describe the plasma flow through the reactor.In Figure1, the chamber of the CVD apparatus is shown. It is modeled as a porous medium, where the porosity is given as the ratio between the void space and the total volume of the ceramic or bulk material:
(2.2)ϕ=VvVb,
where Vv is the void space and Vb is the volume of the ceramic or bulk material.Figure 1
Chamber of the CVD apparatus.For the one-phase model, the velocityv of the underlying fluid is calculated by Darcy's law and given by:
(2.3)v=-κμ∇pϕ,
where κ is the permeability, ∇p is the pressure gradient, ϕ is the porosity and μ the dynamic viscosity.To model a retarded gas concentration in the porous medium, we have to consider multiphases of the underlying fluid. The new model includes immobile and adsorbed phases, where the velocity of the fluid is zero and impermeable layers are given.In the model, we consider both absorption and adsorption taking place simultaneously and with given exchange rates. Therefore we consider the effect of the gas concentrations' being incorporated into the porous medium.We extend the model to two more phases:(i)
immobile phase,(ii)
adsorbed phase.In Figure2, the mobile and immobile phases of the gas concentration are shown in the macroscopic scale of the porous medium. Here the exchange rate between the mobile gas concentration and the immobile gas concentration control the flux to the medium.Figure 2
Mobile and immobile phase.In Figure3, the mobile and adsorbed phases of the gas concentration are shown in the macroscopic scale of the porous medium. To be more detailed in the mobile and immobile phases, where the gas concentrations can be adsorbed or absorbed, we consider a further phase. Here the adsorption in the mobile and immobile phase is treated as a retardation and given by a permeability in such layers.Mobile-adsorbed phase and immobile-adsorbed phase.
(a)(b)The model equations for the multiple phase equations are(2.4)ϕ∂tci+∇⋅Fi=g(-ci+ci,im)+kα(-ci+ci,ad)-λi,iϕci+∑k=k(i)λi,kϕck+Q̃i,inΩ×[0,t],(2.5)Fi=vci-De(i)∇ci+ηEci,(2.6)ϕ∂tci,im=g(ci-ci,im)+kα(ci,im,ad-ci,im)-λi,iϕci,im+∑k=k(i)λi,kϕck,im+Q̃i,im,inΩ×[0,t],(2.7)ϕ∂tci,ad=kα(ci-ci,ad)-λi,iϕci,ad+∑k=k(i)λi,kϕck,ad+Q̃i,ad,inΩ×[0,t],(2.8)ϕ∂tci,im,ad=kα(ci,im-ci,im,ad)-λi,iϕci,im,ad+∑k=k(i)λi,kϕck,im,ad+Q̃i,im,ad,inΩ×[0,t],(2.9)ci(x,t)=ci,0(x),ci,ad(x,t)=0,ci,im(x,t)=0,ci,im,ad(x,t)=0,onΩ,(2.10)ci(x,t)=ci,1(x,t),ci,ad(x,t)=0,ci,im(x,t)=0,ci,im,ad(x,t)=0,on∂Ω×[0,t],
where the  initial value is given as ci,0 and we assume a Dirichlet boundary conditions with the function ci,1(x,t) sufficiently smooth, all other initial and boundary conditions of the other phases are zero.The parameters are given as:ϕ: effective porosity [-],ci: concentration of the ith gaseous species in the plasma chamber,ci,im: concentration of the ith gaseous species in the immobile zones of the plasma chamber phase [mol/cm3],v: velocity through the chamber and porous substrate [20] [cm/nsec],De(i): element-specific diffusions-dispersions tensor [cm2/nsec],λi,i: decay constant of the ith species [1/nsec],Q̃i: source term of the ith species [mol/(cm3nsec)],g: exchange rate between the mobile and immobile concentration [1/nsec],kα: exchange rate between the mobile and adsorbed concentration or immobile and immobile adsorbed concentration (kinetic controlled sorption) [1/nsec],E: stationary electricfield in the apparatus [V/cm],η: the mobility rate in the electricfield; see [13] [cm2/nsec],withi=1,…,M and M denotes the number of components.The parameters in (2.4) are further described; see also [12].The four phases are treated in the full domain, such that we have a full coupling in time and space.The effective porosity is denoted byϕ and declares the portion of the porosities of the aquifer that is filled with plasma, and we assume a nearly fluid phase. The transport term is indicated by the Darcy velocity v, that presents the flow direction and the absolute value of the plasma flux. The velocity field is divergence free. The decay constant of the ith species is denoted by λi. Thereby, k(i) denotes the indices of the other species which react with species i.


## 2.3. Electric Field (Distribution)

In addition, in order to find the distribution of the electric field, it is necessary to solve the Poisson equation which relates the divergence of the local electric fields to the charge density [13], and is given by
(2.11)∇⋅E=eϵ0∑i=1nZici,
where e is the charge of the electron, ϵ0 the permittivity of the gas mixture, ci the particle density of species i and Zi the relative charge of species i; see [21].A general model for PECVD reactors is very complicated and we therefore simplify the model with some assumptions to obtain tractable problems.(i)
Predefined electric fieldE respecting different predefined areas Ωj, we apply:
(2.12)∇⋅E≈eϵ0∑i=1nZicconst,i|Ωj,
where we assume that ci≈cconst,iforΩj is a constant particle density of species i in domain Ωj with j∈1,…,m and m is the number of domains.(ii)
The plasma chemistry in the reactor is neglected.


## 2.4. Multiphysics Model

We derive the multiphysics model based on the multiphase and electric field model.We start by developing the multiphysics model in the following steps:(i)
transport system (multiphase model),(ii)
electrical field is approximated by a permeable layer (field model),(iii)
multiphysics model with multiphase transport and embedded electrical field.We assume that the number of gas particles in the plasma atmosphere Ar (Argon) is low and the density of the particles can be treated as nearly constant.Therefore, we can derive a constant field in each sector of the apparatus(2.13)E≈EΩj,
where EΩj is constant and Ωj⊂Ω are sectors in the domain.Such an electric field, see Figure4, can be assumed in near-field experiments to be nearly homogeneous,  where the field intensity is decreasing and can be assumed to be constant in small areas; see [4].Figure 4
Electrical sectors in the domain of the apparatus.Next the idea is to embed the electric field into a porous medium model. A simplified model can be derived with given domain-dependent parameters and a retardation factor that includes the influence of the electric field.First, the mobile phase equation (2.4) can be reformulated with retardation factor given by: (2.14)Ri∂∂tci+∇F̃i-RiRg,i=0,inΩ×[0,t],F̃i=vci-De(i)∇ci,(2.15)ci(x,t)=ci,0(x),onΩ,ci(x,t)=ci,1(x,t),on∂Ω×[0,t],
where ci is the particle density and F̃i the flux of species i. Rg,i is the reaction term of species i, meaning the right-hand side of (2.4), all other parameters as in (2.4)–(2.8). Ri is defined as a specific and constant Henry isotherm, called the retardation factor; see [22]. Such factors are known as sorption coefficients and model the rate between the mobile and immobile concentration, they are derived from experiments; see [23].In a second step, we calculate the retardation factor, which is based on the influence of the electric field.We subtract (2.4) from (2.14) and obtain (2.16)∇Fi=∇F̃iRi.
We include the assumption of the electric field dependence of the domain parts and consider (2.17)(1Ri-1)∇(v-De(i)∇)ci≈eϵ0∑k=1nZkcconst,k|Ωj,forj=1,…,m,
so we obtain the derivation of Ri for a specific domain Ωj⊂Ω and we have: (2.18)Ri=(eϵ0∑k=1nZkcconst,k|Ωj1∇(v-De(i)∇)ci+1)-1,forΩj.
In the third step we add the equations (2.6)–(2.8) to the modified mobile phase equation (2.14) and we obtain the full coupled system.Here the novel contribution is to derive a full coupled system of a multiphase and multiflow problem to include all the physical processes in the chamber. Such ideas have been developed and considered in fluid dynamical models for many years [24].Remark 2.1.
For the flow through the chamber, for which we assume a homogeneous medium and nonreactive plasma, we have considered a constant flow [1]. A further simplification is given by the very small porous substrate, for which we can assume the underlying velocity to a first approximation as constant [4].Remark 2.2.
For an in-stationary medium and reactive or ionized plasma, we have to take into account the relations of the electrons in thermal equilibrium. Such spatial variation can be considered by modeling the electron drift. Such modeling of the ionized plasma is done with Boltzmanns relation, [3].


## 3. Physical Experiments

The base of the experimental setup is the plasma reactor chamber of a NISTGEC reference cell. The spiral antenna of a hybrid ICP/CCP-RF plasma source was replaced by a double spiral antenna [25]. This reduces the asymmetry of the magnetic field due to the superposition of the induced fields of both antennas. Also, the power coupling to the plasma increases and enhances the efficiency of the source. A set of MKS mass flow controllers allows any defined mixture of gaseous precursors. Even the flows of liquid precursors with high vapour pressure is controlled by this system. All other liquid and all solid precursors will be directly transported to the chamber by a controlled carrier gas flow. Besides the precursor flow, the density can also be changed by varying the pressure inside the recipient. Controlling the pressure is achieved with a valve between the recipient and vacuum pumps. In addition, a heated and insulated substrate holder was mounted. Thus, a temperature up to 800°C and a bias voltage can be applied to the substrate. While the pressure and RF power determine the undirected particle energy (plasma temperature), the bias voltage adds, only to the charged particles, energy directed at the substrate. Apart from the pressure and RF power control, the degree of ionization and number as well as the size of the molecular fractions can be controlled.Altogether, this setup provides as free process parameters:(i)
pressure (typically10-1–10-2 mbar),(ii)
precursor composition (TMS, TMS + H2, TMS + O2),(iii)
precursor flow rate (ranging from SCCM to SLM),(iv)
RF Power (up to 1100 W),(v)
substrate temperature (RT—800°C),(vi)
bias voltage (DC, unipolar and bipolar pulsed, floating).During all experiments, the process was observed by optical emission spectroscopy (OES) and mass spectroscopy (MS). The stoichiometry of the deposited films was analyzedex situ in a scanning electron microscope (SEM) by energy dispersive X-ray analysis (EDX).


### 3.1. Realization of Physical Experiments

The following parameters were used for the physical experiments. Such a reduction allows of concentrating on the important flow and transport processes in the gas phase. Further,  we apply the underlying mathematical model (far-field model, see Section2.2) so that we can switch between physical and mathematical parameters (see Table 1).(1)
Precursor: Tetramethylsilan (TMS).(2)
Substrate: VA-Steel.(3)
Film at substrate:SiCx.We apply the parameters  in Table2 for interpolation of the substrate temperature.For the substrate temperature and plasma power, we use the parameters in Table3.Table 1
The physical experiments with the real-life apparatus.
TestPR (mbar)ϑS (C)Pplasma (W)ϕ (TMS) (SCCM)ϕ(H2) (SCCM)Ratio (C : Si)Mass (growth) (g)Time (min)080701-01-VA9.7E-240090010.2300.978110.00012120080718-01-VA1.1E-140090010.0001.001740.00050130080718-02-VA4.5E-240090010.0001.248110.00070110080618-01-VA4.3E-240050010.2301.32078—127080716-01-VA1.1E-140050010.0001.425440.00250120080715-02-VA1.1E-140010010.0001.588720.00337122080804-01-VA4.5E-240010010.0002.915450.00356129080630-01-VA9.9E-280090010.2301.091160.00102120080807-01-VA4.5E-280090010.0001.180780.00118120080625-01-VA3.9E-280050010.2301.06373—120080626-01-VA9.3E-280050010.2301.128180.00174130080806-01-VA4.8E-280010010.0001.739130.00219121080715-01-VA1.1E-180010010.0001.624670.00234120081016-01-VA1.0E-160030010.0001.728980.00321123081020-01-VA1.1E-160030010.00501.490750.00249114081028-01-VA1.1E-160030010.00151.535490.00273120081023-01-VA1.1E-160030010.00101.542780.00312127081027-01-VA1.1E-160030010.005.51.558180.00277126081024-01-VA1.1E-160030010.003.51.643670.00299120081022-01-VA1.0E-160030010.002.51.695890.00318127Table 2Temperature (C)Ratio (SiC : C)4002.4 : 16001.5 : 17001.211 : 18001.1 : 1Table 3Temperature (C)Power (W)Ratio (SiC : C)4009001 : 0.974005001.3 : 18009001.18 : 1Remark 3.1.
In the process, the temperature and power of the plasma is important and experiments show that these are significant parameters. Based on these parameters, we initialize the mathematical model and interpolate the flux and reaction parts.


## 3.1. Realization of Physical Experiments

The following parameters were used for the physical experiments. Such a reduction allows of concentrating on the important flow and transport processes in the gas phase. Further,  we apply the underlying mathematical model (far-field model, see Section2.2) so that we can switch between physical and mathematical parameters (see Table 1).(1)
Precursor: Tetramethylsilan (TMS).(2)
Substrate: VA-Steel.(3)
Film at substrate:SiCx.We apply the parameters  in Table2 for interpolation of the substrate temperature.For the substrate temperature and plasma power, we use the parameters in Table3.Table 1
The physical experiments with the real-life apparatus.
TestPR (mbar)ϑS (C)Pplasma (W)ϕ (TMS) (SCCM)ϕ(H2) (SCCM)Ratio (C : Si)Mass (growth) (g)Time (min)080701-01-VA9.7E-240090010.2300.978110.00012120080718-01-VA1.1E-140090010.0001.001740.00050130080718-02-VA4.5E-240090010.0001.248110.00070110080618-01-VA4.3E-240050010.2301.32078—127080716-01-VA1.1E-140050010.0001.425440.00250120080715-02-VA1.1E-140010010.0001.588720.00337122080804-01-VA4.5E-240010010.0002.915450.00356129080630-01-VA9.9E-280090010.2301.091160.00102120080807-01-VA4.5E-280090010.0001.180780.00118120080625-01-VA3.9E-280050010.2301.06373—120080626-01-VA9.3E-280050010.2301.128180.00174130080806-01-VA4.8E-280010010.0001.739130.00219121080715-01-VA1.1E-180010010.0001.624670.00234120081016-01-VA1.0E-160030010.0001.728980.00321123081020-01-VA1.1E-160030010.00501.490750.00249114081028-01-VA1.1E-160030010.00151.535490.00273120081023-01-VA1.1E-160030010.00101.542780.00312127081027-01-VA1.1E-160030010.005.51.558180.00277126081024-01-VA1.1E-160030010.003.51.643670.00299120081022-01-VA1.0E-160030010.002.51.695890.00318127Table 2Temperature (C)Ratio (SiC : C)4002.4 : 16001.5 : 17001.211 : 18001.1 : 1Table 3Temperature (C)Power (W)Ratio (SiC : C)4009001 : 0.974005001.3 : 18009001.18 : 1Remark 3.1.
In the process, the temperature and power of the plasma is important and experiments show that these are significant parameters. Based on these parameters, we initialize the mathematical model and interpolate the flux and reaction parts.


## 4. Discretization and Solver Methods

We distinguish between the mobile and immobile phases. Here the mobile phases are parabolic partial differential equations and the immobile phases are ordinary differential equations. So for the space-discretization of the PDE's we apply finite volume methods as mass-conserved discretization schemes and for the time-discretization of the PDE's and ODE's we apply Runge-Kutta methods.So first we discretize in space, while we then apply a splitting method, taking into account the different matrices.


### 4.1. Spatial Discretization

We deal with the transport part of the mobile phase, see (2.14): (4.1)Ri∂∂tci+∇F̃i=0,inΩ×[0,t],F̃i=vci-De(i)∇ci,(4.2)ci(x,t)=ci,0(x),onΩ,ci(x,t)=ci,1(x,t),on∂Ω×[0,t],
For the convection part, we use a piecewise constant finite volume method with upwind discretization; see [26]. For the diffusion-dispersion part, we also apply a finite volume method and we assume the boundary values are denoted by n·De(i)∇ci(x,t)=0, where x∈Γ is the boundary Γ=∂Ω; compare with [27]. The initial conditions are given by ci(x,0)=ci,0(x).We integrate (4.1) over space and obtain (4.3)∫ΩjRi∂∂tcidx=∫Ωj∇⋅(-vci+De(i)∇ci)dx.
The time integration is done later in the decomposition method with implicit-explicit Runge-Kutta methods. Further the diffusion-dispersion term is lumped; compare with [12]. Equation (4.3) is discretized over space using Green's formula (4.4)VjRi∂∂tcidx=∫Γjn⋅(-vci+De(i)∇ci)dγ,
where Γj is the boundary of the finite volume cell Ωj and Vuj is the volume of the cell j. We use the approximation in space; see [12].The spatial integration for the diffusion part (4.4) is done by the mid-point rule over its finite boundaries and the convection part is done with a flux limiter and we obtain: (4.5)VjRi∂∂tci,j=∑e∈Λjne∇vciedγ+∑e∈Λj∑k∈Λje|Γjke|njke⋅Djke∇ci,jke,
where |Γjke| is the length of the boundary element Γjke. The gradients are calculated with the piecewise finite-element function ϕl.We decide to discretize the ux with an up-winding scheme and obtain the following discretization for the convection part:(4.6)Fj,e={vj,eci,j,ifvj,e≥0,vj,eci,k,ifvj,e<0,
where vj,e=∫ev·nj,eds.We obtain for the diffusion part:(4.7)∇cjke=∑l∈Λecl∇ϕl(xjke).
We get, using difference notation for the neighbour points j and l; compare with [28], the full semi-discretization: (4.8)VjRi∂∂tci,j=∑e∈ΛjFj,e+∑e∈Λj∑l∈Λe∖{j}(∑k∈Λje|Γjke|njke⋅Djke∇ϕl(xjke))(cj-cl),
where j=1,…,m.Remark 4.1.
For higher-order discretization of the convection equation, we apply a reconstruction which is based on Godunov's method. We apply a limiter function that fulfills the local min-max property. The method is explained in [26]. The linear polynomials are reconstructed by the element-wise gradient and are given by
(4.9)u(xj)=cj,∇u|Vj=1Vj∑e=1E∫Te∩Ωj∇cdx,withj=1,…,I.
The piecewise linear functions are denoted by
(4.10)ujk=cj+ψj∇u|Vj(xjk-xj),withj=1,…,I,
where ψj∈(0,1) is the limiter function and it fulfills the discrete minimum-maximum property, as described in [26].


### 4.2. Splitting Method to Couple Mobile, Immobile, and Adsorbed Parts

The motivation of the splitting method is the following observations.(i)
The mobile phase is semidiscretized with fast finite volume methods and can be stored in a stiffness matrix. We achieve large time steps, if we consider implicit Runge-Kutta methods of lower order (e.g., implicit Euler) as a time discretization method.(ii)
The immobile, adsorpted, and immobile-adsorpted phases are pure ordinary differential equations and each is cheap to solve with explicit Runge-Kutta schemes.(iii)
The ODEs can be seen as perturbations and are solved explicitly in the iterative scheme.For the full equation we use the following matrix notation:(4.11)∂tc=A1c+A2c+B1(c-cim)+B2(c-cad)+Q,∂tcim=A2cim+B1(cim-c)+B2(cim-cim,ad)+Qim,∂tcad=A2cad+B2(cad-c)+Qad,∂tcim,ad=A2cim,ad+B2(cim,ad-cim)+Qim,ad,
where c=(c1,…,cm)T is the spatial discretized concentration in the mobile phase, see (2.4), and cim=(c1,im,…,cm,im)T is the concentration in the immobile phase, and similarly for the other phase concentrations. A1 is the stiffness matrix of (2.4), A2 is the reaction matrix of the right-hand side of (2.4), B1 and B2 are diagonal matrices with entries the immobile and kinetic parameters, see (2.7) and (2.8).FurthermoreQ,…,Qim,ad are the spatially discretized source vectors.Now we have the following ordinary differential equation:(4.12)∂tC=(A1+A2+B1+B2-B1-B20-B1A2+B1+B20-B2-B20A2+B200-B20A2+B2)C+Q̃,
where C=(c,cim,cad,cim,ad)T and the right-hand side is Q̃=(Q,Qim,Qad,Qim,ad)T.For such an equation we apply the decomposition of the matrices(4.13)∂tC=ÃC+Q̃,∂tC=Ã1C+Ã2C+Q̃,
where (4.14)Ã1=(A1+A20000A20000A20000A2),Ã2=(B1+B2-B1-B20-B1B1+B20-B2-B20B200-B20B2).
This system of equations is numerically solved by an iterative scheme.Algorithm 4.2.
We divide our time interval[0,T] into subintervals [tn,tn+1], where n=0,1,…N, t0=0 and tN=T.
We start withn=0.
(1) The initial conditions are given byC0(tn+1)=C(tn). We start with k=0.
(2) Compute the fixed-point iteration scheme given by(4.15)∂tCk=Ã1Ck+Ã2Ck-1+Q̃,
where k is the iteration index; see [29]. For the time integration, we apply Runge-Kutta methods as ODE solvers; see [30, 31].
(3) The stop criterion for the time intervaltn,tn+1 is given by
(4.16)‖Ck(tn+1)-Ck-1(tn+1)‖≤err⁡,
where ∥·∥ is the maximum norm over all components of the solution vector. err is a given error tolerance, for example, err=10-4.
If (4.16) is fulfilled, we have the result
(4.17)C(tn+1)=Ck(tn+1).
Ifn=N (which means we reach the maximum iterative steps), then we stop and are done.
If (4.16) is not fulfilled, we do k=k+1 and go to (2).The error analysis of the schemes is given in the following theorem.Theorem 4.3.
LetA,B∈ℒ(X) be given linear bounded operators in a Banach space ℒ(X). We consider the abstract Cauchy problem
(4.18)∂tC(t)=ÃC(t)+B̃C(t),tn≤t≤tn+1,C(tn)=Cn,forn=1,…,N,
where t1=0 and the final time is tN=T∈ℝ+. Then problem (4.18) has a unique solution. For finite steps with time size τn=tn+1-tn, the iteration (4.15) for k=1,2,…,q is consistent with an order of consistency 𝒪(τnq).Proof.
The outline of the proof is given in [15].


## 4.1. Spatial Discretization

We deal with the transport part of the mobile phase, see (2.14): (4.1)Ri∂∂tci+∇F̃i=0,inΩ×[0,t],F̃i=vci-De(i)∇ci,(4.2)ci(x,t)=ci,0(x),onΩ,ci(x,t)=ci,1(x,t),on∂Ω×[0,t],
For the convection part, we use a piecewise constant finite volume method with upwind discretization; see [26]. For the diffusion-dispersion part, we also apply a finite volume method and we assume the boundary values are denoted by n·De(i)∇ci(x,t)=0, where x∈Γ is the boundary Γ=∂Ω; compare with [27]. The initial conditions are given by ci(x,0)=ci,0(x).We integrate (4.1) over space and obtain (4.3)∫ΩjRi∂∂tcidx=∫Ωj∇⋅(-vci+De(i)∇ci)dx.
The time integration is done later in the decomposition method with implicit-explicit Runge-Kutta methods. Further the diffusion-dispersion term is lumped; compare with [12]. Equation (4.3) is discretized over space using Green's formula (4.4)VjRi∂∂tcidx=∫Γjn⋅(-vci+De(i)∇ci)dγ,
where Γj is the boundary of the finite volume cell Ωj and Vuj is the volume of the cell j. We use the approximation in space; see [12].The spatial integration for the diffusion part (4.4) is done by the mid-point rule over its finite boundaries and the convection part is done with a flux limiter and we obtain: (4.5)VjRi∂∂tci,j=∑e∈Λjne∇vciedγ+∑e∈Λj∑k∈Λje|Γjke|njke⋅Djke∇ci,jke,
where |Γjke| is the length of the boundary element Γjke. The gradients are calculated with the piecewise finite-element function ϕl.We decide to discretize the ux with an up-winding scheme and obtain the following discretization for the convection part:(4.6)Fj,e={vj,eci,j,ifvj,e≥0,vj,eci,k,ifvj,e<0,
where vj,e=∫ev·nj,eds.We obtain for the diffusion part:(4.7)∇cjke=∑l∈Λecl∇ϕl(xjke).
We get, using difference notation for the neighbour points j and l; compare with [28], the full semi-discretization: (4.8)VjRi∂∂tci,j=∑e∈ΛjFj,e+∑e∈Λj∑l∈Λe∖{j}(∑k∈Λje|Γjke|njke⋅Djke∇ϕl(xjke))(cj-cl),
where j=1,…,m.Remark 4.1.
For higher-order discretization of the convection equation, we apply a reconstruction which is based on Godunov's method. We apply a limiter function that fulfills the local min-max property. The method is explained in [26]. The linear polynomials are reconstructed by the element-wise gradient and are given by
(4.9)u(xj)=cj,∇u|Vj=1Vj∑e=1E∫Te∩Ωj∇cdx,withj=1,…,I.
The piecewise linear functions are denoted by
(4.10)ujk=cj+ψj∇u|Vj(xjk-xj),withj=1,…,I,
where ψj∈(0,1) is the limiter function and it fulfills the discrete minimum-maximum property, as described in [26].


## 4.2. Splitting Method to Couple Mobile, Immobile, and Adsorbed Parts

The motivation of the splitting method is the following observations.(i)
The mobile phase is semidiscretized with fast finite volume methods and can be stored in a stiffness matrix. We achieve large time steps, if we consider implicit Runge-Kutta methods of lower order (e.g., implicit Euler) as a time discretization method.(ii)
The immobile, adsorpted, and immobile-adsorpted phases are pure ordinary differential equations and each is cheap to solve with explicit Runge-Kutta schemes.(iii)
The ODEs can be seen as perturbations and are solved explicitly in the iterative scheme.For the full equation we use the following matrix notation:(4.11)∂tc=A1c+A2c+B1(c-cim)+B2(c-cad)+Q,∂tcim=A2cim+B1(cim-c)+B2(cim-cim,ad)+Qim,∂tcad=A2cad+B2(cad-c)+Qad,∂tcim,ad=A2cim,ad+B2(cim,ad-cim)+Qim,ad,
where c=(c1,…,cm)T is the spatial discretized concentration in the mobile phase, see (2.4), and cim=(c1,im,…,cm,im)T is the concentration in the immobile phase, and similarly for the other phase concentrations. A1 is the stiffness matrix of (2.4), A2 is the reaction matrix of the right-hand side of (2.4), B1 and B2 are diagonal matrices with entries the immobile and kinetic parameters, see (2.7) and (2.8).FurthermoreQ,…,Qim,ad are the spatially discretized source vectors.Now we have the following ordinary differential equation:(4.12)∂tC=(A1+A2+B1+B2-B1-B20-B1A2+B1+B20-B2-B20A2+B200-B20A2+B2)C+Q̃,
where C=(c,cim,cad,cim,ad)T and the right-hand side is Q̃=(Q,Qim,Qad,Qim,ad)T.For such an equation we apply the decomposition of the matrices(4.13)∂tC=ÃC+Q̃,∂tC=Ã1C+Ã2C+Q̃,
where (4.14)Ã1=(A1+A20000A20000A20000A2),Ã2=(B1+B2-B1-B20-B1B1+B20-B2-B20B200-B20B2).
This system of equations is numerically solved by an iterative scheme.Algorithm 4.2.
We divide our time interval[0,T] into subintervals [tn,tn+1], where n=0,1,…N, t0=0 and tN=T.
We start withn=0.
(1) The initial conditions are given byC0(tn+1)=C(tn). We start with k=0.
(2) Compute the fixed-point iteration scheme given by(4.15)∂tCk=Ã1Ck+Ã2Ck-1+Q̃,
where k is the iteration index; see [29]. For the time integration, we apply Runge-Kutta methods as ODE solvers; see [30, 31].
(3) The stop criterion for the time intervaltn,tn+1 is given by
(4.16)‖Ck(tn+1)-Ck-1(tn+1)‖≤err⁡,
where ∥·∥ is the maximum norm over all components of the solution vector. err is a given error tolerance, for example, err=10-4.
If (4.16) is fulfilled, we have the result
(4.17)C(tn+1)=Ck(tn+1).
Ifn=N (which means we reach the maximum iterative steps), then we stop and are done.
If (4.16) is not fulfilled, we do k=k+1 and go to (2).The error analysis of the schemes is given in the following theorem.Theorem 4.3.
LetA,B∈ℒ(X) be given linear bounded operators in a Banach space ℒ(X). We consider the abstract Cauchy problem
(4.18)∂tC(t)=ÃC(t)+B̃C(t),tn≤t≤tn+1,C(tn)=Cn,forn=1,…,N,
where t1=0 and the final time is tN=T∈ℝ+. Then problem (4.18) has a unique solution. For finite steps with time size τn=tn+1-tn, the iteration (4.15) for k=1,2,…,q is consistent with an order of consistency 𝒪(τnq).Proof.
The outline of the proof is given in [15].


## 5. Experiments for the Multiple Phase Model

In the following subsections, we present our experiments based on the mobile, immobile and adsorbed gaseous phases.We contribute ideas to obtain an optimal layer deposition, which is based on the PE-CVD process, while different additional phases are considered, for example, plasma and precursor media.The main contributions are an optimal collection of point sources, line sources or moving sources to cover the deposition area, with respect to the remainder concentration in the immobile and adsorbed phases.We simulate the deposition process with our boundary value solver algorithms and could deal with many different conditions that might be impossible for physical experiments. Such simulation results may benefit the physical experiment and give new ideas to optimize such deposition problems of a complicated physical process.


### 5.1. Verification of the Code with Analytical Solutions

In the verification of our software code UG; see [32], we apply the multiphase model (2.4) and (2.6) with the following parameters.We use ascending parameters for the retardation factors. The retardation factors are given byR1=16, R2=8, R3=4, R4. The reaction factors are given by λ1=0.4, λ2=0.3, λ3=0.2, λ0·0.The initial conditions for the mobile phase are given by(5.1)ci(x,0)={aix+bi,x∈(x1,x2),0,otherwise,
where the parameters are given by x1=0.0, x2=1.0, a1=1.0, b1=1.0, a2=-1.0, b2=1.0, a3=1.0, b3=0.0, a4=2.0, b4=1.0.The mobile-immobile exchange rate isg=0.5.The initial conditions for the immobile concentrations arecim(0)=(1,1,1,1)t.For the time integration method we apply a fourth-order Runge-Kutta method and we applyk=1,…,4 iterative steps for the solver scheme.The errors are given by(5.2)err⁡i,j=∑k=1N|ui,coupled,j-1(xk,t)-ui,coupled,j(xk,t)|,
where xk=(k-1)Δx, k=0,…,N, Δx=0.01, N=1000, and j=2,3,4,5.The solutions of the species are given in Figures5 and 6.Experiment with the coupled model for the mobile solutionscim, with j=2 and j=4 iteration steps at time point t=8.
(a)(b)(c)(d)Experiment with the coupled model for the immobile solutionscim, with j=2 and j=4 iteration steps at time point t=8.
(a)(b)Remark 5.1.
The iterative scheme reduced the amount of numerical work to couple the mobile and immoble equations. We obtain more accurate results by using 3-4 iterative steps. Such schemes help the flexibility in coupling model equations by reducing the amount of recoding needed. More accuracy can be simply achieved with more iterative steps.


### 5.2. Verification of the Model with Physical Experiments

For the next experiments, we have the following parameters of the model, discretization, and solver methods.We apply interpolation and regression methods to couple the physical parameters that were discussed in Section3, to the mathematical parameters, see Figure 7 and Table 4.Table 4
Physical and mathematical parameters.
Physical parameterMathematical parameterTemperature, pressure, powerVelocity, diffusion, reactionT, p, WV, D, λFigure 7
Coupling of physical and mathematical parameter space.Parameters of (2.4)–(2.8)
In Table5, we list the parameters for our simulation tool UG; see [32]. The software toolbox has a flexible user interface to allow a large number of numerical experiments and approximations using the known physical parameters.Table 5
Model parameters.
Densityρ=1.0Mobile porosityϕ=0.333Immobile porosity0.333DiffusionD=0.0Longitudinal dispersionαL=0.0Transversal dispersionαT=0.00Retardation factorR=100e-4 (Henry rate)Velocity fieldv=(0.0,-4.0·10-8)tDecay rate of species of 1st EXλAB=1·10-68Decay rate of species of 2nd EXλAB=2·10-8, λBNN=1·<10-68Decay rate of species of 3rd EXλAB=0.25·10-8, λCB=0.5·10-8Geometry (2D domain)Ω=[0,100]×[0,100]BoundaryNeumann boundary at top, left, and right boundariesOutflow boundary at bottom boundaryDiscretization Method
Finite-volume method of 2nd order is shown in Table6.Table 6
Spatial discretization parameters.
Spatial step sizeΔxmin=1.56, Δxmax=2.21Refined levels6LimiterSlope limiterTest functionsLinear test function reconstructed with neighbour gradientsTime Discretization Method
Crank-Nicolson method (2nd order).Solver Method
In Tables7 and 8, we deal with test examples which are approximations of physical experiments.Table 7
Spatial discretization parameters.
Initial time stepΔtinit=5·107Controlled time stepΔtmax=1.298·107,Δtmin=1.158·107Number of time steps time step control100, 80, 30, 25 time steps are controlled with the Courant numberCFLmax=1Table 8
Solver methods and their parameters.
SolverBiCGstab (Bi conjugate gradient method)PreconditionerGeometric multigrid methodSmootherGauss-Seidel method as smoothers for the multigrid methodBasic level0Initial gridUniform grid with 2 elementsMaximum level6Finest gridUniform grid with 8192 elementsNumerical Experiment
In the test example, we deal with the following reaction:(5.3)2SiC+4H⟶λSiC+CH4+Si.Here, we have physical experiments and approximate temperature parametersT=400,600,800.We  compute the SiC : C ratio at a given temperatureT=400,600,800 with the UG program and fit to the parameter λ.We use a Lagrangian formula to computeλ for the new temperatures T=500,700 and apply the ratio of the simulated new parameters (Table 9). These numerical results are used to initialise the physical experiments  in Table 12.Table 9
Computed and experimental fitted parameters from UG simulations.
Tλ (fitted)λ (interpolated)Ratio (SiC : C) (computed with UG)4001/2·10-8—2.4 : 1500—0.35·10-81.85 : 16001/4·10-8—1.5 : 1700—0.171·10-81.211 : 18001/8·10-8—1.1 : 1One Source
See Table10.
In Figure8, we present the concentration of one point source at (50,20) with number of time steps equal to 25.Table 10
Parameters of source concentration.
Point source at position(x,y)=(50,20)Starting point of source concentrationtstart=0.0End point of source concentrationtend=1108Amount of permanent source concentrationcsource=1.0Number of time steps25Table 11
Rate of concentration.
RateSiCsource,max:SiCtarget,max9·106:6.5·106=1.38Table 12
L2-norm of the solution with refinements and CPU time.
Max·levelL2·normMax·normElementsCPU·timeTime steps61.6232+097.5026+0781921.1 min1971.6222+097.7068+07327683.55 min1981.6155+097.6753+0713107215.7 min1991.6137+097.6267+0752428859.55 min19One point source at(50,20) with 25 time steps, full concentration is given in red and zero concentration with blue.
(a)(b)In Figure9, we show the deposition rates of one point source at (50,20), with number of time steps equal to 25.Figure 9
Deposition rates for one point source at(50,20), with 25 time steps (x-axis: time scale in (nsec), y-axis: deposition rate in (mol/nsec)).Remark 5.2.
In the numerical experiments, we could approximate the physical experiments and obtain the same deposition rates for the Si: C deposition. Therefore, the model allows of making predictions about future deposition rates with various parameters. Such numerical experiments can replace expensive physical experiments and allow of restricting them to more special cases.


### 5.3. Experiments with Full Model(α=4·10-14 and Immobile Rate g=8·10-14)

In this experiment, we have taken into account the full model and verify with the previous simpler model equations that are verified by physial experiments.Such novel experiments allow of giving prognoses of which direction of model modification can be important for future diagnostics of physical experiments.We have the following outline of the experiment.The exchange in the following experiments of the carbon (C) species between the mobile and immobile concentrations is very low, it is aboutg=8·10-14, we assume less activities in the plasma environment. Further the exchange between the mobile and adsorbed mobile concentrations is also very low, it is about α=4·10-14, also the exchange rates between the immobile and adsorbed immobile concentration is the same as in the mobile and adsorbed mobile phases.In previous  experiments; see [33], we obtained optimal deposition results by combining multiple point sources which can be moved in the spatial directions (moving sources). Further we could approximate the physical results, while using mobile and immobile phase models.In this experiment we increase the model to four phases in order to predict the delicate adsorped regions. Physicists called such regions lost deposition areas, where the deposition material vanishes into the apparatus; see [25].We taken into account 11 point sources at the positionsY=20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, and these sources are moving in the X direction in equal spatial steps with 15.The concentration of each source in each step is equal to 1. The starting movement ofX is given with 50→35→20→35→50→65→80→65→50.In Figure10, we present the experiment with 11 moving sources at the positions in four phases. We obtain high depositions in the mobile phase, while we have only marginal depositions in the immobile and adsorbed phases.Figure 10
mobile (lowermost), immobile (lower-middle), mobile-adsorbed (upper-middle) and immobile-adsorbed (uppermost) case of 11 moving sources moving with number of time step equal to 15 (left figures) and 25 (right figures), where full concentration is given in red and zero concentration in blue.Here the important observations to give a measurement of the deposition are the various deposition rates in the different phases.In Figure11, we present the deposition rate of mobile concentration of 11 moving sources moving in the X-direction in steps of 15: X moves from 50→35→20→35→50→65→80→65→50, with the number of time steps equal to 25.Deposition rates in four phases: mobile (a), immobile (b), mobile-adsorbed (c), immobile-adsorbed (d) concentration of 11 moving sources moving, (x-axis: time scale in (nsec), y-axis: deposition rate in (mol/nsec)).
(a)(b)(c)(d)In see Table12, we discuss the CPU time for the spatial refinements. We consider the L2-norm of the solution and obtain convergent results of the solutions.Remark 5.3.
With moving sources we gain improved deposition rates, this was also observed in [33]. Nevertheless the remaining concentration in the immobile and adsorbed phases are important for seeing the lost areas. In Figure 11, the maximum deposition rate of carbon in the mobile phase is 1.2·108[mol], while the maximum lost deposition rate in the other phases are at least 18000[mol]. So in the immobile and adsorbed phases we lost only 0.00018 percent concentration of the deposition material. This may seem to be neglectable, but we have considered only 1[sec] deposition time, while a full cycle of deposition could be 1000[sec]. Therefore we lost about 10% of the concentration in the apparatus, which is immense. Due to this fact, such models are important for foreseeing the percentage of lost concentration. One of the simulation results is the usage of higher amount of depositon material, that is, the higher amount of the underlying precursor gases are used to overcome the lost of concentrations in the apparatus. The second result is to consider more detailed models to predict the influence of the multiphases of the underlying gaseous species.


## 5.1. Verification of the Code with Analytical Solutions

In the verification of our software code UG; see [32], we apply the multiphase model (2.4) and (2.6) with the following parameters.We use ascending parameters for the retardation factors. The retardation factors are given byR1=16, R2=8, R3=4, R4. The reaction factors are given by λ1=0.4, λ2=0.3, λ3=0.2, λ0·0.The initial conditions for the mobile phase are given by(5.1)ci(x,0)={aix+bi,x∈(x1,x2),0,otherwise,
where the parameters are given by x1=0.0, x2=1.0, a1=1.0, b1=1.0, a2=-1.0, b2=1.0, a3=1.0, b3=0.0, a4=2.0, b4=1.0.The mobile-immobile exchange rate isg=0.5.The initial conditions for the immobile concentrations arecim(0)=(1,1,1,1)t.For the time integration method we apply a fourth-order Runge-Kutta method and we applyk=1,…,4 iterative steps for the solver scheme.The errors are given by(5.2)err⁡i,j=∑k=1N|ui,coupled,j-1(xk,t)-ui,coupled,j(xk,t)|,
where xk=(k-1)Δx, k=0,…,N, Δx=0.01, N=1000, and j=2,3,4,5.The solutions of the species are given in Figures5 and 6.Experiment with the coupled model for the mobile solutionscim, with j=2 and j=4 iteration steps at time point t=8.
(a)(b)(c)(d)Experiment with the coupled model for the immobile solutionscim, with j=2 and j=4 iteration steps at time point t=8.
(a)(b)Remark 5.1.
The iterative scheme reduced the amount of numerical work to couple the mobile and immoble equations. We obtain more accurate results by using 3-4 iterative steps. Such schemes help the flexibility in coupling model equations by reducing the amount of recoding needed. More accuracy can be simply achieved with more iterative steps.


## 5.2. Verification of the Model with Physical Experiments

For the next experiments, we have the following parameters of the model, discretization, and solver methods.We apply interpolation and regression methods to couple the physical parameters that were discussed in Section3, to the mathematical parameters, see Figure 7 and Table 4.Table 4
Physical and mathematical parameters.
Physical parameterMathematical parameterTemperature, pressure, powerVelocity, diffusion, reactionT, p, WV, D, λFigure 7
Coupling of physical and mathematical parameter space.Parameters of (2.4)–(2.8)
In Table5, we list the parameters for our simulation tool UG; see [32]. The software toolbox has a flexible user interface to allow a large number of numerical experiments and approximations using the known physical parameters.Table 5
Model parameters.
Densityρ=1.0Mobile porosityϕ=0.333Immobile porosity0.333DiffusionD=0.0Longitudinal dispersionαL=0.0Transversal dispersionαT=0.00Retardation factorR=100e-4 (Henry rate)Velocity fieldv=(0.0,-4.0·10-8)tDecay rate of species of 1st EXλAB=1·10-68Decay rate of species of 2nd EXλAB=2·10-8, λBNN=1·<10-68Decay rate of species of 3rd EXλAB=0.25·10-8, λCB=0.5·10-8Geometry (2D domain)Ω=[0,100]×[0,100]BoundaryNeumann boundary at top, left, and right boundariesOutflow boundary at bottom boundaryDiscretization Method
Finite-volume method of 2nd order is shown in Table6.Table 6
Spatial discretization parameters.
Spatial step sizeΔxmin=1.56, Δxmax=2.21Refined levels6LimiterSlope limiterTest functionsLinear test function reconstructed with neighbour gradientsTime Discretization Method
Crank-Nicolson method (2nd order).Solver Method
In Tables7 and 8, we deal with test examples which are approximations of physical experiments.Table 7
Spatial discretization parameters.
Initial time stepΔtinit=5·107Controlled time stepΔtmax=1.298·107,Δtmin=1.158·107Number of time steps time step control100, 80, 30, 25 time steps are controlled with the Courant numberCFLmax=1Table 8
Solver methods and their parameters.
SolverBiCGstab (Bi conjugate gradient method)PreconditionerGeometric multigrid methodSmootherGauss-Seidel method as smoothers for the multigrid methodBasic level0Initial gridUniform grid with 2 elementsMaximum level6Finest gridUniform grid with 8192 elementsNumerical Experiment
In the test example, we deal with the following reaction:(5.3)2SiC+4H⟶λSiC+CH4+Si.Here, we have physical experiments and approximate temperature parametersT=400,600,800.We  compute the SiC : C ratio at a given temperatureT=400,600,800 with the UG program and fit to the parameter λ.We use a Lagrangian formula to computeλ for the new temperatures T=500,700 and apply the ratio of the simulated new parameters (Table 9). These numerical results are used to initialise the physical experiments  in Table 12.Table 9
Computed and experimental fitted parameters from UG simulations.
Tλ (fitted)λ (interpolated)Ratio (SiC : C) (computed with UG)4001/2·10-8—2.4 : 1500—0.35·10-81.85 : 16001/4·10-8—1.5 : 1700—0.171·10-81.211 : 18001/8·10-8—1.1 : 1One Source
See Table10.
In Figure8, we present the concentration of one point source at (50,20) with number of time steps equal to 25.Table 10
Parameters of source concentration.
Point source at position(x,y)=(50,20)Starting point of source concentrationtstart=0.0End point of source concentrationtend=1108Amount of permanent source concentrationcsource=1.0Number of time steps25Table 11
Rate of concentration.
RateSiCsource,max:SiCtarget,max9·106:6.5·106=1.38Table 12
L2-norm of the solution with refinements and CPU time.
Max·levelL2·normMax·normElementsCPU·timeTime steps61.6232+097.5026+0781921.1 min1971.6222+097.7068+07327683.55 min1981.6155+097.6753+0713107215.7 min1991.6137+097.6267+0752428859.55 min19One point source at(50,20) with 25 time steps, full concentration is given in red and zero concentration with blue.
(a)(b)In Figure9, we show the deposition rates of one point source at (50,20), with number of time steps equal to 25.Figure 9
Deposition rates for one point source at(50,20), with 25 time steps (x-axis: time scale in (nsec), y-axis: deposition rate in (mol/nsec)).Remark 5.2.
In the numerical experiments, we could approximate the physical experiments and obtain the same deposition rates for the Si: C deposition. Therefore, the model allows of making predictions about future deposition rates with various parameters. Such numerical experiments can replace expensive physical experiments and allow of restricting them to more special cases.


## 5.3. Experiments with Full Model(α=4·10-14 and Immobile Rate g=8·10-14)

In this experiment, we have taken into account the full model and verify with the previous simpler model equations that are verified by physial experiments.Such novel experiments allow of giving prognoses of which direction of model modification can be important for future diagnostics of physical experiments.We have the following outline of the experiment.The exchange in the following experiments of the carbon (C) species between the mobile and immobile concentrations is very low, it is aboutg=8·10-14, we assume less activities in the plasma environment. Further the exchange between the mobile and adsorbed mobile concentrations is also very low, it is about α=4·10-14, also the exchange rates between the immobile and adsorbed immobile concentration is the same as in the mobile and adsorbed mobile phases.In previous  experiments; see [33], we obtained optimal deposition results by combining multiple point sources which can be moved in the spatial directions (moving sources). Further we could approximate the physical results, while using mobile and immobile phase models.In this experiment we increase the model to four phases in order to predict the delicate adsorped regions. Physicists called such regions lost deposition areas, where the deposition material vanishes into the apparatus; see [25].We taken into account 11 point sources at the positionsY=20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, and these sources are moving in the X direction in equal spatial steps with 15.The concentration of each source in each step is equal to 1. The starting movement ofX is given with 50→35→20→35→50→65→80→65→50.In Figure10, we present the experiment with 11 moving sources at the positions in four phases. We obtain high depositions in the mobile phase, while we have only marginal depositions in the immobile and adsorbed phases.Figure 10
mobile (lowermost), immobile (lower-middle), mobile-adsorbed (upper-middle) and immobile-adsorbed (uppermost) case of 11 moving sources moving with number of time step equal to 15 (left figures) and 25 (right figures), where full concentration is given in red and zero concentration in blue.Here the important observations to give a measurement of the deposition are the various deposition rates in the different phases.In Figure11, we present the deposition rate of mobile concentration of 11 moving sources moving in the X-direction in steps of 15: X moves from 50→35→20→35→50→65→80→65→50, with the number of time steps equal to 25.Deposition rates in four phases: mobile (a), immobile (b), mobile-adsorbed (c), immobile-adsorbed (d) concentration of 11 moving sources moving, (x-axis: time scale in (nsec), y-axis: deposition rate in (mol/nsec)).
(a)(b)(c)(d)In see Table12, we discuss the CPU time for the spatial refinements. We consider the L2-norm of the solution and obtain convergent results of the solutions.Remark 5.3.
With moving sources we gain improved deposition rates, this was also observed in [33]. Nevertheless the remaining concentration in the immobile and adsorbed phases are important for seeing the lost areas. In Figure 11, the maximum deposition rate of carbon in the mobile phase is 1.2·108[mol], while the maximum lost deposition rate in the other phases are at least 18000[mol]. So in the immobile and adsorbed phases we lost only 0.00018 percent concentration of the deposition material. This may seem to be neglectable, but we have considered only 1[sec] deposition time, while a full cycle of deposition could be 1000[sec]. Therefore we lost about 10% of the concentration in the apparatus, which is immense. Due to this fact, such models are important for foreseeing the percentage of lost concentration. One of the simulation results is the usage of higher amount of depositon material, that is, the higher amount of the underlying precursor gases are used to overcome the lost of concentrations in the apparatus. The second result is to consider more detailed models to predict the influence of the multiphases of the underlying gaseous species.


## 6. Conclusions and Discussions

We have presented a continuous model for the multiple phases, we assumed gaseous behaviour with exchange rates to adsorbed and immobile phases at very low-pressure and low temperature while dealing with catalyst processes, for example, plasma environment and precursor gases. We have to take into acount the remaining gas concentrations in each processes. Such detailed modelling allows of seeing the delicate retardation and sorption processes of the underlying plasma medium. From the methodology side of the numerical simulations, the contributions were to decouple the multiphase problem into single-phase problems, where each single problem can be solved with more accuracy. The iterative schemes allows of coupling the simpler equations and for each additional iterative step, we could reduce the splitting error. Such iterative methods allow of accelerating the solver process of multiphase problems.We can see in the numerical experiments a loss of the deposition concentration, which means that we have to consider a higher inflow rate of the species. Further, we could verify the underlying model with physical experiments. In the numerical experiments, we presented various ideas to control an optimal deposition process. While numerical experiments are cheap and all possible parameter variations of the model can be done in less than a few weeks, the help in reducing the need for expensive physical experiments is enormous. Such simulations reduce the need for physical experiments and allow of foreseeing new directions of helpful optimizations.In the future we are interested on analysing such fast processes due to very small time scales, for example, in a nanoscaled models (molecular dynamics model).



---

*Source: 641920-2011-03-09.xml*