UNIVERSITÀ DEGLI STUDI DEL SANNIO

FACOLTÀ DI INGEGNERIA

CORSO DI LAUREA MAGISTRALE IN

INGEGNERIA ENERGETICA

ESAME DI “DINAMICA E CONTROLLO DEI PROCESSI ENERGETICI”

Studio di un CSTR in regime transitorio non isotermo che scambia con aria a temperatura ambiente o con un jacket di raffreddamento

Studenti:

Aldo Mottola

Gianmarco Aversano

Nicola Minutella

Stefano Guerriero

1: CSTR IN REGIME TRANSITORIO NON ADIABATICO CHE SCAMBIA CON ARIA A TEMPERATURA AMBIENTE

Si consideri un CSTR nel quale ha luogo la reazione irreversibile esotermica del primo ordine:

A Br(C, T) = - k(T)CA

Il sistema nella forma adimensionale che descrive la dinamica delle variabili di stato è il seguente:

In cui:

β = Entalpia di reazione adimensionale = ;

= Energia di attivazione adimensionale = ;

Da = numero di Damköhler = τkoe-γ = τkrif ;

ϕo = coefficiente di scambio termico adimensionale = ;

= temperatura dell’aria esterna adimensionale;

= temperatura adimensionale del fluido della corrente di alimentazione del CSTR;

= concentrazione adimensionale nella corrente di alimentazione al CSTR.

Lo scopo di quest’analisi è quello di osservare la dinamica nell’intorno del punto di regime che si ha per valori di:

β = 0.81;

= 16;

Da = 0.16;

ϕo = 13;

= 0;

= 0;

= 0.

Si è visto che per tali parametri nasce un regime oscillante, il che mostra la presenza di un ciclo limite. Utilizzando il codice “pplane8” implementato in Matlab è stato ricavato il diagramma delle fasi e le serie temporali delle variabili di stato del sistema per diversi set di condizioni iniziali.

Fig. 1.1: Diagramma delle fasi

Da questo diagramma si evince la presenza di un ciclo limite e di un punto di regime.

Tale punto compare per valori di:

Gli autovalori della matrice jacobiana sono complessi e coniugati con parte reale maggiore di 0, ciò indica che questo punto è un fuoco instabile.

Il ciclo limite, che ha un periodo di oscillazione 5.25τ, è invece un attrattore. Ciò si può vedere anche analizzando le traiettorie ottenute sul diagramma delle fasi per diverse condizioni iniziali.

Fig. 1.2: Traiettorie sul diagramma delle fasi

Notiamo come tutte le traiettorie convergano al ciclo limite, il che conferma la sua stabilità. Questo può essere apprezzato meglio rappresentando le serie temporali. In particolare, notiamo che se il punto iniziale è preso fuori da ciclo limite allora l’ampiezza delle oscillazioni si smorza durante il transitorio fino a raggiungere l’ampiezza del ciclo limite stesso. Se il punto è invece scelto all’interno del ciclo allora l’ampiezza delle oscillazioni si amplifica fino al valore di regime.

Fig. 1.3: Serie temporale variabili di stato con condizioni iniziali esterne al ciclo

Fig. 1.4: Serie temporale variabili di stato con condizioni iniziali interne al ciclo

2: CSTR IN REGIME TRANSITORIO CON JACKET DI RAFFREDDAMENTO

Si consideri un CSTR nel quale ha luogo la reazione irreversibile esotermica del primo ordine:

A Br(C,T) = - k(T)CA

Il calore generato dalla reazione è rimosso da una camicia di raffreddamento (jacket) di volume costante W percorsa da acqua di raffreddamento a flusso volumetrico con temperatura in ingresso Twin e di uscita Tw.

Nella trattazione si fanno le seguenti ipotesi:

· Continuità delle portate volumetriche nel reattore e nel jacket;

· Jacket perfettamente miscelato;

· Calori specifici costanti con la temperatura;

· Volume del reattore costante;

· Densità nel reattore e nel jacket costanti;

· Si trascurano gli effetti termici della parete metallica della camicia.

Iniziamo scrivendo l’equazione di bilancio di materia per un CSTR non isotermo con singola reazione in regime transitorio:

In cui

= Portata volumetrica in ingresso e uscita dal reattore;

= Concentrazione di A in ingresso al reattore;

= Concentrazione nel reattore = Concentrazione in uscita;

V = Volume del reattore;

= tempo dimensionale.

Utilizzando l’espressione di Arrhenius per la costante cinetica di reazione e portando al primo membro il termine differenziale si ottiene:

Abbiamo bisogno ora di scrivere altre 2 equazioni di bilancio, cioè l’equazione di bilancio di energia sul reattore e sul jacket. L’equazione di bilancio di energia sul reattore è:

In cui:

Ctot = Concentrazione totale;

= Calore specifico molare;

Tin = Temperatura della corrente in ingresso al reattore;

T = Temperatura nel reattore = Temperatura corrente in uscita dal reattore;

U = Coefficiente di scambio termico globale fra reattore e jacket;

S = Area superficie esterna del reattore;

Tw = Temperatura nel jacket (T water);

V(∆Hr)r(C,T) = Potenza termica liberata dalla reazione;

L’ultima equazione che ci resta da scrivere è il bilancio di energia sul jacket, che può essere scritto come:

In cui;

ρ = densità del fluido nel jacket;

cp,w = calore specifico del fluido;

Temperatura fluido in ingresso nel jacket;

Uest = Coefficiente di scambio termico globale fra jacket e ambiente;

= Area della superficie esterna del jacket;

= Temperatura ambiente;

= Portata di fluido in ingresso e in uscita al jacket;

W = volume del jacket;

Trascuriamo in via preliminare lo scambio fra jacket e ambiente esterno (scambio gas-solido, molto meno efficace dello scambio liquido-solido). Portando al primo membro il termine differenziale e definendo anche per il jacket un tempo di residenza τw = W / , ottengo:

A questo punto il sistema di equazioni dimensionali da risolvere è:

Essendo a conoscenza dei vantaggi dell’adimensionalizzazione, quali la riduzione del numero di parametri da controllare, quindi la maggior semplicità di analisi che ne consegue, nonché la universalità che il modello adimensionale assume, si procede con la ricerca di un sistema di equazioni adimensionali immagine di quello sopra riportato. Per fare ciò si definiscono delle grandezze di riferimento Crif e Trif scelte in maniera arbitraria. Definiamo quindi:

· ;

· = Temperatura nel reattore adimensionale =

in cui ∆Tad rappresenta il massimo incremento di temperatura raggiungibile nel reattore, in condizioni adiabatiche e per conversione totale a partire da una concentrazione di reagente uguale a quella di riferimento (caso ideale);

· = Temperatura nel jacket adimensionale = ;

·

Sostituendo nell’equazioni e raggruppando i termini si ottiene:

In cui:

β = Entalpia di reazione adimensionale = ;

= Energia di attivazione adimensionale = ;

Da = numero di Damköhler = τkoe-γ = τkrif ;

ϕo = coefficiente di scambio termico adimensionale = ;

Ω = Rapporto tra i tempi di residenza = ;

λ = Rapporto tra le capacità termiche dei due fluidi = .

ANALISI DEI RISULTATI

Utilizzando Matcont si procede al calcolo delle soluzioni di regime per le diverse variabili di stato variando i parametri adimensionali Ω e , perché sono questi due parametri che contengono informazioni riguardo la portata e la temperatura del fluido in ingresso al jacket.

Nelle simulazioni effettuate sono stati considerati valori costanti per gli altri parametri, perché non si vogliono variare né le condizioni di alimentazione del fluido nel reattore né il tipo di reazione.

Iniziamo l’analisi con il calcolo delle soluzioni di regime al variare di Ω, nei diagrammi che seguono riportiamo i risultati.

Fig 2.1. Diagramma delle soluzioni per x al variare di Ω. Con Da=0.16 , β=0.81 , γ = 16 , ϕo = 13 , λ= 0.01 , xin=0, θin=0 , =0.

In Fig. 2.1 si nota una diminuzione del valore di x all’aumentare di Ω , perché un aumento di Ω può rappresentare un aumento di portata alimentata al jacket e questo permette di aumentare la quantità di energia termica asportata al reattore con conseguente riduzione della temperatura nel CSTR e quindi riduzione della velocità di reazione, perché è chiaro che all’aumentare di Ω si limita la possibilità che la temperatura del jacket aumenti inseguendo quella del CSTR.

Nel diagramma viene evidenziata la presenza di un punto di biforcazione di Hopf, da cui nascono dei cicli limite che si conservano all’aumentare di Ω.

Fig 2.2. Diagramma delle soluzioni per θ al variare di Ω. Con Da=0.16 , β=0.81 , γ = 16 , ϕo = 13 , λ= 0.01 , xin=0, θin=0 , =0.

Il diagramma delle soluzioni di θ contro Ω si ripete in maniera analoga a x, in particolare la diminuzione di θ all’aumentare di Ω è sempre dovuta all’aumento dell’energia termica asportata al reattore.

Fig 2.3. Diagramma delle soluzioni per θw al variare di Ω. Con Da=0.16 , β=0.81 , γ = 16 , ϕo = 13 , λ= 0.01 , xin=0, θin=0 , =0.

Anche in questo caso ritroviamo un andamento simile ai precedenti, ma qui ci sono oscillazioni meno evidenti perché si considera un λ=0.01, ossia la capacità termica del fluido nel jacket è scelta circa 100 volte superiore la capacità termica del fluido nel reattore. La maggiore capacità termica consente di contenere in un range molto minore le oscillazioni di temperatura del fluido nel jacket.

Una volta visto che le variabili di stato assumono essenzialmente 2 comportamenti distinti al variare di Ω, cioè quello di punto fisso (per Ω<0.192) e quello di regime periodico (per Ω >0.192), passiamo ad analizzare il comportamento delle soluzioni di regime al variare di per 2 distinti valori di Ω.

Figura 2.4. Diagramma delle soluzioni per x al variare di . Con Da=0.16 , β=0.81 , γ = 16 , ϕo = 13 , λ= 0.01 , xin=0, θin=0 , Ω=0.01.

Nel diagramma possiamo individuare 3 zone distinte: una prima zona in cui ho un unico punto fisso stabile che si mantiene su valori di x bassi; una seconda zona in cui ho tre soluzioni di regime, due delle quali sono stabili e una terza, posta per valori intermedi di x, instabile. Su questa troviamo la nascita di un ciclo limite ‘half-stable’. Nella terza ed ultima zona abbiamo un unico punto fisso con valori di x elevati.

Tale andamento è conseguenza del fatto che il diagramma delle soluzioni presenta 2 biforcazioni sella-nodo che comportano la nascita di un ciclo di isteresi.

Figura 2.5. Diagramma delle soluzioni per θ al variare di . Con Da=0.16 , β=0.81 , γ = 16 , ϕo = 13 , λ= 0.01 , xin=0, θin=0 , Ω=0.01.

Nel grafico l’andamento di θ contro si ripete in maniera qualitativamente analoga a quello visto in figura 4.

Figura 2.6. Diagramma delle soluzioni per θw al variare di . Con Da=0.16 , β=0.81 , γ = 16 , ϕo = 13 , λ= 0.01 , xin=0, θin=0 , Ω=0.01.

Anche l’andamento di θw in funzione di segue lo stesso andamento del grafico in figura 4, in questo caso le oscillazioni sono piccole perché valgono le stesse osservazioni fatte per la figura 3.

Conclusa l’analisi per Ω =0.01 si propone una nuovo studio delle variabili di stato al variare di per un valore di Ω = 1 (ossia maggior del valore critico), che tra l’altro è il caso con τ = τw.

Figura 2.7. Diagramma delle soluzioni per x al variare di . Con Da=0.16 , β=0.81 , γ = 16 , ϕo = 13 , λ= 0.01 , xin=0, θin=0 , Ω=1.

La differenza che notiamo in questo caso rispetto al caso analogo per Ω = 0.01 descritto in Fig. 4, è che scompaiono le biforcazioni sella-nodo. Resta però la biforcazione di Hopf che ha come conseguenza la nascita di un ciclo limite, stavolta stabile.

L’andamento di θ e θw è qualitativamente lo stesso.

Figura 2.8. Diagramma delle soluzioni per θ al variare di . Con Da=0.16 , β=0.81 , γ = 16 , ϕo = 13 , λ= 0.01 , xin=0, θin=0 , Ω=1.

Figura 2.9. Diagramma delle soluzioni per θw al variare di . Con Da=0.16 , β=0.81 , γ = 16 , ϕo = 13 , λ= 0.01 , xin=0, θin=0 , Ω=1.

In particolare in Fig. 9 notiamo che l’andamento di θw con è lineare e rappresenta la bisettrice del primo quadrante perché il fluido permane nel jacket un tempo insufficiente affinché esso possa subire variazioni di temperatura significative e quindi ne esce alla temperatura a cui entra.

ANALISI DEI REGIMI OSCILLANTI

Abbiamo ritenuto utile analizzare a questo punto le variazioni di periodo e ampiezza dei diversi regimi periodici per capire come variano i range di funzionamento delle diverse variabili di stato al variare dei parametri e Ω.

I risultati ottenuti con un codice scritto in Matlab, riportato in allegato alla presente relazione, sono i seguenti:

Figura 2.10. Serie temporale delle 3 variabili di stato. Con Da=0.16 , β=0.81 , γ = 16 , ϕo = 13 , λ= 0.01 , xin=0, θin=0 ,=0, Ω=1.

Il diagramma rappresenta le serie temporali centrate rispetto al proprio valore medio. Per tutte le variabili di stato è stato calcolato il periodo di oscillazione e il ritardo tra le diverse grandezze. Il periodo risulta essere uguale per tutte le grandezze e pari a 4.9 volte τ. Mentre il ritardo tra le grandezze mostra che la variazione θw è in anticipo ed è seguita in ordine da una variazione di x e una di θ. In particolare x è in ritardo rispetto a θw di 0.06 τ; mentre il ritardo θ rispetto ad x è di 0.09 τ.

Figura 2.11.Diagramma Periodo di oscillazione in funzione di Ω. Con Da=0.16 , β=0.81 , γ = 16 , ϕo = 13 , λ= 0.01 , xin=0, θin=0 =0.

La prima zona rappresenta soluzioni di regime a punto fisso per cui non ha alcun senso definire il periodo di oscillazione. Il punto di minimo della curva si registra nel punto di biforcazione di Hopf. Si può osservare come il periodo di oscillazione sia crescente all’aumentare di Ω.

ANALISI DELLE SERIE TEMPORALI E DEI TRANSITORI

Analizzate le possibili soluzioni di regime raggiungibili del sistema CSTR-jacket per il considerato range di variazione dei parametri e Ω, si vogliono adesso simulare i regimi transitori e, quindi, le evoluzioni nel tempo del sistema, con cui evidentemente questo si porta a regime. Nelle immagini che seguono, la scala in ascissa è stata limitata, nonostante le grandezze di stato stiano ancora crescendo, per ragioni di visibilità delle curve.

Fig. 2.12 Andamento della variabile x nel caso di CSTR con jacket, per un valore di Ω lontano dal valore critico.

Fig. 2.13 Andamento della variabile ϑ nel caso di CSTR con jacket, per un valore di Ω lontano dal valore critico.

Fig. 2.14 Andamento della variabile ϑw nel caso di CSTR con jacket, per un valore di Ω lontano dal valore critico.

Nelle tre figure sovrastanti si può infatti vedere come, per il set di valori di parametri riportato in legenda, il sistema evolve da condizioni iniziali ‘spente’ fino alla soluzione di regime di punto fisso determinata anche dalla analisi parametrica delle soluzioni di regime esposta nei paragrafi precedenti. Il fatto che il sistema si porti, in generale, ad una condizione di funzionamento ‘acceso’ pur partendo da condizioni iniziali ‘spente’ (il sistema reagisce fortemente già a temperatura ambiente, ϑ=0) è causa dei valori impostati per i parametri Da, γ, β, la cui variazione non è interesse di questo studio.

È possibile notare un transitorio oscillante dovuto probabilmente ad un conflitto tra la temperatura dell’acqua nel jacket, bassa, e il calore sviluppato dalla reazione nel CSTR. Come evidente, le grandezze di stato tendono, comunque, già durante le oscillazioni verso il regime stazionario finale, a prova del fatto che è l’energia termica della reazione a governare maggiormente l’evoluzione del sistema, nonostante il raffreddamento del jacket, che comunque rappresenta un ‘ostacolo’ da superare negli istanti iniziali. Osservando l’andamento delle variabili di stato (Fig. successiva) nel caso di CSTR adiabatico (x che si porta con un gradino, come la temperatura, alla soluzione di regime stazionario), la presenza, invece, di un transitorio di questo tipo nel caso esaminato non può che confermare questa analisi.

Fig. 2.15 Andamento della variabile x nel caso di CSTR adiabatico.

Come già detto in precedenza, aumentando il valore del parametro Ω, il sistema passa dall’avere una soluzione stabile di punto fisso all’avere una soluzione dinamica oscillante stabile. È interessante allora vedere, con l’avvicinarsi di Ω al valore critico, come e se il transitorio o in generale la serie temporale della soluzione vari. Sono state a tal fine condotte delle simulazioni del sistema per un valore di Ω prossimo a quello critico, pari a circa 0.192, da cui si è notato che il sistema in prossimità di quel valore inizi ad assumere un transitorio anziché grossomodo crescente, come nel caso precedente, oscillante smorzato intorno al regime di punto fisso finale, quasi come se esso fosse ‘pronto’ ad accogliere quanto accade per Ω ≥ 0.192.

Fig. 2.16 Andamento della variabile di stato x nel caso di CSTR con jacket, per un valore di Ω prossimo a quello critico.

Nel passaggio, infatti, di Ω a ridosso del valore critico, quello che accade è che il tempo di smorzamento di questa soluzione oscillante inizia ad essere sempre più lungo, fino a diventare infinito.

Fig. 2.17 Andamento della variabile di stato x nel caso di CSTR con jacket, per un valore di Ω leggermente superiore di quello critico.

Superato il valore di Ω critico la soluzione di regime diventa periodica, e all’aumentare di Ω l’ampiezza dell’oscillazione del ciclo limite aumenta. Si noti infatti che all’aumentare di Ω la temperatura del jacket inizia a diventare costante (e sempre più pari a quella in ingresso al jacket) e pertanto il comportamento del sistema inizia, ovviamente, a imitare quello di un CSTR che scambia con l’aria ambiente, in quanto in entrambi i casi la temperatura della sorgente esterna è ininfluenzabile.

Fig. 2.18 Andamento delle variabili di stato x e ϑ nel caso di CSTR con jacket, per un valore di Ω molto superiore rispetto a quello critico.

Fig. 2.19 Andamento della variabile di stato ϑw nel caso di CSTR con jacket, per un valore di Ω molto superiore rispetto a quello critico.

clear all; close all; clc;

%%Time history of a CSTR in terms of conversion X, reactor temperature Theta and

%%jacket temperature Thetaw

%%

%%Parameters:

par.gamma=16; %16 Ea/RTrif

par.beta=0.81; %0.81 Crif*deltaHr/cp*Trif*Ctot

par.phi=13; %13 US/V*Ctot*cp*krif

par.omega=1; % tau/tauw

par.lambda=0.01; % V*Ctot*cp/W*rho*cw

par.Da=0.16; %0.16 tau*krif with krif=k0*exp(-gamma)

%%

%B.C.:

par.x=0;

par.theta=0;

par.thetaw=0;

%I.C.:

ic.x=0.6;

ic.theta=0;

ic.thetaw=0;

%%

%Solver tools:

par.tend=200; n=100;

%tspan=linspace(0,par.tend,n); %dimensionless time vector

tspan=[0 par.tend];

Y=[ic.x, ic.theta, ic.thetaw];

%%

%Solver:

options=odeset('RelTol',1e-6,'AbsTol',1e-12);

sol=ode45(@ODE,tspan,Y,options,par);

x=sol.y(1,:);

theta=sol.y(2,:);

thetaw=sol.y(3,:);

%%

%PLOT:

%figure(1);plot(sol.x,x);xlabel('t'); ylabel('x');

%figure(2);plot(sol.x,theta);xlabel('t'); ylabel('theta');

%figure(3);plot(sol.x,thetaw);xlabel('t'); ylabel('thetaw');

plot(sol.x,x,sol.x,theta,sol.x,thetaw); legend('x','theta','thetaw'); xlabel('time');

% %%

% %Post-processing:

% %evaluating means of outputsignals:

% ave.x=x(1,end/2:end)-mean(x(1,end/2:end));

% ave.theta=theta(1,end/2:end)-mean(theta(1,end/2:end));

% ave.thetaw=thetaw(1,end/2:end)-mean(thetaw(1,end/2:end));

%

% sol.x2=sol.x(1,end/2:end);

% plot(sol.x2,ave.x,sol.x2,ave.theta,sol.x2,ave.thetaw); legend('x','theta','thetaw'); xlabel('time');

%

% %Cheking periods of the output signals:

% % counter=0;

% % for i=1:length(ave.x)-1

% % if (ave.x(i)*ave.x(i+1)<0)

% % counter=counter+1;

% % end

% % end

% % counterx=counter; counter=0;

% %

% % for i=1:length(ave.x)-1

% % if (ave.theta(i)*ave.theta(i+1)<0)

% % counter=counter+1;

% % end

% % end

% % countertheta=counter; counter=0;

% %

% % for i=1:length(ave.x)-1

% % if (ave.thetaw(i)*ave.thetaw(i+1)<0)

% % counter=counter+1;

% % end

% % end

% % counterthetaw=counter; counter=0;

%

% %disp(counterx); disp(countertheta); disp(counterthetaw);

%

% %Finding the period:

% check=(ave.x(1:end-1)).*(ave.x(2:end));

% ind=find(check<0);

% stepsx=sol.x2(ind(:));

%

% check=(ave.theta(1:end-1)).*(ave.theta(2:end));

% ind=find(check<0);

% stepstheta=sol.x2(ind(:));

%

% check=(ave.thetaw(1:end-1)).*(ave.thetaw(2:end));

% ind=find(check<0);

% stepsthetaw=sol.x2(ind(:));

%

% period.x=round(mean(stepsx(3:end)-stepsx(1:end-2))*10)/10;

% period.theta=round(mean(stepstheta(3:end)-stepstheta(1:end-2))*10)/10;

% period.thetaw=round(mean(stepsthetaw(3:end)-stepsthetaw(1:end-2))*10)/10;

%

% %Evaluating the periods lag:

% lag.x_theta=mean(stepsx-stepstheta)/(period.x);

% lag.x_thetaw=mean(stepsx-stepsthetaw)/(period.x);

% lag.thetaw_theta=mean(stepsthetaw-stepstheta)/(period.x);

function ydot = ODE(t, Y, par)

% system dimensionless of cstr with jacket

ydot=zeros(3,1);

%disp(t/par.tend)

ydot(1,1)=(par.x-Y(1))+(1-Y(1))*par.Da*exp(par.gamma*par.beta*Y(2)/(1+par.beta*Y(2)));

ydot(2,1)=(par.theta-Y(2))+(1-Y(1))*par.Da*exp(par.gamma*par.beta*Y(2)/(1+par.beta*Y(2)))-par.Da*par.phi*(Y(2)-Y(3));

ydot(3,1)=(par.thetaw-Y(3))*par.omega+par.phi*par.Da*par.lambda*(Y(2)-Y(3));

end