Page 213 - Numerical methods for chemical engineering

Differential-algebraic equation (DAE) systems 199

Example. Dynamics on the 2-D unit circle

We demonstrate using ode15s to solve a DAE-IVP for a simple system in which an algebraic
equation constrains a 2-D trajectory to the unit circle,

dx 1
= f 1 =−[x 1 − cos θ end ] + [x 2 − sin θ end ]
dt
2
2
0 = f 2 = x + x − 1 (4.217)
1 2
The trajectory is a simple relaxation to (cos θ end , sin θ end ) along the unit circle. The following
routine returns the function vector:
function f = calc f(t,x);
f = zeros(2,1);
*
theta end = 10/180 pi;
f(1) = -(x(1) - cos(theta end)) + (x(2) - sin(theta end));
f(2) = x(1)ˆ2 + x(2)ˆ2-1;
return;

[0]
We solve this DAE-IVP from an initial guess with x 1 = 0by
x 0 = [0; 0.8]; % set initial guess
M=[10;00]; % set mass matrix
Options = odeset(‘Mass’,M,‘MassSingular’,‘yes’, . . .
‘MState Dependence’,‘none’); % inform ode15s of mass matrix
[t traj,x traj] = ode15s(@calc f, [0 10], x 0, Options);
Here, we have a constant mass matrix, but we also could have provided a function that returns
M(t, x) were the mass matrix state-dependent (see help odeset for further details). Note
[0]
that for x 1 = 0, the two possible values of x 2 that satisfy the algebraic equation f 2 = 0
[0] [0]
are x =±1, but here the initial guess is x = 0.8. The resulting trajectory (Figure 4.14)
2 2
shows that the trajectory starts on the unit circle, satisfying f 2 = 0 and demonstrates that
ode15s first computes a consistent initial state.

Example. Heterogeneous catalysis in a packed bed reactor

We demonstrate the DAE-IVP use of ode15s for a more realistic chemical engineering
problem of modeling an isothermal packed bed reactor with a solid-catalyzed gas-phase
reaction A ⇔ B + C. First, A adsorbs reversibly to a vacant site S on the catalyst. The
adsorbed species A · S then reacts to form an adsorbed B · S and a free C molecule in the
gas phase. The adsorbed product B · S must then desorb, freeing again the active site.

A + S ⇔ A · S ˆ r aA = k aA p A c v − k dA c A·S

c B·S p C
A · S ⇔ B · S + C ˆ r s = k s c A·S − (4.218)
K s
B · S ⇔ B + S ˆ r dB = k dB c B·S − k aB p B c v