Cinetica Chimica
Applicata
Applicazione N. 1
(D. Dalle Nogare, P. Canu; a.a. 2006-2007)
Dati del problema
La produzione di P da A e B appare avvenire secondo il
meccanismo:
k1
CD
A B
k
1
2C
k2
k
P
2
k3
C A
k
R
3
Sono date
le cinetiche kj+ e kj+
la composizione in ingresso CIN [mol/L]
la portata volumetrica V [L/s]
Richieste (i)
i profili di concentrazione di tutte le specie C(z)
lungo un reattore tubolare di diametro 5 cm,
l’evoluzione della concentrazione di tutte le
specie C(t) in un reattore batch,
la lunghezza ottimale del reattore tubolare per
avere:
la massima concentrazione di P in uscita,
la massima selettività a P (= P prodotto sul totale dei
prodotti)
la massima resa in P di A (= P prodotto per reagente A
consumato)
la massima conversione di B
Richieste (ii)
concentrazione di P, selettività a P, resa in P di
A e conversione di B (variabili usate ai punti
A..D precedenti) in uscita da un reattore CSTR
di volume pari al PFR dei casi A, B, C, D
precedenti e confrontare con i corrispondenti
valori dati dal PFR.
[facoltativo] comparare VPFR e VCSTR a cui
concentrazione di P, selettività a P, resa in P di
A e conversione di B (variabili usate ai punti
A..D precedenti) sono massime.
Riassunto della teoria
Reattore differenziale:
PFR
dC
z
Batch
r,
,t
d
v
Reattore CSTR
C C
IN
r
Velocità di produzione
r R
Significato della formula (I)
Velocità di produzione:
Esempio per C:
(vC,1=1, vC,2=-2, vC,3=-1)
rC=R1-2R2-R3
r R
k1
CD
A B
k
1
2C
R1 = k1 ·CA·CB - k1 ·CC·CD
R2 = k2+ ·CC2 - k2- ·CP
R3 = k3+ ·CC·CA - k3- ·CR
+
Specie con c.s.<0
-
k2
k
P
2
k3
C A
k
Specie con c.s.>0
3
R
Applicazione numerica (I)
Velocità di produzione:
function r = erre(C,nu,k)
r R
kd
ki
k Nr
[Nc,Nr] = size(nu);
nur = (nu-abs(nu))/2; % coeff. stech dei reagenti
nup = (nu+abs(nu))/2; % coeff. stech dei prodotti
Nc
Nr
for j = 1:Nr
Rd = k(j,1)*prod(C.^abs(nur(:,j))) ; % vel. di reazione diretta
Ri = k(j,2)*prod(C.^nup(:,j)) ;
% vel. di reazione inversa
R(j) = Rd-Ri;
end
r = nu*R';
% velocita' di produzione delle singole specie
Significato della formula (II)
Reattore differenziale:
dC
z PFR Batch
r,
,t
d
v
Noto Vdot si conosce v.
Noto CIN si calcola r ed è quindi nota la derivata
della composizione in τ=0: dC/dτ|0.
Si può integrare numericamente il sistema di
equazioni differenziali.
Applicazione numerica (II)
Reattore differenziale:
dC
z PFR Batch
r,
,t
d
v
% nel programma principale, si chiama l’integratore “ode45”
options=[];
[tau,C] = ode45(@BMidiff,[0 tauf],C0,options,nu,k);
%-----------------------------------------------------------------------% bilancio materiale DIFF
function Cprimo = BMidiff(tau,C,nu,k)
r = erre(C,nu,k);
% velocita' di produzione delle singole specie
Cprimo = r;
% BMi = derivate nel tempo delle concentrazioni
Significato della formula (III)
Reattore CSTR:
C C
IN
r
Noto Vdot si conosce θ.
Il bilancio è macroscopico: si possono calcolare
solo le composizioni in uscita.
Poiché r=r(C), si tratta di risolvere un sistema
di equazioni algebriche non lineari: richiede una
risoluzione numerica.
Applicazione numerica (III)
Reattore CSTR:
C C
IN
r
% nel programma principale,
si chiama il risolutore di equazioni
algebriche non lineari “fsolve”
options = optimset('Display','off');
Cguess = CIN./2; % stima della possibile soluzione
C=fsolve(@BMicstr,Cguess,options,CIN,thi,nu,k)
%----------------------------------------------------------------------% bilanci materiali CSTR
function err = BMicstr(C,CIN,th,nu,k)
r = erre(C,nu,k); % velocita' di produzione delle singole specie
DC = C-CIN;
% differenza OUT-IN
err = r*th - DC; % 0 = BMi
Struttura di un file matlab
Gli “script” di matlab, con estensione “.m”,
si strutturano in questo modo:
function applicazione1
% Cinetica Chimica Applicata – Appl. 1 – Programma principale
|||| Corpo del programma principale:
|||| vengono eseguiti tutti i comandi in succesione
|||| e tutte le chiamate alle altre funzioni.
|||| Solo alla fine del programma principale si posizionano le
|||| funzioni chiamate dal programma stesso
function Cprimo = BMidiff(tau,C,nu,k)
|||| Corpo della funzione
function err = BMicstr(C,CIN,th,nu,k)
|||| Corpo della funzione
function r = erre(C,nu,k)
|||| Corpo della funzione
Nota Bene
Fare
uso il più possibile dell’help di
Matlab, per capire come agiscono le
varie funzioni e soprattutto per
controllare la sintassi dei vari
comandi
Mettere “;” alla fine di ogni comando
per evitare che vengano stampati tutti i
risultati. Eliminarlo solo quando serve
monitorare un calcolo.
Algoritmo
Definire tutti i dati
Chiamare l’integratore ODE
Graficare i risultati (pto 1)
Rielaborare i risultati per ottenere il pto 2
Calcolare le proprietà richieste al pto 3
Calcolare i 4 CSTR relativi al pto 4
[facoltativo] Calcolare una serie di CSTR di
dimensioni crescenti, analizzare i risultati
come al pto 3 e confrontare
Algoritmo – dettagli (1)
Definire tutti i dati
function gruppo1
% Cinetica Chimica Applicata – Appl. 1 – Programma principale
% Definizione dati cinetici
nu=[ .. .. ..
.. .. ..];
k=[.. ..
.. ..];
% Definizione dati di ingresso
C0=[..; ..];
Vdot=..;
Algoritmo – dettagli (2)
Chiamare l’integratore ODE
% Scegliere un tempo di permanenza totale
tauf=..;
% Chiamare l’integratore di eq. differenziali “ode45”
(vedi “reattore differenziale”)
% Salvare i dati in una variabile dedicata
CPFR=C; tauPFR=tau;
N.B. A priori non si conosce il tempo di
permanenza totale necessario al completarsi
della reazione. Si fa un primo tentativo con
60s e si aumenta o diminuisce a seconda dei
risultati. In questo ci si aiuta con il grafico.
Algoritmo – dettagli (3)
Graficare i risultati (pto 1)
Poichè il testo richiede i profili di
concentrazione, il modo più intuitivo di darli è
attraverso un grafico. Attenzione: sono
richiesti i profili C(z), e non C(τ)!
% Grafico profili di concentrazione nel reattore tubolare
figure
plot(zPFR,CPFR)
xlabel(‘z [cm]’),ylabel(‘C [mol/L]’),axis tight
legend('A','B','C','D','P','R',-1)
title('PFR')
Algoritmo – dettagli (4)
Rielaborare i risultati per ottenere il pto 2
Anche in questo caso si tratta di graficare i
profili di concentrazione C(t).
% Grafico profili di concentrazione nel reattore Batch
figure
plot(t,C)
xlabel(‘t [s]’),ylabel(‘C [mol/L]’),axis tight
legend('A','B','C','D','P','R',-1)
title(‘Batch’)
Algoritmo – dettagli (5)
Calcolare le proprietà richieste al pto 3
Massime CP SP YPA XB
% Trovare il massimo in un vettore
% Definire il vettore
CP=CPFR(:,5); % SP=...; YPA=...; XB=...;
% Individuare il massimo e l’indice corrispondente nel vettore
% con la funzione “max” di matlab
[CPmax,ICPmax]=max(CP);
% Trovare il corrispondente valore in z
zCPmax=z(ICPmax);
% Ripetere per SP, YPA e XB, ottenendo SPmax e zSPmax,
% YPAmax e zYPAmax, XBmax e zXBmax
Algoritmo – dettagli (6)
Calcolare i 4 CSTR relativi al pto 4
% CSTR 1:
thi=tauPFR(ICPmax);
% Chiama il risolutore di eq. algebriche non lineari “fsolve”
(vedi reattore CSTR)
% Calcola CP,SP,YPA,XB per il CSTR
CP1cstr=...; SP1cstr=...; YPA1cstr=...; XB1cstr=...;
% Valuta CP,SP,YPA,XB per il PFR
CP1pfr=CP(ICPmax); SP1pfr=SP(ICPmax);
YPA1pfr=YP(ICPmax); XB1pfr=XBP(ICPmax);
% CSTR 2:
thi=tauPFR(ISPmax);
...
Algoritmo – dettagli (7) [facoltativo]
Calcolare una serie di CSTR di dimensioni
crescenti, analizzare i risultati come al pto
3 e confrontare con il PFR
% Ciclo CSTR:
thf=...; % Massimo tempo di permanenza. A scelta, come nel PFR
dth=...; % Incremento del theta da un CSTR al successivo
th=1:dth:thf;
Ccstr=[];
for i=1:length(th)
thi=th(i);
% Chiama il risolutore di eq. algebriche non lineari “fsolve”
(vedi reattore CSTR)
Ccstr=[Ccstr C];
end
% Analizzare i risultati come in dettagli(5) e confrontare con PFR
Osservazioni sulla relazione
Stendere un sintetico resoconto, senza
dilungarsi nella teoria
Preferire le figure alle tabelle, laddove
possibile, e ridurre al minimo la quantità
di figure raggruppando il più possibile i
risultati
Allegare il programma
Scarica
