Flow Reactor

In this section we solve the equations for a flow reactor, first for an isothermal case and then for a non-isothermal case. In many flow reactors the effect of radial and axial dispersion is neglibible, and we neglect them here; they are treated elsewhere. (link)

The first problem is very simple, and is patterned after a problem on the California Professional Engineers License Examination, according to Fogler (ref). Here we modify it. We take a reactor in which components A and C are fed in equimolar amounts, and the following reaction takes place:

2 A -> B

The equations are (ref)

where F_j is the molar flow rate of the j-th species, V_R is the reactor volume (cumulative from the inlet), and r_j is the rate of reaction of the j-th species, in moles/volume time. Here the rate of reaction is taken as second order.

where the units on k are volume/mole time. The equations for all three species are then

At the inlet we take

Solve for the case of k = 0.3 m³/kmol s, and integrate until V_R = 1.2 m³.

The equations look a bit funny. There is a derivative of the molar flow rate, but the reaction rate depends on the concentration. It is possible to derive the concentration of each species from the molar flow rates, and this section illustrates how to handle this complication in Matlab. The calculation of concentration depends on some assumptions, and here we choose the case in which the total concentration remains constant. This is appropriate to liquids and dilute gas streams at constant temperature and pressure. Then the concentration of a species is the mole fraction times the total concentration, e.g.

The total concentration is the total concentration at the inlet, which is taken as 5 kmol/m³. Thus, the flow rate is

The problem is summarized here.

Keep in mind that the reactor volume is the cross section of the reactor times the length. Thus, we can write

and the corresponding equations would be

Here we solve in terms of the accumulated reactor volume, V_R. We recognize that this problem can be solved analytically, but we do it numerically.

The Matlab program requires that we write a function which defines the right-hand side, and allows the right-hand side to be calculated, given the input, which are the flow rates of all species. Thus, to solve the problem we use the variables

The function will have available to it the accumulated reactor volume V_R, and the local value of y_j, j=1, 2, 3 at the same reactor volume. We then calculate the concentrations at that reactor volume

using these equations.

The rates of reaction are then evaluated, and the function returns the numerical value of the right-hand side. The code for the function is

% rhs1.m
% This function gives the right-hand side for a simple
% reactor problem which is isothermal.
function ydot=rhs1(z,y)
rate = 0.3*y(1)*y(1);
ydot(1) = - 2.*rate;
ydot(2) = + rate;
ydot(3) = 0.;
ydot = ydot';

We test this by calling it with specific values for y_j, j=1, 2, 3 to insure that it is correct. Using

we get

which agrees with hand calculations. This is a most important step, because this is where we are adding value. Matlab will integrate whatever equations we give it, right or wrong, and only we can insure that we've solved the right equations.

Next we write the code as the driver. This code must set the constants (note that one of them was just put into the function ydot for simplicity), set the initial conditions and total reactor volume, and then call the ode solver.

% run_rhs1.m
% This is the driver program to solve the simple
% flow reactor example.

% set the initial conditions
y0 = [2 0 2]

% set the total volume of the reactor
vrspan = [0 1.2]

% call the ode solver
[vr y] = ode45('rhs1',vrspan,y0)

% plot the result
plot(vr,y)
xlabel('Reactor Volume')
ylabel('Molar Flow Rates')
legend('A', 'B', 'C')

The result is shown in the figure. Now we can vary the parameters and see how the solution changes with changes in the parameters.

A more complicated, non-isothermal case is also treated.

Take Home Message: For reactor problems, we need to write a function that will evaluate the right-hand side of the set of differential equations, given the solution at that point. This function defines the problem we are solving, and we must check it carefully. The ode solvers in Matlab do the rest.