Calcolo dell'energia libera di Gibbs in funzione della temperatura e della pressione

L'esercitazione che segue utilizza un modulo python (gibb_tp.py) che fa parte di una suite di programmi scritti per scopi professionali nell'ambito del modeling termodinamico dei sistemi geochimici e petrologici. Se ne illustreranno solo alcune caratteristiche e, precisamente, quelle utili ai fini di questa esercitazione.

La caratteristica principale del programma consiste nell'uso di una classe python (la classe mineral) per caricare da un database l'insieme completo delle proprietร  termodinamiche di un certo numero di fasi minerali, insieme con tutte le funzioni necessarie per calcolarne l'energia libera o altre grandezze di interesse.

Come vedremo, ogni fase minerale trattata รจ associata a una specifica istanza della classe.

Il codice che definisce la classe mineral รจ:

class mineral:
    def __init__(self,name,nick):
        self.name=name
        self.nick=nick
        self.eos='m'
        self.cp=None
        self.al=None
        self.k0=None
        self.kp=None
        self.dkt=None
        self.v0=None
        self.g0=None
        self.s0=None

    def info(self):
        print("Mineral: %s\n" % self.name)
        print("K0: %4.2f GPa, Kp: %4.2f, dK0/dT: %4.4f GPa/K, V0: %6.4f J/bar" \
              % (self.k0, self.kp, self.dkt, self.v0))
        print("G0: %8.2f J/mol, S0: %6.2f J/mol K\n" % (self.g0, self.s0))
        print("Cp coefficients and powers:")
        for ci in self.cp:
            print('{:>+8.4e}{:>+6.1f}'.format(ci[0],ci[1]))

        print("\nAlpha coefficients and powers:")
        for ai in self.al:
            print('{:>+8.4e}{:>+6.1f}'.format(ai[0],ai[1]))

    def load_ref(self,v0,g0,s0):
        self.v0=v0
        self.g0=g0
        self.s0=s0

    def load_bulk(self,k0, kp, dkt):
        self.k0=k0
        self.kp=kp
        self.dkt=dkt

    def load_cp(self,cpc,cpp):
        cl=list(zip(cpc,cpp))
        item=0
        self.cp=np.array([])
        for ic in cl:
            self.cp=np.append(self.cp,[cl[item][0], cl[item][1]])
            item=item+1
        self.cp=self.cp.reshape(item,2)

    def load_alpha(self, alc, alp):
        cl=list(zip(alc,alp))
        item=0
        self.al=np.array([])
        for ic in cl:
            self.al=np.append(self.al,[cl[item][0], cl[item][1]])
            item=item+1
        self.al=self.al.reshape(item,2) 

    def cp_t(self,tt):
        cpt=0.
        iterm=0
        for it in self.cp:
            cpt=cpt+self.cp[iterm,0]*(tt**self.cp[iterm,1])
            iterm=iterm+1
        return cpt 

    def alpha_t(self,tt):
        alt=0.
        iterm=0
        for it in self.al:
            alt=alt+self.al[iterm,0]*(tt**self.al[iterm,1])
            iterm=iterm+1
        return alt

    def kt(self,tt):
        return self.k0+(tt-298.15)*self.dkt

    def entropy(self,tt):
        fc=lambda ti: (self.cp_t(ti))/ti
        integ, err=scipy.integrate.quad(fc,298.15,tt)
        return integ+self.s0

    def volume_t(self,tt):
        fc= lambda ti: self.alpha_t(ti)
        integ,err=scipy.integrate.quad(fc,298.15,tt)
        return (self.v0)*np.exp(integ)

    def volume_p(self,tt,pp):
        k0t=self.kt(tt)
        vt=self.volume_t(tt)
        if self.eos=='m':
           fact=(1.-(pp*self.kp)/(pp*self.kp+k0t))**(1./self.kp)
           return fact*vt
        elif self.eos=='bm':
           pf=lambda f: (3*k0t*f*(1+2*f)**(5/2)*(1+3*f*(self.kp-4)/3)-pp)**2
           ff=scipy.optimize.minimize(pf,1,tol=0.00001)            
           return vt/((2*ff.x[0]+1)**(3/2))

    def volume_fix_t_p(self,tt):
        return lambda pp: self.volume_p(tt,pp) 

    def vdp(self,tt,pp):
        fv=self.volume_fix_t_p(tt)
        integ,err=scipy.integrate.quad(fv,0.0001, pp)
        return integ*1e4

    def g_t(self,tt):
        integ,err=scipy.integrate.quad(self.entropy, 298.15, tt)
        return integ

    def g_tp(self,tt,pp):
        return self.g0+self.vdp(tt,pp)-self.g_t(tt) 

    def alpha_p(self, tt, pp):
        v=self.volume_p(tt,pp)
        t_list=np.linspace(tt-10, tt+10, 5)
        vt_list=np.array([])
        for ti in t_list:
            vi=self.volume_p(ti,pp)
            vt_list=np.append(vt_list,vi)
        fitpv=np.polyfit(t_list,vt_list,2)
        fitder1=np.polyder(fitpv,1)
        altp=np.polyval(fitder1,tt)
        return 1*altp/v   

    def s_tp(self,tt,pp):
        gtp=lambda tf: self.g_tp(tf,pp)
        t_list=np.linspace(tt-5, tt+5, 5)
        g_list=np.array([])
        for ti in t_list:
            gi=self.g_tp(ti,pp)
            g_list=np.append(g_list,gi)
        fit=np.polyfit(t_list,g_list,2)
        fitder=np.polyder(fit,1)
        return -1*np.polyval(fitder,tt)

    def h_tp(self,tt,pp):
        g=self.g_tp(tt,pp)
        s=self.s_tp(tt,pp)
        return g+tt*s

Cerchiamo di interpretarlo nelle caratteristiche essenziali con qualche esempio.

Tanto per cominciare, lanciamo i due comandi che seguono per avere la grafica inline e caricare il programma

In [238]:
%matplotlib inline
%run gibbs_tp.py

La lista dei minerali che sono trattabili in questo sistema รจ contenuta nella variabile name_list.mineral_names che possiamo visualizzare scrivendo semplicemente:

In [239]:
name_list.mineral_names
Out[239]:
['coe', 'q', 'fo', 'ky', 'sill', 'andal', 'per', 'sp', 'py', 'ens', 'cor']

Le sigle delle diverse fasi minerali sono conformi (con qualche eccezione) a quelle usate internazionalmente ('coe', coesite; 'q', quarzo; 'andal', andalusite; 'per', periclasio; 'sp', spinello, ecc...)

Come scrivevo sopra, ogni fase minerale รจ associata a una istanza della classe; per esempio:

py=mineral("pyrope","py")

crea una variabile py che contiene tutti i dati termodinamici del piropo e le funzioni che li utilizzano. Questi dati possono essere visualizzati con la funzione py.info():

In [240]:
py.info()

Mineral: pyrope

K0: 173.70 GPa, Kp: 4.00, dK0/dT: -0.0261 GPa/K, V0: 11.3180 J/bar
G0: -5934105.00 J/mol, S0: 266.30 J/mol K

Cp coefficients and powers:
+6.3350e+02  +0.0
-5.1961e+06  -2.0
-4.3152e+03  -0.5

Alpha coefficients and powers:
+4.3600e-05  +0.0
-4.3600e-04  -0.5

Abbiamo quindi il bulk modulus $K_0$, il $K'$ (Kp), il $dK_0/dT$... e poi le espressioni che esprimono il calore specifico a pressione costante e l'espansione termica in funzione di $T$. Per esempio, per il piropo abbiamo un'equazione del tipo

$$C_P(T)=633.5 - 5.1961\cdot 10^{6}T^{-2} - 4.3152\cdot 10^3T^{-0.5}$$

per il calore specifico (si noti che questa espressione non va bene per ottenere un $C_P$ valido per le temperature molto basse... perchรจ?); volendo sapere il $C_P$ del piropo a T=500 K non dovremo far altro che usare una delle funzioni della classe e, precisamente, cp_t che accetta, come argomento, la temperatura a cui si voglia fare il calcolo:

In [241]:
py.cp_t(300).round(2)
Out[241]:
326.63

py.cp_t(T) calcola il $C_P$ alla temperatura $T$; il round(2) postfisso ha solo un effetto cosmetico limitando a 2 il numero di cifre decimali. Il dato รจ calcolato in J/mole K.

Possiamo anche mettere in grafico i valori del CP in un dato intervallo di temperature, scrivendo il semplice codice che segue:

In [242]:
t_list=np.linspace(300,1000,40)   # lista di 40 valori di temperatura tra 300 e 1000 K
cp_list=py.cp_t(t_list)           # lista dei corrispondenti valori di CP

plt.figure()
plt.plot(t_list,cp_list)
plt.xlabel("T (K)")
plt.ylabel("CP (J/mole K)")
plt.title("Calore specifico del piropo")
plt.show()

Volendo il $C_P$, che so, della coesite alla stessa temperatura, basterร  richiedere:

In [243]:
coe.cp_t(500).round(2)
Out[243]:
58.74

Per quanto riguarda il piropo, il suo volume molare รจ conservato in py.v0:

In [244]:
py.v0
Out[244]:
11.318

Il suo valore รจ dato in J/bar... unitร  di misura abbastanza strana per un volume! A cosa corrisponde? Vediamo:


1 J = 1 N m; 1 bar 10^5 Pa = 10^5 N/m^2, da cui 1 J/bar = 10^-5 N m m^2/N = 10^-5 m^3 = 10 cm^3

Quindi il volume molare del piropo รจ 113.18 cm^3. Ma perchรจ si usa un'unitร  di misura cosรฌ inusuale?
Perchรจ nel modeling termodinamico tutto viene definito a partire dall'energia libera di Gibbs, e

$$V=\left(\frac{\partial G}{\partial P}\right)_T$$

Poichรจ $G$ รจ dato in J (รจ un'energia) e la pressione รจ data in bar, allora $V$ รจ dato in J/bar.

Siamo curiosi adesso di sapere quale sia il volume molare del piropo a T=500 K, e P=5 GPa (notare che la pressione viene qui data in GPa e non in bar): usiamo la funzione volume_p(T,P):

In [245]:
py.volume_p(500,5).round(3)
Out[245]:
11.053

che utilizza i dati dell'espansione termica e dell'equazione di stato (con il bulk modulus e la sua dipendenza dalla temperatura) per fare il calcolo.

La funzione g_tp(T,P) ci calcola l'energia libera $G$ in date condizioni $T$, $P$

In [246]:
print('{:.4e}'.format(py.g_tp(500,5)))

-5.4481e+06

OK, d'accordo, l'istruzione sopra sempra arabo, ma in realtร  รจ solo un modo per ottenere un numero scritto con notazione scientifica con quattro cifre decimali. Per esempio, se la variabile che vogliamo scrivere fosse:

In [247]:
a=2564.566791

scriveremo un comando di print con incorporate le istruzioni per la formattazione: {:.3e}, che scriveranno la variabile in notazione scientifica (e) con 3 cifre decimali:

In [248]:
print('{:.3e}'.format(a))

2.565e+03

La funzione plot_g_t esegue il plot dell'energia libera a npoint valori di temperatura, un dato range tmin, tmax, a una pressione fissata P: per esempio

plot_g_t(300, 700, 12, 3)

con tmin=300 K; tmax=700 K, npoint=12, P=3 GPa

In [249]:
plot_g_t(300,700,12,3)

All'interno della funzione sono giร  specificate le fasi su cui si vuole fare il calcolo. Qui abbiamo 3 curve:

  • piropo: Mg3Al2Si3O12 (curva in nero)
  • enstatite + corindone: 3/2 Mg2Si2O6 + Al2O3 (curva in rosso)
  • periclasio + corindone + quartzo: 3 MgO + Al2O3 + 3 SiO2 (curva in blu)

Vediamo che (a 3 GPa) l'insieme per+cor+q non รจ mai stabile rispetto ad altre fasi (la sua energia libera รจ sempre la piรน alta a qualunque T, il che รจ vero anche a pressioni piรน alte). Enstatite + corindone (ens+cor) รจ piรน stabile rispetto al piropo a bassa temperatura, mentre ad alta T รจ piรน stabile il piropo (la T di transizione รจ intorno ai 500 K).

Abbiamo una funzione analoga per fare un calcolo di $G$ in funzione della pressione, a temperatura fissata T: plot_g_p(pmin,pmax,npoint,T); qui abbiamo Pmin=0, Pmax=10 GPa (12 punti), T=500 K:

In [250]:
plot_g_p(0,10,12,500)

E' molto semplice modificare queste funzioni di plot per considerare fasi minerali diverse. Per esempio, plot_g_t รจ riscritta nella cella seguente (cella che puรฒ essere eseguita, caricando la nuova funzione in luogo della precedente). Per esempio, vogliamo considerare l'equilibrio cianite (ky), andalusite (andal):

In [251]:
def plot_g_t(tmin,tmax,npoint,pres,ret=False):
    t_list=np.linspace(tmin,tmax,npoint)
      
    g_1=np.array([])
    g_2=np.array([])
    for it in t_list:
        ig=ky.g_tp(it,pres)        # specificare qui la cianite:   ky.g_tp 
        ig2=andal.g_tp(it,pres)    # specificare qui l'andalusite: andal.g_tp
        g_1=np.append(g_1,ig)
        g_2=np.append(g_2,ig2)
        
    fig=plt.figure()
    ax=fig.add_subplot(111) 
    ax.title.set_text("G in funzione di T, a P costante = "+str(pres)+" GPa")
    ax.plot(t_list,g_1,"k-",label="ky")                     # cambiare la label
    ax.plot(t_list,g_2,"r-",label="andal")                  # cambiare la label
    ax.yaxis.set_major_locator(plt.MaxNLocator(5))
    ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%4.2e'))
    ax.xaxis.set_major_locator(plt.MaxNLocator(8))
    ax.xaxis.set_major_formatter(mtick.FormatStrFormatter('%.1f'))
    ax.set_xlabel("Temperatura (K)")
    ax.set_ylabel("G (J/mole)")
    ax.legend(frameon=False)
    plt.show()
    
    if ret:
        return t_list, g_1, g_2   # variabili in uscita se ret=True

... e chiamare la funzione con gli argomenti desiderati:

In [252]:
plot_g_t(500,900,12,0.3)

Vogliamo adesso determinare la temperatura precisa a cui abbiamo la transizione da cianite ad andalusite per una pressione di 0.2 GPa. Richiamiamo la funzione plot_g_t specificando l'argomento opzionale plot=True: La funzione restituisce in tal caso, le liste dei valori di temperatura, e le corrispondenti liste delle energie libere di cianite e andalusite

In [253]:
t, g_ky, g_andal=plot_g_t(300,800,200,0.2,ret=True)

Calcoliamo adesso, temperatura per temperatura, la funzione $dg = (G_{ky}-G_{andal})^2$ e cerchiamone il minimo; il minimo (idealmente 0) รจ dove le energie libere delle due fasi sono uguali e rappresenta dunque la pressione di equilibrio:

In [254]:
dg=(g_ky-g_andal)**2
In [255]:
plt.figure()
plt.plot(t,dg)
plt.xlabel("T (K)")
plt.ylabel("dg")
plt.show()

La posizione nella lista dg a cui la stessa dg รจ minima รจ estratta dalla funzione argmin di numpy

In [256]:
np.argmin(dg)
Out[256]:
132

Enfatizzo che questo numero (65) non รจ il minimo di dg, ma รจ la posizione nella lista alla quale รจ presente quel valore minimo. Il valore minimo di dg si ottiene invece con la funzione min:

In [257]:
min(dg).round(2)
Out[257]:
90.33

Non รจ esattamente 0, a causa della risoluzione con cui abbiamo costruito la lista del valori di temperatura (200 punti in un intervallo di 500 K)

Andiamo a vedere il corrispondente valore di temperatura nella posizione 132 della lista delle temperature:

In [258]:
t[132].round(1)
Out[258]:
631.7

Se vogliamo fare una ricerca piรน accurata con un metodo piรน sofisticato, conoscendo la posizione approssimata del minimo in T, costruiamo le liste t e dg nei dintorni del minimo, poi interpoliamo i risultati con una funzione

$$dg = at^2 + bt + c$$

usando la funzione polyfit di numpy. Calcoliamo poi la posizione del minimo uguagliando a 0 la derivata analitica della funzione dg(t) interpolata: $t_{min}=-b/2a$

In [259]:
tmin=620
tmax=640
t, g_ky, g_andal=plot_g_t(tmin,tmax,10,0.2,ret=True)
dg=(g_ky-g_andal)**2
reg=np.polyfit(t,dg,2)
tt=-1*reg[1]/(2*reg[0])
print("Temperatura di transizione: %5.2f K" % tt)

# Notare l'uso della funzione polyval che accetta i parametri
# a, b e c (conservati in reg e determinati da polyfit)
# insieme a un valore di t, per restituire il valore del
# polinomio appunto al valore di t
print("Minimo di dg: %3.2f" % np.polyval(reg,tt))  

# Parte grafica (opzionale)
new_tlist=np.linspace(tmin,tmax,100)
new_dg=np.polyval(reg,new_tlist)
plt.figure()
plt.plot(t,dg,"k*",label="punti calcolati")
plt.plot(new_tlist,new_dg,"b-",label="funzione interpolante")
plt.legend(frameon=False)
plt.xlabel("T (K)")
plt.ylabel("dg")
plt.show()

Temperatura di transizione: 630.63 K
Minimo di dg: -0.31

Questa รจ la temperatura a cui si ha l'equilibrio cianite/andalusite a una pressione di 0.2 GPa; a T piรน bassa รจ stabile la cianite; a T piรน alta รจ stabile l'andalusite.

Facciamo lo stesso giochetto con l'equilibrio quarzo/coesite, alla temperatura di 300 K. Modifichiamo la funzione plot_g_p

In [260]:
def plot_g_p(pmin,pmax,npoint,temp,ret=False):
    p_list=np.linspace(pmin,pmax,npoint)
      
    g_1=np.array([])
    g_2=np.array([])
    for ip in p_list:
        ig=coe.g_tp(temp,ip)           # specifichiamo qui la coesite
        ig2=q.g_tp(temp,ip)            # specifichiamo qui il quarzo
        g_1=np.append(g_1,ig)
        g_2=np.append(g_2,ig2)
        
    fig=plt.figure()
    ax=fig.add_subplot(111) 
    ax.title.set_text("G in funzione di P, a T costante = "+str(temp)+" K")
    ax.plot(p_list,g_1,"k-",label="coe")                  # cambiamo l'etichetta
    ax.plot(p_list,g_2,"r-",label="q")                    # cambiamo l'etichetta
    ax.yaxis.set_major_locator(plt.MaxNLocator(5))
    ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%4.2e'))
    ax.xaxis.set_major_locator(plt.MaxNLocator(8))
    ax.xaxis.set_major_formatter(mtick.FormatStrFormatter('%.1f'))
    ax.set_xlabel("Pressione (GPa)")
    ax.set_ylabel("G (J/mole)")
    ax.legend(frameon=False)
    plt.show()  
    
    if ret:
         return p_list, g_1, g_2
In [261]:
plot_g_p(0, 6, 12, 300)
In [262]:
p, gq, gc=plot_g_p(2,4,50,300,ret=True)
In [263]:
dg=(gq-gc)**2
print("Pressione di transizione: %3.2f GPa" % p[np.argmin(dg)])

Pressione di transizione: 2.90 GPa

Abbiamo poi una funzione piรน evoluta che calcola la pressione di equilibrio tra diverse fasi minerali, in un dato intervallo di temperatura. Si tratta della funzione equilib che ha come argomenti:

equilib(Tmin, Tmax, npoint, prod, rea)

prod e rea sono liste di prodotti e reagenti con i rispettivi coefficienti stechiometrici. Per esempio, per la reazione

3/2 enstatite + 1 corindone <---> 1 piropo

diamo prod=['py',1] e rea=['ens', 1.5, 'cor',1]

In [264]:
equilib(300,700,12,prod=['py',1], rea=['ens',1.5,'cor',1])

 T (K)  P (GPa)  DH(J/mol)  DS (J/mol K)  DV (J/bar)  Slope (bar/K)
 300.0   3.67     5304.500     17.682       -0.527       -33.54    
 336.4   3.55     5971.184     17.752       -0.531       -33.43    
 372.7   3.43     6650.804     17.844       -0.535       -33.34    
 409.1   3.31     7325.018     17.906       -0.539       -33.20    
 445.5   3.19     7981.999     17.919       -0.544       -32.97    
 481.8   3.07     8613.967     17.878       -0.548       -32.63    
 518.2   2.95     9215.754     17.785       -0.552       -32.20    
 554.5   2.83     9783.931     17.643       -0.557       -31.69    
 590.9   2.72    10316.250     17.458       -0.561       -31.10    
 627.3   2.61    10811.291     17.235       -0.566       -30.45    
 663.6   2.50    11268.213     16.979       -0.570       -29.76    
 700.0   2.39    11686.593     16.695       -0.575       -29.03    



Average Clapeyron Slope (from Delta S/Delta V): -31.94 bar/K
Clapeyron slope (from a linear fit of the P/T curve): -32.13 bar/K

E' noto che lo spinello MgAl2O4, in alta pressione, si decompone in periclasio (MgO) + corindone (Al2O3). Vogliamo indagare quale sia il campo P/T a cui avviene la reazione. Cominciamo con l'esplorare gli andamenti dell'energia libera delle fasi coinvolte. In particolare, modifichiamo la funzione plot_g_p ponendo nelle due variabili g_1 e g_2 in modo che:

$g_1=G_{sp}$
$g_2=G_{cor}+G_{per}$

In [265]:
def plot_g_p(pmin,pmax,npoint,temp,ret=False):
    p_list=np.linspace(pmin,pmax,npoint)
      
    g_1=np.array([])
    g_2=np.array([])
    for ip in p_list:
        ig=sp.g_tp(temp,ip)                       # spinello
        ig2=cor.g_tp(temp,ip)+per.g_tp(temp,ip)   # corindone + periclasio         
        g_1=np.append(g_1,ig)
        g_2=np.append(g_2,ig2)
        
    fig=plt.figure()
    ax=fig.add_subplot(111) 
    ax.title.set_text("G in funzione di P, a T costante = "+str(temp)+" K")
    ax.plot(p_list,g_1,"k-",label="sp")                 
    ax.plot(p_list,g_2,"r-",label="per+cor")
    ax.yaxis.set_major_locator(plt.MaxNLocator(5))
    ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%4.2e'))
    ax.xaxis.set_major_locator(plt.MaxNLocator(8))
    ax.xaxis.set_major_formatter(mtick.FormatStrFormatter('%.1f'))
    ax.set_xlabel("Pressione (GPa)")
    ax.set_ylabel("G (J/mole)")
    ax.legend(frameon=False)
    plt.show()  
    
    if ret:
         return p_list, g_1, g_2

Vediamo quindi l'andamento delle energie libere (tra 10 e 22 GPa, in 12 punti, a T=1700 K)

In [269]:
plot_g_p(10,22,12,1700)

Quindi usiamo equilib tra T=1500 K e 2000 K (in 12 punti):

In [270]:
equilib(1500,2000,12,prod=['per',1,'cor',1], rea=['sp',1])

  T (K)  P (GPa)  DH(J/mol)  DS (J/mol K)  DV (J/bar)  Slope (bar/K)
 1500.0   12.86  -18043.761    -12.029       -0.247        48.73    
 1545.5   13.09  -18840.691    -12.191       -0.246        49.61    
 1590.9   13.31  -19636.057    -12.343       -0.245        50.46    
 1636.4   13.55  -20429.066    -12.484       -0.244        51.28    
 1681.8   13.78  -21218.984    -12.617       -0.242        52.07    
 1727.3   14.02  -22005.123    -12.740       -0.241        52.84    
 1772.7   14.26  -22786.844    -12.854       -0.240        53.59    
 1818.2   14.51  -23563.554    -12.960       -0.239        54.30    
 1863.6   14.75  -24334.703    -13.058       -0.237        55.00    
 1909.1   15.01  -25099.781    -13.148       -0.236        55.67    
 1954.5   15.26  -25858.318    -13.230       -0.235        56.33    
 2000.0   15.52  -26609.880    -13.305       -0.234        56.96    



Average Clapeyron Slope (from Delta S/Delta V):  53.07 bar/K
Clapeyron slope (from a linear fit of the P/T curve):  53.13 bar/K

Il risultato รจ in sostanziale accordo con le osservazioni sperimentali.