Adaptive Step-Size Control in Simulation of Diffusive CVD Processes

Abstract

We present control strategies of a diffusion process for chemical vapor deposition for metallic bipolar plates. In the models, we discuss the application of different models to simulate the plasma-transport of chemical reactants in the gas-chamber. The contribution are an optimal control problem based on a PID control to obtain a homogeneous layering. We have taken into account one- and two-dimensional problems that are given with constraints and control functions. A finite-element formulation with adaptive feedback control for time-step selection has been developed for the diffusion process. The optimization is presented with efficient algorithms. Numerical experiments are discussed with respect to the diffusion processes of the macroscopic model.

1. Introduction

We motivate our studying on simulating a low-temperature low-pressure plasma that can be found in chemical vapor deposition (CVD) processes. In the last years the research and optimization in producing high-temperature films by depositing low-pressure processes have increased by using simulation tools, see [1, 2]. Theoretical models exist for deposition processes and can be modeled by coupled transport and flow equations, see [3, 4]. Further interest on standard applications to deposit titanium-nitrogen (TiN) and titanium-carbon (TiC) on metallic layers are immense, see [5]. Recently more and more focus on deposition with new material classes known as MAX-phases are becoming important, see [6, 7]. The MAX-phase are nanolayered terniar metal-carbides or -nitrids, where 𝑀 is a transition metal, 𝐴 is an A-group element (e.g., Al, Ga, In, Si, etc.) and 𝑋 is 𝐶 (carbon) or 𝑁 (nitride), see [8]. Such materials combine ceramic and metallic behavior and can be implanted in the metallic bipolar plates to obtain a new material with noncorrosive and good metallic conductivity behavior.

We discuss a model for low-temperature and low-pressure plasma that can be used to implant or deposit thin layers of important materials, see [9]. This model is used for the implantation process. First, we present the process in the plasma-reactor that transport the contaminants to the wafer surface, see [10]. We deal with a continuous flow model, while we assume a vacuum- and a diffusion-dominated processes. Second, the process at the wafer-surface is modeled by the heavy particles problem with underlying drift. This model deals more with the atomic behavior and we do not allow 𝑝=0, see [10].

To solve such optimization problems, we present a PID-controller (proportional, integral, differential) to control our deposition process, see [11]. We improved heuristic methods of deriving the PID parameters, while we discuss the posteriori error estimates respecting the time-step size control. Our contribution is a modified automatically step-size control and a best approximation is obtained with the time-dependent control method based on the Chien-Hrones-Reswick algorithm, see [12].

Numerical methods are described in the context of time- and spatial-discretization methods for the mesoscopic-scale model. We discussed different experiments and their convergence rates.

For the simulations we apply analytical and also numerical methods to obtain results to control the grow of thin layers.

The paper is outlined as follows. In Section 2 we present our mathematical model and a possible reduced model for the further approximations. In Section 3 we discuss the theoretical background for the simulation of CVD processes. The optimal control and their control paths based on the PID-control approach are discussed in Section 4. The software and program-tools are discussed in Section 5. The numerical experiments are given in Section 6. In the contents that are given in Section 7, we summarize our results.

2. Mathematical Model

In the following, the models are discussed in two directions of far-field and near-field problems:

(1)reaction-diffusion equations, see [13] (far-field problem);(2)Boltzmann-Lattice equations, see [9] (near-field problem). The modeling is considered by the Knudsen Number (Kn), which is the ratio of the mean free path 𝜆 over the typical domain size 𝐿. For small Knudsen numbers Kn≈0.01−1.0, we deal with a Navier-Stokes equation or with the convection-diffusion equation, see [5, 14], whereas for large Knudsen numbers Kn≥1.0, we deal with a Boltzmann equation, see [4].

2.1. Modeling with Partial Differential Equations

Dynamic processes with modifications in time and space will be reshaped by partial differential equations. There is (i) the PDE-formula itself which describes the physical laws of nature that influence the process and (ii) initial and boundary conditions in which specific characteristics of the process, like boundary behavior, can be coded.

There are two types of boundary conditions, namely, Dirichlet and Neumann boundary. With the Dirichlet type the exact value of the boundary is known, however, with Neumann boundaries the spatial derivation of the boundary values in normal direction is known, see an example of the boundary conditions in Figure 1.

(a)

(b)

Figure 1 

Dirichlet and Neumann boundary conditions.

(a)

(b)

Figure 2 

Spatial-discretization.

2.2. Model for Optimal Control of the Layer

We will concentrate us on a continuum model of mass transportation and assume that the energy and momentum is conserved, see [13]. Therefore, the continuum flow of the mass can be described as diffusion reaction equation given as𝜕_𝑡𝑐−∇𝐷𝑐(∇𝑥,𝑐0−)=𝑅_𝑔𝑐₀=0,inΩ×[0,𝑇],(2.1)𝜕𝑐(𝑥,𝑡)(𝑥),onΩ,(2.2)𝜕𝑛=𝑐¹(𝑥,𝑡),on𝜕Ω×[0,𝑇],(2.3)where 𝑐 is the molar concentration, 𝐷 is the diffusion parameter, and 𝑅_𝑔 is the reaction and source term.

We modify our model equation (2.1) to a control problem with an additionally right-hand side source:𝜕_𝑡𝑐−∇𝑐𝐷(𝑥∇,0𝑐)==𝑐𝑐_source0,inΩ×[0,𝑇],𝜕𝑐(𝑥,𝑡)(𝑥),onΩ,𝜕𝑛=𝑐¹(𝑥,𝑡),on𝜕Ω×[0,𝑇],(2.4)where 𝑐_source(𝑥,𝑡)=𝑅_𝑔 is the discontinuous or continuous source term of the concentration 𝑐 and we neglect a reaction term of this concentration.

We assume an optimal concentration at the layer𝑐_opt(𝑥,𝑡),(2.5) where the layer is given as 𝑥∈Ω_layer and our constraints are given as𝑐_source,min≤𝑐_source≤𝑐_source,max.(2.6)

Additionally, we have to solve the minimization problem:min(𝐽(𝑐,𝑐_source))∶=12_𝑇_Ωlayer||𝑐(𝑥,𝑡)−𝑐_opt(𝑥,𝑡)||²𝑑𝑥𝑑𝑡+𝜆2_𝑇_Ω||𝑐_source(𝑥,𝑡)||²𝑑𝑥𝑑𝑡,(2.7)where 𝑇 is the time period of the process.

Remark 2.1. We choose the 𝐿₂-error to control our minimization problem. In literature, see [15, 16], there exists further control-errors, which respect the time behavior.

In a first part, we only solve the transport equation with UG software-tool (unstructed grid software, see [17]) and try to find out the optimal control of the sources to obtain the best homogeneous layer.

In a second part, we consider the optimal control problem and solve also the backward problem.

3. Theoretical Background for Simulation of Diffusive CVD Processes

In what follows we discuss the approximation methods and errors for the simulation of the CVD processes.

3.1. Approximation and Discretization

For the numerical solutions we need to apply approximation methods, for example, finite-difference methods and iterative solver methods for the nonlinear differential equations, see [18, 19].

The finite-element discretization is based on Ω_ℎ the variational boundary value problem reduces to find 𝑢_ℎ∈𝑉_ℎ satisfying the initial condition 𝑢_ℎ(0) such that_Ωℎ𝜕𝜕𝑢𝑡^ℎ𝑣^ℎ+𝐷∇𝑢_ℎ⋅∇𝑣_ℎ𝑑𝑥=0,∀𝑣^ℎ∈𝑉_ℎ.(3.1)We define the minimal length of triangle which we get from the spatial-discretization with Δ𝑥.

This leads to the following linear semidiscretized system of ordinary differential equations:𝑀𝑑𝑢^∗𝑑𝑡+𝐴𝑢^∗=0,(3.2)where 𝑀 is the mass and 𝐴 the M-matrix.

Here we have taken into account the Courant-Friedrichs-Levy- (CFL-) condition, which is given asCFL=2𝐷_maxΔ𝑡min_𝑒∈𝐸Δ𝑥^2𝑒,(3.3)where 𝐷_max is the maximal diffusion parameter, 𝐸 is the set of the edges of the discretization. We restrict the CFL-condition to 1, if we use an explicit time-discretization and can lower the condition, if we use an implicit discretization.

For the explicit time-discretization, we apply explicit Euler or Runge-Kutta methods.

We use the explicit lower-order Runge-Kutta methods:0121201(3.4)

Furthermore we use the following Heun method (third-order):0131233023014034(3.5)

The implicit time-discretization is done with implicit Euler or Runge-Kutta methods.

Here, we use the implicit trapezoidal rule:0112211221(3.6)

Furthermore we use the following Gauss Runge-Kutta method:12−√361414−√3162+√3614+√36411221(3.7)

Remark 3.1. We apply implicit time-discretization methods for the pure diffusion part, where we apply explicit time-discretization methods for the pure convection part. Here we have to respect the CLF-condition, see [20].

3.2. Errors and Convergence Rate

For studying the errors and the convergence-rates in our test example, we have to define the following norm in two space-dimensions:

(i)discrete 𝐿_max-norm:err_{𝐿max,Δ𝑥,Δ𝑡}=m^𝑝ax_𝑖=1||𝑐_num(𝑥_𝑖,𝑇)−𝑐_ref(𝑥_𝑖,𝑇)||,(3.8)(ii) discrete 𝐿₁-norm: err_{𝐿1,Δ𝑥,Δ𝑡}=^𝑝_𝑖=1Δ𝑥²||𝑐_num(𝑥_𝑖,𝑇)−𝑐_ref(𝑥_𝑖,𝑇)||,(3.9)(iii) discrete 𝐿₂-norm: err_{𝐿2,Δ𝑥,Δ𝑡}=⎷^𝑝_𝑖=1Δ𝑥²||𝑐_num(𝑥_𝑖,𝑇)−𝑐_ref(𝑥_𝑖,𝑇)||²,(3.10) where Δ𝑥 is the spatial-step of the discretization, Δ𝑡 is the time-step of the discretization, and 𝑇 is the end-time of the computation. 𝑝 is the number of grid points in the discretization method. 𝑐_num is the numerical solution and 𝑐_ref is the reference solution, computed at fine spatial- and time-grids.

The numerical convergence rate are given as follows.

(i)For the spatial error, we define 𝜌_{𝐿2,Δ𝑥1,Δ𝑥2,Δ𝑡}=logerr^{𝐿2,Δ𝑥1,Δ𝑡}/err_{𝐿2,Δ𝑥2,Δ𝑡}log(Δ𝑥₁/Δ𝑥₂),(3.11)where Δ𝑥₁ is the coarse, Δ𝑥₂ is the fine spatial grid-step, and Δ𝑡 is the time-grid step for both results.(ii) For the time error, we define 𝜌_{𝐿2,Δ𝑥,Δ𝑡1,Δ𝑡2}=logerr^{𝐿2,Δ𝑥,Δ𝑡1}/err_{𝐿2,Δ𝑥,Δ𝑡2}log(Δ𝑡₁/Δ𝑡₂),(3.12) where Δ𝑡₁ is the coarse, Δ𝑡₂ is the fine-time-step, and Δ𝑥 is the spatial-grid step for both results.

We often use Δ𝑥₂=Δ𝑥₁/2. In this case, we have 𝜌_{𝐿2,Δ𝑥1,Δ𝑥2,Δ𝑡}=𝜌_{𝐿2,Δ𝑥1,Δ𝑡}. Further we have to choose Δ𝑥₁ with respect to the Δ𝑡∈𝐼=[0,Δ𝑡_max], which have maximal 𝜌. Thus we define ArgMax(Δ𝑥):ArgMax(Δ𝑥)∶=argmax_{Δ𝑡∈𝐼}𝜌_{𝐿2,Δ𝑥,Δ𝑥/2,Δ𝑡}.(3.13)

4. Optimal Control Methods

Here we discuss the control of a diffusion equation with a feedback based on a PID-controller.

4.1. Forward Controller (Simple P-Controller)

The first controller we discuss is the simple P-controller, see [11]. A first idea is to control linearly the error of the solved PDE.

In Figure 3, we present the P-controller.

Figure 3 

P-controller for the solution 𝐶.

Our control problem is given with the control of the error to the optimal concentration of the layer and correct the source-flux:𝜕_𝑡𝑐−∇𝑐𝐷(𝑥∇,0𝑐)==𝑐𝑐_source0,inΩ×[0,𝑇],𝜕𝑐(𝑥,𝑡)(𝑥),onΩ,𝜕𝑛=𝑐¹(𝑥,𝑡),on𝜕Ω×[0,𝑇],(4.1)where 𝑐_source(𝑥,𝑡) is a discontinuous source flow of the concentration 𝑐.

We assume an optimal concentration at the layer with the concentration 𝑐_opt(𝑥,𝑡), where the layer is given as 𝑥∈Ω_layer and 0 elsewhere.

Our constraints are bounded as𝑐_source,min≤𝑐_source≤𝑐_source,max.(4.2)

Remark 4.1. Taken into account the hysteresis of the deposition process, we apply a linear increase of 𝑐_source,max in the optimal control with respect to time, see Figure 4.

Figure 4 

Linear constraint 𝑐opt for the deposition process. x-axis: Time, y-axis: c (concentration).

4.2. PID-Controller

The PID-controller is used to control temperature, motion, and flow. The controller is available in analog and digital forms, see [16]. The aim of the controller is to get the output (velocity, temperature, position) in the area of the constraint output, in a short time, with minimal overshoot, and with small errors.

We have three elements in the PID-control, where 𝑃 is the proportional part, 𝐼 the integral part, and 𝐷 is the derivative part of the controller, see [16].

These terms describe three basic mathematical functions applied to the error signal, error=𝐶_optimal−𝐶_computed.

The errors represented the difference between constraint (optimal set) and computed results in the simulation.

To accelerate a PID-controller means to adjust the three multipliers 𝐾_𝑃, 𝐾_𝐼, and 𝐾_𝐷 adding in various amounts of these functions to get the system to behave the way you want, see [11].

Figure 5 summarizes the PID terms and their effect on a control system.

Figure 5 

PID-control. Effect on control system: the main influence in a control loop, 𝐾𝑃 reduces a large part of the overall error. 𝐾𝐼 reduces the final error in a system. Summing even a small error over time produces a drive signal large enough to move the system toward a smaller error. 𝐾𝐷 counteracts the 𝐾𝑃 and 𝐾𝐼 terms when the output changes quickly. This helps reduce overshoot and ringing.

Initialization of the PID-Controller
The algorithm of the initialization of the PID-control (i.e., search 𝐾_𝑃,𝐾_𝐼,𝐾_𝐷) is given as in Algorithm 4.2 (see [15]).

Algorithm 4.2. (1) We initialize the P-controller: 𝐾_𝐼=0.0,𝐾_𝐷=0.0.
(2) The amplifying factor 𝐾_𝑃 is increased till we reached the permanent oscillations as a stability boundary of the closed control system.
(3) We obtain for 𝐾_𝑃 the critical value 𝐾_𝑃,crit..
(4) The period-length of the permanent oscillation is given as 𝑇_crit.
(5) We obtain the parameters from Table 1.

Table 1 

Heuristic derivation of the control parameters (Nichols-Ziegler).

Further we compute the rest parameters as 𝐾_𝐼=𝐾_𝑃/𝑇_𝑛,𝐾_𝐷=𝐾_𝑃∗𝑇_𝑣, see [21].

4.3. Adaptive Time-Control

Often the heuristic assumptions of the PID-parameters are too coarse.

One can improve the method by applying an adaptive step-size control.

We discuss the step-size control with respect to our underlying error, that is, given by the computed and optimal output of our differential equation.

Based on the adaptive control, we can benefit to accelerate the control problem.

According to Hairer and Wanner [22], we apply the automatic control problem with a PID-controller.

The automatically step-size is given as (see [11])Δ𝑡_𝑛+1=𝑒^𝑛−1𝑒_𝑛^𝐾𝑃t𝑒o_𝑛l^𝐾𝐼𝑒𝑒_𝑛𝑒^2𝑛−1_𝑛−2^𝐾𝐷Δ𝑡_𝑛,(4.3)where tol is the tolerance, 𝑒_𝑛 is the error of the quantities of interest in time-step Δ𝑡.

We can control the step-size with respect to our heuristically computed 𝐾_𝑃, 𝐾_𝐼, and 𝐾_𝐷 parameters. Initialization of the adaptive control can be seen in Algorithm 4.3.

Algorithm 4.3. (1) Define Tolerance, Min and Max of the concentration.
(2) Apply the parameters: 𝐾_𝑃,𝐾_𝐷,𝐾_𝐼 form a first run.
(3) Optimize the computations with a first feedback.

4.4. Identification of the Control Path

In the forward problem, we computed a PDE to simulate the CVD process in the underlying domain. For our control problem, we have the behavior of two points in the underlying domain, 𝐶_source our source and 𝐶_𝑥 the source restricted to the response of the controlled system. To analyze the differences between given source the 𝐶_source and the computed source in the control process 𝐶_𝑥, we study the step response from our system:𝑥_𝑎(𝑡)∶=𝐶_𝑥(𝑡)∣𝐶_Source(𝑡)=𝑥_𝑒0.(4.4)

Algorithm 4.4. (1) We determine the model of the control path, for example, PT1 (Proportional time 1), PT2 (Proportional time 2), see [16]. For that purpose we investigate the step response 𝑥_𝑎, his first and second derivate [16, pages 117 and 331].
(2) We determine the parameters of our model: Kp, T1, T2.
(3) Our goal is to control the system with a controller. Also we have to determine the control-parameter [16, page 405].

Identification of the Control Path
With the step response 𝑥_𝑎, we identify the control path as a PT2 element.
Because lim_𝑡→∞𝑥_𝑎(𝑡)=constant (also ≠0 nor ∞), the basic behavior of our element is 𝑃 (and neither 𝐷 nor 𝐼). Now we have to determine the number of time delays in the control path.

Identification of the Parameters
For such control path which has PT2 behavior, there exists a model as ordinary differential equation. Here a generalization can be obtained with four parameters and be divided by 𝑎₀ a normed form with three parameters 𝜔₀= angular frequency, 𝐴 = attenuation, and 𝐾_𝑆 = proportional factor:𝑎₂1⋅𝑥̈_𝑎(𝑡)+𝑎₁⋅𝑥̇_𝑎(𝑡)+𝑎₀⋅𝑥_𝑎(𝑡)=𝑏₀⋅𝑥_𝑒0,𝜔²⁰⋅𝑥̈_𝑎(𝑡)+2𝐴𝜔₀⋅𝑥̇_𝑎(𝑡)+𝑥_𝑎(𝑡)=𝐾_𝑆⋅𝑥_𝑒0.(4.5)The solution of the differential equation results in the characteristic equation:𝛼²𝜔²⁰+2𝐴⋅𝛼𝜔₀+1=0with𝛼_1,2=−𝐴⋅𝜔₀±𝜔₀⋅√𝐴²−1.(4.6)
So that we have the analytical solution for the normed step response:𝑥_𝑎(𝑡)𝑥_𝑒0=𝐾_𝑆⋅1−𝑇¹𝑇₁−𝑇₂⋅𝑒^{−𝑡/𝑇1}+𝑇²𝑇₁−𝑇₂with𝛼₁⋅𝑒^{−𝑡/𝑇2},=−1𝑇₁,𝛼₂=−1𝑇₂.(4.7)We go now the inverse way to identify the system. We start with the step response. Using the inflected tangent at (𝑡_𝑤,𝑥_𝑤,𝑥̇_𝑎(𝑡_𝑤)) and the limit 𝑥_𝑎(𝑡→∞), we obtain our parameters.
Firstly, we have 𝐾_𝑆=𝑥_𝑎(𝑡→∞)/𝑥_𝑒0, 𝑡_𝑔=𝐾_𝑆⋅(𝑥_𝑒0/𝑥̇_𝑎(𝑡_𝑤)), and 𝑡_𝑢=𝑡_𝑤−𝑡_𝑔⋅(𝑥_𝑎(𝑡_𝑤)/𝑥_𝑎(𝑡→∞)), see Figure 6. In the second step, we obtain 𝑇₁ and 𝑇₂ by the following formulas:𝑡_𝑢𝑡_𝑔=𝛼^{𝛼/(1−𝛼)}⋅(𝛼⋅ln(𝛼)+𝛼²−1)𝛼−1−1,𝑇𝑡^𝑔₁=𝛼^{𝛼/(1−𝛼)},𝛼=𝑇²𝑇₁.(4.8)

Figure 6 

Step response and inflection tangent.

Figure 7 

Empirical determination of the model of the control path. The three main cases: (1) 𝑥𝑎(𝑡=0)≠0∶ P,PDT1,PPT1,PIDT1, (2) 𝑥𝑎(𝑡=0)=0,𝑥̇𝑎(𝑡=0)≠0∶ I,PT1,PDT2, (3) 𝑥𝑎(𝑡=0),𝑥̇𝑎(𝑡=0)=0, 𝑥̈𝑎(𝑡=0)≠0∶ IT1,PT2,…., (1) 𝑥𝑎(0)=0, (2) 𝑥̇𝑎(0)=0, (3) 𝑥̈𝑎(0)≠0.

Control Parameters
The parameter occurrence also in the transfer function 𝐺(𝑠), which is the Laplace transformation of the differential equation. Our first choice to determine the control parameter is here the PT2 element. Alternatively, we use an approximation to a PT1 element, but here we have additional a scan-time parameter:𝐺(𝑠)=𝐾^𝑆𝐺(𝑠)≈𝐾(1_𝑆+𝑇₁𝑠)(1+𝑇₂𝑠)PT2element,⋅𝑒^{−𝑇𝑢𝑠}1+𝑠⋅𝑇_𝑔PT1element.(4.9)Independent from the behavior which is being wished, we obtain values for the parameters 𝐾_𝑃, 𝐾_𝐼, and 𝐾_𝐷, respectively. 𝐾_𝑅=𝐾_𝑃, 𝑇_𝑛 and 𝑇_𝑣 are of the PID-controller. In Table 2 are given the values in dependence to the percentage overshoot (PO). There are the values of PO=0, respectively, PO=20. The rest can be determined by linear interpolation. An alternative schema is given by Takahashi. Here we have in addition the parameter scan-time 𝑇. This is a general form of the Nichols-Ziegler method (see Table 1). Further are 𝐾_𝑃=𝐾_𝑅, 𝐾_𝐼=𝐾_𝑃/𝑇_𝑛, and 𝐾_𝐷=𝐾_𝑃∗𝑇_𝑣.

Table 2 

Derivation of the control parameters (Chien-Hrones-Reswick and Takahashi).

Table 3: MATLAB-toolbox pdetool.

5. Software and Program-Tools

MATLAB Functions
We use the MATLAB-toolbox pdetool for the time- and spatial-discretizations, where we have a finite-element method with 𝑃1-elements and an implicit Euler method for the time-discretization.

DEPOSIT-PID Toolbox
The PID-controller is also programmed in MATLAB. Our combined code is given in the DEPOSIT-PID toolbox and described in what follows. The DEPOSIT-PID toolbox is manipulated by a graphical user interface, by which simulation- and control-models are chosen and corresponding parameters can be manually adjusted.
The objective is to simulate the diffusion and deposition of the vapor in the apparatus and to obtain the optimal vapor concentration at the measuring point.
We can divide the process into three phases, namely, forward step, control step, backward step, which have a cyclic repetition. In the forward step, the diffusion takes place, which can be simulated by a time-step of the heat equation.
After that, in the control step, the actual concentration at the under boundary (0,−1) can be measured (shown by computed in the graph) and compared with the optimal value (shown by optimal).
The control step is followed by the backward step: from the error in control step and the control model, the optimal alteration of the source can be computed. The vapor flows through the source point (0,0) in the apparatus and this value is shown by the SourceOutput in Figure 12. This will be simulated through the addition of SourceOutput to actual concentration at the source point (see Figure 8). Further, as a simplification, we set the concentration at the under boundary to zero, because here the gas transforms into solid matter.

(a)

(b)

Figure 6 

Backward step.

Figure 9 

Layout of the DEPOSIT-PID Gui.

Figure 10 

Convergence diagram tool.

In Section 4 we introduced some fundamentals in order to control the apparatus. The models and corresponding parameters can be altered by the Gui.

The Layout of the DEPOSIT-PID Gui
This Gui contains the following:
(1)short-time plot (2D) of computed, optimal, and SourceOutput; (2) long-time plot (2D);(3)listbox with names parameters:(a)textbox with actual value of the parameter chosen in [12], (b)textbox to change the value of the parameter chosen in [12]; (4)listbox with names parameters:(a)textbox with actual value of the parameter chosen in [23], (b)textbox to change the value of the parameter chosen in [23]; (5)3D plot of distribution;(6)3D grid plot of distribution;(7)listbox with names of parameters: (a)checkbox with actual value of the parameter chosen in [24];(8) push button: save;(9) push button: reset;(10)radio button: start.

KONTOOL
Another software-tool, KONTOOL, is programmed to compute the numerical convergence rates of the applications. The software-tool has implemented the errors and convergence rates defined in Section 3.2. An error-analysis based on successive refinement of space and time is done and the resulting errors and convergence rates are computed. Optimal convergence rates with respect to balance the time- and spatial-grids are calculated, see Algorithm 5.2.

Remark 5.1. The software-tool can be modified and applied to arbitrary spatial- and time-discretization methods. The interface of KONTOOL needs at least the parameters of the spatial- and time-grid and the starting parameters of the underlying methods.

Algorithm 5.2 (The algorithm of the computation of the numerical convergence tableau). (1) We compute reference solutions: (a) numerically: fine time and spatial steps or (b) analytically (if there exists an analytical solution).
(2) We apply one spatial discretization of step Δ𝑥 and apply all time discretization with steps Δ𝑡, where the coarsest Δ𝑡 is given by the CFL condition or till first non-numerical results as oscillations. We compute the error 𝑐_num−𝑐_ref in the 𝐿₂-norm.
(3) We continue the next fine spatial steps, for example, Δ𝑥/2.
(4) We compute the convergence tableau with time and space.

In the next section, we discuss the numerical experiments.

6. Experiment for the Plasma Reactor

In this section, we present our numerical experiments for the CVD processes in a plasma reactor.

6.1. Simulation of a Diffusion Equation with Analytical Solution (Neumann Boundary Conditions)

Here we simulate a diffusion equation with Neumann boundary conditions and right-hand side 0. Our control problem has only the forward problem to solve and we consider the accuracy of our simulations.

We have the following equation:𝜕𝑐_𝑡(𝑐𝑥−𝜕,𝑦𝛽,²0()𝜕=_𝑥𝑥c+os𝜕(2_𝑦𝑦)𝑥𝑐)=+0co,s(2in𝑦)Ω,×o[n0,Ω𝑇],_𝑛,𝑐(𝑥,𝑦,𝑡)=0,on𝜕Ω×[0,𝑇],(6.1)where 𝑐 is the molar concentration, Ω=[−1,1]×[−1,1], and 𝑡∈(0,𝑇). 𝐷=𝛽² is the diffusion parameter of the diffusion equation.

We have the following analytical solution:𝑐_ana(𝑥,𝑦,𝑡)=sin(2)+^∞_𝑛=1𝐴_𝑛exp(−𝛽²𝑛²𝜋²𝑡)(cos(𝑛𝜋𝑥)+cos(𝑛𝜋𝑦)),(6.2)where 𝐴_𝑛=(−1)^𝑛(−4sin(2)/(𝑛²𝜋²−4)). We apply the diffusion coefficient 𝐷=0.01, respectively, 𝛽=0.1. We obtain the following result after time 𝑇=12.0, see Figure 11.

(a)

(b)

(c)

(d)

Figure 8 

2D experiment of the diffusion equation at the end-time 𝑇=12.0, 𝛽=0.1. x-axis: x-coordinate, y-axis: y-coordinate and z-axis: c (concentration).

(a)

(b)

(c)

(d)

Figure 9 

2D experiment of the diffusion equation and control of a single point.

We see that for large Δ𝑥 the numerical solution converge faster to the stable constant endsolution (offset) than solution for smaller Δ𝑥 and the analytical solution. The error between the constant analytical endsolution (offset, sin(2)=0.9093) is also greater than for smaller Δ𝑥. The 𝐿₂-error is given in Table 4.

Table 4 

Offset convergence (Δ𝑡=0.1, time-step=100, 𝛽=0.61644).

Table 5 

Numerical results for the P-controller for different spatial-steps with 𝐷=0.1 and 𝜆=1.0 as P-value for the controller.

Table 6 

The associated Δ𝑡, which have maximal convergence rate (ArgMax), for Δ𝑥=0.1, 0.05, 0.025, 0.0125.

Table 7 

The associated Δ𝑡, which have maximal convergence rate (ArgMax), for Δ𝑥=0.1, 0.05, 0.025, 0.0125.

Table 8 

Convergence of 𝜆 with bisection.

Remark 6.1. We test for the pure diffusion equation our underlying discretization methods and apply finite elements for the spatial-discretization and implicit Runge-Kutta methods for the time-discretization. In the results, we obtain decreasing errors for the different time- and spatial-steps.

6.2. Simulation of an Optimal Control of a Diffusion Equation with Heuristic Choise of the Control Parameters

Here we simulate a first example of a diffusion equation and control the concentrations in the deposition process.

We have the following equation:𝜕_𝑡𝑐−∇𝑐(𝐷𝑥,∇0𝑐)==𝑓𝑐₀(𝑡),inΩ×[0,𝑇],𝜕𝑐(𝑥,𝑡)(𝑥),onΩ,𝜕𝑛=𝑐¹(𝑥,𝑡),on𝜕Ω×[0,𝑇],(6.3)where 𝑐 is the molar concentration, 𝐷 is the diffusion parameter of the diffusion equation, and 𝑓(𝑡) is the right-hand side or diffusion source. We have the following constraint: 𝑐_optimal(𝑥_point,𝑦_point)=0.5, where (𝑥_point,𝑦_point) is the control point in our domain. The parameters are given as 𝐷=0.01 and 𝑓(𝑡)=𝑎𝑡+𝑏, so we deal with a linear source, 𝑎=0.2 and 𝑏=0.1 are constants. In the following tests, we propose the 3 possibilities to control the optimal temperature:

(i)P-control with constant optimal constraint,(ii)PID-control with constant optimal constraint,(iii)PID-control with Linear optimal constraint.

The results for the control methods are given in Figure 12.

6.3. Preliminary Remark for Simulations for the Convergence Order

To determine the function ArgMax, firstly, we have to determine a practicable interval 𝐼 forArgMax(Δ𝑥)∶=argmax_{Δ𝑡∈𝐼}𝜌_{𝐿2,Δ𝑥,Δ𝑥/2,Δ𝑡}.(6.4)One possibility is the interval𝐼_CFL∶=(0,Δ𝑡_CLF],(6.5)with Δ𝑡_CFL∶=Δ𝑠²/2𝐷_max from the CFL-condition. In the experiments, we find another interval, where the convergence rate function is relative stable convex:𝐼_stable∶=(0,Δ𝑡_stable],(6.6)whereΔ𝑡_stable∶=max_Δ𝑡{𝜌∶(0,Δ𝑡)⟶ℝisconvex}.(6.7)It is clear that we can in 𝐼_stable take some restrictions to a subinterval 𝐼_sub∶=[Δ𝑡_𝑠min,Δ𝑡_𝑠max]⊂𝐼_stable, when we know that some Δ𝑡∈𝐼_sub with 𝜌(Δ𝑡)>𝜌(Δ𝑡_𝑠min),𝜌(Δ𝑡_𝑠max).

The function of the numerical convergence rate is discrete in the spatial-discretization variable Δ𝑥 since we get a finer discretization with a bisection of Δ𝑥. With a finer discretization, a triangle is replaced by four subtriangles.

In the temporal discretization variable Δ𝑡, we are in contrast not restricted to such conditions and Δ𝑥 could be chosen to any arbitrary value above 0. To get a first glance, we have selected the methods of bisection. Subsequently, we consider finer discretizations in intervals which are of special interest.

6.4. Simulations for the Convergence Order

We consider the convergence rate 𝜌 as a function in dependent of Δ𝑥 and Δ𝑡. This function also depends on some parameters, for instance, the diffusion coefficient 𝐷 and the control parameter 𝜆. We now present the results of two chosen experiments. For the first experiment, the parameters are 𝐷=0.1 and 𝜆=1, while we increase the propagation velocity in such a way that we increase 𝐷 to 1 for the second experiment.

We observe that, for instance, in the case of spatial-discretization Δ𝑥=0.05 and 𝑥=0.25, the maximal convergence rate lies between Δ𝑡=0.05 and Δ𝑡=0.0125. The maximum itself lies close to 0.025. To determine the precise value of Δ𝑡=ArgMax(Δ𝑥), we refine our method in these interesting intervals.

In Figure 13, we observe that there is a relatively stable convex area starting at 0. This area then switches over to an area with strong fluctuations. We now compare the stable area with the CFL-condition. The area where CFL<1 lies in the stable convex area of the function (Δ𝑡_CFL<Δ𝑡_stable). For Δ𝑥=0.05,𝐷=0.1, the CFL-area ends at Δ𝑡_CFL=0.0125, the stable convex area range to about Δ𝑡_stable=0.2. For our convergence diagram, we try to find for every Δ𝑥 the value Δ𝑡, where the convergence rate becomes maximal (ArgMax). It is clear that we only use Δ𝑡 at which the convergence rate function is stable: Δ𝑡<Δ𝑡_stable.

(a)

(b)

Figure 10 

𝜌 for a P-controller, 𝐷=0.1.

Figure 14 

Convergence diagram (KONTOOL). We can see in the loglog-plot ArgMax(Δ𝑥) linear dependence.

Figure 15 

Convergence diagram (KONTOOL). The loglog-plot ArgMax(Δ𝑥) with linear dependence.

(a)

(b)

(c)

(d)

Figure 11 

𝜌 for a P-controller, 𝐷=1.

Remark 6.2. The experiment shows the linear convergence rate of the P-controller with different 𝜆 values. So we obtain a stable method with respect to the P-controller. In the examples, we apply heuristic methods to derive the control parameters for the P- and PID-controller. We show that we have reached the linear order of the underlying finite element discretization method. We have higher control errors if we did not compute the correct control parameters and the numerical errors are smaller than our control error. To prohibit this problem, we have to compute in the next example the control parameters by a feedback equation, see [25].

6.5. Simulation of an Optimal Control of a Diffusion Equation with Adaptive Control

In the second example, we simulate the diffusion equation and control the temperature with and adaptive control based on a PID-controller, see [11].

We have the following equation:𝜕_𝑡𝑐−∇𝑐(𝐷𝑥∇,𝑡𝑐)==𝑐𝑓₀(𝑡),inΩ×[0,𝑇],𝜕𝑐(𝑥,𝑡)(𝑥),onΩ,𝜕𝑛=𝑐¹(𝑥,𝑡),on𝜕Ω×[0,𝑇],(6.8)where 𝑐 is the molar concentration, 𝐷 is the diffusion parameter of the diffusion equation, and 𝑓(𝑡) is the right hand side or source.

We have the following constraint: 𝑐_optimal(𝑥_point,𝑦_point)=0.5,(6.9) where (𝑥_point,𝑦_point)=(0,−1) is the control point in our domain.

The automatically step-size is given as (see [11])Δ𝑡_𝑛+1=𝑒^𝑛−1𝑒_𝑛^𝐾𝑃t𝑒o_𝑛l^𝐾𝐼𝑒𝑒_𝑛𝑒^2𝑛−1_𝑛−2^𝐾𝐷Δ𝑡_𝑛,(6.10)where tol=1 is the tolerance, 𝑒_𝑛 is the error of the quantities of interest in time-step Δ𝑡.

The errors are given as𝑒_𝑛=||𝑢^𝑛−𝑢_𝑛−1||||𝑢_𝑛||,(6.11)where 𝑢_𝑛 is the result at time-step 𝑡^𝑛.

The parameters are given as 𝐷=0.1, 𝑃_krit=15, 𝑇_krit=5, Δ𝑡 (PID-control), tol=1 (adaptive PID-control), Δ𝑡_max=0.1 (adaptive PID-control), Δ𝑡_min=0.01 (adaptive PID-control) (see Figure 17).

Figure 17 

2D experiment with and without the adaptive time-step control.

Furthermore, we change the parameter tol=1 to tol=0.1 and Δ𝑡_min=0.01 to Δ𝑡_min=0.0001 (see Figure 18).

(a)

(b)

Figure 13 

Adaptive PID with modified parameters.

Remark 6.3. In Figures 17 and 18, we see an oscillating time interval at the beginning of the automatically step-size control. In the first experiments, we had only taken into account a previous heuristically computation of the control parameters 𝐾_𝑃, 𝐾_𝐼, and 𝐾_𝐷 before the step-size control for the whole time-interval 0,𝑇. The optimal control parameters are given to the whole time-interval 0,𝑇. A modified algorithm to compute the control parameters for the initialization time-interval 0,𝑡_begin and the whole time-interval 0,𝑇 overcome the oscillation problems.

A modified automatically step-size control, which minimize the oscillations is given in the following algorithm.

Algorithm 6.4. (1) We compute the reference control parameters 𝐾_𝑃,global, 𝐾_𝐼,global and 𝐾_𝐷,global for the time-interval 0,𝑇.
(2) We apply the automatically step-size control for the global control parameters with tol_global, Δ𝑡_max,global and Δ𝑡_min,global, which can by chosen large.
(3) We stop the computation till we reach the optimal solution and mark remember the time 𝑡_oszill.
(4) We compute the local control parameters 𝐾_𝑃,local, 𝐾_𝐼,local and 𝐾_𝐷,local for the time-interval 0,𝑡_oszill.
(5) We restart the computation with the local control parameters and smaller step-size parameters tol_local, Δ𝑡_max,local and Δ𝑡_min,local till we reach 𝑡_oszill and continue the computation with the global parameters.
(6) We stop the computation if we reach 𝑡=𝑇. If we obtain also high oscillation with the local parameters, we refine the local interval and go to step (3).

Remark 6.5. The modified automatically step-size control had taken into account the local behavior of the control problem. We could adapt the control parameters 𝐾_𝑃, 𝐾_𝐼, and 𝐾_𝐷 with respect to the local time-intervals. This modified algorithm considers a local time behavior more accurate and reduces oscillations at the initialization process.

6.6. Simulation of an Optimal Control of a Diffusion Equation with Adaptive Control and Recovering of the Control Parameters

In what follows, we present the adaptive control based on Algorithm 6.4.

To be more precise the computation of global control parameters 𝐾_𝑃,global=𝐾_𝑃(Δ𝑡), 𝐾_𝐼,global=𝐾_𝐼(Δ𝑡), and 𝐾_𝐷,global=𝐾_𝐷(Δ𝑡) can be automatized in an interval Δ𝑡∈(0,𝑇).

We improve the automatically step-size, which is given in Section 4.3 to(Δ𝑡)_𝑛+1=𝑒^𝑛−1𝑒_𝑛^{𝐾𝑃(Δ𝑡)}t𝑒o_𝑛l^{𝐾𝐼(Δ𝑡)}𝑒𝑒_𝑛𝑒^2𝑛−1_𝑛−2^{𝐾𝐷(Δ𝑡)}(Δ𝑡)_𝑛,(6.12)where tol is the tolerance, 𝑒_𝑛 is the error of the quantities of interest in time-step 𝑛.

With approximations to each subintervals of (0,𝑇), we derive the control parameters 𝐾_𝑃(Δ𝑡), 𝐾_𝐼(Δ𝑡), and 𝐾_𝑃(Δ𝑡).

The improved globalized control parameters are computed in the following algorithm.

Algorithm 6.6. Computing of the control parameters 𝐾_𝑃,global=𝐾_𝑃(𝑡), 𝐾_𝐼,global=𝐾_𝐼(𝑡) and 𝐾_𝐷,global=𝐾_𝐷(𝑡) for the time-interval 0,𝑇.
(1) Compute the critical control function with the minima and maxima. Based on the optimal oscillated function, we can derive the 𝐾_𝑃(𝑇/3) for the interval (0,𝑇/3).
(2) Redo the step (1) for the intervals (𝑇/3,2𝑇/3) and 2𝑇/3,𝑇 and approximate the function 𝐾_𝑃(𝑡). Based on the 𝐾_𝑃(𝑡) we can derive the 𝐾_𝐼(𝑡) and 𝐾_𝐷(𝑡) function.
(3) The optimal control parameters are used in Algorithm 6.4.

Remark 6.7. Further ideas to determine the time-dependent parameters 𝐾_𝑃(Δ𝑡), 𝐾_𝐼(Δ𝑡), and 𝐾_𝐷(Δ𝑡) can be done by using our identification method from Algorithm 4.4 in Section 4.4. There we have discussed two methods to determine the parameters from the control path, namely, Chien-Hrones-Reswick and Takahashi. Alternatively, we have the method from Algorithm 4.2 (Nichols-Ziegler) in Section 4.2. For the identification of the control path, we have only look to the step response. The alternative method demanded to find the critical 𝑃-value.

6.7. Determinate of Critical P with Exponential Regression

To initialize the PID-controller, we have to derive the 𝐾_𝑃, 𝐾_𝐷, and 𝐾_𝐼 parameters.

An idea is derived with the step function response, see [15].

This idea has also be included into the time-step control of the PID-controller.

In Figure 19, we can see that the maximum and minimum of the error function, done with the step function response, have a exponential behavior.

Figure 19 

Error of the time-step control.

(a)

(b)

Figure 15 

Exponential regression of min and max. Behavior of 𝜆.

So we use the exponential regression𝑦=𝑏⋅exp(𝑎⋅𝑡)+𝜖,(6.13)with 𝜖∼𝑁(0,𝜎²) normal distributed. Via log the exponential regression can be formed in a linear regression:log𝑦(𝑦=)=𝑏⋅leoxgp(𝑏(𝑎)+⋅𝑡𝑎),⋅𝑡.(6.14)Before we have to convert the minima into the positive area.

Let 𝜆(𝑛) be the actual lambda value, 𝜆(𝑛)⁺ the next greater 𝜆, and 𝜆(𝑛)⁻ the next smaller one. Firstly, we use 𝜆(𝑛+1)=⎧⎪⎨⎪⎩(𝜆(𝑛)+𝜆(𝑛)⁺)(𝜆(𝑛)+2𝜆(𝑛)₋∶𝑎(max,𝑛)<0,𝑎(min,𝑛)>0(𝜆tosmall),)2∶𝑎(max,𝑛)>0,𝑎(min,𝑛)<0(𝜆togreat).(6.15) In the first step and in limit case, we use 𝜆(𝑛+1)=2⋅𝜆(𝑛), respectively, 𝜆(𝑛+1)=𝜆(𝑛)/2. In the first experiment, we start with 𝜆=359, which we have found with manual tests. The automatic found that the value 396.729 the 𝑎-value of 𝜆=3.59 have 3 zeros, the new less than five.

Remark 6.8. The adaptive step-size control of the PID-controller had taken into account the optimization in each local time-interval. We consider a larger time-interval and derive the optimal control parameters 𝐾_𝑃,global, 𝐾_𝐼,global, and 𝐾_𝐷,global. Based on this parameterization, we apply the adaptation in time to be optimal in such a large time-interval.

6.8. Results of Our Experiments to Find the Time-Dependent Parameters

To automatize the adaptive time-step algorithm, we have to compute the step response and derive the 𝐾_𝑃, 𝐾_𝐼, and 𝐾_𝐷 parameters of this process. Here, we present the experiments in Figure 21, which had taken into account the critical values 𝑃_crit and 𝑇_crit to derive the optimal 𝐾_𝑃, 𝐾_𝐷, and 𝐾_𝐼 parameters based on Algorithm 6.6.

(a)

(b)

(c)

Figure 16 

Δ𝑡 dependence parameters determined by Nichols-Ziegler.

(a)

(b)

(c)

Figure 17 

Δ𝑡 dependence parameters determined by Chien-Hrones-Reswick with PO=0,…,20.

A further algorithm to derive the control parameters is done by Chien et al., see [12].

The results are given in Figure 23.

(a)

(b)

(c)

Figure 18 

Δ𝑡 dependence parameters determined by Takahashi.

The 𝐾_𝑃,𝐾_𝐷, and 𝐾_𝐼 parameters are derived with Algorithm 6.6.

6.9. Experiments with Time-Dependent Parameters

In these experiments, we concentrate on the Chien-Hrones-Reswick method, see [12].

To automatize the adaptive time-step algorithm, we have to compute the step response and derive the 𝐾_𝑃, 𝐾_𝐼, and 𝐾_𝐷 parameters of this process. The minimization problem is presented in Figure 24, which had taken into account the critical parameters.

(a)

(b)

(c)

(d)

Figure 19 

Parameters chosen by Chien-Hrones-Reswick: Summands of the minimization problem.

In Figure 25, we see the step function response and the best fit of this function with the adaptive time control method. Because of the sensitivity to the initial problem, we can overcome the oscillation with the Chien-Hrones-Reswick method.

Figure 25 

Parameters chosen by Chien-Hrones-Reswick: Actual value.

In Figure 26, we see the comparison of our different adaptive methods. Here the best method is the benefit of the time adaptive control method based on Algorithm 6.6.

Figure 26 

Parameters chosen by Chien-Hrones-Reswick: Source.

Remark 6.9. The benefit of the time-dependent control parameters 𝐾_𝑃(Δ𝑡), 𝐾_𝐷(Δ𝑡), and 𝐾_𝐼(Δ𝑡) are important in the initialization of the control problem. Based on the fine-time scales globally chosen, parameters neglect the local time behavior. In the experiments, we found out the accelerated relaxation behavior and the fast control to the optimal value. Such adaptive algorithms to optimize the control parameters can improve the control methods and had taken into account the local-time scales.

7. Conclusions and Discussions

We present a continuous or kinetic model, due to the fare-field or near-field effect of our deposition process. We discuss the PID-controller to automatize our deposition process. Due to heuristic methods of deriving the PID parameters, we discuss aposteriori error estimates to automatize the time-stepping methods. A modified automatically step-size control is discussed and the best approximations are obtained with the time-dependent control method based on the Chien-Hrones-Reswick algorithm. For the mesoscopic-scale model, we discussed different experiments and their convergence rates. In future, we will analyze the validity of the models with physical experiments.