# Model of PE-CVD Apparatus: Verification and Simulations

Authors: J. Geiser; V. Buck; M. Arab
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

In this paper, we present a simulation of chemical vapor deposition with metallic bipolar plates. In chemical vapor deposition, a delicate optimization between temperature, pressure and plasma power is important to obtain homogeneous deposition. The aim is to reduce the number of real-life experiments in a given CVD plasma reactor. Based on the large physical parameter space, there are a hugh number of possible experiments. A detailed study of the physical experiments in a CVD plasma reactor allows to reduce the problem to an approximate mathematical model, which is the underlying transport-reaction model. Significant regions of the CVD apparatus are approximated and physical parameters are transferred to the mathematical parameters. Such an approximation reduces the mathematical parameter space to a realistic number of numerical experiments. The numerical results are discussed with physical experiments to give a valid model for the assumed growth and we could reduce expensive physical experiments.


---

## Body

## 1. Introduction

We motivate our study by simulating the growth of a thin film by PE-CVD (plasma-enhanced chemical vapor deposition) processes; see [1, 2]. Such technical processes are very complex and real-life experiments are enormously extensive and expensive. Based on a large physical parameter space, the number of possible experiments involves at least the variation of all possible parameters. Such a large number of experiments can be reduced with the help of numerical experiments based on a mathematical model. Such modelling results are based on interdisciplinary work by engineers, mathematicians, and physicists. We derive a multiphysics model, that includes a simplification of the dominant physical processes, that is, transport of the reactive species in the gas phase and their deposition rates at the target layer.The approximation between the mathematical parameters and physical parametes can be done by a regression method, such that we can verify the physical experiments. Such approximations help to study the physical experiments with simulation tools which are cheaper and can predict the more appropriate experiments which should be done to understand and control the physical processes.In the following, we introduce the PE-CVD process and the important trends in modelling it.A gas exposed to an electric field under low pressure conditions(<5Torr) results in a nonequilibrium plasma; see [3, 4]. Such ionized media, known as “cold” plasma or glow discharges, are powerful surface-modification tools in Material Science and Technology. Low-pressure plasmas allow to modify the surface chemistry and properties of materials compatible with a low-medium vacuum, through a PE-CVD process; see the applications in [4, 5].Here, PE-CVD processes are attractive methods, because of their reproducible chemical processes that can be controlled by pressure, temperature, and by additional precursor gases. Such methods have been developed in recent years and there is an interest in producing high-temperature films; see [2].We consider models that are related to mesoscopic scales [6], with respect to flows close to the wafer surface, where a wafer is the target material (e.g., metal or ceramic) for deposition, [2]. We assume that the wafer is a homogeneous medium and that the surface can be modelled as a porous medium [7].The physical experiments are used to measure the influence of temperature, pressure and plasma power on deposition rates; see [8]. Here the plasma reactor chamber of an NIST GEC reference cell is used and for the hybrid ICP/CCP-RF plasma source, a double spiral antenna is used; see [8]. Such experiments are important but the variation of all possible parameters is very extensive.Mathematically, we apply an interpolation between physical and mathematical parameters to verify our simulation model. Based on the smaller mathematical parameter space, we can allow much more experiments and obtain via the regression function the resulting parameters to the physical experiments. Such switching between numerical experiments and physical experiments reduces the experiments to a possible amount and we can optimize the deposition process.The numerical results are discussed and applied to validation problems and real-life problems. We discuss an application involving the deposition ofSiC thin films on a metallic plate.The paper is organized as follows.In Section2, we present our mathematical model and a possible reduced model for further approximations. In Section 3, we discuss the physical experiments of the CVD process. The numerical methods of the transport-reaction equation and their parameter approximation to the physical model are described in Section 4. The numerical experiments are given in Section 5. In Section 6, we briefly summarize our results.


## 2. Mathematical Model

In the next, we discuss the derivation of our model.We start by developing a multiphase model in the following steps:(i)
standard transport model (one phase),(ii)
flow model (flow field of plasma medium),(iii)
multiphase model with mobile and immobile zones.In each part of the model, we can refine the transport processes of the deposition gaseous species or reaction gaseous species depending on the influence of flow field, plasma zones and precursor gases.A schematic test geometry of the CVD reactor is given in Figure1.Figure 1
Vertically impinging CVD reactor.


### 2.1. Standard Transport Model

In the following, the models are discussed in terms of far-field and near-field problems, which take into account the scales of the model.Two different types of model can be discussed:(1)
convection-diffusion-reaction equations [6] (far-field problem),(2)
Boltzmann-Lattice equations [9] (near-field problem).The modelling is governed by the Knudsen Number, whereby the Knudsen number is a dimensionless number and defines the ratio of the molecular mean free path length to a representative physical length scale.(2.1)Kn=λL,
where λ is the mean free path and L is the representative physical length scale. This length scale could be, for example, the radius of a body in a fluid. Here, we deal with small Knudsen Numbers Kn≈0.01-1.0 for a convection-diffusion-reaction equation and a constant velocity field, whereas for, large Knudsen Numbers Kn≥1.0, we deal with a Boltzmann equation [2]. From modelling the gaseous transport of the deposition species, we consider the pure far-field model and assume a continuum flow field; see [10].Such assumptions lead to transport equations that can be treated with a convection-diffusion-reaction equation owing to a constant velocity field, thus:(2.2)∂∂tc+∇F-Rg=q(x,t),inΩ×[0,t],F=vc-D∇c,(2.3)c(x,t)=c0(x),onΩ,c(x,t)=c1(x,t),on∂Ω×[0,t],
where c is the molar concentration of the reaction gases (called species) and F is the flux of the species. v is the flux velocity through the chamber and porous substrate [7]. D is the diffusion matrix and Rg is the reaction term. The initial value is given as c0 and we assume a Dirichlet boundary with the function c1(x,t) sufficiently smooth. q(x,t) is a source function depending on time and space and represents the inflow of the species.The parameters of the equation are derived as follows. Diffusion in the modified CVD process is given by Knudsen diffusion [11]. We consider the overall pressure in the reactor as 200 Pa and the substrate temperature (or wafer surface temperature) is about 600–900 K. The pore size in the homogeneous substrate is assumed to be 80 nm. The homogeneous substrate can be either a porous medium, for example, a ceramic material, see [11], or a dense plasma, assumed to be very dense and stationary; see [1]. For such media, we can derive diffusion based on Knudsen diffusion.The diffusion is described by(2.4)D=2ϵμKνr3RT,
where ϵ is the porosity, μK is the shape factor of Knudsen diffusion, r is the average pore radius, R and T are the gas constant and temperature, respectively, and ν is the mean molecular speed, given by(2.5)ν=8RTπW,
where W is the molar mass of the diffusive gas.For homogeneous reactions during the CVD process, we consider a constant reaction ofSi, Ti, and C species given as(2.6)3Ti+Si+2C→Ti3SiC2,
where Si3TiC2 is a MAX-phase material, see [12], which deposits on the wafer surface. For simplicity, we do not consider the intermediate reaction with the precursor gases [1] and we assume we are dealing with a compound gas 3Ti+Si+2C; see [13]. Therefore we can concentrate on the transport of one species.The reaction rate is then given by(2.7)λ=kr[3Si]M[Ti]N[2C]O[Si3TiC2]L,
where kr is the apparent reaction constant and L,M,N,O are the reaction orders of the reactants.A schematic overview of the one-phase model is presented in Figure2. Here, the gas chamber of the CVD apparatus is shown, which is modelled by a homogeneous medium.Figure 2
Gas chamber of the CVD apparatus.


### 2.2. Flow Field

The flow field is derived in which the velocity is used for the transport of species. The velocity in the homogeneous substrate is modelled by a porous medium [14, 15]. We assume a stationary or low reactivity medium, for example, nonionized or low-ionized plasma or less reactive precursor gas. Further, the pressure can be assumed to have the Maxwell distribution [1]:(2.8)p=ρbT,
where ρ is the density, b is Boltzmann's constant, and T is the temperature.The model equations are based on mass and momentum conservation equations, where we assume conservation of energy. Because of the low temperature and low pressure environment, we assume the gaseous flow has a nearly liquid behaviour. Therefore, the derivation of the velocity can be given by Darcy's law:(2.9)v=-kμ(∇p-ρg),
where  v  is the velocity of the fluid, k is the permeability tensor, μ is the dynamic viscosity, p is the pressure, g is the vector of gravity, and ρ is the density of the fluid.We use the continuum equation of the particle density and obtain the equation of the system, which is given as a flow equation:(2.10)∂t(ϕρ)+∇·(ρv)=Q,
where ρ is the unknown particle density, ϕ is the effective porosity, and Q is the source-term of the fluid. We assume a stationary fluid and consider only divergence-free velocity fields, that is,(2.11)∇·v(x)=0,x∈Ω.
The boundary conditions for the flow equation are given as(2.12)p=pr(t,γ),t>0,γ∈∂Ω,n·v=mf(t,γ),t>0,γ∈∂Ω,
where n is the normal unit vector with respect to ∂Ω, where we assume that the pressure pr and flow concentration mf are prescribed by Dirichlet boundary conditions [15].From the nearly stationary fluids, we assume that the conservation of momentum for velocityv  is given [15, 16]. Therefore, we can neglect the computation of the momentum for the velocity.Remark 2.1.
For flow through a gas chamber, for which we assume a homogeneous medium and nonreactive plasma, we have considered a constant flow [17]. A further simplification is given by the very small porous substrate, for which we can assume the underlying velocity in a first approximation as constant [2].Remark 2.2.
For a nonstationary media and reactive or ionized plasma, we have to take into account the relations for electrons in thermal equilibrium. Such a spatial variation can be considered by modelling electron drift. Such modelling of the ionized plasma is done with the Boltzman relation [1].


### 2.3. Multiphase Model: Mobile and Immobile Zones

More complicated processes such as retardation, adsorption, and dissipation processes, the gaseous species are modelled with multiphase equations. We take into account the fact that the concentration of species can be given in mobile and immobile versions, depending on their different reactive states; see [7]. From these behaviours, we have to model the transport and adsorbed state of species; see also Figure 3. Here, the mobile and immobile phases of the gas concentration are shown at a macroscopic scale in the porous medium.Figure 3
Mobile and immobile phases.The model equations are given as combinations of transport and reaction equations (coupled partial and ordinary differential equations) as(2.13)ϕ∂tciL+∇·(vciL-De(i)∇ciL)=g(-ciL+ci,imL)-λi,iϕciL+∑k=k(i)λi,kϕckL+Q̃i,(2.14)ϕ∂tci,imL=g(ciL-ci,imL)-λi,iϕci,imL+∑k=k(i)λi,kϕck,imL+Q̃i,im,
where i=1,…,M and M denotes the number of components.The parameters in (2.13) are further described; see also [18].The effective porosity is denoted byϕ and describes the portion of the porosity 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, which presents the flow direction and absolute value of the plasma flux. The velocity field is divergence-free. The decay constant of the ith species is denoted by λi. k(i) therefore denotes the indices of the other species.Remark 2.3.
The concentrations in the mobile zones is modelled with convection-diffusion-reaction equations, see also Section2.1, while the concentration in the immobile zones are modelled with reaction equations. These two phases represent the mobility of the gaseous species through homogeneous media, where the concentrations in the immobile zones are at least the lost amounts of depositable gases.


### 2.4. Simplified Model: Far-Field Model

We concentrate on a far-field model and assume continuum flow and that the transport equations can be treated by a convection-diffusion-reaction equation with a constant velocity field:(2.15)∂∂tc+∇F-Rg=0,inΩ×[0,t],F=vc-D∇c,c(x,t)=c0(x),onΩ,c(x,t)=c1(x,t),on∂Ω×[0,t],
where c is the molar concentration and F is the flux of the species. v is the flux velocity through the chamber and porous substrate [7]. D is the diffusion matrix and Rg is the reaction term. The initial value is given as c0 and we assume a Dirichlet boundary with the function c1(x,t) sufficiently smooth.Remark 2.4.
The focus on the dominant far field processes in the gas phase involving the reactive species reduces enormously the physical parameter space. Such a realistic reduction with respect to the experiments can also simplify the underlying mathematical model and concentrate on a defined number of experiments. Such experiments can validate the switching between physical and mathematical parameter space and allow to foresee the important processes in the gas phase.


## 2.1. Standard Transport Model

In the following, the models are discussed in terms of far-field and near-field problems, which take into account the scales of the model.Two different types of model can be discussed:(1)
convection-diffusion-reaction equations [6] (far-field problem),(2)
Boltzmann-Lattice equations [9] (near-field problem).The modelling is governed by the Knudsen Number, whereby the Knudsen number is a dimensionless number and defines the ratio of the molecular mean free path length to a representative physical length scale.(2.1)Kn=λL,
where λ is the mean free path and L is the representative physical length scale. This length scale could be, for example, the radius of a body in a fluid. Here, we deal with small Knudsen Numbers Kn≈0.01-1.0 for a convection-diffusion-reaction equation and a constant velocity field, whereas for, large Knudsen Numbers Kn≥1.0, we deal with a Boltzmann equation [2]. From modelling the gaseous transport of the deposition species, we consider the pure far-field model and assume a continuum flow field; see [10].Such assumptions lead to transport equations that can be treated with a convection-diffusion-reaction equation owing to a constant velocity field, thus:(2.2)∂∂tc+∇F-Rg=q(x,t),inΩ×[0,t],F=vc-D∇c,(2.3)c(x,t)=c0(x),onΩ,c(x,t)=c1(x,t),on∂Ω×[0,t],
where c is the molar concentration of the reaction gases (called species) and F is the flux of the species. v is the flux velocity through the chamber and porous substrate [7]. D is the diffusion matrix and Rg is the reaction term. The initial value is given as c0 and we assume a Dirichlet boundary with the function c1(x,t) sufficiently smooth. q(x,t) is a source function depending on time and space and represents the inflow of the species.The parameters of the equation are derived as follows. Diffusion in the modified CVD process is given by Knudsen diffusion [11]. We consider the overall pressure in the reactor as 200 Pa and the substrate temperature (or wafer surface temperature) is about 600–900 K. The pore size in the homogeneous substrate is assumed to be 80 nm. The homogeneous substrate can be either a porous medium, for example, a ceramic material, see [11], or a dense plasma, assumed to be very dense and stationary; see [1]. For such media, we can derive diffusion based on Knudsen diffusion.The diffusion is described by(2.4)D=2ϵμKνr3RT,
where ϵ is the porosity, μK is the shape factor of Knudsen diffusion, r is the average pore radius, R and T are the gas constant and temperature, respectively, and ν is the mean molecular speed, given by(2.5)ν=8RTπW,
where W is the molar mass of the diffusive gas.For homogeneous reactions during the CVD process, we consider a constant reaction ofSi, Ti, and C species given as(2.6)3Ti+Si+2C→Ti3SiC2,
where Si3TiC2 is a MAX-phase material, see [12], which deposits on the wafer surface. For simplicity, we do not consider the intermediate reaction with the precursor gases [1] and we assume we are dealing with a compound gas 3Ti+Si+2C; see [13]. Therefore we can concentrate on the transport of one species.The reaction rate is then given by(2.7)λ=kr[3Si]M[Ti]N[2C]O[Si3TiC2]L,
where kr is the apparent reaction constant and L,M,N,O are the reaction orders of the reactants.A schematic overview of the one-phase model is presented in Figure2. Here, the gas chamber of the CVD apparatus is shown, which is modelled by a homogeneous medium.Figure 2
Gas chamber of the CVD apparatus.


## 2.2. Flow Field

The flow field is derived in which the velocity is used for the transport of species. The velocity in the homogeneous substrate is modelled by a porous medium [14, 15]. We assume a stationary or low reactivity medium, for example, nonionized or low-ionized plasma or less reactive precursor gas. Further, the pressure can be assumed to have the Maxwell distribution [1]:(2.8)p=ρbT,
where ρ is the density, b is Boltzmann's constant, and T is the temperature.The model equations are based on mass and momentum conservation equations, where we assume conservation of energy. Because of the low temperature and low pressure environment, we assume the gaseous flow has a nearly liquid behaviour. Therefore, the derivation of the velocity can be given by Darcy's law:(2.9)v=-kμ(∇p-ρg),
where  v  is the velocity of the fluid, k is the permeability tensor, μ is the dynamic viscosity, p is the pressure, g is the vector of gravity, and ρ is the density of the fluid.We use the continuum equation of the particle density and obtain the equation of the system, which is given as a flow equation:(2.10)∂t(ϕρ)+∇·(ρv)=Q,
where ρ is the unknown particle density, ϕ is the effective porosity, and Q is the source-term of the fluid. We assume a stationary fluid and consider only divergence-free velocity fields, that is,(2.11)∇·v(x)=0,x∈Ω.
The boundary conditions for the flow equation are given as(2.12)p=pr(t,γ),t>0,γ∈∂Ω,n·v=mf(t,γ),t>0,γ∈∂Ω,
where n is the normal unit vector with respect to ∂Ω, where we assume that the pressure pr and flow concentration mf are prescribed by Dirichlet boundary conditions [15].From the nearly stationary fluids, we assume that the conservation of momentum for velocityv  is given [15, 16]. Therefore, we can neglect the computation of the momentum for the velocity.Remark 2.1.
For flow through a gas chamber, for which we assume a homogeneous medium and nonreactive plasma, we have considered a constant flow [17]. A further simplification is given by the very small porous substrate, for which we can assume the underlying velocity in a first approximation as constant [2].Remark 2.2.
For a nonstationary media and reactive or ionized plasma, we have to take into account the relations for electrons in thermal equilibrium. Such a spatial variation can be considered by modelling electron drift. Such modelling of the ionized plasma is done with the Boltzman relation [1].


## 2.3. Multiphase Model: Mobile and Immobile Zones

More complicated processes such as retardation, adsorption, and dissipation processes, the gaseous species are modelled with multiphase equations. We take into account the fact that the concentration of species can be given in mobile and immobile versions, depending on their different reactive states; see [7]. From these behaviours, we have to model the transport and adsorbed state of species; see also Figure 3. Here, the mobile and immobile phases of the gas concentration are shown at a macroscopic scale in the porous medium.Figure 3
Mobile and immobile phases.The model equations are given as combinations of transport and reaction equations (coupled partial and ordinary differential equations) as(2.13)ϕ∂tciL+∇·(vciL-De(i)∇ciL)=g(-ciL+ci,imL)-λi,iϕciL+∑k=k(i)λi,kϕckL+Q̃i,(2.14)ϕ∂tci,imL=g(ciL-ci,imL)-λi,iϕci,imL+∑k=k(i)λi,kϕck,imL+Q̃i,im,
where i=1,…,M and M denotes the number of components.The parameters in (2.13) are further described; see also [18].The effective porosity is denoted byϕ and describes the portion of the porosity 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, which presents the flow direction and absolute value of the plasma flux. The velocity field is divergence-free. The decay constant of the ith species is denoted by λi. k(i) therefore denotes the indices of the other species.Remark 2.3.
The concentrations in the mobile zones is modelled with convection-diffusion-reaction equations, see also Section2.1, while the concentration in the immobile zones are modelled with reaction equations. These two phases represent the mobility of the gaseous species through homogeneous media, where the concentrations in the immobile zones are at least the lost amounts of depositable gases.


## 2.4. Simplified Model: Far-Field Model

We concentrate on a far-field model and assume continuum flow and that the transport equations can be treated by a convection-diffusion-reaction equation with a constant velocity field:(2.15)∂∂tc+∇F-Rg=0,inΩ×[0,t],F=vc-D∇c,c(x,t)=c0(x),onΩ,c(x,t)=c1(x,t),on∂Ω×[0,t],
where c is the molar concentration and F is the flux of the species. v is the flux velocity through the chamber and porous substrate [7]. D is the diffusion matrix and Rg is the reaction term. The initial value is given as c0 and we assume a Dirichlet boundary with the function c1(x,t) sufficiently smooth.Remark 2.4.
The focus on the dominant far field processes in the gas phase involving the reactive species reduces enormously the physical parameter space. Such a realistic reduction with respect to the experiments can also simplify the underlying mathematical model and concentrate on a defined number of experiments. Such experiments can validate the switching between physical and mathematical parameter space and allow to foresee the important processes in the gas phase.


## 3. Physical Experiments

The base of the experimental setup is the plasma reactor chamber of an NISTGEC reference cell. The spiral antenna of a hybrid ICP/CCP-RF plasma source was replaced by a double spiral antenna [8]. 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 allow any defined mixture of gaseous precursors. Even the flows of liquid precursors with high vapor pressure are 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 size of 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).Realisation of Physical Experiments
The following parameters are used for the physical experiments. Such a reduction allows to concentrate on the important flow and transport processes in the gas phase. Further, we apply the underlying mathematical model (far field model, see Section2.4) such that we can switch between physical and mathematical parameters.Precursor: Tetramethylsilan (TMS),Substrate: VA-Steel,Film at substrate:SiCx(see Table 1).Table 1TestPR [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.32078127080716-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.06373120080626-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.00318127We apply the parameters shown in Table2 for interpolation of the substrate temperature.Table 2TemperatureRatio (SiC:C)4002.4: 16001.5: 17001.211: 18001.1: 1For the substrate temperature and plasma power, see Table3.Table 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 are 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. Numerical Methods

In this section, we discuss the numerical methods. To accelerate our numerical methods, we combined the numerical and analytical parts in the solver processes.


### 4.1. Discretization and Solver Methods

For space-discretization of the PDEs, we apply finite-volume methods in mass conserved discretization schemes and for time-discretization of the resulting ODEs we apply Runge-Kutta methods or BDF methods. To accelerate the solver process, we combine the numerical and analytical parts of the solutions.


#### 4.1.1. Discretization Method for Convection Equation

We deal with the following convection equation(4.1)∂tRc-v·∇c=0,
where R is the retardation factor and presents the retention of the concentration; see also (2.2). v is the velocity. We have a simple boundary condition c=0 for the inflow and outflow boundary and the initial values are given as c(xj,0)=cj0(x). We use a piecewise constant discretization method with the upwind discretization done in  [19] and get(4.2)VjRcjn+1=VjRcjn-τn∑k∈out(j)vjkcjn+τn∑k∈in(j)cknvkj,VjRcjn+1=cjn(RVj-τnνj)+τn∑k∈in(j)cknvkj.
The explicit time discretization has to fulfill the discrete minimum-maximum property  [19] and we get the following restriction on time steps:(4.3)τj=RVjνj,τn≤minj=1,…,Iτj.
To obtain improved spatial discretization methods and apply larger time steps, we introduce a reconstruction with linear polynomials as a higher test-function in the next subsection.


#### 4.1.2. Discretization Method for Convection-Reaction Equation Based on Embedded One-Dimensional Analytical Solutions

We apply Godunov's discretization method, compare for example,  [20], and extend the formulation to the analytical solution of convection-reaction equations. We reduce the multidimensional equation to one-dimensional equations and solve each equation exactly. The one-dimensional solution is multiplied by the underlying volume to give the mass formulation. The one-dimensional mass is embedded in the multidimensional mass formulation and we obtain the discretization of the multidimensional equation.The algorithm is given in the following manner:(4.4)∂tcl+∇·vlcl=-λlcl+λl-1cl-1,withl=1,…,m.
The velocity vector v is divided by Rl. The initial conditions are given by c10=c1(x,0), else cl0=0forl=2,…,m and the boundary conditions are trivial cl=0 for l=1,…,m.We first calculate the maximum time step for cellj and concentration i with the use of the total outward fluxes:(4.5)τi,j=VjRiνj,νj=∑k∈out(j)vjk.
We get the restricted time step with the local time steps of cells and their components(4.6)τn≤mini=1,…,mj=1,…,Iτi,j.
The velocity of the discrete equation is given by(4.7)vi,j=1τi,j.
We calculate the analytical solution of the mass, compare for example,  [18], and we get(4.8)mi,jk,outn=mi,out(a,b,τn,v1,j,…,vi,j,R1,…,Ri,λ1,…,λi),mi,j,restn=mi,jnf(τn,v1,j,…,vi,j,R1,…,Ri,λ1,…,λi),
where a=VjRi(ci,jkn-ci,jk'n), b=VjRici,jk′n, and mi,jn=VjRici,jn. Furthermore, ci,jk′n is the concentration at the inflow and ci,jkn is the concentration at the outflow boundary of cell j.The discretization with the embedded analytical mass is calculated by:(4.9)mi,jn+1-mi,restn=-∑k∈out(j)vjkνjmi,jk,out+∑l∈in(j)vljνlmi,lj,out,
where vjk/νj is the retransformation of the total mass mi,jk,out into the partial mass mi,jk. In the next time step, the mass is given as mi,jn+1=Vjci,jn+1 and in the old time step it is the rest mass for concentration i. The proof is given in  [18]. In the next section, we derive an analytical solution for the benchmark problem, compare for example,  [21, 22].In the next subsection, we introduce the discretization of the diffusion-dispersion-equation.


#### 4.1.3. Discretization of the Diffusion-Dispersion-Equation

We discretize the diffusion-dispersion-equation with implicit time-discretization and finite-volume method for the following equation:(4.10)∂tRc-∇·(D∇c)=0,
where c=c(x,t) with x∈Ω and t≥0. The diffusion-dispersion tensor D=D(x,v) is given by the Scheidegger approach, compare for example,  (4.15). The velocity is given as v. The retardation factor is R>0.0.The boundary-values are denoted byn·D∇c(x,t)=0, where x∈Γ is the boundary Γ=∂Ω, compare for example,  [23]. The initial conditions are given by c(x,0)=c0(x).We integrate (4.10) over space and time and derive(4.11)∫Ωj∫tntn+1∂tR(c)dtdx=∫Ωj∫tntn+1∇·(D∇c)dtdx.
The time-integration is done by the backward-Euler method and the diffusion-dispersion term is lumped, compare for example,  [18].(4.12)∫Ωj(R(cn+1)-R(cn))dx=τn∫Ωj∇·(D∇cn+1)dx.
Equation (4.12) is discretized over space by using Green's formula.(4.13)∫Ωj(R(cn+1)-R(cn))dx=τn∫ΓjDn·∇cn+1dγ,
where Γj is the boundary of the finite-volume cell Ωj. We use an approximation in space, compare for example,  [18].The spatial integration of (4.13) is done by the mid-point rule over finite boundaries and is given by(4.14)VjR(cjn+1)-VjR(cjn)=τn∑e∈Λj∑k∈Λje|Γjke|njke·Djke∇cjke,n+1,
where |Γjke| is the length of the boundary-element Γjke. The gradients are calculated with the piecewise finite-element function ϕl, see [24], and we obtain(4.15)∇cjke,n+1=∑l∈Λecln+1∇ϕl(xjke).
We use the difference notation for the neighbour points j and l, compare for example,  [25], and obtain the discretized equation:(4.16)VjR(cjn+1)-VjR(cjn)=τn∑e∈Λj‍∑l∈Λe∖{j}(∑k∈Λje|Γjke|njke·Djke∇ϕl(xjke))(cjn+1-cln+1),
where j=1,…,m.


### 4.2. Interpolation and Regression of Experimental Dates

To simulate the physical experiments with the assumed model, we have to approximate the parameters of the numerical model. We apply interpolation and regression schemes to approximate between the mathematical and physical parameters.Here, we concentrate on the reaction rates of speciesSi, C, and H.The physical data of temperature and pressure are used and validation simulations are done to obtain the rate of deposition.Next we have to interpolate the parameters of the numerical model.(1) Lagrangian Interpolation
We assume an interpolation atΩ=[a1,b1]×⋯×[ad,bd].
(4.17)T=∑ν∈Kf(xνt)Lνt,
where the Lagrangian function is given as
(4.18)Lνt(x)=πi=1dπμ=0,μ≠νimxi-xμ[ai,bi]xνi[ai,bi]-xμ[ai,bi].(2) Linear Regression (Least squares approximation)
Here, we have points with values and we assume we have the best approximation by minimization:(4.19)S=∑k=1m(yk-Ln(xk))2,
where m≥n and Ln is a function that is constructed by the least squares algorithm; see [26].Remark 4.1.
To apply larger parameter spaces, we can generalise to multivariate regression methods see [27]. Here we compute approximations between higher dimensional matrix spaces.


## 4.1. Discretization and Solver Methods

For space-discretization of the PDEs, we apply finite-volume methods in mass conserved discretization schemes and for time-discretization of the resulting ODEs we apply Runge-Kutta methods or BDF methods. To accelerate the solver process, we combine the numerical and analytical parts of the solutions.


### 4.1.1. Discretization Method for Convection Equation

We deal with the following convection equation(4.1)∂tRc-v·∇c=0,
where R is the retardation factor and presents the retention of the concentration; see also (2.2). v is the velocity. We have a simple boundary condition c=0 for the inflow and outflow boundary and the initial values are given as c(xj,0)=cj0(x). We use a piecewise constant discretization method with the upwind discretization done in  [19] and get(4.2)VjRcjn+1=VjRcjn-τn∑k∈out(j)vjkcjn+τn∑k∈in(j)cknvkj,VjRcjn+1=cjn(RVj-τnνj)+τn∑k∈in(j)cknvkj.
The explicit time discretization has to fulfill the discrete minimum-maximum property  [19] and we get the following restriction on time steps:(4.3)τj=RVjνj,τn≤minj=1,…,Iτj.
To obtain improved spatial discretization methods and apply larger time steps, we introduce a reconstruction with linear polynomials as a higher test-function in the next subsection.


### 4.1.2. Discretization Method for Convection-Reaction Equation Based on Embedded One-Dimensional Analytical Solutions

We apply Godunov's discretization method, compare for example,  [20], and extend the formulation to the analytical solution of convection-reaction equations. We reduce the multidimensional equation to one-dimensional equations and solve each equation exactly. The one-dimensional solution is multiplied by the underlying volume to give the mass formulation. The one-dimensional mass is embedded in the multidimensional mass formulation and we obtain the discretization of the multidimensional equation.The algorithm is given in the following manner:(4.4)∂tcl+∇·vlcl=-λlcl+λl-1cl-1,withl=1,…,m.
The velocity vector v is divided by Rl. The initial conditions are given by c10=c1(x,0), else cl0=0forl=2,…,m and the boundary conditions are trivial cl=0 for l=1,…,m.We first calculate the maximum time step for cellj and concentration i with the use of the total outward fluxes:(4.5)τi,j=VjRiνj,νj=∑k∈out(j)vjk.
We get the restricted time step with the local time steps of cells and their components(4.6)τn≤mini=1,…,mj=1,…,Iτi,j.
The velocity of the discrete equation is given by(4.7)vi,j=1τi,j.
We calculate the analytical solution of the mass, compare for example,  [18], and we get(4.8)mi,jk,outn=mi,out(a,b,τn,v1,j,…,vi,j,R1,…,Ri,λ1,…,λi),mi,j,restn=mi,jnf(τn,v1,j,…,vi,j,R1,…,Ri,λ1,…,λi),
where a=VjRi(ci,jkn-ci,jk'n), b=VjRici,jk′n, and mi,jn=VjRici,jn. Furthermore, ci,jk′n is the concentration at the inflow and ci,jkn is the concentration at the outflow boundary of cell j.The discretization with the embedded analytical mass is calculated by:(4.9)mi,jn+1-mi,restn=-∑k∈out(j)vjkνjmi,jk,out+∑l∈in(j)vljνlmi,lj,out,
where vjk/νj is the retransformation of the total mass mi,jk,out into the partial mass mi,jk. In the next time step, the mass is given as mi,jn+1=Vjci,jn+1 and in the old time step it is the rest mass for concentration i. The proof is given in  [18]. In the next section, we derive an analytical solution for the benchmark problem, compare for example,  [21, 22].In the next subsection, we introduce the discretization of the diffusion-dispersion-equation.


### 4.1.3. Discretization of the Diffusion-Dispersion-Equation

We discretize the diffusion-dispersion-equation with implicit time-discretization and finite-volume method for the following equation:(4.10)∂tRc-∇·(D∇c)=0,
where c=c(x,t) with x∈Ω and t≥0. The diffusion-dispersion tensor D=D(x,v) is given by the Scheidegger approach, compare for example,  (4.15). The velocity is given as v. The retardation factor is R>0.0.The boundary-values are denoted byn·D∇c(x,t)=0, where x∈Γ is the boundary Γ=∂Ω, compare for example,  [23]. The initial conditions are given by c(x,0)=c0(x).We integrate (4.10) over space and time and derive(4.11)∫Ωj∫tntn+1∂tR(c)dtdx=∫Ωj∫tntn+1∇·(D∇c)dtdx.
The time-integration is done by the backward-Euler method and the diffusion-dispersion term is lumped, compare for example,  [18].(4.12)∫Ωj(R(cn+1)-R(cn))dx=τn∫Ωj∇·(D∇cn+1)dx.
Equation (4.12) is discretized over space by using Green's formula.(4.13)∫Ωj(R(cn+1)-R(cn))dx=τn∫ΓjDn·∇cn+1dγ,
where Γj is the boundary of the finite-volume cell Ωj. We use an approximation in space, compare for example,  [18].The spatial integration of (4.13) is done by the mid-point rule over finite boundaries and is given by(4.14)VjR(cjn+1)-VjR(cjn)=τn∑e∈Λj∑k∈Λje|Γjke|njke·Djke∇cjke,n+1,
where |Γjke| is the length of the boundary-element Γjke. The gradients are calculated with the piecewise finite-element function ϕl, see [24], and we obtain(4.15)∇cjke,n+1=∑l∈Λecln+1∇ϕl(xjke).
We use the difference notation for the neighbour points j and l, compare for example,  [25], and obtain the discretized equation:(4.16)VjR(cjn+1)-VjR(cjn)=τn∑e∈Λj‍∑l∈Λe∖{j}(∑k∈Λje|Γjke|njke·Djke∇ϕl(xjke))(cjn+1-cln+1),
where j=1,…,m.


## 4.1.1. Discretization Method for Convection Equation

We deal with the following convection equation(4.1)∂tRc-v·∇c=0,
where R is the retardation factor and presents the retention of the concentration; see also (2.2). v is the velocity. We have a simple boundary condition c=0 for the inflow and outflow boundary and the initial values are given as c(xj,0)=cj0(x). We use a piecewise constant discretization method with the upwind discretization done in  [19] and get(4.2)VjRcjn+1=VjRcjn-τn∑k∈out(j)vjkcjn+τn∑k∈in(j)cknvkj,VjRcjn+1=cjn(RVj-τnνj)+τn∑k∈in(j)cknvkj.
The explicit time discretization has to fulfill the discrete minimum-maximum property  [19] and we get the following restriction on time steps:(4.3)τj=RVjνj,τn≤minj=1,…,Iτj.
To obtain improved spatial discretization methods and apply larger time steps, we introduce a reconstruction with linear polynomials as a higher test-function in the next subsection.


## 4.1.2. Discretization Method for Convection-Reaction Equation Based on Embedded One-Dimensional Analytical Solutions

We apply Godunov's discretization method, compare for example,  [20], and extend the formulation to the analytical solution of convection-reaction equations. We reduce the multidimensional equation to one-dimensional equations and solve each equation exactly. The one-dimensional solution is multiplied by the underlying volume to give the mass formulation. The one-dimensional mass is embedded in the multidimensional mass formulation and we obtain the discretization of the multidimensional equation.The algorithm is given in the following manner:(4.4)∂tcl+∇·vlcl=-λlcl+λl-1cl-1,withl=1,…,m.
The velocity vector v is divided by Rl. The initial conditions are given by c10=c1(x,0), else cl0=0forl=2,…,m and the boundary conditions are trivial cl=0 for l=1,…,m.We first calculate the maximum time step for cellj and concentration i with the use of the total outward fluxes:(4.5)τi,j=VjRiνj,νj=∑k∈out(j)vjk.
We get the restricted time step with the local time steps of cells and their components(4.6)τn≤mini=1,…,mj=1,…,Iτi,j.
The velocity of the discrete equation is given by(4.7)vi,j=1τi,j.
We calculate the analytical solution of the mass, compare for example,  [18], and we get(4.8)mi,jk,outn=mi,out(a,b,τn,v1,j,…,vi,j,R1,…,Ri,λ1,…,λi),mi,j,restn=mi,jnf(τn,v1,j,…,vi,j,R1,…,Ri,λ1,…,λi),
where a=VjRi(ci,jkn-ci,jk'n), b=VjRici,jk′n, and mi,jn=VjRici,jn. Furthermore, ci,jk′n is the concentration at the inflow and ci,jkn is the concentration at the outflow boundary of cell j.The discretization with the embedded analytical mass is calculated by:(4.9)mi,jn+1-mi,restn=-∑k∈out(j)vjkνjmi,jk,out+∑l∈in(j)vljνlmi,lj,out,
where vjk/νj is the retransformation of the total mass mi,jk,out into the partial mass mi,jk. In the next time step, the mass is given as mi,jn+1=Vjci,jn+1 and in the old time step it is the rest mass for concentration i. The proof is given in  [18]. In the next section, we derive an analytical solution for the benchmark problem, compare for example,  [21, 22].In the next subsection, we introduce the discretization of the diffusion-dispersion-equation.


## 4.1.3. Discretization of the Diffusion-Dispersion-Equation

We discretize the diffusion-dispersion-equation with implicit time-discretization and finite-volume method for the following equation:(4.10)∂tRc-∇·(D∇c)=0,
where c=c(x,t) with x∈Ω and t≥0. The diffusion-dispersion tensor D=D(x,v) is given by the Scheidegger approach, compare for example,  (4.15). The velocity is given as v. The retardation factor is R>0.0.The boundary-values are denoted byn·D∇c(x,t)=0, where x∈Γ is the boundary Γ=∂Ω, compare for example,  [23]. The initial conditions are given by c(x,0)=c0(x).We integrate (4.10) over space and time and derive(4.11)∫Ωj∫tntn+1∂tR(c)dtdx=∫Ωj∫tntn+1∇·(D∇c)dtdx.
The time-integration is done by the backward-Euler method and the diffusion-dispersion term is lumped, compare for example,  [18].(4.12)∫Ωj(R(cn+1)-R(cn))dx=τn∫Ωj∇·(D∇cn+1)dx.
Equation (4.12) is discretized over space by using Green's formula.(4.13)∫Ωj(R(cn+1)-R(cn))dx=τn∫ΓjDn·∇cn+1dγ,
where Γj is the boundary of the finite-volume cell Ωj. We use an approximation in space, compare for example,  [18].The spatial integration of (4.13) is done by the mid-point rule over finite boundaries and is given by(4.14)VjR(cjn+1)-VjR(cjn)=τn∑e∈Λj∑k∈Λje|Γjke|njke·Djke∇cjke,n+1,
where |Γjke| is the length of the boundary-element Γjke. The gradients are calculated with the piecewise finite-element function ϕl, see [24], and we obtain(4.15)∇cjke,n+1=∑l∈Λecln+1∇ϕl(xjke).
We use the difference notation for the neighbour points j and l, compare for example,  [25], and obtain the discretized equation:(4.16)VjR(cjn+1)-VjR(cjn)=τn∑e∈Λj‍∑l∈Λe∖{j}(∑k∈Λje|Γjke|njke·Djke∇ϕl(xjke))(cjn+1-cln+1),
where j=1,…,m.


## 4.2. Interpolation and Regression of Experimental Dates

To simulate the physical experiments with the assumed model, we have to approximate the parameters of the numerical model. We apply interpolation and regression schemes to approximate between the mathematical and physical parameters.Here, we concentrate on the reaction rates of speciesSi, C, and H.The physical data of temperature and pressure are used and validation simulations are done to obtain the rate of deposition.Next we have to interpolate the parameters of the numerical model.(1) Lagrangian Interpolation
We assume an interpolation atΩ=[a1,b1]×⋯×[ad,bd].
(4.17)T=∑ν∈Kf(xνt)Lνt,
where the Lagrangian function is given as
(4.18)Lνt(x)=πi=1dπμ=0,μ≠νimxi-xμ[ai,bi]xνi[ai,bi]-xμ[ai,bi].(2) Linear Regression (Least squares approximation)
Here, we have points with values and we assume we have the best approximation by minimization:(4.19)S=∑k=1m(yk-Ln(xk))2,
where m≥n and Ln is a function that is constructed by the least squares algorithm; see [26].Remark 4.1.
To apply larger parameter spaces, we can generalise to multivariate regression methods see [27]. Here we compute approximations between higher dimensional matrix spaces.


## 5. Numerical Experiments

For all the experiments, we have the following parameters of the model, discretization and solver methods.We apply interpolation and regression methods to couple the physical parameters to the mathematical parameters; see Figure4 and Table 4.Table 4
Physical and mathematical parameters.
Physical parameterMathematical parameterTemperature, pressure, powerVelocity, diffusion, reactionT,p,WV,D,λFigure 4
Coupling of physical and mathematical parameter space.Parameters of the Equation
In the following, we list the parameters for our simulation tool UG; see [28]. The software toolbox has a flexible user interface to allow a large number of numerical experiments and approximations to the known physical parameters. The model parameters are given in Table 5.Table 5
Model-Parameters.
Densityρ=1.0Mobile porosityϕ=0.333Immobile porosity0.333DiffusionD=0.0Longitudinal dispersionαL=0.0Transversal dispersionαT=0.00Retardation factorR=10.0e-4 (Henry rate).Velocity fieldv=(0.0,-4.0·10-8)t.Decay rate of species of 1st EXλAB=1·10-68.Decay rate of species of 2nd EXλAB=2·10-8, λBNN=1·10-68.Decay rate of species of 3rd EXλAB=0.25·10-8, λCB=0.5·10-8.Geometry (2d domain)Ω=[0,100]×[0,100].BoundaryNeumann boundary at top, left and right boundaries.Outflow boundary at bottom boundaryDiscretization Method
Finite-volume method of 2nd order are given with theparameters in Table6.Table 6
Time discretization parameters.
Spatial step sizeΔxmin=1.56,Δxmax=2.21Refined levels6LimiterSlope limiterTest functionsLinear test function reconstructed with neighbor gradientsTime Discretization Method
Crank-Nicolson method of 2nd order are given with theparameters in Table7.Table 7
Spatial discretization parameters.
Initial time stepΔtinit=5·107Controlled time stepΔtmax=1.298·107,Δtmin=1.158·107Number of time steps100,80,30,25Time step controlTime steps are controlled with the Courant-NumberCFLmax=1Solver Method
In the following, we deal with test examples which are approximations of physical experiments. The parameters of the solvers are given in Table8.Table 8
Solver methods and their parameters.
SolverBiCGstab (Bi conjugate gradient method)PreconditionerGeometric Multi-grid methodSmootherGauss-Seidel method as smoothers for the Multi-grid methodBasic level0Initial gridUniform grid with2 elementsMaximum Level6Finest gridUniform grid with8192 elements


### 5.1. Test Experiment 1: Interpolation with Temperature

In the test example, we deal with the following reaction:(5.1)2SiC+4H→λSiC+CH4+Si.
Here, we have physical experiments and approximate temperature parameters T=400,600,800.We compute theSiC:C ratio at a given temperature T=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. These values can be given back to physical experiments; see Table 9.Table 9
Computed and experimental fitted parameters from UG simulations.
Tλ (fitted)λ (interpolated)Ratio(SiC:C) (computed with UG)4001/2·10-82.4: 15000.35·10-81.85: 16001/4·10-81.5: 17000.171·10-81.211: 18001/8·10-81.1: 1One Source
The parameters of the source concentration are givenin Table10. In Figure 5, we present the concentration of one point source at (50,20) with number of time steps being 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=1·108Amount of permanent source concentrationcsource=1.0Number of time steps25One point source at (50,20) with number of time steps being equal to 25.
(a)(b)In Figure6, we show the deposition rates of one point source at (50,20), with number of time steps being equal to 25. The rate of concentrations is given in Table 12.Figure 6
Deposition rates in case of one point source at (50,20) with number of time steps being equal to 25.Nine Point Sources
In this experiment, we apply nine point sources.
The parameters of the source concentration are givenin Table11. In Figure 7, we present the concentration with nine point sources in a short time.
In Figure8, we show the deposition rate with the nine point sources, with number of time steps being equal to 25. The rate of concentrations is given in Table 13.Table 11
Parameters of source concentration.
9 point sources at positionX=10,…,90,Y=20Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationcsource=1.0Number of time steps25Table 12
Rate of concentration.
RATESiCsource,max:SiCtarget,max9·106:6.5·106=1.38Table 13
Rate of concentration.
RATESiCsource,max:SiCtarget,max9·106:6.7·106=1.34Nine point sources with number of time steps being equal to 25.
(a)(b)Figure 8
Deposition rates in case of nine point sources, with number of time steps being equal to 25.81-Point Sources
In this experiment, we apply 81 point sources.
The parameters of the source concentration are givenin Table14. In Figure 9, we present the concentration with 81 point sources.
In Figure10, we show the deposition rates with 81 point sources, with number of time steps equal to 80. The rate of concentrations is given in Table 15.Table 14
Parameters of source concentration.
81 point sources at positionX=10,11,12,…,90,Y=20Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationcsource=1.0Number of time steps80Table 15
Rate of concentration.
RATESiCsource,max:SiCtarget,max1.5·107:1.5·107=181 point sources with number of time steps equal to 80.
(a)(b)Figure 10
Deposition rates in case of 81 point sources, with number of time steps equal to 80.Line Source
In this part, we will perform experiments with a line source. The parameters of the line source are given in Table16.
In Figure11, we present the results with a line source, x∈[5,95],y∈[20,25] with number of time steps equal to 25.
In Figure12, we show the deposition rates with (a) the line source, x is from 5 to 95, and y is from 20 to 25. The rate of concentrations is given in Table 17.Table 16
Parameters the source concentration.
Line source at positionx∈[5,95],y∈[20,25]Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationcsource=1.0Number of time steps25Table 17
Rate of concentration.
RATESiCsource,max:SiCtarget,max4.7·107:4·107=1.17Line source with number of time steps being equal to 25.
(a)(b)Figure 12
Deposition rates in case of line source,x∈[5,95],y∈[20,25].


### 5.2. Test Experiment 2: Interpolation with Temperature and Power

In the next experiment, we fit the mathematical parameters to the temperature and power of the physical experiments.We deal with the reaction:(5.2)2SiC+4H→λSiC+CH4+Si.
In this case, we have a table which has values for the temperature and power of the plasma and for the ratio between the sources.We have to interpolateλ to the physical parameters temperature T and power of plasma P. In Table 18, the interpolated parameters are given.Table 18
Computed (C) and experimentally fitted (F) parameters from UG simulations.
TPλRatio (SiC:C)Computed Ratio400900S.1/10·10-8F.1: 0.971.01400500S.1/5·10-8F.1.3: 11.334001001/2·10-8C.2.4: 16003001/4·10-8C.1.5: 18005001/8·10-8C.1: 1800900S.1/5.7·10-8F.1.18: 11.252One Source
The parameters of the source concentration are givenin Table19. In Figure 13, we present the concentration of one point source at (50,20) with number of time steps being equal to 25.   In Figure 14, we show the deposition rates with one point source at (50,20), with number of time steps being equal to 25. The rate of concentrations is given in Table 20.Table 19
Parameters of source concentration.
Point source at position(x,y)=(50,20)Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationcsource=1.0Number of time steps25Table 20
Rate of concentration.
RATESiCsource,max:Ctarget,max3·106:3·106=1One point source at (50,20) with number of timesteps being equal to 25.
(a)(b)(c)(d)Figure 14
Deposition rates in case of one point source at (50,20), with number of time steps being equal to 25.Nine Point Sources
In this experiment, we apply nine point sources.
The parameters of the source concentration are given in Table21. In Figure 15, we present the concentration of the nine point sources in a short time.
In Figure16, we show the deposition rates of nine point sources, with number of time steps being equal to 25.Table 21
Parameters of source concentration.
Nine point sources at position(x=10,20,30,40,50,60,70,80,90;y=20)Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationcsource=1.0Number of time steps25Nine point sources with number of time steps being equal to 25.
(a)(b)(c)(d)Figure 16
Deposition rates in case of nine point sources, with number of time steps being equal to 25.81-Point Sources
In this experiment, we apply 81 point sources.
The parameters of the source concentration are givenin Table23. In Figure 17, we present the concentration with 81 point sources.
In Figure18, we show the deposition rates with 81 point sources, with number of time steps equal to 100. The rate of concentrations is given in Table 22.Table 22
Rate of concentration.
RATESiCsource,max:Ctarget,max3.106:3.106=1Table 23
Parameters of source concentration.
81 point sources at positionX=10,11,12,000,90;Y=20.Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationcsource=1.0Number of time steps10081 point sources with number of time steps equal to 100.
(a)(b)(c)(d)Figure 18
Deposition rates in case of 81 point sources, with number of time steps equal to 100.Line Source
In this part, we will perform an experiments with a line source,x∈[5,95],y∈[20,25].
In Figure19, we present the result of the line source, x∈[5,95],y∈[20,25], with number of time being steps equal to 30.
In Figure20, we show the deposition rates with a line source, with x from 5 to 95, and y from 20 to 30.Line source with number of time steps being equal to 25.
(a)(b)(c)(d)Figure 20
Deposition rates in case of line source,x∈[5,95],y∈[20,25].


### 5.3. Test Experiment 3: Regression with Temperature and Power

In the next experiment, we apply a more flexible approximation method to obtain the parameters of the mathematical method. We apply a regression and fit to all the physical parameters, because we are not restricted to a given interpolation grid.The reaction is given asA→B and B→C and we apply to 2SiC+4H→SiC+CH4+Si.We computed the ratioSiC:C for temperatures T=400,600,800 and plasma powers 100,300,500,900 and fit the simulated ratio given by the UG program to the mathematical model with a reaction parameter λ.We use linear regression, see Section4, and compute λ for the new temperatures T=450,500,800 and apply the ratio of the simulated new parameters. These values can be given back to the physical experiments; see Table 24.Table 24
Rate of concentration.
RATESiCsource,max:Ctarget,max7.5·106:7·106=1.07Table 25
Parameters of source concentration.
Line source at positionx∈[5,95],y∈[20,25]Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationcsource=1.0Number of time steps25Table 26
Rate of concentration.
RATESiCsource,max:Ctarget,max1.8·107:2.2·107=0.81Table 27
Parameters of source concentration.
TPExactλRegressionλExact ratio (SiC:C)Regression ratio (SiC:C)4009001e-091.703e-091:0.970.8354005000.2e-082.903e-091.3:11.6164001000.5e-084.103e-092.4:12.0116003000.25e-083.303e-091.5:11.7748005000.125e-82.503e-091:11.1928009000.175e-81.303e-091.2:11.1325005002.803e-091.586006002.4030e-091.4338008001.603e-091.2064004003.203e-091.7154504502.703e-091.578001003.703e-091.93Table 28
Parameters of source concentration.
Point source ofSiC at position(x,y)=(50,20)Point source ofH at position(x,y)=(50,20)Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationcsource=1.0Number of time steps100Table 29
Rate of concentration.
RATECsource,max:SiCtarget,max1.8·107:1·107=1.8Table 30
Parameters of source concentration.
Nine point sources ofSiC at position(x=10,20,30,40,50,60,70,80,90;y=20)Nine point sources ofSiC at position(x=10,20,30,40,50,60,70,80,90;y=20)Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationcsource=1.0Number of time steps25Table 31
Rate of concentration.
RATECsource,max:SiCtarget,max5·106:4.4·106=1.13Table 32
Parameters of source concentration.
81 point sources ofSiC at positionX=10,11,12,000,90;Y=20Line source ofH at positionx∈[5,95],y∈[20,25]Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationSicsource=1.0,Hsource=0.20Number of time steps100Table 33
Rate of concentration.
RATECsource,max:SiCtarget,max1.8·107:0.75·107=2.4One Source
Here, we take points sources.
In Figure21, we present the concentration in the one-point source experiment.
In Figure22, we show the deposition rates in the one-point source experiment.One point source experiment.
(a)(b)(c)(d)(e)(f)Figure 22
Deposition rates in case of one-point source experiment.Nine Point Sources
Here, we take nine point sources of both concentrations.
In Figure23, we present the concentration of the nine point sources experiment.
In Figure24, we show the deposition rates of the nine point sources experiment.Nine-point source experiment.
(a)(b)(c)Figure 24
Deposition rates in case of nine-point source experiment.81-Point Sources, thalf=2·108
In this first experiment, the value of temperature is 400 deg C andλ is 0.5·10-8.   Here, we take the concentration of SiC as a point source, and the concentration of H as a line source.
In Figure25, we present the concentration of the 81-point source experiment.
In Figure26, we show the deposition rates in the 81 point source experiment.81-point source experiment.
(a)(b)(c)(d)(e)(f)Figure 26
Deposition rates in case of 81 point source experiment.Remark 5.1.
The regression method is more flexible for approximating to the physical parameters. We obtain numerical results for different parameter studies, that are fitted to the physical experiments. The first test examples with multiple sources and temperature regions which are interesting to physicists are simulated. Here, we have coupled a mathematical model to a physical experiment and studied the near region of the deposition process.


## 5.1. Test Experiment 1: Interpolation with Temperature

In the test example, we deal with the following reaction:(5.1)2SiC+4H→λSiC+CH4+Si.
Here, we have physical experiments and approximate temperature parameters T=400,600,800.We compute theSiC:C ratio at a given temperature T=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. These values can be given back to physical experiments; see Table 9.Table 9
Computed and experimental fitted parameters from UG simulations.
Tλ (fitted)λ (interpolated)Ratio(SiC:C) (computed with UG)4001/2·10-82.4: 15000.35·10-81.85: 16001/4·10-81.5: 17000.171·10-81.211: 18001/8·10-81.1: 1One Source
The parameters of the source concentration are givenin Table10. In Figure 5, we present the concentration of one point source at (50,20) with number of time steps being 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=1·108Amount of permanent source concentrationcsource=1.0Number of time steps25One point source at (50,20) with number of time steps being equal to 25.
(a)(b)In Figure6, we show the deposition rates of one point source at (50,20), with number of time steps being equal to 25. The rate of concentrations is given in Table 12.Figure 6
Deposition rates in case of one point source at (50,20) with number of time steps being equal to 25.Nine Point Sources
In this experiment, we apply nine point sources.
The parameters of the source concentration are givenin Table11. In Figure 7, we present the concentration with nine point sources in a short time.
In Figure8, we show the deposition rate with the nine point sources, with number of time steps being equal to 25. The rate of concentrations is given in Table 13.Table 11
Parameters of source concentration.
9 point sources at positionX=10,…,90,Y=20Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationcsource=1.0Number of time steps25Table 12
Rate of concentration.
RATESiCsource,max:SiCtarget,max9·106:6.5·106=1.38Table 13
Rate of concentration.
RATESiCsource,max:SiCtarget,max9·106:6.7·106=1.34Nine point sources with number of time steps being equal to 25.
(a)(b)Figure 8
Deposition rates in case of nine point sources, with number of time steps being equal to 25.81-Point Sources
In this experiment, we apply 81 point sources.
The parameters of the source concentration are givenin Table14. In Figure 9, we present the concentration with 81 point sources.
In Figure10, we show the deposition rates with 81 point sources, with number of time steps equal to 80. The rate of concentrations is given in Table 15.Table 14
Parameters of source concentration.
81 point sources at positionX=10,11,12,…,90,Y=20Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationcsource=1.0Number of time steps80Table 15
Rate of concentration.
RATESiCsource,max:SiCtarget,max1.5·107:1.5·107=181 point sources with number of time steps equal to 80.
(a)(b)Figure 10
Deposition rates in case of 81 point sources, with number of time steps equal to 80.Line Source
In this part, we will perform experiments with a line source. The parameters of the line source are given in Table16.
In Figure11, we present the results with a line source, x∈[5,95],y∈[20,25] with number of time steps equal to 25.
In Figure12, we show the deposition rates with (a) the line source, x is from 5 to 95, and y is from 20 to 25. The rate of concentrations is given in Table 17.Table 16
Parameters the source concentration.
Line source at positionx∈[5,95],y∈[20,25]Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationcsource=1.0Number of time steps25Table 17
Rate of concentration.
RATESiCsource,max:SiCtarget,max4.7·107:4·107=1.17Line source with number of time steps being equal to 25.
(a)(b)Figure 12
Deposition rates in case of line source,x∈[5,95],y∈[20,25].


## 5.2. Test Experiment 2: Interpolation with Temperature and Power

In the next experiment, we fit the mathematical parameters to the temperature and power of the physical experiments.We deal with the reaction:(5.2)2SiC+4H→λSiC+CH4+Si.
In this case, we have a table which has values for the temperature and power of the plasma and for the ratio between the sources.We have to interpolateλ to the physical parameters temperature T and power of plasma P. In Table 18, the interpolated parameters are given.Table 18
Computed (C) and experimentally fitted (F) parameters from UG simulations.
TPλRatio (SiC:C)Computed Ratio400900S.1/10·10-8F.1: 0.971.01400500S.1/5·10-8F.1.3: 11.334001001/2·10-8C.2.4: 16003001/4·10-8C.1.5: 18005001/8·10-8C.1: 1800900S.1/5.7·10-8F.1.18: 11.252One Source
The parameters of the source concentration are givenin Table19. In Figure 13, we present the concentration of one point source at (50,20) with number of time steps being equal to 25.   In Figure 14, we show the deposition rates with one point source at (50,20), with number of time steps being equal to 25. The rate of concentrations is given in Table 20.Table 19
Parameters of source concentration.
Point source at position(x,y)=(50,20)Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationcsource=1.0Number of time steps25Table 20
Rate of concentration.
RATESiCsource,max:Ctarget,max3·106:3·106=1One point source at (50,20) with number of timesteps being equal to 25.
(a)(b)(c)(d)Figure 14
Deposition rates in case of one point source at (50,20), with number of time steps being equal to 25.Nine Point Sources
In this experiment, we apply nine point sources.
The parameters of the source concentration are given in Table21. In Figure 15, we present the concentration of the nine point sources in a short time.
In Figure16, we show the deposition rates of nine point sources, with number of time steps being equal to 25.Table 21
Parameters of source concentration.
Nine point sources at position(x=10,20,30,40,50,60,70,80,90;y=20)Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationcsource=1.0Number of time steps25Nine point sources with number of time steps being equal to 25.
(a)(b)(c)(d)Figure 16
Deposition rates in case of nine point sources, with number of time steps being equal to 25.81-Point Sources
In this experiment, we apply 81 point sources.
The parameters of the source concentration are givenin Table23. In Figure 17, we present the concentration with 81 point sources.
In Figure18, we show the deposition rates with 81 point sources, with number of time steps equal to 100. The rate of concentrations is given in Table 22.Table 22
Rate of concentration.
RATESiCsource,max:Ctarget,max3.106:3.106=1Table 23
Parameters of source concentration.
81 point sources at positionX=10,11,12,000,90;Y=20.Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationcsource=1.0Number of time steps10081 point sources with number of time steps equal to 100.
(a)(b)(c)(d)Figure 18
Deposition rates in case of 81 point sources, with number of time steps equal to 100.Line Source
In this part, we will perform an experiments with a line source,x∈[5,95],y∈[20,25].
In Figure19, we present the result of the line source, x∈[5,95],y∈[20,25], with number of time being steps equal to 30.
In Figure20, we show the deposition rates with a line source, with x from 5 to 95, and y from 20 to 30.Line source with number of time steps being equal to 25.
(a)(b)(c)(d)Figure 20
Deposition rates in case of line source,x∈[5,95],y∈[20,25].


## 5.3. Test Experiment 3: Regression with Temperature and Power

In the next experiment, we apply a more flexible approximation method to obtain the parameters of the mathematical method. We apply a regression and fit to all the physical parameters, because we are not restricted to a given interpolation grid.The reaction is given asA→B and B→C and we apply to 2SiC+4H→SiC+CH4+Si.We computed the ratioSiC:C for temperatures T=400,600,800 and plasma powers 100,300,500,900 and fit the simulated ratio given by the UG program to the mathematical model with a reaction parameter λ.We use linear regression, see Section4, and compute λ for the new temperatures T=450,500,800 and apply the ratio of the simulated new parameters. These values can be given back to the physical experiments; see Table 24.Table 24
Rate of concentration.
RATESiCsource,max:Ctarget,max7.5·106:7·106=1.07Table 25
Parameters of source concentration.
Line source at positionx∈[5,95],y∈[20,25]Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationcsource=1.0Number of time steps25Table 26
Rate of concentration.
RATESiCsource,max:Ctarget,max1.8·107:2.2·107=0.81Table 27
Parameters of source concentration.
TPExactλRegressionλExact ratio (SiC:C)Regression ratio (SiC:C)4009001e-091.703e-091:0.970.8354005000.2e-082.903e-091.3:11.6164001000.5e-084.103e-092.4:12.0116003000.25e-083.303e-091.5:11.7748005000.125e-82.503e-091:11.1928009000.175e-81.303e-091.2:11.1325005002.803e-091.586006002.4030e-091.4338008001.603e-091.2064004003.203e-091.7154504502.703e-091.578001003.703e-091.93Table 28
Parameters of source concentration.
Point source ofSiC at position(x,y)=(50,20)Point source ofH at position(x,y)=(50,20)Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationcsource=1.0Number of time steps100Table 29
Rate of concentration.
RATECsource,max:SiCtarget,max1.8·107:1·107=1.8Table 30
Parameters of source concentration.
Nine point sources ofSiC at position(x=10,20,30,40,50,60,70,80,90;y=20)Nine point sources ofSiC at position(x=10,20,30,40,50,60,70,80,90;y=20)Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationcsource=1.0Number of time steps25Table 31
Rate of concentration.
RATECsource,max:SiCtarget,max5·106:4.4·106=1.13Table 32
Parameters of source concentration.
81 point sources ofSiC at positionX=10,11,12,000,90;Y=20Line source ofH at positionx∈[5,95],y∈[20,25]Starting point of source concentrationtstart=0.0End point of source concentrationtend=1·108Amount of permanent source concentrationSicsource=1.0,Hsource=0.20Number of time steps100Table 33
Rate of concentration.
RATECsource,max:SiCtarget,max1.8·107:0.75·107=2.4One Source
Here, we take points sources.
In Figure21, we present the concentration in the one-point source experiment.
In Figure22, we show the deposition rates in the one-point source experiment.One point source experiment.
(a)(b)(c)(d)(e)(f)Figure 22
Deposition rates in case of one-point source experiment.Nine Point Sources
Here, we take nine point sources of both concentrations.
In Figure23, we present the concentration of the nine point sources experiment.
In Figure24, we show the deposition rates of the nine point sources experiment.Nine-point source experiment.
(a)(b)(c)Figure 24
Deposition rates in case of nine-point source experiment.81-Point Sources, thalf=2·108
In this first experiment, the value of temperature is 400 deg C andλ is 0.5·10-8.   Here, we take the concentration of SiC as a point source, and the concentration of H as a line source.
In Figure25, we present the concentration of the 81-point source experiment.
In Figure26, we show the deposition rates in the 81 point source experiment.81-point source experiment.
(a)(b)(c)(d)(e)(f)Figure 26
Deposition rates in case of 81 point source experiment.Remark 5.1.
The regression method is more flexible for approximating to the physical parameters. We obtain numerical results for different parameter studies, that are fitted to the physical experiments. The first test examples with multiple sources and temperature regions which are interesting to physicists are simulated. Here, we have coupled a mathematical model to a physical experiment and studied the near region of the deposition process.


## 6. Conclusions

We present a numerical simulation of a CVD process for depositingSiC films. Based on the different scales of physical and mathematical experiments, we apply a parameter approximation to fit the physical experiment to the mathematical experiment. Numerical approximations to the experimental data include the new parameters of the mathematical model. Such experiments allow to reduce the number of physical experiments to an acceptable number and give engineers and experimentalists a mathematical tool for predicting complex physical processes.The first numerical results show predictions of physical experiments with a transport-reaction equation of the deposition process.The temperature of the target and power of the plasma are chosen in such manner that the simulation results can help find the optimal deposition conditions. Furthermore, multiple sources obtain the best results with homogeneous layer deposition.Such numerical simulations help to predict the deposition rates of the underlying film, for example,SiC. In future, we will analyze the validity of the models with more complicated precursor gases. Here the outstanding multivariate analysis will be important for approximating a large number of parameters.



---

*Source: 407561-2010-06-29.xml*