(This question comes from a Chemical Engineering background, I hope it still falls in the scope of engineering.stackexchange.com If not, please move.)

I am trying to calculate and plot fractional yield as a function of temperature for a parallel reaction of first-order reactions forming 3 products R, S and P from reactant A.

enter image description here

I know the highest possible reaction temperature is 500k. S is the desirable product and I am trying to find the reaction temperature to maximize its production and calculate reactor outlet concentration of S for XA = 100%. Thus far, I have tried to construct the respective rate equations and thus corresponding Arrhenius equations (which I have data for), what I'm not sure is if I am plotting this equation against temperature, is this a valid solution to evaluate the fractional yield?

Here is my MATLAB code and corresponding graph:

CA0 = 100; % Initial concentration of A (mol/m^3)

V_reactor = 1; % Reactor volume (for example, 1 m^3)

T_range = linspace(0, 1000, 10); % Temperature range from 0 K to 500 K

YieldS = zeros(size(T_range));

for i = 1:length(T_range)
    T = T_range(i); % Temperature at this iteration
    R = 8.314; % Universal gas constant (J/(mol*K))

k1 = 5*10^8 * exp(-65000 / (R * T)); % Adjust activation energy (kJ/mol) as needed
k2 = 10^9 * exp(-50000 / (R * T)); % Adjust activation energy (kJ/mol) as needed
k3 = 10^7 * exp(-45000 / (R * T)); % Adjust activation energy (kJ/mol) as needed

dC = @(V, C) [-k1 * C(1) - k2 * C(1) - k3 * C(1); % dCA/dV
               k1 * C(1); % dCR/dV
               k2 * C(1); % dCS/dV
               k3 * C(1)]; % dCP/dV

% Initial conditions (at V = 0)
C0 = [CA0; 0; 0; 0]; % Initial concentrations of A, R, S, and P

[V, C] = ode45(dC, [0, V_reactor], C0);

CA_outlet = C(end, 1);
CS_outlet = C(end, 3);

YieldS(i) = (CS_outlet / CA0) * 100;

enter image description here

  • $\begingroup$ I don't fully understand your code but the result does not seem 'physical'. Is not the yield of S just $k_2/(k_1+k_2+k_3)$ where the $k$'s depend on temperature via the Arrhenius eqn. $\endgroup$
    – porphyrin
    Commented Oct 9, 2023 at 6:55

1 Answer 1


I would like to point three things.

1. Mistake in mole balance As we frequently do in chemical engineering, we write the balances too fast and make some mistakes. In the following equations $F_j$ denotes the molar flow rate of species $j$ in $\pu{mol s^-1}$

The PFR mole balance, for a general species $j$ in the case for $k$ chemical reactions, yields \begin{align} \frac{\mathrm{d}F_j}{\mathrm{d}V} &= \sum_k \nu_{j,k}r_{k} \\ \frac{\mathrm{d}(v_0C_j)}{\mathrm{d}V} &= \sum_k \nu_{j,k}r_{k} \rightarrow \boxed{\frac{\mathrm{d}C_j}{\mathrm{d}V} = \frac{1}{v_0}\sum_k \nu_{j,k}r_{k}} \tag{1} \end{align} where $v_0$ is the volumetric flow rate in $\pu{m^3 s^-1}$. You forgot this guy in the mole balance. It is nice to remember that we can only take out $v_0$ from the differential when it is constant. For isothermal operations in liquid-phase reactions, this is an excellent approximation, because we assume that the fluid is incompressible. Make sure the species are liquids.

The four differential equations are in fact \begin{align} \frac{\mathrm{d}C_A}{\mathrm{d}V} &= -\frac{1}{v_0}(k_1C_A + k_2 C_A + k_3 C_A) \tag{2} \\ \frac{\mathrm{d}C_R}{\mathrm{d}V} &= \frac{1}{v_0}k_1C_A \tag{3} \\ \frac{\mathrm{d}C_S}{\mathrm{d}V} &= \frac{1}{v_0}k_2C_A \tag{4} \\ \frac{\mathrm{d}C_P}{\mathrm{d}V} &= \frac{1}{v_0}k_3C_A \tag{5} \\ \end{align}

2. Definition of fractional yield Chemical engineering literature is not homogeneous in the definition of this quantity. For your case you are calculating the overall fractional yield, but you have to change the denominator. The overall fractional yield of species $j$ with respect to reactant $A$, denoted by $\phi_{j,A}$, is \begin{align} \require{cancel} \phi_{j,A} &\equiv \frac{F_{j} - \cancel{F_{j0}}}{F_{A0} - F_{A}} \rightarrow \phi_{j,A} = \frac{\cancel{v}C_{j}}{\cancel{v_0}C_{A0} - \cancel{v}C_{A}} \rightarrow \boxed{\phi_{j,A} = \frac{C_{j}}{C_{A0} - C_{A}}} \tag{6} \end{align} for the case that: (1) no species $j$ enters the reactor, i.e. $F_{j0} = 0$; and (2) the fluid is incompressible, i.e. the volumetric flow rate is constant $v_0 = v$.

In your code you are only dividing by $C_{A0}$.

3. Solving If you are learning how to use MATLAB that is perfectly fine. However, Eqs. (2-5) are a system of ordinary linear differential equations, which are perfectly solved by pencil and paper. With some effort try to find that \begin{align} C_A(V) &= \exp\left[-\left(\frac{k_1 + k_2 + k_3}{v_0}\right)V\right] C_{A0} \tag{7} \\ C_j(V) &=\frac{k_j C_{A0}}{k_1 + k_2 + k_3} \left\{1 - \exp\left[-\left(\frac{k_1 + k_2 + k_3}{v_0}\right)V\right]\right\} \quad j = R,S,P \tag{8} \end{align} Try to manipulate Eqs. (7) and (8) to find what you want. Retaining analytical expressions may give you much more insight of the system, because you can obtain clear dependences of the variables. When integrating numerically, you lose track of this and, you resort to changing the parameters to see how the systems evolves differently along time.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.