Difference between revisions of "Software for integrating ODEs"

From Process Control: 3P4
Jump to navigation Jump to search
(Created page with "== Background == Numerically integrating ODE's is something you should have mastered in your pre-requisite 3E4 course. If you are not familiar with this topic, it is your res...")
 
Line 5: Line 5:
Here is a tutorial a wrote a few years ago when I [http://modelling3e4.connectmv.com/wiki/Software_tutorial/Integration_of_ODEs taught 3E4 in 2010]. The tutorial below is similar, but uses a reactor design example.
Here is a tutorial a wrote a few years ago when I [http://modelling3e4.connectmv.com/wiki/Software_tutorial/Integration_of_ODEs taught 3E4 in 2010]. The tutorial below is similar, but uses a reactor design example.


== Example ==
The example we will work with is a common modelling reaction: a liquid-based stirred tank reactor, with (initially) constant physical properties, a second order chemical reaction where species A is converted to B according to  \({\sf A} \rightarrow {\sf B} \), with reaction rate \(r = k C_{\sf A}^2\).  One can find the time-dependent mass balance for this system to be:


Consider the pair of differential equations which you have or will cover in your Reactor Design course:
\[ \frac{dC_{\sf A}(t)}{dt} = \frac{F(t)}{V} \bigg( C_{{\sf A},\text{in}} - C_{\sf A} \bigg) - k C_{\sf A}^2 \]
<rst>
<rst-options: 'toc' = False/>
<rst-options: 'reset-figures' = False/>
.. math::


\dfrac{dX}{dW} &= \dfrac{-r_A'}{F_{A0}} \\ \dfrac{dy}{dW} &= -\dfrac{\alpha}{2y}\left(1 + \varepsilon X\right)
where \(C_{{\sf A},\text{in}} = 5.5\) mol/L, we will initially assume constant volume of \(V = 100\) L and constant inlet flow of \(F(t) = 20.1 \) L/min.  The reaction rate constant is 0.15 \(\frac{\text{L}}{\text{mol}.\text{min}}\).  We must specify an initial condition for every differential equation: we will let \( C_{\sf A}(t=0) = 0.5\) mol/L.
Some terminology (recap from your pre-requisite math courses)


* The independent variable is \(W\)
In the code below we will integrate the ODE from \(t_\text{start} = 0.0\) minutes up to \(t_\text{final} = 10.0\) minutes and plot the time-varying trajectory of concentration in the tank.
* The two dependent variables (variables being integrated with respect to the independent variable) are \(X\) and \(y\)
Since there are two dependent variables, we require initial conditions for each variable. In this case:
 
* :math:`X(0) = 0.0`
* :math:`y(0) = 1.0`
We also need to specify initial **and final conditions** for the independent variable:
 
* :math:`W` initially is 0.0 at the reactor entrance
* :math:`W` finally is 28.0 kg at the reactor exit
But there are a few other unknowns we must first specify:
 
* :math:`r_A' = -k' \dfrac{(1-X)y}{(1+\varepsilon X)}`
* :math:`q = \dfrac{(1 + \varepsilon X)}{y}`
* :math:`F_{A0} = 0.1362\,\text{mol.s}^{-1}`
* :math:`k' = 0.0074\,\text{mol.s}^{-1}\text{.(kg catalyst)}^{-1}`
* :math:`\alpha = 0.0367\,\text{kg}^{-1}`
* :math:`\varepsilon = -0.15`
</rst>
 
With MATLAB we create two files: a file containing the model equation, and another that drives the ODE integrator. Once you have both files on your computer, call <tt>ODEdriver</tt> from the MATLAB command line.
 
=== Model file: <tt>pbr.m</tt>===
 
<syntaxhighlight lang="MATLAB">
function d_depnt__d_indep = pbr(indep, depnt)
% Dynamic balance for the packed bed reactor (PBR); demo problem class 05C
%
%    indep: the independent ODE variable, such as time or length
%    depnt: a vector of dependent variables
%
%    X = depnt(1) = the conversion
%    y = depnt(2) = the pressure ratio = P/P_0 = y
%
%    Returns d(depnt)/d(indep) = a vector of ODEs
% Assign some variables for convenience of notation
X = depnt(1);
y = depnt(2);
 
% Constant(s)
FA0  = 0.1362;  % [mol/s]
kdash = 0.0074;  % [mol/(kg catalyst . s)]
alpha = 0.0367;  % [1/kg]
eps  = -0.15;  % [-]
 
% Algebraic equation(s)
rAdash = -kdash * (1-X)/(1+eps*X) * y;
 
% Output from this ODE function must be a COLUMN vector, with n rows
n = numel(depnt);
d_depnt__d_indep = zeros(n,1);
d_depnt__d_indep(1) = -rAdash / FA0;
d_depnt__d_indep(2) = -alpha/(2*y) * (1+eps*X);
</syntaxhighlight>
 
=== ODE integrator: <tt>ODEdriver.m</tt> ===
<syntaxhighlight lang="MATLAB">
% Integrate the ODE
% -----------------
% The independent variable always requires an initial and final value:
indep_start = 0.0;  % kg
indep_final = 28.0;  % kg
% Set initial condition(s): for integrating variables (dependent variables)
X_depnt_zero = 0.0;  % i.e. X(W=0) = 0.0
y_depnt_zero = 1.0;  % i.e. y(W=0) = 1.0
% Integrate the ODE(s):
[indep, depnt] = ode45(@pbr, [indep_start, indep_final], [X_depnt_zero, y_depnt_zero]);
% Plot the results:
clf;
plot(indep, depnt(:,1), 'b')
grid('on')
hold('on')
plot(indep, depnt(:,2), 'g')
xlabel('Catalyst weight, W [kg]')
ylabel('X and y')
legend('X', 'y')
</syntaxhighlight>

Revision as of 16:23, 13 January 2014

Background

Numerically integrating ODE's is something you should have mastered in your pre-requisite 3E4 course. If you are not familiar with this topic, it is your responsibility to quickly catch up, because we will use this intensively for the rest of the course.

Here is a tutorial a wrote a few years ago when I taught 3E4 in 2010. The tutorial below is similar, but uses a reactor design example.

The example we will work with is a common modelling reaction: a liquid-based stirred tank reactor, with (initially) constant physical properties, a second order chemical reaction where species A is converted to B according to \({\sf A} \rightarrow {\sf B} \), with reaction rate \(r = k C_{\sf A}^2\). One can find the time-dependent mass balance for this system to be:

\[ \frac{dC_{\sf A}(t)}{dt} = \frac{F(t)}{V} \bigg( C_{{\sf A},\text{in}} - C_{\sf A} \bigg) - k C_{\sf A}^2 \]

where \(C_{{\sf A},\text{in}} = 5.5\) mol/L, we will initially assume constant volume of \(V = 100\) L and constant inlet flow of \(F(t) = 20.1 \) L/min. The reaction rate constant is 0.15 \(\frac{\text{L}}{\text{mol}.\text{min}}\). We must specify an initial condition for every differential equation: we will let \( C_{\sf A}(t=0) = 0.5\) mol/L.

In the code below we will integrate the ODE from \(t_\text{start} = 0.0\) minutes up to \(t_\text{final} = 10.0\) minutes and plot the time-varying trajectory of concentration in the tank.

Retrieved from "https://learnche.org/wiki_3P4/index.php?title=Software_for_integrating_ODEs&oldid=147"