verticale

Misure sperimentali su un sistema microcogenerativo basato su motore Stirling

Lo scopo principale del progetto è stato quello di eseguire delle prove su di una unità Stirling cogenerativa (WhisperGen EU1) al variare delle condizione di pressione del fluido di lavoro (azoto) e confrontare i risultati sperimentali con il modello numerico del motore Stirling elaborato. Inoltre è stato sviluppato un modello numerico per il calcolo delle proprietà termodinamiche dell’acqua e dei principali composti dell’aria e dei combustibili in fase gassosa comunemente utilizzati. Il modello è in grado di trattare sia le singole specie, che una miscela costituita da esse (monofase allo stato gassoso o bifase con presenza di acqua liquida). I risultati sperimentali concordano con le previsioni del modello numerico dello Stirling, e l’utilizzo del modello numerico per le proprietà termodinamiche a portato ad una riduzione sull’errore della chiusura dei bilanci.

Scarica il PDF Scarica il PDF
Aggiungi ai preferiti Aggiungi ai preferiti


Articoli tecnico scientifici o articoli contenenti case history
Tesi di Laurea, Politecnico di Milano, Anno Accademico 2012-2013

Pubblicato
da Alessia De Giosa
VerticaleSegui aziendaSegui




Settori: 

Parole chiave: 


Estratto del testo
POLITECNICO DI MILANO Facoltà di Ingegneria Industriale Corso di Laurea in Ingegneria Meccanica




MISURE SPERIMENTALI SU UN SISTEMA MICROCOGENERATIVO BASATO SU MOTORE STIRLING AL VARIARE DELLE CONDIZIONI DEL FLUIDO DI LAVORO



Relatore: Ing. Gianluca VALENTI
Co-relatore: Ing. Nicola FERGNANI



Tesi di Laurea di:

Luca BANFI Matr. 781314


Anno Accademico 2012 - 2013 3 Ringraziamenti

Desidero innanzitutto ringraziare i miei genitori, Sergio e
Lucia, per avermi permesso di frequentare l''università senza
farmi mai mancare nulla, la mia fidanzata, Elena, per avermi
sostenuto durante tutti questi anni di studi, e mia sorella, Laura.

Inoltre desidero ringraziare tutte le persone che mi hanno
aiutato nel presente lavoro di tesi: l''Ing. Valenti, che mi ha
assistito durante tutte le fasi dello sviluppo del modello
numerico, l''Ing. Ravidà per il prezioso aiuto durante le fasi di
smontaggio e rimontaggio dello WhsiperGen, l''Ing. Fergnani
per il supporto sul modello numerico dell''unità Stirling, l''Ing.
Lazzari e l''Ing. Di Marcoberardino per il loro aiuto durante le
fasi sperimentali e di stesura del codice VBA, e il laureando
Sosio, anch''esso tesista nel Laboratorio di Micro
Cogenerazione, per l''ottimo rapporto di collaborazione che si è
venuto a creare.
4
5 Indice Generale




Introduzione

1
Motore Stirling Cogenerativo ............................................................. 13 1.1 Ciclo Stirling ............................................................................... 15 1.1.1 Ciclo termodinamico ideale ........................................... 16
1.1.2 Ciclo termodinamico reale ............................................. 21
1.1.3 Fluido di lavoro .............................................................. 23 1.2 Configurazione tecnologica ......................................................... 25 1.2.1 Configurazione Alpha .................................................... 26
1.2.2 Configurazione Beta ...................................................... 29
1.2.3 Configurazione Gamma ................................................. 33 1.3 WhisperGen EU1 ........................................................................ 34

2 Modello numerico per fluidi puri e miscele monofase / bifase .......... 39 2.1 Visual Basic for Application (VBA) ........................................... 40 2.2 Modulo ''Water' .......................................................................... 40 2.2.1 Regione 1 ....................................................................... 43
2.2.2 Regione 2 ....................................................................... 45
2.2.3 Regione 3 ....................................................................... 49
2.2.4 Regione 4 ....................................................................... 52
2.2.5 Regione 5 ....................................................................... 53
2.2.6 Struttura modulo ''Water' .............................................. 56
2.2.7 Energia libera di Gibbs .................................................. 57
2.2.8 Energia libera di Helmholtz ........................................... 59 2.3 Modulo ''Gas'.............................................................................. 60 2.3.1 Struttura modulo ''Gas' .................................................. 62 2.4 Modulo ''Mixture' ....................................................................... 63 2.4.1 VLE: Vapor-Liquid Equilibrium ................................... 64
2.4.2 Struttura modulo ''Mixture' ........................................... 66

3 Prove sperimentali al variare della pressione del fluido di
lavoro
...................................................................................................... 69 3.1 Breve cenni sul modello numerico del motore Stirling ............... 70 3.1.1 Modello adiabatico di Urieli .......................................... 71 3.1.2 Analisi semplificata di secondo ordine .......................... 75 6
3.2 Prove Sperimentali ...................................................................... 77 3.2.1 Prove in stato stazionario ............................................... 79 3.2.2 Procedura carico / scarico dell''azoto ............................. 81 3.2.3 Bilancio energetico ........................................................ 84 3.2.4 Risultati prove sperimentali ........................................... 88 3.3 Confronto tra prove sperimentali e modello numerico ............... 89 3.4 Propagazione degli errori ............................................................ 95


4 Conclusione e sviluppi futuri ............................................................... 97 4.1 Sviluppi futuri.............................................................................. 98

Bibliografia ..................................................................................................... 101


Allegato 1 ......................................................................................................... 103

Allegato 2 ......................................................................................................... 105

Allegato 3 ......................................................................................................... 107

Allegato 4 ......................................................................................................... 109

Allegato 5 ......................................................................................................... 111

7 Elenco delle Figure




1.1 Confronto tra produzione in cogenerazione e produzione separata ........ 14 1.2 Brevetto originale del motore Stirling del 1816 ..................................... 16 1.3 Diagramma p-v (pressione '' volume specifico) del ciclo Stirling ideale ....................................................................................................... 17 1.4 Diagramma T-s (Temperatura '' entropia) del ciclo Stirling ideale .... ....18 1.5 Confronto tra ciclo di Carnot e Stirling a pari differenza di volume ...... 20 1.6 Configurazione Alpha - Fase di riscaldamento isocoro. ......................... 26 1.7 Configurazione Alpha - Fase di espansione ............................................ 27 1.8 Configurazione Alpha - Fase di raffreddamento isocoro ........................ 27 1.9 Configurazione Alpha - Fase di compressione ....................................... 28 1.10 Configurazione Alpha a doppio effetto ................................................... 28 1.11 Configurazione Beta. ............................................................................... 29 1.12 Configurazione Beta - Fase di riscaldamento isocoro. ............................ 30 1.13 Configurazione Beta - Fase di espansione. ............................................. 31 1.14 Configurazione Beta - Fase di raffreddamento isocoro. ......................... 31 1.15 Configurazione Beta - Fase di compressione .......................................... 32 1.16 Configurazione Gamma .......................................................................... 33 1.17 Esempio di applicazione del WhisperGen EU1. ..................................... 34 1.18 Schema d''impianto WhisperGen EU1. ................................................... 35 1.19 Principali componenti dell''unità WhisperGen EU1. .............................. 37 1.20 Vista esterna dell''unità WhisperGen (a) e cinematismo Wobble Yoke (b). ................................................................................................. 37


8
2.1 Suddivisione in regioni della IAPWS ..................................................... 42 2.2 Suddivisione in sotto-regione per la Regione 3 ...................................... 50 2.3 Ingrandimento delle sotto-regioni da 3c a 3r .......................................... 51

3.1 Rendimento elettrico e termico riferito al PCS ....................................... 69 3.2 Schema semplificato del motore Stirling ................................................ 71 3.3 Schema di impostazione del modello adiabatico e profilo di temperatura all''interno del motore .......................................................... 72 3.4 Sistema di equazioni del modello ideale adiabatico ............................... 74 3.5 Schema di misurazione ........................................................................... 78 3.6 Andamento della temperatura dell''acqua di ritorno (8 bar) .................... 80 3.7 Manometro .............................................................................................. 82 3.8 Punto di ricarica dell''unità Stirling ......................................................... 83 3.9 Flussi energetici dell''unità WhisperGen EU1 ........................................ 90 3.10 Confronto potenze prove sperimentali - modello ................................... 92 3.11 Confronto potenza elettrica prove sperimentali - modello...................... 93 3.12 Calore entrante nel ciclo Stirling ............................................................ 93 3.13 Confronto rendimento elettrico ............................................................... 94

9 Elenco delle Tabelle




1.1 Dati indicativi per diversi fluidi di lavoro. .............................................. 24 1.2 Caratteristiche nominali del motore WhisperGen EU1 .......................... 36

2.1 Valori numerici per i coefficienti delle equazioni (2.1) e (2.2) .............. 43 2.2 Relazioni proprietà termodinamiche '' Regione 1................................... 44 2.3 Derivate dell''energia libera di Gibbs ...................................................... 45 2.4 Relazioni proprietà termodinamiche '' Regione 2................................... 46 2.5 Derivate della parte di gas ideale dell''energia libera di Gibbs ............... 47 2.6 Derivate della parte residua dell''energia libera di Gibbs ........................ 48 2.7 Derivate della parte di gas ideale dell''energia libera di Gibbs ............... 54 2.8 Derivate della parte residua dell''energia libera di Gibbs ........................ 55 2.9 Relazioni proprietà termodinamiche '' Regione 5................................... 56 2.10 Unità di misura VBA .............................................................................. 57

3.1 Coefficienti correttivi Bi per il misuratore di portata Bronkhorst ........... 85 3.2 Risultati prove sperimentali .................................................................... 88 3.3 Confronto prove sperimentali '' modello numerico ................................ 91
4.1 Principali caratteristiche scambiatore LU-VE ........................................ 99 10
11 Sommario

La presente Tesi di Laurea verte sull''attività svolta presso il Laboratorio di
Micro Cogenerazione (LMC) del Politecnico di Milano. Lo scopo
principale è stato quello di eseguire delle prove su di una unità Stirling
cogenerativa (WhisperGen EU1) al variare delle condizione di pressione
del fluido di lavoro (azoto) e confrontare i risultati sperimentali con il
modello numerico del motore Stirling elaborato. Inoltre è stato sviluppato
un modello numerico per il calcolo delle proprietà termodinamiche
dell''acqua e dei principali composti dell''aria e dei combustibili in fase
gassosa comunemente utilizzati. Il modello è in grado di trattare sia le
singole specie, che una miscela costituita da esse (monofase allo stato
gassoso o bifase con presenza di acqua liquida). I risultati sperimentali
concordano con le previsioni del modello numerico dello Stirling, e
l''utilizzo del modello numerico per le proprietà termodinamiche a portato
ad una riduzione sull''errore della chiusura dei bilanci. ' stata inoltre
sviluppata una procedura interna per le operazione di carica / scarica del
fluido di lavoro all''interno dell''unità Stirling.

Parole chiave: Motore Stirling, Cogenerazione, Modello numerico

Abstract

This Thesis focuses on activities carried out at the Laboratory of Micro
Cogeneration (LMC) of the Polytechnic of Milan. The main purpose has
been to perform the tests on a Stirling cogeneration unit (WhisperGen
EU1) changing working fluid pressure (nitrogen) and compare the
experimental results with the numerical model of the Stirling engine
developed at the laboratory. Furthermore, a numerical model has been
developed for the calculation of the thermodynamic properties of the water
and of the main compounds of the air and commonly used gas fuel. The
model is able to treat both the individual species, that a mixture consisting
of these (single-phase or two-phase gaseous state with the presence of
liquid water). Experimental results are in agreement with the numerical
model for Stilring engine, and the numerical model for thermodynamic
properties has allowed a reduction of the energy balance error. Furthermore
has been developed a procedure for the charge / discharge of working fluid
inside the Stirling engine.

Keywords: Stirling engine, Cogeneration, Numerical Model 12
13 Capitolo 1

Motore Stirling cogenerativo

Con il termine cogenerazione s''intende la produzione combinata di energia
elettrica/meccanica e di energia termica (calore) ottenute in appositi
impianti utilizzanti la stessa energia primaria. Per produrre la sola energia
elettrica si utilizzano generalmente centrali termoelettriche che disperdono
parte dell'energia primaria nell'ambiente: questa è energia termica di scarso
valore termodinamico essendo a bassa temperatura. Invece, per produrre la
sola energia termica, tradizionalmente si usano delle caldaie che
convertono l'energia primaria contenuta nei combustibili, di elevato valore
termodinamico, in energia termica di ridotto valore termodinamico. Se
un'utenza richiede energia elettrica ed energia termica, anziché installare
una caldaia e acquistare energia elettrica dalla rete, si può pensare di
realizzare un ciclo termodinamico per produrre energia elettrica sfruttando
i livelli termici più alti, e cedendo il calore residuo a più bassa temperatura
per soddisfare le esigenze termiche. Da questo punto di vista la
cogenerazione permette un notevole risparmio energetico, che però non è
scontato: si tratta allora di valutare quando è davvero vantaggiosa e rispetto
a quale alternativa. L'obiettivo fondamentale che si vuole perseguire con la
cogenerazione è quello di sfruttare al meglio l'energia contenuta nel
combustibile: a ciò consegue un minor consumo di combustibile e di
conseguenza anche un minor impatto ambientale.

La cogenerazione, quindi, permette di risparmiare energia e assicura
benefici oggettivi, misurabili e quantificabili. Su tale principio guida si
basa anche la trigenerazione, cioè la produzione simultanea di energia
termica, elettrica e frigorifera da un'unica fonte energetica. Cogenerazione
e trigenerazione rientrano nelle scelte strategiche delle aziende che vedono
nell'efficienza energetica un'opportunità essenziale per ridurre i costi e
aumentare la loro competitività. Altrettanto significativi sono i vantaggi a
livello di impatto ambientale, in quanto vengono ridotte drasticamente le
emissioni di CO2 grazie al minor consumo di combustibile fossile. In
sintesi, rispetto alla produzione separata delle stesse quantità di energia
elettrica e calore, la produzione combinata, se efficace, comporta:
' un risparmio economico conseguente al minor consumo di combustibile; Capitolo 1
14
' una riduzione dell'impatto ambientale, conseguente sia alla riduzione delle emissioni sia al minor rilascio di calore residuo nell'ambiente
(minor inquinamento atmosferico e minor inquinamento termico); ' minori perdite di trasmissione e distribuzione per il sistema elettrico nazionale, conseguenti alla localizzazione degli impianti in prossimità
dei bacini di utenza o all'autoconsumo dell'energia prodotta; ' la sostituzione di modalità di fornitura del calore meno efficienti e più inquinanti (caldaie, caratterizzate da bassi livelli di efficienza, elevato
impatto ambientale e scarsa flessibilità).
Per chiarire il significato di risparmio energetico connesso ad un impianto
cogenerativo rispetto alla produzione separata delle medesime quantità di
energia utile, si veda l'esempio riportato nella Figura 1.1. Supponendo che
un impianto cogenerativo, per produrre 35 unità di energia elettrica e 50
unità di calore utile, consumi 100 unità di combustibile, il rendimento
termodinamico complessivo di conversione, inteso come rapporto tra
l'energia utile prodotta e l'energia primaria del combustibile utilizzato,
risulta dell'85%. Se si considera invece il caso di produzione separata,
supponendo di produrre 35 unità di energia elettrica con una centrale
termoelettrica avente un rendimento elettrico di circa il 40% e 50 unità di
calore utile con una caldaia avente un rendimento termico pari a circa 80%,
si avrebbe un consumo di combustibile pari a 140 unità di combustibile.
Nel caso di produzione separata delle stesse quantità di energia elettrica e
calore, risulterebbe quindi un consumo di 140 unità di combustibile
anziché le 100 richieste dall'impianto di cogenerazione. Il risparmio di
energia primaria conseguibile con la cogenerazione è dunque pari al 28%.
Figura 1.1 Confronto tra produzione in cogenerazione e produzione separata [1]. Motore Stirling cogenerativo 15 1.1 Ciclo Stirling

Il motore Stirling fu brevettato da Robert Stirling, sacerdote scozzese, nel
1816 come evoluzione dei ''motori ad aria calda', con l''obiettivo di
competere con i più diffusi motori a vapore inventati durante la prima
rivoluzione industriale.

Recentemente si ha avuto maggiore interesse verso questo tipo di
tecnologia, e ciò è dovuto alle sue caratteristiche favorevoli in termini di:
' potenziali bassi consumi di combustibili (elevati rendimenti); ' possibilità di lavorare con diverse fonti di combustibile, in quanto la fonte di energia può essere di qualsiasi tipo purché ad un sufficiente
livello di temperatura. Sono stati progettati motori Stirling alimentati
da energia solare così come da una varietà di combustibili liquidi e
solidi; ' ragionevole potenza specifica; ' basse emissioni;
' semplicità meccanica; ' il più alto lavoro specifico in uscita rispetto ogni altro ciclo chiuso; ' silenziosità e basse vibrazione in esercizio; ' curva di coppia favorevole per applicazioni in ambito automotive;
Nonostante queste caratteristiche favorevoli, il motore Stirling, non ha
avuto molto successo all''epoca, principalmente a causa di limitazioni
tecnologiche, soprattutto in ambito metallurgico (lo sviluppo del primo
motore Stirling fu ostacolato dal problema del surriscaldamento e
successiva ossidazione dei materiali utilizzati), e di difficoltà nel prevedere
con un''adeguata precisione le sue prestazioni. Molti aspetti del motore
Stirling sono unici comparati con i cicli termodinamici del tempo, e ne
fanno una delle più incredibili invenzioni del suo genere, di molti anni in
avanti rispetto le conoscenze tecnologiche del tempo. Per esempio non si fa
alcun uso di valvole, opera secondo un ciclo termodinamico chiuso, in cui
il fluido di lavoro non viene a contatto con l''ambiente esterno permettendo
di lavorare ad una pressione minima del ciclo al di sotto di quella
atmosferica. Un ulteriore vantaggio è dato dalla combustione continua,
infatti, in questo caso, è molto più semplice far avvenire un processo di
combustione pulito rispetto ad una combustione intermittente. Ma
l''innovazione più importante, e che sostanzialmente determina un
incremento dell''efficienza del ciclo termodinamico, fu l''introduzione di un
rigeneratore, avente il compito di ''immagazzinare' calore subito prima che
il fluido raggiungesse il pozzo freddo, per poi restituirlo nelle fasi Capitolo 1
16
successive, ovvero effettuando un pre-riscaldo del fluido di lavoro. Si
capisce fin da subito come questa soluzione permetta un notevole
miglioramento dell''efficienza del ciclo termodinamico rispetto ad un
generico motore tradizionale in cui i gas combusti, ad elevata temperatura,
vengono scaricati direttamente nell''ambiente.

Le conoscenze moderna, e l''uso dei computer, hanno fatto in modo che la
tecnologia del motore Stirling possa ora essere di largo interesse,
soprattutto in vista della sua elevata efficienza teorica raggiungibile in
assetto cogenerativo per basse potenze.

Figura 1.2 Brevetto originale del motore Stirling del 1816 [2].

1.1.1 Ciclo Termodinamico ideale
Come tutti i motori termici, il motore Stirling si basa sulle quattro
trasformazioni ideali di compressione, riscaldamento, espansione e
raffreddamento. Il ciclo ideale Stirling, è quindi così composto:
' una compressione isoterma, con cessione di una quantità di calore al pozzo freddo, equivalente al lavoro di compressione compiuto (1-2); Motore Stirling cogenerativo 17 ' un riscaldamento isocore del fluido di lavoro per opera del rigeneratore, che cede il calore precedentemente ricevuto nella
precedente isocora (2-3); ' un''espansione isoterma del fluido di lavoro, con introduzione di una quantità di calore dall''esterno nel pozzo caldo del ciclo, equivalente al
lavoro di espansione compiuto (3-4); ' una cessione del calore isocora dal fluido di lavoro ad un generatore ad accumulo (4-1);
Nel ciclo Stirling ideale gli unici scambi di calore con l''esterno sono
isotermi e avvengono alle temperature delle sorgenti. Invece i passaggi
dalla temperatura massima alla minima (e viceversa) avvengono a volume
costante; ciò è possibile grazie al fatto che, quando il fluido di lavoro si
raffredda, cede calore ad una matrice solida porosa, il rigeneratore. Questo
calore viene poi restituito dal rigeneratore al fluido stesso durante la fase di
riscaldamento.
Figura 1.3 Diagramma p-v (pressione '' vol. specifico) del ciclo Stirling ideale [3]. Capitolo 1
18
Figura 1.4 Diagramma T-s (Temperatura '' entropia) del ciclo Stirling ideale [3].

Sopra sono riportati i diagrammi termodinamici relativi al ciclo Stirling
ideale (Figura 1.3. e Figura 1.4.), da cui è possibile valutare
rispettivamente il lavoro specifico prodotto dal ciclo (equivalente all''area
interna al ciclo diagramma p-v), ed il calore assorbito (equivalente all''area
sottesa alle trasformazioni ad entropia crescente nel diagramma T-s). Tale
ciclo è comunemente detto ''ciclo quadrato', per la sua particolare forma
nei diagrammi p-v (pressione '' volume specifico) e T-s (temperatura ''
entropia).

Come è noto, il ciclo di Carnot ideale è il ciclo termodinamico che
permette di ottenere il miglior rendimento data una temperatura massima
ed una temperatura minima. Nel 1873 Reitlinger dimostrò che è possibile
raggiungere il medesimo rendimento del ciclo di Carnot in tutte quelle
macchine che si basano su un ciclo termodinamico composto da due
isoterme e altre due trasformazioni omologhe reversibili, operanti nel
medesimo intervallo di temperature. Quindi se si parte da un ciclo Motore Stirling cogenerativo 19 costituito da due isoterme e si realizzano le altre due trasformazioni
mediante isocore, politropiche o isobare reversibili, si otterranno tre cicli
termodinamici con lo stesso rendimento ideale.

I tre cicli risultanti, quindi, sono i seguenti:
' Il ciclo Stirling, caratterizzato da due isoterme e due isocore; ' Il ciclo Ericsson, caratterizzato da due isoterme e due isobare; ' Il ciclo di Reitlinger, caratterizzato da due isoterme e due politropiche;
Il rendimento di un ciclo Stirling ideale è quindi così definito:


(1.1)
nella quale compaiono le quantità di calore scambiate con l''esterno Q34 = Qin e Q12 = Qout, che possono essere calcolati come: ' '
(1.2)
' '
(1.3)
Sostituendo quindi le espressioni (1.2) e (1.3) nella (1.1), si ottiene:



(1.4)
poiché

, come si vede dalla figura 1.3. Ovvero il rendimento del ciclo Stirling è pari al rendimento di Carnot, ovviamente per lo stesso
rapporto tra la temperatura massima Tmax e minima Tmin.

(1.5) Il calore che viene assorbito e poi ceduto dal ciclo durante le
trasformazioni isocore, ovvero il calore di rigenerazione, è pari a: Capitolo 1
20
( ) (1.6)
Dove si è indicato con ''Cv' il calore specifico del fluido a volume costante, e con ''R' la costante universale dei gas.

Poiché le trasformazioni isocore sono rigenerative, per esse si ha che la
somma algebrica della quantità di calore assorbita e ceduta è nulla. Il
lavoro complessivamente sviluppato in un ciclo motore ideale, come già
indicato nella formula (1.1), si ottiene per differenza tra l''energia termica
in ingresso e in uscita.


( ) (1.7)
Si confronta ora il ciclo Stirling ideale (1,2,3,4) con quello di Carnot ideale
(1'',2'',3,4). Come si vede dalla Fig. 1.4., il ciclo Stirling ideale presenta due
trasformazioni isocore al posto delle due adiabatiche che caratterizzano il
ciclo di Carnot ideale, quindi a pari differenza di volumi, l''area del ciclo
aumenta sensibilmente. In questo modo si ottiene un lavoro ragionevole,
dal ciclo Stirling, senza che le pressioni e i volumi siano troppo grandi,
cosa che invece avviene nel caso del ciclo di Carnot.
Figura 1.5 Confronto tra ciclo di Carnot e Stirling a pari differenza di volume [3]. Motore Stirling cogenerativo 21 1.1.2 Ciclo Termodinamico reale
Durante la precedente trattazione del ciclo Stirling ideale, sono state
considerate delle importanti ipotesi di idealità, tuttavia è impossibile che le
macchine reali possano operare seguendo questo tipo di ciclo
termodinamico. La prima ipotesi su cui si basa il ciclo Stirling ideale è che
le trasformazioni siano reversibili: viene trascurato ogni attrito
fluidodinamico o meccanico, e il fluido si trova, in ogni istante, in una
condizione di equilibrio per quanto riguarda la pressione e la temperatura.
Ciò significa che le trasformazioni sono quasi statiche e non vi è scambio
di calore trasversale nel fluido. Inoltre le equazioni sono state ricavate
considerando il fluido completamente contenuto in uno dei due cilindri, di
compressione o espansione, ipotizzando quindi che non ve ne sia nel
rigeneratore, nelle zone di connessione dei cilindri stessi. Infine si è
considerato un rigeneratore perfetto, caratterizzato da una capacità termica
infinita.

L''efficienza reale viene valutata mediante un coefficiente detto efficienza
relativa ηrel che è il rapporto tra l''efficienza reale e quella ideale dello Stirling, e che varia da valori di 0,4 a 0,7:

(1.8)
Le cause di scostamento dall''idealità vengono classificate di seguito:
' Perdite legate allo scambio termico: Compressione ed espansione non isoterme
Data l''impossibilità pratica di realizzare un moto infinitamente lento
da approssimare una trasformazione isoterma, le macchine reali
operano secondo trasformazioni adiabatiche. La frequenza di lavoro
molto elevata (migliaia di cicli al minuto), la rapidità con cui il fluido
evolve all''interno del cilindro, non permettono al gas di mantenere
costante la sua temperatura nelle fasi di scambio termico durante
l''espansione e la compressione, pertanto piuttosto che trasformazioni
isoterme si riscontra un comportamento più simile a trasformazioni
adiabatiche.

Perdite termiche per convenzione, conduzione ed irraggiamento
Poiché le superfici di scambio termico durante le fasi di introduzione e
rimozione del calore sono limitate, la temperatura del gas nel cilindro
non eguaglia i valori di temperatura del pozzo caldo e del pozzo Capitolo 1
22
freddo, come invece succederebbe in un caso ideale. Questo significa
che la stima del rendimento eseguita come dalla formula (1.4) e (1.5),
in realtà è una sovrastima, poiché tali temperature non sono mai
raggiunte dal fluido di lavoro. L''efficienza è inoltre ridotta da
fenomeni naturali dissipativi con l''ambiente esterno: il fluido di lavoro
riscaldato cede parte del calore per convenzione alle pareti metalliche
dell''involucro, da cui per conduzione attraversa l''intero spessore delle
pareti e quindi per convezione ed irraggiamento avviene la dispersione
in ambiente.

Perdite termiche nei gas combusti
La perdita è associata al contenuto entalpico dei gas combusti allo
scarico, e può essere ridotta adottando sistemi di recupero termico per
il preriscaldo dell''aria comburente. Nei motori cogenerativi questo
calore viene in parte recuperato tramite la potenza termica utile.

Rigenerazione imperfetta
Per rispettare il modello ideale, gli scambi termici tra fluido di lavoro e
rigeneratore dovrebbero essere istantanei, ma questo sarebbe possibile
solo se il fluido di lavoro avesse capacità termica nulla e viceversa il
rigeneratore capacità termica infinita. Nella realtà, durante la fase di
riscaldamento del fluido di lavoro, il rigeneratore non cede tutta
l''energia termica accumulata durante la precedente fase di
raffreddamento.
' Perdite fluidodinamiche:
Presenza dei volumi morti
La realizzazione del ciclo ideale comporta contenere il fluido di
lavoro, nelle fasi di espansione e compressione, interamente nei
rispettivi volumi di espansione e compressione, ed ottenere quindi un
riscaldamento e raffreddamento omogeno del fluido. Ciò comporta che
la corsa dei pistoni sia tale da annullare il volume all''interno del
cilindro in corrispondenza del punto morto superiore, e disporre di
elementi quali condotti di collegamento fra le camere, riscaldatore,
refrigeratore e rigeneratore di volume nullo. Poiché tutto questo non è
fisicamente attuabile, la presenza dei volumi morti non può essere
evitata.

Perdite di carico e trafilamenti
Questo tipo di perdite sono essenzialmente attribuibili alla viscosità
del fluido di lavoro, e quindi ad attriti di tipo viscoso tra pareti del
cilindro e fluido di lavoro, alle quali si aggiungono le perdite dovute ai Motore Stirling cogenerativo 23 trasferimenti del fluido lungo i condotti. Occorre tenere in
considerazione anche le perdite concentrate dovute ad irregolarità del
circuito percorso dal fluido di lavoro, come variazioni di sezione o di
direzione.
' Perdite dovute al cinematismo:
Il cinematismo, qualunque esso sia, non è in grado di eseguire
perfettamente le trasformazioni isocore, possono solo essere
approssimate.
' Perdite meccaniche:
In generale con esse sono intese le dissipazioni energetiche dovute al
moto relativo di superfici a contatto tra loro, come ad esempio la forza
di attrito sui cuscinetti o sulle pareti del pistone.
' Perdite elettriche:
Dovute alla conversione di energia meccanica in energia elettrica.

1.1.3 Fluido di lavoro
Il fluido di lavoro, sottoposto alle sollecitazioni termiche, dà origine al
movimento degli organi mobili del motore, producendo quindi lavoro
meccanico. I primi motori Stirling utilizzavano aria come fluido di lavoro,
favorita soprattutto per la sua elevata disponibilità e il basso costo, e le
condizioni operative non si discostavano molte dalla pressione atmosferica.
La presenza di aria e di altri liquidi lubrificanti derivanti da idrocarburi,
poteva però dare origine a miscele esplosive. In seguito sono stati adottati
gas riducenti come l''idrogeno, che tuttavia non eliminavano il problema
delle miscele esplosive, o gas inerti come elio e azoto, più sicuri ma più
costosi e meno accessibili.

In generale il fluido di lavoro di un ciclo termodinamico dovrebbe
possedere le seguenti proprietà:
1. Buon coefficiente di scambio termico.
2. Elevata capacità termica.
3. Stabilità alle alte temperature
4. Bassa viscosità. Capitolo 1
24
Inoltre non è da trascurare la disponibilità del fluido, il suo costo, la
capacità di stoccaggio e la sicurezza durante l''impiego. Utile ai fini di una
selezione preliminare del fluido di lavoro è un indice definito da Martini e
Clarke come ''capability factor' [16], che consiste nel rapporto tra la
conduttività termica ed il prodotto della capacità termica per la densità.
Nella pratica si evidenzia come siano preferibili fluidi a molecola leggera,
come idrogeno o elio, pe gli elevati coefficienti di scambio termico,
tuttavia essi, proprio a causa delle piccole dimensioni delle loro molecole,
sono difficilmente contenibili dalle tenute, e inoltre l''idrogeno è altamente
infiammabile. Fluidi a molecola più pesante, come aria e azoto, sono meno
adatti per il minor coefficiente di scambio termico, inoltre l''ossigeno
contenuto nell''aria causa fenomeni di corrosione alle alte temperature,
tuttavia è più semplice realizzare le tenute. Per individuare il fluido di
lavoro migliore occorrerebbe analizzare le prestazioni dell''intero sistema
variando di volta in volta il fluido di lavoro, tuttavia un''indagine
sperimentale di questo tipo è troppo onerosa sia in termini di tempo sia in
termini di denaro. Walker suggerisce un approccio più semplice basato su
un''analisi del flusso in regime stazionario: basandosi sull''analogia di
Reynolds, si ricava una relazione tra calore scambiato e forze di attrito che
risulta essere:

'' (1.9)
Confrontando tra loro vari potenziali fluidi di lavoro sulla base del
capability factor e del calore trasferito al fluido di lavoro, si nota che
l''utilizzo dell''idrogeno e dell''elio non è pienamente giustificato, ma che
per esempio si ottengono prestazioni migliori usando una miscela di sodio
e potassio.


Fluido di lavoro Calore trasferito Capability Factor Aria 1.00 1.00 Elio 1.42 0.83 Idrogeno 3.42 0.68 Acqua 1.95 0.39 Sodio-Potassio 32.62 1.32 Tabella 1.1 Dati indicativi per diversi fluidi di lavoro. Motore Stirling cogenerativo 25 Ovviamente i due indici devono essere calcolati con riferimento alle
condizioni prevalenti in cui si trova ad operare il fluido di lavoro, poiché le
proprietà coinvolte sono variabili con pressione e temperatura.
1.2 Configurazione tecnologica

Il motore Stirling, nel corso degli anni, è stato interpretato in numerose
varianti dal punto di vista cinematico e fluidodinamico, mantenendo però
la presenza di una zona calda e di una zona fredda di scambio termico. Il
motore Stirling realizza un''oscillazione ciclica del fluido di lavoro che in
parte viene trasformata in energia meccanica; le oscillazioni del fluido
vengono trasferite a masse connesse a generatori elettrici lineari o
cinematismi (come il classico biella-manovella) preposti alla conversione
del moto lineare dei pistoni in un moto rotazionale che aziona un
generatore elettrico.

' possibile classificare i motori Stirling secondo diversi criteri:
' Il modo in cui il motore scambia il lavoro con l''esterno e sposta il fluido: ' Tramite un pistone, sottoposto ad una importante differenza di pressione, che richiede la presenza di tenute per impedire il
trafilamento; ' Tramite un dislocatore, avente facce sottoposte alla stessa pressione, che sposta il fluido; ' Il numero di effetti: ' Motori a singolo effetto, nei quali i volumi di espansione e compressione dei diversi cilindri sono separati tra loro e sono
presenti uno o due pistoni e un dislocatore; ' Motori a doppio effetto, nei quali lo spazio di espansione di ogni cilindro è connesso tramite uno scambiatore di calore con
lo spazio di compressione del cilindro adiacente; ' La disposizione dei pistoni e dei cilindri: ' La configurazione ''Alpha' prevede che i pistoni siano in cilindri separati; ' La configurazione ''Beta' è caratterizzata da un sistema di dislocatore-pistone con camere di movimento unite; ' La configurazione ''Gamma' presenta invece un sistema dislocatore-pistone con camere di movimento separate. Capitolo 1
26
1.2.1 Configurazione Alpha
La configurazione Alpha è la più semplice dal punto di vista concettuale:
essa prevede l''utilizzo di due pistoni separati, il pistone motore e quello di
spostamento, collocati nei rispettivi cilindri mantenuti a diverse
temperature. Il pistone motore è a contatto con lo scambiatore di calore ad
alta temperatura, mentre il pistone di spostamento è a contatto con quello a
bassa temperatura. I due pistoni sono interconnessi tra loro mediante un
sistema biella-manovella e albero motore con volano. Questa
configurazione permette di ottenere un alto rapporto tra potenza e volume,
ma presenta numerosi problemi tecnici dovuti alla necessità di utilizzare
guarnizioni per contenere il gas. L''alta temperatura di lavoro, infatti, è
causa della scarsa durata delle guarnizioni stesse.

Analizzando il funzionamento del motore Stirling in configurazione Alpha,
si possono distinguere quattro fasi:
' Riscaldamento isocoro; ' Espansione; ' Raffreddamento isocoro; ' Compressione.
Dalle seguenti figure si nota che:
' Riscaldamento isocoro: il pistone di spostamento spinge il gas verso il pistone motore, che si sposta di conseguenza;
Figura.1.6 Configurazione Alpha - Fase di riscaldamento isocoro [4].
' Espansione: il gas scaldandosi si espande, muove il pistone motore e compie lavoro. Successivamente il gas ritorna verso il pistone di Motore Stirling cogenerativo 27 spostamento. I due sistemi biella-manovella fanno sì che il pistone
resti quasi fermo e che quello di spostamento si muova.
Figura 1.7 Configurazione Alpha - Fase di espansione [4].
' Raffreddamento isocoro: il movimento del pistone di spostamento permette l''entrata nel cilindro freddo del gas caldo, che viene a
contatto con il dissipatore di calore raffreddandosi. Nel frattempo il
pistone motore si muove verso il punto morto superiore. Anche in
questo caso è lo sfasamento tra i due sistemi biella-manovella a
regolare il moto relativo tra i due pistoni.
Figura 1.8 Configurazione Alpha - Fase di raffreddamento isocoro [4].
' Compressione: mentre il pistone motore si trova al punto morto superiore, l''inerzia accumulata dal volano fa sì che il sistema biella
manovella muova il pistone di spostamento così da spingere il gas
verso il pistone motore. Il ciclo così ricomincia. Capitolo 1
28
Il volano presente sull''asse del motore accumula energia durante le fasi,
rilasciandola poi tra una fase e la successiva. In questo modo si passa da
uno stadio all''altro senza l''arresto del motore quando i pistoni sono a fine
corsa, rendendo il movimento più omogeneo. Il motore di tipo Alpha può
anche essere utilizzato nella configurazione a cilindri multipli per ottenere
potenze elevate, generalmente i pistoni sono disposti a cerchio e si
muovono con uno sfasamento prestabilito mediante un disco rotante
inclinato.
Figura 1.9 Configurazione Alpha - Fase di compressione [4].
La configurazione Alpha può essere adottata anche in caso di disposizione
di più cilindri in serie in modo che la parte inferiore di ognuno, che funge
da volume di compressione, sia interconnessa alla camera di espansione del
cilindro adiacente. Tale configurazione è denominata Alpha a doppio
effetto
.
Figura 1.10 Configurazione Alpha a doppio effetto [5]. Motore Stirling cogenerativo 29 Il doppio effetto consente un notevole incremento della potenza specifica.
Altri vantaggi sono la riduzione delle perdite di trafilamento (resta
un''unica tenuta verso l''ambiente esterno ed è l''anello che premette lo
scorrimento delle bielle connesse a ciascun pistone) e la maggiore
compattezza. Inoltre, la tipologia di manovellismo consente un''elevata
continuità nel moto con tre o quattro pistoni sfasati rispettivamente di 120°
o 90°, con ridotti effetti inerziali e una buona efficienza meccanica: il
lavoro di compressione è sottratto direttamente al lavoro di espansione in
modo unidirezionale. Per contro gli aspetti critici della configurazione
Alpha a doppio effetto consistono nel maggior numero di cilindri, pistoni,
scambiatori e tenute rispetto alla configurazione Beta (tutti fattori che
contribuiscono ad aumentarne il prezzo). ' anche necessario un
posizionamento verticale del motore per consentire una corretta
lubrificazione del cinematismo e ridurre gli attriti.


1.2.2 Configurazione Beta
La configurazione beta è la classica configurazione del motore Stirling
descritta nel brevetto del 1816. Questo tipo di motore è composto dal
pistone dislocatore (displacer), la cui funzione è trasferire il fluido di
lavoro dal volume di espansione al volume di compressione, e dal pistone
motore contenuti entrambi nello stesso cilindro ed entrambi collegati,
mediante un manovellismo, allo stesso albero motore dotato di volano.

Figura 1.11 Configurazione Beta [5].
La presenza del cinematismo biella-manovella fa sì che i movimenti
compiuti dal pistone seguano i movimenti del dislocatore con un ritardo di
fase di 90° esatti. Il dislocatore non è a tenuta, quindi non fornisce alcuna
potenza aggiuntiva al motore, esso serve esclusivamente a muovere il gas Capitolo 1
30
di lavoro a volume costante. Il pistone motore, dotato di guarnizione,
comprime ed espande periodicamente l''aria nel cilindro. Il gas, per passare
dalla zona di espansione a quella di compressione, attraversa la serie
riscaldatore-rigeneratore-refrigeratore: quando il fluido di lavoro è
nell''estremità calda, si espande spingendo il pistone motore, quando invece
è nell''estremità fredda si contrae. Il momento motore della macchina ed il
volano poi, spingono il pistone motore nel senso inverso per comprimere il
gas.

Il ciclo del motore in configurazione Beta segue le quattro fasi già illustrate
per la configurazione Alpha.
' Riscaldamento isocoro: grazie al moto del dislocatore il fluido di lavoro viene trasferito dalla zona inferiore del cilindro, a contatto con
il refrigeratore, alla zona superiore, prossima allo scambiatore caldo. Il
pistone rimane fermo cosicché il volume totale a disposizione del
fluido rimanga costante e si realizzi la trasformazione isocora. Il
riscaldamento avviene prevalentemente per cessione di calore dal
rigeneratore al fluido e, per la restante parte, ad opera dello
scambiatore caldo fino al raggiungimento della temperatura massima
del fluido.

Figura 1.12 Configurazione Beta - Fase di riscaldamento isocoro [6].
' Espansione: il gas è trasferito totalmente nel volume di espansione dove, dopo essere stato riscaldato, si espande spingendo il pistone
motore verso il punto morto inferiore. Contemporaneamente al Motore Stirling cogenerativo 31 pistone, anche il displacer si sposta verso il basso per far si che il
volume, lato scambiatore freddo, resti costante e non vi sia passaggio
di flusso attraverso il rigeneratore.

Figura 1.13 Configurazione Beta - Fase di espansione [6].
' Raffreddamento isocoro: la maggior parte del gas viene trasferita nel volume di compressione dalla parziale salita del dislocatore. Il
raffreddamento avviene prevalentemente per cessione di calore al
rigeneratore ed è completato ad opera dello scambiatore freddo. In
questa fase si registra una consistente diminuzione di pressione.
Figura 1.14 Configurazione Beta - Fase di raffreddamento isocoro [6]. Capitolo 1
32
' Compressione: si verifica la risalita del pistone che. Riducendo il volume della camera collocata nella zona inferiore del cilindro,
comprime il fluido di lavoro raffreddato nella fase precedente.
Idealmente il displacer rimane fisso al punto morto superiore per fare
in modo che non vi sia passaggio di fluido verso la zona di espansione
a contatto con il riscaldatore.

Figura 1.15 Configurazione Beta - Fase di compressione [6]. La fase di espansione, durante la quale il pistone motore realizza la corsa
verso il punto morto inferiore, è quella in cui è compiuto lavoro. Si nota
che le fasi di riscaldamento e raffreddamento possono essere ritenute
isocore solo nell''ambito dell''analisi del ciclo ideale. In realtà durante la
fase di riscaldamento si verifica via via l''espansione del gas che, grazie
all''azione del dislocatore, passa nel condotto esterno di travaso e spinge
verso il basso il pistone. Analogamente, durante la fase di raffreddamento
si verifica la contrazione del gas che agevola la risalita del pistone.

Tra i vantaggi del motore Stirling in configurazione beta vi è la semplicità
costruttiva e la compattezza del sistema. Le tenute possono essere di
qualità inferiore e sono concentrate esclusivamente sul pistone. I volumi
morti sono piuttosto contenuti e ciò consente benefici in termini di
efficienza. Inoltre la compattezza e la possibilità di pressione medie
elevate, fanno sì che la potenza specifica sia elevata. Questa
configurazione, insieme a quella Alpha a doppio effetto, risulta la più
promettente per applicazioni di micro-cogenerazione. Gli aspetti di
maggiore criticità riguardano la modellazione termodinamica e meccanica.
Motore Stirling cogenerativo 33 1.2.3 Configurazione Gamma
Si tratta di una configurazione simile alla configurazione beta, la differenza
sta nel fatto che il pistone motore è contenuto in un cilindro separato da
quello del pistone dislocatore, nonostante i due pistoni siano collegati alla
stessa trasmissione. Tale soluzione si pone come evoluzione della
configurazione Beta, nella quale durante la fase di espansione si verifica
che parte dell''espansione ha luogo nello spazio di compressione portando a
una riduzione della potenza specifica. Tuttavia, ciò dipende dal moto non
ideale del displacer e avviene anche nella configurazione Gamma con
cinematismo biella-manovella.

Si nota come il volume di compressione sia diviso tra i due cilindri che
costituiscono la macchina: è possibile scegliere di eseguire il
raffreddamento solo sul cilindro di compressione come anche suddividerlo
su entrambi i cilindri. Il vantaggio è la semplicità meccanica che rende tale
soluzione adottabile nei motori a cilindri multipli. Questa soluzione,
rispetto alle altre, presenta un maggiore ingombro (soprattutto se si
mantiene il cinematismo biella-manovella) e i volumi morti risultano
essere sensibilmente maggiori. Altri svantaggi quali i ritardi di scambio
termico, i problemi di trafilamento e l''elevato livello vibratorio dovuto alla
tipologia di cinematismo che muove il sistema, fanno sì che la
configurazione Gamma abbia un utilizzo di nicchia, principalmente
macchine di piccola taglia e modellini didattici.


Figura 1.16 Configurazione Gamma [4].
Capitolo 1
34
1.3 WhisperGen EU1

L''unità cogenerativa montata presso il Laboratorio di Micro
Cogenerazione (LMC) del Politecnico di Milano (gruppo GECOS '' Group
of Energy Conversion Systems), è un WhisperGen EU1 progettato
dall''azienda neozelandese Whisper Tech ltd e prodotto dalla spagnola
EHE. Si tratta di un micro-cogeneratore a ciclo Stirling alimentato a gas
naturale a pressione di rete, progettato per soddisfare i consumi domestici
di riscaldamento ed energia elettrica fino a un massimo di 1000 W elettrici,
con l''intento, quindi, di rimpiazzare la classica caldaia domestica.
Figura 1.17 Esempio di applicazione dello WhisperGen EU1 [7].
Il motore Stirling presente nell''unità WhisperGen è costituito da un
meccanismo Alpha a doppio effetto, illustrato nel paragrafo 1.2.1, a quattro
pistoni (configurazione FCDA Four Cycle Double Acting) il cui moto è
dato dall''espansione/contrazione del fluido di lavoro (azoto a 20 bar @
25°C) dovuta al suo riscaldamento e successivo raffreddamento. Il
meccanismo di conversione del moto da alternato a rotatorio, necessario
per l''accoppiamento con il generatore ad induzione 4 poli singola fase, è
denominato Wobble Yoke.

La macchina è dotata di due bruciatori catalitici premiscelati: il principale
contribuisce alla produzione di potenza elettrica e termica, fino ad una
potenza nominale di circa 8 kW termici, mentre il bruciatore secondario
fornisce esclusivamente potenza termica per un totale di circa 14 kW di
potenza elettrica. L''utilizzo con bassa potenza termica è fortemente
penalizzante dal punto di vista del rendimento e quindi causa di un Motore Stirling cogenerativo 35 maggior tempo di ritorno economico dell''investimento, pertanto sarebbe
opportuno, in questi casi, prevedere l''utilizzo di un sistema di accumulo.

Figura 1.18 Schema d''impianto WhisperGen EU1 [7].
L''unità opera come di seguito:
1. Il gas naturale viene bruciato nella camera di combustione, situata all''estremità superiore del motore. 2. Il calore generato espande l''azoto all''interno del motore, muovendo di conseguenza i pistoni. 3. Il moto dei pistoni aziona il generatore, che produce potenza elettrica. 4. L''acqua di ritorno viene preriscaldata dai gas caldi di combustione. Capitolo 1
36
5. L''acqua a questo punto passa attraverso a delle camicie nel motore, dove viene riscaldata ancora, e contemporaneamente raffredda la zona
''fredda' del motore situata all''estremità inferiore. 6. L''acqua calda viene espulsa verso l''utenza termica finale.

Di seguito sono elencate le principali caratteristiche tecniche del
cogeneratore WhisperGen EU1:

Combustibile: Gas naturale Potenza elettrica prodotta: 1 kW Potenza elettrica assorbita: 11-60 W Efficienza energetica: 96% Dimensioni: 49,1 x 83,8 x 56,3 cm Peso: 142 kg Rumore: 46 dB Potenza termica prodotta a 60-80°C: 7,5 '' 8,4 kW Potenza termica prodotta a 60-80°C con bruciatore
secondario:
13,2 '' 14,5 kW Tabella 1.2 Caratteristiche nominali del motore WhisperGen EU1. Motore Stirling cogenerativo 37

Principali componenti dell''unità
WhisperGen EU1:
1. Generatore elettrico 2. Motore Stirling 3. Bruciatore 4. Bruciatore ausiliario 5. Scambiatore di calore 6. Ventilatore centrifugo







Figura 1.19 Principali componenti dell''unità WhisperGen [7].


Figura 1.20 Vista esterna dell''unità WhisperGen (a) e cinematismo Wobble Yoke (b) [7]. Capitolo 1
38
39 Capitolo 2

Modello numerico per fluidi puri e
miscele monofase / bifase

Nel presente capitolo sarà presentato il modello numerico realizzato per
determinare le caratteristiche termodinamiche di acqua, gas puri e miscele
monofase e bifase. La necessità di realizzare tale modello è nata durante le
prime fasi di sperimentazione del motore Stirling, nelle quali il precedente
foglio di calcolo, basato sull''applicativo per Excel Refprop, era di difficile
condivisione su più computer e più piattaforme. Questo era dovuto al fatto
che per funzionare, l''applicativo, fa riferimento ad indirizzi locali del PC
su qui è installato. Ogni volta che il foglio di calcolo doveva essere
utilizzato da un altro PC era necessario modificare tutti gli indirizzi locali,
o in alternativa inserire nuovamente tutti i nuovi dati anche sul computer
sulla quale si voleva lavorare, rendendo il processo di revisione complicato
e laborioso.

Lo scopo principale di questo modello numerico, quindi, è di facilitare la
condivisione del foglio di calcolo tra diversi computer, e di poter
facilmente calcolare tutte le proprietà termodinamiche di fluidi puri e
miscele con una precisione adeguata. Le due alternative su cui
implementare il codice numerico sono state Matlab, e il linguaggio VBA,
attraverso l''Editor di Visual Basic per Excel. Il primo è un ambiente di
programmazione sofisticato e conosciuto nel settore dell''ingegneria, ma
che richiede necessariamente l''installazione dell''omonimo software per
poter essere utilizzato. Il secondo invece è un linguaggio di
programmazione meno usato, ma che si integra perfettamente in Excel, e
pertanto non richiede nient''altro tranne il famoso software di calcolo della
Microsoft, su qui peraltro è già implementato il foglio di calcolo per il
bilancio energetico del motore Stirling. Poiché lo scopo principale della
realizzazione di questo modello matematico è quello di rendere facilmente
condivisibile il foglio di calcolo, come già precedentemente detto, si è
deciso di utilizzare VBA come linguaggio di programmazione. Inoltre si è
deciso di implementare il codice in lingua inglese, ovvero tutti gli input di
ciascuna funzione dovranno essere in lingua inglese, così come lo sono
tutte le note inserite all''interno del programma che descrivono punto per
punto i vari step del codice.
Capitolo 2
40
Il modello numerico realizzato, è quindi costituito da tre diversi moduli:
''Water', per le proprietà termodinamiche dell''acqua, ''Gas', per le
proprietà termodinamiche dei gas puri (sono stati implementati i principali
gas d''interesse ingegneristico), e ''Mixture' per le proprietà
termodinamiche delle miscele monofase e bifase. Ogni modulo viene
chiamato dal foglio di calcolo tramite un''apposita funzione (LMC_water,
LMC_gas, LMC_mixture, dove ''LMC' è l''acronimo di ''Laboratorio di
Micro Cogenerazione' presso il quale si sono svolte le sperimentazioni).

La principale difficoltà riscontrata nell''implementare il modello numerico,
oltre al dover ''imparare' un nuovo linguaggio di programmazione, è stata
quella di trovare delle correlazioni numeriche facilmente implementabili e
allo stesso tempo adeguate in termini di precisione e stabilità.


2.1 Visual Basic for Application (VBA)

Visual Basic for Application (VBA) è un''implementazione di Visual Basic
inserita all''interno di applicazioni Microsoft, quali la suite Microsoft
Office, e altri programmi come AutoCAD o WordPerfect, che però
contengono solo un''implementazione parziale di VBA. Nonostante il suo
stretto legame con Visual Basic, VBA non può essere usato per eseguire
applicazioni stand-alone, ma è comunque possibile una certa
interoperatività fra applicazioni (ad esempio è possibile creare un report in
Word a partire da dati in Excel).

VBA è un linguaggio di programmazione ad alto livello, molto familiare
alla logica del nostro linguaggio naturale, i cui principali oggetti sono
subroutine e function. La subroutine, chiamata anche procedura macro,
esegue automaticamente un insieme di operazioni, nella cartella, foglio o
cella selezionate al momento del lancio. La function, invece, richiede come
input almeno un parametro numerico o testuale per almeno una variabile
indipendente.


2.2 Modulo ''Water'

Nel modulo Water, come precedentemente accennato, è contenuta la
stringa di codice per la determinazione delle proprietà termodinamiche
dell''acqua, sia in fase liquida che vapore, in funzione di temperatura e
pressione. L''implementazione di tale modello numerico è stata eseguita
seguendo la pubblicazione ''Release on the IAPWS Industrial Formulation Modello numerico per fluidi puri e miscele monofase / bifase 41 1997 for the Thermodynamic Properties of Water and Steam' della
IAPWS (The International Association for the Properties of Water and
Steam).

La IAPWS è un''associazione internazionale no-profit, di cui fa parte anche
l''Italia, con interesse nella determinazione delle proprietà termodinamiche
di acqua e vapore, di particolare importanza nell''ambito dei cicli di
produzione di potenza. I suoi principali obiettivi sono:
' Fornire formulazioni per la determinazione delle proprietà di vapore e acqua riconosciute a livello internazionale; ' Promuovere e coordinare la ricerca sul vapore, acqua e sistemi acquosi importanti per cicli termoelettrici; ' Raccogliere, valutare e diffondere i dati risultanti; ' Fornire un forum internazionale per lo scambio di esperienze, idee e risultati di ricerche sui mezzi acquosi ad alta temperatura.
La pubblicazione utilizzata per l''implementazione del modello numerico è
nata principalmente per un uso industriale, e consiste in un''accurata
approssimazione della "IAPWS Formulation 1995 for the Thermodynamic
Properties of Ordinary Water Substance for General and Scientific Use.",
comunemente detta IAPWS-95, che, come suggerito dal titolo, è adatta ad
un uso prettamente scientifico. Il mondo industriale però ha esigenze
diverse da quello scientifico; le proprietà dei fluidi possono essere
richiamate milioni di volte durante le fasi di progettazione, pertanto la
facilità d''implementazione e la velocità nell''eseguire i calcoli diventano
aspetti decisivi. Da qui in avanti, con il termine IAPWS, si farà riferimento
alla formulazione industriale, se non diversamente specificato.

La IAPWS consiste in un set di equazioni per differenti regioni che copre il
seguente range di pressione e temperatura:
273.15 ' T ' 1073.15 K P ' 100 MPa 1073.15 ' T ' 2273.15 K P ' 50 MPa Capitolo 2
42
Figura 2.1 Suddivisione in regioni della IAPWS [9].

In figura 2.2 è raffigurata la suddivisione nelle 5 regioni in cui l''intero
range di validità è stato suddiviso. L''equazione di base per descrivere le
regioni 1, 2 e 5, è un''equazione fondamentale per l''energia libera di Gibbs
in funzione di temperatura e pressione, mentre per la regione 3 viene
utilizzata un''equazione fondamentale per l''energia libera di Helmholtz,
funzione di densità e pressione, e la curva di saturazione nella forma ps(T). Le proprietà termodinamiche calcolate nativamente dalla IAPWS sono su
base massica, ma all''interno della VBA, tramite un opportuno comando, è
possibile convertirle in base molare.

I limiti di ciascuna regione possono essere direttamente individuati dalla
figura 2.2, ad eccezione del confine tra le regioni 2 e 3; questo confine è
individuato dalla cosiddetta equazione B-23, qui di seguito riportata:

(2.1)
dove ' = p/p* e θ = T/T* con p* = 1 MPa e T* = 1 K. I coefficienti da n1 a n3 sono riportati nella tabella 2.1.
In alternativa l''equazione (2.1) può anche essere espressa come:

[( ) ' ] ' (2.2) Modello numerico per fluidi puri e miscele monofase / bifase 43 L''equazione (2.1) descrive approssimativamente una line isoentropica; i
valori di entropia lungo questa curva di delimitazione variano tra s = 5.047
kJ Kg-1 K-1 e s = 5.261 kJ Kg-1 K-1.

i ni i ni 1 0.348 051 856 289 69 x 103 4 0.572 544 598 627 46 x 103 2 -0.116 718 598 799 75 x 101 5 0.139 544 598 627 46 x 102 3 0.101 929 700 393 26 x 10-2 Tabella 2.1 Valori numerici per i coefficienti delle equazioni (2.1) e (2.2).

2.2.1 Regione 1
Come già accennato in precedenza, l''equazione di partenza per la
determinazione delle proprietà termodinamiche della regione 1, è una
equazione fondamentale per l''energia specifica di Gibbs. Questa equazione
è espressa nella forma adimensionale:

( ) ( ) '' ( ) ( ) (2.2)
dove ' = p/p* e ' = T/T* con p* = 16.53 MPa e T* = 1386 K.

I coefficienti ni e gli esponenti Ii e Ji sono riportati nell''Allegato 1.
Tutte le proprietà termodinamiche possono essere derivate direttamente
dall''equazione (2.2) usando delle appropriate combinazioni dell''energia
libera di Gibbs e delle sue derivate. Tutte le relazioni tra le proprietà
termodinamiche e γ e le sue derivate sono riportate nella Tabella 2.2,
mentre le derivate di γ sono riportate nella Tabella 2.3.

Fin dalla ''5th International Conference on the Properties of Steam' a
Londra nel 1956, l''energia specifica e l''entropia specifica di liquido saturo
nel punto triplo sono state fissate pari a zero. Pertanto i coefficienti ni sono stati sviluppati in tal senso, restituendo un valore di entalpia specifica nel
punto triplo pari a:
Capitolo 2
44
Proprietà Relazione Volume specifico:
( ' ) ( ) Energia interna specifica:
( ' ) ( ' ) ( ) Entropia specifica:
( ' ) ( ) Entalpia specifica:
( ' ) ( ) Calore specifico isobarico:
( ' ) ( ) Calore specifico isocoro:
( ' ) ( ) ( ) Velocità del suono:
[ ( ' ) ] ' ( ) ( ) [
] , [
] , [ ] , [ ] , [ ]
Tabella 2.2 Relazioni proprietà termodinamiche '' Regione 1. Modello numerico per fluidi puri e miscele monofase / bifase 45 '' ( ) ( ) '' ( ) ( ) '' ( )( ) ( ) '' ( ) ( ) '' ( ) ( )( ) '' ( ) ( ) Tabella 2.3 Derivate dell''energia libera di Gibbs.
Dalla Tabella 2.2, si può vedere come tutte le proprietà termodinamiche
d''interesse possono essere calcolate a partire dall''espressione dell''energia
libera di Gibbs e delle sue derivate.

I limiti di validità della regione 1 risultano:
273.15 K ' T ' 623.15 K ps(T) ' p ' 100 MPa

2.2.2 Regione 2
Anche per questa regione l''equazione di partenza è l''equazione
fondamentale per l''energia libera di Gibbs g(p,T), espressa nella sua forma
adimensionale. In questo caso, però, la formulazione adimensionale è la
risultante di due contributi separati: una prima parte relativa ai gas ideali
γ°, e una parte residua γr.

( ) ( ) ( ) ( ) (2.3)
dove ' = p/p* e ' = T/T* con p* = 1 MPa e T* = 540 K. Capitolo 2
46
L''equazione per la parte relativa ai gas ideali dell''energia libera di Gibbs
adimensionalizzata risulta:

( ) '' (2.4)
La parte residua γr, invece, risulta essere:

( ) '' ( ) (2.5)
dove ' = p/p* e ' = T/T* con p* = 1 MPa e T* = 540 K. I coefficienti
presenti nell''equazione (2.4) e (2.5) sono riportati nell''Allegato 2.
Proprietà Relazione Volume specifico:
( ' ) ( ) ( ) Energia interna specifica:
( ' ) ( ' ) ( ) ( ) ( ) Entropia specifica:
( ' ) ( ) ( ) ( ) Entalpia specifica:
( ' ) ( ) ( ) Calore specifico isobarico:
( ' ) ( ) ( ) Calore specifico isocoro:
( ' ) ( ) ( ) ( ) Velocità del suono:
[ ( ' ) ] ' ( ) ( ) ( ) ( ) Tabella 2.4 Relazioni proprietà termodinamiche '' Regione 2. Modello numerico per fluidi puri e miscele monofase / bifase 47 Come per la regione 1, di seguito sono elencate tutte le componenti delle
derivate dell''energia libera di Gibbs nella sua forma adimensionale,
necessarie per il calcolo delle proprietà termodinamiche dell''acqua
all''interno della regione 2. Nella Tabella 2.5 sono presenti le derivate della
componente ideale γo, mentre nella Tabella 2.6 sono presenti le derivate
della componente relativa γr.
'' ' ' '' '' ( ) [ ] , [ ] , [ ] , [ ] , [
]
Tabella 2.5 Derivate della parte di gas ideale dell''energia libera di Gibbs. Capitolo 2
48
'' ( ) '' ( ) '' ( ) ( ) '' ( ) '' ( )( ) '' ( ) [ ] , [ ] , [ ] , [ ] , [ ] Tabella 2.6 Derivate della parte residua dell''energia libera di Gibbs.

L''equazione (2.3), e quindi tutte le proprietà ad essa correlate, vale per il
seguente range di pressione e temperatura:
273.15 K ' T ' 623.15 K 0 ' p ' ps(T)Eq.(2.6) 623.15K < T ' 863.15 K 0 ' p ' ps(T)Eq.(2.1) 863.15K < T ' 1073.15 K 0 ' p ' 100 MPa
I primi due range di pressione sono a loro volta delimitati superiormente da
un valore di pressione che dipende dalla temperatura. L''equazione (2.1) cui
fa riferimento il secondo intervallo, già introdotta nel capitolo 2.2,
rappresenta il confine tra la regione 2 e la regione 3. Modello numerico per fluidi puri e miscele monofase / bifase 49 L''equazione (2.6), invece, rappresenta l''equazione della pressione di
saturazione in funzione della temperatura:


[ ( ) ' ] (2.6)
Dove p*= 1MPa, e:


Dove ' rappresenta:

( ' ) (2.6)
con p*= 1MPa e T*=1 K. I valori dei coefficienti ni sono riportati nell''Allegato 2.


2.2.3 Regione 3
Per la regione 3 l''equazione di partenza è un''equazione fondamentale per
l''energia libera di Helmholtz, anch''essa nella sua forma adimensionale:

( ) ( ) '' (2.7)
dove δ = ρ / ρ*, ' = T* / T con ρ* = 322 kg m-3, T* = 647.096 K.

Come si vede l''equazione fondamentale, e quindi le sue derivate attraverso
le quali calcolare le proprietà termodinamiche, sono funzione della densità,
a sua volta variabile all''interno della Regione 3.

Per la determinazione della densità in funzione della temperatura e
pressione (gli input del modello termodinamico), la IAPWS fornisce una
release supplementare denominata ''Supplementary Release on Backward
Equations for Specific Volume as a Function of Pressure and Temperature Capitolo 2
50
v(p,T) for Region 3 of the IAPWS Industrial formulation 1997 for the
Thermodynamic Properties of Water and Steam'.

La Regione 3 viene a sua volta divisa in 20 sotto-regioni, caratterizzate da
coefficienti e modelli numerici differenti. Figura 2.2 Suddivisione in sotto-regione per la Regione 3 [10]. Modello numerico per fluidi puri e miscele monofase / bifase 51 Figura 2.3 Ingrandimento delle sotto-regioni da 3c a 3r [10].
Capitolo 2
52
' evidente la notevole difficoltà nell''individuare correttamente i confini
delle 20 sotto-regioni, nonché la consistente quantità di dati necessari per
calcolare tutte le proprietà. Considerando l''elevato range di temperatura
cui si rivolge questa regione, si è deciso di non implementala, almeno per
ora, all''interno del modello matematico.


2.2.4 Regione 4
La Regione 4 rappresenta la linea di saturazione dell''acqua, descritta da
un''equazione quadratica che può essere direttamente risolta per quanto
riguarda pressione di saturazione ps e temperatura di saturazione Ts.
(2.8)
Dove,

( ' ) ' (2.8a)
( ' ) (2.8b)
Con p*= 1MPa e T*=1 K. I valori dei coefficienti ni sono riportati nell''Allegato 3.

La soluzione dell''equazione (2.8), in merito alla pressione di saturazione,
risulta essere:


[ ( ) ' ] (2.9)
Nella quale p*= 1MPa, e:


La soluzione in termini di temperatura di saturazione è:

[( ) ( )] ' (2.10) dove T*=1 K, e Modello numerico per fluidi puri e miscele monofase / bifase 53
( ) '
con β in accordo all''equazione (2.8a).

Il range di validità delle equazioni (2.8), (2.9) e (2.10) è:
273.15 K ' T ' 647.096 K 611.213 Pa ' p ' 22.064 MPa

2.2.5 Regione 5
Anche per la regione 5, ad alta temperatura, l''equazione di partenza è
l''equazione fondamentale per l''energia libera di Gibbs g(p,T), espressa
nella sua forma adimensionale. Come nel caso della regione 2, la
formulazione adimensionale è la risultante di due contributi separati: una
prima parte relativa ai gas ideali γ°, e una parte residua γr.

( ) ( ) ( ) ( ) (2.11)
dove ' = p/p* e ' = T/T* con p* = 1 MPa e T* = 1000 K.
L''equazione per la parte relativa ai gas ideali dell''energia libera di Gibbs
adimensionalizzata risulta:

( ) '' (2.12)

La parte residua γr, invece, risulta essere:

( ) '' (2.13) Capitolo 2
54
dove ' = p/p* e ' = T/T* con p* = 1 MPa e T* = 1000 K. I coefficienti
presenti nell''equazione (2.4) e (2.5) sono riportati nell''Allegato 4.

Come per le altre regioni del modulo ''Water', nelle Tabelle 2.7, 2.8 e 2.9
sono riportate le derivate della componente ideale γo e residua γr
dell''equazione (2.11) e le correlazioni per le proprietà termodinamiche.

'' ' ' '' '' ( ) [ ] , [ ] , [ ] , [ ] , [
]
Tabella 2.7 Derivate della parte di gas ideale dell''energia libera di Gibbs.





Modello numerico per fluidi puri e miscele monofase / bifase 55 '' '' '' ( ) '' '' ( ) '' [ ] , [ ] , [ ] , [ ] , [ ] Tabella 2.8 Derivate della parte residua dell''energia libera di Gibbs.

Proprietà Relazione Volume specifico:
( ' ) ( ) ( ) Energia interna specifica:
( ' ) ( ' ) ( ) ( ) ( ) Entropia specifica:
( ' ) ( ) ( ) ( ) Capitolo 2
56
Entalpia specifica:
( ' ) ( ) ( ) Calore specifico isobarico:
( ' ) ( ) ( ) Calore specifico isocoro:
( ' ) ( ) ( ) ( ) Velocità del suono:
[ ( ' ) ] ' ( ) ( ) ( ) ( ) Tabella 2.9 Relazioni proprietà termodinamiche '' Regione 5.
L''equazione (2.11) ha validità nel campo di pressioni e temperature:
1073.15 K ' T ' 2273.15 K 0 ' p ' 50 MPa

2.2.6 Struttura modulo ''Water'
In generale, i tre moduli della VBA, sono strutturati come segue: la prima
parte, ''Declaretions', nella quale sono dichiarate le principali variabili e
costanti, necessarie per il calcolo numerico. Successivamente è presente
una sezione denominata ''Gateways Function', nel quale viene strutturata
la funzione pubblica con cui l''utente richiama il modulo, e che ha il
compito di richiamare le diverse funzioni necessarie per il calcolo della
proprietà richiesta dall''utente (subroutine nel caso si voglia eseguire una
serie di operazioni, o private function nel caso si voglia restituire un
valore). In pratica, per quanto riguarda il modulo Water, a seconda della
proprietà termodinamica richiesta, e della temperatura e pressione
impostata, la sezione Gateways Functions fa in modo di implementare i
calcoli necessari alla sola proprietà termodinamica richiesta, evitando
quindi di fare calcoli inutili e dispendiosi. Per ultima è presente la sezione
''Property Functions', nella quale sono inserite tutte le funzioni richiamate
dalla Gateways Functions.

Il modulo ''Water' è richiamato dalla function ''LMC_water', nella quale
sono implementate le 5 regioni come descritto nei paragrafi precedenti. La Modello numerico per fluidi puri e miscele monofase / bifase 57 funzione necessita come input la proprietà termodinamica desiderata e la
base cui fare riferimento (massica ''mass', che rappresenta l''output nativo
della IAPWS, o molare ''mole'). Come input opzionali (ovvero tutti gli
ingressi necessari solo per il calcolo di alcune proprietà termodinamiche),
la funzione riceve la temperatura in Kelvin ''K', la pressione assoluta in
''MPa' ed una ulteriore variabile, denominata ''Reference', con lo scopo di
indicare se si vuole calcolare l''entalpia (o l''entropia) inglobando anche
l''entalpia (o l''entropia) di formazione oppure no. Nella tabella sottostante
sono presenti le unità di misura di tutte le funzioni del modulo Water:

Base massica Base molare Volume specifico (v) m3 kg-1 m3 kmol-1 Calore specifico (Cp) kJ kg-1 K-1 kJ kmol-1 K-1 Calore specifico (Cv) kJ kg-1 K-1 kJ kmol-1 K-1 Entalpia (h) kJ kg-1 kJ kmol-1 Energia interna (u) kJ kg-1 kJ kmol-1 Entropia (s) kJ kg-1 K-1 kJ kmol-1 K-1 Velocità del suono (w) m s-1 m s-1 Temp. di saturazione (Tsat) K K Press. di saturazione (Psat) MPa MPa Tabella 2.10 Unità di misura VBA.

2.2.7 Energia libera di Gibbs
L'energia libera di Gibbs (o entalpia libera) è una funzione di stato usata in
termodinamica e termochimica per rappresentare l'energia libera nelle
trasformazioni isotermobariche (a cui appartengono la maggior parte delle
reazioni chimiche), e può essere usata per determinare la spontaneità di una
reazione.

L''energia libera di Gibbs ''G' è definita come l''opposto della trasformata
di Legendre dell''entalpia ''H' rispetto all''entropia ''S':

( ) ( ) (2.14)
dove ''T' è la temperatura assoluta. Quest''ultima equivalenza segue dalla
definizione di entropia. In effetti essendo l''entalpia a sua volta l''opposto Capitolo 2
58
della trasformata dell''energia interna ''U' in base al primo principio della
termodinamica, ''G' rappresenta la trasformata inversa dell''energia interna
rispetto all''entropia e al volume.

Il differenziale totale della funzione energia libera di Gibbs (2.14) risulta:

(2.15)
Considerando il differenziale totale della funzione entalpia:

(2.16)
si ottiene quindi:

(2.17)
Dal primo principio della termodinamica si ha che:

(2.18)
quindi sostituendo la (A.5) nella (A.4) si ottiene:

(2.19)
Se il processo è spontaneo, ovvero irreversibile, allora per tale processo è
valida la disuguaglianza di Clausius:
(2.20)
da cui si ricava:

(2.21)
L''equazione (A.6), quindi, può essere riscritta nella disuguaglianza:

(2.22)
ovvero:

(2.23)
Modello numerico per fluidi puri e miscele monofase / bifase 59 Quindi se il processo avviene a temperatura e pressione costante, poiché
risulterebbe dT=0 e dp=0 risulta:

(2.24)
Questo permette di affermare che in un processo spontaneo a temperatura e
pressione costante, l''energia libera ''G' diminuisce. Di conseguenza se
dG > 0, il processo non è spontaneo nel verso indicato, mentre lo sarebbe
nel verso opposto. Infine se dG = 0 il sistema si trova all''equilibrio e non
ha alcuna tendenza ad evolvere né in un senso né nel verso opposto.

Dall''espressione (2.14) si possono valutare come i due termini concorrano
alla spontaneità della reazione:
' Se la reazione è esotermica ('H < 0) e avviene con un aumento di entropia ('S > 0), allora essa sarà sempre spontanea, poiché 'G < 0
sempre. ' Se la reazione è endotermica ('H > 0) e avviene con una diminuzione di entropia ('S < 0), entrambi i termini saranno sfavorevoli alla
spontaneità del processo, poiché 'G sarà sempre positivo.
Quando 'H e 'S hanno lo stesso segno, la temperatura ''T' diventa
decisiva per determinare la spontaneità o meno del processo: se T'S > 'H
allora 'G < 0 e il processo sarà spontaneo, mentre se 'H e 'S sono
entrambi negativi sarà necessario che la temperatura sia sufficientemente
bassa affinché la loro differenza abbia valore negativo e la reazioni risulti
spontanea. Esisterà inoltre una temperatura tale per cui 'H e T'S saranno
numericamente identici e quindi 'G sarà pari a zero, ovvero il sistema sarà
all''equilibrio.


2.2.8 Energia libera di Helmholtz
L''energia libera di Helmholtz è una funzione di stato utilizzata in
termodinamica per rappresentare l''energia libera in una trasformazione
isotermocora, ovvero a temperatura e volume costante.

L''energia libera di Helmholtz ''F' è definita come l''opposto della
trasformata di Legendre dell''energia interna ''U' rispetto all''entropia ''S'.

( ) ( ) (2.25) Capitolo 2
60
dove ''T' è la temperatura assoluta. Quest''ultima equivalenza segue dalla
definizione di entropia.

La disuguaglianza di Clausius impone che:

(2.26)
Nel caso in cui la sola forma di lavoro sia quella di pressione-volume
(espansione o compressione), si ha che il calore diventa il differenziale
esatto dQ = dU, e quindi la disuguaglianza può anche essere scritta nel
seguente modo:

(2.27)
Introducendo l''equazione (B.1) nella (B.3), quest''ultima può essere
riscritta come:

(2.28)
Questa relazione indica che nelle trasformazioni isotermocore l''energia
libera di Helmholtz diminuisce per un processo spontaneo (differenziale
negativo), mentre è nullo ad un valore minimo (differenziale nullo) per un
processo reversibile, cioè in condizioni di equilibrio.


2.3 Modulo ''Gas'

Nel modulo ''Gas' sono state implementate le proprietà termodinamiche di
tutte quelle specie chimiche d''interesse ingegneristico per quanto riguarda
un bilancio energetico di un ciclo termodinamico generico, ovvero tutte
quelle specie solitamente presenti in una combustione. In questo caso sono
stati inseriti i principali composti del gas naturale (Metano, Etano,
Propano, Normal-Butano, Iso-Butano, Normal-Pentano, Iso-Pentano ed
Esano), e i principali componenti presenti nell''aria ambiente (azoto,
ossigeno, argon, biossido di carbonio, monossido di carbonio), ed in
aggiunta l''idrogeno. Ovviamente non è stata inclusa l''acqua poiché per
quest''ultima vi è il modulo dedicato.

In questo modulo è possibile calcolare entalpia, entropia, calore specifico a
pressione costante, massa molare, potere calorifico inferiore e superiore
delle specie precedentemente elencate, supposte tutte nello stato aeriforme. Modello numerico per fluidi puri e miscele monofase / bifase 61 La stima delle proprietà termodinamiche dei singoli composti, in funzione
questa volta della sola temperatura, viene eseguita implementando
l''equazione di Shomate [11]. Si tratta di un metodo ormai validato per la
determinazione delle proprietà termodinamica attraverso l''uso di
un''equazione polinomiale. In letteratura esistono diverse costanti da
utilizzare in funzione del tipo di polinomio utilizzato per la regressione, di
queste l''equazione di Shomate è una delle più accurate e usate. I
coefficienti implementati all''interno della VBA sono stati presi tramite il
WebBook della NIST (''National Institute of Standards and Technology'),
il quale riporta i coefficienti dell''equazione di Shomate per ogni composto,
divisi per classi di temperatura.

Il National Institute of Standards and Technology (NIST) è un'agenzia del
governo degli Stati Uniti d'America che si occupa della gestione delle
tecnologie. Fa parte del Dipartimento del Commercio e il suo compito è la
promozione dell'economia Americana attraverso il lavoro con l'industria
per sviluppare standard, tecnologie e metodologie che favoriscano la
produzione e il commercio. Inizialmente si chiamava National Bureau of
Standards (NBS).

Di seguito sono riportati i polinomi interpolatori, come indicati nel
WebBook della NIST, per la stima di calore specifico a pressione costante,
entalpia ed entropia dei composti puri inseriti nel modulo ''Gas', tutti
riferiti alla pressione costante di 1 bar assoluto:

' (2.29)
' ' ' ' (2.30)
' ' ( ) ' (2.31)
Come è logico aspettarsi, i polinomi riferiti al delta entalpico e all''entropia
delle singole specie, sono ricavabili a partire dall''equazione del calore
specifico, a meno di termini costanti, rispettivamente derivando ed
integrando l''equazione (2.14).

Dalla pagina WebBook, oltre alle costanti da utilizzare per l''equazione di
Shomate su base molare, sono disponibili anche i dati relativi a massa
molare, entalpia di formazione ed entropia di formazione; tali informazioni
sono stati utilizzati per calcolare, in analogia al modulo Water, le proprietà
specifiche alla massa e per il calcolo di entalpia ed entropia inglobando le
grandezze di riferimento. In questo modo l''equazione (2.15) viene in realtà
implementata nella forma seguente: Capitolo 2
62
' ' ' ' (2.32)
Poiché in realtà la costante H dell''equazione (2.15) è proprio il valore
dell''entalpia di formazione della specie considerata, all''interno del modulo
''Gas' non si è fatto altro che porre uguale a zero tale costante per ogni
specie chimica.

Per quanto riguarda l''equazione (2.16), è stato deciso di affinarla
introducendo anche la pressione come variabile indipendente, pertanto il
polinomio presente nel modulo ''Gas' risulta essere:

' ' ( ) ' ( ) (2.33)
dove R = 8.314 kJ/kmol K, ovvero la costante universale dei gas, e Prif risulta essere la pressione di riferimento, che in questo caso è pari ad 1 bar.

Il NIST non fornisce per tutti i composti presenti nella VBA le costanti da
implementare nella Shomate; per questi composti (Etano, Propano,
Normal-Butano, Iso-Butano, Normal-Pentano, Iso-Pentano ed Esano) sono
però tabulati i valori di Cp(T). Partendo da questi valori, attraverso un''interpolazione polinomiale del terzo ordine, è stato possibile ricavare i
valori delle costanti da applicare all''equazione seguente:

' ' ' ' (2.34)
analogamente viene calcolata l''entropia in funzione della temperatura.

Il potere calorifico superiore e inferiore sono calcolati ipotizzando una
combustione stechiometrica:

( ) (2.35)
dove sono rispettivamente l''entalpia di formazione dell''idrocarburo considerato, l''entalpia di formazione dell''acqua (allo stato
aeriforme per il potere calorifico inferiore e allo stato liquido per il potere
calorifico superiore) e l''entalpia di formazione per l''anidride carbonica.


2.3.1 Struttura modulo ''Gas'
Il modulo ''Gas' è strutturato nel medesimo modo del modulo ''Water': è
presente una prima sezione ''Declaretions' per la dichiarazione delle
variabili e delle costanti, una ''Gateways Function' per ottimizzare il Modello numerico per fluidi puri e miscele monofase / bifase 63 processo di calcolo ed infine una ''Property Functions' in cui sono
presenti le funzioni descritte sopra per il calcolo delle proprietà richieste.

Analogamente al modulo ''Water', il modulo ''Gas' è richiamato attraverso
la function ''LMC_gas', la cui sintassi è analoga a quella del modulo
''Water', tranne per il fatto che l''input iniziale della funzione riguarda la
specie chimica di cui si vuole calcolare la proprietà termodinamica. Come
già accennato nei capitoli precedenti, si è scelto di implementare il codice
numerico in lingua inglese, pertanto anche l''input delle specie dovrà essere
in inglese.


2.4 Modulo ''Mixture'

Il modulo ''Mixture', come già anticipato consente di calcolare le proprietà
termodinamiche di miscele monofase, allo stato gassoso, e bifase liquido-
gas. In quest''ultimo caso, si è semplificato il modello ipotizzando che la
fase liquida sia composta solo da acqua, tutti gli altri composti sono stati
trattati come incondensabili. In entrambi i casi, per il calcolo delle
proprietà termodinamiche, si è considerata una miscela ideale di gas ideali,
ovvero una miscela nella quale tutte le proprietà parziali molari di un
componente sono uguali alle corrispondenti proprietà molari del composto
puro come gas perfetto, calcolate alla stessa temperatura della miscela e ad
una pressione pari alla pressione parziale del componente. In questo modo
si possono trattare le particelle di ogni singolo gas che compone la miscela
come masse puntiformi, trascurando le forze intermolecolari che
dipendono dalla loro natura chimica. Viceversa, se gli effetti di
miscelamento non possono essere trascurati, la miscela sarà detta ''reale'.

In questo modo è possibile applicare quella che è comunemente chiamata
legge di Gibbs-Dalton, che può essere riassunta come segue:

''La pressione totale di una miscela ideale di gas è uguale alla somma
delle pressioni parziali di ogni singolo componente; si definisce pressione
parziale di un componente la pressione che assumerebbe tale componente
qualora occupasse da solo, alla medesima temperatura e nella quantità in
cui è presente nella miscela, lo stesso volume occupato dalla miscela.
L''energia interna totale di una miscela ideale di gas è uguale alla somma
dell''energia interna di ogni singolo componente; si definisce energia
interna di un componente l''energia interna propria alla quantità del
singolo componente alla medesima temperatura della miscela.'
Capitolo 2
64
Si ha pertanto:

'' '' (2.36)
Tenendo conto della definizione di entalpia, è facile dedurre come anche
l''entalpia H della miscela sia data dalla somma delle entalpie Hi dei singoli componenti.

'' (2.37)
Analogamente si possono calcolare il calore specifico e l''entropia della
miscela considerata:

'' '' (2.38)
Il codice numerico implementato nella VBA, quindi, non fa altro che
calcolare le proprietà di ogni singolo composto utilizzando le funzioni
dedicate dei moduli Water e Gas, e quindi restituire la proprietà
termodinamica della miscela.

Fino ad ora si è trattato la miscela come costituita dalla sola fase gassosa,
ma nel caso di presenza di acqua tra i composti, è necessario effettuare una
verifica sulla presenza o meno di condense. Questo viene fatto attraverso
l''esecuzione di una VLE (Vapor-Liquid Equilibrium); il procedimento
adottato viene descritto nel capitolo seguente.


2.4.1 VLE: Vapor-Liquid Equilibrium
L'equilibrio liquido-vapore (in inglese Vapor-Liquid Equilibrium o VLE) è
la condizione in cui due fasi, una fase liquida e una fase vapore, sono in
equilibrio termodinamico tra loro. L'equilibrio che si viene a creare è di
tipo dinamico, ovvero la velocità di evaporazione del liquido eguaglia la
velocità di condensazione del vapore. Infatti l'interfaccia liquido-vapore è
interessata da continui scambi di materia tra le due fasi.

Come precedentemente accennato, il modello numerico è stato
semplificando considerando solo l''acqua come gas condensabile, pertanto
la VLE viene eseguita solamente in presenza di quest''ultima tra i composti
della miscela. Nel caso non sia presente acqua, la VLE non viene eseguita,
e la miscela viene trattata come monofase allo stato gassoso, ed il calcolo Modello numerico per fluidi puri e miscele monofase / bifase 65 delle proprietà termodinamiche è eseguito secondo quanto descritto nel
capitolo 2.4. Viceversa in caso di presenza di acqua, si esegue per prima
cosa un controllo sulla sua quantità presente nella miscela; se la pressione
parziale dell''acqua in miscela è inferiore alla sua pressione di saturazione,
alla temperatura della miscela, allora non è presente una frazione liquida e
ancora una volta si utilizzano le equazione descritte precedentemente per il
calcolo delle proprietà termodinamiche della miscela. Se invece la
pressione parziale risulta maggiore della pressione di saturazione, allora un
parte di acqua si trova allo stato liquido; in questo caso è necessario
eseguire la VLE per determinare la composizione del vapore, la frazione di
liquido e la frazione di vapore (la composizione, per quanto detto in
precedenza è nota e si tratta di sola acqua).

Di seguito viene riportato il procedimento adottato per il calcolo
dell''equilibrio liquido-vapore (il pedice ''L' indica la frazione liquida,
mentre il pedice ''V' indica la frazione di vapore, ''F' indica la somma
della frazione liquida e di vapore):


Equilibrio Meccanico (2.39)
Equilibrio Termico (2.40)
Si indica con xi e yi rispettivamente le frazioni di liquido e vapore del componente i-esimo della miscela, mentre con zi indica la frazione totale presente nella miscela.
Risulterà quindi:

Bilancio miscela (2.41)

(2.42)
Pertanto valutando la frazione ''α' di liquido è possibile calcolare la
frazione ''β' di vapore; per il calcolo della frazione liquida si suppone che
la fase vapore sia alla saturazione, pertanto tutta quantità d''acqua
eccedente condenserà:

(2.43)
dove rappresenta la massima quantità di acqua presente alla saturazione nella fase vapore della miscela, ed è ricavabile a partire dalla Capitolo 2
66
pressione di saturazione dell''acqua alla medesima temperatura della
miscela, e dalla pressione della miscela stessa.

(2.44)
Combinando l''equazione (2.26) con la (2.27), si ottiene:

( ) (2.45)
Come detto più volte, l''unico composto in fase liquido può essere solo
l''acqua, per il modello considerato, pertanto risulta:

{ (2.46)
A questo punto è possibile calcolare facilmente le frazioni molari yi che compongono la fase vapore della miscela:

(2.47)

Una volta note le frazioni molari della fase liquida e della fase vapore, le
proprietà termodinamiche della miscela bifase possono essere calcolate
facilmente nel seguente modo nel seguente modo:

'' ( ) '' ( ) (2.48)
dove ''T' è la temperatura della miscela, '' ' la sua pressione e '' ' la pressione parziale della specie i-esima.

L''equazione (2.33) può essere usate anche per le altre grandezze
termodinamiche calcolabili.


2.4.2 Struttura modulo ''Mixture'
La struttura del modulo ''Mixture' riprende quella dei moduli ''Water' e
''Gas' precedentemente descritti. Anche in questo caso il modulo è diviso
in ''Declaretions', ''Gateways Function' e ''Property Functions', che
svolgono la medesima funzione. Modello numerico per fluidi puri e miscele monofase / bifase 67 Il modulo è chiamato dalla function ''LMC_mixture', che ha la medesima
sintassi della funzione per richiamare il modulo gas tranne per il fatto che il
secondo input riguarda le frazioni molari che compongono la miscela.
Capitolo 2
68
69 Capitolo 3

Prove sperimentali al variare della
pressione del fluido di lavoro

Nel presente capitolo verranno illustrate le attività sperimentali svolte
durante il periodo di tesi. In precedenza sono state svolte sperimentazione
sullo WhisperGen EU1 per valutarne le prestazioni in termini di
rendimento elettrico e rendimento termico, riferiti al potere calorifico
inferiore e superiore (PCI e PCS). I test sono stati eseguiti variando la
portata di acqua cogenerativa, entro i limiti imposti dal costruttore, e la
temperatura dell''acqua di ritorno (30°C, 50°C e 70°C) per verificare
l''andamento delle prestazioni dell''unità cogenerativa al variare dei
principali parametri operativi. Il rendimento termico ha mostrato una forte
dipendenza dalla temperatura dell''acqua di ritorno; i migliori risultati sono
stati rilevati nelle prove con temperatura dell''acqua di ritorno da 50 a
30 °C, dove la temperatura dei fumi risulta pari a circa 30 °C, nella quale si
verifica una notevole condensazione dei fumi e quindi una notevole
potenza termica scambia trasferita all''acqua, raggiungendo rendimenti di
primo principio quasi unitari rispetto al PCS. Per quanto riguarda il
rendimento elettrico si è rilevato un leggero incremento all''aumentare della
temperatura dell''acqua di ritorno, con un massimo a 70°C.
Figura 3.1 Rendimento elettrico e termico riferito al PCS. 0 20 40 60 80 100 0.140 0.194 0.247 0.140 0.194 0.247 0.194 0.247 R e n d im e n to η [% ] Portata acqua [kg/s] ηel ηth 30°C 50°C 70°C Capitolo 3
70
In tutte queste prove, però, non si era intervenuti sulla pressione del fluido
di lavoro (azoto), che è rimasta al suo valore nominale pari a 20 bar alla
temperatura di 25°C; si è quindi deciso di compiere delle prove al variare
della pressione dell''azoto. Scendere al di sotto del valore di progetto non
rappresenta alcun rischio per l''integrità della macchina, mentre non si
conosce il limite di pressione oltre la quale è possibile danneggiare lo
Stirling; per questo motivo si è deciso di effettuare le prove solo al di sotto
del valore nominale, per evitare di rendere inutilizzabile la macchina e
quindi non poter effettuare tutte le prove programmate per il futuro.

Gli scopi di questi test sono molteplici: da un lato si vogliono misurare le
prestazioni della macchina al variare della pressione del fluido di lavoro, in
termini di potenza elettrica e termica, mentre dall''altro lato si vuole
verificare il corretto funzionamento del modello numerico descritto nel
Capitolo 2. L''obiettivo principale, comunque, è quello di confrontare le
prestazioni rilevate sperimentalmente con quelle valutate da un modello
numerico dello WhisperGen EU1 implementato dall''Ing. Nicola Fergnani.
Il confronto, tra il modello numerico e i rilievi sperimentali, sarà di estrema
utilità per una eventuale affinazione del modello numerico stesso.

Di seguito saranno descritti brevemente il modello numerico del motore
Stirling, la tipologia di test effettuati ed i risultati delle prove sperimentali,
sia in termini di confronto con il modello numerico, che in termini di
potenza elettrica e termica sviluppati.


3.1 Brevi cenni sul modello numerico del motore Stirling

Il motore Stirling è caratterizzato da un funzionamento piuttosto complesso
e dalla presenza di una serie di trasformazioni termodinamiche che, come
descritto nel Capitolo 1, si discostano fortemente dal ciclo ideale.

Tra i diversi metodi numerici riportati in letteratura per la modellizzazione
e la simulazione di motori Stirling di diversa tipologia, il primo tentativo di
analisi è riconosciuto a Finkelstein. Seguirono numerosi miglioramenti di
questo primo modello, tra i quali il più importante è dovuto ad Urieli nel
1973, il cui approccio risulta tuttora un riferimento per la modellizzazione
del motore Stirling nelle sue varianti, e sul quale si basa il modello di
simulazione dell''Ing. Fergnani. Il modello numerico proposto da Urieli,
specialmente nella sua variante simple analysis, rappresenta infatti un buon
compromesso tra semplicità di calcolo e accuratezza della soluzione.
Prove sperimentali al variare della pressione del fluido di lavoro 71 Nel seguito viene descritto il modello di simulazione basato sull''analisi
adiabatica di Urieli, che costituisce la base per il successivo sviluppo
dell''analisi semplificata di secondo ordine, trattata nel par. 3.1.2 e
finalizzata alla valutazione dei fattori di non idealità che intervengono nel
ciclo.

3.1.1 Modello adiabatico di Urieli
Urieli condusse a partire dal 1977 numerosi studi relativi al ciclo Stirling,
partendo da modelli semplificati per giungere, negli anni successivi e con
la collaborazione di Berchowitz, allo sviluppo di analisi più complesse.
Partendo dal ciclo ideale (analisi isoterma - first order), essi sostituirono le
trasformazioni isoterme proprie del ciclo ideale con trasformazioni
adiabatiche, più rappresentative della reale evoluzione del gas all''interno
del motore (analisi adiabatica '' second order), introducendo in seguito un
ulteriore livello di analisi (''simple analysis') con il proposito di valutare
gli effetti delle non idealità presenti nel ciclo e delle relative conseguenze
in termini di prestazioni dello stesso.
Tutti i livelli di analisi si basano su un approccio a parametri concentrati,
che prevede la schematizzazione del motore mediante cinque elementi
connessi in serie.

Figura 3.2 Schema semplificato del motore Stirling [12].

In figura risultano visibili il volume di compressione (c), il refrigeratore o
cooler (k), il rigeneratore (r), il riscaldatore o heater (h) e il volume di
espansione (e). Tale approccio implica che ogni componente sia
considerato un''entità omogenea, in cui il fluido ivi contenuto è
caratterizzato istantaneamente dalla massa complessiva interna alla cella
m, dalla temperatura assoluta T e dal volume della cella V, che risulta
variabile nel caso delle celle di compressione ed espansione. La pressione,
invece, viene considerata uniforme in tutto il motore, trascurando in prima
analisi gli effetti delle perdite di carico. Capitolo 3
72
Il modello adiabatico di Urieli si basa sull''assunzione che nelle celle di
compressione ed espansione abbiano luogo trasformazioni perfettamente
adiabatiche, e che gli scambi di calore verso l''esterno intervengano per
effetto dalla transizione del gas da una cella all''altra. Si trascurano quindi
le irreversibilità di scambio termico, imponendo quindi che la temperatura
del gas all''interno delle celle di cooler ed heater sia la medesima della
temperatura di parete delle suddette celle. Un altro aspetto importante
riguarda l''assegnazione di temperature condizionali di ingresso del fluido
nelle diverse celle, come schematizzato in Figura 3.3.
Figura 3.3 Schema di impostazione del modello adiabatico e profilo di temperatura all''interno del motore [12].

Le frecce rappresentano la direzione positiva del fluido di lavoro, definita
arbitrariamente dal volume di compressione a quello di espansione. Le
temperature Tck, Tkr, Trh, e The sono definite condizionali poiché
dipendono dalla direzione del flusso. Per la direzione rappresentata in
figura, esse indicano che il gas che fluisce dal volume di compressione a
quello del cooler verrà a trovarsi ''istantaneamente' alla temperatura di
parete del cooler, che coincide anche con la temperatura dell''estremo
freddo del rigeneratore. Il gas che dal rigeneratore fluisce nell''heater non
subirà variazioni di temperatura in quanto, al presente livello di analisi, il
rigeneratore si comporta in modo ideale, mentre il gas che passa
dall''heater al volume di espansione si troverà alla temperatura di parete Prove sperimentali al variare della pressione del fluido di lavoro 73 dell''heater. Effetti analoghi si ottengono nel caso di passaggio del fluido in
senso opposto.
Secondo tale approccio il calore scambiato rispettivamente nel cooler e
nell''heater deriva unicamente dalla variazione di temperatura che il fluido
di lavoro subisce durante le trasformazioni adiabatiche che si verificano
all''interno delle celle di compressione ed espansione. Il flusso termico
scambiato con l''ambiente risulta pari, istante per istante, alla potenza
necessaria per riscaldare/raffreddare la portata di fluido entrante nella
rispettiva cella fino alla temperatura condizionale specifica, a partire dalla
temperatura antecedente l''ingresso.

Risulta facilmente intuibile come l''analisi adiabatica non consideri gli
importanti effetti relativi alla non idealità del rigeneratore e alle
irreversibilità di scambio termico che si verificano all''interno dell''heater e
del cooler. Essa permette tuttavia di valutare in prima approssimazione
l''andamento di diverse grandezze termodinamiche durante il ciclo, nonché
dei fluissi massici che si instaurano tra le varie celle. La risoluzione del
modello Adiabatico avviene mediante l''applicazione delle equazioni di
conservazione della massa e dell''energia tra un volume di controllo e il
successivo, nel rispetto dell''equazione di continuità relativa all''intero
sistema. Si ottiene così un sistema di 22 variabili in 16 incognite, come
riportato in Figura 3.4.

Le equazioni così ottenute non sono esplicitamente integrabili e il sistema
diviene non lineare, richiedendo il ricorso a metodi di soluzione numerici.
Il ciclo è quindi discretizzato in una serie di intervalli (tipicamente 36),
ciascuno dei quali è individuato da una posizione iniziale e una posizione
finale dell''albero motore. All''interno di ciascun intervallo di
discretizzazione il sistema viene studiato come se si trattasse di una
condizione stazionaria, trascurando gli effetti inerziali e di transitorio legati
alla variazione di velocità e condizioni termodinamiche del fluido. Note
tutte le condizioni al contorno che caratterizzano il motore, la soluzione del
sistema si riconduce alla soluzione di un problema ai valori iniziali,
costituiti dalle temperature del fluido contenuto nelle celle di compressione
e di espansione per l''angolo di manovella iniziale (θ=0°). Il problema si
considera risolto quando al termine di un manovellismo completo (θ
=360°) le temperature del fluido contenuto nelle due celle menzionate
risultano pari a quelle iniziali, a meno di un margine di tolleranza
predefinito. Tale vincolo equivale all''imposizione della condizione di
stazionarietà di funzionamento del motore, ossia all''assenza di fenomeni di
accumulo termico all''interno delle diverse celle da un ciclo motore al
successivo. Capitolo 3
74
Figura 3.4 Sistema di equazioni del modello ideale adiabatico [12].
La simulazione mediante il modello adiabatico del funzionamento di
motori operanti in condizioni standard porta all''individuazione di
rendimenti di secondo principio prossimi al 90%. ' evidente come tale
valore non sia rappresentativo delle reali prestazioni del ciclo, che
corrispondono generalmente a rendimenti di secondo principio compresi
tra il 30 e il 60%. Pertanto il modello adiabatico di Urieli non risulta
soddisfacente ai fini della descrizione del comportamento di un motore
Stirling reale.
Prove sperimentali al variare della pressione del fluido di lavoro 75 Nel seguito sono ricapitolate le principali cause di tale scostamento:
' Ipotesi di comportamento ideale del rigeneratore che risulta in grado di fornire/assorbire dal fluido di lavoro tutto il calore necessario
a realizzare il riscaldamento/raffreddamento dalla temperatura del
cooler a quella dell''heater, e viceversa; ' Ipotesi di idealità dei processi di introduzione e rimozione del calore, con conseguente isotermicità del fluido di lavoro rispetto alla
temperatura di parete di heater e cooler; ' Ipotesi di assenza delle perdite di carico che si originano per effetto degli attriti viscosi all''interno del fluido, fenomeno particolarmente
marcato all''interno del rigeneratore che è caratterizzato da sezioni di
passaggio ridotte e percorsi tortuosi; ' Mancata valutazione delle perdite per conduzione sulle pareti che danno luogo a un flusso dissipato dall''heater verso il cooler, e dal
cooler verso l''ambiente; ' Mancata valutazione delle perdite per attrito meccanico a carico del cinematismo e delle tenute atte a garantire l''assenza di trafilamenti del
fluido di lavoro verso l''ambiente; ' Valutazione di tutti i fenomeni in regime stazionario, trascurando tutti gli effetti dinamici e di transitorio; ' Ipotesi di adiabaticità delle celle di compressione ed espansione, nelle quali si impone che il fluido non scambi calore con le pareti; ' Assegnazione di una pressione istantaneamente uniforme in tutto il motore, determinata in funzione del volume e della temperatura
istantanea di ogni cella, trascurano le disuniformità conseguenti la
presenza delle perdite di carico.

3.1.2 Analisi semplificata di secondo ordine
Come già accennato nei paragrafi precedenti, l''analisi semplificata di
secondo ordine si pone come step successivo rispetto all''analisi adiabatica
di Urieli ed è stata sviluppata al fine di fornire una valutazione dei seguenti
fenomeni:
Capitolo 3
76
' Comportamento termico non ideale del rigeneratore; ' Irreversibilità di scambio termico all''interno delle celle di heater e cooler; ' Perdite di carico concentrate e distribuite all''interno delle celle di heater, cooler e rigeneratore; ' Perdite per conduzione dall'' heater verso il cooler attraverso le pareti e le camicie del motore.
Purtroppo la maggior parte degli studi riguardanti i fenomeni di scambio
termico riportati in letteratura si riferiscono a condizioni di flusso
stazionario, ben lontane dalla condizione di flusso oscillante a velocità
variabile riscontrabile all''interno del motore Stirling. La modellizzazione si
basa quindi, nella maggior parte dei casi, su un approccio di tipo quasi-
stazionario analogo a quello utilizzato relativamente all''analisi adiabatica.
L''analisi semplificata di secondo ordine prevede che per ciascun intervallo
di discretizzazione del ciclo siano assegnati dei valori medi di velocità del
flusso, rappresentativi dell''effettiva distribuzione di velocità all''interno
dell''intervallo stesso. Tali valori medi sono successivamente utilizzati per
la determinazione dei coefficienti di scambio termico e di attrito viscoso
caratteristici di ciascuna sezione del motore e di ciascun intervallo
temporale di discretizzazione.

L''approccio quasi-stazionario rappresenta un''assunzione forte nello studio
del funzionamento del ciclo Stirling poiché tende nella maggior parte dei
casi a sovrastimare i coefficienti di scambio termico e sottostimare quelli di
attrito viscoso, portando all''individuazione di rendimenti ottimistici
rispetto alle reali prestazioni del ciclo. Le possibili alternative a tale
approccio sono due:
' Valutazione dei coefficienti di scambio termico e di attrito viscoso mediante correlazioni derivanti da studi specifici condotti in
condizioni di flusso oscillante,
per geometrie analoghe rispetto a
quella oggetto di studio; ' Impiego di metodi CFD per la determinazione del reale comportamento del fluido all''interno degli scambiatori e, più in
generale, di tutte le sezioni del motore caratterizzate dalla presenza di
fenomeni di attrito viscoso o di scambio termico. Prove sperimentali al variare della pressione del fluido di lavoro 77 Il secondo metodo è più generale rispetto al primo in quanto consente di
simulare motori caratterizzati da architetture anche molto diverse tra loro,
spaziando sia in termini di geometria costruttiva che di tipologia di
materiali e conformazione delle sezioni di scambio termico. Esso risulta
però scarsamente conciliabile con l''approccio a parametri concentrati
scelto nel presente studio e richiede tempi lunghi sia in termini di messa a
punto del modello di analisi, sia di tempi di calcolo, risultando
difficilmente conciliabile con le tempistiche richieste per un confronto
parallelo rispetto allo sviluppo delle prove sperimentali.
Si è quindi optato per la prima alternativa che, in relazione alla specifica
tipologia di rigeneratore presente all''interno del motore oggetto di studio,
risulta attuabile grazie agli studi condotti da Gedeon e Wood nel 1996
presso le strutture della NASA [18] e finalizzati allo studio delle
prestazioni di diverse tipologie di rigeneratore per applicazioni Stirling in
condizioni di flusso oscillante.

3.2 Prove sperimentali

Il WhisperGen EU1 è una macchina in grado di lavorare in diverse
condizioni operative, caratterizzate da diversi valori di portata e
temperatura di ritorno dell''acqua, e come si è già visto i maggiori benefici
si hanno alle basse temperatura dell''acqua, sfruttando la condensazione dei
fumi. Per questi test, si è deciso di rimanere nelle condizione nominali
dell''unità cogenerativa, ovvero 0.194 Kg/s e 50°C per la temperatura
dell''acqua di ritorno, mentre per quanto riguarda il settaggio delle
pressione dell''azoto si sono scelti i valori di 20, 16, 12 e 8 bar.

Il motore Stirling è stato installato nella cella di prova del laboratorio in
modo analogo a quanto sarebbe eseguito se il motore fosse istallato in una
unità abitativa reale. Durante le prove le grandezze misurate e controllate
sono: portata massica e temperatura dell''acqua di ritorno alla macchina; la
prima tramite il numero di giri della pompa e la seconda tramite un
controllore Proporzionale-Derivativo-Integrale (PID) che regola la
dissipazione nel circuito di raffreddamento. Inoltre vengono misurate, ma
non controllate:
' portata massica e composizione del gas naturale; ' temperatura ed umidità dell''aria ambiente; Capitolo 3
78
' temperatura dell''acqua di mandata della macchina; ' potenza elettrica generata; ' temperatura, composizione ed emissioni (CO, NO, NO2, SO2 e tenore di O2) nei fumi; ' quantità di condensa dai fumi.
Figura 3.5 Schema di misurazione [13].
Prove sperimentali al variare della pressione del fluido di lavoro 79 L''apporto di combustibile è garantito attraverso una linea di gas naturale
pressurizzato a 10 bar da un compressore a vite e poi laminata, all''interno
della cella, a 2 bar per permettere una misura di portata massica accurata
per mezzo del Bronkhorst IN-FLOW 112-AI. La composizione è misurata
tramite micro-gascromatografo Pollution VEGA-GC, che viene
periodicamente tarato con miscele preparate gravimetricamente da Sapio
S.p.a. La composizione e le emissioni nei fumi sono misurate anidre a
seguito di una deumidificazione per mezzo di un frigorifero ABB Sample
gas cooler Advance SCC-C. La composizione molare è misurata con il
micro-gascromatografo precedentemente citato. Le condizioni dell''aria
comburente aspirata dall''ambiente sono misurate da una termo resistenza
Scandura Pt100 1/3 DIN e da un semplice e portatile termometro Testo H1
in grado di rilevare anche l''umidità relativa.

La potenza elettrica scambiata dal motore con la rete (assorbita in alcuni
momenti e generata in molti altri), la sua tensione e l''angolo di fase sono
rilevati tramite un wattmetro Fluke Norma 4000 accoppiato con
trasformatori amperometrici collegati a Signaltec MCTS.
La temperatura dell''acqua in ingresso al motore, come anticipato, è
controllata tramite PID implementato sul PLC del laboratorio. La
temperatura in ingresso e quella in uscita sono misurate con termo
resistenze Scandura Pt100 1/3 DIN, la portata massica con Emerson
Micromotion Coriolis Elite CMF025 e, infine, la pressione differenziale tra
le due sezioni tramite Emerson Rosemount 2051 CF.

Le condense dei gas combusti vengono misurate tramite la bilancia Bel
Mark 1700 (bilancia che verrà sostituita in futuro con una cella di carico).

Tutte queste grandezze misurate, vengono successivamente elaborate per
verificare il bilancio energetico delle singole prove; il test si ritiene valido
se la differenza tra energia in ingresso ed energia in uscita è minore di una
tolleranza predefinita.


3.2.1 Prove in stato stazionario
Per lo svolgimento delle prove a diversa pressione dell''azoto, si è deciso di
considerare le misurazioni una volta raggiunta la condizione di regime;
dopo aver impostata la quantità di acqua da far circolare, regolando il
numero di giri della pompa DAB, il sistema si considera in stato
stazionario quando la potenza sviluppata e la temperatura dell''acqua di
ritorno si stabilizza al valore desiderato, ammettendo al più una variazione
di ± 0.3 °C dell''acqua. Capitolo 3
80
Come già accennato nelle pagine precedenti, la tempera dell''acqua di
ritorno viene regolata e stabilizzata al valore desiderato tramite un
controllore PID (Proporzionale Integrale Derivativo). Durante le prime
prove abbiamo riscontrato diversi problemi con il controllore stesso, che ci
impediva di raggiungere una situazione di regime stabile. Le variazioni di
temperatura erano molto maggiori della tolleranza ammessa, pertanto si è
dovuto procedere ad una nuova taratura dei parametri propri del controllore
Kp, Ki e Kd. Questo ha richiesto una notevole quantità di prove e di tempo, ma alla fine si è riuscito a trovare una nuova terna che non solo
permettesse di rientrare nelle tolleranze richieste, ma anche di mantenere
una temperatura dell''acqua di ritorno molto più costante e precisa di
quanto si potesse ottenere in precedenza. Di seguito si riporta a titolo di
esempio il grafico della temperatura dell''acqua di ritorno della prova a 8
bar, in cui si può notare come le variazioni della temperatura rispetto al
valore desiderato sono molto piccole (circa ± 0.1 °C contro i ± 0.3 °C
ammissibili). Si ricorda che il valore di temperatura impostato è di 50°C:

Figura 3.6 Andamento della temperatura dell''acqua di ritorno (8 bar).

All''avvio dell''impianto viene imposto un determinato numero di giri e
quindi una certa portata volumetrica, in modo che la portata massica letta
dal misuratore corrisponda con quella richiesta dalla prova. Il
riscaldamento dell''acqua del circuito ha come conseguenza la diminuzione
della viscosità dell''acqua e quindi delle perdite di carico distribuite lungo 49,75 49,80 49,85 49,90 49,95 50,00 50,05 50,10 50,15 17:38 17:48 17:58 18:08 18:18 Temperatura a regime Temperatura Prove sperimentali al variare della pressione del fluido di lavoro 81 le tubazioni; in condizioni di regime termico è quindi necessario
modificare il numero di giri della pompa per avere la portata massica
desiderata.

Una volta raggiunto la condizione di regime, viene avviata l''acquisizione
dei dati tramite il programma di controllo LabVIEW del laboratorio, dal
quale è possibile inoltre verificare l''andamento nel tempo della
temperatura dell''acqua. Per quanto riguarda la durata della prova, prima si
considerava un intervallo di tempo pari a circa 4 cicli di oscillazione della
temperatura dell''acqua di ritorno, per un massimo di 40-50 minuti.
Considerando la maggiore stabilità ottenuta durante queste ultime prove, si
è scelto di considerare la prova conclusa dopo circa 30 minuti dall''avvio.
Durante i test si è osservata una caratteristica diminuzione del consumo di
gas naturale più o meno ogni 15 minuti, per la durata di circa 10 secondi,
in corrispondenza delle quali è possibile osservare anche la diminuzione
della potenza elettrica prodotta e delle emissioni; la presenza di questo
fenomeno non è citata sul manuale del WhisperGen ma ha una scarsa
influenza sulle prestazioni, in quanto interessa un intervallo di tempo molto
breve.

Per eseguire i bilanci energetici delle varie prove, sono stati usati come
valori indicativi di ogni grandezza misurata il proprio valor medio durante
l''intero periodo di test.


3.2.2 Procedura di carico / scarico dell''azoto
Lo scopo di variare la pressione del fluido di lavoro è quello di modificare
la quantità di azoto all''interno dell''unità Stirling. Fissata una temperatura,
infatti, una abbassamento della pressione comporta un diminuzione di
azoto all''interno del motore, viceversa un aumento di pressione comporta
una maggiore quantità di gas. ' evidente che anche la temperatura giochi
un ruolo decisivo, poiché agisce anch''essa sulla pressione misurata,
pertanto risulta importante scegliere una condizione nominale alla quale
fare riferimento per le variazioni di pressione. In tal senso, il manuale del
WhisperGen [7] riporta come la pressione dell''azoto all''interno del motore
sia di 20 bar @ 20°C (o in alternativa 24 bar @ 70°C); si è scelto quindi di
tenere questa condizione come quella di riferimento.

Durante le operazione di carica / scarica dell''azoto, quindi, è stato
importante mantenere l''unità ad una temperatura più vicina possibile ai
20°C; per fare ciò si è sfruttata l''acqua di torre, che è ad una temperatura di
circa 18-22°C, e non subisce grosse variazioni durante l''anno. In pratica si
procedeva a far circolare acqua all''interno del motore Stirling fino a Capitolo 3
82
quando la sua temperatura in uscita dall''unità fosse sostanzialmente uguale
alla temperatura in ingresso; in questo modo si raffredda l''unità
WhisperGen fino ad una temperatura prossima a quella di riferimento.
Dopodiché è possibile modificare correttamente la pressione dell''azoto.

Per le operazioni di carica / scarica, è stato necessario utilizzare un
manometro fornito direttamente dalla casa madre, che si collega al punto di
ricarica del motore Stirling.

Figura 3.7 Manometro [7].

La procedura di scarica consiste in:
' Rimuovere il tappo dal punto di ricarica sul motore ' Assicurarsi che le valvole B e D siano chiuse ' Collegare il dado A al punto di ricarica Prove sperimentali al variare della pressione del fluido di lavoro 83 ' Aprire la valvola D e rilevare la pressione dal manometro ' Aprire la valvola B lentamente fino al valore di pressione desiderato
Figura 3.8 Punto di ricarica dell''unità Stirling. Capitolo 3
84
Nel caso in cui si debba aumentare la pressione dell''azoto, si può eseguire
la medesima procedura indicata sopra, con l''unica differenza che la
connessione C deve essere collegata ad una bombola di azoto (la linea di
azoto della cella di prova può eventualmente arrivare fino a 18 bar).


3.2.3 Bilancio energetico
Un importante controllo che viene effettuato per validare le prove eseguite
è quello sul bilancio di energia, che idealmente dovrebbe permettere di
osservare l''uguaglianza tra i flussi energetici entranti e quelli uscenti. Il
bilancio è realizzato su base molare e non su base massica, poiché le
composizioni dei campioni analizzati dal micro gascromatografo sono
fornite su base molare.

La portata molare del gas naturale è calcolata a partire dalla portata
massica attraverso i coefficienti C1 e C2:
' ' (3.1)
dove C1 è un coefficiente di conversione da base massica a base molare e da minuti a secondi:

(3.2)
mentre C2 è un coefficiente correttivo per il misuratore di portata Bronkhorst, tarato sulla misura di metano puro. Il coefficiente C2 è calcolato tenendo conto della frazione molare, xi,gn [%], e del coefficiente correttivo, Bi, di ogni componente del gas naturale e del coefficiente correttivo del gas naturale utilizzato per la taratura dello strumento, Btar, (nel nostro caso, essendo il gas di taratura metano, Btar = BCH4) indicati nel manuale del misuratore Bronkhorst e riportati nella tabella di seguito.
Composto Bi CH4 0.76 C2H6 0.49 C3H8 0.34 C4H10 0.25 C5H12 '' C6H14 0.49 Prove sperimentali al variare della pressione del fluido di lavoro 85 CO2 0.74 N2 1 Tabella 3.1 Coefficienti correttivi Bi per il misuratore di portata Bronkhorst.
Il coefficiente C2 risulta quindi:
('' ) (3.3)
Il Potere Calorifico Superiore del gas naturale viene calcolato utilizzando
la funzione implementata nel modello VBA descritto nel capitolo
precedente.

Per considerare i cambiamenti della composizione e della portata dei
prodotti di combustione al camino al variare della portata di condensa,
sono stati denominati ''gas combusti' i prodotti formatasi dalla
combustione e ''fumi' ciò che transita nel camino, al netto della condensa
formatasi. Si definisce β il rapporto molare tra gas combusti e
combustibile:

'' ( ) (3.4)
Nota quindi la portata molare di gas naturale, è possibile ricavare la portata
molare di gas combusti:

' ' (3.5)
pari alla somma della portata molare dei fumi umidi e della portata molare
di condensa:

' ' ' (3.6)
dove ' viene calcolato a partire dalla misura della massa condensata attraverso la bilancia BEL Mark 1700:

' (3.7) Capitolo 3
86
La portata di fumi umidi risulta quindi essere:

' ' ' (3.8)
La portata molare di acqua presente nei gas combusti, formatasi in parte
dalla combustione e in parte proveniente dall''aria comburente, è pari a:

' ' '' ( ( ) ) (3.9)
Avendo calcolato in precedenza la portata molare di condensa è possibile
ricavare la portata molare di acqua nei fumi, pari a:

' ' ' (3.10)
che permette di calcolare la frazione molare di acqua nei fumi:

' ' (3.11)
attraverso cui si può ricavare la composizione molare dei fumi umidi dalla
misura gascromatografica dei fumi secchi:

( ) (3.12)
Definito il rapporto molare stechiometrico aria / combustibile:

'' (3.13)
è quindi possibile calcolare il rapporto molare aria / combustibile:

( ) (3.14)
Nota la portata molare di gas naturale si ricava la portata di aria
comburente:

' ' (3.15) Prove sperimentali al variare della pressione del fluido di lavoro 87 Precedentemente il bilancio energetico veniva effettuato esplicitando le
potenze scambiate e verificando che:

{ (3.16)
In seguito è stato deciso un approccio diverso, valutando i flussi entalpici
in ingresso e uscita e verificando la loro uguaglianza:

{ '' ' '' ' (3.17)
La prova è ritenuta valida quando la differenza tra il flusso entalpico in
ingresso e il flusso entalpico in uscita, risulta inferiore ad una tolleranza
massima, in questo caso pari al 4%.

| | ( ) (3.18)
Una volta costatato che il bilancio energetico rispettasse la tolleranza
massima indicata, è stato possibile valutare effettivamente le potenze
termiche scambiate. La potenza elettrica è misurata direttamente dal
wattmetro Fluke Norma 4000, di conseguenza il valore preso in
considerazione è la media dei dati acquisiti sull''intera durata della prova.
Per quanto riguarda la potenza termica scambiata, anch''essa effetto utile
del ciclo termodinamico, viene calcolata come segue:

' ' ( ( ) ( )) (3.19)
' quindi possibile calcolare il rendimento elettrico, termico e di primo
principio delle diverse prove:

' (3.20)
'
' (3.21) Capitolo 3
88
' ' (3.22)
Dove ' rappresenta la potenza termica entrante attraverso la combustione, riferita al Potere Calorifico Superiore (PCS) poiché si sfrutta
la condensazione dei fumi:

' ' (3.23)

3.2.4 Risultati prove sperimentali
Di seguito sono riportati i risultati delle prove sperimentali al variare della
pressione del fluido di lavoro. Tali valori sono riferiti all''intera unita
cogenerativa, considerando il calore entrante come definito dall''equazione
(3.23). Le stime delle proprietà termodinamiche di tutte le specie coinvolte
sono state eseguite utilizzando il modello numerico descritto nel Capitolo
2, utilizzando le formule descritte precedentemente per il calcolo delle
potenze scambiate.

Si può notare come la potenza elettrica netta sviluppata aumenti
all''aumentare della pressione dell''azoto, viceversa si ha un aumento della
potenza termica al diminuire della pressione dell''azoto.

8 bar 12 bar 16 bar 20 bar ' [kW] 8.68 9.13 8.99 9.52 ' [ ] 4.54 5.10 5.41 6.03 ' [kW] 7.63 7.88 7.57 8.15 [kW] 0.374 0.611 0.745 0.885 [%] 4.31 6.70 8.29 9.29 [%] 87.96 86.32 84.19 85.56 [%] 92.27 93.02 92.47 94.86 Tabella 3.2 Risultati prove sperimentali.

Un abbassamento della pressione dell''azoto corrisponde ad una minor
quantità di fluido di lavoro nel ciclo termodinamico, e, poiché i volumi Prove sperimentali al variare della pressione del fluido di lavoro 89 massimi e minimi sono imposti dal cinematismo dell''unità Stirling, si
avranno anche pressioni massime e minime inferiore rispetto alla
condizione di design, e di conseguenza una diminuzione della potenza
elettrica sviluppata. Si ricorda, a tal proposito, come il lavoro specifico
prodotto dal ciclo termodinamico, sia equivalente all''area interna al ciclo
nel diagramma pressione-volume.

Si nota, inoltra, come al diminuire della pressione dell''azoto, aumenti il
rendimento termico dell''unità cogenerativa; questo perché la capacità
termica del fluido di lavoro diminuisce, a causa di una sua minore massa,
determinando una capacità di assorbimento di calore dall''heater inferiore.
Di conseguenza una maggiore porzione di calore entrante viene disperso
tramite le pareti dell''heater, e trasferito all''acqua cogenerativa. In termini
assoluti la potenza termica dell''unità cogenerativa tende a diminuire con la
pressione del fluido di lavoro, ma, come si vede dalla Tabella 3.2., anche il
calore in ingresso diminuisce, e in maniera maggiore rispetto alla
diminuzione di potenza termica rilevata, determinando, globalmente, un
aumento del rendimento termico.

Si ricorda infine che tutti i rendimenti sono calcolati seguendo la procedura
descritta nel paragrafo precedente, considerando il calore entrante
attraverso il combustibile nell''unità cogenerativa ' .

3.3 Confronto tra prove sperimentali e modello numerico

Per poter effettuare il confronto tra le prove sperimentali ed il modello
numerico sviluppato dall''Ing. Fergnani, è necessario un approccio
leggermente diverso rispetto a quanto descritto nel Capitolo 3.2.3. Difatti il
modello numerico riguarda la sola unità Stirling, non l''intera unità
cogenerativa, pertanto è essenziale poter valutare nella maniera più
accurata possibile, la quantità di calore entrante sulle teste dei cilindri
(heater) rispetto al totale entrante nello WhisperGen. Questo perché la
restante quota parte non contribuisce alla produzione di potenza elettrica,
non essendo in realtà calore entrante nel ciclo termodinamico. La
quantificazione della potenza termica entrante nel ciclo Stirling può quindi
essere effettuata mediante due modalità di calcolo:
' Come somma del calore rimosso al cooler e nel circuito di raffreddamento del generatore (ceduto all''acqua), del calore ceduto
all''ambiente dalle pareti del motore e della potenza elettrica lorda Capitolo 3
90
prodotta, ottenuta a sua volta come somma di potenza elettrica netta
e assorbimento degli ausiliari; ' Per differenza considerando come input energetico la somma del calore entrante sotto forma di potere calorifico superiore del gas
naturale inviato al combustore, del calore sensibile posseduto
dall''aria in ingresso al combustore (a seguito dei pre-riscaldamenti)
e del calore latente legato all''umidità presente nell''aria. A tali
fattori devono essere sottratti tutti i flussi termici ad esclusione di
quelli elencati in precedenza. Figura 3.9 Flussi energetici dell''unità WhisperGen EU1. Prove sperimentali al variare della pressione del fluido di lavoro 91 Come valore indicativo della potenza termica entrante, si è considerato il
valore medio derivante dai due metodi descritti sopra.

Di seguito si riporta il confronto tra il modello sperimentale e il modello
numerico in termini di quantità di calore in ingresso nel ciclo Stirling,
potenza elettrica netta prodotta e rendimento elettrico. Si ricorda che tali
valori di rendimento non sono riferiti al flusso energetico
complessivamente entrante nell''unità cogenerativa sotto forma di PCS del
combustibile, ma alla sola potenza termica entrante nel motore Stirling.

8 bar 12 bar 16 bar 20 bar Prove sperimentali [kW] 8.68 9.13 8.99 9.52 [kW] 4.54 5.10 5.41 6.03 Var. % VS 20 bar -25% -15% -10% 0% [K] 682 639 566 566 [K] 55 56 56 56 [kW] 0.374 0.611 0.745 0.885 Var. % VS 20 bar -44% -18% -6% 0% 8.2% 12.0% 13.8% 14.7% Var. % VS 20 bar -44% -18% -6% 0% Simulazione numerica [kW] 4.31 4.96 5.31 6.15 Var. % VS 20 bar -30% -19% -14% 0% [kW] 0.446 0.689 0.781 0.993 Var. % VS 20 bar -55% -31% -21% 0% 10% 14% 15% 16% Var. % VS 20 bar -36% -14% -9% 0% Tabella 3.3 Confronto prove sperimentali - modello numerico. Capitolo 3
92
Figura 3.10 Confronto potenze prove sperimentali - modello.

Com''è logico aspettarsi, dalla Figura 3.10 si nota come al diminuire della
pressione del fluido di lavoro, si ha un decremento della potenza elettrica
netta prodotta. Difatti ad un abbassamento della pressione dell''azoto
corrisponde ad una minor quantità di fluido di lavoro nel ciclo
termodinamico, e poiché i volumi massimi e minimi sono imposti dal
cinematismo dell''unità Stirling, si avranno anche pressioni massime e
minime inferiore rispetto alla condizione di design, e quindi una
diminuzione della potenza elettrica sviluppata.

Inoltre è interessante notare come ad un abbassamento della pressione
operativa, corrisponda una minore quantità di calore entrante nell''unità
Stirling; questo fenomeno è l''effetto principale della minore capacità
termica di assorbire calore dall''heater già citata in precedenza. A conferma
di ciò, in Figura 3.12, si nota come la temperatura misurata di una delle
teste del motore Stirling (che corrisponde sostanzialmente all''heater),
aumenti al diminuire della pressione. 0 1 2 3 4 5 6 7 8 12 16 20 P ot e n za [ kW ] Pressione azoto [bar] Pel sperimentale Pel simulazione Q ingresso sperimentale Q ingresso simulazione Prove sperimentali al variare della pressione del fluido di lavoro 93 Figura 3.11 Confronto potenza elettrica prove sperimentali - modello.
Figura 3.12 Calore entrante nel ciclo Stirling. 0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 8 12 16 20 P ot e n za [ kW ] Pressione azoto [bar] Pel sperimentale Pel simulazione 4,00 5,00 6,00 7,00 8,00 400 450 500 550 600 650 700 750 800 8 12 16 20 Ca lor e e n tr an te n e l ci cl o Sti rl in g [kW ] Te m p e ra tu ra a le tt e He at e r C] Pressione di carica [bar] Capitolo 3
94
Bisogna inoltre ricordare che la potenza elettrica prodotta, non è
influenzata solo dalla pressione del fluido di lavoro, ma anche dalle
temperature del ciclo stesso. Difatti la diminuzione della potenza elettrica
al diminuire della pressione dell''azoto è in parte mitigata dall''aumento di
temperatura dell''heater; ci si aspetta che questo fenomeno si presenti, in
maniera invertita, all''aumentare della pressione. In questo caso, infatti, si
dovrebbe osservare anche una diminuzione della temperatura dell''heater
dovuta ad una maggiore capacità termica dell''azoto.

In figura 3.13 si riportano i dati relativi ai rendimenti elettrici misurati
sperimentalmente e quelli previsti dal modello numerico

Figura 3.13 Confronto rendimento elettrico.
Si nota coma alle basse pressione operative, si ha una maggiore incidenza
delle perdite fisse (attriti meccanici e perdite per il consumo degli
ausiliari).

La potenza elettrica rilevata sperimentalmente risulta mediamente inferiore
di circa il 12%; tuttavia, la prova effettuata nelle condizioni nominali di 20
bar, ha dato risultati inferiore, in termini di potenza sviluppata, anche a 0,0% 5,0% 10,0% 15,0% 20,0% 25,0% 30,0% 8 12 16 20 R e n d im e n to e le ttr ico [ %] Pressione azoto [bar] Rendimento elettrico
sperimentale Rendimento elettrico
simulazione Prove sperimentali al variare della pressione del fluido di lavoro 95 tutte le precedenti prove svolte a 20 bar. Pertanto si ritiene che l''unità
cogenerativa possa essere stata in qualche modo danneggiata, o rimontata
in malo modo, dalle ripetute operazioni di montaggio e smontaggio svolte
durante questi anni di test, ed è quindi necessario effettuare ulteriori prove
nelle 4 diverse condizioni di pressione per accertarsi dei valori misurati.


3.4 Propagazione degli errori

Poiché l''apparato strumentale per la misurazione di tutti i parametri
necessari all''esecuzione del bilancio energetico dell''unità cogenerativa,
non è cambiato rispetto le precedenti prove, si rimanda, per la
determinazione delle incertezze, al Capitolo 6 della tesi di E. Zattoni ''
Prove Sperimentali su un Sistema Micro-Cogenerativo basato
su ciclo Stirling'.

Si riportano comunque una stima degli errori, relativi alla prova in
condizione nominale a 20 bar:

= 0.885 ± 0.002 [kW] ' = 8.15 ± 0.15 [kW] = 9.29 ± 0.26 [%] = 85.56 ± 2.66 [%]
Capitolo 3
96
97 Capitolo 4

Conclusioni e sviluppi futuri

Si riportano, di seguito, le osservazioni conclusive sulle attività svolte:
' Per quanto riguarda la bontà del modello numerico per il calcolo delle proprietà termodinamiche, descritto nel Capitolo 2, ci si può
ritenere soddisfatti sia in termini di facilità di utilizzo, che in
termini di errori numerici, infatti si è riusciti ad chiudere i bilanci
energetici con uno scarto mediamente inferiore a quanto era
possibile in precedenza. Questo è probabilmente dovuto alla
migliore valutazione delle proprietà termodinamiche, in primo
luogo dell''entalpia, delle miscele bifase.
' ' stata elaborata una nuova procedura per le operazione di carica / scarica dell''azoto, che potrà essere utilizzata, in seguito, anche
come procedure per il cambio del fluido di lavoro all''interno
dell''unità Stirling.
' Le prove sull''unità Stirling cogenerativa hanno dato i risultati ipotizzati: si nota una diminuzione della potenza elettrica e del
rendimento elettrico al diminuire della pressione del fluido di
lavoro. Viceversa il rendimento termico aumenta al diminuire della
pressione del fluido di lavoro. La quantità di combustibile, e quindi
di potenza termica in ingresso all''unità cogenerativa, invece, è
determinata da una logica di controllo propria della macchina.
' Il confronto tra le prove sperimentali ed il modello numerico dell''unità Stirling, ha dato risultati concordi. La previsione sulla
potenza termica entrante nel motore Stirling è coerente con i
risultati sperimentali, mentre la potenza elettrica netta misurata e
quella prevista dal modello numerico, risultano meno concordi
numericamente ma presentano il medesimo andamento al variare
della pressione operativa. La differenza, come già precedentemente
accennato, potrebbe essere causata da un danneggiamento dello
WhisperGen durante le numerosi fasi di montaggio / smontaggio in
questi anni di prove.
Capitolo 4
98
4.1 Sviluppi futuri

Poiché è stato possibile effettuare una singola prova per ogni pressione
operativa, nell''immediato futuro sarà necessario ripetere più volte le prove
al fine di verificarne la ripetibilità e la correttezza dei dati rilevati,
soprattutto per la potenza elettrica sviluppata che come già citato
precedentemente è stata inferiore rispetto alle precedenti prove a pari
condizioni operative. Una volta verificata la ripetibilità, verranno eseguito
prove sperimentali a pressioni maggiori di quella di design.

In secondo luogo si eseguiranno prove sperimentali con l''unità
cogenerativa in una nuova configurazione. L''obiettivo sarà quello di
preriscaldare l''aria comburente, effettuando un ricircolo di energia
posseduta dai fumi verso l''ingresso del motore Stirling. In questo modo, a
parità di energia entrante nel sistema, si diminuisce la frazione di energia
direttamente imputabile al combustibile, aumentando di conseguenza il
rendimento globale del sistema.

' ' ' ' ( ) ' (4.1)
Come si vede dall''equazione (4.1), a pari energia termica entrante nel
ciclo, Qin, un aumento della temperatura dell''aria comburente, e quindi dell''ossigeno, comporta una minore portata di combustibile necessaria , mc, e quindi una minore energia chimica in ingresso. Un ulteriore vantaggio è
dato dalla maggiore temperatura dei gas inerti in camera di combustione,
presenti in gran quantità nell''aria ambiente, e quindi un minor calore
assorbito da questi ultimi per raggiungere la temperatura adiabatica di
fiamma. A livello pratico, il preriscaldo dell''aria comburente, è reso
possibile mediante l''utilizzo di uno scambiatore aria/fumi a fascio tubiero,
posto a monte dello scambiatore secondario e a valle della camera di
combustione. La scelta del preriscaldatore di aria comburente, date le
piccole portate di gas combusti e aria, è stata vincolata dalla disponibilità
del mercato, in quanto poche aziende hanno interesse nel produrre
scambiatori per alte temperature di così piccole dimensioni. La scelta,
quasi obbligata, è quindi ricaduta su di uno scambiatore della ''LU-VE
S.p.a.' a fascio tubiero alettato e flusso incrociato, che in condizioni di
design permette il preriscaldo dell''aria fino ad una temperatura di circa
320°C e di raffreddare i gas combusti fino ad una temperatura stimata di
106 °C. La configurazione impiantistica finale, con l''introduzione dello
scambiatore recuperativo, prevede quest''ultimo a valle della camera di
combustione, ed a monte dello scambiatore secondario, il quale è Conclusioni e sviluppi futuri 99 comunque necessario per raffreddare ulteriormente i fumi e ottenere elevati
rendimenti di primo principio dell''unità cogenerativa. L''inserimento del
preriscaldo dell''aria ha introdotto notevoli problematiche sul layout
impiantistico, le principali delle quali sono state la scelta del collocamento
del miscelatore aria/combustibile e del ventilatore aspirante aria
dall''ambiente. Nella configurazione standard, questi due componenti, sono
posti entrambi a monte della camera di combustione, con il miscelamento
che avviene prima dell''ingresso nella stessa. Mantenere la miscelazione
aria-combustibile immediatamente a valle del ventilatore, comporterebbe il
rischio di autoaccensione della miscela di combustibile all''interno dello
scambiatore recuperativo, poiché la miscela raggiungerebbe una
temperatura di circa 320°C. In realtà, con il rapporto aria/combustibile
usato, la temperatura di autoaccensione è di circa 500°C, quindi più elevata
della temperatura raggiungibile, ma rimane il rischio di contatto con punti
caldi dello scambiatore, pertanto si è deciso di porre il miscelatore a valle.
Il ventilatore, viceversa, avendo la girante in materiale plastico, è stato
lasciato a monte dello scambiatore; inoltre, dal punto di vista energetico,
risulterebbe più dispendioso ''comprimere' la miscela calda piuttosto che la
sola aria alla temperatura ambiente. Di seguito sono riportati i principali
dati di progetto dello scambiatore scelto: Fluido Esterno ai Tubi: Aria Portata Massica: 0.005 kg/s Temperatura (Inlet / Outlet): 25.000 / 323.317 °C Velocità (std / reale): 0.017 / 0.026 m/s Perdita di pressione: 0.014 Pa Fluido Interno ai Tubi: Gas Combusti Portata Massica: 0.005 kg/s Temperatura (Inlet / Outlet): 400.000 / 106.262 °C Velocità (std / reale): 9.634 m/s Perdita di pressione: 1.028 Pa Proprietà: Potenza termica scambiata: 1.521 kW Tabella 4.1 Principali caratteristiche scambiatore LU-VE. Capitolo 4
100
101 Bibliografia

[1] R. Cremonesi, G. Pilati, G. Bergamini, Dossier Microcogenerazione Progetto RES & RUE Dissemination.
[2] Israel Urieli, David M. Berchowitz, Stirling Cycle Engine Analysis. Adam Hilger Ltd, Redcliffe Way,.Bristol (1984).
[3] N. Benedetti, Tesi di Laurea: Il motore Stirling tra tradizione e innovazione: una scommessa per il futuro' Università degli Studi di Udine (2007-2008).
[4] http://www.triz-journal.com/archives/2009/09/02/.

[5] http://solarcellcentral.com/stirling_page.html.

[6] M. Ciavatta, Tesi di Laurea: Modello Numerico e attività di sperimentazione su un micro-cogeneratore Stirling a gas
naturale. Politecnico di Milano (2011).
[7] EHE, WhisperGen® (EU1) micro CHP Installation & Technical Manual, 2009.
[8] http://glassocean.net/office-2010-breaks-backwards- compatibility-with-certain-vba-projects/.
[9] IAPWS, Revised Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam Switzerland (August 2007). http://www.iapws.org/
[10] IAPWS, Supplementary Release on Backward Equations for Specific Volume as a Function of Pressure and Temperature
v(p,T) for Region 3 of the IAPWS Industrial Formulation 1997
for the Thermodynamic Properties of Water and Steam. Greece (July 2005). http://www.iapws.org/.
[11] http://www.nist.gov.

102
[12] E. Macchi, P. Silva, G. Valenti, N. Fergnani, N. Lazzari, C. Sacco, A. Ravidà, Sviluppo e realizzazione di un micro-
cogeneratore a ciclo Stirling alimentato a gas naturale
(progetto MICROGEN), Politecnico di Milano (Giugno 2012).
[13] E. Zattoni, Tesi di Laurea: Prove Sperimentali su un sistema Micro-Cogenerativo basato su Ciclo Stirling, Politecnico di
Milano (2012).
[14] Stephen R. Turns, An Introduction to Combustion.
Published by McGraw-Hill Higher Education Wiley Publishing Inc., Indianapolis (2004)
[15] P. Kimmel, S. Bullen, J. Green, R. Bovey, R. Rosenberg,
Excel 2003 VBA Programmer''s Reference. Published by Wiley Publishing Inc., Indianapolis (2004)
[16] William R, Martini. Stirling Engine Design Manual.
Martini Engineering Publications (1983) Prepared for National Aeronautics and Space Administration, Lewis Research Center, Under Grant NSG - 3194
[17] G. Valenti, P. Silva, N. Fergnani, G. Di Marcoberardino, S. Campanari, E. Macchi, Experimental and numerical study of a
micro-cogeneration Stirling engine for residential application. In 68° Convegno Nazionale Associazione Termotecnica,
Bologna (2013)
[18] D. Gedeon, J.G. Wood, Oscillating-Flow Regenerator Test Rig: Hardware and Theory With Derived Correlations for Screens
and Felts. NASA Contractor Report 198442. (1996) 103 Allegato 1: Tabella coefficienti Regione 1dell'' acqua.
i Ii Ji ni 1 0 -2 0.146 329 712 131 67 2 0 -1 -0.845 481 871 691 14 3 0 0 -0.375 636 036 720 40 x 101 4 0 1 0.338 551 691 683 85 x 101 5 0 2 -0.957 919 633 878 72 6 0 3 0.157 720 385 132 28 7 0 4 -0.166 164 171 995 01 x 10-1 8 1 5 0.812 146 299 835 68 x 10-3 9 1 -9 0.283 190 801 238 04 x 10-3 10 1 -7 -0.607 063 015 658 74 x 10-3 11 1 -1 -0.189 900 682 184 19 x 10-1 12 1 0 -0.325 297 487 705 05 x 10-1 13 1 1 -0.218 417 171 754 14 x 10-1 14 1 3 -0.528 383 579 699 30 x 10-4 15 2 -3 -0.471 843 210 732 67 x 10-3 16 2 0 -0.300 017 807 930 26 x 10-3 17 2 1 0.476 613 939 069 87 x 10-4 18 2 3 -0.441 418 453 308 46 x 10-5 19 2 17 -0.726 949 962 975 94 x 10-15 20 3 -4 -0.316 796 448 450 54 x 10-4 21 3 0 -0.282 707 979 853 12 x 10-5 22 3 6 -0.852 051 281 201 03 x 10-9 23 4 -5 -0.224 252 819 080 00 x 10-5 24 4 -2 -0.651 712 228 956 01 x 10-6 25 4 10 -0.143 417 299 379 24 x 10-12 26 5 -8 -0.405 169 968 601 17 x 10-6 104
i Ii Ji ni 27 8 -11 -0.127 343 017 416 41 x 10-8 28 8 -6 -0.174 248 712 306 34 x 10-9 29 21 -29 -0.687 621 312 995 31 x 10-18 30 23 -31 0.144 783 078 285 21 x 10-19 31 29 -38 0.263 357 816 627 95 x 10-22 32 30 -39 -0.119 476 226 400 71 x 10-22 33 31 -40 0.182 280 945 814 04 x 10-23 34 32 -41 0.182 280 945 814 04 x 10-23 105 Allegato 2: Tabella coefficienti Regione 2 dell'' acqua
i Ii Ji ni 1 1 0 -0.177 317 424 732 13 x 10-2 2 1 1 -0.178 348 622 923 58 x 10-1 3 1 2 -0.459 960 136 963 65 x 10-1 4 1 3 -0.575 812 590 834 32 x 10-1 5 1 6 -0.503 252 787 279 30 x 10-1 6 2 1 -0.330 326 416 702 03 x 10-4 7 2 2 -0.189 489 875 163 15 x 10-3 8 2 4 -0.393 927 772 433 55 x 10-2 9 2 7 -0.437 972 956 505 73 x 10-1 10 2 36 -0.266 745 479 140 87 x 10-4 11 3 0 0.204 817 376 923 09 x 10-7 12 3 1 0.438 706 672 844 35 x 10-6 13 3 3 -0.322 776 772 385 70 x 10-4 14 3 6 -0.150 339 245 421 48 x 10-2 15 3 35 -0.406 682 535 626 49 x 10-1 16 4 1 -0.788 473 095 593 67 x 10-9 17 4 2 0.127 907 178 522 85 x 10-7 18 4 3 0.482 253 727 185 07 x 10-6 19 5 7 0.229 220 763 376 61 x 10-5 20 6 3 -0.167 147 664 510 61 x 10-10 21 6 16 -0.211 714 723 213 55 x 10-2 22 6 35 -0.238 957 419 341 04 x 102 23 7 0 -0.590 595 643 242 70 x 10-17 24 7 11 -0.126 218 088 991 01 x 10-5 25 7 25 -0.389 468 424 357 39 x 10-1 26 8 8 0.112 562 113 604 59 x 10-10 106
i Ii Ji ni 27 8 36 -0.823 113 408 979 98 x 101 28 9 13 0.198 097 128 020 88 x 10-7 29 10 4 0.104 069 652 101 74 x 10-18 30 10 10 -0.102 347 470 959 29 x 10-12 31 10 14 0.100 181 793 795 11 x 10-8 32 16 29 -0.808 829 086 469 85 x 10-10 33 16 50 0.106 930 318 794 09 34 18 57 -0.336 622 505 741 71 35 20 20 0.891 858 453 554 21 x 10-24 36 20 35 0.306 293 168 762 32 x 10-12 37 20 48 -0.420 024 676 982 08 x 10-5 38 21 21 -0.590 560 296 856 39 x 10-25 39 22 53 0.378 269 476 134 57 x 10-5 40 23 39 -0.127 686 089 346 81 x 10-14 41 24 26 0.730 876 105 950 61 x 10-28 42 24 40 0.554 147 153 507 78 x 10-16 43 24 58 -0.943 697 072 412 10 x 10-6 107 Allegato 3: Tabella coefficienti Regione 4 dell''acqua.
i ni 1 0.116 705 214 527 67 x 104 2 -0.724 213 167 032 06 x 106 3 -0.170 738 469 400 92 x 102 4 0.120 208 247 024 70 x 105 5 -0.323 255 503 223 33 x 107 6 0.149 151 086 135 30 x 102 7 -0.482 326 573 615 91 x 104 8 0.405 113 405 420 57 x 106 9 -0.238 555 575 678 49 10 0.650 175 348 447 98 x 103
108
109 Allegato 4: Tabella coefficienti Regione 5 dell''acqua.
i 1 0 -0.131 799 836 742 01 x 102 2 1 0.685 408 416 344 34 x 101 3 -3 -0.248 051 489 344 66 x 10-1 4 -2 0.369 015 349 803 33 5 -2 -0.311 613 182 139 25 x 101 6 2 -0.329 616 265 389 17

i Ii Ji ni 1 1 1 0.157 364 048 552 x 10-2 2 1 2 0.901 537 616 739 44 x 10-3 3 1 3 -0.502 700 776 776 48 x 10-2 4 2 3 0.224 400 374 094 85 x 10-5 5 2 9 -0.411 632 754 534 71 x 10-5 6 3 7 0.379 194 548 229 55 x 10-7
110
111 Allegato 5: Modulo ''Water'.

' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' DECLARATIONS ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '

' ' add reference !!!

' ' Numerical values of the coefficients and exponents of the dimensionless Gibbs
free energy for region 1, Eq. (7) (Table 2, page 7) '''
Private J_R1(1 To 34) As Integer 'J(i) Private I_R1(1 To 34) As Integer 'I(i) Private N_R1(1 To 34) As Double 'n(i)
' ' Numerical values of the coefficients and exponents of the dimensionless Gibbs
free energy for region 2, Eq.(16) (table 10, page13) and Eq.(17) (Table 11, page
14) '''
Private J0_R2(1 To 9) As Integer 'J0(i) Private N0_R2(1 To 9) As Double 'n0(i) Private J_R2(1 To 43) As Integer 'J(i) Private I_R2(1 To 43) As Integer 'I(i) Private N_R2(1 To 43) As Double 'n(i)
' ' Numerical values of the coefficients of the B23-equation, Eqs. (5) and (6), for
defining the boundary between regions 2 and 3 (Table 1, page 6) '''
Private N_R23(1 To 5) As Double 'n(i)
' ' Numerical values of the coefficients and exponents of the dimensionless
Helmholtz free energy for region 3, Eq.(28) (Table 30, page 30) '''
Private J_R3(1 To 40) As Integer 'J(i) Private I_R3(1 To 40) As Integer 'I(i) Private N_R3(1 To 40) As Double 'n(i)
' ' Numerical values of the coefficients of the dimensionless saturation, EQs (29)
to (30) (TAble 34, page 33) '''
Private N_R4(1 To 10) As Double

' ' Numerical values of the coefficients and exponents of the dimensionless Gibbs
free energy for region 5 '''
Private J0_R5(1 To 6) As Integer 'J0(i) Private N0_R5(1 To 6) As Double 'n0(i) Private J_R5(1 To 6) As Integer 'J(i) Private I_R5(1 To 6) As Integer 'I(i) Private N_R5(1 To 6) As Double 'n(i)
Private Tred As Integer '[K] Reducing Temperature 112
Private Pred As Double '[MPa] Reducing Pressure Private i As Integer 'Index
Private Const MM As Double = 18.0153 '[kg/kmol] Water molar mass Private Const Tc As Double = 647.096 '[K] Water Critical Temperature Private Const Pc As Double = 22.064 '[MPa] Water Critical Pressure Private Const roc As Double = 322 '[kg/m3] Water Critical Density Private Const R As Double = 0.461526 '[kJ/kgK] Specific gas constant
Option Explicit
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' END DECLARATIONS ' '
''' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '

' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' GATEWAY FUNCTION ''
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '

Public Function LMC_water(prop As String, basis As String, Optional T As
Double, Optional P As Double, Optional Reference As String) As Variant
' Gateway function to call water property functions

' Local declaration
Dim X As Double ' numerical results from property functions
Dim msg As String ' string result in case of error message
Dim fOK As Integer ' flag for occuring errors: 0 an error occured, 1 all fine
Dim fBS As Integer ' flag for changeable base: 0 not changeable, 1 changeable from mass to mole
' initialization
X = 0
fOK = 1
fBS = 0
msg = ""

' choose property to be evaluated
If fOK = 1 Then
Select Case LCase(prop)
Case "v" ' m3/kg
X = v_water(T, P)
fBS = 1
Case "cv" 'kJ/kg K
X = cv_water(T, P)
fBS = 1
Case "cp" 'kJ/kg K
X = cp_water(T, P)
fBS = 1 113 Case "u" 'kJ/kg
X = u_water(T, P)
fBS = 1
Case "h" 'kJ/kg
X = h_water(T, P)
If LCase(Reference) = "formation" Then
X = X - 15970.8876744642
End If
fBS = 1
Case "s" 'kJ/kg K
X = s_water(T, P)
If LCase(Reference) = "formation" Then
X = X + 3.51557936460857
End If
fBS = 1
Case "w" 'm/s
X = w_water(T, P)
Case "tsat" ' K
X = Tsat_water(P)
Case "psat" ' MPa
X = Psat_water(T)
Case "mm" ' kg/kmol o g/mol
X = MM
Case "lhv"
X = 0
fBS = 1
Case "hhv"
X = 0
fBS = 1
Case Else
fOK = 0
msg = "Undefined property"
End Select
End If

' change basis
If fOK = 1 Then
Select Case LCase(basis)
Case "mass"
' nothing to do
Case "mole"
' change specific properties to kmol
If fBS = 1 Then
X = X * MM
End If
Case Else
fOK = 0 114
msg = "Undefined basis"
End Select
End If

' fOK check
If fOK = 1 Then
LMC_water = X
Else
LMC_water = msg
End If
End Function
' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' END GATEWAY FUNCTION ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ''

' ' ' ' ' ' ' ' ' ' ' ' '
' ' PROPERTY FUNCTIONS ' '
' ' ' ' ' ' ' ' ' ' ' ' '
' This function is called to calculate water specific volume
Private Function v_water(T As Double, P As Double)
'Range of temperature and pression of Region 1
If T >= 273.15 And T <= 623.15 And P <= 100 And P >= Psat_water(T) Then
Pred = 16.53
Tred = 1386
v_water = R * T * GammaPR1(T, P, Tred, Pred) / Pred / 1000 'Table 3, page 8 'Range of temperature and pression of Region 2
ElseIf (T >= 273.15 And T <= 623.15 And P > 0 And P <= Psat_water(T)) Or (T > 623.15 And T <= 863.15 And P > 0 And P <= P_R23(T)) Or (T > 863.15
And T <= 1073.15 And P > 0 And P <= 100) Then Pred = 1
Tred = 540
v_water = P / Pred / P / 1000 * R * T * (GammaRPR2(T, P, Tred, Pred) + Gamma0PR2(T, P, Tred, Pred)) 'Table 12, page 15 'Range of temperature and pression of Region 3
'ElseIf (T >= 623.15 And T <= T_R23(P) And P >= P_R23(T) And P <= 100) Then 'v_water = 5555
'Range of temperature and pression of Region 5
ElseIf (T >= 1073.15 And T <= 2273.15 And P >= 0 And P <= 50) Then
Pred = 1
Tred = 1000
v_water = P / Pred * R * T / P * (Gamma0PR5(T, P, Tred, Pred) + GammaRPR5(T, P, Tred, Pred)) / 1000 Else
v_water = "Out of range"
End If 115 End Function

' This function is called to calculate water specific isochoric heat capacity
Private Function cv_water(T As Double, P As Double)
'Range of temperature and pression of Region 1
If T <= 623.15 And T >= 273.15 And P <= 100 And P >= Psat_water(T) Then
Pred = 16.53
Tred = 1386
cv_water = R * (-(Tred / T) ^ 2 * GammaTTR1(T, P, Tred, Pred) + ((GammaPR1(T, P, Tred, Pred) - Tred / T * GammaPTR1(T, P, Tred,
Pred)) ^ 2) / GammaPPR1(T, P, Tred, Pred)) 'Table 3, page 8 'Range of temperature and pression of Region 2
ElseIf (T >= 273.15 And T <= 623.15 And P > 0 And P <= Psat_water(T)) Or (T > 623.15 And T <= 863.15 And P > 0 And P <= P_R23(T)) Or (T > 863.15
And T <= 1073.15 And P > 0 And P <= 100) Then Pred = 1
Tred = 540
cv_water = (-(Tred / T) ^ 2 * (Gamma0TTR2(T, P, Tred, Pred) + GammaRTTR2(T, P, Tred, Pred)) - (1 + P / Pred * GammaRPR2(T, P,
Tred, Pred) - Tred / T * P / Pred * GammaRPTR2(T, P, Tred, Pred)) ^ 2 /
(1 - (P / Pred) ^ 2 * GammaRPPR2(T, P, Tred, Pred))) * R 'Table 12, page
15 'Range of temperature and pression of Region 3
'ElseIf (T >= 623.15 And T <= T_R23(P) And P >= P_R23(T) And P <= 100) Then 'cv_water = 5555
'Range of temperature and pression of Region 5
ElseIf (T >= 1073.15 And T <= 2273.15 And P >= 0 And P <= 50) Then
Pred = 1
Tred = 1000
cv_water = (-(Tred / T) ^ 2 * (Gamma0TTR5(T, P, Tred, Pred) + GammaRTTR5(T, P, Tred, Pred)) - (1 + P / Pred * GammaRPR5(T, P,
Tred, Pred) - Tred / T * P / Pred * GammaRPTR5(T, P, Tred, Pred)) ^ 2 /
(1 - (P / Pred) ^ 2 * GammaRPPR5(T, P, Tred, Pred))) * R Else
cv_water = "Out of range"
End If
End Function

' This function is called to calculate water specific isobaric heat capacity
Private Function cp_water(T As Double, P As Double)
'Range of temperature and pression of Region 1
If T <= 623.15 And T >= 273.15 And P <= 100 And P >= Psat_water(T) Then
Pred = 16.53
Tred = 1386
cp_water = R * -(Tred / T) ^ 2 * GammaTTR1(T, P, Tred, Pred) 'Table 3, page 8 116
'Range of temperature and pression of Region 2
ElseIf (T >= 273.15 And T <= 623.15 And P > 0 And P <= Psat_water(T)) Or (T > 623.15 And T <= 863.15 And P > 0 And P <= P_R23(T)) Or (T > 863.15
And T <= 1073.15 And P > 0 And P <= 100) Then Pred = 1
Tred = 540
cp_water = (-(Tred / T) ^ 2 * (Gamma0TTR2(T, P, Tred, Pred) + GammaRTTR2(T, P, Tred, Pred))) * R 'Table 12, page 15 'Range of temperature and pression of Region 3
'ElseIf (T >= 623.15 And T <= T_R23(P) And P >= P_R23(T) And P <= 100) Then 'cp_water = 5555
'Range of temperature and pression of Region 5
ElseIf (T >= 1073.15 And T <= 2273.15 And P >= 0 And P <= 50) Then
Pred = 1
Tred = 1000
cp_water = (-(Tred / T) ^ 2 * (Gamma0TTR5(T, P, Tred, Pred) + GammaRTTR5(T, P, Tred, Pred))) * R Else
cp_water = "Out of range"
End If
End Function

' This function is called to calculate water specific internal energy
Private Function u_water(T As Double, P As Double)
'Range of temperature and pression of Region 1
If T >= 273.15 And T <= 623.15 And P <= 100 And P >= Psat_water(T) Then
Pred = 16.53
Tred = 1386
u_water = R * T * ((Tred / T) * GammaTR1(T, P, Tred, Pred) - (P / Pred) * GammaPR1(T, P, Tred, Pred)) 'Table 3, page 8 'Range of temperature and pression of Region 2
ElseIf (T >= 273.15 And T <= 623.15 And P > 0 And P <= Psat_water(T)) Or (T > 623.15 And T <= 863.15 And P > 0 And P <= P_R23(T)) Or (T > 863.15
And T <= 1073.15 And P > 0 And P <= 100) Then Pred = 1
Tred = 540
u_water = (Tred / T * (Gamma0TR2(T, P, Tred, Pred) + GammaRTR2(T, P, Tred, Pred)) - P / Pred * (Gamma0PR2(T, P, Tred, Pred) +
GammaRPR2(T, P, Tred, Pred))) * R * T 'Table 12, page 15 'Range of temperature and pression of Region 3
'ElseIf (T >= 623.15 And T <= T_R23(P) And P >= P_R23(T) And P <= 100) Then 'u_water = 5555
'Range of temperature and pression of Region 5
ElseIf (T >= 1073.15 And T <= 2273.15 And P >= 0 And P <= 50) Then
Pred = 1 117 Tred = 1000
u_water = (Tred / T * (Gamma0TR5(T, P, Tred, Pred) + GammaRTR5(T, P, Tred, Pred)) - P / Pred * (Gamma0PR5(T, P, Tred, Pred) +
GammaRPR5(T, P, Tred, Pred))) * R * T Else
u_water = "Out of range"
End If
End Function

' This function is called to calculate water specific enthalpy
Private Function h_water(T As Double, P As Double)
'Range of temperature and pression of Region 1
If T <= 623.15 And T >= 273.15 And P <= 100 And P >= Psat_water(T) Then
Pred = 16.53
Tred = 1386
h_water = Tred * R * GammaTR1(T, P, Tred, Pred) 'Table 3, page 8
'Range of temperature and pression of Region 2
ElseIf (T >= 273.15 And T <= 623.15 And P > 0 And P <= Psat_water(T)) Or (T > 623.15 And T <= 863.15 And P > 0 And P <= P_R23(T)) Or (T > 863.15
And T <= 1073.15 And P > 0 And P <= 100) Then Pred = 1
Tred = 540
h_water = Tred / T * (Gamma0TR2(T, P, Tred, Pred) + GammaRTR2(T, P, Tred, Pred)) * R * T 'Table 12, page 15 'Range of temperature and pression of Region 3
'ElseIf (T >= 623.15 And T <= T_R23(P) And P >= P_R23(T) And P <= 100) Then 'h_water = 5555
'Range of temperature and pression of Region 5
ElseIf (T >= 1073.15 And T <= 2273.15 And P >= 0 And P <= 50) Then
Pred = 1
Tred = 1000
h_water = R * T * Tred / T * (Gamma0TR5(T, P, Tred, Pred) + GammaRTR5(T, P, Tred, Pred)) Else
h_water = "Out of range"
End If
End Function

This function is called to calculate water specific entropy
Private Function s_water(T As Double, P As Double)
'Range of temperature and pression of Region 1
If T >= 273.15 And T <= 623.15 And P <= 100 And P >= Psat_water(T) Then
Pred = 16.53
Tred = 1386
s_water = R * ((Tred / T) * GammaTR1(T, P, Tred, Pred) - GammaR1(T, P, Tred, Pred)) 'Table 3, page 8 118
'Range of temperature and pression of Region 2
ElseIf (T >= 273.15 And T <= 623.15 And P > 0 And P <= Psat_water(T)) Or (T > 623.15 And T <= 863.15 And P > 0 And P <= P_R23(T)) Or (T > 863.15
And T <= 1073.15 And P > 0 And P <= 100) Then Pred = 1
Tred = 540
s_water = (Tred / T * (Gamma0TR2(T, P, Tred, Pred) + GammaRTR2(T, P, Tred, Pred)) - (Gamma0R2(T, P, Tred, Pred) + GammaRR2(T, P, Tred,
Pred))) * R 'Table 12, page 15 'Range of temperature and pression of Region 3
'ElseIf (T >= 623.15 And T <= T_R23(P) And P >= P_R23(T) And P <= 100) Then 's_water = 5555
'Range of temperature and pression of Region 5
ElseIf (T >= 1073.15 And T <= 2273.15 And P >= 0 And P <= 50) Then
Pred = 1
Tred = 1000
s_water = R * (Tred / T * (Gamma0TR5(T, P, Tred, Pred) + GammaRTR5(T, P, Tred, Pred)) - (Gamma0R5(T, P, Tred, Pred) +
GammaRR5(T, P, Tred, Pred))) Else
s_water = "Out of range"
End If
End Function

' This function is called to calculate water speed of sound
Private Function w_water(T As Double, P As Double) 'Function miltiplied by
sqrt(1000) to obtain the correct value
'Range of temperature and pression of Region 1
If T <= 623.15 And T >= 273.15 And P <= 100 And P >= Psat_water(T) Then
Pred = 16.53
Tred = 1386
w_water = (Sqr(R * T * GammaPR1(T, P, Tred, Pred) ^ 2 / ((GammaPR1(T, P, Tred, Pred) - Tred / T * GammaPTR1(T, P, Tred, Pred)) ^ 2 / ((Tred / T)
^ 2 * GammaTTR1(T, P, Tred, Pred)) - GammaPPR1(T, P, Tred, Pred)))) *
Sqr(1000) 'Table 3, page 8 'Range of temperature and pression of Region 2
ElseIf (T >= 273.15 And T <= 623.15 And P > 0 And P <= Psat_water(T)) Or (T > 623.15 And T <= 863.15 And P > 0 And P <= P_R23(T)) Or (T > 863.15
And T <= 1073.15 And P > 0 And P <= 100) Then Pred = 1
Tred = 540
w_water = (Sqr(R * T * ((1 + 2 * P / Pred * GammaRPR2(T, P, Tred, Pred) + (P / Pred) ^ 2 * GammaRPR2(T, P, Tred, Pred) ^ 2) / (1 - (P / Pred) ^ 2 *
GammaRPR2(T, P, Tred, Pred) + (1 + P / Pred * GammaRPR2(T, P, Tred,
Pred) - Tred / T * P / Pred * GammaRPTR2(T, P, Tred, Pred)) ^ 2 / ((Tred / 119 T) ^ 2 * (Gamma0TTR2(T, P, Tred, Pred) + GammaRTTR2(T, P, Tred,
Pred))))))) * Sqr(1000) 'Table 12, page 15 'Range of temperature and pression of Region 3
'ElseIf (T >= 623.15 And T <= T_R23(P) And P >= P_R23(T) And P <= 100) Then 'w_water = 5555
'Range of temperature and pression of Region 5
ElseIf (T >= 1073.15 And T <= 2273.15 And P >= 0 And P <= 50) Then
Pred = 1
Tred = 1000
w_water = (Sqr(R * T * ((1 + 2 * P / Pred * GammaRPR5(T, P, Tred, Pred) + (P / Pred) ^ 2 * GammaRPR5(T, P, Tred, Pred) ^ 2) / (1 - (P / Pred) ^ 2 *
GammaRPR5(T, P, Tred, Pred) + (1 + P / Pred * GammaRPR5(T, P, Tred,
Pred) - Tred / T * P / Pred * GammaRPTR5(T, P, Tred, Pred)) ^ 2 / ((Tred /
T) ^ 2 * (Gamma0TTR5(T, P, Tred, Pred) + GammaRTTR5(T, P, Tred,
Pred))))))) * Sqr(1000) Else
w_water = "Out of range"
End If
End Function

' This function is called to calculate saturation temperature
' Solution of EQ.(29) with regard to saturation temperature
Private Function Tsat_water(P As Double)

Dim beta As Double, d As Double, e As Double, f As Double, g As Double

'Range of pressure of Region 4
If P >= 0.000661213 And P <= 22.064 Then
Call CoefR4
Tred = 1
Pred = 1
beta = (P / Pred) ^ 0.25 'EQ.(29a)
e = beta ^ 2 + N_R4(3) * beta + N_R4(6)
f = N_R4(1) * beta ^ 2 + N_R4(4) * beta + N_R4(7)
g = N_R4(2) * beta ^ 2 + N_R4(5) * beta + N_R4(8)
d = 2 * g / (-f - (f ^ 2 - 4 * e * g) ^ 0.5)
Tsat_water = Tred * (N_R4(10) + d - ((N_R4(10) + d) ^ 2 - 4 * (N_R4(9) + N_R4(10) * d)) ^ 0.5) / 2 Else
Tsat_water = "Out of range"
End If
End Function

' This function is called to calculate saturation pressure
' Solution of EQ.(29) with regard to saturation pressure
Private Function Psat_water(T As Double) 120

Dim teta As Double, A As Double, b As Double, c As Double
'Range of temperature of Region 4
If T >= 273.15 And T <= 647.096 Then
Call CoefR4
Tred = 1
Pred = 1
teta = T / Tred + N_R4(9) / (T / Tred - N_R4(10)) 'EQ.(29b)
A = teta ^ 2 + N_R4(1) * teta + N_R4(2)
b = N_R4(3) * teta ^ 2 + N_R4(4) * teta + N_R4(5)
c = N_R4(6) * teta ^ 2 + N_R4(7) * teta + N_R4(8)
Psat_water = Pred * (2 * c / (-b + Sqr(b ^ 2 - 4 * A * c))) ^ 4
Else
Psat_water = "Out of range"
End If
End Function
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' END PROPERTY FUNCTIONS ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' REGION 1 ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '

' ' The dimensionless Gibbs free energy GAMMA and its derivatives according to
Eq.(7) (Table 4, page 8)

Private Function GammaR1(T As Double, P As Double, Tred As Integer, Pred As
Double)
Call CoefR1
GammaR1 = 0
For i = 1 To 34 GammaR1 = GammaR1 + N_R1(i) * (7.1 - P / Pred) ^ I_R1(i) * (Tred / T -
1.222) ^ J_R1(i) Next i
End Function

Private Function GammaTR1(T As Double, P As Double, Tred As Integer, Pred
As Double)
Call CoefR1
GammaTR1 = 0
For i = 1 To 34 GammaTR1 = GammaTR1 + N_R1(i) * (7.1 - P / Pred) ^ I_R1(i) * J_R1(i)
* (Tred / T - 1.222) ^ (J_R1(i) - 1) Next i
End Function
121 Private Function GammaPR1(T As Double, P As Double, Tred As Integer, Pred
As Double)
Call CoefR1
GammaPR1 = 0
For i = 1 To 34 GammaPR1 = GammaPR1 - N_R1(i) * I_R1(i) * (7.1 - P / Pred) ^ (I_R1(i)
- 1) * (Tred / T - 1.222) ^ J_R1(i) Next i
End Function

Private Function GammaTTR1(T As Double, P As Double, Tred As Integer, Pred
As Double)
Call CoefR1
GammaTTR1 = 0
For i = 1 To 34 GammaTTR1 = GammaTTR1 + N_R1(i) * (7.1 - P / Pred) ^ I_R1(i) *
J_R1(i) * (J_R1(i) - 1) * (Tred / T - 1.222) ^ (J_R1(i) - 2) Next i
End Function

Private Function GammaPPR1(T As Double, P As Double, Tred As Integer, Pred
As Double)
Call CoefR1
GammaPPR1 = 0
For i = 1 To 34 GammaPPR1 = GammaPPR1 + N_R1(i) * I_R1(i) * (I_R1(i) - 1) * (7.1 - P
/ Pred) ^ (I_R1(i) - 2) * (Tred / T - 1.222) ^ J_R1(i) Next i
End Function

Private Function GammaPTR1(T As Double, P As Double, Tred As Integer, Pred
As Double)
Call CoefR1
GammaPTR1 = 0
For i = 1 To 34 GammaPTR1 = GammaPTR1 - N_R1(i) * I_R1(i) * (7.1 - P / Pred) ^
(I_R1(i) - 1) * J_R1(i) * (Tred / T - 1.222) ^ (J_R1(i) - 1) Next i
End Function

' ' Numerical values of the coefficients and exponents of the dimensionless Gibbs
free energy for region 1, Eq. (7) (Table 2, page 7)
Private Sub CoefR1()
For i = 1 To 8
I_R1(i) = 0
Next i
For i = 9 To 14 122
I_R1(i) = 1
Next i
For i = 15 To 19
I_R1(i) = 2
Next i
For i = 20 To 22
I_R1(i) = 3
Next i
For i = 23 To 25
I_R1(i) = 4
Next i
I_R1(26) = 5
I_R1(27) = 8
I_R1(28) = 8
I_R1(29) = 21
I_R1(30) = 23
I_R1(31) = 29
I_R1(32) = 30
I_R1(33) = 31
I_R1(34) = 32

For i = 1 To 8
J_R1(i) = i - 3
Next i
J_R1(9) = -9
J_R1(10) = -7
J_R1(11) = -1
J_R1(12) = 0
J_R1(13) = 1
J_R1(14) = 3
J_R1(15) = -3
J_R1(16) = 0
J_R1(17) = 1
J_R1(18) = 3
J_R1(19) = 17
J_R1(20) = -4
J_R1(21) = 0
J_R1(22) = 6
J_R1(23) = -5
J_R1(24) = -2
J_R1(25) = 10
J_R1(26) = -8
J_R1(27) = -11
J_R1(28) = -6
J_R1(29) = -29
J_R1(30) = -31
J_R1(31) = -38 123 J_R1(32) = -39
J_R1(33) = -40
J_R1(34) = -41

N_R1(1) = 0.14632971213167
N_R1(2) = -0.84548187169114
N_R1(3) = -3.756360367204
N_R1(4) = 3.3855169168385
N_R1(5) = -0.95791963387872
N_R1(6) = 0.15772038513228
N_R1(7) = -0.016616417199501
N_R1(8) = 8.1214629983568E-04
N_R1(9) = 2.8319080123804E-04
N_R1(10) = -6.0706301565874E-04
N_R1(11) = -0.018990068218419
N_R1(12) = -0.032529748770505
N_R1(13) = -0.021841717175414
N_R1(14) = -5.283835796993E-05
N_R1(15) = -4.7184321073267E-04
N_R1(16) = -3.0001780793026E-04
N_R1(17) = 4.7661393906987E-05
N_R1(18) = -4.4141845330846E-06
N_R1(19) = -7.2694996297594E-16
N_R1(20) = -3.1679644845054E-05
N_R1(21) = -2.8270797985312E-06
N_R1(22) = -8.5205128120103E-10
N_R1(23) = -2.2425281908E-06
N_R1(24) = -6.5171222895601E-07
N_R1(25) = -1.4341729937924E-13
N_R1(26) = -4.0516996860117E-07
N_R1(27) = -1.2734301741641E-09
N_R1(28) = -1.7424871230634E-10
N_R1(29) = -6.8762131295531E-19
N_R1(30) = 1.4478307828521E-20
N_R1(31) = 2.6335781662795E-23
N_R1(32) = -1.1947622640071E-23
N_R1(33) = 1.8228094581404E-24
N_R1(34) = -9.3537087292458E-26
End Sub
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' END REGION 1 ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' REGION 2 ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
124
' ' The ideal-gas part GAMMA0 of the dimensionless Gibbs free energy and its
derivatives according to Eq.(16) (Table 13, page 16)

Private Function Gamma0R2(T As Double, P As Double, Tred As Integer, Pred
As Double)
Call CoefR2
Gamma0R2 = Log(P / Pred)
For i = 1 To 9
Gamma0R2 = Gamma0R2 + N0_R2(i) * (Tred / T) ^ J0_R2(i)
Next i
End Function

Private Function Gamma0PR2(T As Double, P As Double, Tred As Integer, Pred
As Double)
Gamma0PR2 = 1 / P * Pred
End Function

Private Function Gamma0PPR2(T As Double, P As Double, Tred As Integer,
Pred As Double)
Gamma0PPR2 = -1 / (P / Pred) ^ 2
End Function

Private Function Gamma0TR2(T As Double, P As Double, Tred As Integer, Pred
As Double)
Call CoefR2
Gamma0TR2 = 0
For i = 1 To 9 Gamma0TR2 = Gamma0TR2 + N0_R2(i) * J0_R2(i) * (Tred / T) ^
(J0_R2(i) - 1) Next i
End Function

Private Function Gamma0TTR2(T As Double, P As Double, Tred As Integer,
Pred As Double)
Call CoefR2
Gamma0TTR2 = 0
For i = 1 To 9 Gamma0TTR2 = Gamma0TTR2 + N0_R2(i) * J0_R2(i) * (J0_R2(i) - 1) *
(Tred / T) ^ (J0_R2(i) - 2) Next i
End Function

Private Function Gamma0PTR2(T As Double, P As Double, Tred As Integer,
Pred As Double)
Gamma0PTR2 = 0
End Function
125 ' ' The residual part GAMMAR of the dimensionless Gibbs free energy and its
derivatives according to Eq.(17) (Table 14, page 16)

Private Function GammaRR2(T As Double, P As Double, Tred As Integer, Pred
As Double)
Call CoefR2
GammaRR2 = 0
For i = 1 To 43 GammaRR2 = GammaRR2 + N_R2(i) * (P / Pred) ^ I_R2(i) * (Tred / T -
0.5) ^ J_R2(i) Next i
End Function

Private Function GammaRPR2(T As Double, P As Double, Tred As Integer, Pred
As Double)
Call CoefR2
GammaRPR2 = 0
For i = 1 To 43 GammaRPR2 = GammaRPR2 + N_R2(i) * I_R2(i) * (P / Pred) ^ (I_R2(i) -
1) * (Tred / T - 0.5) ^ J_R2(i) Next i
End Function

Private Function GammaRTR2(T As Double, P As Double, Tred As Integer, Pred
As Double)
Call CoefR2
GammaRTR2 = 0
For i = 1 To 43 GammaRTR2 = GammaRTR2 + N_R2(i) * (P / Pred) ^ I_R2(i) * J_R2(i) *
(Tred / T - 0.5) ^ (J_R2(i) - 1) Next i
End Function

Private Function GammaRPPR2(T As Double, P As Double, Tred As Integer,
Pred As Double)
Call CoefR2
GammaRPPR2 = 0
For i = 1 To 43 GammaRPPR2 = GammaRPPR2 + N_R2(i) * I_R2(i) * (I_R2(i) - 1) * (P /
Pred) ^ (I_R2(i) - 2) * (Tred / T - 0.5) ^ J_R2(i) Next i
End Function
Private Function GammaRTTR2(T As Double, P As Double, Tred As Integer,
Pred As Double)
Call CoefR2
GammaRTTR2 = 0
For i = 1 To 43 126
GammaRTTR2 = GammaRTTR2 + N_R2(i) * (P / Pred) ^ I_R2(i) *
J_R2(i) * (J_R2(i) - 1) * (Tred / T - 0.5) ^ (J_R2(i) - 2) Next i
End Function

Private Function GammaRPTR2(T As Double, P As Double, Tred As Integer,
Pred As Double)
Call CoefR2
GammaRPTR2 = 0
For i = 1 To 43 GammaRPTR2 = GammaRPTR2 + N_R2(i) * I_R2(i) * (P / Pred) ^
(I_R2(i) - 1) * J_R2(i) * (Tred / T - 0.5) ^ (J_R2(i) - 1) Next i
End Function

' ' Numerical values of the coefficients and exponents of the dimensionless Gibbs
free energy for region 2

Private Sub CoefR2()
'ideal-gas part GAMMA0, Eq. (16)
J0_R2(1) = 0
J0_R2(2) = 1
J0_R2(3) = -5
J0_R2(4) = -4
J0_R2(5) = -3
J0_R2(6) = -2
J0_R2(7) = -1
J0_R2(8) = 2
J0_R2(9) = 3

N0_R2(1) = -9.6927686500217
N0_R2(2) = 10.086655968018
N0_R2(3) = -0.005608791128302
N0_R2(4) = 0.071452738081455
N0_R2(5) = -0.40710498223928
N0_R2(6) = 1.4240819171444
N0_R2(7) = -4.383951131945
N0_R2(8) = -0.28408632460772
N0_R2(9) = 0.021268463753307

'residual part GAMMAR, Eq. (17)
For i = 1 To 5
I_R2(i) = 1
Next i
For i = 6 To 10
I_R2(i) = 2
Next i 127 For i = 11 To 15
I_R2(i) = 3
Next i
For i = 16 To 18
I_R2(i) = 4
Next i
I_R2(19) = 5
For i = 20 To 22
I_R2(i) = 6
Next i
For i = 23 To 25
I_R2(i) = 7
Next i
I_R2(26) = 8
I_R2(27) = 8
I_R2(28) = 9
For i = 29 To 31
I_R2(i) = 10
Next i
I_R2(32) = 16
I_R2(33) = 16
I_R2(34) = 18
For i = 35 To 37
I_R2(i) = 20
Next i
I_R2(38) = 21
I_R2(39) = 22
I_R2(40) = 23
For i = 41 To 43
I_R2(i) = 24
Next i

J_R2(1) = 0
J_R2(2) = 1
J_R2(3) = 2
J_R2(4) = 3
J_R2(5) = 6
J_R2(6) = 1
J_R2(7) = 2
J_R2(8) = 4
J_R2(9) = 7
J_R2(10) = 36
J_R2(11) = 0
J_R2(12) = 1
J_R2(13) = 3
J_R2(14) = 6
J_R2(15) = 35 128
J_R2(16) = 1
J_R2(17) = 2
J_R2(18) = 3
J_R2(19) = 7
J_R2(20) = 3
J_R2(21) = 16
J_R2(22) = 35
J_R2(23) = 0
J_R2(24) = 11
J_R2(25) = 25
J_R2(26) = 8
J_R2(27) = 36
J_R2(28) = 13
J_R2(29) = 4
J_R2(30) = 10
J_R2(31) = 14
J_R2(32) = 29
J_R2(33) = 50
J_R2(34) = 57
J_R2(35) = 20
J_R2(36) = 35
J_R2(37) = 48
J_R2(38) = 21
J_R2(39) = 53
J_R2(40) = 39
J_R2(41) = 26
J_R2(42) = 40
J_R2(43) = 58

N_R2(1) = -1.7731742473213E-03
N_R2(2) = -0.017834862292358
N_R2(3) = -0.045996013696365
N_R2(4) = -0.057581259083432
N_R2(5) = -0.05032527872793
N_R2(6) = -3.3032641670203E-05
N_R2(7) = -1.8948987516315E-04
N_R2(8) = -3.9392777243355E-03
N_R2(9) = -0.043797295650573
N_R2(10) = -2.6674547914087E-05
N_R2(11) = 2.0481737692309E-08
N_R2(12) = 4.3870667284435E-07
N_R2(13) = -3.227767723857E-05
N_R2(14) = -1.5033924542148E-03
N_R2(15) = -0.040668253562649
N_R2(16) = -7.8847309559367E-10
N_R2(17) = 1.2790717852285E-08
N_R2(18) = 4.8225372718507E-07 129 N_R2(19) = 2.2922076337661E-06
N_R2(20) = -1.6714766451061E-11
N_R2(21) = -2.1171472321355E-03
N_R2(22) = -23.895741934104
N_R2(23) = -5.905956432427E-18
N_R2(24) = -1.2621808899101E-06
N_R2(25) = -0.038946842435739
N_R2(26) = 1.1256211360459E-11
N_R2(27) = -8.2311340897998
N_R2(28) = 1.9809712802088E-08
N_R2(29) = 1.0406965210174E-19
N_R2(30) = -1.0234747095929E-13
N_R2(31) = -1.0018179379511E-09
N_R2(32) = -8.0882908646985E-11
N_R2(33) = 0.10693031879409
N_R2(34) = -0.33662250574171
N_R2(35) = 8.9185845355421E-25
N_R2(36) = 3.0629316876232E-13
N_R2(37) = -4.2002467698208E-06
N_R2(38) = -5.9056029685639E-26
N_R2(39) = 3.7826947613457E-06
N_R2(40) = -1.2768608934681E-15
N_R2(41) = 7.3087610595061E-29
N_R2(42) = 5.5414715350778E-17
N_R2(43) = -9.436970724121E-07
End Sub
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' END REGION 2 ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '

' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' BOUNDARY BETWEEN REGIONS 2 AND 3 ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' Auxiliary equation for the boundary between Regions 2 and 3 (page 5)
Function P_R23(T As Double)
Pred = 1
Tred = 1
Call CoefR23 P_R23 = Pred * (N_R23(1) + N_R23(2) * T / Tred + N_R23(3) * (T /
Tred) ^ 2) 'Eq.(5), page 5 End Function

Function T_R23(P As Double)
Pred = 1
Tred = 1
Call CoefR23
If P < 16.5291 Or P > 100 Then 130
T_R23 = 0
Else
T_R23 = (N_R23(4) + ((P / Pred - N_R23(5)) / N_R23(3)) ^ 0.5) * Tred
'Eq.(6), page 6
End If
End Function

' ' Numerical values of the coefficients of the B23-equation, Eqs. (5) and (6), for
defining the boundary between regions 2 and 3 (Table 1, page 6)

Private Sub CoefR23()
N_R23(1) = 348.05185628969
N_R23(2) = -1.1671859879975
N_R23(3) = 1.0192970039326E-03
N_R23(4) = 572.54459862746
N_R23(5) = 13.91883977887
End Sub
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '' ' '
' ' END BOUNDARY BETWEEN REGIONS 2 AND 3 ''
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '

' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' REGION 3 ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '

' 'The dimensionless Helmholtz free energy equation and its derivatives according
to Eq.(28) for Region 3 (Table 32, page 31)

Private Function FI(T As Double, ro As Double)
Call CoefR3
FI = N_R3(1) * Log(ro / roc)
For i = 1 To 40
FI = FI + N_R3(i) * (ro / roc) ^ I_R3(i) * (Tc / T) ^ J_R3(i)
Next
End Function

Private Function FIro(T As Double, ro As Double)
Call CoefR3
FIro = N_R3(1) / (ro / roc)
For i = 1 To 40 FIro = FIro + N_R3(i) * I_R3(i) * (ro / roc) ^ (I_R3(i) - 1) * (Tc / T) ^
J_R3(i) Next i
End Function

Private Function FIroro(T As Double, ro As Double)
Call CoefR3 131 FIroro = -N_R3(1) / (ro / roc) ^ 2
For i = 1 To 40 FIroro = FIroro + N_R3(i) * I_R3(i) * (I_R3(i) - 1) * (ro / roc) ^ (I_R3(i) -
2) * (Tc / T) ^ J_R3(i) Next i
End Function

Private Function FIT(T As Double, ro As Double)
Call CoefR3
FIT = 0
For i = 1 To 40
FIT = FIT + N_R3(i) * (ro / roc) ^ I_R3(i) * J_R3(i) * (Tc / T) ^ (J_R3(i) -
1)
Next i
End Function

Private Function FITT(T As Double, ro As Double)
Call CoefR3
FITT = 0
For i = 1 To 40 FITT = FITT + N_R3(i) * (ro / roc) ^ I_R3(i) * J_R3(i) * (J_R3(i) - 1) *
(Tc / T) ^ (J_R3(i) - 2) Next i
End Function

Private Function FIroT(T As Double, ro As Double)
Call CoefR3
FIroT = 0
For i = 1 To 40 FIroT = FIroT + N_R3(i) * I_R3(i) * (ro / roc) ^ (I_R3(i) - 1) * J_R3(i) *
(Tc / T) ^ (J_R3(i) - 1) Next i
End Function

' ' Numerical values of the coefficients and exponents of the dimensionless
Helmholtz free energy for region 3, Eq.(28) (Table 30, page 30)

Private Sub CoefR3()
For i = 1 To 8
I_R3(i) = 0
Next i
For i = 9 To 12
I_R3(i) = 1
Next i
For i = 13 To 18
I_R3(i) = 2
Next i 132
For i = 19 To 23
I_R3(i) = 3
Next i
For i = 24 To 27
I_R3(i) = 4
Next i
For i = 28 To 30
I_R3(i) = 5
Next i
For i = 31 To 33
I_R3(i) = 6
Next i
I_R3(34) = 7
I_R3(35) = 8
I_R3(36) = 9
I_R3(37) = 9
I_R3(38) = 10
I_R3(39) = 10
I_R3(40) = 11
J_R3(1) = 0
J_R3(2) = 0
J_R3(3) = 1
J_R3(4) = 2
J_R3(5) = 7
J_R3(6) = 10
J_R3(7) = 12
J_R3(8) = 23
J_R3(9) = 2
J_R3(10) = 6
J_R3(11) = 15
J_R3(12) = 17
J_R3(13) = 0
J_R3(14) = 2
J_R3(15) = 6
J_R3(16) = 7
J_R3(17) = 22
J_R3(18) = 26
J_R3(19) = 0
J_R3(20) = 2
J_R3(21) = 4
J_R3(22) = 16
J_R3(23) = 26
J_R3(24) = 0
J_R3(25) = 2
J_R3(26) = 4
J_R3(27) = 26
J_R3(28) = 1 133 J_R3(29) = 3
J_R3(30) = 26
J_R3(31) = 0
J_R3(32) = 2
J_R3(33) = 26
J_R3(34) = 2
J_R3(35) = 26
J_R3(36) = 2
J_R3(37) = 26
J_R3(38) = 0
J_R3(39) = 1
J_R3(40) = 26

N_R3(1) = 1.0658070028513
N_R3(2) = -15.732845290239
N_R3(3) = 20.944396974307
N_R3(4) = -7.6867707878716
N_R3(5) = 2.6185947787954
N_R3(6) = -2.808078114862
N_R3(7) = 1.2053369696517
N_R3(8) = -8.4566812812502E-03
N_R3(9) = -1.2654315477714
N_R3(10) = -1.1524407806681
N_R3(11) = 0.88521043984318
N_R3(12) = -0.64207765181607
N_R3(13) = 0.38493460186671
N_R3(14) = -0.85214708824206
N_R3(15) = 4.8972281541877
N_R3(16) = -3.0502617256965
N_R3(17) = 0.039420536879154
N_R3(18) = 0.12558408424308
N_R3(19) = -0.2799932969871
N_R3(20) = 1.389979956946
N_R3(21) = -2.018991502357
N_R3(22) = -8.2147637173963E-03
N_R3(23) = -0.47596035734923
N_R3(24) = 0.0439840744735
N_R3(25) = -0.44476435428739
N_R3(26) = 0.90572070719733
N_R3(27) = 0.70522450087967
N_R3(28) = 0.10770512626332
N_R3(29) = -0.32913623258954
N_R3(30) = -0.50871062041158
N_R3(31) = -0.022175400873096
N_R3(32) = 0.094260751665092
N_R3(33) = 0.16436278447961
N_R3(34) = -0.013503372241348 134
N_R3(35) = -0.014834345352472
N_R3(36) = 5.7922953628084E-04
N_R3(37) = 3.2308904703711E-03
N_R3(38) = 8.0964802996215E-05
N_R3(39) = -1.6557679795037E-04
N_R3(40) = -4.4923899061815E-05
End Sub
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' END REGION 3 ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '

' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' REGION 4 ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '

' ' Numerical values of the coefficients of the dimensionless saturation, EQs (29)
to (30) (TAble 34, page 33)
Private Sub CoefR4()
N_R4(1) = 1167.0521452767
N_R4(2) = -724213.16703206
N_R4(3) = -17.073846940092
N_R4(4) = 12020.82470247
N_R4(5) = -3232555.0322333
N_R4(6) = 14.91510861353
N_R4(7) = -4823.2657361591
N_R4(8) = 405113.40542057
N_R4(9) = -0.23855557567849
N_R4(10) = 650.17534844798
End Sub
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' END REGION 4 ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '

' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' REGION 5 ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' The ideal-gas part GAMMA0 of the dimensionless Gibbs free energy and its
derivatives according to Eq.(33) (Table 40, page 38)
Private Function Gamma0R5(T As Double, P As Double, Tred As Integer, Pred
As Double)
Call CoefR5
Gamma0R5 = Log(P / Pred)
For i = 1 To 6
Gamma0R5 = Gamma0R5 + N0_R5(i) * (Tred / T) ^ J0_R5(i)
Next i
End Function
135 Private Function Gamma0PR5(T As Double, P As Double, Tred As Integer, Pred
As Double)
Gamma0PR5 = 1 / P * Pred
End Function

Private Function Gamma0PPR5(T As Double, P As Double, Tred As Integer,
Pred As Double)
Gamma0PPR5 = -1 / (P / Pred) ^ 2
End Function

Private Function Gamma0TR5(T As Double, P As Double, Tred As Integer, Pred
As Double)
Call CoefR5
Gamma0TR5 = 0
For i = 1 To 6 Gamma0TR5 = Gamma0TR5 + N0_R5(i) * J0_R5(i) * (Tred / T) ^
(J0_R5(i) - 1) Next i
End Function

Private Function Gamma0TTR5(T As Double, P As Double, Tred As Integer,
Pred As Double)
Call CoefR5
Gamma0TTR5 = 0
For i = 1 To 6 Gamma0TTR5 = Gamma0TTR5 + N0_R5(i) * J0_R5(i) * (J0_R5(i) - 1) *
(Tred / T) ^ (J0_R5(i) - 2) Next i
End Function

Private Function Gamma0PTR5(T As Double, P As Double, Tred As Integer,
Pred As Double)
Gamma0PTR5 = 0
End Function

' ' The residual part GAMMAR of the dimensionless Gibbs free energy and its
derivatives according to Eq.(34) (Table 41, page 38)
Private Function GammaRR5(T As Double, P As Double, Tred As Integer, Pred
As Double)
Call CoefR5
GammaRR5 = 0
For i = 1 To 6 GammaRR5 = GammaRR5 + N_R5(i) * (P / Pred) ^ I_R5(i) * (Tred / T) ^
J_R5(i) Next i
End Function
136
Private Function GammaRPR5(T As Double, P As Double, Tred As Integer, Pred
As Double)
Call CoefR5
GammaRPR5 = 0
For i = 1 To 6 GammaRPR5 = GammaRPR5 + N_R5(i) * I_R5(i) * (P / Pred) ^ (I_R5(i) -
1) * (Tred / T) ^ J_R5(i) Next i
End Function

Private Function GammaRTR5(T As Double, P As Double, Tred As Integer, Pred
As Double)
Call CoefR5
GammaRTR5 = 0
For i = 1 To 6 GammaRTR5 = GammaRTR5 + N_R5(i) * (P / Pred) ^ I_R5(i) * J_R5(i) *
(Tred / T) ^ (J_R5(i) - 1) Next i
End Function

Private Function GammaRPPR5(T As Double, P As Double, Tred As Integer,
Pred As Double)
Call CoefR5
GammaRPPR5 = 0
For i = 1 To 6 GammaRPPR5 = GammaRPPR5 + N_R5(i) * I_R5(i) * (I_R5(i) - 1) * (P /
Pred) ^ (I_R5(i) - 2) * (Tred / T) ^ J_R5(i) Next i
End Function

Private Function GammaRTTR5(T As Double, P As Double, Tred As Integer,
Pred As Double)
Call CoefR5
GammaRTTR5 = 0
For i = 1 To 6 GammaRTTR5 = GammaRTTR5 + N_R5(i) * (P / Pred) ^ I_R5(i) *
J_R5(i) * (J_R5(i) - 1) * (Tred / T) ^ (J_R5(i) - 2) Next i
End Function

Private Function GammaRPTR5(T As Double, P As Double, Tred As Integer,
Pred As Double)
Call CoefR5
GammaRPTR5 = 0
For i = 1 To 6 GammaRPTR5 = GammaRPTR5 + N_R5(i) * I_R5(i) * (P / Pred) ^
(I_R5(i) - 1) * J_R5(i) * (Tred / T) ^ (J_R5(i) - 1) 137 Next i
End Function

' ' Numerical values of the coefficients and exponents of the dimensionless Gibbs
free energy for region 5
Private Sub CoefR5()
'ideal-gas part GAMMA0, Eq. (33) (Table 37, page 36)
J0_R5(1) = 0
J0_R5(2) = 1
J0_R5(3) = -3
J0_R5(4) = -2
J0_R5(5) = -1
J0_R5(6) = 2

N0_R5(1) = -13.179983674201
N0_R5(2) = 6.8540841634434
N0_R5(3) = -0.024805148933466
N0_R5(4) = 0.36901534980333
N0_R5(5) = -3.1161318213925
N0_R5(6) = -0.32961626538917

'residual part GAMMAR, Eq. (34) (Table 38, page 37)
I_R5(1) = 1
I_R5(2) = 1
I_R5(3) = 1
I_R5(4) = 2
I_R5(5) = 2
I_R5(6) = 3

J_R5(1) = 1
J_R5(2) = 2
J_R5(3) = 3
J_R5(4) = 3
J_R5(5) = 9
J_R5(6) = 7

N_R5(1) = 1.5736404855259E-03
N_R5(2) = 9.0153761673944E-04
N_R5(3) = -5.0270077677648E-03
N_R5(4) = 2.2440037409485E-06
N_R5(5) = -4.1163275453471E-06
N_R5(6) = 3.7919454822955E-08
End Sub
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' END REGION 5 ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
138
Allegato 6: Modulo ''Gas'.

' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' DECLARATIONS ' ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '

' ' Reference for constanst: NIST webbook ''
' ' http://webbook.nist.gov/ ''

Private A As Double, b As Double, c As Double, d As Double, e As Double, f As
Double, g As Double, h As Double, MM As Double
Private Const Ru As Double = 8.314 '[kJ/kmolK] Universal gas constant
Private Const Trif As Double = 298.15 '[K] Standard temperature
Private Const Prif As Double = 0.101325 '[MPa] Standard pressure

Option Explicit
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' END DECLARATIONS ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '

' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' GATEWAY FUNCTION ''
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' Gateway function to call (ideal) gas property functions
Public Function LMC_gas(species As String, prop As String, base As String,
Optional T As Double, Optional P As Double) As Variant
'species indicate the selected species:
'prop indicate the property:
'base indicate the selected base: mole or mass
'T is the temperature in [K]
'P is the pressure in MPa

' Local declaration
Dim X As Double ' numerical results from property functions Dim msg As String ' string result in case of error message Dim fOK As Integer ' flag for occuring errors: 0 an error occured, 1 all fine Dim fBS As Integer ' flag for changeable base: 0 not changeable, 1 changeable from mass to mole
' initialization
X = 0
fOK = 1
fBS = 0
msg = ""
139 If LCase(prop) = "mm" Or LCase(prop) = "hhv" Or LCase(prop) = "lhv"
Then T = Trif
End If

' retrieve gas parameters
If fOK = 1 Then
Select Case LCase(species)
Case "argon"
If T >= 273.15 And T < 6000 Then
Call CoefAr_R1
End If
Case "carbondioxide"
If T >= 273.15 And T < 1200 Then
Call CoefCO2_R1
ElseIf T >= 1200 And T < 6000 Then
Call CoefCO2_R2
End If
Case "carbonmonoxide"
If T >= 273.15 And T < 1300 Then
Call CoefCO_R1
ElseIf T >= 1300 And T < 6000 Then
Call CoefCO_R2
End If
Case "hydrogen"
If T >= 273.15 And T < 1000 Then
Call CoefH2_R1
ElseIf T >= 1000 And T < 2500 Then
Call CoefH2_R2
ElseIf T >= 2500 And T <= 6000 Then
Call CoefH2_R3
End If
Case "methane"
If T >= 273.15 And T < 1200 Then
Call CoefCH4_R1
ElseIf T >= 1200 And T < 6000 Then
Call CoefCH4_R2
End If
Case "nitrogen"
If T >= 100 And T < 500 Then
Call CoefN2_R1
ElseIf T >= 500 And T < 2000 Then
Call CoefN2_R2
ElseIf T >= 2000 And T <= 6000 Then
Call CoefN2_R3
End If
Case "oxygen" 140
If T >= 100 And T < 700 Then
Call CoefO2_R1
ElseIf T >= 700 And T < 2000 Then
Call CoefO2_R2
ElseIf T >= 2000 And T <= 6000 Then
Call CoefO2_R3
End If
Case "ethane"
If T >= 200 And T < 700 Then
Call CoefC2H6
End If
Case "propane"
If T >= 200 And T < 700 Then
Call CoefC3H8
End If
Case "isobutane"
If T >= 200 And T < 700 Then
Call CoefiC4H10
End If
Case "normalbutane"
If T >= 200 And T < 700 Then
Call CoefnC4H10
End If
Case "isopentane"
If T >= 200 And T < 700 Then
Call CoefiC5H12
End If
Case "normalpentane"
If T >= 200 And T < 700 Then
Call CoefnC5H12
End If
Case "hexane"
If T >= 200 And T < 700 Then
Call CoefnC6H14
End If
Case Else
fOK = 0
msg = "Undefined species"
End Select
End If

' choose property to be evaluated
If fOK = 1 Then
Select Case LCase(prop)
Case "cp" 'J/mol = kJ/kmol
X = fcp(T / 1000)
fBS = 1 141 Case "h" 'kJ/mol trasformed in kJ/kmol multiplying by 1000
X = fh(T / 1000) * 1000 '
fBS = 1
Case "s" 'J/mol K = kJ/kmol K
X = fs(T / 1000, P)
fBS = 1
Case "mm" ' kg/kmol
X = MM
Case "lhv" 'J/mol = kJ/kmol
X = fHV(species, prop)
fBS = 1
Case "hhv" 'J/mol = kJ/kmol
X = fHV(species, prop)
fBS = 1
Case Else
fOK = 0
msg = "Undefined property"
End Select
End If

' change basis
If fOK = 1 Then
Select Case LCase(base)
Case "mass"
' change specific properties
If fBS = 1 Then
X = X / MM
End If
Case "mole"
' nothing to do
Case Else
fOK = 0
msg = "Undefined basis"
End Select
End If


' fOK check
If fOK = 1 Then
LMC_gas = X
Else
LMC_gas = msg
End If
End Function
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' END GATEWAY FUNCTION ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' 142
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' ANCILLARY PROPERTY FUNCTIONS ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '

' Gas Phase Heat Capacity (Shomate Equation from NIST Webbook)

Private Function fcp(T As Double)
' Return the Cp(T) value of a gas
fcp = A + b * T + c * T ^ 2 + d * T ^ 3 + e / T ^ 2 'Cp in kJ/kmol
End Function

Private Function fh(T As Double)
' Return the h(T) value of a gas
fh = A * T + b / 2 * T ^ 2 + c / 3 * T ^ 3 + d / 4 * T ^ 4 - e / T + f - h 'Enthalpy
in kJ/mol
End Function

Private Function fs(T As Double, P As Double)
' Return the s(T) value of a gas
fs = A * Log(T) + b * T + c / 2 * T ^ 2 + d / 3 * T ^ 3 + e / 2 / T ^ 2 + g - Ru *
Log(P / Prif) 'Entropy in kJ/kmol K
End Function

Private Function fHV(species As String, prop As String) As Double 'LHV or
HHV of the selected species in kJ/Kmol
' declarations
Dim h_water As Double
Dim h_species As Double
Dim h_carbondioxide As Double
Dim P_water As Double
Dim a_water, a_carbondioxide As Integer

' water enthalpy, kJ/kmol at formation conditions
Select Case LCase(prop)
Case "lhv" ' compute saturate vapor
P_water = LMC_water("PSat", "mole", Trif) * (1 - 0.000001)
Case "hhv" ' compute undercooled liquid
P_water = Prif
End Select
h_water = LMC_water("h", "mole", Trif, P_water, "formation")

' species enthalpy, kJ/kmol
h_species = LMC_gas(species, "h", "mole", Trif, Prif)

' carbon dioxide, kJ/kmol
h_carbondioxide = LMC_gas("carbondioxide", "h", "mole", Trif, Prif)
143 'stechiometric coefficients
Select Case LCase(species)
Case "hydrogen"
a_water = 1
a_carbondioxide = 0
Case "methane"
a_water = 2
a_carbondioxide = 1
Case "ethane"
a_water = 3
a_carbondioxide = 2
Case "propane"
a_water = 4
a_carbondioxide = 3
Case "isobutane", "normalbutane"
a_water = 5
a_carbondioxide = 4
Case "isopentane", "normalpentane"
a_water = 6
a_carbondioxide = 5
Case "hexane"
a_water = 7
a_carbondioxide = 6
Case "carbonmonoxide"
a_water = 0
a_carbondioxide = 1
Case Else
h_species = 0
a_water = 0
a_carbondioxide = 0
End Select
fHV = h_species - a_water * h_water - a_carbondioxide * h_carbondioxide
End Function
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' END AUXILIARY PROPERTY FUNCTIONS ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '

' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' ' ' ' COEFFICENT '' ' ' ' ' taken from WebBook, NIST
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
'Argon
Private Sub CoefAr_R1()
A = 20.786
b = 0.0000002825911
c = -0.0000001464191
d = 0.00000001092131
e = -0.00000003661371 144
f = -6.19735
g = 179.999
h = 0
MM = 39.948
End Sub
'Carbondioxide
Private Sub CoefCO2_R1()
A = 24.99735
b = 55.18696
c = -33.69137
d = 7.94839
e = -0.13664
f = -403.6075
g = 228.2431
h = 0
MM = 44.0095
End Sub

Private Sub CoefCO2_R2()
A = 58.16639
b = 2.72007
c = -0.49229
d = 0.03884
e = -6.44729
f = -425.9186
g = 263.6125
h = 0
MM = 44.0095
End Sub
'Carbonmonoxide
Private Sub CoefCO_R1()
A = 25.56759
b = 6.09613
c = 4.054656
d = -2.671301
e = 0.131021
f = -118.0089
g = 227.3665
h = 0
MM = 28.0101
End Sub

Private Sub CoefCO_R2()
A = 35.1507
b = 1.300095
c = -0.205921
d = 0.01355 145 e = -3.28278
f = -127.8375
g = 231.712
h = 0
MM = 28.0101
End Sub
'Hydrogen
Private Sub CoefH2_R1()
A = 33.06618
b = -11.36342
c = 11.43282
d = -2.77287
e = -0.15856
f = -9.9808
g = 172.70979
h = 0
MM = 2.01588
End Sub

Private Sub CoefH2_R2()
A = 18.56308
b = 12.25736
c = -2.85979
d = 0.26824
e = 1.97799
f = -1.14744
g = 156.28813
h = 0
MM = 2.01588
End Sub

Private Sub CoefH2_R3()
A = 43.41356
b = -4.293079
c = 1.272428
d = -0.096876
e = -20.533862
f = -38.515158
g = 162.081354
h = 0
MM = 2.01588
End Sub
'Methane
Private Sub CoefCH4_R1()
A = -0.703029
b = 108.4773
c = -42.52157 146
d = 5.862788
e = 0.678565
f = -76.84376
g = 158.7163
h = 0
MM = 16.0425
End Sub

Private Sub CoefCH4_R2()
A = 85.81217
b = 11.26467
c = -2.114146
d = 0.13819
e = -26.42221
f = -153.5327
g = 224.4143
h = 0
MM = 16.0425
End Sub
'Nitrogen
Private Sub CoefN2_R1()
A = 28.98641
b = 1.853978
c = -9.647459
d = 16.63537
e = 0.000117
f = -8.671914
g = 226.4168
h = 0
MM = 28.0134
End Sub

Private Sub CoefN2_R2()
A = 19.50583
b = 19.88705
c = -8.598535
d = 1.369784
e = 0.527601
f = -4.935202
g = 212.39
h = 0
MM = 28.0134
End Sub

Private Sub CoefN2_R3()
A = 35.51872
b = 1.128728 147 c = -0.196103
d = 0.014662
e = -4.55376
f = -18.97091
g = 224.981
h = 0
MM = 28.0134
End Sub
'Oxygen
Private Sub CoefO2_R1()
A = 31.32234
b = -20.23531
c = 57.86644
d = -36.50624
e = -0.007374
f = -8.903471
g = 246.7945
h = 0
MM = 31.9988
End Sub

Private Sub CoefO2_R2()
A = 30.03235
b = 8.772972
c = -3.988133
d = 0.788313
e = -0.741599
f = -11.32468
g = 236.1663
h = 0
MM = 31.9988
End Sub

Private Sub CoefO2_R3()
A = 20.91111
b = 10.72071
c = -2.020498
d = 0.146449
e = 9.245722
f = 5.337651
g = 237.6185
h = 0
MM = 31.9988
End Sub
'Ethane
Private Sub CoefC2H6()
A = 3.045026 * 10 148
b = 1.393816 * 10
c = 2.649906 * 10 ^ 2
d = -2.070447 * 10 ^ 2
e = 0
f = -95.4302974723231
g = 0
h = 0
MM = 30.069
End Sub
'Propane
Private Sub CoefC3H8()
A = 2.845289 * 10
b = 9.049519 * 10
c = 2.826678 * 10 ^ 2
d = -2.560307 * 10 ^ 2
e = 0
f = -119.196887132583
g = 0
h = 0
MM = 44.0956
End Sub
'i-Butane
Private Sub CoefiC4H10()
A = 2.897929 * 10
b = 1.636562 * 10 ^ 2
c = 3.06985 * 10 ^ 2
d = -3.087608 * 10 ^ 2
e = 0
f = -152.216261634091
g = 0
h = 0
MM = 58.1222
End Sub
'n-Butane
Private Sub CoefnC4H10()
A = 4.358927 * 10
b = 9.459586 * 10
c = 4.118574 * 10 ^ 2
d = -3.640003 * 10 ^ 2
e = 0
f = -145.720093789653
g = 0
h = 0
MM = 58.1222
End Sub
'i-Isopentane
Private Sub CoefiC5H12() 149 A = 1.368775 * 10
b = 3.52611 * 10 ^ 2
c = 4.499795 * 10
d = -1.44554 * 10 ^ 2
e = 0
f = -173.56536956816
g = 0
h = 0
MM = 72.1488
End Sub
'n-Isopentane
Private Sub CoefnC5H12()
A = 5.791811 * 10
b = 7.777861 * 10
c = 5.898819 * 10 ^ 2
d = -4.966314 * 10 ^ 2
e = 0
f = -171.755513729507
g = 0
h = 0
MM = 72.1488
End Sub
'hexane
Private Sub CoefnC6H14()
A = 6.709984 * 10
b = 9.741664 * 10
c = 7.083907 * 10 ^ 2
d = -6.085141 * 10 ^ 2
e = 0
f = -196.491835883083
g = 0
h = 0
MM = 86.1754
End Sub

' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' ' ' END COEFFICENT ' ' ' '
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
150
Allegato 7: Modulo ''Mixture'.

' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
' ' GATEWAY FUNCTION ''
' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' ' '
Option Explicit

Public Function LMC_mixture(ispecies As Range, ifrac As Range, prop As
String, base As String, Optional T As Double, Optional P As Double) As Variant
' Gateway function to call mixture property functions

' Declaration
Dim X As Double ' numerical results from property functions
Dim msg As String ' string result in case of error message Dim fOK As Integer ' flag for occuring errors: 0 an error occured, 1
all fine Dim fBS As Integer ' flag for changeable base: 0 not changeable, 1 changeable from mass to mole Dim pWT As Integer ' position of water if present: 0 not present, otherwise position within species Dim species() As String ' array for species names Dim zfrac() As Double ' array of overall mole fractions Dim xfrac() As Double ' array of liquid mole fractions Dim yfrac() As Double ' array of vapor mole fractions Dim lfrac As Double ' liquid mole fraction Dim vfrac As Double ' vapor mole fraction Dim wfrmx As Double ' max water vapor mole fraction Dim toll As Double ' tolerance
' initialization
X = 0
fOK = 1
fBS = 0
pWT = 0
msg = ""
lfrac = 0
vfrac = 0
wfrmx = 0
toll = 0.000001

' create and check species and zfrac arrays. Manipulated variables: fOK, msg, species, zfrac, xfrac, yfrac, pWT Call spiecefrac(fOK, msg, ispecies, ifrac, species, zfrac, xfrac, yfrac, pWT, toll) If LCase(prop) = "lfrac" Or LCase(prop) = "vfrac" Or LCase(prop) = "cp" Or LCase(prop) = "h" Or LCase(prop) = "s" Or LCase(prop) = "wfrmx" Then 151 ' computes vapor-liquid equilibrium. Manipulated variables: lfrac, vfrac, wfrmx, xfrac, yfrac Call VLE(fOK, msg, pWT, zfrac, lfrac, vfrac, wfrmx, xfrac, yfrac, T, P, toll) End If

' choose property to be evaluated
If fOK = 1 Then
Select Case LCase(prop)
Case "lfrac" 'kmol/kmol
X = lfrac
fBS = 0 ' think on how to change from mole to mass fractions
Case "vfrac" 'kmol/kmol
X = vfrac
fBS = 0 ' think on how to change from mole to mass fractions
Case "cp", "h", "s" 'J/mol = kJ/kmol
X = f_mixture1(species, prop, xfrac, yfrac, lfrac, vfrac, T, P)
fBS = 1
Case "mm" ' kg/kmol
X = f_mixture2(species, prop, zfrac)
Case "lhv", "hhv"
X = f_mixture2(species, prop, zfrac)
fBS = 1
Case "wfrmx"
X = wfrmx
fBS = 1
Case Else
fOK = 0
msg = "Undefined property"
End Select
End If

' change basis
If fOK = 1 Then
Select Case LCase(base)
Case "mass"
' change specific properties
If fBS = 1 Then
'X = X / MM
End If
Case "mole"
' nothing to do
Case Else
fOK = 0
msg = "Errore in LMC_mixture: Undefined basis"
End Select
End If
152
' fOK check
If fOK = 1 Then
LMC_mixture = X
Else
LMC_mixture = msg
End If
End Function

Private Sub VLE(fOK As Integer, msg As String, pWT As Integer, zfrac() As Double, lfrac As Double, vfrac As Double, wfrmx As Double, xfrac() As
Double, yfrac() As Double, T As Double, P As Double, toll As Double) ' computes vapor-liquid equilibrium
' input: pWT, zfrac, T, P, toll
' output: msg, lfrac, vfrac, wfrmx, xfrac, yfrac
' input/output: fOK,

' declarations
Dim i As Integer ' counter
Dim wfrac As Double ' water overall fraction, if present
Dim wPsat As Double ' water saturation pressure, if present
Dim chk As Double ' check variable

' Vapor-Liquid Equilibrium (VLE)
If fOK = 1 Then
If pWT = 0 Then
' no water present, all vapor
lfrac = 0
vfrac = 1
wfrmx = 0
yfrac = zfrac
Else
' water is present
wfrac = zfrac(pWT) ' water overall mole fraction wPsat = LMC_water("Psat", "mole", T) ' water saturation pressure
' check if other components are present
If wfrac < 1 Then ' other components are present, so all vapor or liquid-vapor states are possible while all liquid is impossible wfrmx = wPsat / P
' check maximum water vapor fraction
If wfrac <= wfrmx Then
' all vapor
lfrac = 0
vfrac = 1
yfrac = zfrac
Else 153 ' vapor-liquid
lfrac = (wfrac - wfrmx) / (1 - wfrmx)
vfrac = 1 - lfrac
xfrac(pWT) = 1
For i = 1 To UBound(zfrac)
yfrac(i) = (zfrac(i) - lfrac * xfrac(i)) / (1 - lfrac)
Next
End If
ElseIf wfrac = 1 Then
' only water is present
If P < wPsat Then
' all vapor
lfrac = 0
vfrac = 1
yfrac = zfrac
ElseIf P > wPsat Then
' all liquid
lfrac = 1
vfrac = 0
xfrac = zfrac
Else
fOK = 0 msg = "Error in VLE: only water is present and at saturation, so the state is undefined given T and P" End If
Else
fOK = 0 msg = "Error in VLE: water overall mole fraction greater than 1." End If

End If
End If

' check on sum of lfrac+vfrac and also on xfrac and yfrac
If fOK = 1 Then
' check sum of lfrac and vfrac
If Abs(lfrac + vfrac - 1) > toll Then
fOK = 0
msg = "Error in VLE: lfrac and vfrac do not sum to 1"
End If

' check sum of xfrac if liquid present
If lfrac > 0 Then
chk = 0
For i = 1 To UBound(xfrac)
chk = chk + xfrac(i)
Next 154
If Abs(chk - 1) > toll Then
fOK = 0
msg = "Error in VLE: xfrac do not sum to 1"
End If
End If

' check sum of yfrac if vapor is present
If vfrac > 0 Then
chk = 0
For i = 1 To UBound(yfrac)
chk = chk + yfrac(i)
Next
If Abs(chk - 1) > toll Then
fOK = 0
msg = "Error in VLE: yfrac do not sum to 1"
End If
End If

' check sum of lfrac*xfrac(i) + vfrac*yfrac(i) is equal to zfrac(i)
If fOK = 1 Then
For i = 1 To UBound(yfrac)
If Abs(zfrac(i) - lfrac * xfrac(i) - vfrac * yfrac(i)) > toll Then
fOK = 0 msg = "Error in VLE: sum of liquid and vapor fraction of each phase is different to overall species mole" End If
Next
End If

End If

End Sub

Private Function f_mixture1(species() As String, prop As String, xfrac() As Double, yfrac() As Double, lfrac As Double, vfrac As Double, T As Double, P
As Double) ' compute the value of "cp", "h", "s" of a mixture
' input: species, prop, xfrac, yfrac, lfrac, vfrac, T and P
' input/output:
' output: f_mixture1

' Declarations
Dim i As Integer ' counter
Dim lprop As Double ' liquid phase property
Dim vprop As Double ' vapor phase property
Dim Ppart As Double ' species partial pressure
Dim basis As String ' basis for the properties 155 Dim refer As String ' property reference for water

' init
basis = "mole" ' all calculations must be in molar basis
refer = "formation" ' formation property in use for water
lprop = 0
vprop = 0

' loop over species
For i = 1 To UBound(species)
' vapor phase partial pressure
Ppart = yfrac(i) * P
If species(i) = "water" Then
' liquid phase
If xfrac(i) > 0 Then
lprop = lprop + xfrac(i) * LMC_water(prop, basis, T, P, refer)
End If
' vapor phase
If yfrac(i) > 0 Then vprop = vprop + yfrac(i) * LMC_water(prop, basis, T, Ppart, refer) End If
Else
' species only in the vapor phase
If yfrac(i) > 0 Then vprop = vprop + yfrac(i) * LMC_gas(species(i), prop, basis, T, Ppart) End If
End If
Next
' mixture property as weighted balance between liquid and vapor phases f_mixture1 = lfrac * lprop + vfrac * vprop
End Function

Private Function f_mixture2(species() As String, prop As String, zfrac() As Double)
' compute the value of "mm", "lhv", "hhv" of a mixture
' input: species, prop, xfrac and zfrac
' input/output:
' output: f_mixture2

Dim i As Integer ' counter
Dim base As String

base = "mole"
f_mixture2 = 0
156
For i = 1 To UBound(species)
If zfrac(i) > 0 Then
If species(i) = "water" Then
f_mixture2 = f_mixture2 + zfrac(i) * LMC_water(prop, base)
Else f_mixture2 = f_mixture2 + zfrac(i) * LMC_gas(species(i), prop, base) End If
End If
Next
End Function

Private Sub spiecefrac(fOK As Integer, msg As String, ispecies As Range, ifrac
As Range, species() As String, zfrac() As Double, xfrac() As Double, yfrac() As
Double, pWT As Integer, toll As Double)

' creates and checks species and zfrac arrays
' input: fOK, msg, ispecies, ifrac, prop, xfrac and zfrac
' input/output: species, zfrac, pwt, xfrac, yfrac
' output:

' Declarations
Dim i As Integer ' counter
Dim j As Integer ' counter
Dim index As Integer ' array index
Dim nrow_ispecies As Integer ' auxiliary
Dim ncol_ispecies As Integer ' auxiliary
Dim nrow_ifrac As Integer ' auxiliary
Dim ncol_ifrac As Integer ' auxiliary
Dim fchk As Double ' check sum of zfrac is 1

' dimensions of ispecies and ifrac
If fOK = 1 Then
nrow_ispecies = ispecies.Rows.Count
ncol_ispecies = ispecies.Columns.Count
ncol_ifrac = ifrac.Columns.Count
nrow_ifrac = ifrac.Rows.Count

' check on ispecies and ifrac dimensions If (nrow_ispecies <> nrow_ifrac) Or (ncol_ispecies <> ncol_ifrac) Then
fOK = 0
msg = "ispecies and ifrac have different dimensions"
End If End If

' create and check spieces and zfrac arrays
If fOK = 1 Then
'Creation of spieces array 157 ReDim species(1 To nrow_ispecies * ncol_ispecies) As String
ReDim zfrac(1 To nrow_ifrac * ncol_ifrac) As Double
ReDim xfrac(1 To nrow_ifrac * ncol_ifrac) As Double
ReDim yfrac(1 To nrow_ifrac * ncol_ifrac) As Double

' init
index = 0
fchk = 0
' loop over species
For i = 1 To nrow_ispecies
For j = 1 To ncol_ispecies
index = index + 1
species(index) = ispecies.Cells(i, j).Value
zfrac(index) = ifrac.Cells(i, j).Value
xfrac(index) = 0
yfrac(index) = 0
'check water presence
If LCase(species(index)) = "water" Then
pWT = index
End If
' compute fraction sum
fchk = fchk + ifrac.Cells(i, j).Value
Next
Next

'check on sum of mole fractions
If Abs(fchk - 1) > toll Then
fOK = 0
msg = "sum of molar fractions must be equal to one"
End If
End If
End Sub


© Eiom - All rights Reserved     P.IVA 00850640186