Calcolo della tensione di vapore

Oggetto:Calcolo della tensione di vapore di un componente puro da dati critici utilizzando equazioni di stato cubiche
Difficoltà:media
Software:maxima,scilab

Descrivere il comportamento volumetrico di una sostanza è stato
e continua ad essere un problema di non semplice soluzione,o
meglio se ci sono soluzioni semplici,queste valgono in campi
ristretti ed in certe condizioni ben precise.Come è noto le
equazioni di stato di tipo cubico (Van der Waals (VDW),Redlich
Kwong,Peng-Robinson e tutte le altre derivate) sono di uso comune
per descrivere il comportamento dei fluidi.Il problema
nell’utilizzarle per effettuare previsioni è sempre lo stesso:
i dati sperimentali sono scarsi. Non è noto a tutti che alcune
di esse,anche più complicate della pur sempre venerabilissima
VDW,possono essere utilizzate semplicemente conoscendo due dei
parametri critici fondamentali,parametri facilmente misurabili
e localizzabili in molte tabelle.

In questo esercizio vedremo come arrangiare una di queste equazioni
costitutive per esprimerla in funzione dei parametri critici.
Inoltre calcoleremo la tensione di vapore ad una temperatura
assegnata,facendo uso della condizione per la quale fase liquida
e vapore hanno la stessa energia libera molare,questo scrivendo
un po’ di codice scilab.

Come terzo obiettivo plotteremo un grafico della tensione di vapore
con la temperatura.

Infine calcoleremo il calore latente di vaporizzazione.

Iniziamo.
L’equazione di stato che utilizziamo è la Van der Waals,che didatticamente
è più comoda,ma non deve essere utilizzata per fare alcuna previsione (perchè

non sarebbe accurata):

dove V è il volume molare.
Come detto sopra vogliamo i parametri a e b in funzione di parametri di
più semplice approvvigionamento quali ad esempio Tc e Pc
(temperatura e pressione critiche).
Per far questo è sufficiente ricordare che alla temperatura e
pressione critiche c’è il passaggio gas reale - gas perfetto,e
da un punto di vista analitico la funzione scritta sopra presenta
un flesso a tangente orizzontale,per cui si possono scrivere le
seguenti equazioni:

Esplicitando le 3-4:

Un siffatto sistema può essere risolto al volo con Maxima,clicchiamo
su “equazioni->risolvi sitema algebrico” ed inseriamo le equazioni e
le variabili,diamo ok:

quindi abbiamo trovato come variano a e b con Tc e Vc:

sostituendo le (7) nella (2) otteniamo:

ricordo che vogliamo a(Pc,Tc) e b(Pc,Tc),dobbiamo fare qualche altro
passaggio,intanto ricaviamo il fattore di compressibilità critico:

ottenuta “mescolando” le (7) e (8), questo significa che tutti i
fluidi che seguono la VDW hanno lo stesso fattore di compressibilità
critico,mentre come si sa i fluidi che seguono la legge del gas
perfetto hanno Z=1 per ogni valore di T,quindi la VDW è già un
notevole miglioramento nella descrizione volumetrica.
Ora sosituendo la (9) nelle (7) otteniamo:

dunque abbiamo i valori delle costanti a e b in funzione di Tc e Pc,
e ci siamo.

Possiamo dunque scrivere l’equazione di VDW generalizzata:

cioè abbiamo un’equazione che ci permette di prevedere il comportamento
volumetrico conoscendo solo le proprietà critiche della sostanza e non
i parametri a e b “fittati” da dati sperimentali.

Ora passiamo alla seconda parte dell’esercizio,ricordo che dobbiamo trovare
la tensione di vapore di una sostanza nota la temperatura.Esistono equazioni
(vedi Antoine) che permettono di calcolare la tensione di vapore ad una
determinata temperatura,ma anche in questo caso ci sono dei parametri di
non facile approvvigionamento.

Cominciamo col dire che la tensione di vapore ad una data temperatura,è
la pressione di equilibrio,cioè la pressione tale che la fase vapore e
quella liquida presentano una eguale energia libera molare:

Dobbiamo esprimere questa uguaglianza in termini di P e V per poter utilizzare
la nostra equazione di stato (EOS) generalizzata,ricordiamo che:

all’equilibrio la temperatura e la pressione (di vapore) nelle due fasi è costante
(dT=0,dP=0):

ma noi abbiamo una EOS in cui P è funzione di V e non viceversa,allora:

ora non resta che integrare tra la fase liquida e quella vapore:

siccome P è la tensione o pressione di vapore:

Vv e Vl sono i volumi molari di vapore e liquido che evidentemente sono diversi.
Dunque ci siamo,possiamo scrivere il sistema che ci permette di calcolare Ps:

E’ un sistema di 3 equazioni nelle 3 incognite Ps,Vl,Vv,si potrebbe utilizzare
un metodo numerico per la soluzione di sistemi non lineari,ma non è
consigliabile,a causa della alta non linearità.

Procederemo per bisezione.
Prima scriviamo la terza delle (1 in questo modo:

l’integrale a secondo membro lo risolviamo subito,inseriamo la (1) in maxima
(sostituiamo in seguito ad a e b le loro funzioni di Pc e Tc) ed integriamo,
si ottiene:

Fissiamo due valori adeguati di Ps per determinare un intervallo di ricerca
[Pa,Pb],ci mettiamo in mezzo a questo intervallo per avere un primo valore
di Ps,Ps*.
Ci calcoliamo i volumi Vv e Vl per Pa e Ps* con le formule di Cardano
(già implementate come funzione in “ztools”),e quindi con la (20)
phi_a e phi*.Per calcolare i volumi,dobbiamo riscrivere la (11) come
un polinomio ed esplicitarne i coefficienti:

Valutiamo il prodotto tra phi_a e phi*,se risulta positivo allora la radice
non è contenuta tra Pa e Ps* ma sarà tra Ps* e Pb,si pone Pa=Ps*,si dimezza
e si rivaluta il prodotto e così via,finchè non è soddisfatto il criterio di
arresto:

chiaramente faremo fare tutto in automatico da un programma scritto in Scilab.
La scelta del valore di Ps di primo tentativo è relativamente semplice.
Sappiamo che a questa pressione devono corrispondere due volumi,quello del
liquido e quello del vapore.L’unica zona in cui la funzione P(V,T) ha due
radici (tre,una è senza significato fisico*,più esattamente quella centrale)
è quella nel range Pmin e Pmax che corrispondono a Vmin e Vmax.Quindi sarà
sufficiente calcolarci i valori Vmax e Vmin di P(V,T) e quindi Pmax e Pmin
e scegliere un valore intermedio.

Per calcolarci Vmax e Vmin annulliamo la
derivata di P(V,T) fatta rispetto a V:

Riarrangiamolo con Maxima,”semplifica->scomponi in fattori”:

dividendo per -RT si ha:

è un polinomio cubico,i cui coefficienti sono:

le sue radici si possono calcolare mediante le formule di Cardano.
Risolto il problema di trovare Vmin e Vmax,ci calcoliamo Pmin e
Pmax dalla (11) sostituendo i valori del volume appena trovati.

Non possiamo scegliere Pmax e Pmin come Pa e Pb proprio perchè sono

punti estremali.Possiamo scegliere Pa e Pb in questo modo:

Ma siccome la pressione deve essere positiva,possiamo scegliere come Pb:

La scelta del delta è arbitraria,ne basta uno “sufficientemente piccolo”
,ad esempio come percentuale dell’intervallo Pmax-Pmin:

Il calcolo del calore di vaporizzazione alla temperatura T si può
calcolare sapendo che in condizioni isoterme-isobare
(cioè nella condizione di equilibrio liq-vap) si ha:

Da qui si procede in questo modo,utilizzando le equazioni Maxwell:

A questo punto ci siamo,non rimane altro da fare che scrivere il
codice in scilab.Il nostro programmino avrà la struttura di una
funzione (che chiamo ps_vdw) che avrà come INPUT la temperatura
T alla quale vogliamo la tensione di vapore,e le caratteristiche
critiche della sostanza cioè Pc e Tc e come OUTPUT il valore
della tensione di vapore.

Iniziamo:

//dichiaro la funzione dove il valore ritornato è ps e gli argomenti
//sono quelli tra parentesi

function [ps]=ps_vdw(T,Tc,Pc)
//dichiaro alcune costanti,eps è la tolleranza per il check,R è la costante dei gas
//m è il parametro della (30)

eps=0.00000001;
R=8.31;
m=0.01;

//calcolo a e b dai dati critici con le (10):

a=27*R^2*Tc^2/(64*Pc);
b=R*Tc/(8*Pc);

//dichiarazione della funzione vdw (1):

function [p]=vdw(v,T,a,b)
p=(8.31*T/(v-b))-a/v^2;
endfunction

//vettore coefficienti del polinomio per vmax vmin (26):

c_1=zeros(3);
c_1(1)=-2*a/(R*T);
c_1(2)=4*a*b/(R*T);
c_1(3)=-2*a*b^2/(R*T);

//vettore dei coefficienti del polinomio per vv e vl (21)(dipende da pa,ps o pb):

function [c_2]=cp(T,p,a,b)
R=8.31;
c_2=zeros(3);
c_2(1)=-R*T/p-b;
c_2(2)=a/p;
c_2(3)=-a*b/p;
endfunction

//dichiaro la funzione phi che è quella con la quale bisezioniamo (20):

function [f]=phi(T,p,vl,vv,a,b)
R=8.31;
f=p*(vv-vl)-R*T*log((vv-b)/(vl-b))-a*(1/vv-1/vl);
endfunction

//calcolo i limiti vmax e vmin,per farlo uso la funzione “card” di “ztools”
//cioè il metodo risolutivo di Cardano per polinomi cubici,vm è il vettore
//soluzione (le tre radici),la funzione “sort” ordina le radici dalla più
//grande alla più piccola,a noi interessano le prime 2
//poi calcolo pmax e pmin

vm=card(c_1);
vm=sort(vm);
vmin=vm(2);
vmax=vm(1);
pmin=vdw(vmin,T,a,b);
pmax=vdw(vmax,T,a,b);

//mi calcolo il delta della (28),calcolo pa e verifico che sia positivo con
//un ciclo che dimezza m (e quindi delta) finchè non lo è
//calcolo pb

delta=m*(pmax-pmin);
pa=pmax-delta;
while(pa<0)
m=m/2;
delta=m*(pmax-pmin);
pa=pmax-delta;
end;
pb=0+delta;

//mi calcolo i coefficienti del polinomio a pb
//e ricavo il vettore delle radici
//uso la funzione size() per conoscere la dimensione del vettore v_b
//(che deve essere 3)
//nota:sb è il vettore in cui size() mette due valori,nella prima componente
//c’è la dimensione di v_b

cb=cp(T,pb,a,b);
v_b=card(cb);
sb=size(v_b);

//siccome devo essere sicuro che in corrispondenza di pb ho tre
//radici reali positive di volume,creo un ciclo che fintantochè
//non è verificata la condizione precedente (tre e positive)
//innalza il valore di pb di un delta (lo stesso di prima)
//che non dovrà scavalcare il valore di pa

while(sb(1)==1|positive(v_b)==%F)
printf(”*”);
pb=pb+delta;
if(pb>pa-2*delta)
break;
end;
cb=cp(T,pb,a,b);
v_b=card(cb);
sb=size(v_b);
s=sb(1);
end;

//ora mi calcolo i volumi di pa e pb,sono sicuro che i vettori v_a e v_b
//contengono 3 radici reali positive;come abbiamo visto mi servono la più
//piccola e la più grande;le funzioni max e min di scilab permettono di estrarle
//v_l è il volume del liquido,v_v è quello del vapore

ca=cp(T,pa,a,b);
v_a=card(ca);
va_l=min(v_a);
va_v=max(v_a);
vb_l=min(v_b);
vb_v=max(v_b);

//mi calcolo il valore di phi in corrispondenza dei volumi ricavati:

phi_a=phi(T,pa,va_l,va_v,a,b);
phi_b=phi(T,pb,vb_l,vb_v,a,b);

//inizio il metodo di bisezione che ha come check proprio la (22)

while(abs(phi_a-phi_b)>eps)

//se il check è verificato,seziono l’intervallo [pa,pb] trovando la
//prima approssimazione della tensione di vapore,calcolo i volumi a ps
//e valuto la phi

ps=(pa+pb)/2;
c=cp(T,ps,a,b);
v_s=card(c);
vs_l=min(v_s);
vs_v=max(v_s);
phi_s=phi(T,ps,vs_l,vs_v,a,b);

//verifico che la radice sia all’interno dell’intervallo [phi_a,phi_s]

value=phi_a*phi_s;

//se è verificato che il prodotto precedente è positivo,allora la radice
//è tra pb e ps,e pongo pa=ps

if (value>0)
pa=ps;
c=cp(T,pa,a,b);
v_a=card(c);
va_l=min(v_a);
va_v=max(v_a);
phi_a=phi(T,pa,va_l,va_v,a,b);
end;

//se invece non è verificato,allora la radice si trova tra pa e ps
//allora pongo pb=ps

if (value<0)
pb=ps;
c=cp(T,pb,a,b);
v_b=card(c);
vb_l=min(v_b);
vb_v=max(v_b);
phi_b=phi(T,pb,vb_l,vb_v,a,b);
end;

//fine della bisezione
end;

//calcolo il calore di vaporizzazione:

lambda=T*8.31*log((vs_v-b)/(vs_l-b));

//fine della funzione
endfunction

Questo è il codice base per ottenere il valore della tensione di vapore alla temperatura scelta.
Ora creiamo la funzione che ci permette di plottare il valore della pressione di vapore per un range di temperatura dato:

//N è il numero dei sottointervalli di T,h è la lunghezza di questi sottointervalli
//x è il vettore le cui componenti sono le temperature,in sx è contenuta la dimensione

function ps_vdw_T(Tc,Pc,T1,T2)
N=30;
h=(T2-T1)/N;
x=[T1:h:T2]‘;

//valuto la tensione di vapore per ogni temperatura

for i=1:N+1
printf(”T=%f \n”,x(i));
y(i)=ps_vdw(x(i),Tc,Pc);
end;

//plotta con triangolini rossi con la punta rivolta a destra
//faccio in modo che sul grafico ci siano temperature in celsius
//e pressioni in atm

plot(x-273,y/101300,’r.->’);
endfunction

Metto questo programmino nella libreria Ztools (lievemente modificato)
Ora caricando Ztools.sci dall’editor e “loadandolo” in scilab,possiamo scrivere
ad esempio per l’esano (T in kelvin,P in pascal)

ps_vdw_T(507.4,2969000,300,500) :

il profilo è corretto,infatti la successione è monotona crescente.

E’ possibile far vedere che l’equazione di Van der Waals non predice correttamente il valore della tensione di vapore per temperature lontane dal punto critico,questo è ragionevole se pensiamo che abbiamo considerato costanti i parametri a e b,e che abbiamo ricavato questi parametri dalle condizioni di flesso nel punto critico.

Ci siamo.