$$\require{mhchem}$$
Computer exercise 2
Consider a simplified $\ce{O3}$-$\ce{NO_x}$-VOC system in which the VOC is represented by a single molecule, propene ($\ce{C3H6}$). Writing this as a generalized mechanism (propene written as $\ce{RH}$), reactions and rate constants are shown in the following table (Seinfeld and Pandis, 2006):
\begin{align} \newcommand{\conc}{\text{cm$^{3}$ molec$^{-1}$ s$^{-1}$}} & \text{Reaction} & \text{Rate constant (298 K)} \\ \hline 1\quad & \ce{RH}+\ce{OH} \stackrel{\ce{O2}}{\to} \ce{RO2}+\ce{H2O} & 26.3 \times 10^{-12} \conc\\ 2\quad & \ce{RO2}+\ce{NO} \stackrel{\ce{O2}}{\to} \ce{NO2}+\ce{R^'CHO}+\ce{HO2} & 7.7 \times 10^{-12} \conc\\ 3\quad & \ce{HO2}+\ce{NO} \to \ce{NO2}+\ce{OH} & 8.1 \times 10^{-12} \conc\\ 4\quad & \ce{OH}+\ce{NO2} \stackrel{\ce{M}}{\to} \ce{HNO3} & 1.1 \times 10^{-11} \conc \text{(at 1 atm)}\\ 5\quad & \ce{HO2}+\ce{HO2}\to \ce{H2O2}+\ce{O2} & 2.9 \times 10^{-12} \conc\\ 6\quad & \ce{RO2}+\ce{HO2}\to \ce{ROOH}+\ce{O2} & 5.2 \times 10^{-12} \conc\\ 7\quad & \ce{NO2}+h\nu \stackrel{\ce{O2}}{\to} \ce{NO}+\ce{O3} & \text{typical} \approx 0.015 \text{s$^{-1}$} \\ 8\quad & \ce{O3}+\ce{NO} \to \ce{NO2}+\ce{O2} & 1.9 \times 10^{-14} \conc \end{align}
We will define maximum ozone concentrations from a 10-hour simulation specified from intial values of $\ce{RH}$, $\ce{NO}$, and $\ce{NO2}$ (constant photolysis rate).
Problems:
Conditions:
Assume:
We can write the rate equation above with initial conditions in vector form:
$$ \frac{d \mathbf{x}}{dt} = \boldsymbol{\nu} \mathbf{r} = \boldsymbol{f}\left(\mathbf{x}, \mathbf{k}, \boldsymbol{\nu}\right), \quad \mathbf{x}(t=0) = \mathbf{x}_0 \\ $$
Essentially, the rates of change in concentrations are a functions $\boldsymbol{f}$ of concentrations, rate coefficients $\mathbf{k}$, and stoichiometric numbers $\boldsymbol{\nu}$ relevant for each reaction.
The PSSA (pseudo steady-state assumption) (also referred to as the quasi steady-state assumption, or QSSA) assumes the existence of a set of slow and fast species $\mathbf{x} = \{\mathbf{x}^{(1)}, \mathbf{x}^{(2)}\}$; the dynamics for the fast species drive them toward their equilibrium distribution while the concentration of the slow species vary little over the same timescale. Mathematically, PSSA leads to the assumption that the dynamics of these species can be decoupled and solved by different methods to obtain the same chemical trajectories obtained by solving the full set of equations:
\begin{align} \frac{d \mathbf{x}^{(1)}}{dt} &= \boldsymbol{f}^{(1)}\left(\mathbf{x}, \mathbf{k}, \boldsymbol{\nu}\right), \quad \mathbf{x}^{(1)}(t=0) = \mathbf{x}^{(1)}_0 \\ \frac{d \mathbf{x}^{(2)}}{dt} &= \boldsymbol{f}^{(2)}\left(\mathbf{x}, \mathbf{k}, \boldsymbol{\nu}\right) \approx \boldsymbol{0} \\ \end{align}
Clarifications about PSSA:
Turányi and Tomlin (2014) can be consulted for more information (with historical perspective), if desired.
For this chemical system, $\mathbf{x}^{(1)} = \{\ce{[RH]}, \ce{[NO]}, \ce{[NO2]}, \ce{[O3]}, \ldots\}$, $\mathbf{x}^{(2)} = \{\ce{[OH]}, \ce{[HO2]}, \ce{[RO2]}\}$. The problem statement determines which specific reactions to consider for the steady-state expressions for each of the species in $\mathbf{x}^{(2)}$. $\{\ldots\}$ for $\mathbf{x}^{(1)}$ includes other products which are formed in the reactions that can be tracked but can also be neglected as they do not participate in further reactions leading to ozone production (for the purpose of this exercise). To obtain solutions:
Note that the function of rate equations you will want to integrate with the numerical solver should accept a vector of concentrations, $\mathbf{x}^{(1)}$ and return a vector of rates $d\mathbf{x}^{(1)}/dt$. Concentrations for $\mathbf{x}^{(2)}$ should also be calculated (through the PSSA) within the same function as they are necessary for calculating $d\mathbf{x}^{(1)}/dt$ at every time step.
The numerical solver recommended for GNU Octave is lsode
and for MATLAB is ode45
or ode15s
(for further information, see Beers, 2007).
Good practice:
For instance, functions you will want to have:
Readability:
It is often the case that you have a vector of concentrations. To access individual elements of the vector, it may be convenient to pass around a vector of indices so the species referred to is clear. For example, you can define a structure
like the following in your main script to hold values of vector indices:
{octave}
ind = struct('RH', 1, 'NO', 2, 'NO2', 3, 'O3', 4);
If x
is vector of concentrations at a given time step, you can access the $\ce{NO_2}$ concentrations using the syntax, x(ind.NO2)
.
Rate of change in concentrations. Here is an example of the function to be integrated - function #3 described above. You fill in the {...}
:
{octave}
function xdot = rateeqns(x, t, k, P_HOx, ind)
%% x is concentration
%% t is time
%% k are rate constants
%% P_HOx is HOx production rate
%% ind is the structure of indices
%% PSSA
[OH, HO2, RO2] = pssasoln(x, k, P_HOx, ind);
%% rates
r = zeros(size(k));
r(1) = k(1)*x(ind.RH)*OH;
r(2) = {...}
r(3) = {...}
r(4) = {...}
r(5) = {...}
r(6) = {...}
r(7) = {...}
r(8) = k(8)*x(ind.O3)*x(ind.NO);
%% output
xdot = zeros(size(x)); % {RH, NO, NO2, O3}
xdot(ind.RH) = -r(1);
xdot(ind.NO) = {...}
xdot(ind.NO2) = {...}
xdot(ind.O3) = r(7) - r(8);
end
Plotting contours. For Problem 4a, you should obtain a matrix which contains maximum ozone concentrations for each initial value of $[\ce{RH}]_0$ and $[\ce{NO_x}]_0$. To plot a contour plot, you can find some help on these pages: 1 2 rather than in the help documentation. You will likely have a vector of initial RH and a vector of initial NOx concentrations. If the vectors are called init_RH
and init_NOx
with ozone matrix max_O3
, you will call meshgrid
prior to contour
as shown in the following example:
{octave}
[grid_RH, grid_NOx] = meshgrid(init_RH, init_NOx);
contour(grid_RH, grid_NOx, max_O3)
scale(gca, 'xscale', 'log', 'yscale', 'log');
Beers, K. (2007): Numerical Methods for Chemical Engineering: Applications in MATLAB. Cambridge University Press, New York, NY.
Seinfeld, J. H., and S. N. Pandis (2006): Atmospheric Chemistry and Physics: From Air Pollution to Climate Change. 2nd edition. John Wiley & Sons, New York, NY.
Turányi, T. and A. S. Tomlin (2014): Analysis of Kinetic Reaction Mechanisms. Springer, New York, NY.