(182a) A Combined Graphical / Algebraic Method for Model Reduction and Analysis of Chemical Reaction Networks: Application to Atomic Layer Deposition Process

In many cases, the study of chemical engineering
systems involves mathematical modeling of a potentially complex reaction
network (RN) consisting numerous species and reactions. One example of such
systems is atomic layer deposition (ALD), a thin-film manufacturing technique based
on sequential reactions between gaseous precursors and a solid substrate to conformally
deposit thin-film with accurately controlled composition. An ALD process starts
with precursor adsorption on the substrate surface and proceeds with surface
reactions by different mechanisms such as dissociation, transition state formation,
ligand exchange, and densification; this reaction sequence ends with desorption
and ultimately evacuation of reaction by-products from the reactor [1]. These
reactions may occur with significantly different rates. For example in many
cases the initial adsorption of precursor molecules on surface sites is
barrierless and has a very short timescale; likewise equilibrium reactions result
in pseudo-steady state conditions. Some ligand exchange surface reactions, on
the other hand, may have a large energy barrier making them finite-rate
processes.

One path to formulating the molar-balance system of
equations for the species present in the RN system is to represent each
equilibrium reaction in form of a finite-rate process by introducing an
artificial time constant (1/ɛ), where ɛ can later be set to 0 to recover
the original equilibrium relationship [1,2]. Due to interconnected nature of
the reactions in the network, the system of differential equations obtained
from writing molar balances in this manner sometimes produces a singularly
perturbed system in non-standard form, meaning ɛ cannot be directly set to
0 to recover each equilibrium process; by doing so, information regarding the
dynamics of a subset of chemical species would be lost. Moreover, if the number
of independent reactions in the RN under study is less than number of species, the
RN is expected to have invariant states giving rise to additional algebraic
constraints. These issues can be systematically addressed by finding a
transformation of species molar amounts (original variables) to a new state
space in a process that identifies redundant dynamic modes (invariants) and
decouples reactions with different time scales to produce a well-posed
differential algebraic system of equations (DAE) upon setting ɛ→0 [1,2,3].

We recently proposed a graphical approach to identifying
invariant states associated with an RN [4]. Our approach makes use of
species-reaction (SR) graphs first introduced by Craciun and Feinberg
[5]. We have shown that a RN’s invariant states correspond to paths on the SR
graph and that these paths can be generated by a relatively simple set of rules
(cite our paper with KP, Aisha). In this work, we introduce a computational procedure
to facilitate identification of these paths. This procedure follows the
mathematical algorithm proposed by Schuster and Hofer [6] and is based on
constructing a series of tableaus based on the stoichiometric matrix
constructed from RN. The invariant states obtained from this graphical /
algebraic method are all semi-positive (nonzero and non-trivial linear
combination of species molar amounts); the semi-positive properties of the
invariants will be shown to be crucial to interpreting the physical
significance of each invariant state; furthermore, the properties of the invariants
themselves also can be used to identify potential defects in the RN structure. We
will show that this procedure can be modified to not just obtain invariant
states, but also to decouple reactions with different time scales to complete the
transformation described earlier. Extension of our algorithmic and graphical
techniques to the study of reaction variants also will be presented.  

Fig1. SR graph associated with RN of P+S→A
(e0), A→S+F (f0) giving a simplified picture of one ALD half-reaction
(precursor P, site S, adsorbed species A and film material F) and possessing 2
invariant states.  One of the invariants, associated with conservation of
surface sites, can be formulated from the cyclic path shown in blue (left).
Variant state associated with the effect of finite rate reaction f0 can be
obtained from the P to A path shown in red (right). The node corresponding to
equilibrium reaction e0 is surrounded by two species nodes
with different stoichiometric signs, resulting in removal of e0 effect in the
formulated variant state and leaving the f0 node in decoupled form:

[1] Adomaitis R., 2016. J Vac. Sci. Technol. A 34,
01A104.

[2] Daoutidis, P., 2015. Surveys in diferential-algebraic
equations II., Springer International Publishing, New York, 69.

[3] Rodrigues,
D., Srinivasan, S., Billeter, J., Bonvin, D., 2015. Comput. Chem. Eng. 73, 23.

[4] Salami, H., Ramakrishnasubramanian, K.,
Adomaitis, R., 2017. Phys. Status Solidi B 1700091.

[5] Craciun, G., Feinberg, M., 2006. SIAM J. Appl.
Math. 66(4), 1321.

[6] Schuster, S., Hofer, T., 1991. J. Chem. Soc.
Faraday Trans. 87(16), 2561.