# Introduction to Matlab - PowerPoint by tyndale

VIEWS: 36 PAGES: 16

• pg 1
```									 Differential-Algebraic Equation Systems
1. DAE systems
2. Binary flash unit
3. Exothermic CSTR
4. Matlab example
DAE Systems
   General model structure
dy
 f ( x, y , z ) y ( x0 )  y 0 y  R n
dx
0  g( x, y, z ) 0  g( x0 , y 0 , z 0 ) x  R m

» Differential-algebraic equation (DAE) system
» Comprised of differential (x) and algebraic variables (z)
» Initial condition for algebraic variables must be computed
   Numerical solution
» Considerably more difficult than ODE systems due to
presence of algebraic equations (AEs)
» Need specialized techniques
» Difficulty determined by DAE system index
DAE System Index
   Definition
» The minimum number of times the algebraic equations
must differentiated to recover the algebraic variables
» Most codes developed for index 0 & index 1 systems
   Index 0 DAE systems
0  g( x, y, z )  z  Γ( x, y )

dy
 f ( x, y, z )  f [ x, y, Γ( x, y )]  F( x, y )
dx
Index 2 DAE Systems

g g dy g dz
0  g ( x, y , z )  0         
x y dx z dx
1                             1
dz     g   g g dy      g              g g            
     x y dx    z 
                        f ( x, y , z )   F ( x, y , z )
 x y            
dx     z                                                  

dy
 f ( x, y , z )
dx
Implicit Function Theorem
   Problem
z h ( x , y )
m  n 1
0  g( x, y, z) g : R              R   m
             0  g[ x, y, h( x, y )]

   Assumptions
 g       
g (a, b, c)  0 rank  (a, b, c)  m
 z       

   Result – there is a neighborhood of a & b where:
 h ( a, b )  c
h( x, y )  
 g[ x, y, h( x, y )]  0
   Implication – the function g(x,y,z) can be solved
numerically for z
Binary Flash Unit
           Schematic diagram                             Assumptions
»   Saturated liquid feed
V, y
Vapor          »   Perfect mixing
»   Negligible heat losses
»   Negligible vapor holdup
F, xF                                                  »   Negligible energy
Feed                                                                  accumulation
Q                  »   Constant heat of
H                         Heat
vaporization
»   Constant relative
L, x
Liquid                      volatility

           Vapor-liquid equilibrium
x
y
1  (  1) x
Model Formulation
   Mass balance
dH
 F V  L
dt
   Component balance

d ( Hx)
 FxF  Vy  Lx
dt

d ( Hx)    dH    dx                     dx
x    H     x( F  V  L)  H
dt      dt    dt                     dt

dx
H     ( F  V  L) x  FxF  Vy  Lx  F ( xF  x)  V ( y  x)
dt
Model Formulation cont.
Q
VH v  Q  V 
H v
   Index 0 DAE model

dx                  Q
H     F ( xF  x)        ( y  x)   dx
dt                 H                    f ( x, y )
x                 dt
0 y                             0  g ( x, y )
1  (  1) x
Exothermic CSTR

   Reaction: A  B
   Eliminate assumption of
constant cooling jacket
temperature
   Approximate cooling jacket
as well-mixed system

    Mass and energy balances on tank

dC A
V       q(C Ai  C A )  Vk0 exp( E / RT )C A
dt
dT
VCp     qC p (Ti  T )  (H )Vk0 exp( E / RT )C A  UA(Tc  T )
dt
Exothermic CSTR cont.
   Energy balance on cooling jacket
dTc
cVcC pc       c qcC pc (Tci  Tc )  UA(Tc  T )
dt

   Assume small thermal capacitance
dTc
cVcC pc        0  c qcC pc (Tci  Tc )  UA(Tc  T )
dt
   Index 0 DAE model
dC A
V       q(C Ai  C A )  Vk0 exp( E / RT )c A
dt
dT
VCp     qC p (Ti  T )  (H )Vk0 exp( E / RT )c A  UA(Tc  T )
dt
0   c qcC pc (Tci  Tc )  UA(Tc  T )
Matlab Mass Matrix
   Mass matrix representation
dy
M ( x, y )     f ( x, y )
dx

   CSTR example
1 0 0
M  0 1 0 
      
0 0 0 
      
 dy1 
 dx 
 dy   1
f ( x, y1 , y2 , y3 ) 
( x, y )  f ( x, y )   2    f 2 ( x, y1 , y2 , y3 )
dy
M
dx                          dx                           
 0   f 3 ( x, y1 , y2 , y3 ) 
                        

     

M-file cstr_dae.m

function f = cstr_dae(y,qc)

q = 1;           %   m^3/h
Cai = 10;        %   kmol/m^3
Ti = 298;        %   K
ko = 3.493e7;    %   h^-1
E = 11843;       %   kcal/kmol
H = -5960;       %   kcal/kmol
rhoCp = 500;     %   kcal/m^3/K
UA = 150;        %   kcal/H/K
R = 1.987;       %   kcal/kmol/K
V = 1;           %   m^3
Tci = Ti;

Ca = y(1);   T   = y(2);    Tc = y(3);

f = [
q/V*(Cai - Ca) - ko*exp(-E/R/T)*Ca;
q/V*(Ti-T) - H/rhoCp*ko*exp(-E/R/T)*Ca + UA/V/rhoCp*(Tc-T);
rhoCp*qc*(Tci-Tc)-UA*(Tc-T)
];
CSTR Model Solution
>> qc = 1;
>> f = @(y,qc) cstr_dae(y,qc);
>> yss = fsolve(f,[9 310 300],[],qc)
yss = 8.2931 314.5318 301.8150
 Integrate model equations from steady state

>> M = [1 0 0; 0 1 0; 0 0 0];
>> df = @(t,y,qc) cstr_dae(y,qc);
>> [t1,y1]=ode15s(df,[0,50],yss,[],qc)
>> ax1=plotyy(t1,y1(:,1),t1,y1(:,2));
>> xlabel('Time [h]');
>> ylabel(ax1(1),'Concentration [kmol/m^{3}]');
>> ylabel(ax1(2),'Temperature [K]');
Matlab Results for qc = 1

8.2931                                                           314.5318
Concentration [kmol/m 3]

Temperature [K]
8.2931                                                           314.5318

8.2931                                                           314.5318
0   5    10   15   20      25      30   35   40   45   50
Time [h]
CSTR Model Solution cont.
   Integrate model equations for different qc
>> qc=0.5;
>> [t2,y2]=ode15s(df,[0,50],yss,[],qc)
>> ax2=plotyy(t2,y2(:,1),t2,y2(:,2));
>> xlabel('Time [h]');
>> ylabel(ax2(1),'Concentration [kmol/m^{3}]');
>> ylabel(ax2(2),'Temperature [K]');
Matlab Results for qc = 0.5
10                                                           400

8                                                            380
Concentration [kmol/m 3]

Temperature [K]
6                                                            360

4                                                            340

2                                                            320

0                                                            300
0   5    10   15   20      25      30   35   40   45   50
Time [h]

```
To top