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 , 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 , we deal with a Navier-Stokes equation or with the convection-diffusion equation, see [5, 14], whereas for large Knudsen numbers , 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)
(a)
(b)
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 aswhere 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:where 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 where the layer is given as and our constraints are given as
Additionally, we have to solve the minimization problem: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 such thatWe 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:where is the mass and the M-matrix.
Here we have taken into account the Courant-Friedrichs-Levy- (CFL-) condition, which is given aswhere is the maximal diffusion parameter, is the set of the edges of the discretization. We restrict the CFL-condition to , 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:
Furthermore we use the following Heun method (third-order):
The implicit time-discretization is done with implicit Euler or Runge-Kutta methods.
Here, we use the implicit trapezoidal rule:
Furthermore we use the following Gauss Runge-Kutta method:
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 -norm:(ii) discrete -norm: (iii) discrete -norm: 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. is the numerical solution and 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 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 where is the coarse, is the fine-time-step, and is the spatial-grid step for both results.
We often use . In this case, we have . Further we have to choose with respect to the , which have maximal . Thus we define :
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.
Our control problem is given with the control of the error to the optimal concentration of the layer and correct the source-flux:where is a discontinuous source flow of the concentration .
We assume an optimal concentration at the layer with the concentration , where the layer is given as and elsewhere.
Our constraints are bounded as
Remark 4.1. Taken into account the hysteresis of the deposition process, we apply a linear increase of in the optimal control with respect to time, see Figure 4.
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, .
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.
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: .
(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 .
(4) The period-length of the permanent oscillation is given
as .
(5) We obtain the parameters from Table 1.
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])where 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, our source and the source restricted to the response of the controlled system. To analyze the differences between given source the and the computed source in the control process , we study the step response from our system:
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 (also 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:The solution of the differential
equation results in the characteristic
equation:
So that we have the analytical solution for the normed
step response: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 , , and , see Figure 6. In the second step, we obtain and by the
following formulas:
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: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 , respectively, . 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 .
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 -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 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 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)
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 in the
-norm.
(3) We continue the next fine
spatial steps, for example, .
(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 . Our control problem has only the forward problem to solve and we consider the accuracy of our simulations.
We have the following equation:where is the molar concentration, , and . is the diffusion parameter of the diffusion equation.
We have the following analytical solution:where . We apply the diffusion coefficient , respectively, . We obtain the following result after time , see Figure 11.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
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, ) is also greater than for smaller . The -error is given in Table 4.
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: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: where is the control point in our domain. The parameters are given as and , so we deal with a linear source, and 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 , firstly, we have to determine a practicable interval forOne possibility is the intervalwith from the CFL-condition. In the experiments, we find another interval, where the convergence rate function is relative stable convex:whereIt is clear that we can in take some restrictions to a subinterval , when we know that some with .
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 . 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 and , while we increase the propagation velocity in such a way that we increase to for the second experiment.
We observe that, for instance, in the case of spatial-discretization and the maximal convergence rate lies between and . The maximum itself lies close to . To determine the precise value of , we refine our method in these interesting intervals.
In Figure 13, we observe that there is a relatively stable convex area starting at . This area then switches over to an area with strong fluctuations. We now compare the stable area with the CFL-condition. The area where lies in the stable convex area of the function (). For , the CFL-area ends at , the stable convex area range to about . 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: .
(a)
(b)
(a)
(b)
(c)
(d)
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: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: where is the control point in our domain.
The automatically step-size is given as (see [11])where is the tolerance, is the error of the quantities of interest in time-step .
The errors are given aswhere is the result at time-step .
The parameters are given as , , , (PID-control), (adaptive PID-control), (adaptive PID-control), (adaptive PID-control) (see Figure 17).
Furthermore, we change the parameter to and to (see Figure 18).
(a)
(b)
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 . The optimal control parameters are given to the whole time-interval . A modified algorithm to compute the control parameters for the initialization time-interval and the whole time-interval 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 , and for the
time-interval .
(2) We apply the automatically step-size control for
the global control parameters with , and , which can by chosen large.
(3) We stop the computation till we reach the optimal
solution and mark remember the time .
(4) We compute the local control parameters , and for the
time-interval .
(5) We restart the computation with the local control
parameters and smaller step-size parameters , and till we reach 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 , , and can be automatized in an interval .
We improve the automatically step-size, which is given in Section 4.3 towhere is the tolerance, is the error of the quantities of interest in time-step .
With approximations to each subintervals of , 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 , and for the
time-interval .
(1) Compute the critical control function with the
minima and maxima. Based on the optimal oscillated function, we can derive the for the
interval .
(2) Redo the step (1) for the intervals and 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.
(a)
(b)
So we use the exponential regressionwith normal distributed. Via the exponential regression can be formed in a linear regression: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 In the first step and in limit case, we use , respectively, . In the first experiment, we start with , which we have found with manual tests. The automatic found that the value the -value of 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 , , and . 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 and to derive the optimal , , and parameters based on Algorithm 6.6.
(a)
(b)
(c)
(a)
(b)
(c)
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)
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)
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.
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.
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.