Calcoliamo l'entropia di stato standard (P=1 atm, T=298.15 K) e l'entropia a una data temperatura T (a pressione ambiente) a partire dai dati di calore specifico a pressione costante (misurati a partire da 20 fino a 500K).
Importiamo nel notebook anche la libreria inspect che consente di visualizzare il codice python delle funzioni contenute in entropy.py.
import inspect as ins
%matplotlib inline
%run entropy.py
I valori sperimentali del calore specifico sono salvati nella lista Cp_list: uno per ogni valore di temperatura della lista T_list
print("Temperature (K):\n", T_list)
print("\nCalore specifico (J/K mole):\n", Cp_list)
Volendo visualizzare meglio la lista dei valori T/Cp si può ricorrere alla libreria pandas (che è già importata all'interno di entropy.py con l'alias pd) si procede nel seguente modo:
serie=(T_list,Cp_list)
df=pd.DataFrame(serie, index=['Temp','Cp_exp'])
Si stampa il dataframe eliminando (opzionalmente) la numerazione delle colonne con header=False)
print(df.to_string(header=False))
Facciamo un fit del calore specifico in funzione della temperatura, usando la funzione fit; il fit è fatto sulla base di una polinomiale del tipo:
$$C_P(T)=aT+bT^{-1}+cT^2+dT^{-2}+eT^{0.5}+fT^{-0.5}$$I valori delle potenze sono già contenuti nella funzione Cp che calcola il calore specifico a una data T:
print(ins.getsource(Cp))
print(ins.getsource(fit))
il fit è fatto usando la funzione curve_fit della libreria scipy.optimize qui la documentazione generale. curve_fit vuole come argomenti
I valori ottimizzati dei coefficienti della polinomiale sono restituiti da curve_fit nella lista opt (la lista err contiene le corrispondenti deviazioni standard) e sono salvati nella variabile reg.par (istanza della classe fpar) usando il metodo .set
fit()
Usando il polinomio ottimizzato, ci si può calcolare il calore specifico a qualsiasi temperatura con la funzione Cp passandogli oltre al valore di temperatura, anche la lista dei coefficienti conservata in reg.par:
Cp(300,*reg.par)
Si noti l'uso dell'asterisco che precede il nome reg.par: specifica il passaggio alla funzione Cp di una lista di valori di lunghezza arbitraria.
Usiamo la funzione check_cp per visualizzare i risultati del fit. Il codice della funzione è:
print(ins.getsource(check_cp))
check_cp()
La funzione entropia calcola appunto l'entropia a una data temperatura, attraverso l'integrale
$$S(T)=\int_{0}^{T}\frac{C_P(T)}{T}\,\mathrm dT$$La funzione integranda $C_P/T$ è la funzione integrand del programma. La funzione entropia calcola l'integrale su scritto usando la funzione quad della libreria scipy.integrate (la documentazione la trovate qui).
La sintassi di quad vuole come input:
quad restituisce i risultati del calcolo in una lista di cui il primo elemento (quello con indice 0) è il valore dell'integrale. Quindi salviamo l'output di quad nella variabile ent; la funzione entropia restituisce il valore ent[0].
print(ins.getsource(integrand))
print(ins.getsource(entropia))
Per esempio, calcoliamo l'integrale:
$$I=\int_{0}^{5}x^2\,\mathrm dx$$Scriviamo la funzione integranda (myfunc) in modo che restituisca il valore di $x^{deg}$ dati $x$ e deg; quindi calcoliamo l'integrale e salviamo il risultato nella variabile integ:
def myfunc(x,d):
return x**d
deg=2
integ=quad(myfunc,0,5,args=deg)
integ è una lista di 2 valori: il valore dell'integrale e la stima dell'errore di integrazione; il valore dell'integrale è il primo elemento della lista (python indicizza le liste a partire dal valore 0 e non 1...)
print(integ)
print("\nValore dell'integrale %5.3f:" % integ[0])
Calcoliamo l'entropia a 400 K:
entropia(400)
Possiamo fare un plot dell'entropia in funzione della temperatura, usando la funzione plot_entropy, il cui codice è:
print(ins.getsource(plot_entropy))
plot_entropy(1000)
Tutte le funzioni viste sopra possono essere richiamate in sequenza usando la funzione start:
print(ins.getsource(start))
Come si vede, la funzione start esegue il fit del $C_P$, lo visualizza, fa un plot dell'entropia e calcola l'entropia alla temperatura standard:
start(600)