# How to numerically solve a reaction network with fast equilibrium steps?

For example If I have to solve the following reaction numerically

\begin{align} *\ce{ + A(g) &<=>[K1] [A*]} \tag{1}\\ \ce{[A{*}] &<=>[k1][k-1] [B{*}]} \tag{2}\\ \ce{[B*] &<=>[K2] [B](g) +}* \tag{3} \end{align}

Where the first and third reactions are fast equilibrium steps.

[* ] represents vacnt site and [A*] and [B*] represent adsorbed species.

The reaction is carried out in a closed container, therefore the concentrations of A(g) and B(g) are not constant. Given initial A(g) and B(g) pressure and assuming zero coverage of reactants initially, how do I solve a system like this numerically.

• After you specify all the parameters, it's essentially a differential equation problem. I would suggest use matlab (ode45 or other solver). If it fails, use BZZmath or SUNDIALS solver for stiff problem like the one you mentioned. – mamun Sep 12 '17 at 4:09
• I understand that its a differential equation problem and have solved problems of this nature before. But the rapid equilibrium steps confuse me, I have solved problems which dont have rapid equilibration. I dont have a reverse and backward constant for the first and third reaction, only an equilibrium constant. – ace_101 Sep 12 '17 at 4:32

For A adsorption: $$r_1=0$$ because the reaction is in equilibrium.So, $$K_1=\frac{[A]}{P_A \times [*]}$$ From this equation you can relate $[A]$ to $[*]$.Similarly, using the last reaction you can relate $[B]$ to $[*]$. Your [A] and [B] are connected via reaction 2. $$\frac{d[A]}{dt}=-\frac{d[B]}{dt}=-k_1[A]+k_{-1}[B]$$ This is a first order differential equation, so you need only one boundary condition. You have three unknown and three boundary conditions $([A]_0=[B]_0=0$ and $[*]=1.00)$.

• I think that when A * is in rapid equilibrium with A and * you can not write the rate of change of A as you have written. This expression would be valid if A* had not been in equilibrium with A and *. – ace_101 Sep 12 '17 at 22:28
• Yes, you are right. There should be two solution regime for this kind of problem (I forgot about that). One is for short time scale and another one for long time scale. I'll suggest consulting William Deen's analysis of transport phenomena book (maybe it was discussed in 3rd or 4th chapter) – mamun Sep 13 '17 at 5:16
• See chapter 4 of this link (mat.unimi.it/users/scacchi/didattica_2017/biomat2/asy.pdf). This is exactly what you are looking for. – mamun Sep 13 '17 at 5:18
• Thanks this seems useful, I did not think that you could apply the solution matching approach to a problem of this nature. But they solve the problem analytically, what if I want to do this numerically what do I do? This was just an example system that I chose, the real system that I am interested in solving is rather complicated and definitely doesn't have an analytical solution. Thanks – ace_101 Sep 14 '17 at 3:17
• Practically, I have no idea; however, you can simplify your problem by taking few relevant perturbation order (say up to 3rd order), form your system of differential equation system, solve it using a ODE/DAE/PDE solver. In essence, before you form your system of DE just make sure you have a well defined set, i.e., your number of equations and number of boundary conditions are consistent. – mamun Sep 14 '17 at 5:32