ON THE CALCULATION OF KINETIC CONSTANTS OF CHEMICAL REACTIONS

A. I. RUBAN

(Presented by the scientific seminar on automatic control of the Department of PMAKhP)

On the basis of the use of the linearization method*) a problem is solved of calculating the kinetic constants of homogeneous chemical reactions from experimental data. An example is given.

1. Introduction. For the proper organization of chemical-technological processes, it is first of all necessary to know the mathematical description of chemical reactions, which can be obtained only on the basis of a laborious, long-term experiment. In this case, several hypotheses about reaction mechanisms are formulated from the preliminary experiment. For each mechanism kinetic equations are compiled, which must give a complete representation of the rate of chemical transformation of the components in the process under study. The kinetic equations, which are systems of linear and nonlinear differential equations, include a number of unknown parameters that must be determined from experimental curves of the change in the concentrations of the reacting components over time. In the simplest cases reactions are studied under isothermal conditions, and the unknown constants are the coefficients of the reaction rates. In more complex situations the reaction order is also unknown, and the experiment is carried out under a variable temperature regime. At this stage of research it is necessary to have mathematical methods that make it possible quickly, accurately, and with minimal experimental costs to determine the unknown parameters of the mathematical model of a chemical reaction. But, as is evident from the monographic specialized literature [1—3], chemical kinetics does not yet possess such methods. For this reason, mathematicians have recently [4—5] been intensively developing a number of effective methods as applied to the investigation of chemical kinetics. In the present article the possibility is shown of using the linearization method for calculating the order of chemical reactions.

1. Statement of the problem and method of its solution

Let us consider, for definiteness, a reaction proceeding without a change in volume. Then the variables of the reaction will be the concentrations \(x\) or the degre-

*) A. I. Ruban. Some questions of the mathematical description of dynamic objects. Dissertation. Tomsk State University, 1969.

of transformation. We assume that the reaction proceeds among \(n\) substances, some of which are reaction products and some the initial substances. Then the reaction rate in the general case can be written in the form [6]

\[ r=\kappa_1 \prod_{j=1}^{n} x_j^{\alpha_j}-\kappa_2 \prod_{j=1}^{n} x_j^{\beta_j}, \tag{1,1} \]

where \(\kappa_1,\ \kappa_2,\ \{\alpha_j\},\ \{\beta_j\}\) are unknown constants, with \(\alpha_j,\ \beta_j\) called the orders of the forward and reverse reactions with respect to substance \(A_j\). Since these constants are usually determined from the kinetic curves \(\{x_j(t)\}\) for \(t\in[0,T_1]\), it is necessary to formulate differential equations characterizing the changes of concentrations in time [6]

\[ \frac{dx_j}{dt}=m_j r,\quad j=1,\ldots,n,\quad x_j(0)=x_{j0}, \]

where \(m_j\) are stoichiometric coefficients. In general form this equation can be rewritten as

\[ \frac{d\bar{x}}{dt}=\bar{m}r(\bar{x},\bar{a}),\quad \bar{x}(0)=\bar{x}_0 . \tag{1,3} \]

Here \(\bar{x}(t)\) is an \(n\)-dimensional column vector of concentrations;

\(\bar{a}\) is an \(r\)-dimensional column vector of unknown constants, including all unknown constants \(\kappa_1,\ \kappa_2,\ \{\alpha_j\},\ \{\beta_j\}\);

\(\bar{m}\) is an \(n\)-dimensional column vector of stoichiometric coefficients;

\(r(\bar{x},\bar{a})\) is the reaction rate (1.1).

It is necessary to choose the coefficient vector \(\bar{a}\) so that

\[ \bar{\eta}(t)=\bar{q}(\bar{x},t) \tag{1,3} \]

is close to the quantities \(\eta^*(t)\) measured on the interval \(t\in[0,T_1]\). Here \(\bar{q}\) is an \(m\)-dimensional column vector of known transformations of \(x,t\), corresponding, as a rule, to particular measurement methods and measuring devices. To solve the posed problem we shall apply the method of linearization.

Denote by \(\bar{x}^{\,l}(t)\) the solution of the system of differential equations (1.2) under the initial conditions \(\bar{x}^{\,l}(0)=\bar{x}_0\) and with coefficients \(a^l\), obtained at the \(l\)-th iteration,

\[ \frac{d\bar{x}^{\,l}}{dt}=\bar{m}r(\bar{x}^{\,l},\bar{a}^{\,l}),\quad \bar{x}^{\,l}(0)=\bar{x}_0 . \tag{1,4} \]

This is a Cauchy problem, and its solution can be carried out on a digital computer by one of the standard numerical methods.

Since the quantity \(\bar{x}^{\,l+1}(t)\) satisfies an equation analogous to (1.4), for the concentration increments \(\Delta\bar{x}^{\,l+1}(t)=\bar{x}^{\,l+1}(t)-\bar{x}^{\,l}(t)\) we obtain the following equations:

\[ \frac{d\Delta\bar{x}^{\,l+1}}{dt} = \bar{m}\,[r(\bar{x}^{\,l+1},\bar{a}^{\,l+1})-r(\bar{x}^{\,l},\bar{a}^{\,l})], \quad \Delta\bar{x}^{\,l+1}(0)=0 . \tag{1,5} \]

We expand \(r(\bar{x}^{\,l+1},\bar{a}^{\,l+1})\) in a Taylor series at the point \((\bar{x}^{\,l},\bar{a}^{\,l})\), retaining only the linear terms and assuming the existence and continuity of the first derivatives. Then the differential equation (1.5) takes the form

\[ \frac{d\Delta\bar{x}^{\,l+1}}{dt} = J_1(\bar{x}^{\,l},\bar{a}^{\,l})\,\Delta\bar{x}^{\,l+1} + J_2(\bar{x}^{\,l},\bar{a}^{\,l})\,\Delta\bar{a}^{\,l+1}, \quad \Delta\bar{x}^{\,l+1}(0)=0 . \tag{1,6} \]

Here \(J_1^l(t)\), \(J_2^l(t)\) are the Jacobian matrices at the point \((\bar{x}^l(t), \bar{a}^l)\), which in our case degenerate into a row vector

\[ J_1^l(t)=\left(\frac{\partial r}{\partial x_j}\right)_{\bar{x}^l,\bar{a}^l},\quad J_2^l(t)=\left(\frac{\partial r}{\partial a_k}\right)_{\bar{x}^l,\bar{a}^l}; \tag{1,7} \]

\(\doteq\) denotes approximate equality;

\[ \Delta\bar{a}^{\,l+1}\doteq \bar{a}^{\,l+1}-\bar{a}^{\,l}. \tag{1,8} \]

It is now not difficult to obtain [7] the solution of the linear vector equation (1,6)

\[ \Delta\bar{x}^{\,l+1}(t)=\int_0^t \bar{X}^{\,l}(t)\bigl(\bar{X}^{\,l}(s)\bigr)^{-1}J_2^l(s)\,ds\,\Delta\bar{a}^{\,l+1}, \tag{1,9} \]

where \(\bar{X}^{\,l}(t)\) is the solution of the matrix equation

\[ \frac{d\bar{X}^{\,l}(t)}{dt}=J_1^l(t)\bar{X}^{\,l}(t),\quad \bar{X}^{\,l}(0)=I. \tag{1,10} \]

Here: \(I\) is the identity matrix, \(\bigl(\bar{X}^{\,l}(s)\bigr)^{-1}\) is the matrix inverse to \(\bar{X}^{\,l}(s)\).

It is not difficult to verify that the matrix standing before \(\Delta\bar{a}^{\,l+1}\) consists of the sensitivity functions [8]

\[ \bar{U}^{\,l}(t)=\left(\frac{\partial x_i}{\partial a_k}(t)\right)_{\bar{x}^l,\bar{a}^l}, \tag{1,11} \]

which, as is known [8], in the case under consideration are the solution of the system of ODEs

\[ \frac{d\bar{U}^{\,l}(t)}{dt}=J_1^l(t)\bar{U}^{\,l}(t)+J_2^l(t),\quad \bar{U}^{\,l}(0)=0. \tag{1,12} \]

At each iteration this system should be solved jointly with the original system of ODEs (1,4).

In order finally to linearize the original problem, when computing \(\eta(t)\) by (1,3) we also use the Newton–Raphson iterative procedure [5]

\[ \bar{\eta}^{\,l+1}\doteq \bar{\eta}^{\,l} +\left(\frac{\partial q_i}{\partial x_j}\right)_{\bar{x}^l,\bar{a}^l} \Delta\bar{x}^{\,l+1} =\bar{\eta}^{\,l}+\Gamma_{(t)}^l\Delta\bar{x}^{\,l+1}. \tag{1,13} \]

Substituting into this dependence the values \(\Delta\bar{x}^{\,l+1}\) from (1,9), (1,11), and seeking the minimum of the total squared deviations of the measured \(\bar{\eta}^{\,*}(t)\) and predicted \(\bar{\eta}(t)\) kinetic dependences,

\[ J(\bar{a})=\int_0^{T_1}(\bar{\eta}^{\,*}-\bar{\eta}^{\,l+1})^T \bar{Q}(t)(\bar{\eta}^{\,*}-\bar{\eta}^{\,l+1})\,dt =\min_{\Delta\bar{a}^{\,l+1}}, \tag{1,14} \]

that is,

\[ \Delta\bar{a}^{\,l+1} = \left(\int_0^{T_1}(\Gamma^l\bar{U}^{\,l})^T\bar{Q}\Gamma^l\bar{U}^{\,l}\,dt\right)^{-1} \int_0^{T_1}(\Gamma^l\bar{U}^{\,l})^T\bar{Q}(\bar{\eta}^{\,*}-\bar{\eta}^{\,l})\,dt. \tag{1,15} \]

Then

\[ \bar{a}^{\,l+1}=\bar{a}^{\,l}+\gamma^l\Delta\bar{a}^{\,l+1}, \tag{1,16} \]

where

\[ 0<\gamma^l\leqslant 1. \]

Here the following notation has been used:

\[ \bar{Q}(t) \]

is the matrix of known weight functions,

\[ T \]

is the symbol for matrix transposition.

If the value \(\bar{\eta}^{*}(t)\) is measured at discrete instants of time, then in the preceding formulas (1.14), (1.15) the integration should be replaced by summation over the corresponding discrete values. In this case the system of differential equations (1.4), (1.12) should be integrated in such a way as to obtain the exact solution at the prescribed instants of time.

2. Example. Diphenylchloromethane reacts with ethyl alcohol with the formation of hydrochloric acid and ether. It is assumed that the reaction proceeds according to the scheme

\[ \mathrm{RCl}+\mathrm{R'OH}\underset{k_2}{\stackrel{k_1}{\rightleftarrows}}\mathrm{HCl}+\mathrm{ROR'} \tag{2.1} \]

and experimental data are available for the variable HCl [1]:

\[ \begin{gathered} t_i\;(\text{min})=13;\quad 119;\quad 142;\quad 162;\quad 182;\quad 212;\\ i=1;\quad 2;\quad 3;\quad 4;\quad 5;\quad 6;\\ x_{\mathrm{HCl}}^{*}(t_i)\;(\text{mol})=0.00346;\;0.0268;\;0.0309;\;0.0343;\;0.0375;\;0.0418. \end{gathered} \tag{2.2} \]

Here the amount of alcohol was so large that in the course of the reaction it remained practically constant. The initial amount of diphenylchloromethane was \(\Theta=0.09966\) mole.

If by \(x(t)\) we denote the number of moles of HCl at time \(t\), then the differential equation describing the reaction takes the form

\[ \frac{dx(t)}{dt}=k_1(\Theta-x(t))-k_2x^2(t),\quad x(0)=0, \tag{2.3} \]

where \(k_1\) and \(k_2\) are unknown constants. To compute them we apply the linearization method described above. In this case

\[ \bar{\eta}^{*}(t)=x^{*}(t),\quad \eta(t)=x(t) \tag{2.4} \]

and at the \((l+1)\)-st iteration the parameter increments \(\Delta k_1^{\,l+1}\), \(\Delta k_2^{\,l+1}\), in accordance with criterion (1.14),

\[ J=\sum_{i=1}^{6}\left[ x^{*}(t_i)-x^{l}(t_i) -\frac{\partial x^{l}(t_i)}{\partial k_1}\Delta k_1^{\,l+1} -\frac{\partial x^{l}(t_i)}{\partial k_2}\Delta k_2^{\,l+1} \right]^2 = \min_{\Delta k_1^{\,l+1},\,\Delta k_2^{\,l+1}} \tag{2.5} \]

satisfy the system of linear nonhomogeneous algebraic equations

\[ \sum_{i=1}^{6}\left(\frac{\partial x^{l}(t_i)}{\partial k_1}\right)^2 \Delta k_1^{\,l+1} + \sum_{i=1}^{6} \frac{\partial x^{l}(t_i)}{\partial k_2} \frac{\partial x^{l}(t_i)}{\partial k_1} \Delta k_2^{\,l+1} = \]

\[ = \sum_{i=1}^{6}\bigl(x^{*}(t_i)-x^{l}(t_i)\bigr) \frac{\partial x^{l}(t_i)}{\partial k_1}, \]

\[ \sum_{i=1}^{6} \frac{\partial x^{l}(t_i)}{\partial k_1} \frac{\partial x^{l}(t_i)}{\partial k_2} \Delta k_1^{\,l+1} + \sum_{i=1}^{6}\left(\frac{\partial x^{l}(t_i)}{\partial k_2}\right)^2 \Delta k_2^{\,l+1} = \]

\[ = \sum_{i=1}^{6}\bigl(x^{*}(t_i)-x^{l}(t_i)\bigr) \frac{\partial x^{l}(t_i)}{\partial k_2}, \tag{2.6} \]

where the variables \(x^{l}(t_i)\),

\[ \frac{\partial x^{l}(t_i)}{\partial k_1}=u_1^{l}(t_i),\quad \frac{\partial x^{l}(t_i)}{\partial k_2}=u_2^{l}(t_i) \]

are the solution of the following system of differential equations:

Table 1

\[ \frac{d x^l(t)}{d t}=\kappa_1^l\bigl(\Theta-x^l(t)\bigr)-\kappa_2^l\bigl(x^l(t)\bigr)^2,\quad x^l(0)=0, \]

\[ \frac{d u_1^l(t)}{d t}=\Theta-x^l(t)-\kappa_1^l u_1^l(t)-2\kappa_2^l x^l(t)u_1^l(t),\quad u_1^l(0)=0, \tag{2.7} \]

\[ \frac{d u_2^l(t)}{d t}=-\bigl(x^l(t)\bigr)^2-\kappa_1^l u_2^l(t)-2\kappa_2^l x^l(t)u_2^l(t),\quad u_2^l(0)=0. \]

The algorithm described was implemented on the M-20 digital computer, and the results of adjusting the parameters \(\kappa_1, \kappa_2\) for various initial approximations (which are computed by some approximate method) are given in Table 1. The quantity \(J^l=\sum_{i=1}^{6}|x^*(t_i)-x^l(t_i)|^2\) characterizes the degree of agreement between model (2.3) and the real processes.

The simplicity of implementation, the high rate of convergence, and the possibility of solving identification problems under incomplete observability of the variables of the processes under study make it possible to recommend the linearization method as an effective method for calculating kinetic constants of complex chemical reactions.

In conclusion, the author considers it his duty to thank the students V. Doronin and M. Krasheninnikov, who implemented the linearization algorithm for the example described above on the M-20 digital computer.

REFERENCES

1. L. M. Batchner, M. E. Pozin. Mathematical Methods in Chemical Engineering. Khimiya Publishing House, Moscow, 1968.

2. N. M. Emanuel, D. G. Knorre. A Course in Chemical Kinetics (Homogeneous Reactions). Vysshaya Shkola, Moscow, 1969.

3. V. A. Lutsenko, L. N. Finyakin. Analog Computing Machines in Chemistry and Chemical Technology. Khimiya Publishing House, Moscow, 1969.

4. Applications of Computational Mathematics in Chemical and Physical Kinetics. Nauka, Moscow, 1969.

5. H. Rosenbrock, C. Storey. Computational Methods for Chemical Engineers. Mir, Moscow, 1968.

6. R. Aris. Analysis of Processes in Chemical Reactors. Khimiya, Leningrad, 1967.

7. R. Bellman. Introduction to Matrix Theory. Nauka, Moscow, 1969.

8. Sensitivity of Automatic Systems. Nauka, Moscow, 1968.